Skip to main content
BMC is moving to Springer Nature Link. Visit this journal in its new home.
Download PDF

WUSCHEL-dependent chromatin regulation in maize inflorescence development at single-cell resolution

Genome Biology volume 26, Article number: 329 (2025) Cite this article

Abstract

Background

See More

WUSCHEL (WUS) is a homeodomain transcription factor vital for stem cell proliferation in plant meristems. In maize, ZmWUS1 is expressed in the inflorescence meristem, including the central zone reservoir of stem cells. ZmWUS1 overexpression in the Barren inflorescence3 (Bif3) mutant perturbs inflorescence development due to stem cell over-proliferation.

Results

Single-cell Assay for Transposase Accessible Chromatin sequencing (scATAC-seq) shows that Bif3 alters central zone chromatin accessibility compared to normal inflorescence meristems. The CAATAATGC motif, a known homeodomain recognition site, is enriched within regions with increased chromatin accessibility in Bif3, suggesting ZmWUS1 could function as a transcriptional activator in the central zone. This motif differs from the TGAATGAA motif identified by DNA Affinity Purification sequencing (DAP-seq) of ZmWUS1, which showed low enrichment in the central zone. Conversely, regions with decreased chromatin accessibility in Bif3 are instead adjacent to AUXIN RESPONSE FACTOR genes, suggesting possible reduced auxin signaling in the Bif3 central zone.

Conclusions

This study characterized how Bif3 overexpression of ZmWUS1 influences chromatin accessibility in the central zone, reducing auxin signaling, while raising questions about differential ZmWUS1 motif usage in distinct cellular contexts.

Peer Review reports

Background

In inflorescence meristems, the delicate balance between self-renewal and stem cell differentiation is controlled by the transcription factor (TF) WUSCHEL (WUS) [1]. Stem cell maintenance involves a negative feedback loop where WUS promotes CLAVATA3 (CLV3) expression, and CLV3 represses WUS [2]. In Arabidopsis, mutant studies revealed the functional interplay between these two components: in clv3, WUS expression expands, whereas in wus, CLV3 expression markedly decreases [3]. This feedback loop ensures a balanced stem cell population, where an overabundance of WUS increases stem cells and its scarcity depletes them [4]. The CLV-WUS regulatory pathway is conserved among plant species, including tomato, maize, and rice [5,6,7], even though WUS function is not [8, 9].

In maize, two paralogs, ZmWUS1 and ZmWUS2, are co-orthologs of Arabidopsis WUS [10]. ZmWUS1 is expressed within the female inflorescence meristems, which eventually give rise to mature ears [11]. ZmWUS1 is expressed in both the apical inflorescence meristem and the reproductive axillary meristems (spikelet-pair meristems, spikelet meristems, and floral meristems) that eventually form kernels [12]. In Barren inflorescence3 (Bif3) mutants, a tandem ZmWUS1 duplication perturbs the CLV3-WUS feedback loop and elevates ZmWUS1 expression [11, 13]. Bif3 ZmWUS1 overexpression drastically rearranges the functional domains of ear inflorescence meristems, causing small, spherical ears with a reduced kernel count [11].

In maize, similar to Arabidopsis, a functional homolog of CLV3, ZmCLE7, is implicated in the regulation of inflorescence stem cell proliferation [14]. Diminished ZmCLE7 expression increases kernel count, producing kernels that are narrower and less round, whereas ZmCLE7 upregulation lowers kernel yield in ears [15]. In Arabidopsis, CLV3 is expressed in the three uppermost cell layers, L1, L2, and L3, where the stem cells are located [16]. WUS is expressed in cells starting from the L3 layer and extending to more inner layers; the region of WUS expression is known as the organizing center [17]. Although WUS is expressed beneath the stem cell niche, WUS protein accumulates across a wider area including CLV3 expressing stem cells [18]. It is posited that WUS diffuses to activate CLV3 in these stem cells, as evidenced by the observation that a decrease in WUS expression leads to a reduction in CLV3 expression [19]. However, CLV3 is not co-expressed with WUS in the organizing center, and some have suggested that the WUS transcriptional regulation may vary by cellular context; in the central zone, WUS promotes CLV3 expression, but in the organizing center, it could repress CLV3 [20]. HAIRY MERISTEMs (HAMs) are GRAS-family TFs that function as transcriptional co-factors, and genetic interactions suggest a model where in the organizing center WUS-HAM complexes inhibit CLV3 transcription [20]. In contrast, another model proposes that WUS forms homodimers that repress CLV3 in the organizing center, but, as WUS diffuses to nearby cells, its concentration diminishes and favoring WUS monomer formation and CLV3 activation [21]. These findings collectively suggest that WUS’s function depends on cellular context, emphasizing the need for cell-specific data to fully understand its regulatory roles.

WUS is a homeodomain TF, with a helix–turn–helix DNA-binding domain that forms hydrogen bonds with specific DNA sequences called “motifs” [17]. There have been significant efforts to identify WUS-dependent cis-regulation at CLV3 and other targets. For example, Electrophoretic Mobility Shift Assays (EMSA) and Chromatin Immunoprecipitation followed by quantitative PCR (ChIP-qPCR) have detected binding of WUS to TAAT motifs both upstream and downstream of CLV3 [18, 21]. This TAAT motif is bound by other homeodomain TFs [17, 22]. WUS and SHOOT MERISTEMLESS (STM), another stem-cell regulator, dimerize and interact with the TGACA motif upstream of CLV3 [23]. Additionally, ChIP-qPCR has revealed WUS binding upstream of the ARABIDOPSIS RESPONSE REGULATOR 7 (ARR7) [24], whereas EMSA has identified WUS binding downstream of the AGAMOUS (AG), corresponding to a CC(A/T)6GG [25, 26] or TTAATGG motif [27]. Genomic approaches have also identified WUS targets. For example, Arabidopsis ChIP-seq using inducible WUS and WUS DNA affinity purification sequencing (DAP-seq) [28] identified the TGAATGAA motif [28]. Microscale thermophoresis has demonstrated an affinity preference for the TGAATGAA motif over TAAT [29]. Likewise, the crystal structure of the WUS homodimer suggests it preferentially binds a TGAATGAA motif [29]. These diverse, and occasionally conflicting, WUS binding motifs reflect the hierarchical nature of protein:DNA interactions, suggesting that some motifs are favored in certain cellular contexts. These molecular and biochemical data highlight that WUS binding, and thus downstream activity, likely varies based on cell-type-specific interacting partners, driving differential WUS cis-regulatory element binding preferences and transcriptional regulation.

Auxin related genes are candidate WUS targets supported by Arabidopsis ChIP-seq [30]. Auxin is a plant hormone that plays a concentration dependent role in inducing differentiation or maintaining stem cells [31]. Auxin signaling is mediated by AUXIN RESPONSE FACTOR (ARF) family TFs [32]. ARFs are divided into clade A and clade B, that, respectively, activate and repress transcription some in an auxin-dependent fashion [32,33,34]; a third clade, clade C, exists but its function is less clear. WUS expressing stem cells have low auxin concentrations, inhibiting differentiation [35]. However, when auxin levels rise, Aux/IAA proteins are degraded, de-repressing ARF activity and causing auxin transcriptional reprogramming [36], promoting meristematic cell differentiation. Despite the parallel role of auxin in stem-cell maintenance, the extent to which WUS controls meristematic ARF activity, and auxin signaling, remains poorly understood.

Although the significance of ZmWUS1 in maize inflorescence development is known [11], the detailed mechanisms, and particularly the identity of ZmWUS1 targets, are not. The complex regulatory dynamics of ZmWUS1, and its unique cell-type-specific expression pattern, complicate target gene identification. To accurately identify cis-regulatory elements influenced by ZmWUS1 at single-cell resolution, we used single-cell Assay for Transposase-Accessible Chromatin with sequencing (scATAC-seq) comparing wild-type (WT) and Bif3 mutant immature maize ears. Contrasting cell-type resolved chromatin accessibility between WT and Bif3 developing female inflorescences identified cell-type-specific mechanisms of ZmWUS1 cis-regulation. This Bif3-to-WT comparison suggests ZmWUS1 directly promotes transcription via central zone specific cis-regulation, while indirectly repressing stem-cell auxin signaling. This research supports the model whereby ZmWUS1 targets distinct cis-regulatory elements depending on the cell type.

Results

scATAC-seq captures central zone nuclei in the immature maize ear

We hypothesized that Bif3 ZmWUS1 overexpression alters cis-regulatory element activity, and the accessible chromatin landscape, with ZmWUS1 binding cis-regulatory elements that are not targeted in WT. To assess this, we compared the chromatin accessibility landscape across cell types by performing scATAC-seq, with two biological replicates, on developing WT and Bif3 female inflorescences. After filtering nuclei based on Tn5 insertion number, transcription start site enrichment, organellar ratio, Fraction of Reads in Peak (FRiP) scores, and doublet scores, an average of 4426 nuclei per replicate remained, each with an average of ~ 25,000 Tn5 insertions. Separate cell annotations for WT and Bif3 identified 13 and 14 distinct nuclei clusters, respectively, with sub-clustering further refining cell-type annotations (Additional file 1: Fig. S2A; Additional file 1: Fig. S2B; Additional file 2: Table S1-3). Cluster cell-type annotation was determined using marker gene body chromatin accessibility, an established proxy for gene transcription [37] (Fig. 1A).

Fig. 1
figure 1

Annotation of central zone nuclei with chromatin accessibility in WT and Bif3. A UMAP plot for WT and Bif3, colored to represent the 14 cell types. B Illustrations depicting the phenotypes of WT and Bif3, including longitudinal sections of the inflorescence meristem, the spikelet pair meristem, the spikelet meristem, and the floral meristem are colored to correspond with the cell types in the UMAP. Grey colors represent cell types with an unknown annotation. C–G UMAP plots highlighting nuclei exhibiting high gene body chromatin accessibility, depicted in purple, surrounding the marker genes of ZmWUS1 (C), ZmCLE7 (D), ZmOCL4 (E), ZYB15 (F), and RA1 (G). H Genome browser track displaying chromatin accessibility around marker genes for meristem-associated cells in WT and Bif3. The ranges of number of indicate CPM (counts per million) normalized Tn5 insertion number. I Gene body chromatin accessibility patterns of 30 marker genes in WT and Bif3. The dot size represents the percentage of nuclei that have chromatin accessibility around the marker genes. The values are Z-score normalized aggregated number of Tn5 insertions by cell types

Full size image

Unlike Arabidopsis WUS, which has organizing center exclusive expression [38], the expression domain of ZmWUS1 (Zm00001eb067310) in maize female inflorescence spans a broader area than the organizing center [11], covering multiple cell layers 1 through 10 across the meristem (Fig. 1B). Further distinguishing maize from Arabidopsis, the inflorescence central zone stem cells, which express the CLV3 homolog ZmCLE7, also express ZmWUS1. This pattern continues in axillary meristem cells, which include the spikelet meristem, spikelet pair meristem, and floral meristem, with ZmCLE7 and ZmWUS1 exhibiting precisely overlapping expression patterns [11]. In Bif3 mutants, the overlap of ZmWUS1 and ZmCLE7 expression is retained in both axillary and inflorescence meristems. Therefore, we designated nuclei with ZmWUS1 and ZmCLE7 high gene body accessibility as central zone cells (Fig. 1B). ZmWUS1 gene body chromatin accessibility was widespread and enriched in WT Cluster (C) 3–1 & 3–3, and in Bif3 C14 (Fig. 1C). ZmCLE7 exhibited pronounced gene body chromatin accessibility in the clusters that also show the highest gene body chromatin accessibility for ZmWUS1 (Fig. 1D). The other central zone marker gene is LOG7a (Zm00001eb160700), which is highly expressed in the WT central zone, but weaker in Bif3 [11]. Correspondingly, gene body chromatin accessibility for LOG7a was enriched in WT C3-1 & C3-3, but diminished in Bif3 C14 (Additional file 1: Fig. S2C). Next, using the defined clusters, we identified 1,165 genes in the central zone that display enriched gene body chromatin accessibility relative to other cell types (Additional file 2: Table S4).

WT C6 and C11, and Bif3 C8 and C9 were annotated as epidermis, as they exhibited significant gene body chromatin accessibility for OUTER CELL LAYER4 (OCL4, Zm00001eb024680) [39], a TF expressed in the meristem L1 layer and in the epidermis (Fig. 1E). WT C3-2 and C3-4, along with Bif3 C1-1 and C1-3, were annotated as suppressed bract or glume primordia, due to increased gene body chromatin accessibility of ZEA YABBY15 (ZYB15, Zm00001eb248640) [11, 40] (Fig. 1F). WT C10 and Bif3 C1 and C2 were classified as the base of the axillary meristem, as they were characterized by notable gene body chromatin accessibility for RAMOSA1 (RA1, Zm00001eb312340) [41] (Fig. 1G). The remaining clusters were annotated using previously published cell-type-specific markers (Additional file 1: Fig. S3). All clusters were identified in both WT and Bif3, except for a few unannotated (unknown) clusters and nuclei at different stages of the cell cycle (Fig. 2A; Additional file 1: Fig. S1F-J; Additional file 1: Fig. S2D and Additional file 1: Fig.S2E).

Fig. 2
figure 2

The intergenic chromatin accessibility profile of WT and Bif3. A A heatmap depicts the correlation of intergenic ACRs between WT and Bif3. The 2000 most variable intergenic ACRs were selected for analysis. The correlation is based on the aggregated number of Tn5 insertions by cell type between WT and Bif3. B A dot plot displays the number of differential intergenic ACRs across various cell types. The illustration provides examples of differential ACRs with either higher chromatin accessibility in WT or Bif3. The total number of differential ACRs, as well as those exhibiting higher chromatin accessibility in either Bif3 or WT, are distinctly colored. The ratio of differential ACRs is calculated by dividing the number of differential ACRs by the total number of intergenic ACRs

Full size image

Although there are significant morphological differences between WT and Bif3 [11], there was an unexpectedly high consistency in cell-type-specific patterns of chromatin accessibility (Fig. 1H; Additional file 1: Fig. S3). To evaluate the gene body chromatin accessibility, and thus transcription, patterns more comprehensively, we generated cluster-level dot plots for 30 marker genes in WT and Bif3 (Fig. 1I), which revealed similar marker gene chromatin accessibility between WT and Bif3 and concordant cell type annotations. We observed a clear correspondence of de novo markers in gene body chromatin accessibility between WT and Bif3 (Additional file 1: Fig. S4). These findings indicate similarity in the patterns of gene body chromatin accessibility for primary marker genes between WT and Bif3 and corresponding cell types. Histologically, Bif3 mutants exhibited a reduced number of cells in the axillary meristem base, coupled with an increased number in the suppressed bract cells [11]. The distribution of the proportion of annotated cells aligns with previous histological cell-type numbers (Additional file 1: Fig. S5), although some variability exists between replicates.

Intergenic chromatin accessibility specifically differs in the Bif3 central zone

To understand the differences in Accessible Chromatin Regions (ACRs) that underlie the meristem morphological variation, we identified cell-type-level differential ACRs between WT and Bif3. Using the union of all ACRs identified from each cell type, the total number of ACRs was 91,386 in WT and 77,393 in Bif3 (Additional file 2: Table S5), each spanning 500 bp from the peak summit. In WT, approximately 46.6% of the ACRs are gene-overlapping ACRs, 36.6% are distal ACRs—located more than 2 kbp from the closest gene, and 16.7% are proximal ACRs—situated less than 2 kbp from a gene. Similarly, in Bif3, the distribution includes ~ 47.3% gene-overlapping ACRs, ~ 35.6% distal ACRs, and 17% proximal ACRs.

We investigated genome-wide intergenic accessible chromatin changes within WT and Bif3 cell types (Fig. 2A). To capture ACRs with active cis-regulatory regions, we created intergenic ACR sets, defined by ACR midpoint not overlapping with any exons. We aggregated intergenic ACRs across all cell types and genotypes, resulting in 58,578 intergenic ACRs. Indicating high similarity between WT and Bif3, ACR read density correlations between corresponding cell types exceeded 0.95 of Pearson correlation coefficient for all the cell types except the central zone, which displayed a correlation of 0.81. This suggests drastic alterations in chromatin accessibility within the Bif3 central zone, amidst an otherwise conserved chromatin landscape, suggesting that the central zone is the source of the meristem morphological differences observed in Bif3.

To investigate how intergenic chromatin accessibility is changed, we identified differentially ACRs (differential ACRs) between Bif3 and WT samples. For each cell type, we combined intergenic ACRs from both WT and Bif3, covering the full extent of ACRs found in either genotype. The procedure was performed independently for each cell type, resulting in intergenic ACR sets of varying numbers. On average, 40,208 intergenic ACRs were pairwise tested for differential chromatin accessibility by comparing their Tn5 insertion density between WT and Bif3. Most cell types had fewer than 355 differential ACRs, representing less than 1% of the total ACRs per cell type. In contrast, the central zone exhibited 3,194 differential ACRs, accounting for 8.5% of its total ACRs (Fig. 2B). Among the central zone differential ACRs, 1437 increased chromatin accessibility in Bif3 (termed “increased in Bif3”), whereas 1757 decreased chromatin accessibility in Bif3 (termed “decreased in Bif3”) (Additional file 2: Table S5). Consistent with our cell-type correlations, this differential ACR analysis indicates Bif3 impacts the regulatory environment of the central zone, while leaving neighboring cell types relatively unchanged.

Regions with increased chromatin accessibility in Bif3 were enriched for the CAATAATGC motif

To uncover TFs potentially involved in the Bif3 altered chromatin accessibility landscape, we identified enriched motifs using a de novo motif search in differential ACRs with increased and decreased central zone cell accessibility in Bif3 compared to WT samples. The increases in Bif3 ACRs were enriched for a single significant motif (CAATAATGC), whereas the decreases in Bif3 ACRs showed enrichment for five distinct motifs (Fig. 3A; Additional file 2: Table S7). In total, 381 (~ 26.51%) differential ACRs increased in Bif3 had a CAATAATGC motif. The CAATAATGC motif has a core TAAT motif, a previously identified WUS binding motif, making these ACRs good candidates for direct ZmWUS1 binding [17, 22]. To confirm if enrichment of the CAATAATGC motif is specific to the central zone, we performed a known motif search for CAATAATGC across ACRs with increased accessibility in all cell types. In the central zone, 296 (~ 20.6%) of the increased ACRs contained the motif, compared to only 7 (~ 9.3%) on average in other cell types (Additional file 1: Fig. S6). Thus, CAATAATGC preferentially marks altered ACRs in the ZmWUS1 overexpressing central zone.

Fig. 3
figure 3

Characteristics of differentially Accessible Chromatin Regions (differential ACRs) between WT and Bif3 by cell types. A The Position Weight Matrix (PWM) illustrates significant motifs discovered within differential ACRs in the central zone (E-value < 1). The ratio indicates the number of differential ACRs containing the motif divided by the total number of differential ACRs. The color denotes whether the differential ACR sets are increased or decreased in Bif3 in the central zone. B Box plots display the log2 fold change (log2FC) of gene body chromatin accessibility surrounding differential ACRs increased in Bif3 with the CAATAATGC motif. Log2FC was calculated by normalizing gene body chromatin accessibility of differential ACRs in Bif3 to those of WT. Rows represent individual genes, whereas columns represent cell types. Pairwise two-tailed Student’s t-test, ****FDR < 2e-16 and *FDR < 0.05. C A genome browser view illustrates chromatin accessibility in the central zone for the peaks with the CAATAATGC motif

Full size image

Beyond direct ZmWUS1 activities, an alternative explanation for the increased chromatin accessibility of CAATAATGC-containing ACRs is that other homeodomain TFs are altered in the Bif3 central zone as the TAAT motif can be bound by other homeodomain TFs [22]. Many diverse TFs include a homeodomain, including WUSCHEL-related homeobox (WOX) family genes [42,43,44] and they could drive the increases in Bif3 chromatin accessibility. In total, we identified 161 TFs containing homeodomains [45] spanning six families: WOX, homeodomain leucine zipper (HD-ZIP) [46], three-amino-acid-loop-extension (TALE) [47], zinc finger homeodomain (ZF-HD) [48], homeobox (HB)-other [46], and homeobox plant homeodomain finger (HB-PHD) [49]. Among them, 23 showed increased gene body chromatin accessibility in the Bif3 central zone (Additional file 2: Table S6; Additional file 1: Fig. S7B). Other TFs with a homeodomain had decreased or unchanged gene body chromatin accessibility. For example, an ACR near the promoter of ZmWUS2 (Zm00001eb433010) exhibits reduced chromatin accessibility in the Bif3 central zone nuclei (Additional file 1: Fig. S7C) and is downregulated in Bif3 RNA-seq datasets [11]. The 23 TFs with increased accessibility included two WOX members (ZmWUS1 and ZmWOX1, Zm00001eb015500), 13 HD-ZIP, five TALE, one ZF-HD, and two HB-other genes (Additional file 1: Fig. S7D). However, some, such as ZmWOX1, showed fivefold lower in gene body chromatin accessibility compared to ZmWUS1 in the Bif3 central zone, suggesting limited contribution. Only six of the 23 TFs exhibited higher chromatin accessibility than ZmWUS1: five HD-ZIP TFs (Zm00001eb005440, Zm00001eb075230, Zm00001eb319390, Zm00001eb009550, and Zm00001eb402230) and one TALE TF (Zm00001eb202150). Zm00001eb005440 is orthologous to Arabidopsis ATHB20, while Zm00001eb075230 and Zm00001eb319390 are orthologous to ATHB12 [50]. DAP-seq data from Arabidopsis [28] show that ATHB20 and ATHB12 recognize motifs such as (C/T)(C/T)AATAATT(A/G)(A/T) and AAT(C/G)ATT(A/G), respectively. Although ATHB12 shares the TAAT core with the CAATAATGC motif we identified in the central zone, the flanking sequences differ. These HD-ZIP and TALE TFs may thus contribute to the recognition and usage of the CAATAATGC motif, and the observed enrichment could arise indirectly as a downstream consequence of ZmWUS1 overexpression.

We next investigated the change of gene body chromatin accessibility linked to the 381 differential ACRs with increased accessibility in Bif3 in the central zone nuclei with a CAATAATGC motif. We used the fold change of the gene body accessibility in the Bif3 mutant to WT. Most of the genes close to these 381 ACRs increased in central zone gene body accessibility in Bif3 compared to WT (Fig. 3C), displaying an average log2 fold change of ~ 0.38. In contrast, the gene body chromatin accessibility of these genes in other cell types showed smaller changes (Pairwise two-tailed Student’s t-test, FDR < 2e-16), with average log2 fold changes less than 0.05. Of note, we observed that Bif3 increased ACRs containing the CAATAATGC motif were linked to KNOTTED-like homeobox1 (KNOX1, Zm00001eb001720) (Fig. 3C). KNOX1 belongs to the same gene family as KNOTTED1, an important master regulator of maize meristems [51,52,53]. Our findings suggest that KNOX1 could be activated in Bif3 by either ZmWUS1 itself. Furthermore, the scATAC-seq data support widespread elevated chromatin accessibility possessing the CAATAATGC motif in Bif3.

Accessible chromatin regions with decreased accessibility due to overexpressed ZmWUS1 are associated with ARFs

We next examined ACRs with decreased central zone chromatin accessibility in Bif3 vs WT female inflorescence; these ACRs were associated with five different motifs: GCACAGCAGC, GCAGCATGC, CGCGCCGCGCC, GCTAGCTAGC, and AG repeats (Fig. 3A). The multiple motifs present in the decreased ACR set suggest that distinct TF families may be active in the Bif3 central zone, each recognizing their specific motif sequence. Near the decreased ACRs with these five motifs, we observed decreased gene body chromatin accessibility in Bif3 central zone genes (Additional file 1: Fig. S8). This indicates that differential ACRs decreased in Bif3 might be associated with reduced nearby gene transcription.

Gene ontology analysis of the central zone genes linked to ACRs decreased in Bif3 were associated with transcription regulation, developmental regulation, and auxin hormone responses (Fig. 4A). As meristem size regulation is associated with auxin [54], we investigated AUXIN RESPONSE FACTORS (ARFs) genes, the major regulators of the transcriptional auxin response [55]. We discovered 10 ZmARF genes nearby differential ACRs in the central zone, all of which had decreased accessibility in Bif3 compared to WT (Additional file 1: Fig. 4B; Additional file 1: Fig. S9A). ZmARFs linked to 10 of the reduced differential ACRs in the central zone, except ZmARF3, ZmARF18, and ZmARF36, showed reduced gene body chromatin accessibility in the central zone nuclei of Bif3 (Fig. 4B). ZmARF23, ZmARF20, and ZmARF4 had the three lowest logFC in gene body chromatin accessibility (− 0.75, − 0.66, and − 0.43, respectively; Fig. 4C). Thus, ACRs decreased in the Bif3 central zone are likely associated with decreased transcription of some ZmARFs, indicating that the negative regulation of ZmARFs is associated with increased meristem size.

Fig. 4
figure 4

TheARFgenes are associated with decreased differential ACRs by overexpression of ZmWUS1. A Bar graph illustrates significant biological process terms identified using a Fisher’s exact test using the genes closest to the differential ACRs in the central zone (FDR < 0.05). B The left heatmap shows the logFC of differential ACRs around ARF genes, whereas the right heatmap shows the gene body chromatin accessibility of ARFs by cell types. C Chromatin accessibility in the central zone around ZmARF4 and ZmARF23. The blue bar and yellow highlighted region indicate the differential ACRs

Full size image

We further categorized central zone ZmARFs near ACRs decreased in Bif3 into clade A and B [32, 33]. Among the 10 ZmARFs found near differential ACRs, ZmARF4, ZmARF18, ZmARF20, and ZmARF30 belong to clade A, whereas ZmARF10, ZmARF23, ZmARF25, and ZmARF36 are in clade B [56]. The differential ACRs near clade A ZmARFs lacked the five motifs enriched within the ACRs decreased in Bif3, whereas those ACRs near clade B ZmARFs contained at least one of these motifs, with all differential ACRs near clade B ZmARFs containing the CGCGCCGCGCC motif (Additional file 1: Fig. S9B). This CGCGCCGCGCC resembles the CG repeats found in AP2/ERF motifs from the Jaspar database (ERF109: MA1053.1; CRF2:MA0975.1; RAP2-3: MA1051.1). This indicates that the Bif3 negative regulation of ZmARFs involves both clade A and clade B, with clade B being predominantly regulated by a single CGCGCCGCGCC motif.

ZmWUS1 DAP-seq identified the TGAATGAA motif rather than the CAATAATGC motif

To supplement our in vivo data, we used DAP-seq [28] to identify DNA sequences from the maize genome directly bound in vitro by ZmWUS1 (Fig. 5A). In total, 10,382 peaks were identified, which uncovered TGAATGAA as the top-ranked motif consistently observed at the center of the peaks (Fig. 5B). The first occurrence of TGAA shows strong motif bit scores, whereas the second occurrence exhibits variability, especially in the tail of the motif (Fig. 5C). Analysis using Tomtom [57], the JASPAR motif database for plants [58], and the DAP-seq cistrome database [28] revealed that the ZmWUS1 DAP-seq TGAATGAA motif only shares similarity with the Arabidopsis WUS motif (E-value = 0.079), identified previously using ChIP-seq and DAP-seq [11, 28, 30]. This highlights the differences of the WUS TGAATGAA motif from that bound by other plant homeodomain TFs, corresponding to the unique role WUS plays in meristem maintenance.

Fig. 5
figure 5

The ZmWUS1 motif identified using DAP-seq shows different activity within ACRs depending on the cell type. A A genome browser view showing the peaks from ZmWUS1 DAP-seq. The y-axis shows CPM values. B A metaplot displaying the ratio of reads within a 1-kilobase pair region surrounding the ZmWUS1 DAP-seq peaks. C Identification of the ZmWUS1 motif using a k-mer set analysis, along with PWM models of the most enriched and significant motifs. D The motif deviation of ZmOCL1 in WT and Bif3. The motif is known for Arabidopsis HDG1, which is orthologous to ZmOCL1. The box with arrows indicates the epidermis cells. E The motif deviation of ZmWUS1 in WT and Bif3. The box with arrows indicates the central zone cells

Full size image

We hypothesized that the TGAATGAA motif would be accessible for binding in the central zone nuclei where ZmWUS1 is active. First, we identified all intergenic ACRs across all cell types in each genotype, leaving 58,424 and 49,618 ACRs in WT and Bif3, respectively. Next, we used chromVar [59] to calculate motif deviation scores, which represents the likelihood of observing a specific motif with high chromatin accessibility within the intergenic ACRs in each nucleus. We previously indicated that TFs with high gene body chromatin accessibility tend to have high cognate motif deviations [37]. For example, based on orthology to Arabidopsis HOMEODOMAIN GLABROUS1 (HDG1), the ZmOCL1 TF is predicted to bind a TAATTAATG motif in outer cell layers. Motif deviation for ZmOCL1 indicated that TAATTAATG is enriched in intergenic ACRs in epidermal cells for both WT and Bif3 (Fig. 5D). This suggests that the ZmOCL1 TF, expressed in the epidermis [60], binds to ACRs within the same cells where it is expressed. For the ZmWUS1 TGAATGAA motif identified via DAP-seq, we observed low motif deviation scores in the central zone of both WT and Bif3 (Fig. 5E), indicating that TGAATGAA motifs are less prevalent in the central zone ACRs compared to other cell types. In both genotypes, the highest motif deviation scores for the TGAATGAA motif were in the epidermis, proto/meta xylem, or procambial cells. This suggests that there is a spatial separation between where Bif3 alters chromatin accessibility, the central zone, and where the in vitro determined ZmWUS1 motif is most accessible.

Discussion

ZmWUS1 plays a crucial role in inflorescence meristem development as a TF that binds to cis-regulatory elements eventually facilitating kernel formation. ZmWUS1 expression is confined to a population of cells in the central zone of meristems [10]. The restricted domain of ZmWUS1 expression complicates identifying its associated cis-regulatory elements. To address this, we compared cell-type-resolved chromatin accessibility differences in the central zone between WT and Bif3, a mutant genotype with increased ZmWUS1 expression. While the chromatin accessibility changes observed in Bif3 may not fully reflect ZmWUS1’s natural binding behavior, this overexpression-based approach provides a valuable alternative for uncovering candidate binding motifs and inferring TF activity. We identified 1437 regions with increased chromatin accessibility and 1757 regions with decreased chromatin accessibility in Bif3. The ACRs with increased accessibility in the Bif3 central zone nuclei were enriched for the CAATAATGC motif, comprising approximately 26.51% of the ACRs with increased accessibility. Arabidopsis WUS can bind TAAT repeats in in vitro binding assays [17, 23], which is contained in the CAATAATGC motif. The extended flanking sequences in the CAATAATGC motif likely reflect ZmWUS1-specific binding preferences, potentially influenced by central zone-specific cofactors or dimerization.

The increases in accessibility of ACRs containing the CAATAATGC motif in the central zone cells of Bif3 are either driven by direct ZmWUS1 actions or other differentially expressed homeodomain TFs. To explore the latter, we investigated the TFs with homeodomain that have increased gene body chromatin accessibility in the Bif3 central zone compared to WT. We identified five homeodomain TFs whose gene body chromatin accessibility in Bif3 exceeds that of ZmWUS1, making them prime candidates to outcompete ZmWUS1 at shared binding sites. This criterion ensures we focus on TFs likely to exert a stronger regulatory influence than ZmWUS1 itself. We caution that ZmWUS1’s measured gene body chromatin accessibility is probably underestimated. ZmWUS1 expression was increased due to a tandem duplicated copy of the ZmWUS1 gene, ZmWUS1-B, along with a novel 119 bp region in the promoter region [11, 13]. The 119 bp insertion in the ZmWUS1-B disrupts read mapping to the v5 reference genome, likely reducing the gene body chromatin accessibility signal of ZmWUS1. Therefore, the most parsimonious explanation is that the “increased in Bif3” ACRs with CAATAATGC motifs are either driven by direct ZmWUS1 actions in the central zone over CAATAATGC motifs or by together with other homeodomain TFs. As WUS is proposed to be a dual function TF, both activating or repressing gene expression in Arabidopsis depending on context [20], we examined gene body chromatin accessibility changes near ACRs harboring the CAATAATGC motif and discovered that their gene body chromatin accessibility mostly increases. This suggests that ZmWUS1 could act as an activator in the maize central zone.

The TGAATGAA motifs identified by ZmWUS1 DAP-seq were not predominant in the ACRs of both Bif3 and WT central zone nuclei, where, instead, the CAATAATGC motif was prevalent. DAP-seq captures ZmWUS1 binding, either as a monomer or a homodimer [28], but without interacting partners. Therefore, TGAATGAA binding is recognized by either ZmWUS1 monomers or homodimers. Based on structural studies of the Arabidopsis WUS, the TGAATGAA binding likely reflects a homodimer interaction [11]. The enrichment of the TGAATGAA motif in the epidermis suggests that ZmWUS1 homodimers may be active in epidermis cells (Fig. 6B). However, the absence of central zone TGAATGAA motif accessibility change in Bif3 suggests additional factors influence ZmWUS1 in vivo central zone binding preferences. One possibility is that there could be other ZmWUS1-TF protein interactions, similar to the demonstrated HAM-WUS binding in Arabidopsis [20], that may alter DNA affinity. If these cofactors occur in the central zone and not in other cell types, like the epidermis, it may hypothetically explain the variance in motif usage—in the central zone, the cofactors outcompete ZmWUS1 homodimer formation, but movement of ZmWUS1 out of the central zone allows a monomer or homodimer formation, altering motif usage (Fig. 6A). Alternatively, cell-type-specific post-translational modifications of ZmWUS1 could likewise affect motif preference by altering binding partner or direct DNA-binding affinity [61, 62]. Future research into ZmWUS1 protein–protein interactions and post-translational modifications may reveal the molecular mechanisms that explain the discrepancy between the in vitro predicted ZmWUS1 motif and the observed in vivo differences in Bif3 central zone motif chromatin accessibility.

Fig. 6
figure 6

A hypothetical model that could explain central zone specific ZmWUS1 cis-regulation. A Proposed model of motif usage of ZmWUS1 in central zone cells. The CAATAATGC motif was identified in the central zone cells as a potential ZmWUS1 binding site, due to its presence in the differential ACRs of the central zone, where chromatin accessibility increased in Bif3. The cell-type-specific motif usage by ZmWUS1 could be attributed to unique molecular interactions in central zone cells. B A hypothetical epidermal model. In contrast, to the central zone, ACRs in the epidermis were largely unchanged in the Bif3 mutant, suggesting that the effect of overexpressed ZmWUS1 may be minimal or absent in epidermis cells. However, the ZmWUS1 DAP-seq motif, TGAATGAA, was enriched in the epidermal ACRs. It remains unclear how ZmWUS1 might engage this motif, although one possibility is that epidermal cells acquire ZmWUS1 activity through translocation from the central zone. A structural prediction of ZmWUS1 as a homodimer binding the TGAATGAA motif is shown

Full size image

The Arabidopsis WUS and ZmWUS1 TGAATGAA binding motif is unique amongst known plant homeodomain TF motifs; it seems unlikely any other characterized TF uses this motif outside of the central zone. Therefore, given the mobility of Arabidopsis WUS, it is tempting to speculate that ZmWUS1 moves into these cell types enriched for the TGAATGAA motif. However, why then would the Bif3 ZmWUS1 overexpression not impact the chromatin environments in these more distal cell types? One possibility is that the ZmWUS1 transport pathway becomes rate limiting, both restricting the amount of ZmWUS1 that reaches these distal cell types to a level analogous to WT levels and potentially constraining ZmWUS1 overaccumulation to the central zone, where we observe most Bif3 molecular changes. It would be interesting to see how the TGAATGAA motif accessibility changes after reduction or loss of ZmWUS1 expression/function. Furthermore, given that the canonical ZmWUS1 target, ZmCLE7, has overlapping expression patterns with ZmWUS1, it remains possible that, unlike in Arabidopsis, the ZmWUS1 protein is cell autonomous—further inquiry into the scope and potential molecular impacts of ZmWUS1 movement is needed to shed light on these questions.

In Arabidopsis, the central zone and organizing center, defined by localized WUS expression, are fairly distinct domains, and the models proposed for WUS function [20, 63, 64] indicate that WUS functions as a repressor in the organizing center and as an activator in the central zone. In maize, however, ZmWUS1 expression is broader in inflorescence meristems [11, 14] and there is no clear evidence that a domain separation between the central zone and organizing center exists. Indeed, in our data, ZmWUS1 and CLE7 gene body chromatin accessibility shows significant overlap. It is therefore likely that, unlike in Arabidopsis, a clear spatial separation of ZmWUS1 function as an activator or a repressor of transcription is not occurring in maize meristems. The decreased chromatin accessibility in the central zone nuclei due to overexpressed ZmWUS1 raises intriguing questions as the motifs (GCACAGCAGC, GCAGCATGC, CGCGCCGCGCC, GCTAGCTAGC, and AG repeats) enriched within these ACRs do not correspond to known ZmWUS1 binding sequences. ZmWUS1 could potentially interact and alter other TF complexes, producing the reduction in chromatin accessibility and local gene expression. Alternatively, the decreased chromatin accessibility may be from ZmWUS1 indirect effects. For example, AG repeats are one class of polycomb response elements [65, 66], suggesting potential alteration to Polycomb repressive complex 2 silencing in Bif3. Again, finding direct ZmWUS1 targets will help to resolve these possibilities.

The reduced chromatin accessibility near ZmARF genes links ZmWUS1 overexpression and diminished central zone auxin signaling. This observation aligns with the known interactions between Bif3 and auxin-insensitive mutants[11]. In Arabidopsis, the shoot apical meristem exhibits lower auxin signaling in the central zone compared to peripheral regions [67]. Previous research indicates that induced WUS expression further decreases auxin signaling in the central zone, suggesting WUS shields stem cells from auxin-driven differentiation [30]. Our results possibly suggest that this auxin-related regulatory mechanism involves ZmWUS1 negatively affecting auxin signaling via reducing chromatin accessibility of ARFs. In Bif3, this may retard stem cell differentiation, causing increased stem cell divisions, and thus count, in the inflorescence meristems.

Conclusions

Our study presents that overexpression of ZmWUS1 induces domain-specific changes in chromatin accessibility in the maize inflorescence meristem. The enrichment of the CAATAATGC motif suggests ZmWUS1 might recognize context-dependent binding, potentially shaped by central zone-specific cofactors or dimerization. Decreased chromatin accessibility near ZmARF genes suggests that ZmWUS1 overexpression interferes with auxin signaling in the central zone. Overall, these findings highlight the molecular mechanisms through which ZmWUS1 maintains the inflorescence meristem stem cell niche.

Methods

Plant material and growth conditions

Immature ears were collected from the inbred line A619 wild-type, and +/Bif3 plants (back-crossed 12 times in A619) grown in winter 2020 in the Waksman Institute greenhouse (Piscataway, NJ, USA), and in summer 2021 in the Waksman Institute field (Piscataway, NJ, USA). Immature ears were 2–6 mm in size and were collected at the same time for both genotypes. The Bif3 mutant phenotype is completely dominant in A619 ears [29]. Immature ears were harvested (36 and 31 samples for A619, and 41 and 36 for +/Bif3 in 2020 and 2021, respectively) approximately at the same time of the day (between 11 am and 2 pm) and bulked for subsequent analysis. Greenhouse-sown seeds were grown in 5-gallon pots filled with PROMIX BX General Purpose Soil and supplemented with 14–14-14 Osmocote. Soil was saturated with tap water and placed under GE Lucalox 400 W lighting. Seedlings were grown under a photoperiod of 16 h of light, 8 h of dark. The temperature was maintained at ~ 25 °C during light hours with a relative humidity of approximately 30%. Greenhouse growing conditions were monitored and controlled by MicroGrow Greenhouse Systems.

Single-cell ATAC-seq library construction

The protocol for nuclei isolation and purification was adapted from a previously described method with specific modifications to improve nuclei yield and quality [68]. In brief, approximately 5–6 immature maize ears were fine chopped on ice for approximately 1 min using 300 µL of pre-chilled Nuclei Isolation Buffer (NIB, 10 mM MES-KOH at pH 5.4, 10 mM NaCl, 250 mM sucrose, 0.1 mM spermine, 0.5 mM spermidine, 1 mM DTT, 1% BSA, and 0.5% TritonX-100). Following chopping, the mixture was filtered through a 40-µm cell strainer and then centrifuged at 500 rcf for 5 min at 4 °C. The supernatant was carefully removed, and the pellet was washed in 500 µL of NIB wash buffer, comprised of 10 mM MES-KOH at pH 5.4, 10 mM NaCl, 250 mM sucrose, 0.1 mM spermine, 0.5 mM spermidine, 1 mM DTT, and 1% BSA. Then, the sample was filtered again through a 20-µm cell strainer, and then gently layered onto the surface of 500 µL of a 35% Percoll buffer, prepared by mixing 35% Percoll with 65% NIB wash buffer, in a 1.5-mL centrifuge tube. The nuclei were centrifuged at 500 rcf for 8 min at 4 °C. After centrifugation, the supernatant was carefully removed, and the pellets were resuspended in 10 µL of diluted nuclei buffer (DNB, 10 × Genomics Cat# 2,000,207). Approximately 5 µL of nuclei were diluted 10 times, stained with DAPI (Sigma Cat. D9542), and then assessed for quality and density with a hemocytometer under a microscope. The original nuclei were diluted with DNB buffer to a final concentration of 3200 nuclei per µL. Finally, 5 µL of nuclei (16,000 nuclei in total) were used as input for scATAC-seq library preparation. scATAC-seq libraries were generated using the Chromium scATAC v1.1 (Next GEM) kit from 10 × Genomics (Cat# 1,000,175), following the manufacturer’s instructions (10xGenomics, CG000209_Chromium_NextGEM_SingleCell_ATAC_ReagentKits_v1.1_UserGuide_RevE). The libraries were sequenced with Illumina NovaSeq 6000 in dual-index mode with 8 and 16 cycles for i7 and i5 index, respectively.

DAP-seq library construction

DAP-seq libraries for assessing ZmWUS1 binding to maize genomic DNA were created as described in Galli et al. [56]. In brief, genomic DNA was extracted from the leaves of 14-day-old B73 seedlings using a phenol:cholorform:Iso amyl alcohol procedure. Five micrograms of extracted DNA was sheared in a Covaris S2 sonicator to 200 bp fragments. The fragmented DNA was then end repaired using the END-IT DNA repair kit (Lucigen) per the manufacturer’s recommendations and cleaned using a QiaQuick PCR clean-up method (Qiagen). Eluted DNA was then A-tailed using Klenow exo- (NEB) for 30 min at RT. Samples were cleaned using a QiaQuick PCR-clean-up method. Eluted DNA was ligated to truncated Illumina Y-adapters overnight at 16 °C using T4 DNA ligase (NEB). Samples were heat inactivated for 10 min at 70 °C followed by a 1:1 bead clean-up using AmpureXP beads (Beckman-Coulter). Libraries were quantified using the Qubit HS kit. One microgram of this library was then incubated with in vitro expressed HALO-WUS1 protein bound to 10 µl of MagneHALO beads (TNT rabbit reticulocyte expression system and MagneHALO beads, Promega) generated from 1 µg of pIX-HALO-WUS1 plasmid. Unbound DNA was washed away using six washes of phosphate buffered saline (PBS), and bound DNA was eluted from the beads using 30 µl of 10 mM Tris (EB) and heating to 98 °C for 10 min. DNA was barcoded and enriched with 19 cycles and Illumina primers prior to sequencing on a NextSeq500 with 75 bp Single End reads.

High-quality nuclei identification

The reads from scATAC-seq are mapped using the CellRanger-atac count (v2.0.0) by 10 × Genomics to the modified B73 v5 reference genome. In B73 v5, the scaffolds are removed, while Mt (mitochondria) and Pt (plastids) are added from the v3 genome. The multi-mapped reads with BWA added tags of “XA” and “SA” were removed [69]. De-duplication was performed by picard MarkDuplicates (v2.16.0) (https://broadinstitute.github.io/picard/). Single-base pair Tn5 integration was identified using a python script [70].

Socrates, R package (https://github.com/plantformatics/Socrates), was used to filter out nuclei that did not meet quality thresholds [37]. These criteria included the quantification of Tn5 transposase insertions, the proximity ratio of Tn5 insertions to the transcription start site (TSS), the fraction of reads in peaks (FRiP) score, and the organelle DNA ratio. We excluded nuclei falling below the knee point identified in the plot correlating unique Tn5 integration sites per nucleus with barcode rank. Nuclei were filtered if they exhibited TSS ratio or FRiP score beyond three standard deviations from the mean, and if the organelle DNA ratio exceeded 5%. Additionally, nuclei with fewer than 100 accessible chromatin regions, or with an open chromatin peak ratio below 0.01 or above 0.05, were excluded from subsequent analysis.

We used the presence or absence of Tn5 insertions within 500 bp windows as features for dimensionality reduction. A blacklist generated from comparative ATAC-seq to genomic DNA and control ChIP-seq data was applied to remove features. The cleanData function in the Socrates package was used to apply feature filters based on their frequency across the cell population. Features present in more than 5% of cells were retained using the “max.t = 0.05” parameter, while those observed in fewer than 1% of cells were excluded with the “min.t = 0.01” threshold. Feature normalization was performed using the term frequency-inverse document frequency transformation (TF-IDF) [71]. The top third of the features, identified as highly variable, were selected for further dimensionality reduction using the reduceDims function in Socrates.

Singular value decomposition (SVD) is used to reduce dimensionality and compute the principal components [72]. The top 100 principal components were used for Uniform Manifold Approximation Projection (UMAP). Doublets were identified and removed using Scrublet using the software Scrublet as implemented in detectDoublets and filterDoublets functions in Socrates with the option of filterRatio = 1.5 [73]. Finally, to integrate the two replicates, we applied the Harmony algorithm with “theta = 2, sigma = 0.1, max.iter.cluster = 100, and max.iter.harmony = 30” [74]. Cell clustering was performed using Louvain clustering on k = 50 nearest neighborhood graph with a resolution of 0.3. Subsequent cluster assignment by Louvain community detection identified distinct cell clusters.

Gene body chromatin accessibility analysis

Gene body chromatin accessibility is Tn5 insertion counts within gene bodies. Gene body chromatin accessibility served as a basis for cell annotation, leveraging the variability across cells. To capture cis-regulatory elements influencing gene expression, we extended the regions of analysis 500 bp upstream and downstream of gene bodies. Tn5 insertions within gene bodies and these extended regions were counted using the GenomicRanges package and findOverlaps function [75]. We used marker genes referenced in prior scATAC-seq literature [37] and additionally collected marker genes known for their localization in specific cell types, as identified in other literature through in situ hybridization. The cell cycle ratio was evaluated for each cell type using cell-cycle stage marker genes that include both established cell cycle regulators and newly identified in maize [76]. We added the ZmCLE7 gene from reference genome AGPv3.19 into the B73 v5 genome, by finding the sequence match at “chr4:8,531,149–8,531,929” for ZmCLE7 in v5 using BLAST. To ensure complete coverage of the manually added ZmCLE7 region in v5, we extended the region by 1000 bp upstream and downstream of the gene body. Cell annotations were performed using a UMAP plot, depicting gene body accessibility per nucleus. For enhanced clarity in visualization, we applied smoothing to the normalized gene accessibility scores using a diffusion nearest neighbor graph [77]. The other usage of gene body chromatin accessibility is to assess differences in gene body chromatin accessibility between genotypes. We calculated the log2 FC in chromatin accessibility after excluding genes with fewer than 50 Tn5 insertions to mitigate extreme values caused by sparse data. For log2 FC calculation, CPM normalization was performed followed by quantile normalization. When comparing gene body chromatin accessibility between different genes, Transcripts Per Million (TPM) values were used.

To identify de novo cell type marker genes, we used chromatin gene body accessibility for clusters. For each target cluster, we constructed a null model by randomly sampling an equal number of cells from all other clusters, separately for each biological replicate. These null cells were matched to the target cluster cells by total Tn5 insertion counts. Differential analysis for each gene was performed using DESeq2 [78], applying a p-value cutoff of 0.01. The top five de novo marker genes were selected based on the lowest p-values.

Identification of accessible chromatin regions

Accessible chromatin regions were identified by aggregating cells with the same annotation into pseudo bulks. Peak calling was executed on pooled and individual replicates. MACS2 (v 2.2.7.1) was utilized with “–nomoel –keep-dup auto –extsize 150 –shift 75 –qvalue 0.05” [79]. The summits identified in each replicate were expanded by 250 bp on both sides. The peaks with less than 20 Tn5 coverage were filtered. Only the peaks present in all replicates and that overlapped a pooled set peak were retained [70].

Correlation analysis for the ACR by cell types

To compare chromatin accessibility within peaks across cell types, we used peaks called individually for each genotype and corresponding cell type. For the correlation analysis, we established a set of 500 bp union peaks. When using a fixed size of 500 bp peaks, we assessed the p-value of each peak across cell types, utilizing a chromatin accessibility score that was normalized per million reads to identify the most representative peaks [80]. The aggregated Tn5 insertions for each cell type were quantile normalized, and Pearson correlation coefficients were calculated using read density for the 2000 most variable peaks.

Differential ACR analysis

For differential ACR analysis, we created the sets of union peaks per cell type. We merged peaks from both genotypes for each cell type, resulting in distinct sets of union peaks specific to each genotype. Consequently, the lengths of the peaks used for differential ACR analysis vary. Two pseudo replicates were generated for each biological replicate within cell types to enhance the robustness by randomly partitioning cells in one replicate into two groups. To quantify chromatin accessibility in the pseudo bulk replicates, we aggregated the chromatin accessibility within peaks by pseudo replicates to each cell type. We normalized aggregated chromatin accessibility as Trimmed Mean of M values (TMM) to adjust for different library sizes and performed statistical tests under generalized linear model using edgeR (v3.32.1) [81].

Motif analysis in accessible chromatin regions

We used three approaches to find the motif enrichment in ACRs: (1) Motif deviation calculation across cells, (2) de novo motif searches, and (3) known motif searches. To identify the cells where the motif is active, we computed the motif deviation score using chromVar [59]. We used Tn5 insertion counts at intergenic ACRs per cell as input for chromVar. For motif PWM, we used the non-redundant core plant PWM database from JASPAR2022 and PWMs derived from DAP-seq motif discovery. We applied smoothing to the bias-corrected motif deviations for each nucleus, integrating them into UMAP embedding for visualization, which is the same method used for visualizing gene body chromatin accessibility [77]. Additionally, we rescaled the bias-corrected motif deviations to fit a color scale ranging from −1 to 1 across all motifs.

De novo motif searches in differential ACRs were performed using XSTREME version 5.5.3 within the MEME suite package (v5.5.0) [82, 83]. XSTREME leverages STREME, which finds the enriched motif in the test set relative to control sets [84]. While the test set is differential ACRs, control regions were used to determine the significance of motif occurrence compared to the background. To create the control set, we randomly selected the same number and length of ACRs from all ACRs, ensuring that they had a similar GC content ratio to the test set. To account for varied ACR lengths in the test set, we ensured the control set featured a comparable distribution of lengths. The motif search with STREME halts when it encounters a succession of motifs with p-values exceeding 0.05, applying “–thresh” option. We showcased motifs boasting an e-value below 1. Furthermore, utilizing the PWM of identified de novo motifs, FIMO scrutinized specific ACR regions [85], pinpointing motif occurrences with p-values under 0.0001.

DAP-seq analysis and motif discovery

Raw reads were trimmed using trimmomatic [86] with the following settings: ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:2:keepBothReads LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50. Trimmed reads were mapped to the B73v5 genome using bowtie2 [87] with default settings. Mapped reads were filtered to retain those with a MAPQ score greater than 30 using samtools [88] view -b -q 30. To identify peaks and motifs, we employed GEM (v 3.0) [89] with the following options: “–q 5 –k_min 5 –k_max 14”. This option uses significant peaks with a q-value below 5 and specifies a range of k-mer lengths from 5 to 14 for motif discovery.

Data availability

Single-cell ATAC-seq and DAP-seq have been deposited at the NCBI GEO database under accession code GSE266007 [90]. Metadata by barcode for WT and *Bif3* are provided in Additional file 2: Table S2 and Additional file 2: Table S3, respectively. All scripts used in this study are available on GitHub under an MIT open-source license [91, 92]. A version of the source code used in this manuscript has been archived in Zenodo [93].

References

  1. Laux T, Mayer KF, Berger J, Jürgens G. The wuschel gene is required for shoot and floral meristem integrity in Arabidopsis. Development. 1996;122:87–96.

    Article  CAS  PubMed  Google Scholar 

  2. Clark SE, Running MP, Meyerowitz EM. CLAVATA3 is a specific regulator of shoot and floral meristem development affecting the same processes as CLAVATA1. Development. 1995;121:2057–67.

    Article  CAS  Google Scholar 

  3. Schoof H, Lenhard M, Haecker A, Mayer KF, Jürgens G, Laux T. The stem cell population of Arabidopsis shoot meristems is maintained by a regulatory loop between the CLAVATA and WUSCHEL genes. Cell. 2000;100:635–44.

    Article  CAS  PubMed  Google Scholar 

  4. Gaillochet C, Lohmann JU. The never-ending story: from pluripotency to plant developmental plasticity. Development. 2015;142:2237–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Suzaki T, Toriba T, Fujimoto M, Tsutsumi N, Kitano H, Hirano H-Y. Conservation and diversification of meristem maintenance mechanism in Oryza sativa: function of the FLORAL ORGAN NUMBER2 gene. Plant Cell Physiol. 2006;47:1591–602.

    Article  CAS  PubMed  Google Scholar 

  6. Xu C, Liberatore KL, MacAlister CA, Huang Z, Chu Y-H, Jiang K, et al. A cascade of arabinosyltransferases controls shoot meristem size in tomato. Nat Genet. 2015;47:784–92.

    Article  CAS  PubMed  Google Scholar 

  7. Rodríguez-Leal D, Lemmon ZH, Man J, Bartlett ME, Lippman ZB. Engineering quantitative trait variation for crop improvement by genome editing. Cell. 2017;171(470–480):e478.

    Google Scholar 

  8. Suzuki C, Tanaka W, Tsuji H, Hirano H-Y. Tillers ABSENT1, the WUSCHEL ortholog, is not involved in stem cell maintenance in the shoot apical meristem in rice. Plant Signal Behav. 2019;14:1640565.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Tanaka W, Ohmori Y, Ushijima T, Matsusaka H, Matsushita T, Kumamaru T, et al. Axillary meristem formation in rice requires the WUSCHEL ortholog TILLERS ABSENT1. Plant Cell. 2015;27:1173–84.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Nardmann J, Werr W. The shoot stem cell niche in angiosperms: expression patterns of WUS orthologues in rice and maize imply major modifications in the course of mono-and dicot evolution. Mol Biol Evol. 2006;23:2492–504.

    Article  CAS  PubMed  Google Scholar 

  11. Chen Z, Li W, Gaines C, Buck A, Galli M, Gallavotti AJ. Structural variation at the maize WUSCHEL1 locus alters stem cell organization in inflorescences. Nat Commun. 2021;12:1–12.

    Google Scholar 

  12. Chen Z, Gallavotti A. Improving architectural traits of maize inflorescences. Mol Breed. 2021;41:1–13.

    Article  CAS  Google Scholar 

  13. Chen Z, Cortes L, Gallavotti A. Genetic dissection of cis-regulatory control of ZmWUSCHEL1 expression by type B response regulators. Plant Physiol. 2024;194:2240–8.

    Article  CAS  PubMed  Google Scholar 

  14. Je BI, Gruel J, Lee YK, Bommert P, Arevalo ED, Eveland AL, et al. Signaling from maize organ primordia via FASCIATED EAR3 regulates stem cell proliferation and yield traits. Nat Genet. 2016;48:785–91.

    Article  CAS  PubMed  Google Scholar 

  15. Liu L, Gallagher J, Arevalo ED, Chen R, Skopelitis T, Wu Q, et al. Enhancing grain-yield-related traits by CRISPR–Cas9 promoter editing of maize CLE genes. Nat Plants. 2021;7:287–94.

    Article  CAS  PubMed  Google Scholar 

  16. Fletcher JC, Brand U, Running MP, Simon R, Meyerowitz EM. Signaling of cell fate decisions by CLAVATA3 in Arabidopsis shoot meristems. Science. 1999;283:1911–4.

    Article  CAS  PubMed  Google Scholar 

  17. Mayer KF, Schoof H, Haecker A, Lenhard M, Jürgens G, Laux T. Role of WUSCHEL in regulating stem cell fate in the Arabidopsis shoot meristem. Cell. 1998;95:805–15.

    Article  CAS  PubMed  Google Scholar 

  18. Yadav RK, Perales M, Gruel J, Girke T, Jönsson H, Reddy GV. Wuschel protein movement mediates stem cell homeostasis in the Arabidopsis shoot apex. Genes Dev. 2011;25:2025–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Daum G, Medzihradszky A, Suzaki T, Lohmann JU. A mechanistic framework for noncell autonomous stem cell induction in Arabidopsis. Proc Natl Acad Sci U S A. 2014;111:14619–24.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Zhou Y, Yan A, Han H, Li T, Geng Y, Liu X, et al. Hairy meristem with WUSCHEL confines CLAVATA3 expression to the outer apical meristem layers. Science. 2018;361:502–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Perales M, Rodriguez K, Snipes S, Yadav RK, Diaz-Mendoza M, Reddy GV. Threshold-dependent transcriptional discrimination underlies stem cell homeostasis. Proc Natl Acad Sci U S A. 2016;113:E6298-306.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Noyes MB, Christensen RG, Wakabayashi A, Stormo GD, Brodsky MH, Wolfe SA. Analysis of homeodomain specificities allows the family-wide prediction of preferred recognition sites. Cell. 2008;133:1277–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Su YH, Zhou C, Li YJ, Yu Y, Tang LP, Zhang WJ, et al. Integration of pluripotency pathways regulates stem cell maintenance in the Arabidopsis shoot meristem. Proc Natl Acad Sci U S A. 2020;117:22561–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Leibfried A, To JP, Busch W, Stehling S, Kehle A, Demar M, et al. WUSCHEL controls meristem function by direct regulation of cytokinin-inducible response regulators. Nature. 2005;438:1172–5.

    Article  CAS  PubMed  Google Scholar 

  25. Liu X, Kim YJ, Müller R, Yumul RE, Liu C, Pan Y, et al. AGAMOUS terminates floral stem cell maintenance in Arabidopsis by directly repressing WUSCHEL through recruitment of Polycomb Group proteins. Plant Cell. 2011;23:3654–70.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Huang H, Mizukami Y, Hu Y, Ma H. Isolation and characterization of the binding sequences for the product of the Arabidopsis floral homeotic gene AGAMOUS. Nucleic Acids Res. 1993;21:4769–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Lohmann JU, Hong RL, Hobe M, Busch MA, Parcy F, Simon R, et al. A molecular link between stem cell regulation and floral patterning in Arabidopsis. Cell. 2001;105:793–803.

    Article  CAS  PubMed  Google Scholar 

  28. O’Malley RC, Huang S-sC, Song L, Lewsey MG, Bartlett A, Nery JR, et al. Cistrome and epicistrome features shape the regulatory DNA landscape. Cell. 2016;165:1280–92.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Sloan J, Hakenjos JP, Gebert M, Ermakova O, Gumiero A, Stier G, et al. Structural basis for the complex DNA binding behavior of the plant stem cell regulator WUSCHEL. Nat Commun. 2020;11:2223.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Ma Y, Miotk A, Šutiković Z, Ermakova O, Wenzl C, Medzihradszky A, et al. WUSCHEL acts as an auxin response rheostat to maintain apical stem cells in Arabidopsis. Nat Commun. 2019;10:5093.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Vernoux T, Besnard F, Traas J. Auxin at the shoot apical meristem. Cold Spring Harb Perspect Biol. 2010;2:a001487.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Finet C, Berne-Dedieu A, Scutt CP, Marlétaz F. Evolution of the ARF gene family in land plants: old domains, new tricks. Mol Biol Evol. 2013;30:45–56.

    Article  CAS  PubMed  Google Scholar 

  33. Ulmasov T, Hagen G, Guilfoyle TJ. Activation and repression of transcription by auxin-response factors. Proc Natl Acad Sci U S A. 1999;96:5844–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Hernández-García J, Carrillo-Carrasco VP, Rienstra J, Tanaka K, de Roij M, Dipp-Álvarez M, et al. Evolutionary origins and functional diversification of auxin response factors. Nat Commun. 2024;15:10909.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Vernoux T, Brunoud G, Farcot E, Morin V, den Van Daele H, Legrand J, et al. The auxin signalling network translates dynamic input into robust patterning at the shoot apex. Mol Syst Biol. 2011;7:508.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Liscum E, Reed J. Genetics of aux/IAA and ARF action in plant growth and development. Plant Mol Biol. 2002;49:387–400.

    Article  CAS  PubMed  Google Scholar 

  37. Marand AP, Chen Z, Gallavotti A, Schmitz RJ. A cis-regulatory atlas in maize at single-cell resolution. Cell. 2021;184(3041–3055):e3021.

    Google Scholar 

  38. Fuchs M, Lohmann JU. Aiming for the top: non-cell autonomous control of shoot stem cells in Arabidopsis. J Plant Res. 2020;133:297–309.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Vernoud V, Laigle G, Rozier F, Meeley RB, Perez P, Rogowsky PM. The HD-ZIP IV transcription factor OCL4 is necessary for trichome patterning and anther development in maize. Plant J. 2009;59:883–94.

    Article  CAS  PubMed  Google Scholar 

  40. Strable J, Vollbrecht E. Maize YABBY genes drooping leaf1 and drooping leaf2 regulate floret development and floral meristem determinacy. Development. 2019;146:dev171181.

    Article  CAS  PubMed  Google Scholar 

  41. Bortiri E, Chuck G, Vollbrecht E, Rocheford T, Martienssen R, Hake S. Ramosa2 encodes a lateral organ boundary domain protein that determines the fate of stem cells in branch meristems of maize. Plant Cell. 2006;18:574–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. van der Graaff E, Laux T, Rensing SA. The WUS homeobox-containing (WOX) protein family. Genome Biol. 2009;10:1–9.

    Google Scholar 

  43. Kawai T, Shibata K, Akahoshi R, Nishiuchi S, Takahashi H, Nakazono M, et al. Wuschel-related homeobox family genes in rice control lateral root primordium size. Proc Natl Acad Sci U S A. 2022;119:e2101846119.

    Article  CAS  PubMed  Google Scholar 

  44. Nardmann J, Zimmermann R, Durantini D, Kranz E, Werr W. Wox gene phylogeny in Poaceae: a comparative approach addressing leaf and embryo development. Mol Biol Evol. 2007;24:2474–84.

    Article  CAS  PubMed  Google Scholar 

  45. Jin J, Tian F, Yang D-C, Meng Y-Q, Kong L, Luo J, Gao G. PlantTFDB 4.0: toward a central hub for transcription factors and regulatory interactions in plants. Nucleic acids research 2016;45(D1):D1040–5.

  46. Ariel FD, Manavella PA, Dezar CA, Chan RL. The true story of the HD-Zip family. Trends Plant Sci. 2007;12:419–26.

    Article  CAS  PubMed  Google Scholar 

  47. Hamant O, Pautot V. Plant development: a TALE story. Comptes Rendus Biologies. 2010;333:371–81.

    Article  CAS  PubMed  Google Scholar 

  48. Windhövel A, Hein I, Dabrowa R, Stockhaus J. Characterization of a novel class of plant homeodomain proteins that bind to the C4 phosphoenolpyruvate carboxylase gene of Flaveria trinervia. Plant Mol Biol. 2001;45:201–14.

    Article  PubMed  Google Scholar 

  49. Halbach T, Scheer N, Werr W. Transcriptional activation by the PHD finger is inhibited through an adjacent leucine zipper that binds 14-3-3 proteins. Nucleic Acids Res. 2000;28:3542–50.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Portwood JL, Woodhouse MR, Cannon EK, Gardiner JM, Harper LC, Schaeffer ML, et al. MaizeGDB 2018: the maize multi-genome genetics and genomics database. Nucleic Acids Res. 2019;47:D1146–54.

    Article  PubMed  Google Scholar 

  51. Bharathan G, Janssen B-J, Kellogg EA, Sinha N. Phylogenetic relationships and evolution of the KNOTTED class of plant homeodomain proteins. Mol Biol Evol. 1999;16:553–63.

    Article  CAS  PubMed  Google Scholar 

  52. Kerstetter R, Vollbrecht E, Lowe B, Veit B, Yamaguchi J, Hake S. Sequence analysis and expression patterns divide the maize knotted1-like homeobox genes into two classes. Plant Cell. 1994;6:1877–87.

    CAS  PubMed  PubMed Central  Google Scholar 

  53. Vollbrecht E, Veit B, Sinha N, Hake S. The developmental gene Knotted-1 is a member of a maize homeobox gene family. Nature. 1991;350:241–3.

    Article  CAS  PubMed  Google Scholar 

  54. Kitagawa M, Jackson D. Control of meristem size. Annu Rev Plant Biol. 2019;70:269–91.

    Article  CAS  PubMed  Google Scholar 

  55. Guilfoyle TJ, Ulmasov T, Hagen G. The ARF family of transcription factors and their role in plant hormone-responsive transcription. Cell Mol Life Sci. 1998;54:619–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Galli M, Khakhar A, Lu Z, Chen Z, Sen S, Joshi T, et al. The DNA binding landscape of the maize AUXIN RESPONSE FACTOR family. Nat Commun. 2018;9:4526.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Gupta S, Stamatoyannopoulos JA, Bailey TL, Noble WS. Quantifying similarity between motifs. Genome Biol. 2007;8:1–9.

    Article  Google Scholar 

  58. Sandelin A, Alkema W, Engström P, Wasserman WW, Lenhard B. JASPAR: an open-access database for eukaryotic transcription factor binding profiles. Nucleic Acids Res. 2004;32:D91-4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Schep AN, Wu B, Buenrostro JD, Greenleaf WJ. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nat Methods. 2017;14:975–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Ingram GC, Magnard J-L, Vergne P, Dumas C, Rogowsky PM. ZmOCL1, an HDGL2 family homeobox gene, is expressed in the outer cell layer throughout maize development. Plant Mol Biol. 1999;40:343–54.

    Article  CAS  PubMed  Google Scholar 

  61. Barlev NA, Liu L, Chehab NH, Mansfield K, Harris KG, Halazonetis TD, et al. Acetylation of p53 activates transcription through recruitment of coactivators/histone acetyltransferases. Mol Cell. 2001;8:1243–54.

    Article  CAS  PubMed  Google Scholar 

  62. Yuan Z-l, Guan Y-j, Chatterjee D, Chin YE. Stat3 dimerization regulated by reversible acetylation of a single lysine residue. Science. 2005;307:269–73.

    Article  CAS  PubMed  Google Scholar 

  63. Zhou Y, Liu X, Engstrom EM, Nimchuk ZL, Pruneda-Paz JL, Tarr PT, et al. Control of plant stem cell function by conserved interacting transcriptional regulators. Nature. 2015;517:377–80.

    Article  CAS  PubMed  Google Scholar 

  64. Rodriguez K, Do A, Senay-Aras B, Perales M, Alber M, Chen W, et al. Concentration-dependent transcriptional switching through a collective action of cis-elements. Sci Adv. 2022;8:eabo6157.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  65. Xiao J, Jin R, Yu X, Shen M, Wagner JD, Pai A, et al. Cis and trans determinants of epigenetic silencing by Polycomb repressive complex 2 in Arabidopsis. Nat Genet. 2017;49:1546–52.

    Article  CAS  PubMed  Google Scholar 

  66. Yan H, Mendieta JP, Zhang X, Marand AP, Liang Y, Luo Z, Minow MA, Roulé T, Wagner D, Tu X: Evolution of plant cell-type-specific cis-regulatory elements. bioRxiv 2024:2024.2001. 2008.574753.

  67. Shi B, Guo X, Wang Y, Xiong Y, Wang J. Hayashi K-i, Lei J, Zhang L, Jiao Y: Feedback from lateral organs controls shoot apical meristem growth by modulating auxin transport. Dev Cell. 2018;44:204–216.e206.

  68. Zhang X, Marand AP, Yan H, et al. scifi-ATAC-seq: massive-scale single-cell chromatin accessibility sequencing using combinatorial fluidic indexing. Genome Biol. 2024;25:90. https://doi.org/10.1186/s13059-024-03235-5.

  69. Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10:giab008.

  70. Mendieta JP, Tu X, Yan H, Zhang X, Marrand AP, Zhong S, et al. Investigating the cis-regulatory basis of C3 and C4 photosynthesis in grasses at single-cell resolution. bioRxiv. 2024:2024.2001. 2005.574340.

  71. Cusanovich DA, Hill AJ, Aghamirzaie D, Daza RM, Pliner HA, Berletch JB, et al. A single-cell atlas of in vivo mammalian chromatin accessibility. Cell. 2018;174(1309–1324):e1318.

    Google Scholar 

  72. Witten DM, Tibshirani R, Hastie T. A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis. Biostatistics. 2009;10:515–34.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Wolock SL, Lopez R, Klein AM. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 2019;8(281–291):e289.

    Google Scholar 

  74. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with harmony. Nat Methods. 2019;16:1289–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, et al. Software for computing and annotating genomic ranges. PLoS Comput Biol. 2013;9:e1003118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Nelms B, Walbot V. Defining the developmental program leading to meiosis in maize. Science. 2019;364:52–6.

    Article  CAS  PubMed  Google Scholar 

  77. Van Dijk D, Sharma R, Nainys J, Yim K, Kathail P, Carr AJ, et al. Recovering gene interactions from single-cell data using data diffusion. Cell. 2018;174(716–729):e727.

    Google Scholar 

  78. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:1–21.

    Article  Google Scholar 

  79. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008;9:1–9.

    Article  Google Scholar 

  80. Corces MR, Granja JM, Shams S, Louie BH, Seoane JA, Zhou W, et al. The chromatin accessibility landscape of primary human cancers. Science. 2018;362:eaav1898.

    Article  PubMed  PubMed Central  Google Scholar 

  81. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    Article  CAS  PubMed  Google Scholar 

  82. Bailey TL, Johnson J, Grant CE, Noble WS. The MEME suite. Nucleic Acids Res. 2015;43:W39–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Grant CE, Bailey TL: XSTREME: Comprehensive motif analysis of biological sequence datasets. BioRxiv 2021:2021.2009. 2002.458722.

  84. Bailey TL. Streme: accurate and versatile sequence motif discovery. Bioinformatics. 2021;37:2834–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  85. Grant CE, Bailey TL, Noble WS. FIMO: scanning for occurrences of a given motif. Bioinformatics. 2011;27:1017–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  86. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

    Article  PubMed  PubMed Central  Google Scholar 

  89. Guo Y, Mahony S, Gifford DK. High resolution genome wide binding event finding and motif discovery reveals transcription factor spatial binding constraints. PLoS Comput Biol. 2012;8(8):e1002638.

  90. Schmitz, RJ. WUSCHEL-dependent chromatin regulation in maize inflorescence development at single-cell resolution. Datasets. Gene Expression Omnibus. 2025. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE266007.

  91. Bang, S. WUSCHEL maize scATAC-seq. GitHub. 2025. https://github.com/schmitzlab/WUSCHEL_maize_scATAC-seq.

  92. Bang, S. Single-cell genomics WT vs mutant comparison. Github. 2025 https://github.com/sohyunbk/single_cell_genomics_wt_mutant_comparison.

  93. Bang, S. WUSCHEL maize scATAC-seq. Zenodo. 2025. https://doi.org/10.5281/zenodo.16929585.

Download references

Acknowledgements

We are grateful to John Pablo Mendieta and Alexandre P. Marand for their valuable comments and for providing code used in prior analyses.

Peer review information

Wenjing She was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team. The peer-review history is available in the online version of this article.

Funding

This research was supported by the National Science Foundation (IOS-2026554) to RJS and AG as well as the UGA Office of Research to RJS.

Author information

Authors and Affiliations

  1. Institute of Bioinformatics, University of Georgia, Athens, GA, 30602, USA

    Sohyun Bang

  2. Department of Genetics, University of Georgia, Athens, GA, 30602, USA

    Xuan Zhang, Ziliang Luo, Mark A. A. Minow & Robert J. Schmitz

  3. Waksman Institute of Microbiology, Rutgers University, Piscataway, NJ, 08540, USA

    Jason Gregory, Zongliang Chen, Mary Galli & Andrea Gallavotti

  4. Department of Plant Biology, Rutgers University, New Brunswick, NJ, 08901, USA

    Jason Gregory & Andrea Gallavotti

Authors
  1. Sohyun Bang
    View author publications

    Search author on:PubMed Google Scholar

  2. Xuan Zhang
    View author publications

    Search author on:PubMed Google Scholar

  3. Jason Gregory
    View author publications

    Search author on:PubMed Google Scholar

  4. Ziliang Luo
    View author publications

    Search author on:PubMed Google Scholar

  5. Zongliang Chen
    View author publications

    Search author on:PubMed Google Scholar

  6. Mark A. A. Minow
    View author publications

    Search author on:PubMed Google Scholar

  7. Mary Galli
    View author publications

    Search author on:PubMed Google Scholar

  8. Andrea Gallavotti
    View author publications

    Search author on:PubMed Google Scholar

  9. Robert J. Schmitz
    View author publications

    Search author on:PubMed Google Scholar

Contributions

S.B., R.J.S., and A.G. designed the experiments and wrote the paper; S.B. and M.G. analyzed the data; M.M. contributed to paper preparation; X.Z., Z.L., J.G., Z.C., M.G. and A.G. generated materials and/or conducted the experiments. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Andrea Gallavotti or Robert J. Schmitz.

Ethics declarations

Ethics approval and consent to participate

Wenjing She was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team. The peer-review history is available in the online version of this article.

Consent for publication

Not applicable.

Competing interests

R.J.S. is a co-founder of REquest Genomics, LLC, a company that provides epigenomic services. The remaining authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Bang, S., Zhang, X., Gregory, J. et al. WUSCHEL-dependent chromatin regulation in maize inflorescence development at single-cell resolution. Genome Biol 26, 329 (2025). https://doi.org/10.1186/s13059-025-03791-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13059-025-03791-4

Keywords

Genome Biology

ISSN: 1474-760X

Contact us