Spatiotemporal Developmental Trajectories in the Arabidopsis Root Revealed Using High-Throughput Single-Cell RNA Sequencing

•scRNA-seq of Arabidopsis root cells captures precise spatiotemporal information•Defining expression features for cell types identify new developmental regulators•Cluster arrangement reflects developmental time with a centrally localized niche•Intricate waves of gene expression finely resolve developmental trajectories High-throughput single-cell RNA sequencing (scRNA-seq) is becoming a cornerstone of developmental research, providing unprecedented power in understanding dynamic processes. Here, we present a high-resolution scRNA-seq expression atlas of the Arabidopsis root composed of thousands of independently profiled cells. This atlas provides detailed spatiotemporal information, identifying defining expression features for all major cell types, including the scarce cells of the quiescent center. These reveal key developmental regulators and downstream genes that translate cell fate into distinctive cell shapes and functions. Developmental trajectories derived from pseudotime analysis depict a finely resolved cascade of cell progressions from the niche through differentiation that are supported by mirroring expression waves of highly interconnected transcription factors. This study demonstrates the power of applying scRNA-seq to plants and provides an unparalleled spatiotemporal perspective of root cell differentiation. High-throughput single-cell RNA sequencing (scRNA-seq) is becoming a cornerstone of developmental research, providing unprecedented power in understanding dynamic processes. Here, we present a high-resolution scRNA-seq expression atlas of the Arabidopsis root composed of thousands of independently profiled cells. This atlas provides detailed spatiotemporal information, identifying defining expression features for all major cell types, including the scarce cells of the quiescent center. These reveal key developmental regulators and downstream genes that translate cell fate into distinctive cell shapes and functions. Developmental trajectories derived from pseudotime analysis depict a finely resolved cascade of cell progressions from the niche through differentiation that are supported by mirroring expression waves of highly interconnected transcription factors. This study demonstrates the power of applying scRNA-seq to plants and provides an unparalleled spatiotemporal perspective of root cell differentiation. In recent years, high-throughput single-cell transcriptomics has developed to a point of becoming a fundamental, widely used method in mammalian research (Potter, 2018Potter S.S. Single-cell RNA sequencing for the study of development, physiology and disease.Nat. Rev. Nephrol. 2018; 14: 479-492Crossref PubMed Scopus (220) Google Scholar). Thousands of cells can be profiled simultaneously and analyzed accurately, revealing unique insights into developmental progressions, transcriptional pathways, and the molecular heterogeneity of tissues. The increasingly high-throughput nature of single-cell RNA sequencing (scRNA-seq) has been facilitated by the development of droplet technology (Macosko et al., 2015Macosko E.Z. Basu A. Satija R. Nemesh J. Shekhar K. Goldman M. Tirosh I. Bialas A.R. Kamitaki N. Martersteck E.M. et al.Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets.Cell. 2015; 161: 1202-1214Abstract Full Text Full Text PDF PubMed Scopus (3700) Google Scholar, Klein et al., 2015Klein A.M. Mazutis L. Akartuna I. Tallapragada N. Veres A. Li V. Peshkin L. Weitz D.A. Kirschner M.W. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells.Cell. 2015; 161: 1187-1201Abstract Full Text Full Text PDF PubMed Scopus (1885) Google Scholar) and increased automation (Zheng et al., 2017Zheng G.X. Terry J.M. Belgrader P. Ryvkin P. Bent Z.W. Wilson R. Ziraldo S.B. Wheeler T.D. McDermott G.P. Zhu J. et al.Massively parallel digital transcriptional profiling of single cells.Nat. Commun. 2017; 8: 14049Crossref PubMed Scopus (2343) Google Scholar). In brief, a cell is encapsulated within an oil droplet and lysed, and its transcripts reverse transcribed on barcoded beads. Following library production and sequencing, transcripts from individual cells can be identified from the bead-derived barcode and individual transcripts accounted for using unique molecular identifiers (UMIs) (Prakadan et al., 2017Prakadan S.M. Shalek A.K. Weitz D.A. Scaling by shrinking: empowering single-cell ‘omics’ with microfluidic devices.Nat. Rev. Genet. 2017; 18: 345-361Crossref PubMed Scopus (209) Google Scholar). However, while commonly used in animal systems, additional technical demands such as the necessity to break down cell walls (with subsequent transcriptional effects), high osmotic pressure sensitivities, and high cell size variability present potential challenges when applying this technology to plants. The Arabidopsis root provides an ideal tissue for analyzing the promise of scRNA-seq. The transcriptomes of key cell types have been well profiled, and the root shows a strict spatiotemporal organization. Radially, the root is organized in concentric rings of endodermis, cortex, and epidermis that surround a central stele, comprising the pericycle, phloem, and xylem (Figure S1A). These cell types originate from a specialized stem cell niche in which initials, surrounding the quiescent center (QC), divide in a predictable manner, giving rise to long cell files that capture their developmental trajectory along the length of the root (Figure S1B). Several gene expression atlases of the Arabidopsis root have been produced (Birnbaum et al., 2003Birnbaum K. Shasha D.E. Wang J.Y. Jung J.W. Lambert G.M. Galbraith D.W. Benfey P.N. A gene expression map of the Arabidopsis root.Science. 2003; 302: 1956-1960Crossref PubMed Scopus (985) Google Scholar, Brady et al., 2007aBrady S.M. Orlando D.A. Lee J.-Y. Wang J.Y. Koch J. Dinneny J.R. Mace D. Ohler U. Benfey P.N. A high-resolution root spatiotemporal map reveals dominant expression patterns.Science. 2007; 318: 801-806Crossref PubMed Scopus (860) Google Scholar, Li et al., 2016Li S. Yamada M. Han X. Ohler U. Benfey P.N. High-resolution expression map of the Arabidopsis root reveals alternative splicing and lincRNA regulation.Dev. Cell. 2016; 39: 508-522Abstract Full Text Full Text PDF PubMed Scopus (156) Google Scholar). These, however, have focused primarily on describing either radial or temporal expression profiles and typically relied on reporter lines to assess select cell types. scRNA-seq, on the other hand, allows the simultaneous, unbiased sampling of every type of cell at every developmental stage in one experiment. Here, we present a high-resolution scRNA-seq expression atlas of the Arabidopsis root that captures its precise spatiotemporal information, revealing key regulators and defining features for all major cell types. We show how QC cells and meristematic cells are distinguished and resolve intricate developmental trajectories that cells undergo during their transition from stem cell through differentiation. The precise waves of gene expression characterizing this process are mirrored by similar expression changes of highly interconnected transcription factors (TFs). Our atlas offers an unparalleled spatiotemporal perspective of root cell-type differentiation at a resolution not previously achievable. 4,727 Arabidopsis root cells from two biological replicates were isolated and profiled using droplet-based scRNA-seq. At ∼87,000 reads per cell, the median number of genes and transcripts detected per cell was 4,276 and 14,758, respectively (Figure S1C; Table S1). In total, transcripts for 16,975 genes were detected (RPM ≥ 1), which, after correction for read depth, represents ∼90% of genes detected by bulk RNA sequencing (RNA-seq) of protoplasted root tissue. Further, the global gene expression profiles of pooled scRNA-seq and bulk RNA-seq are highly correlated (r = 0.9; Figure S1D), indicating that plant scRNA-seq is highly sensitive. This methodology is also highly reproducible, as demonstrated by the facts that ∼96% of genes expressed (RPM ≥ 1) in one scRNA-seq replicate are detectable in the second and that expression across the two replicates is highly correlated (r = 0.99; Figure S1E). To identify distinct cell populations based on gene expression profiles, an unbiased, graph-based clustering was performed on the 4,727 single-cell transcriptomes using the Seurat software package (Satija et al., 2015Satija R. Farrell J.A. Gennert D. Schier A.F. Regev A. Spatial reconstruction of single-cell gene expression data.Nat. Biotechnol. 2015; 33: 495-502Crossref PubMed Scopus (2101) Google Scholar, Butler et al., 2018Butler A. Hoffman P. Smibert P. Papalexi E. Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species.Nat. Biotechnol. 2018; 36: 411-420Crossref PubMed Scopus (4241) Google Scholar) (Figure 1A). Genes induced by protoplasting (≥2-fold; q < 0.05) were identified by standard RNA-seq and dismissed prior to analysis (Figure S1F; Table S1). 15 distinct clusters were identified, each containing between 81 and 596 cells. These clusters harbored similar numbers of cells from each replicate, and their gene expression profiles were highly correlated across the replicates (r between 0.95 and 1; Table S2), highlighting again the impressive reproducibility of this technique. In order to attribute cell identities to these clusters, expression of cell-type-specific marker genes, either well established or identified from a curated collection of root transcriptomic datasets (Table S2; Efroni et al., 2015Efroni I. Ip P.L. Nawy T. Mello A. Birnbaum K.D. Quantification of cell identity from single-cell gene expression profiles.Genome Biol. 2015; 16: 9Crossref PubMed Scopus (98) Google Scholar), was compared across clusters. This allowed cell identities to be confidently assigned to 8 of the 15 clusters in the cluster cloud (Figure 1B; Table S2). Expression of key root development genes among these markers, such as PLT1, SCR, SHR, APL, COBL9, and GL2, shows high specificity to particular clusters (Figure S2). Cluster identities were confirmed with a complementary approach, whereby transcription profiles of differentially expressed (DE) genes governing the clusters were harvested from microarray datasets (Brady et al., 2007aBrady S.M. Orlando D.A. Lee J.-Y. Wang J.Y. Koch J. Dinneny J.R. Mace D. Ohler U. Benfey P.N. A high-resolution root spatiotemporal map reveals dominant expression patterns.Science. 2007; 318: 801-806Crossref PubMed Scopus (860) Google Scholar) and analyzed for tissue specificity (Figure S1G). Together, these approaches revealed that, with the exception of lateral root cap cells (for which limited marker data are available), all known major tissue types in the root were captured and are represented by identifiable clusters. Clusters 9 and 13 comprise cortex and endodermal cells, respectively (Figures 1B and S1G; Table S2). The identity of the endodermal cluster was further validated by the localized accumulation of GFP transcripts in one of the replicates generated from the pMIR166A:erGFP reporter line (see STAR Methods; Figures 2A and S3A). In addition, when cells were re-clustered incorporating scRNA-seq data of shortroot mutants (shr-3), which lack a defined endodermis (Helariutta et al., 2000Helariutta Y. Fukaki H. Wysocka-Diller J. Nakajima K. Jung J. Sena G. Hauser M.T. Benfey P.N. The SHORT-ROOT gene controls radial patterning of the Arabidopsis root through radial signaling.Cell. 2000; 101: 555-567Abstract Full Text Full Text PDF PubMed Scopus (809) Google Scholar; Figure S3A), otherwise well-dispersed shr-3 cells were absent from a cluster comprising endodermal cells of both wild-type replicates (Figures 2B and S3). This cluster analysis also shows that shr cells, while present in all other clusters, localize on the outskirts of some (Figure 2B). This points to subtler effects of SHR on cell types other than the endodermis; although some of this phenomenon may also be attributable to the fact that shr-3 is in the Ler background. Irrespective, this observation nicely highlights the potential of applying scRNA-seq to identify hidden phenotypic changes, whether stemming from natural variation or mutations. Clusters 10 and 3 comprise trichoblast and atrichoblast cells, respectively (Figures 1B and S1G; Table S2). Cluster 5 also contains trichoblast cells (Figure S1G). Although cells in this cluster show low expression of a number of atrichoblast marker genes, crucially, the trichoblast marker COBL9 is expressed in this cluster, whereas the atrichoblast marker, GL2, is not (Figures 1B and S2). The co-expression of atrichoblast marker genes hints at a degree of commonality between this subset of trichoblasts and its epidermal counterparts, perhaps reflecting a distinction in developmental stage to the trichoblast cells contained in cluster 10. Cluster 4 comprises stele cells while a neighboring cluster (12) comprises maturing xylem cells (Figures 1B and S1G; Table S2). Consistent with the tissue complexity of the stele, subclustering reveals cell heterogeneity within cluster 4. Particularly, phloem and pericycle cells are separated into two discrete subclusters (Figure S4), as highlighted by the highly subcluster-specific expression of genes such as APL (4.2), LBD29, and TIP2-3 (4.1) (Figure 2D; Bonke et al., 2003Bonke M. Thitamadee S. Mähönen A.P. Hauser M.T. Helariutta Y. APL regulates vascular tissue identity in Arabidopsis.Nature. 2003; 426: 181-186Crossref PubMed Scopus (352) Google Scholar, Porco et al., 2016Porco S. Larrieu A. Du Y. Gaudinier A. Goh T. Swarup K. Swarup R. Kuempers B. Bishopp A. Lavenus J. et al.Lateral root emergence in Arabidopsis is dependent on transcription factor LBD29 regulation of auxin influx carrier LAX3.Development. 2016; 143: 3340-3349Crossref PubMed Scopus (89) Google Scholar, Gattolin et al., 2009Gattolin S. Sorieul M. Hunter P.R. Khonsari R.H. Frigerio L. In vivo imaging of the tonoplast intrinsic protein family in Arabidopsis roots.BMC Plant Biol. 2009; 9: 133Crossref PubMed Scopus (74) Google Scholar). Finally, cluster 11 comprises both columella and QC cells (Figures 1B and S1G; Table S2), which can be separated into two subclusters. Subcluster 11.2 contains columella cells that express marker genes such as COBL2, NCED2, and ATL63 (Figure 2E; Brady et al., 2007bBrady S. Song S. Dhugga K. Rafalski A. Benfey P. Combining expression and comparative evolutionary analysis. The COBRA gene family.Plant Phys. 2007; 143: 172-187Crossref PubMed Scopus (105) Google Scholar, Efroni et al., 2015Efroni I. Ip P.L. Nawy T. Mello A. Birnbaum K.D. Quantification of cell identity from single-cell gene expression profiles.Genome Biol. 2015; 16: 9Crossref PubMed Scopus (98) Google Scholar). In contrast, transcripts for the QC-expressed genes AGL42, BBM, and TEL1 are largely limited to cells in subcluster 11.1 (Figure 2E; Nawy et al., 2005Nawy T. Lee J.Y. Colinas J. Wang J.Y. Thongrod S.C. Malamy J.E. Birnbaum K. Benfey P.N. Transcriptional profile of the Arabidopsis root quiescent center.Plant Cell. 2005; 17: 1908-1925Crossref PubMed Scopus (260) Google Scholar, Efroni et al., 2015Efroni I. Ip P.L. Nawy T. Mello A. Birnbaum K.D. Quantification of cell identity from single-cell gene expression profiles.Genome Biol. 2015; 16: 9Crossref PubMed Scopus (98) Google Scholar). Given the small number of QC cells per root, this cluster may well contain other transcriptionally similar cells, perhaps the adjacent initials in the niche. However, importantly, the fact that QC cells are captured illustrates well the possibilities of this methodology for studying rare cell types or elucidating transcriptional subtleties affecting small numbers of cells within a tissue. The identity of cells in the remaining clusters is less obvious. Overall gene expression in cells within clusters 0, 1, and 14 is comparatively low (Figure S5A), likely masking their identity at this level of sequencing resolution. However, expression values extracted from a longitudinal microarray dataset (Brady et al., 2007aBrady S.M. Orlando D.A. Lee J.-Y. Wang J.Y. Koch J. Dinneny J.R. Mace D. Ohler U. Benfey P.N. A high-resolution root spatiotemporal map reveals dominant expression patterns.Science. 2007; 318: 801-806Crossref PubMed Scopus (860) Google Scholar) for the top DE genes defining these clusters suggest that they comprise mature cells of mixed identity (Figure S5B). In contrast, cells in the final four clusters (2, 6, 7, and 8) show markedly meristematic-based expression profiles (Figure S5B). Notable histone and cytokinesis-linked genes, such as KNOLLE, ENODL14, and ENODL15, are among the most prominently DE genes for these clusters (Figure S2; Table S2; Lauber et al., 1997Lauber M.H. Waizenegger I. Steinmann T. Schwarz H. Mayer U. Hwang I. Lukowitz W. Jürgens G. The Arabidopsis KNOLLE protein is a cytokinesis-specific syntaxin.J. Cell Biol. 1997; 139: 1485-1493Crossref PubMed Scopus (400) Google Scholar, Adrian et al., 2015Adrian J. Chang J. Ballenger C.E. Bargmann B.O. Alassimone J. Davies K.A. Lau O.S. Matos J.L. Hachez C. Lanctot A. et al.Transcriptome dynamics of the stomatal lineage: birth, amplification, and termination of a self-renewing population.Dev. Cell. 2015; 33: 107-118Abstract Full Text Full Text PDF PubMed Scopus (94) Google Scholar). Subclustering revealed some cell-type identities, albeit that they are generally less distinct than those of the clusters described above. For example, subcluster 2.4 shows a distinct cortex identity (Figure S4). Curiously, this subcluster is positioned adjacent to the main cortex cell cluster. A comparable pattern is seen in clusters 7 and 8 with trichoblast cell identity apparent in those subclusters (7.3 and 8.3) closest to the adjacent defined trichoblast clusters (Figure S4). It is interesting to note that when comprehending all the clusters together, the meristematic clusters are closely localized in the center of the cluster cloud with the subcluster containing QC cells (11.1) at the heart of this. Meanwhile, those clusters with distinct, mature cell identities span out from the meristematic clusters (Figure 1A; Video S1), suggesting an overall cluster arrangement that reflects developmental time. Subclustering of the meristematic clusters refines this idea, showing a degree of closeness of mature and developing cells of the same eventual fate. This notion is further supported by pseudotime analysis across all cells, which reveals that genes DE in cells of the central clusters describe the beginning of cell fate progressions (Figure 2C). eyJraWQiOiI4ZjUxYWNhY2IzYjhiNjNlNzFlYmIzYWFmYTU5NmZmYyIsImFsZyI6IlJTMjU2In0.eyJzdWIiOiI5NjNlY2RlOGY1NzEyZGEwN2RmNzViYWNkMmExN2FkYyIsImtpZCI6IjhmNTFhY2FjYjNiOGI2M2U3MWViYjNhYWZhNTk2ZmZjIiwiZXhwIjoxNjc5MjAyMDgwfQ.lXQ2SjISSrYbOdXmDwGdypKwzRuF2jH9itgPajZjM2wTolRnf3qQK3SS72k9gzFtEIxo2Jb2d_fuRSSdKaIddE75dHvYKawNCua3fP4q0Db0Q49dA5RTBGZAuvtBdlNq5KQdePE2P-9_tXBbNKBu6cVXMbXS56TxhoG2xBF13KGYx2tyj0R_yC_CWHFnlthmHYauu0-Du-6k31attdaK8xqU4Y64CMd0OMTQ8vgq93uO01RRFZdawTAxlG-oGu3XneTP_VyKk0GfQPqH4jECUS8RRLrP5KObYg-yK_EnB4-79SMtKb8PzVBC84Sh-Yc3WPabp2EWnsbGVCPnobAW-w Download .mp4 (15.93 MB) Help with .mp4 files Video S1. scRNA-Seq of the Arabidopsis Root Reveals Distinct Clusters, Related to Figure 1 Likewise, the cluster cloud reveals an organization that captures the “lineage” relationships between cell and tissue types. For instance, the trichoblast and atrichoblast clusters, as well as the xylem, vasculature, and the cortex and endodermis clusters, are positioned next to each other within the cluster cloud. The position of the columella in a cluster with QC cells indicates a higher degree of transcriptional accord between these cell types than between these cell types individually and others. This is reflected in the fact that key developmental regulators, such as PLT2, PLT3, and PIN4 are co-expressed in the columella and QC (Galinha et al., 2007Galinha C. Hofhuis H. Luijten M. Willemsen V. Blilou I. Heidstra R. Scheres B. PLETHORA proteins as dose-dependent master regulators of Arabidopsis root development.Nature. 2007; 449: 1053-1057Crossref PubMed Scopus (594) Google Scholar, Feraru and Friml, 2008Feraru E. Friml J. PIN polar targeting.Plant Physiol. 2008; 147: 1553-1559Crossref PubMed Scopus (108) Google Scholar). This way of contemplating clusters, along with pseudotime visualization, thus offers a valuable director for early comprehension of developmental trajectories, particularly in the absence of a priori knowledge, such as a reference atlas. Given that detailed reference datasets are available only for select tissue and organ types in very few plant species, we developed an unbiased approach to assign cell type identities to scRNA-seq-generated cell clusters. Genes DE in a given cluster compared to all other clusters (q < 0.01; average log fold change [FC] ≥ 0.25) were identified using “biomod” on Seurat (McDavid et al., 2013McDavid A. Finak G. Chattopadyay P.K. Dominguez M. Lamoreaux L. Ma S.S. Roederer M. Gottardo R. Data exploration, quality control and testing in single-cell qPCR-based gene expression experiments.Bioinformatics. 2013; 29: 461-467Crossref PubMed Scopus (217) Google Scholar). DE genes were further narrowed down by applying the criteria that cluster-specific marker genes must be expressed in ≥10% of cells within the cluster (PCT1), and ≤10% of cells across all other clusters (PCT2). Applying these criteria, we uncovered expected marker genes alongside hundreds of additional genes diagnostic for a given developmental stage or cell type that encompass every cluster (Table S3). The top two cluster-specific genes (based on average log FC) for each cluster are expressed across a substantial proportion of cells specifically within one cluster, with the exception of genes for clusters 0 and 1, which show substantial co-expression in cluster 14 (Figure 3A). In addition, Gene Ontology (GO) overrepresentation analysis on cluster-specific gene sets reveal GO terms appropriate to their biology (Table S3). For example, the meristematic clusters 2 and 8 show an abundance of marker genes implicated in processes related to cell proliferation and DNA replication, respectively. Further, markers for the root-hair-cell cluster 10 are enriched in trichoblast differentiation and maturation terms; for the QC- and columella-containing cluster 11, in root development and starch biosynthesis; and for cluster 12, in xylem development and secondary cell wall biogenesis. Finally, genes required for the formation and suberization of the Casparian strip are among the markers for cluster 13, which comprises endodermal cells. However, a notable outcome of this analysis is the number of marker genes for which a root function has yet to be assigned. This illustrates the potential of scRNA-seq for identifying new developmental regulators. To further validate this strategy for marker gene calling and for assigning cell identity to clusters without other references, we assessed the spatiotemporal patterns of expression for select genes using transcriptional promoter:3xVenus-NLS reporter lines. Prioritizing by a balance of high log-fold change, high PCT1, low PCT2, and a lack of prior biological information relating to cell-type specificity and root development, we selected ten genes from across clusters. Expression for eight of the ten genes tested localized to specific cell types and/or root zones in line with predictions. Specific expression in the cortex (AT1G62510) and maturing trichoblasts (MES15) was observed for marker genes for clusters 9 and 10, respectively (Figures 3B and 3C), while genes selected from cluster 4 revealed highly specific phloem (PME32) and pericycle (ATL75) expression (Figures 3D and 3E). MLP34 is expressed in the atrichoblasts, as expected for a marker for cluster 3 (Figure 3F). However, expression is also seen in cells of the lateral root cap (Figure 3F), a cell type to which a cluster could not be assigned. MLP34 shows expression in some cells in cluster 1 (Figure S2), indicating that this cluster may in fact contain cells of the lateral root cap, although further analysis is needed to confirm this. Finally, expression of genes selected from the meristematic cluster 2 was found to localize to the meristematic cortex and endodermis (AT3G22120), the meristematic cortex (AT1G62500), or the meristematic vasculature (PIP2-8) (Figures 3G–3I). Given the common occurrence of cis-regulatory motifs in the introns of genes, the fact that promoter fusions for eight out of the ten marker genes tested confirm predictions is notable. This unbiased approach for assigning identities to cell clusters could prove invaluable when no reference data are available. Moreover, our results reveal a level of sensitivity beyond that of assigning whole cluster identity. This is typified by PME32 and ATL75 whose promoter fusions show expression in the phloem and pericycle, respectively, in accordance with their expression being predominant in cells of subclusters 4.2 and 4.1, respectively (Figure S2). Furthermore, the meristematic vasculature marker, PIP2-8, is primarily expressed in cells of subcluster 2.3, positioned adjacent to the vasculature cluster, while the meristematic cortex and endodermis marker, AT3G22120, is mainly expressed in cells of subcluster 2.2 positioned near these mature cell types (Figure S2). The expression patterns observed for the latter reporters thus further reinforce the hypothesis that meristematic subclusters share expression features with the mature cell types they are closest to. Moreover, a gene’s expression profile across the cluster cloud is a confident predictor of its localized expression in planta. The fact that QC cells are captured offers a unique opportunity to study this rare cell type, a point quite pertinent given that RNA-seq analysis shows the established marker WOX5 to be induced upon protoplasting (Table S1). Within subcluster 11.1, 36 cells express at least half of 15 proposed QC genes (Figure 4A; Table S2; Efroni et al., 2015Efroni I. Ip P.L. Nawy T. Mello A. Birnbaum K.D. Quantification of cell identity from single-cell gene expression profiles.Genome Biol. 2015; 16: 9Crossref PubMed Scopus (98) Google Scholar, Nawy et al., 2005Nawy T. Lee J.Y. Colinas J. Wang J.Y. Thongrod S.C. Malamy J.E. Birnbaum K. Benfey P.N. Transcriptional profile of the Arabidopsis root quiescent center.Plant Cell. 2005; 17: 1908-1925Crossref PubMed Scopus (260) Google Scholar), which is in line with the sampling depth of scRNA-seq and the relatively low expression of most QC genes. The high number of QC cells captured likely reflects a bias in the methodology toward capturing small cells (see STAR Methods), which may also account for an overrepresentation of meristematic cells. Reinforcing our QC cell calling, it is of note that genes marking initial cells directly neighboring the QC, such as AT3G22120, AT1G62500, and PIP2-8 (Figures 3G–3I), show no (33 cells), or negligible (3 cells) expression in the QC cells. Additionally, cells expressing such genes cluster away from the QC, in localized regions of meristematic cluster 2, adjacent to their mature-cell counterparts (Figure S2). Transcriptomic comparison between the QC cells and undifferentiated cells of the meristem (cells in clusters 2, 6, 7, and 8), identified 254 genes preferentially expressed in the QC (Table S3). While meristematic cells are distinguished by expression of genes involved in cell division and DNA replication, cells of the QC are not. Instead, transcription is an enriched GO term, as is auxin biosynthesis, which is fitting given the role of auxin in QC specification (Sabatini et al., 1999Sabatini S. Beis D. Wolkenfelt H. Murfett J. Guilfoyle T. Malamy J. Benfey P. Leyser O. Bechtold N. Weisbeek P. et al.An auxin-dependent distal organizer of pattern and polarity in the Arabidopsis root.Cell. 1999; 99: 463-472Abstract Full Text Full Text PDF PubMed Scopus (1051) Google Scholar, Galinha et al., 2007Galinha C. Hofhuis H. Luijten M. Willemsen V. Blilou I. Heidstra R. Scheres B. PLETHORA proteins as dose-dependent master regulators of Arabidopsis root development.Nature. 2007; 449: 1053-1057Crossref PubMed Scopus (594) Google Scholar). Further, unexpectedly, genes with functions in glucosinolate biogenesis and callose deposition are overrepresented among those genes DE in the QC (Table S3). This finding, in particular, is intriguing. Both processes are characteristic of a defense response, which seems curious given the QC’s internal location, insulated from external stimuli. Their prominence instead points toward a biology of QC cells not previously appreciated. The recent finding that 3-hydroxypropylglucosinolate acts as a reversible inhibitor of root growth (Malinovsky et al., 2017Malinovsky F.G. Thomsen M.F. Nintemann S.J. Jagd L.M. Bourgine B. Burow M. Kliebenstein D.J. An evolutionarily young defense metabolite influences the root growth of plants via the ancient TOR signaling pathway.Elife. 2017; 6: e29353Crossref PubMed Scopus (58) Google Scholar) is, in this regard, intriguing. Likewise, that small RNAs are prevented from moving in and out of the QC (Skopelitis et al., 2018Skopelitis D.S. Hill K. Klesen S. Marco C.F. von Born P. Chitwood D.H. Timmermans M.C.P. Gating of miRNA movement at defined cell-cell interfaces governs their impact as positional signals.Nat. Commun. 2018; 9: 3107Crossref PubMed Scopus (59) Google Scholar) points to a unique regulation of cell-cell communication via plasmodesmata in the QC. Among the genes DE between the QC and meristematic cells of the root, 47 show a particularly strong expression bias to the QC cells (log FC ≥ 0.25; PCT1 ≥ 10; PCT2 ≤ 10