Regulatory modules controlling maize inflorescence architecture

Molecular signatures of auxin response are detected prior to morphological changes in ramosa mutants

We first characterized morphological features of inflorescence primordia associated with developmental transitions in determinate and indeterminate axillary meristems. We used loss-of-function ramosa (ra) mutant ears to monitor these transitions upon genetic perturbation of the SPM determinacy pathway. Scanning electron microscopy (SEM) showed that branching phenotypes of ra mutant ears are first evident when the ears are ∼2 mm in length. At a length of 1 mm, wild-type and ra mutant ears have initiated SPMs, and are indistinguishable (Fig. 1B,G; Supplemental Fig. S1). By the time they are 2 mm, wild-type ears have initiated SMs, whereas ra1 mutant ears have initiated elongating branches, and lack SMs (Fig. 1C,E,H,I). In contrast, ra2 and ra3 mutants have initiated both SMs and elongating branches by 2 mm (Supplemental Fig. S1).

The phytohormone auxin induces organ formation, and expression of DR5, an auxin response reporter, marks the conversion from SPM to SM identity in maize (Gallavotti et al. 2008). To define the developmental window when SPM determinacy is established and SM identity is initiated, we analyzed DR5-ER∷RFP expression in developing wild-type and ra mutant ears. Upon formation of the SPM, low levels of DR5 mark its central domain in all genotypes (Supplemental Fig. S1). Subsequently, a strong DR5 maximum develops at the lower flank of the SPM, followed by a second maximum at the opposite flank (Fig. 1D,I). One millimeter ears were enriched for SPMs with one to two DR5 maxima (Supplemental Fig. S1). Interestingly, DR5 signal also was detected linking the two maxima in 1-mm ears of ra mutants (Fig. 1I; Supplemental Fig. S1), but not in wild-type (Fig. 1D), indicating that molecular changes are established prior to visible morphological phenotypes. In wild-type tassels, DR5 signal linking the two maxima was evident only in indeterminate BMs (Fig. 1N), but never in determinate SPMs. It is possible that maintaining a low level of response to auxin in the meristem, i.e., linking the two flanking maxima, could inhibit transition to determinate growth. These observations provide reference points for profiling developmental transitions in normal and branching mutant backgrounds.

Genetic perturbation of the RAMOSA branching network reveals converging developmental pathways

We used RNA-seq to profile transcriptional changes upon genetic disruption of the RA network, and compared differentially expressed (DE) genes in ra1, ra2, and ra3 mutant ears, relative to wild-type siblings, and across a 1- to 2-mm developmental transition (Fig. 2A; Supplemental Fig. S1; Supplemental Tables S1, S2). We observed a dynamic shift in the number of DE genes relative to wild-type shared between ra1 and ra2 mutants from 1 to 2 mm; by 2 mm, 89% of DE genes in ra2 were also DE in ra1, whereas at 1 mm, only 13% were DE in ra1 (Fig. 2A; Supplemental Table S2). A similar trend was observed for DE putative TF genes (Fig. 2B). This shift is concomitant with establishment of the branching phenotype and suggests that RA1- and RA2-dependent pathways converge toward a common molecular phenotype by 2 mm.

A shift in functional enrichment among DE genes shared between ra1 and ra2 mutants was also observed from 1 to 2 mm (Fig. 2C), but not for DE genes shared among all three mutants. For the latter, enriched processes among shared DE genes included macromolecular complex subunit organization (P = 4.07 × 10⁻¹⁰), nucleosome assembly (P = 5.12 × 10⁻⁸), and chromatin (P = 1.21 × 10⁻⁷); this functional profile was maintained in 1- and 2-mm ears (Fig. 2C) and may reflect the early, or primary, molecular responses during the developmental transition. At 2 mm, the largely shared molecular profile of ra1 and ra2 likely includes secondary responses downstream from the decision to branch, while the smaller number of DE genes in ra3 at this stage is consistent with a shift back to determinate fate after only a few indeterminate BMs are laid down at the base in this mutant. Genes DE in ra3 also showed unique functional enrichment compared with the other mutants, possibly due to its predicted role in sugar metabolism and signaling (Fig. 2C; Satoh-Nagasawa et al. 2006).

To identify early signatures of transcriptional regulation in axillary meristem fate, we analyzed expression patterns of DE TFs in 1-mm mutants. In general, DE TFs were expressed at lower levels in ra1 mutants relative to wild type, whereas in ra3 mutants, they tended to be up-regulated. Interestingly, those TFs DE in all mutants displayed common trends, either up- or down-regulated in all three (Fig. 2D). This suggests that certain TFs are independently regulated while others show common regulation, possibly functioning in shared pathways among ra genes. Notably, down-regulation of plant-specific TCP genes was observed in all mutants (Fig. 2D,E), consistent with members of this family regulating cell proliferation and branching (Cubas et al. 1999). In contrast, AP2-EREBP, bHLH, and C3H TF families tended to be up-regulated (Fig. 2D), as were genes encoding chromatin- and RNAi-associated factors (Supplemental Fig. S2; Supplemental Table S3). In some cases, differential expression profiles were observed for TFs within a family. For example, 12 members of a MADS-box family, implicated in developmental patterning, could be separated into up- and down-regulated classes across ra mutants (Fig. 2E). Interestingly, some MADS-box genes up-regulated in ra mutants were also up-regulated in tassel primordia compared with ears, suggesting they may underlie differences between ear and tassel architectures.

Co-expression signatures across spatiotemporal and genetic contexts predict a SPM determinacy module

To characterize gene expression signatures associated with developmental transitions in normal inflorescences, we clustered genes based on their transcript profiles across eight wild-type ear and tassel stages. In addition to 1- and 2-mm wild-type ears (above), we performed RNA-seq on samples sectioned from tip, middle, and base along a 10-mm ear, which were enriched for specific meristem types: IM/SPM, SM, and FM, respectively, and tassels at three stages of early development (stage [stg]1–3) (Supplemental Fig. S3; Supplemental Tables S1, S2, S4, S5). Stg1 tassels are enriched for indeterminate BMs, with determinate SPMs and SMs formed at stg2, and FMs formed by stg3. Therefore, profiling tassel development provided additional perspectives on indeterminate vs. determinate fate. From a set of 16,272 genes with dynamic expression profiles, we identified 20 co-expression clusters representing unique developmental signatures using a k-means clustering approach (Fig. 3A; Supplemental Table S6).

Co-expression clusters 8 and 11 were highly enriched for genes DE in 1-mm ra mutant ears, suggesting these signatures are associated with SPM determinacy (Fig. 3A,B,D). Both clusters showed enriched expression in stg2 tassels (Fig. 3B,D), and this is notable since the transition from stg1 to stg2 marks the time when production of BMs shifts to determinate SPMs (Fig. 1L,M). Expression of genes in cluster 8 was also strongly enriched in 1-mm wild-type ears, and those genes DE in all three mutants were coordinately repressed at 1 mm (Fig. 3E). Genes in cluster 11 and DE in all 1-mm mutant ears were also coordinately regulated, either co-induced or co-repressed, and those up-regulated showed strongest changes in ra3 (Fig. 3C; Supplemental Table S7).

The expression signature of cluster 8 is consistent with production of determinate SPMs in 1-mm ears and stg2 tassels. We hypothesized that reduced expression of determinacy factors in this cluster by loss of ra gene function would promote indeterminate fate in SPMs. As predicted, the suite of down-regulated genes in cluster 8 showed strongest reductions in expression in the ra1 mutant, consistent with its extreme phenotype, compared with less severe phenotypes of ra2 and ra3 mutants (Fig. 3F). Among genes with this expression pattern, many function in promoting determinacy and differentiation, and some were previously shown to form modular units in other developmental contexts, e.g., flower development in Arabidopsis and leaf development in maize and other species (Supplemental Table S8). In addition to being co-expressed during normal development, many of these determinacy factors were also co-expressed in mutant backgrounds (Fig. 3G). This module of co-expressed genes suggests conserved functions for determinacy factors co-opted for SPM development in maize and identified grass-specific genes of unknown function that may confer specificity in grass inflorescence morphology (Fig. 3G; Supplemental Tables S6, S8). The latter include GRMZM2G130354, which is shown in Figure 3G, and others predicted to be lineage-specific based on taxonomic dating (Supplemental Table S6; Supplemental Methods).

Among genes in this module, the MYB TF rough sheath 2 (rs2), is expressed in lateral organ primordia and their initials, and acts to promote organogenesis (Timmermans et al. 1999). We identified a co-expressed gene encoding a glutaredoxin with sequence similarity to ROXY1 from Arabidopsis, implicated in petal development (Xing et al. 2005). This gene showed localized expression in the same pattern as rs2 in developing spikelet primordia (Fig. 3G,H; Timmermans et al. 1999). Also co-expressed in this module was a NAC TF, an ortholog of boundary-specifying CUP SHAPED COTYLEDON (CUC) genes from Arabidopsis, which displayed a somewhat different spatial profile, marking boundaries between initiating organs (Fig. 3G,I). This co-expression module represents a signature of determinate SPMs, which is dependent on a functional RA1. We also define a module for SM initiation by identifying clusters that specifically target differences in the 1- to 2-mm transition between wild-type and ra mutant ears (Supplemental Note; Supplemental Fig. S4; Supplemental Tables S9, S10). The SM initiation module included many genes involved in auxin and ethylene biosynthesis and signaling, as well as multiple members of the NON-PHOTOTROPIC HYPOCOTYL (NPH) and LIGHT-DEPENDENT SHORT HYPOCOTYLS (LSH) gene families.

To identify putative upstream regulators of genes in the SPM determinacy module, we used position weight matrices (PWM) to identify cis-regulatory elements enriched in proximal promoters. For 379 co-expressed genes from cluster 8 that were DE in ra1 mutants, we identified 31 significantly enriched motifs (Fig. 3J; Supplemental Table S11). Among these were putative binding sites for TFs that function in integrating environmental cues: i.e., ABF and ABRE motifs, which are bound by stress-responsive bZIP TFs; and LEAFY (LFY), a key factor in modulating the floral transition (Weigel et al. 1992; Bomblies et al. 2003). Motifs recognized by MADS-box TFs were also enriched; e.g., AGAMOUS (AG) and SEPELLATA3 (SEP3) (Fig. 3J; Supplemental Table S11; Kaufmann et al. 2009). Recent work from Arabidopsis showed that bZIPs and MADS-box TFs are likely co-factors of LFY in promoting floral differentiation (Winter et al. 2011; Wu et al. 2012) and thus may play qualitatively similar roles in maize in promoting branch determinacy.

RAMOSA1 functions in activation and repression of co-expressed target genes

To identify targets of RA1 and to distinguish direct vs. indirect interactions, we performed chromatin immunoprecipitation (ChIP)-seq and compared the results to the gene networks described above. Plants expressing complementing RA1 transgenes tagged with HA or YFP were used in parallel experiments (Supplemental Fig. S5; Supplemental Table S12). As expected, YFP-RA1 was expressed in an adaxial domain subtending the SPMs in developing inflorescences (Fig. 4A). HA-RA1 expression was confirmed by detection of a ∼30-kDa fusion protein from immature ear extracts (Supplemental Fig. S5). Ear and tassel primordia were collected and tag-specific antibodies were used to pull down RA1 bound to its target loci. Genome-wide analysis of RA1 occupancy revealed thousands of putative binding sites (i.e., peaks significantly enriched [P < 1 × 10⁻⁵] compared with input DNA). Here we considered 2105 high-confidence peaks, which either overlapped in HA- and YFP-tagged libraries from ear or tassel and/or overlapped in three of the four libraries. High-confidence peaks were strongly enriched between −1.5 and +1.5 kb from the transcriptional start site (TSS) of genes (Fig. 4B); however, RA1 bound in various genomic contexts, frequently introns and exons, and distal intergenic regions (Fig. 4C; Supplemental Figs. S6, S7; Supplemental Tables S13, S14).

Genes within 10 kb of high-confidence peaks were considered putative targets of RA1. Among 1094 putative targets, 22% were DE in ra1 mutant ears and considered directly modulated (Supplemental Table S15). Of these, 70% were down-regulated, suggesting positive regulation by RA1. This was unexpected since RA1 has two ethylene-responsive element binding factor-associated amphiphilic repression (EAR) motifs, which can confer strong repressive activity to TFs (Vollbrecht et al. 2005; Gallavotti et al. 2010; Kagale and Rozwadowski 2011). However, our data also indicate that RA1 acts as a repressor, since 73 bound genes were significantly up-regulated in ra1 mutants. Among repressed targets of RA1 were a number of genes associated with signal transduction and nucleic acid-related processes, such as chromatin and RNAi (Fig. 4D; Supplemental Table S15).

RA1 also bound and modulated genes with known effects on inflorescence development including tasselseed2, encoding an enzyme involved in sex determinacy (DeLong et al. 1993), and compact plant2, the alpha subunit of a heterotrimeric G-protein that functions in IM maintenance (Fig. 4E; Bommert et al. 2013). TFs implicated in cell specification were generally activated by RA1; however, the SBP TF liguleless1 (lg1) was repressed (Fig. 4E,F). lg1 regulates leaf architecture (Moreno et al. 1997; Tian et al. 2011), and loss-of-function lg1 mutants display upright tassel branches (Bai et al. 2012). Although a branch number phenotype has not been reported for lg1 mutants, a recent genome-wide association study identified lg1 as a candidate locus controlling tassel branch number (Brown et al. 2011). Here we show that LG1 protein accumulated in a domain at the base of indeterminate long branches in developing tassel primordia (Fig. 4I), but not determinate SPMs of tassels or ears (Fig. 4G; Supplemental Fig. S8). The LG1 expression domain overlaps that of RA1; however, lg1 is expressed only in the absence of RA1, as shown in ra1 mutant ears (Fig. 4H), supporting direct repression of lg1 by RA1.

RA1 appears to maintain control of its targets during development since 91% of modulated targets were expressed or repressed more differentially at 2 mm compared with 1 mm (Fig. 4J). Strikingly, targets of RA1 were co-expressed in distinct spatiotemporal clusters (from Fig. 3A) depending on whether they were activated or repressed (Fig. 4K,L; Supplemental Table S15; Supplemental Fig. S9). This suggests that RA1 acts to promote or repress gene expression in a manner dependent on spatiotemporal context, possibly through interactions with multiple cofactors.

RA1 is characterized by a single zinc finger domain, implicated in binding DNA through a consensus AGT sequence (Dathan et al. 2002). To determine whether additional sequences were associated with RA1 binding, we tested for enrichment of de novo motifs under high-confidence peaks in target gene promoters. Most notable was a GAGA repeat, strongly enriched adjacent to the mean peak center (Fig. 4M,N). Drosophila GAGA factor (dGAF), which activates and represses gene expression primarily through association with chromatin-associated proteins, shares several structural and functional features with RA1 (Omichinski et al. 1997; Dathan et al. 2002; Vollbrecht et al. 2005). One hypothesis is that RA1 may have been co-opted in Panicoid grasses to control developmental transitions using a mechanism analogous to that of dGAFs along with specific cofactors.

Additional de novo motifs enriched in regions surrounding peak centers, i.e., potential binding sites for RA1 cofactors, were associated with genes involved in flower development and floral transition (Fig. 4M,N; Supplemental Fig. S10). One motif was significantly similar (P = 5.6 × 10⁻⁹) to the in vitro validated binding site for indeterminate1 (id1), a regulator of the floral transition in maize (Colasanti et al. 1998). Little is known about the other INDETERMINATE DOMAIN (IDD) genes in maize; however, it was shown that at least some can also bind the id1 site (Kozaki et al. 2004; Colasanti et al. 2006). RA1 binds and positively modulates three IDD genes (Supplemental Table S15).

Common targets of RA1 and KNOTTED1 link meristem determinacy and maintenance networks

While RA1 imposes determinacy on SPMs, KNOTTED1 (KN1) is a master regulator of indeterminacy in all shoot meristems. Experimental and genetic evidence suggest that RA1 and KN1 proteins interact, and their expression domains overlap in a boundary at the base of SPMs (E Vollbrecht, unpubl.). We used the KN1 data set to determine the extent of RA1 and KN1 co-occupancy and modulation of shared target genes (Bolduc et al. 2012). We identified 189 target loci associated with regions where high-confidence RA1 and KN1 peaks overlapped. An additional 292 putative targets were bound by RA1 and KN1, but at different sites. Total shared targets were significantly greater than expected by chance (P < 2.2 × 10⁻¹⁶) (Fig. 5A; Supplemental Table S16). Of the 176 (37%) targets that were modulated, 79% and 39% were DE in ra1 and kn1 loss-of-function mutants, respectively (Fig. 5A,C). This suggests a stronger dependence on RA1 for normal expression of shared targets, possibly due to compensation for KN1 by other knox genes (Fig. 5C; Farnham 2009; Bolduc et al. 2012). RA1 itself may also be involved in kn1 regulation, since it bound the 5′ region of kn1's third intron (Fig. 5B), which is rich in conserved noncoding sequences involved in kn1 repression (Inada et al. 2003). KN1 was shown to bind and auto-regulate itself through the same intron (Tsuda et al. 2011; Bolduc et al. 2012). Although expression of kn1 was not significantly altered in ra1 mutants in our experiments, regulation could be restricted to a small number of cells, since the ra1 expression domain is much smaller than that of kn1.

Shared targets of RA1 and KN1 were enriched for genes encoding TFs (GO:0045449; P = 4 × 10⁻⁴), and displayed several unique expression profiles (Fig. 5D). Among co-modulated TF family members were three HD-Zip class I genes (Fig. 5D,E), which were positively modulated by RA1 and/or KN1. Members of this class have been implicated in integrating responses to the environment with development (Ariel et al. 2007; Whipple et al. 2011). Several other TFs co-bound by RA1 and KN1 function in circadian response and modulation of floral transition (Fig. 5D; Supplemental Fig. S11). For example, RA1 and KN1 co-occupy the promoter of a MYB TF orthologous to LATE ELONGATED HYPOCOTYL/CIRCADIAN CLOCK ASSOCIATED1 (LHY/CCA1), and the intron of bZIP TF liguleless2 (lg2), two genes that function in the transition from vegetative to reproductive development (Walsh et al. 1998; Mizoguchi et al. 2002). Although ra1 is expressed after the transition to inflorescence development, this binding could reveal feedback regulation of the transition. Co-occupancy of RA1 and KN1 at lg2 is also notable since loss-of-function mutants in lg2 display a decreased tassel branching phenotype (Walsh and Freeling 1999), further implicating the liguleless module in inflorescence architecture.

Two of the IDD TFs modulated by RA1 were also bound by KN1 (Fig. 5F). Both RA1 and KN1 bound the promoter of ZmIDD1-p1 (Colasanti et al. 2006), but at distinct sites, and modulated the gene antagonistically. RA1 and KN1 co-occupied the same site in the promoter of another IDD gene, an ortholog of LOOSE PLANT ARCHITECTURE 1, which regulates tiller angle in rice (Wu et al. 2013), but loss of kn1 function did not alter its expression. RA1 also independently bound and modulated an ortholog of ENHYDROUS, an IDD TF implicated in integration of hormone signals in Arabidopsis (Feurtado et al. 2011). Evidence for combinatorial action of IDD genes in developmental patterning has recently emerged (Welch et al. 2007; Reinhart et al. 2013; Slewinski 2013), including novel roles in lateral organ development via regulation of auxin biosynthesis and transport (Cui et al. 2013). One possibility is that RA1 may modulate auxin levels and distribution in the developing inflorescence indirectly via regulation of IDD genes.

KN1 was shown to play a major role in integration of various hormone signaling networks (Bolduc et al. 2012). While only 3% of genes bound and modulated by RA1 were classified as hormone related, functional analysis showed significant enrichment of gibberellic acid (GA)-mediated signaling (P = 5.9 × 10⁻⁴). RA1 and KN1 co-bound a ga2-oxidase (ga2ox1, GRMZM2G427618) and an ortholog of the Arabidopsis SPINDLY (SPY) gene, suggesting they interface at both GA biosynthesis and signaling (Fig. 6A). In addition, a number of auxin-related genes were co-targeted by RA1 and KN1 (Supplemental Fig. S11). Our analysis of DR5 expression suggested that auxin response in the meristem marks indeterminate fate. Therefore, RA1 and KN1 could work cooperatively and/or antagonistically to modulate auxin signals during development.

Plant material for RNA-seq experiments

Segregating families (1:1) of ra1-R, ra2-R, and ra3-fea1 mutant alleles, all introgressed at least six times into the B73 inbred background, were grown at CSHL Uplands Farm. Field-grown plants were genotyped and collected 6–7 wk after germination (V7–V8 stage). First and second ear primordia were immediately hand-dissected, measured, and frozen in liquid nitrogen. For ra1, ra2, and ra3 mutants and wild-type controls, ears were pooled into two size classes: (1) 1-mm class included a range from 0.7- to 1.5-mm sized ears, and nine ears were pooled for each of two biological replicates; (2) 2-mm class included a range from 1.8- to 2.5-mm sized ears, and six ears were pooled for each of three biological replicates. Wild-type samples were proportional mixtures of heterozygote siblings segregating in ra1, ra2, and ra3 populations. Variability factors (e.g., ear size within class, ear rank on the plant, and time of collection) were distributed evenly across pooled samples.

For the wild-type ear and tassel developmental series, greenhouse-grown B73 inbred plants were used. Ten millimeter ears were collected and sectioned as follows from tip to base along the developmental gradient: tip 1-mm sampled (tip; IM/SPM), next 1-mm discarded, next 2-mm sampled (mid; SM), next 2-mm discarded, next 2-mm sampled (base; FM) (Supplemental Fig. S3), and immediately frozen in liquid nitrogen. Sections from ∼30 sampled ears were pooled for each of two biological replicates to represent tip, mid, and base stages. Tassels were hand-dissected, measured, and separated by stage: 1–2 mm (stg1), 3–4 mm (stg2), and 5–7 mm (stg3); and immediately frozen in liquid N. For each stage, ∼20–30 tassels were pooled for each of two biological replicates.

Microscopy

For scanning electron microscopy (SEM), fresh immature ears and tassels were dissected and immediately imaged without fixation using a Hitachi (S-3500N) SEM. For fluorescence microscopy, maize HiII plants were transformed with a DR5-ER∷RFP reporter construct (Gallavotti et al. 2008) and backcrossed to ra mutants. Tassel and ear primordia were hand-dissected at 4–5 wk (V4–V5) and 6–7 wk (V7–V8) after germination, respectively. Optical sections were taken using the Carl Zeiss LSM 710 Laser Scanning Confocal Microscope. Z-stacks of the optical sections were reconstructed into three-dimensional images and the RFP channel signal was converted into the LUT (mapped color) “Fire” using Bitplane Imaris 7 software. Further details on sample preparation and parameters used are found in Supplemental Methods.

RNA-seq library construction, sequencing, and analysis

RNA-seq libraries were generated from 2 to 5 µg total RNA (RNAeasy kit; Qiagen) using methods adapted from Li et al. ( 2010). Libraries were size-selected for a 250- to 300-bp insert for paired-end (PE) sequencing using standard Illumina protocols (Illumina, Inc.). Libraries were quantified on an Agilent bioanalyzer (Agilent) using a DNA 1000 chip, and sequenced using the Illumina GAII platform at the CSHL Genome Center.

The Tuxedo suite (Trapnell et al. 2010) was used for mapping and analysis of RNA-seq data. Further details are provided in Supplemental Methods. TopHat (version 1.2.0) was used to align reads to the maize reference genome (AGPv2) based on an a priori set of 110,028 predicted maize gene models (Working Gene Set v5b.60; maizesequence.org). Cuffdiff (version 1.0.2) was then used to analyze differential expression using a high-confidence subset of 39,656 maize gene models, the Filtered Gene Set (FGS v5b.60; maizesequence.org). Gene-level expression values are represented by fragments per kilobase exon per million reads mapped (FPKM) and a consensus FPKM was determined for each gene based on its representation across biological replicates (Trapnell et al. 2010). Biological replicates showed strong correlations (r > 0.95) in gene expression. We used a corrected P-value of <0.05 to call differentially expressed (DE) genes between each set of compared samples (Supplemental Table S2). Additional information about gene annotation and functional enrichment analyses is provided in Supplemental Methods.

Analyses were performed the same for previously published RNA-seq data sets from knotted1 loss-of-function ear and tassel samples and B73 wild-type controls (Bolduc et al. 2012), two biological replicates from each, after raw data were downloaded from the Gene Expression Omnibus (GEO; GSE38487).

Clustering of gene expression profiles

A k-means approach was used (R Bioconductor package Mfuzz) to cluster 16,272 dynamically expressed genes based on their expression profiles across wild-type libraries. These genes were selected based on the following criteria: (1) showed a twofold change in expression or were identified as DE between at least two of the stages being compared, (2) were covered by at least 50 reads among every two stages being tested, and (3) collectively had an expression value of at least 1 FPKM. With these criteria, we reduced noise in the clusters by omitting genes that did not change expression during development or showed very low expression levels with little confidence. FPKM values were normalized to a Z-score scale prior to clustering. We evaluated fuzzy k-means results based on 10–35 clusters using 1000 iterations and finally chose 20 clusters based on optimal results (Supplemental Methods). Genes were grouped with their best-fit cluster.

ChIP-seq sample collection, library construction, and peak calling

We created two native translational fusion constructs to drive the expression of tagged RA1 proteins in the endogenous expression domain using 2.9 kb of the RA1 promoter. We fused the YFP and HA-FLAG tags in frame with the RA1 coding sequence at the N terminus. Constructs were transformed into the HiII genetic background at the Iowa State University Plant Transformation Facility. T0 generation transformed plants were crossed to the ra1-R mutant. T1 plants were then backcrossed to create a T2 generation segregating 1:1 for the transgene and for ra1. Tassel primordia were harvested ∼4 wk after planting, and immature ears were harvested after 6 wk. Analysis of plant phenotypes in F1BC2 families segregating for the transgene and ra1 mutants, showed that transgenic constructs were capable of complementing the mutant. In populations fixed for the ra1-R mutation, but segregating for the YFP- or HA-tagged transgene, the wild-type phenotype segregated perfectly with the presence of the transgene (n = 160).

ChIP and ChIP-seq library preparation were performed as previously described (Bolduc et al. 2012; Morohashi et al. 2012) with minor modifications. Further details on ChIP, library preparation and sequencing, and ChIP-seq analysis are provided in Supplemental Methods. ChIP-seq reads were aligned to the maize reference genome (AGPv2) and peak calling was performed with MACS version 1.4.0rc2 using only uniquely mapped reads (Zhang et al. 2008). Peaks were identified as significantly enriched (P < 1 × 10⁻⁵) in each of the ChIP-seq libraries compared with input DNA. Significant peaks from individual libraries were considered overlapping if their summits were positioned within 300 bp of each other. Coordinates for KN1 peaks were used from the published data set (Bolduc et al. 2012); GEO (GSE39161).

Analysis of cis-regulatory motif enrichment

To determine enrichment of cis-regulatory elements within promoter regions spanning 1 kb upstream to 500 bp downstream from the TSS of co-expressed genes, we used a computational prediction pipeline, which leveraged the search tool for occurrences of regulatory motifs (STORM) from the comprehensive regulatory element analysis and detection (CREAD) suite of tools (Smith et al. 2006). We identified putative TF binding sites based on 128 experimentally derived position weight matrices (PWMs) from various sources (Supplemental Methods). We considered only those motifs that were over-represented (P-value <0.001) in promoter sequences of protein-coding genes as compared with a background set of the same number of random genomic sequences. STORM uses this P-value to assign a PWM-specific quality score to each predicted cis-regulatory motif (Schones et al. 2007). Predictions were further filtered based on a PWM-specific median score threshold (i.e., quality score greater than or equal to the median score passed the filter) and a motif occurrence frequency of three or more per promoter. A parallel analysis was carried out for all FGS genes in the genome to compare their overall percent occurrence of enriched motifs in the maize genome to determine the extent of cluster specificity for a given motif.

To identify de novo motifs associated with RA1 ChIP-seq peaks, we used the Promzea pipeline (Liseron-Monfils et al. 2013; Supplemental Methods). To optimize signal-to-noise, we used only those high-confidence RA1 ChIP-seq peaks that were associated with proximal promoters of FGS models within 1 kb upstream of TSSs. Enriched motifs were identified and tested for significance by comparing their presence in the input data set (sequences spanning RA1 peaks) with 5000 random maize promoter sequences (500 bp upstream of predicted TSSs). Additional details are provided in Supplemental Methods.

In situ hybridizations

Immature B73 ears (2–5 mm) were sampled as above. Details on sample preparation and in situ hybridizations are provided in Supplemental Methods. The glutaredoxin gene (GRZM2G442791) probe included the entire coding sequence and 241 bp of the 3′UTR. The CUC-like gene (GRMZM2G393433) probe contained 700 bp of nucleotide sequence, including 115 bp of the 3′UTR. Probe sequences were cloned into the TOPO PCRII dual promoter system (Invitrogen) and linearized for probe transcription (Roche).

LG1 immunolocalization

B73 wild-type ears and tassels, and ra1 mutant ears (all 2–5 mm) were sampled as described above and fixed in PFA under vacuum infiltration, dehydrated through an ethanol series into Histoclear, and embedded in paraplast plus. Tissue was sectioned to 10 microns using a Leica microtome and mounted on Probe-on Plus slides. Further details for the immunolocalization are provided in Supplemental Methods. The antibody against full-length LG1 expressed as a HIS-tagged N-terminal protein fusion in E. coli was created as previously described (Chuck et al. 2010) (Cocalico Biologicals, Inc.) and purified against full-length LG1 expressed as a GST-tagged N-terminal fusion protein. Primary α-LG1 (guinea pig) was used at a 1:500 dilution, and α-guinea pig alkaline phosphatase conjugated secondary antibody was used at a 1:5000 dilution (Abcam).