Abstract

The shoot apical meristem (SAM) acts as a reservoir for stem cells. The central zone (CZ) harbors stem cells. The stem cell progenitors differentiate in the adjacent peripheral zone and in the rib meristem located just beneath the CZ. The SAM is further divided into distinct clonal layers: the L1 epidermal, L2 sub-epidermal and L3 layers. Collectively, SAMs are complex structures that consist of cells of different clonal origins that are organized into functional domains. By employing fluorescence-activated cell sorting, we have generated gene expression profiles of ten cell populations that belong to different clonal layers as well as domains along the central and peripheral axis. Our work reveals that cells in distinct clonal layers exhibit greater diversity in gene expression and greater transcriptional complexity than clonally related cell types in the central and peripheral axis. Assessment of molecular functions and biological processes reveals that epidermal cells express genes involved in pathogen defense: the L2 layer cells express genes involved in DNA repair pathways and telomere maintenance, and the L3 layers express transcripts involved in ion balance and salt tolerance besides photosynthesis. Strikingly, the stem cell-enriched transcriptome comprises very few hormone-responsive transcripts. In addition to providing insights into the expression profiles of hundreds of transcripts, the data presented here will act as a resource for reverse genetic analysis and will be useful in deciphering molecular pathways involved in cell type specification and their functions.

INTRODUCTION

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.

Central and peripheral organization

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.

Identification of cell population elevated gene expressions without outliers

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).

Pathways with elevated gene expressions in various cell populations. (A) Layered organization of SAM. (B) A schematic representation of cell types along the central and peripheral dimension on a 3D reconstructed SAM template showing the expression of PLAS::LAS-GFP. CLV3 and FIL expression domains are marked green and yellow, respectively. (C) SAM top view showing PKAN1::KAN1-GFP expression in the outer edges of the PZ. (D) GO term enrichment analysis of DEGs enriched in individual cell layers and cell populations along the central and peripheral zones in SAMs. Gene expression patterns are depicted in the form of heat maps for genes associated with enriched GO terms.

Relaxed CPEGs

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.

Validation of microarray data

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).

The effect of the fold-change cutoffs on CPEG set sizes

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.

CPEGs in SAM cell layers

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.

L1 layer CPEGs

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).

L2 layer CPEGs

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).

L3 layers/corpus CPEGs

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 Ca2+-dependent protein with peroxidase activity; a vacuolar Ca2+-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.

KAN1/outer PZ

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).

Similarities among cell types based on expression profiles

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.

A comparison of genes expressed in shoot and root cell type. (A) Venn diagram showing the intersection of genes detected in shoot and root cell populations. (B) Summary of GO term enrichment results highlighting some of the most dominant classifications. (C) Sample-wise hierarchical clustering of shoot microarray replicates. (D) A comparison of shoot- and root-derived microarray replicates to assess their relative closeness in terms of their transcriptome profile.

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-N6-(Δ2-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).

Orchestration of transcriptional programs

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.

Identification of cis-element and hormone-responsive transcripts in SAM cell populations. (A,B) Enrichment of cis elements in cell layers (A) and in the central and peripheral zones (B). (C) Number of hormone-responsive CPEGs plotted as a percentage against the total number of CPEGs for each cell population. (D) Relative contributions of hormone-responsive CPEGs depicted as a percentage against the number of all hormone-responsive CPEGs.

Hormonal pathways

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.

DISCUSSION

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.

MATERIALS AND METHODS

DNA constructs and transgenic lines

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 and microarray hybridization

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.

Analysis of microarray data

Raw data analysis and quality assessment

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 clustering

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.

DEG analysis

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).

Identification of regulatory motifs

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).

Acknowledgements

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.

Footnotes

Competing interests

Authors declare that they have no competing interests.

Author contributions

R.K.Y. and G.V.R. conceived research and designed experiments. R.K.Y. and M.T. performed the experiments. M.X. contributed H2B:YFP pGreen plasmid. T.G. performed microarray analysis. R.K.Y., T.G. and G.V.R. analyzed the data and wrote the paper.

Funding

This work was supported by NSF 2010 grant (IOS-0820842) to G.V.R. R.K.Y. is Ramalingaswami fellow and recipient of Innovative Young Biotechnologist Award from the Department of Biotechnology, Govt. of India and acknowledges DBT for current funding.

(1996). Identification of a meristem L1 layer-specific gene in Arabidopsis that is expressed during embryonic pattern formation and defines a new class of homeobox genes. Plant Cell8, 2155-2168.doi:10.1105/tpc.8.12.2155

Other journals from The Company of Biologists

When D'Arcy Wentworth Thompson’s On Growth and Form was published 100 years ago, it raised the question of how biological forms arise during development and across evolution. In light of the advances in molecular and cellular biology since then, a succinct modern view of the question states: how do genes encode geometry? Our new special issue is packed with articles that use mathematical and physical approaches to gain insights into cell and tissue patterning, morphogenesis and dynamics, and that provide a physical framework to capture these processes operating across scales.

Read the Editorial by guest editors Thomas Lecuit and L. Mahadevan, as they provide a perspective on the influence of D'Arcy Thompson's work and an overview of the articles in this issue.

The Node currently features a number of posts to accompany the special issue. Authors of articles in the special issue talk about what On Growth and Form means to them. The creators of the commissioned cover of the special issue, Jen Ma and Matthew Spremulli, share how the cover was conceived, designed and completed. And the authors of a special issue Research Article on leaf epidermis topology reveal the extent of D’Arcy Thompson’s influence in current plant development research.

The resources pages of the Node and the BSDB have been combined and refurbished, with the new page designed to provide a range of useful links, including information on advocacy and outreach, new teaching resources for schools and databases for a wide range of species.