
==== Front
Plant Commun
Plant Commun
Plant Communications
2590-3462
Elsevier

S2590-3462(24)00286-4
10.1016/j.xplc.2024.100978
100978
Research Article
Single-cell network analysis reveals gene expression programs for Arabidopsis root development and metabolism
Han Ershang 1
Geng Zhenxing 1
Qin Yue 1
Wang Yuewei 1
Ma Shisong sma@ustc.edu.cn
12∗
1 MOE Key Laboratory for Cellular Dynamics, School of Life Sciences, Division of Life Sciences and Medicine, University of Science and Technology of China, Innovation Academy for Seed Design, Chinese Academy of Sciences, Hefei 230027, China
2 School of Data Science, University of Science and Technology of China, Hefei 230027, China
∗ Corresponding author sma@ustc.edu.cn
22 5 2024
12 8 2024
22 5 2024
5 8 1009783 1 2024
24 3 2024
20 5 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Single-cell RNA-sequencing datasets of Arabidopsis roots have been generated, but related comprehensive gene co-expression network analyses are lacking. We conducted a single-cell gene co-expression network analysis with publicly available scRNA-seq datasets of Arabidopsis roots using a SingleCellGGM algorithm. The analysis identified 149 gene co-expression modules, which we considered to be gene expression programs (GEPs). By examining their spatiotemporal expression, we identified GEPs specifically expressed in major root cell types along their developmental trajectories. These GEPs define gene programs regulating root cell development at different stages and are enriched with relevant developmental regulators. As examples, a GEP specific for the quiescent center (QC) contains 20 genes regulating QC and stem cell niche homeostasis, and four GEPs are expressed in sieve elements (SEs) from early to late developmental stages, with the early-stage GEP containing 17 known SE developmental regulators. We also identified GEPs for metabolic pathways with cell-type-specific expression, suggesting the existence of cell-type-specific metabolism in roots. Using the GEPs, we discovered and verified a columella-specific gene, NRL27, as a regulator of the auxin-related root gravitropism response. Our analysis thus systematically reveals GEPs that regulate Arabidopsis root development and metabolism and provides ample resources for root biology studies.

This study reports a single-cell gene co-expression network analysis performed with single-cell transcriptome datasets from Arabidopsis roots, leading to the identification of 149 gene expression programs (GEPs) that regulate the development of major root cell types and participate in cell-type-specific metabolism. Using these GEPs, the authors identified the columella-specific gene NRL27 and confirmed its role as a regulator of the auxin-related root gravitropism response.

Key words

single-cell RNA sequencing
single-cell gene co-expression network analysis
gene expression program
Arabidopsis root development
NRL27
gravitropism
Published: May 22, 2024
==== Body
pmcIntroduction

Plant roots function in water and nutrient absorption and provide mechanical support for the shoots. Root development is critical for forming the plant body and adapting to the environment. Biologists have long sought to systematically elucidate the gene networks that control root development. In Arabidopsis, roots display two developmental axes (Dolan et al., 1993). Transversely, cell types are arranged in circular layers surrounding a central vasculature, thus forming a radial developmental axis. Longitudinally, cell lineages are positioned along a temporal developmental axis, in which the youngest cells are located near the root-tip stem cell niche and the oldest cells are closest to the shoot. The simple cell-type and structural composition of Arabidopsis roots has made them an ideal model for studies of root systems biology.

Transcriptomic profiling is essential in root research. Combining cell-type-specific reporter lines with cell sorting, previous studies have used bulk transcriptomic approaches to characterize gene expression patterns of 14 root cell types in Arabidopsis (Brady et al., 2007; Li et al., 2016). These studies also revealed developmental-stage-specific expression patterns by separately analyzing the meristem, elongation, and maturation zones of primary roots. However, these assays analyzed gene transcription using mixtures of cells at the tissue level, which masked the molecular events that regulate development and occur at the cellular level. Advances in single-cell RNA-sequencing (scRNA-seq) technology in the past decade have provided unprecedented opportunities to characterize gene expression at the single-cell level. A number of studies have reported scRNA-seq datasets for Arabidopsis root tips (Denyer et al., 2019; Jean-Baptiste et al., 2019; Ryu et al., 2019; Shulse et al., 2019; Zhang et al., 2019; Wendrich et al., 2020; Shahan et al., 2022; Nolan et al., 2023). Focused analyses have also been performed on lateral roots and root phloem cells (Gala et al., 2021; Serrano-Ron et al., 2021; Otero et al., 2022). These studies identified the major cell types in Arabidopsis roots and inferred their developmental trajectories. Cell-type identification and trajectory inference can be categorized as cell-level analyses of single-cell datasets. The findings of cell-level analyses have greatly enhanced our understanding of root cellular heterogeneity.

In addition to cell-level analysis, gene-level analysis can be conducted on scRNA-seq datasets. Gene-level analysis aims to identify the genes or gene modules that drive cell development or regulate cell functions. However, scRNA-seq datasets are usually sparse and have many dropout expression values, posing great challenges for gene-level analysis (Lahnemann et al., 2020). For example, gene co-expression network analysis has been successfully applied to bulk plant transcriptome datasets to identify gene modules that function in various biological processes (Mentzen and Wurtele, 2008; Zhang et al., 2022); however, similar analysis of single-cell datasets is difficult because co-expression signals are extremely weak in single-cell transcriptomes owing to excessive dropouts (Crow and Gillis, 2018). Multiple methods have been developed to address the dropout problem, including grouping cells to form metacells or imputing missing values via statistical models (Li and Li, 2018; Baran et al., 2019). We recently developed a single-cell gene co-expression network analysis algorithm based on the graphical Gaussian model, named SingleCellGGM, which was modified from our previous method for graphical Gaussian model network analysis of bulk transcriptomes (Ma et al., 2007; Xu et al., 2023). SingleCellGGM uses an approach based on random gene sampling to calculate partial correlation coefficients (pcors) between gene pairs to evaluate their co-expression strength. The pcor is the correlation between two genes after removing the effects of other genes and is considered a better choice for gene network analysis than the commonly used Pearson correlation coefficient (Wille et al., 2004; Schafer and Strimmer, 2005). When applied to two mouse single-cell datasets (Han et al., 2018; La Manno et al., 2021), SingleCellGGM achieved satisfactory network results despite these two datasets being highly sparse (Xu et al., 2023).

In the present study, we performed a single-cell gene co-expression network analysis on publicly available Arabidopsis root scRNA-seq datasets via SingleCellGGM. The analysis identified 149 gene co-expression modules, which we considered to be gene expression programs (GEPs). We identified a series of GEPs with specific expression in all major root cell types at different developmental stages. These GEPs are enriched with developmental regulators and provide numerous candidate genes for future functional studies. As an example, GEP M20 is specifically expressed in early sieve elements (SEs) and contains 17 known SE developmental regulators out of 116 genes; other GEPs for middle and late SE development were also identified. We also identified GEPs for metabolic pathways with cell-type-specific expression, indicating the existence of cell-type-specific metabolism in Arabidopsis roots. Using the GEPs, we identified the columella-specific gene NRL27 and verified its role as an auxin-related regulator of root gravitropism response. Our analysis thus provides ample resources for the study of Arabidopsis root development and metabolism.

Results

Single-cell gene co-expression network analysis of Arabidopsis roots

Although scRNA-seq transcriptomes of Arabidopsis roots have been reported previously, no comprehensive gene co-expression network analyses have been conducted. To explore the discovery potential of the available datasets, we performed a single-cell gene co-expression network analysis to identify GEPs regulating root development and metabolism (Figure 1A). Given that co-expression network analysis typically benefits from a larger number of transcriptomes, we collected three publicly available scRNA-seq datasets of Arabidopsis roots from Jean-Baptiste et al., Wendrich et al., and Zhang et al. for use in the analysis (Jean-Baptiste et al., 2019; Zhang et al., 2019; Wendrich et al., 2020). We integrated these three datasets using Harmony and annotated their cells using a root scRNA-seq reference dataset from Shahan et al. via cell-label transfer (Korsunsky et al., 2019; Stuart et al., 2019; Shahan et al., 2022). In addition to cells profiled in their own study, the Shahan dataset contains root cells derived from two other studies by Ryu et al. and Denyer et al. (Denyer et al., 2019; Ryu et al., 2019). As in previous studies, cells belonging to different root cell types were identified in our integrated analysis, including quiescent center (QC), columella, lateral root cap (LRC), atrichoblast (nonhair cell), trichoblast (hair cell), cortex, endodermis, procambium, phloem, xylem, and pericycle cells, which were visualized in a uniform manifold approximation and projection (UMAP) plot (Figure 1B). Cell-type annotations were supported by the expression of known marker genes (Figure 1C). Developmental stages of the cells, corresponding to the meristem, elongation, and maturation zones of the primary roots or distal/proximal columella/LRC, were also identified; the meristem-stage cells were located in the center of the UMAP plot, whereas those from the elongation and maturation stages or proximal and distal columella/LRC spread outward consecutively (Figure 1D).Figure 1 Single-cell co-expression network analysis of Arabidopsis roots reveals gene expression programs (GEPs).

(A) The analysis workflow.

(B) UMAP plot of root cells after integration of three root scRNA-seq datasets via Harmony. Cells are colored according to their cell-type annotations. QC, quiescent center; LRC, lateral root cap; CC, companion cell; PPP, phloem pole pericycle; XPP, xylem pole pericycle.

(C) Expression patterns of known root-cell marker genes.

(D) UMAP plot of the root cells colored according to their developmental stages.

(E) The AtRootGGM gene co-expression network. Dots represent genes, which are colored according to their GEP IDs, and connections between genes indicate that they have co-expression patterns. Only the top 15 largest GEPs are shown owing to space limitations.

(F) A subnetwork for GEP M82 extracted from AtRootGGM. Genes highlighted in red have the GO term “cell fate commitment.”

(G) A subnetwork extracted from AtRootGGM2 for genes of GEP M82 of AtRootGGM. AtRootGGM2 was constructed using the Shahan et al. (2022) dataset. Of the 38 genes in GEP M82, 30 were found in AtRootGGM2, but only 24 of them were connected and formed a connected component.

(H) A histogram showing the robustness of the 149 GEPs in AtRootGGM, as measured by the percentage of genes from a GEP that also formed a module within AtRootGGM2. Those with more than 50% were considered robust.

The Harmony algorithm provided only cell embeddings after integration; therefore, we also used the variance stabilizing transformation (VST) algorithm to process the three scRNA-seq datasets and generate an integrated gene expression matrix (Hafemeister and Satija, 2019). The matrix, which contained 21 626 genes and 22 499 cells, was subsequently used for single-cell gene co-expression network analysis via our SingleCellGGM algorithm (Xu et al., 2023). The algorithm calculates pcors between gene pairs to evaluate their co-expression strength. In total, 158 672 co-expressed gene pairs with pcor ≥ 0.03 and co-expressed in ≥10 cells were identified (Supplemental Table 1). The selected pcor cutoff was stringent, as it landed at the far-right end of the pcor distribution curve (Supplemental Figure 1), and a permutation-based analysis estimated an associated false discovery rate (FDR) of 0.0078. The gene pairs were used to construct a gene co-expression network named AtRootGGM. By comparison, because of the reduced number of cells, separate analyses of the Jean-Baptiste, Wendrich, and Zhang datasets, required higher pcor cutoff values of 0.038, 0.041, and 0.041, respectively, to achieve an FDR level of 0.05, resulting in the identification of only 49 766, 48 333, and 54 004 co-expressed gene pairs. These results underscore the improvement in co-expression network coverage achieved by integrating single-cell datasets.

The AtRootGGM network was then clustered into 149 gene co-expression modules with ≥15 genes via the Markov cluster algorithm (Van Dongen, 2008), and we considered these modules to be GEPs because they contained co-expressed genes that might be regulated by the same transcriptional regulators (Figure 1E; Supplemental Table 2). A gene ontology (GO) enrichment analysis revealed that 125 of these 149 GEPs had enriched GO terms (Benjamini–Hochberg-adjusted P ≤ 0.05), highlighting their biological significance (Supplemental Table 3).

As a comparison, we analyzed the same integrated Arabidopsis root single-cell dataset using hdWGCNA, a method adapted from WGCNA for construction of single-cell gene co-expression networks (Langfelder and Horvath, 2008; Morabito et al., 2023). Previously, this method was used to analyze a root single-cell dataset of Medicago truncatula (Pereira et al., 2024). Using hdWGCNA, we identified 58 gene co-expression modules, 55 of which exhibited enriched GO terms (P ≤ 0.05) (Supplemental Table 4), fewer than the 125 GEPs identified by SingleCellGGM. At more stringent P value cutoffs of 1E−05, 1E−10, and 1E−50, hdWGCNA identified 31, 21, and 3 modules with enriched GO terms, respectively, whereas SingleCellGGM identified 70, 39, and 10 GEPs. Moreover, certain GO terms exhibited higher enrichment in GEPs identified by SingleCellGGM than in the modules identified by hdWGCNA. For example, for the GO term “root hair cell development” (GO:0080147), GEP M2 had the highest enrichment among all SingleCellGGM GEPs, with a P value of 3.62E−16. By contrast, the highest enriched module identified by hdWGCNA, the “red” module, had a P value of 1.09E−08. Thus, SingleCellGGM demonstrated superior capability in grouping genes with this particular GO term compared with hdWGCNA. Of all 1416 GO terms with P ≤ 1E−05 in any of the GEPs or modules (a stringent cutoff was adopted to focus on highly enriched GO terms for comparison), 820 exhibited higher enrichment levels in the SingleCellGGM GEPs, whereas 596 demonstrated higher enrichment in the hdWGCNA modules (Supplemental Table 5). These findings indicate that SingleCellGGM outperformed hdWGCNA in identifying gene modules with significantly enriched GO terms for Arabidopsis roots.

To evaluate the robustness of the identified GEPs, we also used the Shahan dataset, which contains 110 427 cells, to construct a second gene co-expression network named AtRootGGM2 via SingleCellGGM. A total of 135 415 gene pairs with pcor ≥ 0.03 and co-expressed in ≥10 cells were chosen for construction of AtRootGGM2 (Supplemental Table 6) with an FDR less than 6.35E−06. We then checked whether the genes from the GEPs of AtRootGGM still formed modules in AtRootGGM2. For example, GEP M82 contains 38 genes that function mainly in “cell fate commitment” (P = 1.30E−04) in the root epidermis (Figure 1F), including CPC, MYB23, TRY, and TTG2, which regulate root epidermis cell differentiation (Bruex et al., 2012). We extracted a subnetwork for these 38 genes from AtRootGGM2 and found that its largest connected component contained 24 genes, indicating that 63% of the genes in GEP M82 from AtRootGGM still formed a module in AtRootGGM2 (Figure 1G). Among the 149 GEPs, there were 108 in which at least 50% of their genes still formed modules in AtRootGGM2, indicating that these modules were robust and were repeatedly identified from different root datasets (Figure 1H; Supplemental Table 7).

By combining the aforementioned analyses, we identified 135 GEPs from AtRootGGM that contained enriched GO terms (P ≤ 0.05) and/or at least 50% of whose genes still formed modules in AtRootGGM2 (Supplemental Table 7). We further analyzed these GEPs by examining their cell-type- and developmental-stage-specific expression patterns and identified a series of GEPs regulating root development and metabolism. Some of these GEPs were expressed in specific cell types, providing insights into the regulatory programs of root development; some had enriched functions in metabolic pathways with cell-type-specific expression; and some showed shared expression across tissues and represented common pathways among different cell types. Below, we describe these GEPs in detail.

GEPs for QC, columella, and LRC development

The QC is the source of root stem cells and functions directly in root development. We identified a GEP, M26, involved in regulating QC and stem cell niche homeostasis (Figure 2A). M26 contains 96 genes, among which WIP4, BRAVO, and RGF1 are known markers of the QC (Lee et al., 2006; Matsuzaki et al., 2010). To obtain an overview of the gene expression pattern of M26, we mapped its expression, defined as the average expression value of all genes within the GEP in every cell, and found that it was specifically expressed in the QC (Figure 2A). We then checked whether the individual genes within M26 had similar expression patterns. We identified genes with the highest intra-GEP connections from M26 and considered them to be hub genes. Hub genes share co-expression patterns with the greatest number of other genes, and their expression can be considered the most representative within the GEP. The top 20 hub genes in M26, including WIP4, AT3G02245, and AT2G28790, also showed specific expression in the QC (Figure 2B).Figure 2 GEPs for QC, columella, and LRC development.

(A) GEP M26 for QC and stem cell niche homeostasis regulation. A subnetwork for GEP M26 from AtRootGGM is shown. Dots represent genes, and connections between genes indicate that they have co-expression patterns. Names are shown for only a portion of genes owing to space limitations. Gene names in red or boldface indicate that they have a specific GO annotation or belong to a specific gene group. The map on the top right shows the GEP M26 expression pattern, which is defined as the average expression level of all the genes within the GEP in every cell. The same applies to other plots hereafter.

(B) A dot plot showing cell-type-specific expression of the hub genes in the GEPs. The top 20 hub genes were selected from each GEP to draw the plot. The genes are ordered according to their intra-GEP connections, and the gene with the greatest number of connections is located on the left.

(C) GEP M14 for proximal (young) columella. Two genes highlighted by underlining were selected for functional investigation (see the main text).

(D) GEP M1 for the distal columella and LRC.

GO analysis revealed that M26 is enriched with 17 genes related to “meristem maintenance” (P = 1.59E−13) (Supplemental Table 3). It contains at least 20 genes with known functions in the QC and the stem cell niche (Figure 2A). Among them, WOX5 is a key regulator of QC organization (Sarkar et al., 2007), BRAVO encodes a transcriptional repressor that forms a complex with WOX5 to modulate gene expression in the QC (Betegón-Putze et al., 2021), PLT1/2/3/4 regulate stem cell proliferation and differentiation (Aida et al., 2004), WIP2/4/5 act redundantly to determine columella stem cell fate (Crawford et al., 2015), RGF1/5/8 and RAFL34 encode small peptides that regulate stem cell size and primary root growth (Matsuzaki et al., 2010; Nikonorova et al., 2021), CIK4 encodes a kinase that modulates both root and shoot apical meristem signaling (Zhu et al., 2021), and LAX2, TAA1, IAA33, PIN1, and ARF5 are auxin-related genes with functions in the QC (Sabatini et al., 1999; Vidaurre et al., 2007; Zhang et al., 2013; Brumos et al., 2018; Lv et al., 2020). This module also contains genes with known functions in the shoot apical meristem, such as OBO1 and PI, but their roles in the QC have not been characterized (Bouhidel and Irish, 1996; Cho and Zambryski, 2011). The enrichment of known QC-related genes indicates that other uncharacterized genes within M26 are ideal candidates for the discovery of novel QC regulators.

In contrast to M26, M14 is expressed in both the QC and the proximal (young) columella, and its hub genes, such as AT1G56680 and AT5G02070, also display similar expression patterns (Figure 2B and 2C). M14 contains RGF2 and RGF3, two marker genes of the proximal columella, and ADF9, a marker gene of the columella (Matsuzaki et al., 2010; Taniguchi et al., 2017). M14 also contains genes that regulate auxin homeostasis and redistribution and the gravitropism response, such as YUC3, PDK1, PIN4, PIN7, DRO1, ARL2, and SGR9, consistent with the role of the columella in related biological processes (Blilou et al., 2005; Harrison and Masson, 2008; Nakamura et al., 2011; Chen et al., 2014; Taniguchi et al., 2017; Tan et al., 2020).

GEPs for more mature root-tip cells were also identified. M1 is expressed in both the distal columella and the proximal and distal LRC (Figure 2B and 2D). It is enriched with 42 genes that regulate “root development” (P = 4.38E−08), including ARF10, BRN1, and BRN2, three transcription factors (TFs) that modulate columella and LRC development (Wang et al., 2005; Bennett et al., 2010). M1 also contains nine genes involved in “starch metabolic process” (P = 9.72E−06), such as DPE1, DPE2, and APL4, consistent with the role of the columella in statolith production and gravitropism response. Moreover, GEPs with specific or preferred expression in the proximal LRC (M65), distal LRC (M41), and both the distal LRC and the distal columella (M51) were also identified (Supplemental Figure 2). Our analysis thus revealed a series of GEPs that function in development of the QC and root-tip cells at different developmental stages.

GEPs for phloem development

The root phloem poles consist of SEs, companion cells (CCs), and phloem pole pericycle (PPP) cells, which function in metabolite and nutrient transport. We identified GEPs with specific expression in all these cell types (Figure 3A). Among them, M20 is a GEP for the early stage, M22 for the middle stage, and M77 and M76 for the late stage of SE development.Figure 3 GEPs for phloem development.

(A) Expression patterns of GEPs for sieve element (SE), companion cell (CC), and phloem pole pericycle (PPP) development. The phloem portion of the UMAP (rectangular area) is shown and enlarged. The GEPs are arranged according to the developmental stages of the cells in which they are expressed, with arrows pointing toward more mature stages. The same applies to similar figures hereafter.

(B–D) GEPs M20, M22, and M77 for SEs.

M20 contains 116 genes with specific expression in the meristem and elongation protophloem (we will refer to the meristem zone of the protophloem as the meristem protophloem, and the same applies to other zones and cell types) (Figure 3B; Supplemental Figure 3A). At least 17 of these genes are known regulators of early protophloem SE development, including NAC020, an early TF regulator of SE differentiation (Kondo et al., 2016); PEAR1, PEAR2, DOF6, and TMO6, which encode C2C2-Dof TFs that regulate phloem formation (Miyashima et al., 2019); CLE25/26/45, which form a circuit with C2C2-Dof TFs to control phloem organization (Qian et al., 2022); BRX and PAX, encoding two membrane-associated proteins that regulate SE differentiation by adjusting auxin flux (Marhava et al., 2018); and OPS, CVP2, CRN, JUL2, and SMXL3/4/5, also regulators of phloem development (Truernit et al., 2012; Rodriguez-Villalon et al., 2015; Hazak et al., 2017; Wallner et al., 2017; Cho et al., 2018). M20 contains three C2C2-GATA TF genes, GATA18/19/20, but their roles in SEs have not been characterized. M20 also includes FKD1 and its homologs FL1 and FL3. FKD1 regulates PIN1 localization and auxin canalization in leaf veins (Prabhakaran Mariyamma et al., 2018), but whether FKD1, FL1, and FL3 also regulate SE development by modulating auxin flux remains to be tested. Because of its enrichment in relevant known regulators, other uncharacterized genes within M20 could also serve as candidates for functional studies on early SE development.

Genes in M22 are mainly expressed in the elongation protophloem, with some also expressed in the meristem and maturation protophloem (Figure 3C; Supplemental Figure 3A). M22 contains SE marker genes such as SEOR1, SEOR2, ENODL9, CALS7, and SUS5 (Khan et al., 2007; Xie et al., 2011; Anstead et al., 2012; Cayla et al., 2019; Yao et al., 2020), as well as the SE developmental regulators BAM3 and HCA2 (Depuydt et al., 2013; Miyashima et al., 2019). By contrast, genes in M77 and M76 are mainly expressed in the maturation protophloem, with M76 expressed at a later stage (Figure 3D; Supplemental Figure 3A). M77 contains NAC45 and NAC86 and their target gene NEN2, which function together in enucleation of SEs (Furuta et al., 2014). M76 contains NEN4, another target gene of NAC45/86 that functions in enucleation (Supplemental Figure 3B). Previous studies have shown that NEN2 is expressed at the early stage of enucleation, whereas NEN4 is expressed at a later stage (Furuta et al., 2014), consistent with the expression timing of M77 and M76. Our analysis thus revealed GEPs that are expressed chronologically along the trajectory of phloem SE development.

We also identified a series of GEPs for CC development, with M129 and M34 expressed at the early stage and M49 and M31 expressed at a later stage (Figure 3A; Supplemental Figure 3). The genes in M34 are mainly expressed in the elongation stage, with some also expressed in the meristem stage. M34 contains three CC marker genes, NAKR1, AHA3, and SUC2 (Ivashikina et al., 2003; Tian et al., 2010). It also contains APL, which encodes a G2-like TF required for phloem development (Bonke et al., 2003). Interestingly, M34 also contains four uncharacterized G2-like TF genes, HHO4, HHO5, HHO6, and SAPL. In contrast to those of M34, genes in M31 are mainly expressed in maturation CCs. M31 also contains marker genes of CCs, such as PP2-A1, SUS1, and AAP2, which respectively encode a phloem lectin protein, a sucrose synthase, and an amino acid transporter (Dinant et al., 2003; Zhang et al., 2010; Yao et al., 2020). M31 contains other transporter genes, such as AAP4, UMAMIT26, UMAMIT27, SWEET14, and PROT1, whose functions in CC require further investigation.

Another two GEPs, M16 and M39, were identified for the early and late stages of PPP cell development. Genes in M16 are mainly expressed in the elongation PPP and thus represent an early stage of PPP development (Figure 3A; Supplemental Figure 3A). GEP M16 contains two PPP marker genes, CALS8 and AT3G27030 (Ross-Elliott et al., 2017; Otero et al., 2022), as well as four small-peptide-encoding genes, CLE4/5/6/7 (Supplemental Figure 3B), but their functions in phloem development remain uncharacterized. By contrast, GEP M39 contains genes mainly expressed in both the maturation and elongation PPP, thus representing a more mature stage of PPP development (Figure 3A; Supplemental Figure 3A). M39 is enriched in 18 genes with “transporter activity” (P = 9.04E−05), such as SWEET3/11/12, UMAMIT10/12/20/21/30/31, and ZIP5 (Supplemental Figure 3B). It also contains TF genes such as NAC056, bZIP9, and SOM. Interestingly, SWEET11 and SWEET12 are marker genes of phloem parenchyma (PP) cells in Arabidopsis leaves, and they encode sugar transporters that efflux sucrose from PP cells into the apoplast as a prerequisite for phloem loading (Chen et al., 2012; Cayla et al., 2019). In an scRNA-seq study focused on leaf phloem cells, transcripts of SWEET11/12, UMAMIT12/20/21/30, NAC056, and bZIP9 were also found to be enriched in PP cells (Kim et al., 2021). Thus, at least a subset of mature PPP cells in roots and PP cells in leaves express the same set of transporters and TF regulators, indicating that these PPP cells may also function in metabolite transport.

GEPs for development of the endodermis, cortex, and other tissues

We also identified GEPs regulating endodermis and cortex development (Figure 4A). In roots, the endodermis and cortex are derived from the same cortex/endodermis initial cells (CEIs). Among the GEPs regulating endodermis and cortex development, M55 is expressed in both the meristem endodermis and meristem cortex, which include the CEIs (Figure 4B; Supplemental Figure 4A). On the basis of its expression pattern, we hypothesized that M55 might contain regulators of early ground tissue (endodermis and cortex) development. Indeed, M55 contains SCZ, which is expressed in the CEIs, cortex, and endodermis and functions as a key regulator of ground tissue formation (Pernas et al., 2010). However, most other genes in M55 remain to be tested for their potential involvement in ground tissue development.Figure 4 GEPs for endodermis and cortex development.

(A) Expression patterns of GEPs for endodermis and cortex development.

(B–E) GEPs M55, M32, M50, and M93 expressed in the endodermis at different developmental stages.

Another five GEPs, M32, M56, M12, M50, and M93, were identified for endodermis development. M32 and M56 are mainly expressed in the elongation endodermis, with some of their genes also showing expression in the mature endodermis (Figure 4C; Supplemental Figure 4). Interestingly, M32 contains three TF regulators of endodermis development, SCR, BLJ, and MYB36, highlighting the role of this GEP in endodermis formation (Di Laurenzio et al., 1996; Liberman et al., 2015; Moreno-Risueno et al., 2015). M12, M50, and M93 are specifically expressed in the mature endodermis. M50 is enriched with seven genes involved in “suberin biosynthesis” (P = 2.38E−12), including CYP86B1, FAR4, and KCS2 (Figure 4D; Supplemental Figure 4). M50 also contains known TF regulators of suberin biosynthesis, such as MYB93, MYB39, MYB9, and MYB107 (Lashbrooke et al., 2016; Wang et al., 2020; Shukla et al., 2021). M93 mainly functions in Casparian strip formation. It is enriched with six genes with the GO term “Casparian strip” (P = 1.50E−11), such as CASP1–4 and ESB1, and also contains other genes required for Casparian strip formation, such as PRX9, PRX64, PRX72, and SGN1 (Figure 4E) (Alassimone et al., 2016; Rojas-Murcia et al., 2020). The identification of these two GEPs from endodermal cells is consistent, in that both suberin biosynthesis and Casparian strip formation are specialized pathways of the endodermis.

GEPs M98, M9, and M48 were mainly expressed in the cortex at different developmental stages, providing candidate genes for the study of cortex development and function (Figure 4A; Supplemental Figure 4). We also identified GEPs regulating the development of the procambium, xylem, atrichoblasts, and trichoblasts (Supplemental Figures 5and 6). Thus, from the AtRootGGM gene co-expression network, we were able to identify GEPs regulating the development of different root tissues at different developmental stages. These GEPs paint a comprehensive picture of root cell development and provide numerous candidate genes for further studies.

GEPs for metabolism and shared functions across cell types

We also identified GEPs that participate in metabolic pathways, a number of which exhibit cell-type-specific expression (Figure 5A and 5B; Supplemental Figure 7A). For example, M83 is enriched with 18 genes involved in “glucosinolate biosynthetic process” (P = 9.19E−31); this GEP is specifically expressed in a subset of maturation xylem pole pericycle (XPP) and PPP cells, indicating that the glucosinolate biosynthesis pathway is active in these cells. M40 is enriched with 32 genes for “lipid biosynthetic process” (P = 1.24E−07) and is specifically expressed in distal LRC cells. M70 is enriched with 16 genes for “flavonoid biosynthetic process” (P = 2.09E−24) and is mainly expressed in cortex cells. M74 is enriched with 12 genes for “tryptophan metabolic process” (P = 1.73E−20) and is expressed across multiple tissues, including QC and pericycle cells. M130 is enriched with 9 genes for “terpenoid metabolic process” (P = 1.83E−11), and it is expressed in atrichoblasts and the LRC. M24 is enriched with 19 genes for “phenylpropanoid metabolic process” (p = 9.75E−21) and is expressed in the procambium, XPP, and PPP. Thus, our analysis revealed cell-type-specific metabolic pathways in Arabidopsis roots at the transcriptional level, and further investigation is warranted to determine whether related metabolites have similar specific distributions.Figure 5 GEPs for metabolic pathways and GEPs shared across different cell types.

(A and B) Expression patterns and enriched GO terms of GEPs involved in different metabolic pathways.

(C and D) Expression patterns and enriched GO terms of GEPs shared across different cell types.

We identified GEPs with shared expression across different root cell types, which may represent common pathways active in roots (Figure 5C and 5D; Supplemental Figure 7B). These GEPs include M6, with 78 genes for “cell cycle” (P = 1.32E−66); M13, with 42 genes for “DNA replication” (P = 6.35E−55); M58, with 12 genes for “chromatin organization” (P = 4.99E−11); M11, with 53 genes for “ribonucleoprotein complex biogenesis” (P = 6.20E−49); M119, with 3 genes for “sulfate assimilation” (P = 7.89E−05); and M54, with 29 genes for “response to endoplasmic reticulum stress” (P = 9.73E−50). Thus, our algorithm identified both GEPs with specific expression in particular cell types and GEPs with shared expression across different cell types.

Validation of GEP expression patterns using reporter lines

As indicated in the aforementioned analysis, a subset of GEPs exhibited specific spatiotemporal expression patterns. To substantiate these findings, we used transgenic reporter lines. Fifteen hub genes were chosen from five GEPs, and Arabidopsis transgenic lines expressing the Histone2B-sGFP (H2B-sGFP) reporter gene were generated under the control of the corresponding gene promoters (1.5–2 kb upstream of ATG). The expression of these genes in roots was monitored using the reporter lines, and the results are summarized below.

As anticipated, three genes from the QC-specific GEP M26 demonstrated QC-specific expression patterns (Figure 6A–6C). Among these, AT3G02245 and AT2G28790 are novel, whereas SAUR51 has previously shown leaf and root primordium-specific expression (van Mourik et al., 2017). Three genes from GEP M20 exhibited early SE-specific expression patterns (Figure 6D–6F); GATA19 and GATA20 remain uncharacterized, whereas JUL2 has been recognized as a key regulator of phloem differentiation, although its expression pattern has not been reported (Cho et al., 2018). Two genes from GEP M55, AT4G16447 and AT3G04520, were expressed in cortex and endodermis initials, consistent with expression of this GEP in the meristematic cortex and endodermis. Notably, their expression extended into the early endodermis or the QC area (Figure 6G and 6H). Four genes from the CC-specific GEP M34, FTIP1, AT1G10380, AT1G02705, and AT2G02000, exhibited similar expression patterns in CCs (Figure 6I–6L). Among these genes, only FTIP1 has previously been associated with CC-specific expression (Liu et al., 2012). Finally, three uncharacterized genes from M2, a GEP specific for early trichoblast development (AT1G53680, AT1G27140, and AT4G22217), also exhibited specific expression patterns in developing root hair cells (Figure 6M–6O).Figure 6 Validation of GEP expression patterns in roots.

Promoter-driven histone2B-sGFP (H2B-sGFP, green) reporters were used to characterize the expression patterns of 15 hub genes selected from GEPs M26 (A–C), M20 (D–F), M55 (G and H), M34 (I–L), and M2 (M–O). Samples were stained with propidium iodide (red) before observation. In (I)–(L), arrows indicate cells in which the GFP signal initially appears, and stars in (M)–(O) represent trichoblasts. Scale bars, 50 μm.

In summary, the majority of genes from the same GEPs exhibited similar patterns in roots, which also aligned with those of the respective GEPs derived from single-cell datasets. These findings underscore the potential of these genes to serve as cell-type-specific markers for Arabidopsis roots, thereby establishing the GEPs as a valuable resource for identifying novel root marker genes.

GEPs enable identification of a regulator gene of root gravitropism

The GEPs revealed by our analysis can be used to identify candidate genes for root biology studies. As an illustration, we investigated the functions of two uncharacterized genes, AT5G02070 and AT5G48130, from GEP M14 (Figure 2C). Given that GEP M14 exhibits specific expression in proximal columella cells and includes known gravitropism genes, we hypothesized that these two genes might also function in the root gravitropism response. Initial screening revealed that mutation in at5g02070 did not affect the root gravitropism response (see below). Consequently, our focus turned to AT5G48130, a member of the NPH3/RPT2-like (NRL) gene family. The NRL family was named after two founding members, NONPHOTOTROPIC HYPOCOTYL3 (NPH3) and ROOT PHOTOTROPISM2 (RPT2), that function in hypocotyl and root phototropism (Motchoulski and Liscum, 1999; Sakai et al., 2000). Five other NRL genes (NRL6, NRL7, NRL20, NRL21, and NRL30) have been implicated in the regulation of auxin movements crucial for organogenesis and root gravitropism (Furutani et al., 2007; Li et al., 2011; Glanc et al., 2021). Triple, quadruple, or quintuple mutants of these genes exhibit various degrees of gravitropic defects (Li et al., 2011). Previously named NRL27, AT5G48130 remains functionally uncharacterized (Christie et al., 2018), and we will refer to this gene as NRL27 hereafter.

We first examined whether the expression of NRL27 matched that of GEP M14. We generated Arabidopsis transgenic lines expressing the β-D-glucuronidase (GUS) gene driven by the NRL27 promoter. GUS staining indicated that NRL27 was specifically expressed in both primary and lateral root tips (Figure 7A and 7B). Another transgenic line expressing the NRL27–GFP fusion protein driven by the NRL27 promoter showed that NRL27 was specifically expressed in proximal columella cells (Figure 7C). Thus, the expression pattern of NRL27 is consistent with that of GEP M14.Figure 7 NRL27 contributes to root gravitropism.

(A and B) Expression patterns of NRL27. Eight-day-old Arabidopsis NRL27::GUS transgenic plants were stained for GUS expression. Scale bars, 1 mm in (A) and 50 μm in (B).

(C) NRL27 is localized exclusively in columella cells. Arabidopsis NRL27::NRL27-GFP transgenic plants were observed under a confocal microscope. Scale bars, 50 μm.

(D) Schematic representation of the root gravitropism response assay. Arabidopsis seedlings were grown vertically for 5 days in the direction of G1 gravity and for another 24 h in the direction of G2 gravity, and the root angle α (green arrow) was measured as an indicator of the gravitropism response.

(E and F) Gravitropism response of WT, nrl27 mutant (nrl27-1, nrl27-2, and nrl27-3), NRL27-OE, and complementation (#6 and #8) lines (E) and their root angle measurements (F). The experiment and measurements were performed as described in (D). Scale bars in (E), 1 cm. Bars represent means ± SD in (F); n = 50–91; ns, not significantly different; ∗∗∗∗P < 0.001; Student’s t-test.

(G and H) Confocal images (G) and measurements (H) of green fluorescence in DR5rev::GFP/WT and DR5rev::GFP/nrl27-1 lines, as an indicator of local auxin level. Scale bars in (G), 20 μm. Bars represent means ± SD in (H); n = 37–39; ∗∗P < 0.05; Student’s t-test.

(I and J) Gravitropism response of WT, nrl27, and NRL27-OE lines growing on 1/2 MS agar plates supplemented with 10 nM IAA (I) and their root angle measurements (J). Scale bars in (I), 1 cm. Bars represent means ± SD in (J); n = 31–53; ns, not significantly different; Student’s t-test.

To study the physiological functions of NRL27, we obtained a transfer DNA (T-DNA) mutant line (CS872834) harboring a T-DNA insertion in the third exon of NRL27, which we named nrl27-1 (Supplemental Figure 8A). Another two independent mutant lines, nrl27-2 and nrl27-3, were generated using CRISPR technology. These two lines had a 1-bp insertion and a 1-bp deletion in the first exon of NRL27 near the single-guide RNA site, respectively, both of which resulted in early termination of NRL27 protein biosynthesis (Supplemental Figure 8B and 8C). We grew the mutant and wild-type (WT) plants on vertical 1/2 MS agar plates for 5 days, changed the direction of gravity force by rotating the plates 90°, and observed the root gravitropism response after 24 h (Figure 7D and 7E). In this assay, roots with normal gravitropism form a nearly 90° angle, whereas those without gravitropism have an angle close to 180°. Compared with the WT, all three mutant lines displayed significant root gravitropism defects: the average angle of the WT roots was 94.9°, whereas those of nrl27-1, nrl27-2, and nrl27-3 were 124.4°, 136.3°, and 146.0°, respectively, all significantly larger than that of the WT (Figure 7F). By contrast, a similar analysis of a T-DNA insertion mutant of at5g02070 (SALK_043598C), also from GEP M14, revealed that at5g02070 mutation did not affect root gravitropism (Supplemental Figure 9).

To confirm that the gravitropism defects in nrl27 were caused by the NRL27 mutation, we generated complementation lines expressing the NRL27 coding sequence driven by the NRL27 promoter in the nrl27-1 background. Two complementation lines (#6 and #8) rescued the gravitropic phenotype, and their root angles were similar to those of WT plants after the change in gravitropic force direction (Figure 7E and 7F). We also generated an NRL27 overexpression line, NRL27-OE, in which NRL27 was overexpressed 65 fold (Supplemental Figure 10A). Interestingly, NRL27-OE also displayed gravitropism defects (Figure 7E and 7F). These results indicated that the expression of NRL27 must be maintained at a proper level for the root gravitropism response, and both the deletion and the overexpression of this gene led to gravitropism defects. We examined the subcellular localization of the NRL27 protein and found that it mainly localized to the plasma membrane (Supplemental Figure 10B), similar to other NRL proteins (Motchoulski and Liscum, 1999; Sakai et al., 2000; Furutani et al., 2011).

We then investigated whether the gravitropism phenotypes of the nrl27 mutants were related to auxin. By crossing nrl27-1 with the auxin reporter line DR5rev::GFP (Friml et al., 2003), we found that the auxin level of columella cells was significantly lower in nrl27-1 than in WT plants (Figure 7G and 7H). Interestingly, growth of the mutants and the overexpression line in 1/2 MS agar plates supplemented with 10 nM indole-3-acetic acid (IAA) rescued the gravitropism defects to the WT level (Figure 7I and 7J). These results indicated that NRL27 regulates root gravitropism response via an auxin-related pathway. Furthermore, the nrl27 mutants displayed longer primary roots than WT plants when grown on 1/2 MS agar plates, suggesting a role for NRL27 in root length regulation as well (Supplemental Figure 10C and 10D). This case study demonstrates that our single-cell gene co-expression network can be used to pinpoint candidate genes for investigations of root biology.

Discussion

Multiple scRNA-seq studies have identified the major cell types in Arabidopsis roots and characterized their developmental trajectories, but the GEPs underlying the related developmental processes remain to be fully characterized. We integrated three Arabidopsis root scRNA-seq datasets, performed a gene co-expression network analysis using the SingleCellGGM algorithm, and identified 149 GEPs. Notably, although single-cell co-expression network analysis can be performed on individual datasets, relying solely on individual datasets may limit the total number of cells and result in identification of fewer co-expressed gene pairs. By examining their spatiotemporal expression patterns, we identified a series of GEPs with specific expression in different root cells along their developmental trajectories. For example, GEPs M20, M22, M77, and M76 are expressed in SEs at different developmental stages; M3, M25, M59, M62, and M29 in protoxylem; and M2, M28, and M17 in trichoblasts (Figure 3, Supplemental Figures 5 and 6). Moreover, five GEPs—M32, M56, M12, M50, and M93—are implicated in endodermis development, revealing the heterogeneity of this intricate tissue (Figure 4). We next used reporter lines to validate the cell-type-specific expression patterns of hub genes from the identified GEPs, demonstrating their reliability as a source for novel marker gene discovery (Figure 6). In addition, these GEPs are enriched in relevant developmental regulators and provide numerous candidates for future studies. As an example, M20, which is expressed in the meristem and elongation SEs, is enriched with at least 17 regulators of SE development (Figure 3B), covering most regulators known to date (Hardtke, 2023). The uncharacterized genes in M20 thus represent promising candidates for identification of novel SE developmental regulators.

Our analysis also revealed GEPs involved in metabolic pathways, some of which exhibit cell-type-specific expression in roots. Among them, M83 for glucosinolate biosynthesis is specifically expressed in a subset of maturation XPP and PPP cells (Figure 5A). In Arabidopsis leaves, a group of glucosinolate-rich cells (S cells) is found in floral stems and in the PP of the phloem, and they participate in defense against herbivores and pathogens (Koroleva et al., 2010). However, glucosinolate is not produced in S cells because biosynthetic enzymes cannot be detected in isolated S-cell extracts. The cells that produce glucosinolate in leaves remain to be identified. It will be interesting to investigate whether glucosinolate is produced in the root XPP and PPP cells that express M83, whether the resulting glucosinolate participates in the defense response, and whether cells similar to S cells are present in Arabidopsis roots. The relevance of other metabolism-related GEPs, such as M130 for terpenoid biosynthesis in atrichoblasts and LRCs and M70 for flavonoid biosynthesis in the cortex, also warrants further investigation.

Our analysis also revealed GEPs with shared expression across different root cell types, associated with crucial pathways for root functions such as DNA replication, cell-cycle regulation, and sulfate assimilation (Figure 5). These shared GEPs may have been overlooked in previous single-cell analyses owing to the absence of an efficient algorithm for identifying them. Our analysis thus facilitates the study of gene programs shared across various cell types in roots.

Finally, we illustrated the practical utility of the GEPs for root biology studies with a specific example. We pinpointed two candidate genes from the columella-specific GEP M14 and verified that one of the two candidates, NRL27, functions in root gravitropism. NRL27 exhibits specific expression in proximal columella cells (Figure 7A–7C), consistent with the expression pattern of M14. Given that GEP M14 includes several known root gravitropism genes, we hypothesized that NRL27 might play a similar role. Our results confirmed that NRL27 indeed regulates the auxin-related root gravitropism response (Figure 7D–7J). Notably, a single mutation of NRL27 resulted in significant root gravitropic defects, a phenotype observed only in triple or higher-order mutants of the five other NRL genes (NRL6, NRL7, NRL20, NRL21, and NRL30) previously implicated in gravitropism (Li et al., 2011). This uniqueness suggests that NRL27 has a distinct function that cannot be replaced by the other NRL genes, possibly owing to its unique expression pattern (Supplemental Figure 11) or its distinct phylogenetic position (NRL27 and the other five NRLs belong to different clades in a phylogenetic tree of the NRL protein family in Arabidopsis) (Christie et al., 2018).

In conclusion, our analysis systematically identified the GEPs that function in Arabidopsis roots. As more scRNA-seq datasets become available, similar analyses can be extended to cover other plant tissues and species, which will provide a systematic overview of the GEPs that regulate plant development and metabolism.

Methods

scRNA-seq dataset processing and integration

The raw sequence data of three Arabidopsis root scRNA-seq datasets from Jean-Baptiste et al., Wendrich et al., and Zhang et al. (Jean-Baptiste et al., 2019; Zhang et al., 2019; Wendrich et al., 2020) were downloaded from the NCBI GEO and SRA databases (GEO: GSE121619 and GSE141730; SRA: PRJNA517021). Six scRNA-seq samples, all generated by the 10× Chromium platform, were used for the analysis, including one (SAMN10818242) from Zhang et al., two (GSM4466787 and GSM4466788) from Jean-Baptiste et al., and three (GSM4212550, GSM4212551, and GSM4212552) from Wendrich et al. Sequencing reads were mapped to the Arabidopsis genome (araport11) to obtain gene expression count matrices using CellRanger (v.5.0.1). Cells with <1000 or ≥5000 expressed genes or with ≥5% of reads mapped to mitochondrial genes were filtered out. The gene expression counts of the 22 499 remaining cells were used to create a merged single-cell dataset using Seurat (v.4.2.0) in R (Stuart et al., 2019).

The merged dataset was normalized and scaled while regressing out the effects of the cell cycle on gene expression using a list of cell-cycle genes obtained from Zhang et al. (2021). Three thousand highly variable genes (HVGs) were identified from the dataset. To mitigate the effects of protoplasting during scRNA-seq sample preparation, genes affected by protoplasting (|Log2FC| ≥ 2 and FDR < 0.05 according to Denyer et al., 2019) were removed from the HVG list if present. The remaining 1975 HVGs were used as variable features to perform a principal-component analysis of the dataset. The dataset was then integrated via Harmony (v.0.1.1) (Korsunsky et al., 2019), using the original sample IDs of the cells as batch variables. A UMAP plot of the cells after integration was generated on the basis of the Harmony reduction using the RunUMAP function in Seurat with the parameters “dims=1:50, min.dist=0.6”.

Cell-type and developmental-stage annotations

We used the cell label transfer method provided by Seurat to annotate the cells in the merged dataset (Stuart et al., 2019). A reference dataset (GSE152766_Root_Atlas_seu4.rds.gz) produced by Shahan et al. was obtained from GEO (GEO: GSE152766) (Shahan et al., 2022). This dataset contains expert curated annotations with both developmental stage and cell type information for 110 427 Arabidopsis root cells, for example, “Meristem_Cortex” and “Elongation_Trichoblast.” We mapped the merged dataset to the Shahan et al. dataset and transferred the annotations from the Shahan et al. dataset to the merged dataset using the TransferData function in Seurat, with the parameters “dims=1:50, k.weight=10”. After obtaining developmental stage–cell type annotations for the cells, we separated the annotations into cell type annotations and developmental stage annotations for downstream analysis.

Single-cell gene co-expression network analysis

We processed the three scRNA-seq datasets using the VST algorithm provided by the SCTransform R package (Hafemeister and Satija, 2019) to obtain an integrated gene expression matrix. The merged gene count matrix from the three scRNA-seq datasets was processed using the vst function in SCTransform with the parameters “batch_var = ‘orig.ident’, min_cells = 20”, where ‘orgi.ident’ represents the original sample IDs of the cells. An integrated gene expression count matrix corrected for sequencing depth and batch effects was subsequently obtained. The integrated count matrix was log-transformed via log1p and used for gene co-expression network analysis via the SingleCellGGM algorithm (https://github.com/MaShisongLab/SingleCellGGM) (Xu et al., 2023). The algorithm calculates pcors between all gene pairs. For each gene pair, it also checks the number of cells in which both genes are expressed. It also provides an “fdr_control” function to estimate the FDR at different pcor cutoff levels. We chose 158 672 gene pairs with pcor ≥ 0.03 that were co-expressed in ≥10 cells to construct the AtRootGGM gene co-expression network, with an FDR rate of 0.0078.

Identification and annotation of GEPs

The AtRootGGM network was clustered via the Markov cluster algorithm with the parameters “-I 1.7 -scheme 7 -te 20” to obtain gene co-expression modules (Van Dongen, 2008), which we considered to be GEPs. We extracted subnetworks for the GEPs from AtRootGGM and visualized AtRootGGM and the subnetworks in Cytoscape (v.3.8.3) (Shannon et al., 2003). Hub genes were identified from the GEPs according to their number of intra-GEP connections, and the top 20 hub genes in each GEP were used for dot plot visualization. The average expression values of all genes within every GEP in each cell were calculated and considered to be the GEP expression value. The expression patterns of the GEPs were visualized to identify GEPs with cell-type- and developmental-stage-specific expression. GO annotations of Arabidopsis genes were obtained from TAIR as of January 1, 2023 (Berardini et al., 2015) and used to perform GO enrichment analysis for GEPs based on the hypergeometric distribution.

hdWGCNA gene co-expression network analysis

The hdWGCNA package was employed to analyze the same integrated dataset used to construct the AtRootGGM network (Morabito et al., 2023). The analysis began with the integrated and batch-corrected gene expression count matrix, as described in the preceding steps. Default parameters of hdWGCNA were applied, with the following adjustments: (1) all genes and all cells in the matrix were used for network construction; (2) metacells were constructed using the aforementioned Harmony cell dimension reduction and cell type annotations, with the parameters “min_cells = 50, max_shared = 10”; and (3) the minimal module size was set to 15, consistent with the size limit applied for the GEPs identified by SingleCellGGM. The modules identified by hdWGCNA were then analyzed for GO enrichment and compared with the GEPs.

Plant materials and CRISPR mutant lines

Arabidopsis thaliana lines used in this study were derived from the Col-0 background. To confirm the expression patterns of hub genes from selected GEPs, Arabidopsis H2B-sGFP reporter lines driven by their respective promoters (1.5–2 kb) were generated for the following genes using the TQ662 vector (Zhang et al., 2021): AT3G02245, AT2G28790, SAUR51 (AT1G75580), JUL2 (AT5G25490), GATA19 (AT4G36620), GATA20 (AT2G18380), AT4G16447, AT3G04520, FTIP1 (AT5G06850), AT1G10380, AT1G02705, AT2G02000, AT1G53680, AT1G27140, and AT4G22217. The primers used to construct the transformation vectors are listed in Supplemental Table 8. The DR5rev::GFP marker line was described previously (Friml et al., 2003). The nrl27-1 (CS872834) and at5g02070 (SALK_043598C) mutant lines were obtained from TAIR. The NRL27::GUS transgenic line was generated using the pCAMBIA1391Z vector containing a native NRL27 gene promoter (1380 bp). The NRL27::NRL27-GFP transgenic line and the NRL27::NRL27/nrl27-1 complementation lines were generated using the pCAMBIA1300-nos vector with the same NRL27 promoter. The NRL27-OE (35S::NRL27) line was generated using the pCAMBIA1300 vector.

We used a pHY01 vector developed in our laboratory to generate CRISPR-based NRL27 mutants. The pHY01 vector is similar to another vector, pHY07, which was described in Geng et al. (2023). These two vectors were designed to construct CRISPR gene knockouts in Arabidopsis, and both vectors use a promoter from the NUC-1 gene to drive expression of the zCas9 gene (Wang et al., 2015) to improve gene knockout efficiency. pHY01 is an earlier version of pHY07, and it differs from pHY07 only in having lower efficiency of mCherry marker gene expression for visual identification of transgenic seeds. A single-guide RNA (CAAGGATTGTTCGTCCG) was designed via CRISPR-P 2.0 (Liu et al., 2017) and cloned into the pHY01 vector to target NRL27. The vector was transformed into Agrobacterium tumefaciens GV3101, which was subsequently used to transform WT Arabidopsis via the floral dip method (Clough and Bent, 1998). The T1 transgenic plants were sequenced to identify lines harboring mutations in NRL27, and these plants were self-crossed to obtain homozygous mutants.

Gravitropism response assay

Seeds were surface sterilized and stratified at 4°C for 2 days and subsequently grown on vertical 1/2 MS medium plates supplemented with 1% (w/v) sucrose (pH 5.7) in a growth chamber at 21°C under long-day conditions (16-h light/8-h dark). After 5 days, the plates were rotated 90°, and the root angles were measured after another 24 h. For the IAA treatment, the plants were grown on plates supplemented with 10 nM IAA sodium salt (Sigma, I5148) for 5 days and then subjected to gravitropism response assays.

Confocal microscopy

Fluorescence imaging was performed using the Zeiss LSM 710 and LSM 880 confocal microscopes. The manufacturer’s default settings were used to image proteins tagged with GFP (excitation, 488 nm; emission, 495–545 nm). For the H2B-sGFP reporter lines, roots from 5-day-old seedlings were mounted in 10 μg/mL propidium iodide (Sigma) and imaged for both GFP and propidium iodide (excitation, 543 nm; emission, 589–781 nm). Images were analyzed using ImageJ software.

Data and code availability

All single-cell datasets of Arabidopsis roots used in this study are publicly available and were downloaded from the NCBI GEO and SRA databases as cited above.

Funding

This work was supported by grants from the Strategic Priority Research Program of the Chinese Academy of Science (XDA24010303 ), the 10.13039/501100001809 National Natural Science Foundation of China (31770268 ), the Fundamental Research Funds for the Central Universities (WK2070000091 ), and the 10.13039/501100009076 University of Science and Technology of China (start-up fund to S.M.).

Author contributions

S.M. designed and supervised the project. E.H., Z.G., Y.Q., Y.W., and S.M. performed the experiments and analysis. E.H. and S.M. wrote the manuscript. All authors reviewed and approved the final manuscript.

Supplemental information

Document S1. Supplemental Figures 1–11

Supplemental Tables 1–8

Document S2. Article plus supplemental information

Acknowledgments

We thank the USTC Supercomputing Center and USTC School of Life Sciences Bioinformatics Center for providing the computing resources. No conflict of interest is declared.

Published by the Plant Communications Shanghai Editorial Office in association with Cell Press, an imprint of Elsevier Inc., on behalf of CSPB and CEMPS, CAS.

Supplemental information is available at Plant Communications Online.
==== Refs
References

Aida M. Beis D. Heidstra R. Willemsen V. Blilou I. Galinha C. Nussaume L. Noh Y.S. Amasino R. Scheres B. The PLETHORA genes mediate patterning of the Arabidopsis root stem cell niche Cell 119 2004 109 120 10.1016/j.cell.2004.09.018 15454085
Alassimone J. Fujita S. Doblas V.G. van Dop M. Barberon M. Kalmbach L. Vermeer J.E.M. Rojas-Murcia N. Santuari L. Hardtke C.S. Geldner N. Polarly localized kinase SGN1 is required for Casparian strip integrity and positioning Nat. Plants 2 2016 16113 10.1038/nplants.2016.113 27455051
Anstead J.A. Froelich D.R. Knoblauch M. Thompson G.A. Arabidopsis P-protein filament formation requires both AtSEOR1 and AtSEOR2 Plant Cell Physiol. 53 2012 1033 1042 10.1093/pcp/pcs046 22470058
Baran Y. Bercovich A. Sebe-Pedros A. Lubling Y. Giladi A. Chomsky E. Meir Z. Hoichman M. Lifshitz A. Tanay A. MetaCell: analysis of single-cell RNA-seq data using K-nn graph partitions Genome Biol. 20 2019 206 10.1186/s13059-019-1812-2 31604482
Bennett T. van den Toorn A. Sanchez-Perez G.F. Campilho A. Willemsen V. Snel B. Scheres B. SOMBRERO, BEARSKIN1, and BEARSKIN2 regulate root cap maturation in Arabidopsis Plant Cell 22 2010 640 654 10.1105/tpc.109.072272 20197506
Berardini T.Z. Reiser L. Li D. Mezheritsky Y. Muller R. Strait E. Huala E. The Arabidopsis information resource: Making and mining the "gold standard" annotated reference plant genome Genesis 53 2015 474 485 10.1002/dvg.22877 26201819
Betegón-Putze I. Mercadal J. Bosch N. Planas-Riverola A. Marquès-Bueno M. Vilarrasa-Blasi J. Frigola D. Burkart R.C. Martínez C. Conesa A. Precise transcriptional control of cellular quiescence by BRAVO/WOX5 complex in Arabidopsis roots Mol. Syst. Biol. 17 2021 e9864 10.15252/msb.20209864 34132490
Blilou I. Xu J. Wildwater M. Willemsen V. Paponov I. Friml J. Heidstra R. Aida M. Palme K. Scheres B. The PIN auxin efflux facilitator network controls growth and patterning in Arabidopsis roots Nature 433 2005 39 44 10.1038/nature03184 15635403
Bonke M. Thitamadee S. Mähönen A.P. Hauser M.T. Helariutta Y. APL regulates vascular tissue identity in Arabidopsis Nature 426 2003 181 186 10.1038/nature02100 14614507
Bouhidel K. Irish V.F. Cellular interactions mediated by the homeotic PISTILLATA gene determine cell fate in the Arabidopsis flower Dev. Biol. 174 1996 22 31 10.1006/dbio.1996.0048 8626018
Brady 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 318 2007 801 806 10.1126/science.1146265 17975066
Bruex A. Kainkaryam R.M. Wieckowski Y. Kang Y.H. Bernhardt C. Xia Y. Zheng X. Wang J.Y. Lee M.M. Benfey P. A Gene Regulatory Network for Root Epidermis Cell Differentiation in Arabidopsis PLoS Genet. 8 2012 e1002446 10.1371/journal.pgen.1002446
Brumos J. Robles L.M. Yun J. Vu T.C. Jackson S. Alonso J.M. Stepanova A.N. Local Auxin Biosynthesis Is a Key Regulator of Plant Development Dev. Cell 47 2018 306 318.e5 10.1016/j.devcel.2018.09.022 30415657
Cayla T. Le Hir R. Dinant S. Live-Cell Imaging of Fluorescently Tagged Phloem Proteins with Confocal Microscopy Methods Mol. Biol. 2014 2019 95 108 10.1007/978-1-4939-9562-2_8 31197789
Chen L.Q. Qu X.Q. Hou B.H. Sosso D. Osorio S. Fernie A.R. Frommer W.B. Sucrose efflux mediated by SWEET proteins as a key step for phloem transport Science 335 2012 207 211 10.1126/science.1213351 22157085
Chen Q. Dai X. De-Paoli H. Cheng Y. Takebayashi Y. Kasahara H. Kamiya Y. Zhao Y. Auxin Overproduction in Shoots Cannot Rescue Auxin Deficiencies in Arabidopsis Roots Plant Cell Physiol. 55 2014 1072 1079 10.1093/pcp/pcu039 24562917
Cho E. Zambryski P.C. ORGAN BOUNDARY1 defines a gene expressed at the junction between the shoot apical meristem and lateral organs Proc. Natl. Acad. Sci. USA 108 2011 2154 2159 10.1073/pnas.1018542108 21245300
Cho H. Cho H.S. Nam H. Jo H. Yoon J. Park C. Dang T.V.T. Kim E. Jeong J. Park S. Translational control of phloem development by RNA G-quadruplex–JULGI determines plant sink strength Nat. Plants 4 2018 376 390 10.1038/s41477-018-0157-2 29808026
Christie J.M. Suetsugu N. Sullivan S. Wada M. Shining Light on the Function of NPH3/RPT2-Like Proteins in Phototropin Signaling Plant Physiol. 176 2018 1015 1024 10.1104/pp.17.00835 28720608
Clough S.J. Bent A.F. Floral dip: a simplified method for Agrobacterium-mediated transformation of Arabidopsis thaliana Plant J. 16 1998 735 743 10.1046/j.1365-313x.1998.00343.x 10069079
Crawford B.C.W. Sewell J. Golembeski G. Roshan C. Long J.A. Yanofsky M.F. Plant development. Genetic control of distal stem cell fate within root and embryonic meristems Science 347 2015 655 659 10.1126/science.aaa0196 25612610
Crow M. Gillis J. Co-expression in Single-Cell Analysis: Saving Grace or Original Sin? Trends Genet. 34 2018 823 831 10.1016/j.tig.2018.07.007 30146183
Denyer T. Ma X. Klesen S. Scacchi E. Nieselt K. Timmermans M.C.P. Spatiotemporal Developmental Trajectories in the Arabidopsis Root Revealed Using High-Throughput Single-Cell RNA Sequencing Dev. Cell 48 2019 840 852.e5 10.1016/j.devcel.2019.02.022 30913408
Depuydt S. Rodriguez-Villalon A. Santuari L. Wyser-Rmili C. Ragni L. Hardtke C.S. Suppression of Arabidopsis protophloem differentiation and root meristem growth by CLE45 requires the receptor-like kinase BAM3 Proc. Natl. Acad. Sci. USA 110 2013 7074 7079 10.1073/pnas.1222314110 23569225
Di Laurenzio L. Wysocka-Diller J. Malamy J.E. Pysh L. Helariutta Y. Freshour G. Hahn M.G. Feldmann K.A. Benfey P.N. The SCARECROW gene regulates an asymmetric cell division that is essential for generating the radial organization of the Arabidopsis root Cell 86 1996 423 433 10.1016/s0092-8674(00)80115-4 8756724
Dinant S. Clark A.M. Zhu Y. Vilaine F. Palauqui J.C. Kusiak C. Thompson G.A. Diversity of the superfamily of phloem lectins (phloem protein 2) in angiosperms Plant Physiol. 131 2003 114 128 10.1104/pp.013086 12529520
Dolan L. Janmaat K. Willemsen V. Linstead P. Poethig S. Roberts K. Scheres B. Cellular organisation of the Arabidopsis thaliana root Development 119 1993 71 84 8275865
Friml J. Vieten A. Sauer M. Weijers D. Schwarz H. Hamann T. Offringa R. Jürgens G. Efflux-dependent auxin gradients establish the apical–basal axis of Arabidopsis Nature 426 2003 147 153 10.1038/nature02085 14614497
Furuta K.M. Yadav S.R. Lehesranta S. Belevich I. Miyashima S. Heo J.O. Vatén A. Lindgren O. De Rybel B. Van Isterdael G. Plant development. Arabidopsis NAC45/86 direct sieve element morphogenesis culminating in enucleation Science 345 2014 933 937 10.1126/science.1253736 25081480
Furutani M. Kajiwara T. Kato T. Treml B.S. Stockum C. Torres-Ruiz R.A. Tasaka M. The gene MACCHI-BOU 4/ENHANCER OF PINOID encodes a NPH3-like protein and reveals similarities between organogenesis and phototropism at the molecular level Development 134 2007 3849 3859 10.1242/dev.009654 17913786
Furutani M. Sakamoto N. Yoshida S. Kajiwara T. Robert H.S. Friml J. Tasaka M. Polar-localized NPH3-like proteins regulate polarity and endocytosis of PIN-FORMED auxin efflux carriers Development 138 2011 2069 2078 10.1242/dev.057745 21490067
Gala H.P. Lanctot A. Jean-Baptiste K. Guiziou S. Chu J.C. Zemke J.E. George W. Queitsch C. Cuperus J.T. Nemhauser J.L. A single-cell view of the transcriptome during lateral root initiation in Arabidopsis thaliana Plant Cell 33 2021 2197 2220 10.1093/plcell/koab101 33822225
Geng H. Wang Y. Xu Y. Zhang Y. Han E. Peng Y. Geng Z. Liu Y. Qin Y. Ma S. Data-driven optimization yielded a highly-efficient CRISPR/Cas9 system for gene editing in Arabidopsis Preprint at bioRxiv 2023 10.1101/2023.10.09.561629
Glanc M. Van Gelderen K. Hoermayer L. Tan S. Naramoto S. Zhang X. Domjan D. Včelařová L. Hauschild R. Johnson A. AGC kinases and MAB4/MEL proteins maintain PIN polarity by limiting lateral diffusion in plant cells Curr. Biol. 31 2021 1918 1930.e5 10.1016/j.cub.2021.02.028 33705718
Hafemeister C. Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression Genome Biol. 20 2019 296 10.1186/s13059-019-1874-1 31870423
Han X. Wang R. Zhou Y. Fei L. Sun H. Lai S. Saadatpour A. Zhou Z. Chen H. Ye F. Mapping the Mouse Cell Atlas by Microwell-Seq Cell 172 2018 1091 1107.e17 10.1016/j.cell.2018.02.001 29474909
Hardtke C.S. Phloem development New Phytol. 239 2023 852 867 10.1111/nph.19003 37243530
Harrison B.R. Masson P.H. ARL2, ARG1 and PIN3 define a gravity signal transduction pathway in root statocytes Plant J. 53 2008 380 392 10.1111/j.1365-313X.2007.03351.x 18047472
Hazak O. Brandt B. Cattaneo P. Santiago J. Rodriguez-Villalon A. Hothorn M. Hardtke C.S. Perception of root-active CLE peptides requires CORYNE function in the phloem vasculature EMBO Rep. 18 2017 1367 1381 10.15252/embr.201643535 28607033
Ivashikina N. Deeken R. Ache P. Kranz E. Pommerrenig B. Sauer N. Hedrich R. Isolation of AtSUC2 promoter-GFP-marked companion cells for patch-clamp studies and expression profiling Plant J. 36 2003 931 945 10.1046/j.1365-313x.2003.01931.x 14675456
Jean-Baptiste K. McFaline-Figueroa J.L. Alexandre C.M. Dorrity M.W. Saunders L. Bubb K.L. Trapnell C. Fields S. Queitsch C. Cuperus J.T. Dynamics of Gene Expression in Single Root Cells of Arabidopsis thaliana Plant Cell 31 2019 993 1011 10.1105/tpc.18.00785 30923229
Khan J.A. Wang Q. Sjölund R.D. Schulz A. Thompson G.A. An early nodulin-like protein accumulates in the sieve element plasma membrane of Arabidopsis Plant Physiol. 143 2007 1576 1589 10.1104/pp.106.092296 17293437
Kim J.-Y. Symeonidi E. Pang T.Y. Denyer T. Weidauer D. Bezrutczyk M. Miras M. Zöllner N. Hartwig T. Wudick M.M. Distinct identities of leaf phloem cells revealed by single cell transcriptomics Plant Cell 33 2021 511 530 10.1093/plcell/koaa060 33955487
Kondo Y. Nurani A.M. Saito C. Ichihashi Y. Saito M. Yamazaki K. Mitsuda N. Ohme-Takagi M. Fukuda H. Vascular Cell Induction Culture System Using Arabidopsis Leaves (VISUAL) Reveals the Sequential Differentiation of Sieve Element-Like Cells Plant Cell 28 2016 1250 1262 10.1105/tpc.16.00027 27194709
Koroleva O.A. Gibson T.M. Cramer R. Stain C. Glucosinolate-accumulating S-cells in Arabidopsis leaves and flower stalks undergo programmed cell death at early stages of differentiation Plant J. 64 2010 456 469 10.1111/j.1365-313X.2010.04339.x 20815819
Korsunsky I. Millard N. Fan J. Slowikowski K. Zhang F. Wei K. Baglaenko Y. Brenner M. Loh P.-r. Raychaudhuri S. Fast, sensitive and accurate integration of single-cell data with Harmony Nat. Methods 16 2019 1289 1296 10.1038/s41592-019-0619-0 31740819
La Manno G. Siletti K. Furlan A. Gyllborg D. Vinsland E. Mossi Albiach A. Mattsson Langseth C. Khven I. Lederer A.R. Dratva L.M. Molecular architecture of the developing mouse brain Nature 596 2021 92 96 10.1038/s41586-021-03775-x 34321664
Lahnemann D. Koster J. Szczurek E. McCarthy D.J. Hicks S.C. Robinson M.D. Vallejos C.A. Campbell K.R. Beerenwinkel N. Mahfouz A. Eleven grand challenges in single-cell data science Genome Biol. 2110 2020 1186/s13059-020-1926-6
Langfelder P. Horvath S. WGCNA: an R package for weighted correlation network analysis BMC Bioinf. 9 2008 559 10.1186/1471-2105-9-559
Lashbrooke J. Cohen H. Levy-Samocha D. Tzfadia O. Panizel I. Zeisler V. Massalha H. Stern A. Trainotti L. Schreiber L. MYB107 and MYB9 Homologs Regulate Suberin Deposition in Angiosperms Plant Cell 28 2016 2097 2116 10.1105/tpc.16.00490 27604696
Lee J.Y. Colinas J. Wang J.Y. Mace D. Ohler U. Benfey P.N. Transcriptional and posttranscriptional regulation of transcription factor expression in Arabidopsis roots Proc. Natl. Acad. Sci. USA 103 2006 6055 6060 10.1073/pnas.0510607103 16581911
Li 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 39 2016 508 522 10.1016/j.devcel.2016.10.012 27840108
Li W.V. Li J.J. An accurate and robust imputation method scImpute for single-cell RNA-seq data Nat. Commun. 9 2018 997 10.1038/s41467-018-03405-7 29520097
Li Y. Dai X. Cheng Y. Zhao Y. NPY genes play an essential role in root gravitropic responses in Arabidopsis Mol. Plant 4 2011 171 179 10.1093/mp/ssq052 20833732
Liberman L.M. Sparks E.E. Moreno-Risueno M.A. Petricka J.J. Benfey P.N. MYB36 regulates the transition from proliferation to differentiation in the Arabidopsis root Proc. Natl. Acad. Sci. USA 112 2015 12099 12104 10.1073/pnas.1515576112 26371322
Liu H. Ding Y. Zhou Y. Jin W. Xie K. Chen L.L. CRISPR-P 2.0: An Improved CRISPR-Cas9 Tool for Genome Editing in Plants Mol. Plant 10 2017 530 532 10.1016/j.molp.2017.01.003 28089950
Liu L. Liu C. Hou X. Xi W. Shen L. Tao Z. Wang Y. Yu H. FTIP1 is an essential regulator required for florigen transport PLoS Biol. 10 2012 e1001313 10.1371/journal.pbio.1001313 22529749
Lv B. Yu Q. Liu J. Wen X. Yan Z. Hu K. Li H. Kong X. Li C. Tian H. Non-canonical AUX/IAA protein IAA33 competes with canonical AUX/IAA repressor IAA5 to negatively regulate auxin signaling EMBO J. 39 2020 e101515 10.15252/embj.2019101515 31617603
Ma S. Gong Q. Bohnert H.J. An Arabidopsis gene network based on the graphical Gaussian model Genome Res. 17 2007 1614 1625 10.1101/gr.6911207 17921353
Marhava P. Bassukas A.E.L. Zourelidou M. Kolb M. Moret B. Fastner A. Schulze W.X. Cattaneo P. Hammes U.Z. Schwechheimer C. Hardtke C.S. A molecular rheostat adjusts auxin flux to promote root protophloem differentiation Nature 558 2018 297 300 10.1038/s41586-018-0186-z 29875411
Matsuzaki Y. Ogawa-Ohnishi M. Mori A. Matsubayashi Y. Secreted peptide signals required for maintenance of root stem cell niche in Arabidopsis Science 329 2010 1065 1067 10.1126/science.1191132 20798316
Mentzen W.I. Wurtele E.S. Regulon organization of Arabidopsis BMC Plant Biol. 8 2008 99 10.1186/1471-2229-8-99 18826618
Miyashima S. Roszak P. Sevilem I. Toyokura K. Blob B. Heo J.O. Mellor N. Help-Rinta-Rahko H. Otero S. Smet W. Mobile PEAR transcription factors integrate positional cues to prime cambial growth Nature 565 2019 490 494 10.1038/s41586-018-0839-y 30626969
Morabito S. Reese F. Rahimzadeh N. Miyoshi E. Swarup V. hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data Cell Rep. Methods 3 2023 100498 10.1016/j.crmeth.2023.100498 37426759
Moreno-Risueno M.A. Sozzani R. Yardımcı G.G. Petricka J.J. Vernoux T. Blilou I. Alonso J. Winter C.M. Ohler U. Scheres B. Benfey P.N. Transcriptional control of tissue formation throughout root development Science 350 2015 426 430 10.1126/science.aad1171 26494755
Motchoulski A. Liscum E. Arabidopsis NPH3: A NPH1 photoreceptor-interacting protein essential for phototropism Science 286 1999 961 964 10.1126/science.286.5441.961 10542152
Nakamura M. Toyota M. Tasaka M. Morita M.T. An Arabidopsis E3 ligase, SHOOT GRAVITROPISM9, modulates the interaction between statoliths and F-actin in gravity sensing Plant Cell 23 2011 1830 1848 10.1105/tpc.110.079442 21602290
Nikonorova N. Murphy E. Fonseca de Lima C.F. Zhu S. van de Cotte B. Vu L.D. Balcerowicz D. Li L. Kong X. De Rop G. The Arabidopsis Root Tip (Phospho)Proteomes at Growth-Promoting versus Growth-Repressing Conditions Reveal Novel Root Growth Regulators Cells 2021
Nolan T.M. Vukašinović N. Hsu C.W. Zhang J. Vanhoutte I. Shahan R. Taylor I.W. Greenstreet L. Heitz M. Afanassiev A. Brassinosteroid gene regulatory networks at cellular resolution in the Arabidopsis root Science 379 2023 eadf4721 10.1126/science.adf4721
Otero S. Gildea I. Roszak P. Lu Y. Di Vittori V. Bourdon M. Kalmbach L. Blob B. Heo J.-o. Peruzzo F. A root phloem pole cell atlas reveals common transcriptional states in protophloem-adjacent cells Nat. Plants 8 2022 954 970 10.1038/s41477-022-01178-y 35927456
Pereira W.J. Boyd J. Conde D. Triozzi P.M. Balmant K.M. Dervinis C. Schmidt H.W. Boaventura-Novaes C. Chakraborty S. Knaack S.A. The single-cell transcriptome program of nodule development cellular lineages in Medicago truncatula Cell Rep. 43 2024 113747 10.1016/j.celrep.2024.113747 38329875
Pernas M. Ryan E. Dolan L. SCHIZORIZA controls tissue system complexity in plants Curr. Biol. 20 2010 818 823 10.1016/j.cub.2010.02.062 20417101
Prabhakaran Mariyamma N. Clarke K.J. Yu H. Wilton E.E. Van Dyk J. Hou H. Schultz E.A. Members of the Arabidopsis FORKED1-LIKE gene family act to localize PIN1 in developing veins J. Exp. Bot. 69 2018 4773 4790 10.1093/jxb/ery248 29982821
Qian P. Song W. Zaizen-Iida M. Kume S. Wang G. Zhang Y. Kinoshita-Tsujimura K. Chai J. Kakimoto T. A Dof-CLE circuit controls phloem organization Nat. Plants 8 2022 817 827 10.1038/s41477-022-01176-0 35817820
Rodriguez-Villalon A. Gujas B. van Wijk R. Munnik T. Hardtke C.S. Primary root protophloem differentiation requires balanced phosphatidylinositol-4,5-biphosphate levels and systemically affects root branching Development 142 2015 1437 1446 10.1242/dev.118364 25813544
Rojas-Murcia N. Hématy K. Lee Y. Emonet A. Ursache R. Fujita S. De Bellis D. Geldner N. High-order mutants reveal an essential requirement for peroxidases but not laccases in Casparian strip lignification Proc. Natl. Acad. Sci. USA 117 2020 29166 29177 10.1073/pnas.2012728117 33139576
Ross-Elliott T.J. Jensen K.H. Haaning K.S. Wager B.M. Knoblauch J. Howell A.H. Mullendore D.L. Monteith A.G. Paultre D. Yan D. Phloem unloading in Arabidopsis roots is convective and regulated by the phloem-pole pericycle Elife 6 2017 10.7554/eLife.24125
Ryu K.H. Huang L. Kang H.M. Schiefelbein J. Single-Cell RNA Sequencing Resolves Molecular Relationships Among Individual Plant Cells Plant Physiol. 179 2019 1444 1456 10.1104/pp.18.01482 30718350
Sabatini S. Beis D. Wolkenfelt H. Murfett J. Guilfoyle T. Malamy J. Benfey P. Leyser O. Bechtold N. Weisbeek P. Scheres B. An Auxin-Dependent Distal Organizer of Pattern and Polarity in the Arabidopsis Root Cell 99 1999 463 472 10.1016/S0092-8674(00)81535-4 10589675
Sakai T. Wada T. Ishiguro S. Okada K. RPT2. A signal transducer of the phototropic response in Arabidopsis Plant Cell 12 2000 225 236 10.1105/tpc.12.2.225 10662859
Sarkar A.K. Luijten M. Miyashima S. Lenhard M. Hashimoto T. Nakajima K. Scheres B. Heidstra R. Laux T. Conserved factors regulate signalling in Arabidopsis thaliana shoot and root stem cell organizers Nature 446 2007 811 814 10.1038/nature05703 17429400
Schafer J. Strimmer K. A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics Stat. Appl. Genet. Mol. Biol. 4 2005 Article32 10.2202/1544-6115.1175
Serrano-Ron L. Perez-Garcia P. Sanchez-Corrionero A. Gude I. Cabrera J. Ip P.L. Birnbaum K.D. Moreno-Risueno M.A. Reconstruction of lateral root formation through single-cell RNA sequencing reveals order of tissue initiation Mol. Plant 14 2021 1362 1378 10.1016/j.molp.2021.05.028 34062316
Shahan R. Hsu C.W. Nolan T.M. Cole B.J. Taylor I.W. Greenstreet L. Zhang S. Afanassiev A. Vlot A.H.C. Schiebinger G. A single-cell Arabidopsis root atlas reveals developmental trajectories in wild-type and cell identity mutants Dev. Cell 57 2022 543 560.e9 10.1016/j.devcel.2022.01.008 35134336
Shannon P. Markiel A. Ozier O. Baliga N.S. Wang J.T. Ramage D. Amin N. Schwikowski B. Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks Genome Res. 13 2003 2498 2504 10.1101/gr.1239303 14597658
Shukla V. Han J.P. Cléard F. Lefebvre-Legendre L. Gully K. Flis P. Berhin A. Andersen T.G. Salt D.E. Nawrath C. Barberon M. Suberin plasticity to developmental and exogenous cues is regulated by a set of MYB transcription factors Proc Natl Acad Sci USA 118 2021 10.1073/pnas.2101730118
Shulse C.N. Cole B.J. Ciobanu D. Lin J. Yoshinaga Y. Gouran M. Turco G.M. Zhu Y. O'Malley R.C. Brady S.M. Dickel D.E. High-Throughput Single-Cell Transcriptome Profiling of Plant Cell Types Cell Rep. 27 2019 2241 2247.e4 10.1016/j.celrep.2019.04.054 31091459
Stuart T. Butler A. Hoffman P. Hafemeister C. Papalexi E. Mauck W.M. Hao Y. Stoeckius M. Smibert P. Satija R. Comprehensive Integration of Single-Cell Data Cell 177 2019 1888 1902.e21 10.1016/j.cell.2019.05.031 31178118
Tan S. Zhang X. Kong W. Yang X.L. Molnár G. Vondráková Z. Filepová R. Petrášek J. Friml J. Xue H.W. The lipid code-dependent phosphoswitch PDK1-D6PK activates PIN-mediated auxin efflux in Arabidopsis Nat. Plants 6 2020 556 569 10.1038/s41477-020-0648-9 32393881
Taniguchi M. Furutani M. Nishimura T. Nakamura M. Fushita T. Iijima K. Baba K. Tanaka H. Toyota M. Tasaka M. Morita M.T. The Arabidopsis LAZY1 Family Plays a Key Role in Gravity Signaling within Statocytes and in Branch Angle Control of Roots and Shoots Plant Cell 29 2017 1984 1999 10.1105/tpc.16.00575 28765510
Tian H. Baxter I.R. Lahner B. Reinders A. Salt D.E. Ward J.M. Arabidopsis NPCC6/NaKR1 is a phloem mobile metal binding protein necessary for phloem function and root meristem maintenance Plant Cell 22 2010 3963 3979 10.1105/tpc.110.080010 21193571
Truernit E. Bauby H. Belcram K. Barthélémy J. Palauqui J.C. OCTOPUS, a polarly localised membrane-associated protein, regulates phloem differentiation entry in Arabidopsis thaliana Development 139 2012 1306 1315 10.1242/dev.072629 22395740
Van Dongen S. Graph Clustering Via a Discrete Uncoupling Process SIAM J. Matrix Anal. Appl. 30 2008 121 141 10.1137/040608635
van Mourik H. van Dijk A.D.J. Stortenbeker N. Angenent G.C. Bemer M. Divergent regulation of Arabidopsis SAUR genes: a focus on the SAUR10-clade BMC Plant Biol. 17 2017 245 10.1186/s12870-017-1210-4 29258424
Vidaurre D.P. Ploense S. Krogan N.T. Berleth T. AMP1 and MP antagonistically regulate embryo and meristem development in Arabidopsis Development 134 2007 2561 2567 10.1242/dev.006759 17553903
Wallner E.S. López-Salmerón V. Belevich I. Poschet G. Jung I. Grünwald K. Sevilem I. Jokitalo E. Hell R. Helariutta Y. Strigolactone- and Karrikin-Independent SMXL Proteins Are Central Regulators of Phloem Formation Curr. Biol. 27 2017 1241 1247 10.1016/j.cub.2017.03.014 28392107
Wang C. Wang H. Li P. Li H. Xu C. Cohen H. Aharoni A. Wu S. Developmental programs interact with abscisic acid to coordinate root suberization in Arabidopsis Plant J. 104 2020 241 251 10.1111/tpj.14920 32645747
Wang J.W. Wang L.J. Mao Y.B. Cai W.J. Xue H.W. Chen X.Y. Control of root cap formation by MicroRNA-targeted auxin response factors in Arabidopsis Plant Cell 17 2005 2204 2216 10.1105/tpc.105.033076 16006581
Wang Z.P. Xing H.L. Dong L. Zhang H.Y. Han C.Y. Wang X.C. Chen Q.J. Egg cell-specific promoter-controlled CRISPR/Cas9 efficiently generates homozygous mutants for multiple target genes in Arabidopsis in a single generation Genome Biol. 16 2015 144 10.1186/s13059-015-0715-0 26193878
Wendrich J.R. Yang B. Vandamme N. Verstaen K. Smet W. Van de Velde C. Minne M. Wybouw B. Mor E. Arents H.E. Vascular transcription factors guide plant epidermal responses to limiting phosphate conditions Science 370 2020 10.1126/science.aay4970
Wille A. Zimmermann P. Vranova E. Furholz A. Laule O. Bleuler S. Hennig L. Prelic A. von Rohr P. Thiele L. Sparse graphical Gaussian modeling of the isoprenoid gene network in Arabidopsis thaliana Genome Biol. 5 2004 10.1186/gb-2004-5-11-r92 [pii]
Xie B. Wang X. Zhu M. Zhang Z. Hong Z. CalS7 encodes a callose synthase responsible for callose deposition in the phloem Plant J. 65 2011 1 14 10.1111/j.1365-313X.2010.04399.x 21175885
Xu Y. Wang Y. Ma S. SingleCellGGM enables gene expression program identification from single-cell transcriptomes and facilitates universal cell label transfer Preprint at bioRxiv 2023 10.1101/2023.02.05.526424
Yao D. Gonzales-Vigil E. Mansfield S.D. Arabidopsis sucrose synthase localization indicates a primary role in sucrose translocation in phloem J. Exp. Bot. 71 2020 1858 1869 10.1093/jxb/erz539 31805187
Zhang L. Tan Q. Lee R. Trethewy A. Lee Y.-H. Tegeder M. Altered Xylem-Phloem Transfer of Amino Acids Affects Metabolism and Leads to Increased Seed Yield and Oil Content in Arabidopsis Plant Cell 22 2010 3603 3620 10.1105/tpc.110.073833 21075769
Zhang T.Q. Chen Y. Wang J.W. A single-cell analysis of the Arabidopsis vegetative shoot apex Dev. Cell 56 2021 1056 1074.e8 10.1016/j.devcel.2021.02.021 33725481
Zhang T.Q. Xu Z.G. Shang G.D. Wang J.W. A Single-Cell RNA Sequencing Profiles the Developmental Landscape of Arabidopsis Root Mol. Plant 12 2019 648 660 10.1016/j.molp.2019.04.004 31004836
Zhang W. Swarup R. Bennett M. Schaller G.E. Kieber J.J. Cytokinin Induces Cell Division in the Quiescent Center of the Arabidopsis Root Apical Meristem Curr. Biol. 23 2013 1979 1989 10.1016/j.cub.2013.08.008 24120642
Zhang Y. Han E. Peng Y. Wang Y. Wang Y. Geng Z. Xu Y. Geng H. Qian Y. Ma S. Rice co-expression network analysis identifies gene modules associated with agronomic traits Plant Physiol. 190 2022 1526 1542 10.1093/plphys/kiac339 35866684
Zhu Y. Hu C. Cui Y. Zeng L. Li S. Zhu M. Meng F. Huang S. Long L. Yi J. Conserved and differentiated functions of CIK receptor kinases in modulating stem cell signaling in Arabidopsis Mol. Plant 14 2021 1119 1134 10.1016/j.molp.2021.04.001 33823234
