An environmentally regulated dichotomy in cell lineages of the root

The root of Arabidopsis thaliana is a valuable model for studying cell differentiation in plants. Cell types are arranged in concentric layers around a central vascular bundle and arise from distinct stem cell populations maintained near the root tip in the meristem¹. As stem cells divide, their descendants are displaced shootward and progressively differentiate, which creates developmental trajectories in which the spatial arrangement of cells recapitulates the temporal dynamics of differentiation (Fig. 1a).

Single-cell RNA sequencing (scRNA-seq) has refined the transcriptional signatures of Arabidopsis root cell types and enabled high-resolution reconstruction of these developmental trajectories^{2,3,4,5,6,7,8,9}. Although these studies have identified transient cell states in individual lineages, they have largely analysed each lineage independently, which may overlook transcriptional states shared across multiple cell types irrespective of developmental identity. Whether such shared transcriptional states exist and how they contribute to root cell identity remain largely unexplored.

To characterize cell states during root development, here we generated a high-resolution scRNA-seq (10x Genomics) atlas of the A. thaliana Col-0 root tip (Fig. 1b). After filtering low-quality cells, doublets and genes significantly affected by protoplast isolation, batch integration generated a reference atlas comprising 87,056 cells from 9 replicates, with a total of 25,784 detected genes (Fig. 1b, Supplementary Tables 1 and 2 and Supplementary Fig. 1).

We identified 27 transcriptionally distinct clusters and assigned their identities by combining differential expression with optimized whole-mount hybridization chain reaction (HCR)¹⁰ and tissue clearing¹¹ to spatially map the expression of the 15 most specific marker genes (Fig. 1b–d). This process resolved all known cell types, including distinct root cap subpopulations (early, tip and lateral) (Fig. 1b–d). One cluster located between the root cap and epidermis could not be assigned directly by marker expression. However, on the basis of its position and developmental origin of these tissues, it was annotated as ‘epidermis–meristem’ (Fig. 1b).

Cells were grouped according to their developmental lineages (epidermis, root cap, cortex, endodermis, xylem pole pericycle (XPP), phloem pole pericycle (PPP), procambium, phloem and xylem), which confirmed that ontogenic programs remain the dominant determinant of root cell identity (Fig. 1b,e). Cells from each lineage were analysed independently to reconstruct developmental trajectories. A meristematic score, derived from genes enriched in the meristematic epidermis cluster, was calculated for each lineage, and the cell with the highest score was defined as the trajectory starting point (t = 0). Trajectory topology was inferred using partition-based graph abstraction (PAGA), and developmental progression was ordered by diffusion pseudotime from t = 0 (Fig. 1e–g).

The root cap lineage exhibited the most complex topology and recapitulated known developmental events. That is, the split of lateral and tip populations and their convergence with columella cells to form the sloughing border-like cells¹² (Fig. 1g). These trajectories were supported by stage-specific marker genes^13,14,15 (Extended Data Fig. 1a), which further validated the inferred developmental trajectories.

Although some lineages are anticipated to be composed of distinct cell subtypes, such as phloem and xylem, which are made of metaphloem and protophloem and metaxylem and protoxylem, respectively, and the epidermis, which subdivides into root hairs and non-root hairs, each is expected to follow a single developmental trajectory from undifferentiated to differentiated states. Notably, PAGA revealed that all cell types and subcell types instead exhibit two branches (branch 1 and branch 2) that diverge at an early pseudotime point (Fig. 1g). Only the epidermis contained cells preceding this bifurcation and a clearly identifiable meristematic population. In most other lineages, branching occurred directly from t = 0, which probably reflects the difficulty in capturing the rarer inner-layer meristematic cells during protoplasting.

To independently validate the inferred topology, whole-mount HCR was used to detect early and late markers in epidermis, endodermis and XPP (Extended Data Fig. 1b–e). In all three lineages, genes predicted to be expressed near t = 0 were detected closer to the meristem than predicted late markers. Consistent with this finding, lineage-specific meristem, elongation and maturation gene signatures obtained from a previous study⁸ progressed as expected along pseudotime. This result provides support for the robustness of our t = 0 assignment and inferred developmental ordering (Extended Data Fig. 2a).

Independent trajectory inference using Monocle also recovered the bifid topology (Extended Data Fig. 2b–d). Cell ordering in both branches was highly correlated with diffusion pseudotime and supported by RNA velocity analyses (Extended Data Fig. 2e–g).

Together, these data show that in addition to developmental progression, cells in most cell types exhibit two alternative transcriptional states. This result therefore reveals the existence of an additional, previously unrecognized, level of complexity in cell identity.

To identify transcriptional programs associated with cell type, developmental progression and branch identity, we performed high-dimensional weighted gene co-expression network analysis (hdWGCNA). Of the 6,649 genes assigned to 25 identified co-expression modules (Fig. 2a,b and Supplementary Table 3), most were enriched in specific cell types, groups of cell types or developmental stages. By contrast, modules 20 and 21 were specifically enriched in branch 1 across all lineages. This result demonstrates that there is a transcriptional program shared across multiple cell types independently of developmental identity.

Across all cell types, branch 1 cells consistently contained more detected genes and higher unique molecular identifier (UMI) counts than branch 2 cells (Fig. 2c), which indicated greater overall transcriptional activity. Modules 20 and 21 were also among the largest co-expression modules identified (320 and 615 genes, respectively; Fig. 2d), which showed that this transcriptional asymmetry involves a broad coordinated gene expression program. Although the PAGA analysis did not resolve two branches in non-hair cells, enrichment of the branch-specific modules and increased transcript abundance suggest that both populations are also present in this lineage.

Given the unexpected transcriptional differences between the two branches of the same cell type, we next sought to exclude possible experimental or analytical artefacts that might have created spurious trajectory bifurcations. Branch separation was not explained by cell cycle phase, experimental batches or genes affected by protoplast isolation (Extended Data Fig. 3a–e). Because roots grown on the surface of agar experience distinct agar-facing and air-facing microenvironments, we also profiled root tips from plants grown fully embedded in agar. The bifid trajectories remained clearly detectable across all cell types, which therefore ruled out these growth microenvironments as a contributing factor (Extended Data Fig. 3f–h).

Another possibility was that the lower UMI branch (branch 2) reflected damaged cells from which fewer transcripts were captured, despite these cells not exhibiting elevated mitochondrial or chloroplastic transcript content (Extended Data Fig. 3i–j). If so, lowly expressed genes would be preferentially detected in the higher UMI branch 1, whereas highly expressed genes would be shared between both branches. Using the XPP, we found that depletion of branch-1-specific genes in branch 2 substantially exceeded that expected from expression level alone (Extended Data Fig. 4a–c). Likewise, branch separation persisted after randomly downsampling branch 1 cells to match the UMI distribution of their branch 2 counterparts and reprocessing the dataset in three independent replicates (Extended Data Fig. 4d–g). Together, these analyses indicate that the lack of activation of the branch 1 transcriptional program in branch 2 cells cannot be explained solely by their lower transcript abundance. Finally, both branches and the branch-1-specific modules were also detected in a previously published Arabidopsis root tip single-nucleus RNA sequencing (snRNA-seq) dataset¹⁶ (Extended Data Fig. 5). This result rules out technical differences between cells and nuclei, protoplast isolation or cell damage as explanations for the observed bifurcation.

To further demonstrate the robustness of the observed branches, we re-analysed two independent previously published Arabidopsis root tip scRNA-seq datasets^7,8 (Extended Data Figs. 6 and 7). In all three cases (one snRNA-seq dataset and two scRNA-seq datasets), we recovered both branches in all cell types, together with the higher UMI and gene content of branch 1 cells and enrichment of modules 20 and 21. Indeed, branch assignment remained robust across a range of label-transfer parameters (Extended Data Figs. 5–7). The first scRNA-seq study⁷ did not initially identify the branches because many branch 2 cells were excluded owing to their lower UMI content. Conversely, the second study⁸, which integrates two additional published root scRNA-seq datasets^3,4, clearly showed two branches in their uniform manifold approximation and projection (UMAP) plots and noted their presence in the epidermis, cortex and endodermis. However, these findings were not explored further.

Some of our libraries were generated from an equimolar mix of cells from Col-0 and a second ecotype, which enabled experimental determination of cell multiplet rates while simultaneously generating root tip atlases for four additional Arabidopsis ecotypes (C24, Cvi-0, Ler and Ws-2; Methods, Extended Data Fig. 8a and Supplementary Table 2). Both branches were detected across all cell lineages in all ecotypes, with the same co-expression modules enriched in branch 1 (Extended Data Fig. 8b,c). Moreover, fewer than 4% of branch-1-specific module genes were differentially expressed in any ecotype, which demonstrated that the transcriptional program underlying the branching is highly conserved (Extended Data Fig. 8e,f). Among the ecotypes analysed, only Cvi-0 showed a significantly lower proportion of branch 1 cells (relative to branch 2) compared with Col-0 and was the ecotype with the greatest overall transcriptional divergence (Extended Data Fig. 8d–f)

A previous snRNA-seq study¹⁶ also generated a snATAC–seq atlas of the root tip. Using this chromatin accessibility data alone, we identified cells corresponding to each branch, with branch 1 cells displaying elevated fractions of reads in peaks and increased transcription start site enrichment, consistent with greater transcriptional activity (Extended Data Fig. 9a–c). Notably, this increase was specific to genes belonging to the branch-1-specific modules. Through this orthogonal method, this finding confirms the transcriptional asymmetry underlying the dichotomy identified by the scRNA-seq trajectories (Extended Data Fig. 9d,e).

To independently validate this transcriptional state in vivo, we measured the transcripts of five branch-1-specific genes by whole-mount HCR in intact root tips: AT2G41660 (also known as MIZ1), AT3G04010, AT4G26080 (also known as ABI1), AT5G56550 (also known as OXS3) and AT4G24060 (also known as VDOF1) (Fig. 3a and Extended Data Fig. 10a,b). All five genes showed heterogeneous expression in the cell types in which they are expressed, with some cells exhibiting strong transcript enrichment (positive) and others showing depletion (negative). Volumetric fluorescence quantification of segmented epidermal, cortex and endodermal cells confirmed up to 20-fold differences in fluorescence intensity between cells in the same tissue layer (Fig. 3a). Canonical cell-type markers remained detectable in negative cells, which indicated that reduced branch 1 marker detection was unlikely to reflect a general failure of HCR detection (Extended Data Fig. 10c). Although branch 2 lacks a distinct transcriptional program, analyses of epidermal gene expression identified AT3G01720 as largely complementary to the branch 1 marker AT5G56550, a pattern confirmed by HCR in vivo (Fig. 3b–h). Volumetric fluorescence quantification of segmented cells across three independent roots, together with Gaussian mixture modelling of both scRNA-seq scores and HCR fluorescence ratios, identified two discrete epidermal populations corresponding to branch 1 and branch 2 cells. This result provides direct in vivo evidence for the two transcriptionally distinct cell populations in the same tissue layer (Fig. 3e–h).

Together, these independent lines of evidence demonstrate that the bifid topology reflects a reproducible biological state characterized by the activation of a coordinated transcriptional program in a subset of cells across multiple root cell types.

To investigate the biological characteristics and potential regulators of this branch-specific transcriptional program, we analysed gene ontology (GO) term enrichment in the two branch-1-specific modules 20 and 21 (Fig. 4a and Supplementary Table 4). These modules were enriched for responses to environmental stimuli, including nutrients, water and oxygen. Terms linked to positive regulation of RNA synthesis were also enriched, consistent with the observed increase in transcriptional activity of branch 1 cells. Among the signalling pathways analysed, only ABA, a central regulator of environmental responses and the circadian clock (an integrator of multiple environmental signals), were enriched (Fig. 4a and Supplementary Table 4). Consequently, this cellular state will henceforth be referred to as the environment-regulated state (ERS), with branch 1 and branch 2 corresponding to ERS-positive (ERS⁺) and ERS-negative (ERS^–) cells, respectively. These terms reflect the presence or absence of this transcriptional state at the single-cell level rather than any intrinsic difference in the capacity of individual cells to respond to environmental signals.

To complement this enrichment analysis, we investigated transcription factors (TFs) that could potentially regulate the ERS. Using an Arabidopsis DAP-seq dataset¹⁷ and predicted binding similarities from Cis-BP2 (ref. ¹⁸), we predicted putative targets for 800 TFs, filtered using hdWGCNA with the XGBoost algorithm. This process produced putative regulons for 186 TFs expressed in our atlas (average expression scores per cell identity are shown in Fig. 4b). Notably, 84 of these TFs (>45%) showed regulon enrichment in ERS⁺ cells, a result consistent with the elevated activity of this state (Fig. 4b). These TFs included all the profiled ones linked to ABA signalling (ABF1, ABF2, ATHB5, ATHB6, ATHB7 and HAT22)^{19,20,21,22,23,24,25} or known to be circadian-regulated, such as members of the REVEILLE family (RVE1, RVE2, RVE6, EPR1 and RVE8)^26,27,28 as well as cycling DOF factors (CDF2–CDF5). Consistently, two ERS markers validated by HCR (Fig. 3a), AT4G26080 and AT2G41660, are known direct ABA targets^29,30,31. This result provides further support that ABA signalling has a role in ERS regulation.

Because ABA was the only signalling pathway identified by GO enrichment, we functionally tested its role in ERS formation by profiling an ABA biosynthesis mutant (aba2-3)³² and a quintuple ABA receptor mutant (pyl11458)³³ (Fig. 4c–e). Both mutants showed a reduced proportion of ERS⁺ cells across the whole root (Fig. 4d) and in almost all cell types independently, albeit not always significantly (Extended Data Fig. 11a). Consistently, the ERS-specific modules 20 and 21 were the most strongly downregulated co-expression modules, which confirmed that ABA signalling is involved in the establishment of the ERS (Fig. 4e).

The ERS is present in only a subset of cells regardless of lineage, which suggests that it is superimposed on, and independent from, developmental programs. To test this idea, we performed scRNA-seq to profile the root tips of the scr-3 mutant, which exhibits major defects in root cell identity through loss of SCARECROW (SCR) function³⁴ (Extended Data Fig. 11b–d). Using our Col-0 reference atlas for cell-type assignment, we observed the expected depletion of endodermis and protoxylem, together with additional alterations affecting the columella and pericycle. This result is consistent with shifts towards root cap and procambium identities, respectively, and uncovers broader cell identity changes than previously described^8,35,36,37 (Extended Data Fig. 11d). Despite these major developmental perturbations, ERS⁺ cells remained detectable across all cell populations, including those with altered identities (Extended Data Fig. 11b,c). These observations support the conclusion that the ERS represents a transcriptional state superimposed on multiple cell lineage identities regardless of cell type or developmental stage.

Our atlas incorporated a set of experiments designed to test the effect of the composition of the growth medium on cell identities. In these experiments, seedlings from three independent replicates were grown in parallel on either full-strength (1×) or half-strength (0.5×) Murashige and Skoog (MS) medium, which resulted in six out of the nine samples included in the atlas. Although this range had little effect on root growth, reducing the medium concentration induced physiological responses, including increased root hair density (Extended Data Fig. 12a). Compared with full MS, roots grown on half MS contained fewer ERS⁺ cells across most cell types, reaching significance in hair cells, cortex and protophloem (Fig. 5a,b and Extended Data Fig. 12b). Consistently, ERS-specific gene modules were significantly downregulated (Fig. 5c).

Given the well-established role of ABA in water-related stresses and the enrichment of GO terms linked to responses to water and nutrients, we speculated that larger reductions in water potential would further increase the proportion of ERS⁺ cells. To that end, we used a hard-agar model of drought to explore this idea. In this model, reduced water content induces a drought response that more closely resembles soil drying than conventional osmotic treatments³⁸.

Seedlings were grown on half MS containing 100%, 50% or 20% water (two replicates each) and profiled by scRNA-seq (Fig. 5d). Cells were assigned an identity on the basis of our reference atlas. Fewer cells were recovered from the 20% water condition (Supplementary Table 2), together with an enrichment of undifferentiated cells and an underrepresentation of larger differentiated cells (Extended Data Fig. 12c). This result is consistent with their increased fragility under severe water deficit.

We then statistically assessed changes in the proportion of ERS⁺ cells between the conditions, taking into account cells that can either be in ERS⁺ or ERS^– state, using water concentration as a continuous covariate (Fig. 5e). The proportion of ERS⁺ cells increased substantially when the water potential was reduced, from 0.38 in 100% water to 0.66 in 20% water. This result demonstrates that water deficit increases the prevalence of ERS⁺ cells (Fig. 5e). This change was significant both across the whole-root level and in individual cell types (Extended Data Fig. 12d), and was unrelated to differences in cell filtering (here and across the study) (Extended Data Fig. 12h). We observed that cells grown under water-deficit conditions contained more detected transcripts than those grown under control conditions (Extended Data Fig. 12e), which suggests that the increased prevalence of ERS⁺ cells might reflect higher transcriptional output. However, ERS-specific module scores remained elevated in water-deficit cells with comparable transcript abundance, and this difference remained significant after adjusting for log-transformed UMI counts (Extended Data Fig. 12f). Likewise, logistic regression showed that the increased odds of a cell being classified as ERS⁺ remained significant after accounting for transcript abundance (Extended Data Fig. 12g), which indicated that the increased prevalence of ERS⁺ cells cannot be fully explained by higher transcriptional output alone. Whether this result reflects changes in the probability of cells transitioning between ERS⁺ and ERS⁻ states, or more persistent occupancy of the ERS⁺ state, cannot be determined from the current data.

To further validate this finding, we examined the expression of two ERS-specific markers that we had already profiled (Fig. 3a)—AT3G04010 and AT4G26080 (which are ABA-responsive)^29,30—in vivo by whole-mount HCR in intact root tips under both the control (100% water) and water-deficit (20% water) conditions (Fig. 5f). We observed a marked increase in the number of cells expressing both these genes under the water-deficit condition (Fig. 5f). Quantifications revealed an approximate 1.6-fold increase in overall signal intensity and a 2.2–2.5-fold increase in areas showing a positive signal, independent of differences in root length (Fig. 5g and Extended Data Fig. 12i). Consistent with these measurements, empirical cumulative distribution functions showed a global shift of per-cell fluorescence towards higher values under water-deficit conditions for both gene markers (Fig. 5h). Together, these HCR analyses provide independent in planta evidence that more cells express ERS markers under water-deficit conditions, consistent with the increased prevalence of ERS⁺ cells detected by scRNA-seq.

Overall, these findings demonstrate that environmental conditions modulate the prevalence of the ERS, thereby establishing this transcriptional state as a dynamic layer of root cell identity that may contribute to adaptation to water-deficit conditions.

Through deep characterization of root cell transcriptional states coupled with HCR spatial mapping, we identified an alternative trajectory for all cell types above the root meristematic region. This trajectory, which we term the ERS, corresponds to a strong and specific transcriptional signature superimposed on a subset of cells across multiple lineages simultaneously. Cell identity in the root tip can therefore be defined by three factors: (1) the lineage that defines the cell type, which remains stable during development; (2) the developmental stage; and (3) the presence or absence of this additional transcriptional state, the prevalence of which is modulated by environmental conditions (Supplementary Fig. 2). Although detected in parallel with developmental trajectories, the ERS seems to be independent of developmental origin, as it was detected across all cell types throughout developmental progression. Instead, it represents an additional layer of identity strongly associated with ABA signalling. Our analyses also suggest that the circadian clock may influence ERS formation. This aspect will be interesting to explore further, particularly given its known feedback with ABA signalling and stress responses³⁹. We note that although an ERS can be detected across all lineages, it is not entirely equivalent in all cell types. Some ERS-specific genes were enriched in all ERS⁺ cells, whereas others were restricted to specific lineages, which suggests that there is a degree of cell-type specificity in the occurrence of this transcriptional state, which may be shaped by the underlying regulatory networks of each cell type. Differences in ABA sensitivity or signalling capacity among cells within and between cell types could equally contribute to this variability and represent an important avenue for future investigation.

A notable feature of this discovery is the scale of transcriptional heterogeneity observed in cell types, which was largely undetectable using bulk methods before single-cell technologies. Recent imaging approaches have revealed transcriptional heterogeneity in plant tissues^{40,41,42,43,44}. However, to our knowledge, none has reported the clear and consistent dichotomy observed here. The robust presence of two populations across independent datasets and detection methods, together with consistent detection of the same ERS-specific co-expressed gene modules, strongly supports the interpretation that the ERS reflects a regulated transcriptional program. Nonetheless, important questions remain regarding the mechanisms that establish and maintain this state. Exploring the temporal dynamics of the ERS will be essential, including how rapidly it is induced and how stable it remains under constant or changing environments. Moreover, it will be important to determine whether the proportions of ERS⁺ and ERS^– cells detected by scRNA-seq reflect stochastic transitions between states, the probabilities of which are modulated by environmental conditions. Live imaging of ERS marker genes will be an important avenue to address these questions. Equally important will be determining the signals that govern entry into and exit from the ERS, including whether ABA directly influences the probability that cells adopt this state or whether it permissively sets the proportion of cells competent to do so.

Another open question is the spatial arrangement and coordination of ERS⁺ cells. We showed that the ERS can occur in all root layers, independent of cell type, but the signal determining which cells enter this state is unclear. The increase in ERS prevalence under water-deficit conditions suggests that this state is part of a physiological response to the environment. If this response is linked to water and nutrient availability, it would make sense that it spans from the epidermis, which senses and absorbs these resources, to the vasculature, which distributes them throughout the plant. As such, ERS formation may contribute to, or result from, a coordinated response across the tissue. The persistence of the ERS across environmental and genetic perturbations suggests that it represents an important physiological state. Correlating ERS occurrence with known pathways or functional markers will be crucial to clarify its role.

The marked increase in ERS prevalence under water-deficit conditions demonstrates that this transcriptional state is dynamically regulated by environmental conditions. Determining which environmental cues regulate ERS prevalence, how these signals are integrated and how the ERS influences cellular physiology will be important next steps. More broadly, our findings establish the ERS as a previously unrecognized layer of root cell identity. We anticipate that the identification of the ERS will facilitate the discovery and investigation of additional environmentally regulated transcriptional states across plant tissues and species. As such, it provides a framework for understanding how environmental regulation shapes cell identity across plant tissues and how plants dynamically remodel cellular states in response to their environment, which will ultimately advance our understanding of plant physiology.

Plant material and growth conditions

Seeds from wild-type A. thaliana ecotypes (Col-0, Ws-2, Cvi-0, c24 and Ler) or the mutants aba2-3 (Col-0 background, ABRC stock number CS3834), pyl11458 (Col-0 background, ABRC stock number CS72390) and scr-3 (Col-0 background, NASC stock code N3997) were sterilized and densely sown in rows on 1× (or 0.5× when testing the effect of nutrient concentration) MS basal medium with vitamins (Austratec) on 0.9% agar plates supplemented with 1% sucrose (pH 5.7) and stratified in the dark at 4 °C for 2 days. Plates were subsequently grown in long day conditions (16 h of light, 8 h of dark) at 22 °C. At 6–7 days after germination (d.a.g.), root tips (up to around 3 mm) were collected for protoplast isolation using sharp tweezers. For all experiments, seedlings were grown under identical conditions and collected at the same time of day. For all comparative experiments (MS concentrations or water-deficit), seedlings from the different conditions were sterilized, plated, grown and harvested simultaneously.

Hard-agar model of drought

For the hard-agar water-deficit assay, seedlings were grown on 0.5× MS medium under control conditions (100% water content) or on 0.5× MS medium with the water content reduced to 50% or 20% of the control level (corresponding to 1× MS with 1.8% agar and 2.5× MS with 4.5% agar, respectively). Seedlings were otherwise grown and collected under the same conditions as described above.

Roots growing through agar

In standard plate assays, roots are grown along the agar surface and are simultaneously exposed to two microenvironments: air on one side and agar on the other. To test whether this heterogeneous exposure influences the formation of the ERS described in this study, roots were instead grown entirely in the agar to provide a homogeneous environment. To do so, square plates were filled with 1× MS medium as described above. After solidification, the medium was cut in half with a razor blade under sterile conditions. Seeds were placed at the top edge of the cut agar and plates were positioned vertically in the growth chamber, allowing cotyledons to develop at the surface while roots penetrated and grew through the agar. Root tips were collected with fine tweezers at the same developmental stage and time point as described above.

Protoplast isolation and scRNA-seq

Root tips were placed in 15 ml protoplasting solution⁴⁵ (1.25% w/v cellulase, 0.3% w/v macerozyme, 0.4 M d-mannitol, 20 mM MES (pH 5.7), 20 mM KCl, 10 mM CaCl₂, 0.5% w/v BSA and 5 mM β-mercaptoethanol in ultrapure water) and incubated at room temperature for 1 h on a tube rotator. The solution was subsequently filtered through a 70 μm strainer and a 40 μm cell strainer (Falcon) to separate protoplasts from undissociated tissue. Cells were centrifuged for 5 min at 500g in a swinging bucket centrifuge. The pellet was resuspended in 15 ml protoplast resuspension buffer (0.4 M d-mannitol, 20 mM MES (pH 5.7), 20 mM KCl, 10 mM CaCl₂ and 0.5% w/v BSA in ultrapure water) and centrifuged at 500g for 5 min. Cells were resuspended in the protoplast resuspension buffer and filtered through a 35 μm cell strainer Snap Cap (Falcon). When preparations involved two Arabidopsis ecotypes, protoplasts from each background were mixed in equal proportions. The final concentration was adjusted to 900–1,200 cells per μl and barcoded using a Chromium Controller (10x Genomics). Libraries were generated with reagents from a Single Cell Gene Expression kit (10x Genomics, Single Cell 3′ Library kit) according to the manufacturer’s instructions. Libraries were sequenced either on an Illumina NextSeq 500 or a NovaSeq 6000, following recommendations of 10x Genomics.

Preprocessing

All the data preprocessing steps are summarized in Supplementary Fig. 1.

Mapping

Reads were mapped using the 10x Genomics CellRanger pipeline (v.7.1.0) to the TAIR10 Arabidopsis genome with the –expect-cells=10000 parameter.

If the cell preparation contained a mix of two ecotypes, reads were mapped to a version of the Arabidopsis (TAIR10) genome with ambiguous base codes placed at single-nucleotide polymorphism (SNP) positions for the genotypes sequenced, as identified by the Arabidopsis 1001 genomes project⁴⁶.

Low-quality cell filtering

Recent versions of CellRanger identify substantially more cells than earlier releases. This enables the detection of low-UMI root cap cells (which were discarded in earlier releases), but also results in the identification of more cells than expected based on cell loading in the assay. Many of these seem to be low-UMI, low-quality cells, which are probably derived from damaged protoplasts that nonetheless pass the CellRanger threshold. To identify and remove such cells, we developed a filtering approach based on clustering CellRanger-called ‘cell’ barcodes together with barcodes classified as background. The method was applied independently to each library and manually adjusted, which enabled the proportion of filtered cells to vary according to the intrinsic quality of each preparation. The procedure (illustrated in Supplementary Fig. 1) is described below.

A Seurat (v.5.0.0)⁴⁷ object was created from the N barcodes identified as cells by CellRanger, excluding genes expressed in fewer than three cells. Data were log-normalized and scaled. Highly variable genes (HVGs) were identified using FindVariableFeatures (nfeatures=15,000). Genes affected by protoplast isolation (see the section ‘Identification of genes affected by protoplast isolation’), as well as mitochondrial and chloroplastic genes, were excluded from the HVG set. The first 40 principal components were used for UMAP dimensionality reduction and clustering. Major root cell types (root cap, epidermis, cortex, endodermis, pericycle, phloem, xylem and procambium) were annotated.

The following steps were applied to each cell type individually:

• A new Seurat object was created containing the ‘cell type’ cells and the first top N barcodes (ranked by UMI content) that were classified as background by CellRanger. The object was processed as above, with the exception that FindVariableFeatures was run with the parameter nfeatures=5000 and that the first 20 principal components were used for UMAP and cell clustering.

• The object was subsetted to retain only the cluster (or clusters) containing the ‘cell type’ cells, which typically also includes some background barcodes. This subsetted object was then re-processed as in the previous step.

• In this subset, background barcodes typically form distinct clusters, often together with low-quality cells of the given type. All cells clustering with background were considered more similar to background than to high-quality cells and were reclassified as background and discarded.

• Cells retained after filtering across all types were merged into a new Seurat object.

Using this approach, between 8 and 30% of CellRanger-identified cells were typically removed. Preparations processed using 10x Genomics NextGem chemistry tended to generate a higher proportion of discarded cells.

Doublet identification

Doublets were identified using scDblFinder^48,49 following the recommendations of the Single-cell Best Practices guide (www.sc-best-practices.org). The analysis was performed on the Seurat object obtained after cell filtering, with the parameters artificialDoublets=25000, dims=30, dbr.sd=1 and propRandom=1. All cells identified as doublets were subsequently removed from the dataset.

Cell genotyping

After doublet removal, preparations composed of a mix of two Arabidopsis ecotypes were subjected to cell genotyping based on SNPs. SNP UMIs were counted per cell using the countsnps command in the sctools package (https://github.com/timoast/sctools). To separate the two ecotype populations, density-based clustering was performed on the SNP UMI counts for each ecotype using the dbscan⁵⁰ package in R (Supplementary Fig. 1), with the eps and minPts parameters adjusted for each dataset. The few cells that did not cluster with either population were discarded.

Data integration

Each condition included at least three replicates, except for the water-deficit experiments, which had two replicates each. Replicates were integrated using Seurat (v.5.0.0) to generate batch-corrected atlases. In brief, replicates were merged into a single Seurat object. HVGs were detected using the FindVariableFeatures function with nfeatures=25000. Genes affected by protoplasting, as well as mitochondrial and chloroplastic genes, were excluded from the HVG dataset. Data were scaled, and the first 50 principal components were computed using RunPCA. Replicates were then integrated using the anchor-based RPCA method in Seurat (v.5.0.0) with IntegrateLayers, using method=RPCAIntegration and dims=1:50.

For the Col-0 reference atlas, the replicates RL650, RL970 and RL971, which exhibited the highest data complexity, were used as reference during integration with reference=1:3. The resulting 50-dimensional integrated RPCA reduction was used for UMAP and clustering.

In some cases, after the integration step, a few diffuse groups of low-UMI cells could be identified associated with the epidermis and cortex lineages, which probably represent low-quality cells that had passed the initial filtering step. To remove these, separate ‘epidermis–root cap’ and ‘cortex’ objects were created, re-integrated and processed as described above. After removing the low-quality cells, the atlas was re-processed with RPCA integration, and the first 45 dimensions were used for UMAP and clustering.

In other cases beyond the Col-0 atlas, small clusters of low-quality cells were also removed; see the integration scripts for further details.

Identification of genes affected by protoplast isolation

Plants were prepared as for scRNA-seq, as described above. Four replicates each of dissociated and non-dissociated root tips were compared. For dissociated samples, root tips were placed in the protoplasting solution for 1 h at room temperature, as for scRNA-seq. After this incubation step, samples were centrifuged for 5 min at 500g, the supernatant was removed and the tissue samples were snap-frozen in liquid nitrogen and ground into a powder. For non-dissociated samples, root tips were collected and directly snap-frozen and ground into powder in liquid nitrogen. Bulk RNA-seq libraries were generated using a TruSeq stranded mRNA kit (Illumina). Bulk RNA-seq libraries were sequenced on an Illumina HiSeq 1500 (single end, 101 cycles). RNA-seq read quality was assessed using FastQC and MultiQC^51,52. RNA-seq reads were trimmed using Trimmomatic (v.0.36) with the following parameters: ILLUMINACLIP:TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 (ref. ⁵³). Reads were aligned to the TAIR10 genome using STAR⁵⁴ with the following parameters: --alignIntronMax 5000, --alignIntronMin 10, --readFilesCommand zcat, --quantMode GeneCounts, --outSAMtype BAM SortedByCoordinate. Differentially expressed genes resulting from the cell dissociation step were identified using the DESeq2 R package⁵⁵, with a Benjamini–Hochberg adjusted P < 0.001 and a log₂ fold change in expression of >2. This experiment was repeated twice, and genes that were identified as differentially expressed in both experiments were taken as the set of genes affected by protoplast isolation.

Whole-mount in situ HCR

Probe design

Candidate cluster markers for our reference atlas were generated using the FindAllMarkers function in Seurat, with the only.pos = TRUE parameter. To select the best target genes, cluster markers were filtered on the basis of the following criteria: the gene should be detected almost only in the cluster (or clusters) it is enriched in; and the gene should have a large portion of its coding sequence that does not share homology with other Arabidopsis genes, which will enable the design of highly specific probes for these sequences. Split probes for our target genes were designed and synthesized by Molecular Instruments (www.molecularinstruments.com).

Sample fixation and dehydration

Experiments were performed as previously described⁵⁶. Roots were collected at 6 d.a.g., as for scRNA-seq, and immediately placed in FAA (7.4% formaldehyde, 5% acetic acid and 50% ethanol in ultrapure water) for 1 h at room temperature. Samples were then dehydrated in successive 10 min washes at room temperature: once in 70% ethanol, once in 90% ethanol and twice in 100% ethanol. Samples were washed twice in 100% methanol and left in the second wash at −20 °C overnight (or up to a few weeks).

Permeabilization and protease digestion

Roots were rehydrated in a series of 5 min washes at room temperature in 75%, 50%, 25% methanol and in DPBS-T (0.1% Tween 20 in DPBS) and incubated in 1× cell wall digestion solution (10× stock: 250 mg macerozyme, 250 mg cellulase, 125 mg pectolyase and 375 mg pectinase in 50 ml ultrapure water) in DPBS-T for 5 min. After a wash in DPBS-T, samples were fixed in 10% (v/v) formaldehyde (37% solution, Sigma-Aldrich) in DPBS-T for 15 min at room temperature and washed twice in DPBS-T. Roots were incubated in 1 ml proteinase K buffer (0.1 M Tris-HCl (pH 8) and 50 mM EDTA (pH 8) in ultrapure water) with 8 μl proteinase K (10 mg ml^–1) for 15 min at 37 °C. After two washes in DPBS-T, samples were fixed in 10% (v/v) formaldehyde (37% solution, Sigma-Aldrich) in DPBS-T for 15 min at room temperature and washed twice in DPBS-T.

Probe hybridization and amplification

We followed the third-generation in situ HCR procedure developed by Molecular Instruments, and as previously described¹⁰, but the dextran sulfate used in the amplification buffer had to be replaced by a low-molecular-weight dextran sulfate to preserve the root structure. Roots were pre-hybridized in the 30% probe hybridization buffer for 30 min at 37 °C, then placed in the probe solution (2–3 μl of probe mixture in 500 μl of 30% probe hybridization buffer) and incubated at 37 °C overnight.

Samples were washed twice in the 30% probe wash buffer for 30 min at 37 °C, then washed twice in 5× SSCT (5× SSC and 0.1% Tween 20 in ultrapure water) for 10 min at room temperature. Samples were placed in the amplification buffer for 30 min at room temperature. In the meantime, fluorescently labelled hairpins were prepared by snap-cooling 10 μl of 3 μM stock in hairpin storage buffer (heat at 95 °C for 90 s in a thermocycler and cool to room temperature in the dark for 30 min). We predominantly used HCR amplifiers coupled to Alexa 647 fluorophores in this study, with Alexa 488-coupled amplifiers used for experiments targeting two genes simultaneously. The hairpin solution was prepared by adding the snap-cooled hairpins to 500 μl amplification buffer. Samples were placed in the hairpin solution and incubated overnight at room temperature in the dark. They were subsequently washed 3 times in 5× SSCT for 20 min at room temperature.

Clearing, cell wall staining and sample mounting

Samples were incubated in ClearSee¹¹ in the dark for 3 days.

Cell walls were stained by incubating samples in 1 ml ClearSee and 2 μl calcofluor solution (50 mg ml^–1 in DMSO) for 2 h in the dark. Samples were washed in ClearSee for 1 h, in the dark, and mounted in ClearSee in a chamber made with a square coverslip and double-sided tape. The chamber was sealed using top coat nail polish, and samples were imaged on the same day.

Imaging

Mounted roots were imaged using a Nikon Ti-E inverted motorized microscope with a Nikon A1Si spectral detector confocal system or a Nikon A1RMP confocal microscope, controlled by NIS-Elements software, and a ×60 water-immersion objective. Calcofluor white and the Alex 488- and the Alexa 647-coupled amplifiers were excited at 405 nm, 488 nm and 640 nm, respectively. Seven to eight images were taken per root along the differentiation axis. Images were processed using Fiji⁵⁷. Images of the same roots were stitched using the Grid/Collection stitching plugin⁵⁸, and the brightness and contrast were adjusted for each image based on the signal detected. Optical sections were generated using the Orthogonal Views method in Fiji.

Cell annotation

All datasets other than the Col-0 atlas were annotated using the Col-0 atlas as the reference. Anchors between the reference atlas and the other datasets were identified using FindTransferAnchors in Seurat (v.5), specifying the Col-0 RPCA-integrated data as the reference assay. TransferData in Seurat (v.5) was then used to calculate, for each cell, a prediction score for every identity, based on the identified anchors and the atlas annotations. Each cell was assigned the identity with the highest prediction score.

Cell cycle annotation

We used the CellCycleScoring function in Seurat to assign a cell cycle score to each cell of our reference atlas based on the expression of S phase marker genes (E2Fc (also known as AT1G47870), E2Fb (AT5G22220), DPa (also known as AT5G02470), CDKG;2 (also known as AT1G67580), E2Fa (also known as AT2G36010), DPb (also known as AT5G03415), CYCA3;1 (also known as AT5G43080), CYCA3;2 (also known as AT1G47210), KRP5 (also known as AT3G24810), KRP2 (also known as AT3G50630), KRP6 (also known as AT3G19150), KRP3 (also known as AT5G48820), WEE (also known as AT1G02970), KRP4 (also known as AT2G32710), CYCD4;2 (also known as AT5G10440)) and of G2/M phase marker genes (CYCB1;1 (also known as AT4G37490), CYCB2;1 (also known as AT2G17620), CDKB2;2 (also known as AT1G20930), CYCB3;1 (also known as AT1G16330), CYCA2;2 (also known as AT5G11300), CYCA2;4 (also known as AT1G80370), CYCA1;1 (also known as AT1G44110), CYCA2;1 (also known as AT5G25380), CDKD;1 (also known as AT1G73690), CDKB2;1 (also known as AT1G76540), CKS2 (also known as AT2G27970), CYCB1;2 (also known as AT5G06150), CYCD3;1 (also known as AT4G34160), CCS52B (also known as AT5G13840), CYCA2;3 (also known as AT1G15570), CCS52A2 (also known as AT4G11920), CYCB2;3 (also known as AT1G20610) and CYCB1;3 (also known as AT3G11520)) obtained from a previous study⁵⁹. This score was used to provide a predicted classification of the cell cycle phase and to regress out cell cycle effects using vars.to.regress = c(“S.Score”, “G2M.Score”) in the ScaleData function. After regression, the top 50 principal components were calculated and the replicates were integrated again. UMAP dimensionality reduction was then performed, and cell identities were visualized in the embedding to confirm that both branches remained present.

Pseudotime trajectory analyses

PAGA

The Col-0 atlas was subsetted into nine lineage-specific objects: epidermis, root cap, cortex, endodermis, PPP, XPP, procambium, phloem and xylem. Each lineage dataset was then processed independently.

Subclusters in each lineage were identified using FindNeighbors and FindClusters in Seurat (v.5), with the values of dims, k.param and resolution adapted to each lineage.

For each lineage, a ‘root cell’ corresponding to t = 0 in the trajectory was defined, as required by the diffusion pseudotime analysis underlying the PAGA trajectory inference. In the case of the epidermis and root cap lineages, developmental progression is known to originate at their junction. Therefore, for the epidermal lineage, the root cell was defined as the cell with the lowest coordinate along the UMAP y axis. For the root cap lineage, the root cell was defined as the cell with the highest coordinate along the UMAP y axis. For all other lineages, for which such anatomical landmarks are less obvious, the root cell was defined as the cell with the highest meristematic score (see below).

To perform trajectory inference, we used the SCP package⁶⁰, which implements PAGA^61,62,63 analysis directly in R, thereby avoiding conversion of Seurat objects to AnnData format to run scanpy⁶³. Seurat v.5 assays were first converted to v.3 assays to enable SCP. Trajectories were inferred with the RunPAGA function, grouping cells by the subclusters defined above, with the following parameters: inear_reduction = “integrated.rpca”, nonlinear_reduction = “umap.rpca”, infer_pseudotime=TRUE, n_pcs and n_neighbours matching the values used for subclustering, and root_cell containing the barcode of the designated root cell. The trajectory topology was visualized with the PAGAPlot function, using reduction = “umap.rpca” and an edge_threshold adapted to each lineage. The pseudotime value assigned to each cell was then transferred to the metadata of the original Seurat (v.5) object.

Meristematic score calculation

The diffusion pseudotime analysis required the identification of a cell that would serve as t = 0; that is, the starting point of each trajectory. To do so, we developed an approach for identifying the ‘most meristematic’ cell in each lineage.

The epidermis is the only lineage for which the meristematic region can be unambiguously located based on both the topology of its clusters and the known biology of the tissue. It lies at the junction with the root cap, before differentiation into the two epidermal subtypes: hair cells and non-hair cells. This cluster was therefore designated as epidermis–meristem. To identify marker genes for this region, we ran the FindMarkers function on the subsetted Seurat object containing only epidermal clusters, with ident.1 = 5 (corresponding to the ‘epidermis–meristem’ cluster), only.pos=TRUE, and test.use = ‘t’. The resulting gene list was then used with the AddModuleScore function in Seurat (v.5) to assign a meristematic score to every cell in the atlas.

We reasoned that at least a subset of the genes enriched in the epidermal meristem would also be enriched in meristematic regions of other lineages. Thus, the cell with the highest meristematic score in each lineage was taken to represent the most likely earliest developmental state and used as the root cell (t = 0) for pseudotime analysis.

Monocle analysis

Monocle 3 (v.1.3.7)^64,65,66 was used as an alternative method to confirm the developmental trajectories obtained using PAGA. Our Col-0 atlas was converted to a cds object using SeuratWrappers (v.0.4.0) (https://github.com/satijalab/seurat-wrappers). Cells were clustered and partitioned with cluster_cells using resolution=0.0002, and trajectories were computed using learn_graph with learn_graph_control=list(minimal_branch_len=20)). As for the PAGA trajectories, root cells were defined as those with the highest meristematic score per lineage, and these were used in order_cells to assign a pseudotime to each cell in the trajectory. All plots were generated with plot_cells. Concordance between (scaled) Monocle and diffusion pseudotime assignments was assessed using Spearman’s correlation coefficients computed independently for each lineage.

scVelo analysis

RNA velocity analysis was performed using the dynamical model implemented in the R package velociraptor⁶⁷ (v.1.10.0) (wrapping scVelo⁶⁸) on a per-lineage basis. Spliced and unspliced count matrices were generated independently for each library from velocyto⁶⁹ (v.0.17.17) loom files and merged after restricting to genes detected across all samples. Because the intronic signal in Arabidopsis scRNA-seq datasets is limited, only genes with sufficient spliced and unspliced coverage were retained for velocity inference (minimum 50 cells with unspliced counts and 300 total unspliced counts; minimum 50 cells with spliced counts and 1,000 total spliced counts globally across the dataset). The top 4,000 variable genes were identified independently for each lineage, and protoplasting-sensitive genes, mitochondrial genes and chloroplast genes were excluded before analysis. The final gene set used for velocity inference corresponded to the intersection between lineage variable features and genes passing global spliced or unspliced filtering criteria. For each lineage, a SingleCellExperiment object containing spliced, unspliced and total RNA counts was generated and analysed using the scVelo dynamical model with the RPCA integrated embedding used as the neighbourhood graph for kinetic inference. Root cells for latent time estimation were defined independently for each lineage as the top 1% of cells with the highest meristematic score. Given the relatively low proportion of unspliced reads in our dataset (about 3% median unspliced fraction), RNA velocity was used primarily as an orthogonal validation of trajectory structure rather than as the primary trajectory inference approach. Concordance between scVelo latent time and diffusion pseudotime was assessed using Spearman’s correlation coefficients computed independently for each lineage.

Expression-level analysis of ERS-associated genes

To test whether the apparent absence of ERS-associated genes in ERS⁻ cells could be explained by reduced transcript capture, we compared the expression behaviour of genes from ERS⁺-specific modules (modules 20 and 21) and modules expressed in both branches (modules 13 and 14) in the XPP lineage, the largest lineage exhibiting the ERS dichotomy. Average log-normalized expression was calculated independently in ERS⁺ and ERS⁻ cells for each gene. To compare genes across a shared expression range, overlapping expression intervals between ERS⁺-specific and shared-module genes were defined using the 10–90th percentile range of average ERS⁺ expression for each group. A linear model describing the relationship between ERS⁺ and ERS⁻ expression was then fitted using genes from modules expressed in both branches. Residual ERS⁻ expression values (observed minus expected expression) were subsequently calculated for all genes relative to this model. Residual distributions were compared between ERS⁺-specific genes and genes expressed in both branches in the overlapping expression interval using a two-sided Wilcoxon rank-sum test. This analysis tested whether ERS⁺-specific genes exhibited systematically lower ERS⁻ expression than expected based solely on their expression level in ERS⁺ cells.

UMI downsampling

Because the number of transcripts detected in branch 1 was systematically higher than in branch 2, we tested whether the observed branching events could result from coverage artefacts. To control for this factor, UMI counts were selectively downsampled in ERS⁺ cells to match the UMI distribution of ERS⁻ cells in each cell type. Specifically, for each ERS⁺ cell, a target UMI count was assigned based on quantile matching to the corresponding ERS⁻ cell-type distribution, and counts were downsampled using the downsampleMatrix function from DropletUtils⁷⁰ (v1.20.0). ERS⁻ cells were left unchanged. Three independent downsampled count matrices were generated using different random seeds. Each downsampled matrix was used to generate a new Seurat object, which was processed independently using the same workflow as the reference atlas, including log-normalization, variable feature selection, scaling, PCA and UMAP embedding. Protoplasting-sensitive genes, mitochondrial genes and chloroplast genes were excluded from the variable feature set before dimensionality reduction. A meristematic score was then calculated for each cell using the same gene set as in the reference atlas.

Identification of gene co-expression modules

We used hdWGCNA (v0.4.1)^71,72,73 to identify modules of co-expressed genes across our Col-0 reference atlas. To determine the list of genes that would be used in the analysis, we first applied the SetupForWGCNA function to obtain all genes expressed in at least 10% of cells in at least one lineage. This list was intersected with the VariableFeatures of the atlas, which contains all genes detected in every replicate, excluding those affected by protoplast isolation as well as mitochondrial and chloroplastic genes. The resulting set of 12,923 genes was supplied to SetupForWGCNA with the parameter gene_select = “custom”.

To reduce data sparsity before WGCNA, we generated metacells (aggregates of small groups of similar cells per biological sample) and their associated expression matrices, which served as inputs for the network analysis. Metacells were generated using the MetacellsByGroups function using group.by=c(“lineage”,“orig.ident”), reduction = “integrated.rpca”, ident.group = “lineage”, k = 20 and max_shared=5, followed by normalization with NormalizeMetacells.

We then set up the expression dataset for consensus network analysis across all replicates using SetMultiExpr with multi.group.by = “orig.ident”. After determining the soft power threshold for each batch with the TestSoftPowersConsensus function, the consensus co-expression network was constructed using ConstructNetworkusing minModuleSize=30, consensus=TRUE and deepSplit=2.

Module scores were computed using ModuleExprScore using n_genes = “all” and method = “Seurat”, and projected onto the atlas UMAP with ModuleFeaturePlot using features = “scores”.

GO enrichment analysis for ERS-specific modules

GO enrichment analysis was performed on the gene sets from modules 20 and 21 using ShinyGO (v.0.82; https://bioinformatics.sdstate.edu/go/). The analysis was run with the GO Biological Process pathway database, an FDR cut-off value of 0.01 and the number of pathways to display set to 50. Pathways were ranked by fold enrichment, and the results were visualized as a hierarchical clustering tree.

TF regulon analysis

Putative TF target genes were inferred from the DAP-seq dataset¹⁷. The ‘Motifs in peaks (P ≤ 1 × 10^–4)’ BED file (http://neomorph.salk.edu/dap_web/pages/browse_table_aj.php) was used as a DNA-binding database for the profiled TFs of the DAP-seq dataset.

The 500 bp promoter coordinates for Arabidopsis genes were extracted using the promoters function in GenomicRanges⁷⁴, genes(TxDb.Athaliana.BioMart.plantsmart28) as an annotation for Arabidopsis genes and the arguments upstream = 500, downstream = 50. For each TF, the binding coordinates were extracted from the Motif in Peaks.bed file and intersected with the promoter coordinates using the subsetByOverlaps function in GenomicRanges. Every gene that had a binding motif for the TF in the 500 bp region upstream of its transcription start site (TSS) was considered a putative target gene. For TFs that were not profiled in the DAP-seq studies but, according to the cisBP2 study¹⁸, have a strong similarity in their DNA-binding domains with another TF that was studied by DAP-seq, their target genes were considered to be the same as the TF profiled by DAP-seq.

We then converted the TF target predictions into data to be processed with hdWGCNA. A sparse matrix was generated whereby rows are all genes in the dataset, columns are TFs with putative targets and values are set to 1 if the gene was a target of the TF and 0 otherwise. This was input into the Seurat object using the SetMotifMatrix function of hdWGCNA (v.0.4.1). The list of TFs and their target genes was set with SetMotifTargets, and a dataframe containing the target numbers was input with SetMotifs. Because we were interested in what drives the expression of co-expression modules in the dataset, we restricted the analysis to all genes that were assigned to a co-expression module as well as all TFs with target information. A network of TFs and target genes was built using ConstructTFNetwork, in which the XGBoost regression was used to evaluate how important each TF is for predicting each gene based on single-cell data^72,75. The network was then pruned to keep only the strongest TF–gene connections. Specifically, we selected the top ten TFs above a regulatory score threshold as regulators for each gene, using the AssignTFRegulons function with strategy = “A”, reg_thresh=0.01 and n_tfs=10. A positive regulon score was then computed for each cell using RegulonScores, and the scaled average regulon score per cell identity was plotted using pheatmap⁷⁶.

Analysis of published single-cell and single-nucleus datasets

Three published Arabidopsis root single-cell datasets were used to verify that the ERS and associated ERS-specific genes described in this study are also present in other datasets^7,8,16. For ref. ⁸, we directly used the processed Seurat object provided by the authors. For ref. ⁷, we reprocessed the raw fastq files using the same pre-processing steps as described for our Col-0 atlas. For ref. ¹⁶, we used the count matrices provided by the authors and performed batch integration.

Once the final Seurat objects were generated, cells were annotated using our reference atlas following the same procedure described in the section ‘Cell annotation’ above. Cells were also projected onto the atlas UMAP embedding using the MapQuery function in Seurat (v.5.0.0). Co-expression gene modules were projected onto these datasets using the ProjectModules function in hdWGCNA (v.0.4.1), and module scores were calculated with ModuleExprScore using n_genes = “all” and method = “Seurat”.

Label-transfer stability analysis

To assess the robustness of ERS assignments in external datasets, label transfer was repeated across a range of Seurat parameter settings using the integrated Col-0 atlas as reference. Specifically, projections were performed using 30, 40 or 50 principal components and k.anchor values of 5, 10 or 20, which generated nine parameter combinations in total. For each parameter set, transfer anchors were identified using FindTransferAnchors with the RPCA integrated embedding, and cell identities were assigned using TransferData. To focus analyses on cells robustly assigned to the ERS bifurcation, only cells classified as ERS⁺ or ERS⁻ in the baseline parameter set (50 principal components, k.anchor = 5) were retained for downstream comparisons. Predicted identities were separated into cell-type and ERS-state components, and stability was assessed independently for both. Cell-type stability was quantified as the fraction of parameter sets in which each cell retained its baseline cell-type assignment, whereas ERS stability was calculated as the consistency of ERS⁺ and ERS⁻ assignment in the baseline cell type across parameter sets. Global agreement between parameter sets was further quantified using pairwise adjusted rand index (ARI) comparisons across all transferred labels. Furthermore, the proportion of ERS⁺ cells in each cell type, the fraction of cells retaining their baseline cell type assignment and the frequency of ERS state switching across parameter sets were calculated independently for each external dataset.

Analysis of a previously published snATAC-seq dataset

The fragment files generated with Cell Ranger ATAC (v.1.2) and mapped to TAIR10/Araport11 reference were obtained from a previous study¹⁶, imported to snapatac2 (v.2.6.4) and paired-insertions calculated⁷⁷ in 500 bp bins⁷⁸. The cells were initially defined as barcodes associated with more than 3,500 and fewer than 50,000 fragments per barcode and with a TSS enrichment score higher than 2.5. Low-quality cells and outliers were also identified per-sample basis using the 3× median absolute deviation threshold and further removed from analysis⁷⁹. Following peak calling using macs3 (ref. ⁸⁰), peak and gene activity matrices were calculated using the paired-insertion counting strategy. For ATAC–seq data integration into our scRNA-seq dataset, the peak matrices were imported to the Signac package⁸¹ and the dimensions of ATAC–seq were reduced using latent semantic indexing. ATAC–seq nuclei were assigned identities using our reference atlas and the FindTransferAnchors function with the parameters reduction = ‘cca’ and dims = 1:30, followed by TransferData with weight.reduction = obj[[“lsi”]] and dims = 2:30. The TransferData function and previously computed anchors were also used to generate an ‘imputed’ scRNA-seq matrix for the ATAC nuclei. Using this imputed matrix, the snATAC–seq and scRNA-seq reference atlases were merged into a single object, processed with ScaleData and RunPCA, and the first 30 principal components were used to produce a UMAP co-embedding of the two datasets.

Differential module score analysis

To perform differential expression analysis of our co-expression modules, a Seurat object was created with all conditions to compare, and the data were scaled using ScaleData (Seurat v.5.0.0). The modules were projected into the new dataset using ProjectModules (hdWGCNA v.0.4.1). A module score was then calculated using ModuleExprScore. The differential module scores were determined using the FindDMEs function with features = “ModuleScores” and test.use = “wilcox”. The results were plotted using PlotDMEsLollipop.

Differential composition analysis

All changes in cell composition were tested with sccomp (v.1.99.18) in R⁸². Specifically, when evaluating changes in ERS⁺ (and consequently ERS⁻) proportions, the Seurat object was subsetted to retain only cells that had reached the developmental stage at which they could be assigned to either state to avoid taking into account sampling effects that might affect the total number of mature cells captured. This process excluded the root cap, ‘early’ identities and the non-hair cells, for which the trajectories could determine the point of branching. The test was then run on the subsetted data using sccomp_estimate and sccomp_test, and results were plotted with sccomp_boxplot or with plot in the case of the water-deficit experiment for which the covariate is continuous.

Counting of hair cells

Seedlings were grown vertically on full (1×) or half (0.5×) MS plates, as for scRNA-seq experiments. Pictures of the roots were taken at 6 d.a.g., using a stereomicroscope (n = 16 for full MS, n = 23 for half MS, from two different plates each). A portion corresponding to approximately 2.5 mm from the tip was selected for each image of the root, and root hairs were manually counted using the Cell Counter plugin in Fiji. A significant difference in root hair density between the two conditions was determined using a Wilcoxon rank-sum test (wilcox.test function) in R.

Note that the seedlings used for this quantification were in the top row of the plate (seeds were sown in two rows for scRNA-seq experiments). Roots of seedlings grown in the bottom row tended to have fewer root hairs (in both conditions).

Ecotype-specific DEGs

A new Seurat object was created using the Col-0 cells that were grown on full MS (which is all replicates except the three half MS replicates), as well as the cells of the C24, Cvi-0, Ler and Ws-2 atlases. The expression data were then pseudobulked per genotype, cell identity and batch using the AggregateExpression in Seurat (v.5.0.0). For each possible inter-accession replicate pair, a differential expression test was performed on the pseudobulk data using FindMarkers with DESeq2. Genes were considered differentially expressed if the adjusted P < 0.05 and the absolute log₂ fold-change was >0.5. For a gene to be called an accession-specific DEG, it had to be systematically differentially expressed in a given accession compared with the other four, but not be differentially expressed in any pairwise comparison between the remaining four accessions. The proportion of ecotype-specific DEGs per identity or per module were plotted using pheatmap.

Adjustment for transcript abundance in water-deficit analyses

Module scores for hdWGCNA modules 20 and 21 were calculated using AddModuleScore in Seurat (v.5.0.0). Analyses were restricted to cells classified as ERS⁺ or ERS⁻ from the 100% and 20% water conditions. Total UMI counts were log₁₀-transformed after addition of a pseudocount. Module scores were analysed using linear mixed-effects models implemented in the lme4 (v.1.1.36) and lmerTest (v.3.2.1) R packages, with water condition as a fixed effect, log₁₀ UMI count modelled using a natural spline with four degrees of freedom (splines v.4.3.1 package) and biological replicate included as a random intercept. For visualization, cells were divided into ten equal-frequency bins according to log₁₀ UMI count, and mean module scores ± s.e.m. were calculated in each bin. Binning was used for visualization only; all statistical inference was performed on the unbinned cell-level data. The association between water condition and ERS identity was assessed using binomial generalized linear mixed-effects models (glmer; lme4) with biological replicate included as a random intercept. An unadjusted model included water condition alone, whereas the adjusted model also included log₁₀ UMI count as a fixed-effect covariate. Odds ratios and 95% confidence intervals for the effect of water condition were obtained using broom.mixed (v.0.2.9.6).

Quantification of HCR signals in the water-deficit experiment

To compare roots grown under different conditions using whole-mount HCR, all roots were grown, processed and imaged simultaneously using identical confocal settings.

First, a mask defining the cytoplasmic area of the root was created in Fiji. A threshold-based mask of the whole root was generated from the calcofluor channel using Image > Adjust > Threshold, followed by Process > Binary > Close and Process > Binary > Fill Holes. A second threshold-based mask was created from the calcofluor channel to define cell walls. Subtracting the cell wall mask from the whole-root mask (Process > Math > Subtract) produced a mask containing only the cytoplasmic regions. This mask was then multiplied by a mask selecting the region from the elongation zone upwards, corresponding to the area where ERS expression begins in single-cell data.

The resulting cytoplasmic mask was multiplied by the channel containing the HCR signal to retain only signal in cells from the elongation zone upwards. A separate threshold-based mask was created from the HCR signal channel to identify areas above a defined signal threshold. Multiplying this mask by the cytoplasmic mask ensured that only signal inside cells from the elongation zone upwards was considered.

Custom Fiji macros were used to automate measurements across z-stacks. One macro measured the surface area covered by a mask in each slice, which enabled calculation of the total cytoplasmic area for normalization and the area of positive signal. Another macro measured the sum of all pixel intensities per slice, which was applied to the cytoplasmic signal image.

Differences in signal between conditions were assessed using the Wilcoxon rank-sum test.

Cell segmentation and volumetric quantification of single-cell fluorescence

Single-cell segmentation was performed in 3D using the cell-wall staining channel. Individual cells were first manually seeded by placing markers in each cell lumen and storing these positions as regions of interest. A binary marker image was then generated from the seed positions. The cell wall channel was processed using the MorphoLibJ⁸³ Morphological Filters (3D) function with the Gradient filter to enhance cell boundaries. The processed cell-wall image and binary marker image were subsequently used for marker-controlled watershed segmentation using MorphoLibJ. Segmentation outputs were manually curated to correct improper segmentations caused by discontinuities or holes in cell-wall staining, either by removing incorrectly segmented cells or merging oversegmented regions. Before fluorescence quantification, signal channels were background-corrected to minimize the contribution of imaging noise and diffuse background fluorescence. In brief, pixel intensities were measured in a region outside the root tissue, and the 99th percentile background intensity value was subtracted from all pixels in the image using the Subtract function in Fiji. The resulting filtered fluorescence images were combined with the curated segmentation masks to quantify volumetric single-cell mean fluorescence intensity using the Intensity Measurements 2D/3D function in MorphoLibJ. Measured fluorescence values were then mapped back onto segmented cells using the Assign Measure to Labels function in MorphoLibJ to generate colour-coded single-cell fluorescence representations.

For the water-deficit HCR experiments, segmentation was restricted to epidermal and cortical cells. Fluorescence quantification was performed on the original, non-background-subtracted images so that identical image processing was applied across control and water-deficit samples, which enabled direct comparison of fluorescence intensities between conditions. Mean fluorescence intensity was measured for each segmented cell. To account for batch-to-batch differences in HCR signal intensity, the mean fluorescence intensity for each segmented cell was normalized to the mean per-cell fluorescence of the corresponding 100% water control sample processed in parallel. Empirical cumulative distribution functions were then generated to compare the distributions of normalized single-cell fluorescence intensities between conditions.

Detection of bimodal cell populations in fluorescence imaging and scRNA-seq data

To test whether mutually exclusive marker expression defined two separable cell populations, Gaussian mixture modelling was applied independently to fluorescence quantification data and scRNA-seq expression data using the R package mclust⁸⁴ (v.6.0.0). For fluorescence imaging data, mean single-cell fluorescence intensities for AT3G01720 and AT5G56550 were extracted from segmented cells following background correction and volumetric quantification. A marker ratio was calculated for each cell as AT3G01720/(AT3G01720 + AT5G56550), and Gaussian mixture modelling was performed on the resulting distribution using Mclust. Cells were assigned to the inferred populations based on maximum posterior probability. The fitted Gaussian components and their relative proportions were visualized over the observed ratio distribution. For scRNA-seq data, normalized expression values for AT3G01720 and AT5G56550 were extracted from hair cells in the reference atlas. Cells with zero combined expression for both markers were excluded from the analysis. A marker score was then calculated as the difference between AT3G01720 and AT5G56550 expression values, and Gaussian mixture modelling was applied to the score distribution using Mclust. Fitted Gaussian distributions were visualized together with the observed score distribution to assess the presence of two transcriptionally distinct populations.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

The raw FASTQ files have been deposited into the European Nucleotide Archive under project accession PRJEB100815. Processed Seurat objects have been deposited into Zenodo (https://doi.org/10.5281/zenodo.17411120 and https://doi.org/10.5281/zenodo.17411065)^85,86. Previously published datasets re-analysed in this study are available from the Gene Expression Omnibus under accession numbers GSE155304 (ref. ¹⁶), GSE141730 (ref. ⁷) and GSE152766 (ref. ⁸). Source data are provided with this paper.

1. Petricka, J. J., Winter, C. M. & Benfey, P. N. Control of Arabidopsis root development. Annu. Rev. Plant Biol. 63, 563–590 (2012).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

2. Jean-Baptiste, K. et al. Dynamics of gene expression in single root cells of Arabidopsis thaliana. Plant Cell 31, 993–1011 (2019).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

3. Denyer, T. et al. Spatiotemporal developmental trajectories in the Arabidopsis root revealed using high-throughput single-cell RNA sequencing. Dev. Cell 48, 840–852 (2019).

   Article  CAS  PubMed  Google Scholar 

4. Ryu, K. H., Huang, L., Kang, H. M. & Schiefelbein, J. Single-cell RNA sequencing resolves molecular relationships among individual plant cells. Plant Physiol. 179, 1444–1456 (2019).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

5. Shulse, C. N. et al. High-throughput single-cell transcriptome profiling of plant cell types. Cell Rep. 27, 2241–2247 (2019).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

6. 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, 648–660 (2019).

   Article  CAS  PubMed  Google Scholar 

7. Wendrich, J. R. et al. Vascular transcription factors guide plant epidermal responses to limiting phosphate conditions. Science 370, eaay4970 (2020).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

8. Shahan, R. et al. A single-cell Arabidopsis root atlas reveals developmental trajectories in wild-type and cell identity mutants. Dev. Cell 57, 543–560 (2022).

   Article  CAS  PubMed  PubMed Central  Google Scholar 

9. Otero, S. et al. A root phloem pole cell atlas reveals common transcriptional states in protophloem-adjacent cells. Nat. Plants 8, 954–970 (2022).

   Article  CAS  PubMed  Google Scholar 

10. Choi, H. M. T. et al. Third-generation in situ hybridization chain reaction: multiplexed, quantitative, sensitive, versatile, robust. Development 145, dev165753 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

11. Kurihara, D., Mizuta, Y., Sato, Y. & Higashiyama, T. ClearSee: a rapid optical clearing reagent for whole-plant fluorescence imaging. Development 142, 4168–4179 (2015).

    CAS  PubMed  PubMed Central  Google Scholar 

12. Kumpf, R. P. & Nowack, M. K. The root cap: a short story of life and death. J. Exp. Bot. 66, 5651–5662 (2015).

    Article  CAS  PubMed  Google Scholar 

13. Pi, L. et al. Organizer-derived WOX5 signal maintains root columella stem cells through chromatin-mediated repression of CDF4 expression. Dev. Cell 33, 576–588 (2015).

    Article  CAS  PubMed  Google Scholar 

14. Willemsen, V. et al. The NAC domain transcription factors FEZ and SOMBRERO control the orientation of cell division plane in Arabidopsis root stem cells. Dev. Cell 15, 913–922 (2008).

    Article  CAS  PubMed  Google Scholar 

15. Shi, C.-L. et al. The dynamics of root cap sloughing in Arabidopsis is regulated by peptide signalling. Nat. Plants 4, 596–604 (2018).

    Article  CAS  PubMed  Google Scholar 

16. Farmer, A., Thibivilliers, S., Ryu, K. H., Schiefelbein, J. & Libault, M. Single-nucleus RNA and ATAC sequencing reveals the impact of chromatin accessibility on gene expression in Arabidopsis roots at the single-cell level. Mol. Plant 14, 372–383 (2021).

    Article  CAS  PubMed  Google Scholar 

17. O’Malley, R. C. et al. Cistrome and epicistrome features shape the regulatory DNA landscape. Cell 165, 1280–1292 (2016).

    Article  PubMed  PubMed Central  Google Scholar 

18. Lambert, S. A. et al. Similarity regression predicts evolution of transcription factor sequence specificity. Nat. Genet. 51, 981–989 (2019).

    Article  CAS  PubMed  ADS  Google Scholar 

19. Choi, H., Hong, J., Ha, J., Kang, J. & Kim, S. Y. ABFs, a family of ABA-responsive element binding factors. J. Biol. Chem. 275, 1723–1730 (2000).

    Article  CAS  PubMed  Google Scholar 

20. Uno, Y. et al. Arabidopsis basic leucine zipper transcription factors involved in an abscisic acid-dependent signal transduction pathway under drought and high-salinity conditions. Proc. Natl Acad. Sci. USA 97, 11632–11637 (2000).

    Article  CAS  PubMed  PubMed Central  ADS  Google Scholar 

21. Kang, J.-Y., Choi, H.-I., Im, M.-Y. & Kim, S. Y. Arabidopsis basic leucine zipper proteins that mediate stress-responsive abscisic acid signaling. Plant Cell 14, 343–357 (2002).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

22. Himmelbach, A., Hoffmann, T., Leube, M., Höhener, B. & Grill, E. Homeodomain protein ATHB6 is a target of the protein phosphatase ABI1 and regulates hormone responses in Arabidopsis. EMBO J. 21, 3029–3038 (2002).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

23. Johannesson, H., Wang, Y., Hanson, J. & Engström, P. The Arabidopsis thaliana homeobox gene ATHB5 is a potential regulator of abscisic acid responsiveness in developing seedlings. Plant Mol. Biol. 51, 719–729 (2003).

    Article  CAS  PubMed  Google Scholar 

24. Söderman, E., Mattsson, J. & Engström, P. The Arabidopsis homeobox gene ATHB-7 is induced by water deficit and by abscisic acid. Plant J. 10, 375–381 (1996).

    Article  PubMed  ADS  Google Scholar 

25. Liu, T., Longhurst, A. D., Talavera-Rauh, F., Hokin, S. A. & Barton, M. K. The Arabidopsis transcription factor ABIG1 relays ABA signaled growth inhibition and drought induced senescence. eLife 5, e13768 (2016).

    Article  PubMed  PubMed Central  Google Scholar 

26. Kuno, N. et al. The novel MYB protein EARLY-PHYTOCHROME-RESPONSIVE1 is a component of a slave circadian oscillator in Arabidopsis. Plant Cell 15, 2476–2488 (2003).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

27. Zhang, X. et al. Constitutive expression of CIR1 (RVE2) affects several circadian-regulated processes and seed germination in Arabidopsis. Plant J. 51, 512–525 (2007).

    Article  CAS  PubMed  ADS  Google Scholar 

28. Hsu, P. Y., Devisetty, U. K. & Harmer, S. L. Accurate timekeeping is controlled by a cycling activator in Arabidopsis. eLife 2, e00473 (2013).

    Article  PubMed  PubMed Central  Google Scholar 

29. Leung, J., Merlot, S. & Giraudat, J. The Arabidopsis ABSCISIC ACID-INSENSITIVE2 (ABI2) and ABI1 genes encode homologous protein phosphatases 2C involved in abscisic acid signal transduction. Plant Cell 9, 759–771 (1997).

    CAS  PubMed  PubMed Central  Google Scholar 

30. Wang, X. et al. ABRE-BINDING FACTORS play a role in the feedback regulation of ABA signaling by mediating rapid ABA induction of ABA co-receptor genes. New Phytol. 221, 341–355 (2019).

    Article  CAS  PubMed  Google Scholar 

31. Moriwaki, T., Miyazawa, Y., Fujii, N. & Takahashi, H. Light and abscisic acid signalling are integrated by MIZ1 gene expression and regulate hydrotropic response in roots of Arabidopsis thaliana: light and ABA regulate MIZ1-dependent hydrotropism. Plant Cell Environ. 35, 1359–1368 (2012).

    Article  CAS  PubMed  Google Scholar 

32. Laby, R. J., Kincaid, M. S., Kim, D. & Gibson, S. I. The Arabidopsis sugar-insensitive mutants sis4 and sis5 are defective in abscisic acid synthesis and response. Plant J. 23, 587–596 (2000).

    Article  CAS  PubMed  ADS  Google Scholar 

33. Antoni, R. et al. PYRABACTIN RESISTANCE1-LIKE8 plays an important role for the regulation of abscisic acid signaling in root. Plant Physiol. 161, 931–941 (2013).

    Article  CAS  PubMed  Google Scholar 

34. Fukaki, H. et al. Genetic evidence that the endodermis is essential for shoot gravitropism in Arabidopsis thaliana. Plant J. 14, 425–430 (1998).

    Article  CAS  PubMed  ADS  Google Scholar 

35. Scheres, B. et al. Mutations affecting the radial organisation of the Arabidopsis root display specific defects throughout the embryonic axis. Development 121, 53–62 (1995).

    Article  CAS  Google Scholar 

36. Di Laurenzio, L. et al. The SCARECROW gene regulates an asymmetric cell division that is essential for generating the radial organization of the Arabidopsis root. Cell 86, 423–433 (1996).

    Article  PubMed  Google Scholar 

37. Carlsbecker, A. et al. Cell signalling by microRNA165/6 directs gene dose-dependent root cell fate. Nature 465, 316–321 (2010).

    Article  CAS  PubMed  PubMed Central  ADS  Google Scholar 

38. Gonzalez, S. et al. Arabidopsis transcriptome responses to low water potential using high-throughput plate assays. eLife 12, RP84747 (2024).

    Article  PubMed  PubMed Central  Google Scholar 

39. Seo, P. J. & Mas, P. STRESSing the role of the plant circadian clock. Trends Plant Sci. 20, 230–237 (2015).

    Article  CAS  PubMed  Google Scholar 

40. Araújo, I. S. et al. Stochastic gene expression in Arabidopsis thaliana. Nat. Commun. 8, 2132 (2017).

    Article  PubMed  PubMed Central  ADS  Google Scholar 

41. Alamos, S., Reimer, A., Niyogi, K. K. & Garcia, H. G. Quantitative imaging of RNA polymerase II activity in plants reveals the single-cell basis of tissue-wide transcriptional dynamics. Nat. Plants 7, 1037–1049 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

42. Hani, S. et al. Live single-cell transcriptional dynamics via RNA labelling during the phosphate response in plants. Nat. Plants 7, 1050–1064 (2021).

    Article  CAS  PubMed  Google Scholar 

43. Abley, K. & Locke, J. C. W. Noisy transcription under the spotlight. Nat. Plants 7, 996–997 (2021).

    Article  CAS  PubMed  Google Scholar 

44. Zhao, L., Fonseca, A., Meschichi, A., Sicard, A. & Rosa, S. Whole-mount smFISH allows combining RNA and protein quantification at cellular and subcellular resolution. Nat. Plants 9, 1094–1102 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

45. Bargmann, B. O. R. & Birnbaum, K. D. Fluorescence activated cell sorting of plant protoplasts. J. Vis. Exp. https://doi.org/10.3791/1673 (2010).

46. 1001 Genomes Consortium. 1,135 genomes reveal the global pattern of polymorphism in Arabidopsis thaliana. Cell 166, 481–491 (2016).

    Article  Google Scholar 

47. Hao, Y. et al. Dictionary learning for integrative, multimodal and scalable single-cell analysis. Nat. Biotechnol. 42, 293–304 (2024).

    Article  CAS  PubMed  ADS  Google Scholar 

48. Germain, P.-L. scDblFinder. R package version 1.14.0 https://doi.org/10.18129/B9.BIOC.SCDBLFINDER (2019).

49. Germain, P.-L., Lun, A., Garcia Meixide, C., Macnair, W. & Robinson, M. D. Doublet identification in single-cell sequencing data using scDblFinder. F1000Research 10, 979 (2021).

    Article  PubMed  PubMed Central  Google Scholar 

50. Hahsler, M., Piekenbrock, M. & Doran, D. dbscan: fast density-based clustering with R. J. Stat. Softw. https://doi.org/10.18637/jss.v091.i01 (2019).

51. Andrews, S. FastQC: a quality control tool for high throughput sequence data. Babraham Bioinformatics http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (2010).

52. Ewels, P., Magnusson, M., Lundin, S. & Käller, M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics 32, 3047–3048 (2016).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

53. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

54. Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).

    Article  CAS  PubMed  Google Scholar 

55. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).

    Article  PubMed  PubMed Central  Google Scholar 

56. Rozier, F., Mirabet, V., Vernoux, T. & Das, P. Analysis of 3D gene expression patterns in plants using whole-mount RNA in situ hybridization. Nat. Protoc. 9, 2464–2475 (2014).

    Article  CAS  PubMed  Google Scholar 

57. Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. Nat. Methods 9, 676–682 (2012).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

58. Preibisch, S., Saalfeld, S. & Tomancak, P. Globally optimal stitching of tiled 3D microscopic image acquisitions. Bioinformatics 25, 1463–1465 (2009).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

59. Boruc, J. et al. Functional modules in the Arabidopsis core cell cycle binary protein–protein interaction network. Plant Cell 22, 1264–1280 (2010).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

60. Zhang, H. zhanghao-njmu/SCP. GitHub https://github.com/zhanghao-njmu/SCP (2023).

61. Wolf, F. A. et al. PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells. Genome Biol. 20, 59 (2019).

    Article  PubMed  PubMed Central  Google Scholar 

62. Haghverdi, L., Büttner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly reconstructs lineage branching. Nat. Methods 13, 845–848 (2016).

    Article  CAS  PubMed  Google Scholar 

63. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 19, 15 (2018).

    Article  PubMed  PubMed Central  Google Scholar 

64. Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32, 381–386 (2014).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

65. Qiu, X. et al. Reversed graph embedding resolves complex single-cell trajectories. Nat. Methods 14, 979–982 (2017).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

66. Cao, J. et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496–502 (2019).

    Article  CAS  PubMed  PubMed Central  ADS  Google Scholar 

67. Rue-Albrecht, K., Lun, A., Soneson, C. & Stadler, M. velociraptor. R package version 1.10.0 https://bioconductor.org/packages/velociraptor (2026).

68. Bergen, V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat. Biotechnol. 38, 1408–1414 (2020).

    Article  CAS  PubMed  ADS  Google Scholar 

69. La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018).

    Article  PubMed  PubMed Central  ADS  Google Scholar 

70. Lun, A. DropletUtils. R package version 1.20.0 https://doi.org/10.18129/B9.BIOC.DROPLETUTILS (2018).

71. Morabito, S. et al. Single-nucleus chromatin accessibility and transcriptomic characterization of Alzheimer’s disease. Nat. Genet. 53, 1143–1155 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

72. Morabito, S., Reese, F., Rahimzadeh, N., Miyoshi, E. & Swarup, V. hdWGCNA identifies co-expression networks in high-dimensional transcriptomics data. Cell Rep. Methods 3, 100498 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

73. Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008).

    Article  PubMed  PubMed Central  Google Scholar 

74. Aboyoun, P., Pagès, H. & Lawrence, M. GenomicRanges. R package 1.52.0 https://doi.org/10.18129/B9.BIOC.GENOMICRANGES (2017).

75. Childs, J. E. et al. Relapse to cocaine seeking is regulated by medial habenula NR4A2/NURR1 in mice. Cell Rep. 43, 113956 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

76. Kolde, R. raivokolde / pheatmap. GitHub https://github.com/raivokolde/pheatmap (2025).

77. Miao, Z. & Kim, J. Uniform quantification of single-nucleus ATAC-seq data with Paired-Insertion Counting (PIC) and a model-based insertion rate estimator. Nat. Methods 21, 32–36 (2024).

    Article  CAS  PubMed  Google Scholar 

78. Zhang, K., Zemke, N. R., Armand, E. J. & Ren, B. A fast, scalable and versatile tool for analysis of single-cell omics data. Nat. Methods 21, 217–227 (2024).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

79. You, Y. et al. Benchmarking UMI-based single-cell RNA-seq preprocessing workflows. Genome Biol. 22, 339 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

80. Zhang, Y. et al. Model-based analysis of ChIP–seq (MACS). Genome Biol. 9, R137 (2008).

    Article  PubMed  PubMed Central  Google Scholar 

81. Stuart, T., Srivastava, A., Madad, S., Lareau, C. A. & Satija, R. Single-cell chromatin state analysis with Signac. Nat. Methods 18, 1333–1341 (2021).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

82. Mangiola, S. et al. sccomp: robust differential composition and variability analysis for single-cell data. Proc. Natl Acad. Sci. USA 120, e2203828120 (2023).

    Article  CAS  PubMed  PubMed Central  Google Scholar 

83. Legland, D., Arganda-Carreras, I. & Andrey, P. MorphoLibJ: integrated library and plugins for mathematical morphology with ImageJ. Bioinformatics 32, 3532–3534 (2016).

    Article  CAS  PubMed  Google Scholar 

84. Scrucca, L., Fraley, C., Murphy, T. B. & Adrian E. R. Model-Based Clustering, Classification, and Density Estimation Using mclust in R (Chapman and Hall/CRC, 2023).

85. Oliva, M. & Lister, R. Single-cell atlas of the Arabidopsis root tip - part 1 [Dataset]. Zenodo https://doi.org/10.5281/zenodo.17411120 (2025).

86. Oliva, M. & Lister, R. Single-cell atlas of the Arabidopsis root tip - part 2 [Dataset]. Zenodo https://doi.org/10.5281/zenodo.17411065 (2025).

87. Infrastructure and platforms. KAYA Research Compute Platform https://doi.org/10.26182/DH5S-V482 (The University of Western Australia, 2025).

Download references

Authors and Affiliations

1. ARC Centre of Excellence in Plant Energy Biology, School of Molecular Sciences, University of Western Australia, Perth, Western Australia, Australia

   Marina Oliva, Tim Stuart, Dave Tang, Daniel Poppe, Jahnvi Pflueger & Ryan Lister

2. ARC Centre of Excellence in Plants for Space, School of Molecular Sciences, University of Western Australia, Perth, Western Australia, Australia

   Marina Oliva, Adil Khan, Milly Petterson & Ryan Lister

3. Genome Institute of Singapore, Agency for Science Technology and Research (A*STAR), Singapore, Singapore

   Tim Stuart

4. The Harry Perkins Institute of Medical Research, QEII Medical Centre and Centre for Medical Research, The University of Western Australia, Perth, Western Australia, Australia

   Olivier Clement, Dawid Makosa, Daniel Poppe, Jahnvi Pflueger & Ryan Lister

5. The Walter and Eliza Hall Institute of Medical Research, Melbourne, Victoria, Australia

   Jafar S. Jabbari & Scott Gigante

6. Australian Genome Research Facility, Melbourne, Victoria, Australia

   Jafar S. Jabbari

7. La Trobe Institute for Sustainable Agriculture and Food, Department of Animal, Plant and Soil Sciences, School of Life Sciences, La Trobe University, Bundoora, Victoria, Australia

   Jonathan Michael Dragwidge & Mathew G. Lewsey

8. State Key Laboratory of Plant Environmental Resilience, College of Life Science, Zhejiang University, Hangzhou, P. R. China

   James Whelan

9. Australian Research Council Research Hub for Protected Cropping, La Trobe University, Bundoora, Victoria, Australia

   Mathew G. Lewsey

10. ARC Centre of Excellence in Plants for Space, La Trobe University, Bundoora, Victoria, Australia

    Mathew G. Lewsey

Contributions

M.O. and R.L. conceived the study and designed the experiments and analyses. M.O. performed the experiments and analysed the data. T.S. developed the cell genotyping and data integration methods, and analysed the bulk RNA-seq libraries. D.T. performed preliminary analyses of Arabidopsis accessions and assisted in establishing various analytical tools. A.K. and M.P. conducted some of the HCR experiments. O.C. contributed to testing spatial transcriptomics methods and generating single-cell libraries. D.M. analysed the snATAC–seq data. S.G. contributed to the development of the cell genotyping method. D.P. and J.P. sequenced the libraries. J.M.D., J.S.J., J.W. and M.G.L. assisted with pilot experiments. M.O. and R.L. wrote the manuscript. All authors reviewed and approved the final version of the manuscript.

Corresponding authors

Correspondence to Marina Oliva or Ryan Lister.

Competing interests

The authors declare no competing interests.

Peer review information

Nature thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions