A high-resolution gene expression map of the Arabidopsis shoot meristem stem cell niche

A spatiotemporal orchestration of gene expression and cell behaviors is crucial for cell fate specification (Baurle and Laux, 2003; Reddy, 2008). The stem cells in shoot apical meristems (SAMs) of higher plants differentiate into various above-ground organs (Steeves and Sussex, 1989). Traditionally, SAMs have been divided into distinct domains on the basis of their function and cytological criteria (Lyndon, 1998; Steeves and Sussex, 1989; Xie et al., 2009). The central zone (CZ) harbors a set of pluripotent stem cells (Reddy and Meyerowitz, 2005). The stem cell progeny differentiate into leaves or flowers within the adjacent peripheral zone (PZ). Upon specification of leaf or flower primordia, specific cellular behaviors and differentiation events lead to the formation of boundary regions that separate developing organs from SAMs, whereas cells in the rib meristem (RM), situated just below the CZ, differentiate to give rise to the stem. Collectively, stem cell daughters progress through sequential differentiation programs and acquire distinct cellular identities (Reddy et al., 2004; Heisler et al., 2005). SAMs of higher plants are multilayered structures. In Arabidopsis, the tunica consists of the outer epidermal (L1 layer) and an inner sub-epidermal (L2 layer), while the corpus forms a multilayered structure beneath the L2 layer (Satina and Blackeslee, 1941; Satina et al., 1940; Tilney-Basset, 1986). Thus, the SAM contains cell types that are clonally distinct and cells that are at different stages of differentiation within the same cell layer. Both cell type-specific factors and ubiquitously used signals such as plant hormones have been implicated in SAM development (Xie et al., 2009; Perales and Reddy, 2012; Barton, 2010; Shani et al., 2006; Yadav et al., 2011). Elucidating gene networks underlying cell fate specification requires methods that facilitate the gene expression analysis of individual cell types.

An earlier study describes gene expression profiles of the CZ, the PZ and the RM (Yadav et al., 2009); however, it does not resolve gene expression programs associated with cell types in different clonal layers or with specific specialized cell types. To gain further insights into transcriptional programs of SAM cell types, we identified novel promoters that are expressed in discrete clonal layers and in different cell types of CZ, PZ and RM. Fluorescence-activated cell sorting (FACS) method was employed to collect seven new cell populations and Affymetrix GeneChips (ATH1) were used to profile their transcriptomes. These data were combined with expression data of three cell populations from an earlier study (Yadav et al., 2009) to generate a high-resolution gene expression map. Besides the description of molecular pathways enriched in various cell populations, our work reveals that cells in the distinct clonal layers exhibit greater diversity in gene expression and transcriptional complexity than cell populations within the CZ and the PZ.

The HIGH MOBILITY GROUP (PHMG; AT1G04880) promoter was used to label cells of the L1 layer of SAMs (Table 1 and Fig. 1A; HMG/L1-meristematic). The Arabidopsis thaliana MERISTEM LAYER1 (AtML1) promoter was used to mark cells in the L1 layer of both the meristematic region and differentiating organs (Fig. 1B,F; AtML1/L1-ubiquitous) (Lu et al., 1996). The promoter of the HOMEODOMAIN GLABROUS4 (PHDG4, At4g17710) gene was used to label cells in the L2 layer (Table 1; Fig. 1C,G; HDG2/L2-meristematic) (Yadav et al., 2009). The microarray data of WUSCHEL-expressing cells (Yadav et al., 2009) were used to represent cells of the corpus/L3 layers/the RM (Table 1; WUS/RM/L3/corpus) (Mayer et al., 1998). The phloem pole pericycle reporter line S17 (At2g22850) was used to label the phloem (Fig. 1J; see supplementary material Fig. S1B) (Lee et al., 2006), and the S6/AtHB8 reporter was used to mark procambial cells (Fig. 1K; see supplementary material Fig. S1A) (Baima et al., 1995).

The CZ (PCLV3::mGFP5-ER) gene expression dataset was obtained from an earlier study, referred to as CLV3/CZ (Yadav et al., 2009). The KANADI1 (KAN1) promoter was used to label cells in the outer edges of the PZ, referred to as KAN1/outer PZ (Fig. 1E,I; green fluorescent signal) (Yadav et al., 2013; Kerstetter et al., 2001). The organ primordia gene expression data set (PFIL::dsRED-N7) has been described in an earlier study, referred to as FIL/organ primordia (Table 1; Fig. 1E,I; red fluorescent signal) (Yadav et al., 2009). The LATERAL SUPPRESSOR (LAS) promoter was used to isolate cells in the organ boundaries, referred to as LAS/organ boundaries (Table 1; Fig. 1D,H) (Goldshmidt et al., 2008; Greb et al., 2003). Seven new fluorescent reporters were introduced into apetala1;cauliflower1-1 (ap1-1;cal1-1) genetic background and then we performed protoplasting to disrupt the SAM, employed FACS to isolate specific cell types and used the Affymetrix GeneChip ATH1 to perform microarray experiments. The new data sets were collated with the expression profiles of three cell populations that have been reported previously (Yadav et al., 2009) and were analyzed to identify gene expression programs within the shoot apex.

Based on results from our previous study (Yadav et al., 2009), the protoplasting method influences the expression of ∼300 genes. Thus, these genes were excluded from all subsequent analysis steps. Cell layers in the SAM traverse through the CZ and PZ, whereas vasculature cell types partially overlap with deeper cell layers of the RM. Consequently, the reporters used to isolate specific cell populations have overlapping expressions. To minimize gene expression redundancies arising from these overlaps in our analyses of differentially expressed genes (DEGs), we assigned the ten cell population samples to four non-overlapping comparison groups. Group A included the HMG/L1-meristematic, HDG4/L2-meristematic and WUS/L3/corpus cell populations from the SAM cell layers. Group B included the CLV3/CZ, FIL/organ primordia, LAS/organ boundaries and KAN1/outer PZ. Group C included AtHB8shoot, S17shoot and HDG4/L2-meristematic cell populations, and is expected to identify genes expressed in vascular cell populations in comparison with the L2 cell layer. Group D included AtML1 (L1 ubiquitous), HMG (L1 meristematic) and HDG4 (L2 meristematic) cell populations, and is expected to identify genes expressed in the L1 layer of meristematic and non-meristematic cells in comparison with the L2 cell layer. Within those comparison groups we identified DEGs for all possible sample comparisons using the LIMMA package (Smyth, 2004). The DEG results were then queried for genes that showed an elevated expression in one cell type (target) compared with all other cell types (reference) in a comparison group. The resulting ‘cell population differentially expressed genes’ (CPEGs) needed to show a ≥1.5-fold higher expression in the target cell type compared with the reference cell types in a comparison group. In addition, all fold-changes needed to be supported by an adjusted P value from LIMMA of ≤0.05. The four groups A, B, C and D showed 1456, 472, 1031 and 298 CPEGs, respectively (Fig. 2A; see supplementary material Table S1).

Our previously reported experiments involving three cell populations, CLV3, WUS and FIL, revealed a greater number of CPEGs in comparison with the current work (Yadav et al., 2009). This could be due to overlapping expression domains of reporters used in this study, as seen in the case of FIL and KAN1 (supplementary material Fig. S2). A more relaxed version of the CPEGs analysis used the same threshold values as above but allowed one outlier in the target to reference comparisons. For example, a gene was considered a relaxed CPEG if it met those criteria relative to at least three out of a total of four reference cell populations: in this study, the CZ and the PZ domain datasets. In the relaxed CPEG analysis, the CZ and the PZ domain dataset yielded 2182 genes as opposed to the 472 detected without an outlier (see supplementary material Table S1). Many of these newly detected CPEGs were part of PZ cell types, thereby confirming the utility of the outlier analysis.

To validate the cell population transcriptome, we compared the cell type-elevated expression of 91 genes in cell populations with RNA in situ hybridization patterns determined by Yadav et al. (2009) and other studies (see supplementary material Fig. S3 and Table S8). For example, 14 L1 CPEGs were found to be expressed in the L1 layer/epidermis (see supplementary material Table S1). Moreover, TRYPTOPHAN AMINOTRANSFERASE OF ARABIDOPSIS1 (AT1G70560) and UNKNOWN PROTEIN (AT5G06270) were CPEGs in both the L1/HMG (layer-specific analysis data set) and the CLV3 cell population, which is in agreement with the relatively higher levels of expression of these genes in the central part of the L1 layer and to a lesser degree in the L2 sub-epidermal layer (see supplementary material Fig. S4). Among five L2 layer expressed genes, four CPEGs were identified (see supplementary material Table S1). The absence of MEI2 C-TERMINAL RRM LIKE2 (MCT2) (AT5G07930) in the L2 CPEGs may be attributed to the dynamic expression of MCT2 in the flower meristem in addition to its enrichment in the L2 layer of SAMs (Aggarwal et al., 2010). Furthermore, the expression patterns of 12 RM genes were correctly assigned to the WUS cell population (see supplementary material Table S1). Collectively, the layer-specific analysis correctly identified 29 CPEGs out of 30 known candidate genes considered for analysis.

A similar analysis was performed with datasets of the CZ and the PZ domains. Of the 15 genes expressed in the CZ, 13 were identified as CPEGs in CLV3 cell population. One of the genes excluded by the analysis, TERMINAL EAR LIKE2 (AT1G67770), is expressed in the CZ and by the PZ of the early stage flower buds (Anderson et al., 2004) (see supplementary material Table S1). Of the 31 genes that are expressed in the PZ, our CPEG analysis assigned 24 to the FIL cell population. Exceptions include the expression patterns of AUXIN RESPONSE FACTOR3 (AT2G33860) and LATE MERISTEM IDENTITY2 (AT3G61250), which are expressed more diffusely in the PZ (see supplementary material Table S1). Of the seven genes expressed in the outer edge of PZ, six were assigned to the KAN1 cell population. Similarly, nine out of 11 genes showed elevated expression in boundary regions in the LAS cell population (see supplementary material Table S1). In summary, 41 of the 46 CPEGs agreed with the RNA in situ hybridization data. The disagreeing five candidates may have a too dynamic expression pattern to be detectable by RNA in situ hybridization experiments. Last, the expression of all ten genes that were considered for the vascular cell populations were correctly assigned to their respective cell populations in the CPEG analysis. Taken together, our microarray experiments correctly identified 94% of the chosen CPEGs (see supplementary material Tables S1 and S8).

To estimate the impact of the fold change cutoffs on the number of CPEGs identified, we tested four different fold-change cutoffs, including 0, 1.5, 2 and 4, while using in all cases a constant adjusted P value cutoff of ≤0.05. As expected, larger fold change cutoffs result in smaller CPEG sets (supplementary material Table S11). Examples of genes that were CPEGs at the 1.5-fold but not at 2-fold cutoff anymore include HDG11, SPT, MCT1, MCT2, CLV1, BAM1, ARF3, LFY, BOP1, LMI1, PRS, KAN1 and KAN2 (supplementary material Tables S11 and S12). One reason for this trend is that many of the weaker expressed genes are unlikely to show a pronounced fold change in Affymetrix experiments, mainly because of limitations related to the sensitivity and the signal-to-background ratio of the technology. Thus, a relatively small fold change of 1.5 can be biologically meaningful, especially for weakly expressed genes (for example KAN1 and KAN2). This is also supported by the fact that the CPEG sets for the different fold change cutoffs shared many statistically enriched Gene Ontology (GO) terms (see below). Additional reasons could be the relatively broad expression pattern of candidate genes (e.g. CLV1, BAM1), which does not perfectly overlap with the expression domains of reporters used for cell sorting.

To analyze the CPEG sets functionally, we performed enrichment analysis of Gene Ontology (GO) terms (Fig. 2D; see supplementary material Tables S2 and S12). The GO term enrichment tests were performed on the CPEG sets generated for each of the four-fold change cutoffs (see supplementary material Table S12) (see above for details), while the following discussion is restricted to the GO terms enriched (adjusted P≤0.05) in both result sets, the one with the 1.5- and the 2-fold cutoff.

The CPEG datasets from the SAM L1 layer contained several over-represented GO terms related to L1 layer development (GO:0032502), defense (GO:0006952) and responses to bacterial and fungal pathogens (GO:0009607) (Fig. 2D; see supplementary material Table S2 and Fig. S5). This included: genes involved in wax biosynthesis, such as 3-KETOACYL-COA SYNTHASE1 (KCS1), KCS5, KCS6/CUTICULAR1, KCS10/FIDDLEHEAD and WAX2 (Javelle et al., 2011); genes implicated in L1 layer development, such as HOMEODOMAIN GLABROUS2 and 12, ARABIDOPSIS CRINKLY4, PROTODERMAL FACTOR 1 and 2 and AtML1 (Abe et al., 2003); genes that have been shown to function in the fine-tuning of the basal resistance against Pseudomonas syringae-WRKY11 and WRKY17 (Journot-Catalino et al., 2006); and three members of the Toll-like/interleukin-1 (TIR-NB) gene family involved in eliciting responses to fungal pathogens. The cells of the L1 layer also contained CPEGs annotated as auxin influx carrier AUX1 and its paralogs, LAX1 and LAX2. Conversely, the CPEGs of the non-meristematic cell populations of the epidermal layer were enriched in pathways related to photosynthesis (GO:0015979) and flavonoid biosynthesis (GO:0009813) (see supplementary material Table S2) (Zhao et al., 2007).

The developmental roles of the L2 layer are not well understood. The CPEG set of this sample was enriched in GO terms related to the metabolic processes of nucleic acids (GO:0006139) (Fig. 2D; see supplementary material Table S2). The corresponding genes include: POLYMERASE DELTA4 (POLD4), a DNA-directed DNA polymerase; DNA glycosylase family protein (At1g80850); URACIL DNA GLYCOSYLASE enzymes that are involved in DNA excision repair; AtMND1, which plays a role in double-strand break repair during meiotic recombination (Cappadocia et al., 2010); and genes encoding TELOMERASE REVERSE TRANSCRIPTASE, which is involved in the maintenance of chromosomal ends (Ogrocká et al., 2012), and WHIRLY1, which is involved in the negative regulation of telomere length. Furthermore, the CPEGs identified for the L2 layer cells include: DOMAINS REARRANGED METHYLASE1 (DRM1), which is involved in de novo DNA methylation and the maintenance of asymmetric DNA methylation; and INCREASE IN BONSAI1 METHYLATION, a jumonji domain-containing protein involved in H3mK9 demethylation (Stroud et al., 2013). However, at twofold enrichment cutoff, POLD4, DNA GLYCOSYLASE and WHIRLY1 were no longer identified as CPEGs. The preferential enrichment of pathways involved in maintaining DNA fidelity and telomere length within the L2 layer suggests that they may provide necessary safeguards for germ cells that descend from this cell layer (Scott et al., 2004).

The CPEGs identified for the L2 layer sample contained many genes involved in microtubule organization. They included: AtCDC48, which encodes ATP binding microtubule motor protein; two genes that encode P-loop containing triphosphate hydrolase superfamily proteins; and members of the kinesin family, such as phragmoplast-associated kinesin-related protein (PAKRP2) and ARABIDOPSIS THALIANA KINESIN4 (ATK4). PAKRP2, a phragmoplast-localized protein has been implicated in transport of cargo from the Golgi to the division site during cell-plate formation (Lee et al., 2001). Whereas ATK4 is a kinesin-like protein that binds microtubules in ATP-dependent manner to facilitate microtubule-based movement of vesicles. At the twofold cutoff, PAKRP2 and ATK4 were not detected as CPEGs. The enrichment of GO terms related to pathways controlling microtubule organization and dynamics among cells within the L2 sub-epidermal layer is intriguing, considering their ubiquitous requirement (Zhong et al., 2002).

The CPEG set of the cells of the L3 layer was enriched in GO terms related to photosynthesis (GO:0015979) (Fig. 2D; see supplementary material Table S2). Many of these CPEGs encode members of photosystem I light-harvesting complex A (LHCA1, 3, 4 and 5), photosystem II components (LHCB3, 4, 4.2 and 5), proteins of the photosystem I subunits PSAH1, PSAE1, PASF, PSAA and PSAO, and proteins of photosystem II subunits PSBA and PSBO2 (Friso et al., 2004).

The GO categories related to responses to ions (GO:0010038) and salt stress (GO:0009651) were also enriched in the CPEG set of the L3 cell layer. This includes genes encoding FERRETIN1 (an iron-storage protein); AILP1 (which responds to aluminium ions), CAX1-LIKE CATION EXCHANGER 1 and 3, a Ca2+/H+ exchanger; ANNEXIN1, a Ca²⁺-dependent protein with peroxidase activity; a vacuolar Ca²⁺-binding protein that responds to cadmium ions and salt stress (Sudre et al., 2013); protein kinases of the SNF1 family (SNRK2.1, SNRK2.2, SNRK2.6); ABI FIVE BINDING PROTEIN4 (AFP4); and MYB15 (Fujii et al., 2011). At the twofold cutoff, AILP1 and AFP4 were not identified as CPEGs. In addition, the CPEG set of the L3 cells contained GO terms associated with hormonal responses (GO:0009725) that included type B Arabidopsis response regulators (ARR2, ARR4 and ARR10), which may explain enhanced cytokinin sensitivity of RM cells observed in an earlier study (Gordon et al., 2009). In addition, when we applied the BINGO tool on the L3 layer CPEGs using Cytoscape, we found that the L3 cell population responds to the majority of plant hormones, including jasmonic acid, ethylene, cytokinin, abscisic acid and gibberellins (Fig. 3; see supplementary material Fig. S5) (Maere et al., 2005). Collectively, the cells in the L3 layers express genes implicated in maintaining ionic homeostasis, mediating the water stress, the salt stress and the response to multiple plant hormones.

CZ cells primarily express genes involved in stem cell maintenance (GO:0019827), cell differentiation (GO:0030154) and floral organ development (GO:0048437) (Fig. 2D; see supplementary material Table S2). This includes CLAVATA3, WUSCHEL, SAPTULA (SPT), CORONA and PHABULOSA (Barton, 2010). However, at the twofold enrichment cutoff, SPT was not detected as a CPEG anymore.

The GO terms related to organ development (GO:0048513), specification of organ polarity (GO:0010084) and responses to hormonal stimuli (GO:0009725) were enriched in the CPEG set of the organ primordia (see supplementary material Table S2). This CPEG set includes FILAMENTOUS FLOWER, ASYMMETRIC LEAVES1, GROWTH REGULATING FACTOR7, ROXY1, BLADE ON PETIOLE2 and MONOPTEROS/ARF5. It also includes the hormone-responsive genes TONOPLAST INTRINSIC PROTEIN1, AUXIN RESISTANT2, ARF5, ARABIDOPSIS RESPONSE REGULATOR11, AtbZIP12, HOMEOBOX-LEUCINE ZIPPER2, AtMYB13 and BASIC CHITINASE (Barton, 2010).

The CPEG set of the KAN1 cell population was enriched in GO terms related to auxin homeostasis (GO:0010252) (Staswick et al., 2005), wax biosynthesis (GO:0010025), fungal responses (GO:0050832) and defense responses (GO:0006952) (Fig. 2D; see supplementary material Table S2) (Consonni et al., 2006). The CPEGs that encode two auxin-responsive-GH3 family proteins and SLOW MOTION, a F-box protein implicated in the timing of lateral organ development were found in the KAN1 cell population (Lohmann et al., 2010).

Cells of the organ boundary exhibit unique cell behaviors, including infrequent cell division and asymmetrical expansion that create a region of reduced growth (Reddy et al., 2004). GO terms related to hormonal responses (GO:0009725), including brassinosteroid (BR) (GO:0009741), were associated with this cell population (Fig. 2B; see supplementary material Table S2) (Hu et al., 2004). This included CPEGs that encode CYCLIN D3;1, a D-type CYCLIN, and Xyloglucan Endotransglycosylase (XET), which mediates cell wall expansion by hydrolyzing xyloglucan chains. At twofold enrichment cutoff, CYCLIN D3;1 was not retained as a CPEG.

The GO terms related to auxin response and signaling pathways (GO:0009733), and polar auxin transport (GO:0009926) were enriched in the CPEG set of the xylem cell population (Fig. 2D; see supplementary material Table S2) (Vernoux et al., 2011). This includes auxin-induced genes [BODENLOS (BDL), MONOPTEROS, IAA8, IAA13 and IAA30], genes encoding GH3.3 family proteins [DFL1, AUX1, AUXIN INDUCED IN ROOT CULTURES3 (AIR3), AIR12, AXR4, FLAVODOXIN-LIKE QUINONE REDUCTASE1, PINOID-BINDING PROTEIN1 and BRI1-LIKE2] and genes involved in polar auxin transport [AtPIN1, AtPIN6, ENHANCER OF PINOID (ENP), AUX1, ACAULIS5 (ACL5) and MORE AXILLARY BRANCHES2 (MAX2)]. At the twofold and higher cutoffs, BDL, IAA8, AXR4, ENP, ACL5 and MAX2 were no longer detected as CPEGs. An enrichment of the auxin pathway elements in cells that represent the early stages of xylem differentiation is consistent with the documented role of auxin in vascular differentiation (Sachs, 1969,, 1981).

The CPEG set of the phloem cell population contained enriched GO terms related to glucosinolate biosynthesis (GO:0019761) (Piotrowski et al., 2004), to gibberellin- and bacterium-induced responses (GO:0009739; GO:0009617) (Cheng et al., 2004), to cell wall modification (GO:0042545), and to ion transport (GO:0006811) (Fig. 2D; see supplementary material Table S2) (Grotz et al., 1998). Genes involved in cation homeostasis, such as CYCLIC NUCLEOTIDE-GATED ION CHANNEL 5 and 19 (CNGC5 and CNGC19), K+ UPTAKE PERMEASE9 (AtKUP9), COPPER TRANSPORTER5 (COPT5), ZINC TRANSPORTER 1 and 5 PRECURSOR (ZIP), iron-regulated transporters IREG1 and IRT3, nitrate transporters NRT1.5 and NRT1.7, the phosphate transporter PHO1 and sulfate transporter-SULTR2;1, were identified as CPEGs in the phloem sample (Lin et al., 2008). However, at the twofold and higher cutoffs, CNGC5, AtKUP9 and ZIP1 were no longer detected as CPEGs.

We performed the sample-wise hierarchical clustering of transcriptome to determine similarities in the gene expression profiles among cell types. Data from the ten SAM cell populations as well as from 16 cell populations of the Arabidopsis root were considered (Birnbaum et al., 2003; Brady et al., 2007). The distance matrix required for this step was generated from the distance-transformed Spearman correlation coefficients obtained from all pair-wise sample comparisons using the average among their normalized replicates (see supplementary material Table S6). The resulting clustering dendrogram (Fig. 4C) indicates that the expression profiles of adjacent cell types in the CZ and the PZ are more similar than those that are spatially distant. For example, FIL cell population share greater similarities with KAN1 cell population compared with CLV3 (CZ) cell population (Fig. 4C). Similarly, the two vascular cell types, the xylem and phloem, clustered together; however, they showed a lower degree of relatedness to stem cell population of the CZ. Notably, the extent of this similarity was not always related to the spatial proximity of cell types within the SAMs. For example, the clonally distinct but adjacent L1-meristematic (HMG) and L2-meristematic (HDG4) cells showed a lesser degree of relatedness in their gene expression profiles. Similarly, cell populations of the root system and the SAM, including vascular cells, formed separate clusters (Fig. 4D). Clonal separation of SAM cell layers at early stages of embryonic development might have contributed to the greater diversity of gene expression profiles.

The estimation of the amount of expressed (detected) genes in the ten cell populations of the SAM presented in this study and in root samples (S4, J2501, xylem, S17, RM1000, JO121, SUC2, phloem, S32, cortex, COBL9, src5, J0571, wol, LRC, gl2) (Brady et al., 2007) (see Materials and methods for details) revealed 958 and 2187 genes that are expressed exclusively in the shoot and the root, respectively (Fig. 4A; see supplementary material Table S3). GO terms associated with flower development (GO:0009908), stomata development (GO:0010103) and chlorophyll metabolic processes (GO:0015994) were overrepresented in the genes expressed in shoot cell populations (Fig. 4B; see supplementary material Table S4). On the contrary, the gene set of the root cell-types contained GO terms associated with responses to toxins (GO:0009404) (Fig. 4B) and hormones (GO:0042446), which included genes involved in auxin biosynthesis (YUCCA3 and YUCCA5), cytokinin biosynthesis [ISOPENTENYL TRANSFERASE (IPT) 1, 3, 5, 7 and 8] (Zhao, 2010) and ethylene biosynthesis (1-AMINOCYCLOPROPANE-1-CARBOXYLATE SYNTHASE 7, 8 and 11) (Tsuchisaka et al., 2009). The IPT class genes were not expressed in SAM cell populations, which suggests that cytokinin precursor-N⁶-(Δ²-isopentenyl)adenine (iP) may be transported into meristematic cells where it is converted into active forms of cytokinin through the action of locally synthesized LONELYGUY class of proteins (Kurakawa et al., 2007; Kuroha et al., 2009).

To understand the orchestration of the gene regulatory networks (GRNs) mediated by transcription factors (TFs), we mined our data sets for the distribution expression patterns of TFs (Jiao and Meyerowitz, 2010). Of the 1270 TF genes considered for this analysis, 1225 were detected in SAM cell populations. Among these, only 296 TFs were CPEGs (see supplementary material Table S1). About 39% CPEGs encoding for TFs were found in clonally distinct cell layers of the SAMs, whereas only 22% CPEGs encoding for TFs were present in distinct cell populations located in the CZ and the PZ. Thirty-five percent of CPEGs encoding for TFs were present in vasculature cell populations. The greater diversity of TF representation among the distinct layers of SAMs may explain lesser similarity in the transcriptome among cell layers. AP2_EREBP, C2H2, MADS, MY and WRKY family TFs are highly expressed in the meristematic L1 layer; C2C2_DOF members are predominant in the L2 meristematic layer; and bHLH, GRAS, HB, NAC, SBP and ZF_HD are enriched in L3 layer (see supplementary material Fig. S6). Furthermore, members of the AP2_EBERP family, which are involved in mediating stress responses, and members of the WRKY family, which are involved in defense responses, are expressed in the L1 epidermal layer. The CPEGs of the CLV3/stem cell population contained genes encoding HOMEOBOX and MADS TFs, whereas the FIL/differentiating cell population contained genes that code for C2C2_YABBY, MYB and SBP TFs (see supplementary material Fig. S6 and Table S1). Last, the genes that encode AP2_EREBP, AUX_IAA and WKRY TFs were among xylem cell type CEPEGs, whereas genes that encode members of the bHLH, C2C2_DOF, C2H2, GRAS, MYB, NAC and SBP TF families were found in phloem CPEGs. The analysis of the spatial expression of TFs and identification of cis-motifs to which they bind may provide clues to the GRNs that underlie specification of different cell types. We identified the TF-binding motifs by studying the enrichment of upstream regulatory elements within the promoter regions of CPEGs identified in the L1, L2, L3 cell layers, and the CZ and the PZ domains (Fig. 5A,B; see supplementary material Table S7 and see Materials and Methods). A greater degree of diversity was observed among the TF-binding motifs present in the promoters of distinct cell layers CPEGs relative to the CZ and the PZ domain CPEGs, suggesting that a highly complex GRN specifies cell types in the SAM layers.

An earlier study has shown that 1638 genes respond to abscisic acid, indole-3 acetic acid, gibberellin, cytokinin, jasmonate, brassinosteroid and ethylene treatments (Nemhauser et al., 2006). Of the 1638 hormone-responsive genes, 1298 were detected in SAM cell populations, out of which 469 were identified as CPEGs. A majority of these genes were CPEGs in one of the three cell layers of SAMs (see supplementary material Table S1). Notably, stem cells or CZ CPEGs contained very few hormone-responsive genes when compared with the CPEGs in the PZ and the RM (Fig. 5D; see supplementary material Table S9). In addition, we observed ethylene-responsive genes in the L1-meristematic layer, the L3 layer and vascular cell population CPEGs (Fig. 5C; see supplementary material Table S1; Fig. 3A,B,D,E) which suggests that ethylene may be broadly perceived in SAMs.

The SAM is a complex developmental field consisting of cell types that form clonal layers within the CZ and the PZ that are closely related and clonal layers that are distantly related. Our analysis indicates cell types in clonal layers that are distantly related by birth show a lesser degree of relatedness in gene expression profiles than those that are closely related. A greater diversity of TF classes and cis-elements detected in promoters of genes that are differentially expressed in cell layers suggests that a highly complex GRN may underlie cell type specification in the SAM cell layers relative to the GRNs that specify cell types in the CZ and the PZ. Furthermore, small fractions of the total number of genes expressed (∼13% in the root system and ∼5.7% in the SAM) are organ specific, suggesting that the network topology of interacting components within a given cell type may play a crucial role in cell type specification.

Our work shows that CPEGs of the L1 cell populations are implicated in pathways related to biosynthesis and secretion of epicuticular wax. Genetic studies have implicated two functionally redundant HD-ZIP IV classes of TFs (AtML1 and PROTODERMAL FACTOR2) in epidermal cell differentiation (Lu et al., 1996; Abe et al., 2003). Upon differentiation, the L1 epidermal cells asymmetrically secrete a hydrophobic lipid layer consisting of cutin and waxes that forms a protective film (Jeffree, 2006). Furthermore, the AP2/EREBP and MYB TFs have been shown to regulate the L1 layer expression of genes involved in cuticle biosynthesis (Javelle et al., 2011). Analysis presented here may be used in future studies to decipher the mechanisms by which developmental TFs, such as AtML1 and PDF2, function in tandem with the AP2/EREBP and MYB family TFs in controlling epidermal cell type development and its function. Besides protective function, the L1 layer has been shown to play a role in stem cell maintenance (Knauer et al., 2013), regulation of cell division patterns in sub-epidermal layers (Reinhardt et al., 2003a,b; Savaldi-Goldstein et al., 2007), lateral organ primordia specification (Reinhardt et al., 2005, 2003a,b) and establishment of lateral organ polarity (Reinhardt et al., 2005). The L1 layer transcriptome presented here should facilitate functional analysis of these genes, leading to better molecular insights into the involvement of epidermal layer in plant development.

Until recently, it was widely believed that the SAM is heterotrophic and contains only proplastids that lack thylakoids and chlorophyll-binding proteins (Fleming, 2006). However, recent work has revealed that plastids in SAMs contain PSII-associated chlorophyll and thylakoid membranes (Charuvi et al., 2012). Our data show elevated expression of genes encoding photosystem I and II light-harvesting complexes, along with other accessory proteins of the photosynthetic machinery, in the L3 cell population; this provides molecular support for the presence of the functional chloroplast machinery in SAMs. Several developmental processes in SAMs, including leaf initiation and positioning, depend on light (Yoshida et al., 2011). Given that chlorophyll biosynthesis in higher plants is light dependent, the light-dependent developmental regulation may be intricately linked to the chloroplast maturation and metabolic status of cells of SAMs.

Earlier studies have shown that a higher level of auxin signaling is required for specification and growth of the lateral organ primordia (Reinhardt et al., 2003a,b). Furthermore, GA has been implicated in promoting cell differentiation, and members of the class I KNOX family have been shown to repress GA biosynthesis in the meristem center (Sakamoto et al., 2001; Bolduc and Hake, 2009). In agreement with earlier studies, the genome level analysis presented here shows poor representation of auxin and GA-responsive genes in the CZ. The cell population transcriptome presented here can be used to reconstruct cellular responsiveness to various plant hormones, as shown for stem cell promoting TF-WUS (Yadav et al., 2013), which may lead to new insights into the hormonal regulation of stem cell maintenance and organ differentiation.

PAtML1::mGFP5-ER was generated by fusing ∼3.3 kb promoter fragment to mGFP5-ER. A ∼2 kb upstream sequence including the first exon of HMG (locus At1g04880) was fused to H2B-YFP to generate PHMG::H2B-YFP. A 3 kb HDG4 promoter was fused with H2B-YFP to generate PHDG4::H2B-YFP. A ∼2 kb AtHB8 promoter region was fused with H2B-YFP to generate PAtHB8::H2B-YFP. PLAS::LAS-mGFP5 (Goldshmidt et al., 2008), PKAN1::KAN1-mGFP5 (Yadav et al., 2013), and PCLV3::mGFP-ER, PWUS::mGFP-ER and PFIL::ds-Red (Yadav et al., 2009) have been described previously. The oligonucleotides used to amplify promoter fragments of PHMG, PHDG4 and PAtHB8 are in supplementary material Table S10. All transgenic lines were generated in ap1-1;cal1-1 genetic background using the floral dip method (Clough and Bent, 1998).

Protoplasting, RNA extraction and ATH1 gene chip hybridization have been described previously (Yadav et al., 2009). All microarray data generated in this study are deposited in GEO under accession number GSE28109.

Microarray data analyses were performed in the statistical environment R using Bioconductor packages (Gentleman et al., 2005; R Development Core Team, 2011). The probe set to locus mappings for the ATH1 chip were obtained from TAIR (V-2010-12-20). Controls and probe sets matching no or several loci in the Arabidopsis genome were ignored in the downstream analysis steps. In addition, redundant probe sets that represent the same locus several times were counted only once. The quality of the Affymetrix GeneChips was assessed with analysis routines provided by the affyPLM and arrayQualityMetrics packages (Bolstad et al., 2005; Kauffmann et al., 2009). The normalization of the raw data Cel files was performed with the MAS5 algorithm using the default settings of the corresponding R function of the Affymetrix package (Gautier et al., 2004).

Sample-wise hierarchical clustering was used to evaluate global expression similarities among cell types. The distance matrix required for this step was generated from the distance-transformed Spearman correlation coefficients obtained from all pair-wise comparisons of the sample expression vectors using the average among their normalized replicates. The actual clustering was performed with the hclust function in R, by choosing complete linkage as cluster joining rule.

Analysis of differentially expressed genes was performed with the LIMMA package using the normalized expression values as input. The Benjamini and Hochberg method was selected to adjust P values (P) for multiple testing and controlling false discovery (FDR) rates (Benjamini and Hochberg, 1995). As confidence thresholds, we used an adjusted P≤0.05 and a minimum fold change of 1.5. A gene was considered to be a DEG if it met both threshold criteria in at least one of all possible sample comparisons (contrasts) among the ten cell types. Cell population expressed genes (CPEGs) were identified by filtering the DEG results for candidate genes with an N-fold higher expression in the target cell type compared to all of the chosen reference cell types, where N (fold change cutoff) was set to values of ≤0, ≤1.5, ≤2 and ≤4. At the same time, all fold changes needed to be supported by an adjusted P value of ≤0.05 from LIMMA (see columns N-R in supplementary material Table S1). A more relaxed version of the CPEG analysis allowed one outlier. For example, if a gene was higher expressed in cell type ‘A’ compared with at least three out of a total of four reference cell types (e.g. ‘BCDE’), then it was considered to be a relaxed CPEG of cell type ‘A’ (see columns S-U in supplementary material Table S1). To estimate the amount of expressed (detected) and unexpressed genes in the SAM and root samples (S4, J2501, xylem, S17, RM1000, JO121, SUC2, phloem, S32, cortex, COBL9, src5, J0571, wol, LRC, gl2) (Brady et al., 2007), the present call information of the non-parametric Wilcoxon signed rank test was computed. Genes that received present calls in all replicates of a given cell type were considered to be expressed (see supplementary material Tables S5 and S6) (McClintick and Edenberg, 2006).

Enrichments analysis of GO terms was performed as described in an earlier study (Horan et al., 2008) (see supplementary material Tables S2 and S12). The CPEG sets from each cell population was imported into Cytoscape 2.8.1 plug-in BINGO 2.4.4 (Maere et al., 2005; Smoot et al., 2011). The network of enriched categories was organized using first the force directed layout, and subsequently the hierarchical layout (see supplementary material Fig. S5).

The sequence regions 1000 bp upstream of each gene model were considered. For the motif finding, 471 cis-regulatory elements available in the Agris database were downloaded (Yilmaz et al., 2011) and mapped to the upstream sequences of each gene set while allowing one mismatch per ten nucleotides (see supplementary material Table S7). The hypergeometric distribution was used to compute P values for estimating the significance of this enrichment and the Bonferroni method to adjust the P values for multiple testing (Sinha et al., 2006).

FACS, microarray and imaging were carried out at Institute for Integrative Genome Biology and Center for Plant Cell Biology, UCR. We thank Yuval Eshed for sharing reporters and K. S. Sandhu for help with data analysis.