Brassinosteroid gene regulatory networks at cellular resolution

Brassinosteroids (BRs) are plant steroid hormones that regulate diverse processes such as cell division and cell elongation. BRs control thousands of genes through gene regulatory networks that vary in space and time. By using time-series single-cell RNA-sequencing to identify BR-responsive gene expression specific to different cell types and developmental stages of the Arabidopsis root, we uncovered the elongating cortex as a site where BRs trigger a shift from proliferation to elongation associated with increased expression of cell wall-related genes. Our analysis revealed HAT7 and GTL1 as BR-responsive transcription factors that regulate cortex cell elongation. These results establish the cortex as an important site for BR-mediated growth and unveil a BR signaling network regulating the transition from proliferation to elongation, illuminating new aspects of spatiotemporal hormone response.


Introduction
During development, cells pass through different states as they acquire identities and progress towards end-stage differentiation (1). Gene regulatory networks (GRNs) control this progression and must be tuned according to developmental stage, cell identity, and environmental conditions (2,3). Signaling molecules such as hormones are central players in coordinating these networks, but it has been challenging to disentangle how cell identities, developmental states, and hormone responses influence one another. Recent technological advances in single-cell RNA-sequencing (scRNA-seq) (2,4) and tissue-specific gene manipulations (5) make it possible to address this challenge using the Arabidopsis root as a model system. Brassinosteroids (BRs) are a group of plant steroid hormones that affect cell division and cell elongation during root growth (6)(7)(8). BRs are sensed at the plasma membrane by BRI1 family receptors (9)(10)(11), initiating signal transduction events that activate BES1 and BZR1 family transcription factors (TFs) to control thousands of genes (12)(13)(14)(15). The BR GRN is typically represented singularly without consideration of cell specificity (15,16), but BRs lead to different responses depending on the developmental context (17)(18)(19). By profiling BR responses across the cell types and developmental stages of the root using scRNA-seq, we discovered that BRs strongly affect gene expression in the elongating cortex. Reconstruction of cortex trajectories over a scRNAseq time course showed that BRs trigger a shift from proliferation to elongation, which is associated with up-regulation of cell wall-related genes. Loss of BR signaling in the cortex using tissue-specific CRISPR reduced cell elongation. Our time-course data allowed us to infer BR-responsive GRNs, which led to the identification of HAT7 and GTL1 as validated regulators of BR response in the elongating cortex. These datasets represent more than 180,000 single-cell transcriptomes, providing a view of BR-mediated GRNs at unprecedented resolution.

Results
scRNA-seq reveals differential brassinosteroid response in the Arabidopsis root To investigate spatiotemporal BR responses in the root, we used a sensitized system, which involved inhibiting BR biosynthesis using Brassinazole (BRZ) (20), then reactivating signaling with Brassinolide (BL), the most active BR (8,16) (fig. S1). We treated 7-day-old primary roots for 2 hours with BL or a corresponding mock BRZ control and performed scRNA-seq on protoplasts isolated from three biological replicates of 0.5cm root tips (containing meristem, elongation and early differentiation zones) using the 10X Genomics Chromium system (Methods). To annotate cell types and developmental stages, we performed label transfer based on our single-cell atlas of the Arabidopsis root (21). We distinguished between two domains of the meristem: the proliferation domain, where cells have a high probability to divide, and the transition domain, where cells divide less frequently but have not yet begun rapid expansion ( fig. S2A-D and Data S2) (22). After data integration, the 11 major cell types and eight devel- opmental stages identified were logically arranged in 2D uniform manifold approximation and projection (UMAP) space as previously described for root datasets ( fig. S3A-B) (21,23). Marker genes characteristic of cell types and developmental stages remained enriched, suggesting that although BRs alter the expression of thousands of genes, cell identities can be successfully aligned through integration ( fig. S3C). scRNA-seq captures spatiotemporal patterns of BRresponsive gene expression Previous studies have profiled BR-responsive gene expression in bulk tissue or in a handful of cell types, conflating cell type and developmental stage (24)(25)(26). To obtain better spatiotemporal resolution, we performed differential expression analysis for each combination of cell type and developmental stage using pseudobulk expression profiles (see methods). We identified over 8,000 differentially expressed genes (DEGs; Fold-change >1.5, False discovery rate <0.05; Fig.  1A and fig. S4A-B), which were enriched in BES1 and BZR1 targets and had significant overlap with previously identified BR-regulated genes ( fig. S4C and Data S3). Strikingly, we found that 37% of DEGs were significantly altered in a single cell type/developmental stage and more than 82% were differentially expressed in 5 or fewer cell type/developmental stage combinations ( fig. S4B). This indicates that although BRs broadly influence gene expression, they modulate distinct sets of genes in different spatiotemporal contexts. Among the tissues with many DEGs was the epidermis, as previously described (7,17,24,25,27). Atrichoblasts, or non-hair cells in the epidermis were particularly affected, showing marked changes across both the meristem and elongation zone. Unexpectedly, our data also indicated that BRs strongly influence gene expression in the cortex, especially in the elongation zone ( Fig. 1A and fig. S4A). The cortex has been linked to plant environmental interactions, including response to water limitation (28) and hydrotropism (29), but it is not known how BRs modulate gene expression in this tissue nor what processes are affected. To address these questions, we focused on BR-mediated gene expression in the elongating cortex. Cell wall-related genes are up-regulated by BRs in the elongating cortex We found that BR treatment led to approximately 1,000 up-regulated genes and about the same number of downregulated genes in the elongating cortex ( Fig. 1B and  fig. S4A). Gene ontology (GO) analysis indicated the upregulated genes were strongly enriched for genes related to "cell wall organization or biogenesis", which is intriguing given the role of BRs in promoting cell elongation ( fig.  S4D). The cell wall-related DEGs included CELLULOSE SYNTHASES (CESAs), CELLULOSE SYNTHASE INTER-ACTIVE1 (CSI1), and cell-wall loosening enzymes such as EXPANSINS and XYLOGLUCAN ENDOTRANSGLUCOSY-LASES. Cell-wall-related genes such as CESAs have been demonstrated to be direct targets of BES1 and BZR1 (14,15,30), but their spatiotemporal regulation, especially in the cortex, has not been reported.
To monitor their responsiveness to BRs, we generated transcriptional reporters for three of the DEGs with distinct spatiotemporal patterns (Fig. 1C and fig. S5A-E). For example, CELL WALL / VACUOLAR INHIBITOR OF FRUCTOSI-DASE 2 (C/VIF2) was enriched in the transition domain and elongation zone of the cortex and induced by BL (Fig. 1C). These results confirm that our differential expression analysis captures spatiotemporal BR responses and raise the possibility that BR induction of cell-wall-related genes is associated with cortex cell elongation. Waddington optimal transport identifies BR induction of cell wall genes in the cortex associated with the switch to elongation A recently developed analytical approach for analyzing expression trends over a scRNA-seq time course is Waddington-OT (WOT), which connects snapshots of gene expression to facilitate trajectory reconstruction (31). WOT identifies putative ancestors for a given set of cells at earlier time points and descendants at later time points (31). To better understand how BRs influence cell wall-related gene expression we performed scRNA-seq at six time points beginning with BRZ treatment (time 0) and BL treatments for 30 minutes, 1 hour, 2 hours, 4 hours or 8 hours (Fig. 1, D and E). These time points capture the rapid root elongation triggered by re-addition of BRs (24). To examine the trajectories leading to activation of cell wallrelated genes in the elongating cortex, we applied WOT (31) and created a cell wall gene signature using 107 cell wallrelated genes that were induced by BL in the elongating cortex (Data S4). We monitored the relative expression of these genes, resulting in a "cell wall score" for each cell in the time course (see methods). Cortex cells had a higher cell wall score compared to other cell types, which increased with BL treatment (Fig. 1F-G and fig. S6A-B), confirming that the cell wall score represents a BR-responsive module in the cortex. At the 2 hour BL time point, more than 20% of cortex cells had a cell wall score greater than 1, whereas only 5% or fewer cells in other cell types exhibited scores this high (Fig. 1F). We therefore designated cells with a cell wall score of at least 1 as "cortex cell wall+" to indicate their exceptional BR response (Fig. 1, F and G). An advantage of WOT analysis is that it does not rely on pre-specified boundaries between developmental zones. We used this property to examine the relationship between developmental stage annotation and cell wall score. Under BRZ treatment, cortex cell wall+ cells were sparse and predominately annotated as transition domain. Upon BL treatment, the annotation of cortex cell wall+ cells shifted to the elongation zone (Fig. 1G). These results suggest that BRs are involved in initiating elongation of cortex cells via activation of cell wall genes. Differential expression along WOT trajectories identifies BR-responsive transcription factors Using the cells at the 2 hour time point as a reference, we looked at the probability of cells being ancestors or descendants of cortex cell wall+ cells (Fig. 1H). We also constructed a similar trajectory for the remaining cortex cells, which were designated "cortex cell wall-". To reveal potential regulators of cell wall-related genes in the cortex, we performed probabilistic differential expression analysis along WOT trajectories, contrasting cells assigned to cortex cell wall+ versus cortex cell wall-states at each time point (see methods). Among the DEGs identified were known TFs in the BR pathway including BES1 (12), BIM1 (32), and IBL1 (33); Fig.  1I and Data S4). We also identified additional TFs including JACKDAW (JKD), which is involved in ground tissue specification (34), the class I HD-ZIP TF HAT7 (35, 36) and GTL1. Because JKD was uniquely identified in this analysis, we confirmed its BR-responsiveness (Figures S6, C and D). These results indicate that WOT trajectories can identify BR-responsive TFs that may be involved in regulating cell wall-related genes in the cortex. Analysis of the triple receptor mutant bri1-T reveals changes in cortex expression Since our results indicated that exogenous BRs lead to activation of cell wall-related genes in the elongating cortex, we asked if this is also the case for endogenous BRs. A gradient of BRs is present along the longitudinal axis of the root, with low BR levels in the proliferation domain (37). BR biosynthesis increases as cells enter the transition domain and peaks in the elongation zone, shootward of which is a BR signaling maximum (24,37). Interpretation of this endogenous BR gradient requires receptor BRI1 and its close homologs BRL1 and BRL3 (9,38,39). To identify differentially expressed genes, we performed two replicates of scRNA-seq on the BR-blind bri1brl1brl3 triple mutant (bri1-T) along with paired wild-type controls ( Fig.  2A-C). A previous study profiled single cells from bri1-T, suggesting potential BR-responsiveness of the cortex (40), but these data were from a single replicate, were compared to a wild-type sample from a different study, and could not resolve developmental stage-specificity. In contrast, our analysis across both cell types and developmental stages identified the elongating cortex as exhibiting substantial differential gene expression (Fig. 2,D and E, fig. S7A and Data S3). The genes down-regulated in the elongating cortex of bri1-T were enriched for the GO term "cell wall organization or biogenesis" ( Figure S7B). These data indicate that endogenous BR signaling promotes the expression of cell wallrelated genes in the elongating cortex. The epidermis is widely described as the major site for BRpromoted gene expression in the root (7,8,17,24,25,41). Previous studies showed that epidermal expression of BRI1 was sufficient to rescue morphological phenotypes including meristem size of bri1-T (7,17,42). To determine the extent to which BR-regulated gene expression is restored, we performed scRNA-seq on pGL2-BRI1-GFP/bri1-T -a line in which BRI1 is expressed in atrichoblast cells of the epidermis of bri1-T (7,25,42). We identified over 8,000 DEGs in comparison with wild type (Fig. 2F) and in comparison with bri1-T (Fig. 2G-I and fig. S8A-F), indicating that gene expression remains dramatically perturbed and that this is far from a complete rescue of the bri1-T phenotype. Tissue-specific CRISPR of BRI1 confirms a role for the cortex in BR-mediated cell expansion Characterization of cell-type-specific BR signaling has relied on tissue-specific complementation lines, which led to conflicting results and overlooked the role of BR signaling in the cortex (7,17,24,25,40,42,43). To selectively block BR signaling in cell types of interest we performed tissue-specific CRISPR (5) of BRI1. We used a bri1 mutant complemented with pBRI1-BRI1-mCitrine (Fig. 3A) into which we introduced Cas9 driven by tissue-specific promoters to knock out BRI1 either in the epidermis and lateral root cap (pWER-BRI1-CRISPR) or in the cortex (pCO2-BRI1-CRISPR). mCitrine signals were absent in the expected locations of the tissue-specific CRISPR lines, confirming their efficacy and specificity ( Since our scRNA-seq data indicated that BRs promote the expression of cell wall-related genes in the elongating cortex, we hypothesized that loss of BR signaling in the cortex would affect final cell size. Indeed, pCO2-BRI1-CRISPR lines displayed significantly shorter mature cortex cells, while meristematic cortex cell length was relatively unaffected (Fig. 3C-E). In contrast, epidermal knockout of BRI1 in pWER-BRI1-CRISPR lines resulted in both reduced meristem cell size and reduced mature cortex cell length (Fig. 3C-E), which is consistent with the reported role of epidermal BR signaling (7,18,24,25,41) and BR-responsiveness across developmental zones of the epidermis in our scRNA-seq data. These results indicate that in addition to the epidermis, BR signaling in the cortex is required to promote cell expansion in the elongation zone. The cortex could instruct anisotropic growth through its physical connection with the epidermis, but as the outermost tissue, relaxation of the epidermis appears to be required to allow for cell elongation (44,45). This may explain the apparent widening of cortex cells in pWER-BRI1-CRISPR lines. BRI1 was also reported to rescue bri1-T morphology when expressed in the developing phloem using the CVP2 promoter (40,42,43), but gene expression was not fully restored to wild-type levels in either epidermal or phloem rescue lines. Our scRNA-seq of epidermal pGL2-BRI1-GFP/bri1-T lines showed patterns of gene expression distinct from either wild type or bri1-T. Similarly, scRNA-seq of pCVP2-BRI1-CITRINE/bri1-T indicated an intermediate state between wild type and bri1-T (40). Notably, BRI1 driven by its native promoter was still present in the stele of our tissuespecific CRISPR lines when we observed phenotypic defects, suggesting that, unlike pCVP2-BRI1, native expression of BRI1 in the stele is not sufficient for BR-induced cell elongation and root growth. These results confirm the role of the epidermis in BR-regulated root growth and reveal the function of cortex in BR-mediated cell expansion, demonstrating how scRNA-seq can identify a new spatiotemporal context for hormone signaling. HAT7 and GTL1 are BR-responsive regulators along cortex trajectories To define a core set of genes associated with BR response along cortex trajectories we first compared genes induced in the cortex by BL treatment with those down-regulated in the cortex of bri1-T. Of the 768 genes in common, we then asked which vary along developmental time in wild-type cortex trajectories (21). The intersection of these three lists identified a core set of 163 BR responsive DEGs ( Fig. 4A and Data S3). Consistent with regulation by BRs, 69% of the core DEGs are BES1 and BZR1 direct targets from ChIP experiments (14,15,46). Expression along cortex pseudotime illustrates induction by BL treatment and down-regulation in bri1-T (Fig.  4B). HAT7 and GTL1 were induced along these trajectories, suggesting a potential role for these TFs in controlling BRregulated gene expression in the cortex (Fig. 4C-E and fig.  S10A).
To gain insight into their roles, we generated translational reporter lines for HAT7 and GTL1. Under control conditions, pHAT7-HAT7-mCitrine lines showed expression in the transition domain and elongation zone of the cortex (Fig. 4F). We also observed HAT7 signals in the epidermis and endodermis, in line with expression patterns in our wild-type scRNA-seq atlas (21). HAT7 expression was decreased when BR biosynthesis was inhibited with BRZ, and restored upon BL treatment ( Fig. 4F and fig. S10A). pGTL1-GTL1-mCitrine was more broadly expressed, with increasing levels in the cortex and epidermis as cells progress from the transition domain to the elongation zone ( Fig. 4G and S10A). GTL1-mCitrine expression was reduced by BRZ and increased by BL treatment ( Figure 4G). These results confirm that BRs promote the expression of HAT7 and GTL1 coinciding with the onset of cell elongation. Furthermore, HAT7 and GTL1 are direct targets of BES1 and BZR1 (14,15,46), suggesting that they may be part of the BR-directed GRN activated as cells progress from proliferation to elongation. Previous studies inferred global (14,15,26,47) or temporally resolved GRNs (16) for BR response, but they lacked cell type and developmental-stage specificity. To infer GRN configurations across our BR time series we used CellOracle (see methods; Data S8) focusing on BL DEGs and associated TFs. Analysis of network importance scores such as centrality measures is a powerful approach to prioritize candidate regulators among DEGs (48). Since the cell wall signature peaked at 2 hours after BL treatment, we prioritized TFs with high network centrality scores in the elongating cortex at this time point. HAT7 was the top-ranked TF in terms of degree centrality and three close homologs: HB13, HB20 and HB23 were also among the top 10 TFs (Fig. 5A-B and Data S8). Together HAT7, HB13, HB20 and HB23 make up the alpha clade HD-ZIP I TFs (35,36). We used CRISPR to generate hat7 loss-of-function mutants but did not observe strong phenotypes in terms of cortex cell elongation ( Fig. S11A-C). Since HB13, HB20 and HB23 are induced by BRs and predicted to regulate cell wall-related genes in our GRNs (Fig. 5B, fig. S10A-B and Data S10), we next generated hat7 hb13 hb20 hb23 quadruple mutants via multiplex CRISPR. Mature cortex cell length was reduced by approximately 25% in two independent quadruple mutants ( Fig. 5C-E and fig. S11A-C), providing strong evidence that HAT7 and its homologs are required for cell elongation. Despite the decrease in final cell length, the root length of the quadruple mutant was not dramatically reduced (Fig. S11A), suggesting that the decrease in cell length is at least partially compensated for by increased cell production. We next investigated GTL1, which was the 5th highestranked TF in the BL 2 hour elongating cortex GRN (Fig. 5A-B). Given that GTL1 was shown to function redundantly with its close homolog DF1 in terminating root hair growth (49), we examined gtl1 df1 double mutants, finding significantly shorter mature cortex cell lengths and shorter roots compared to wild-type (Fig. 5C-E and fig. S11A-C). DF1 was challenging to detect in scRNA-seq (fig. S10A) due to its low expression level (49), but we observed increasing trends of DF1 expression along WOT trajectories in the BL time course, especially in cortex cell wall+ cells ( fig. S10B) which was verified using a pDF1-DF1-GFP reporter (fig. S10C). Together, our genetic analysis of HAT7 and GTL1 family TFs illustrates the power of GRN-mediated discovery of regulatory factors in spatiotemporal BR response. BES1 and GTL1 physically interact and share a common set of target genes Since BES1 is known to interface with other TFs in controlling BR-regulated gene expression, we compared target genes for BES1 and BZR1 (14,15,46) to ChIP targets of GTL1 and DF1 (49). BES1 and BZR1 share 3,020 common targets with GTL1, significantly more than expected by chance (fig. S12A, P < 0.001, Fisher's exact test). Similarly, BES1 and BZR1 share 2,490 common targets with DF1 (Figure S12A, P < 0.001, Fisher's exact test). When compared to BR-regulated genes from scRNA-seq, BES1 and GTL1 targets showed the strongest enrichment in genes up-regulated by BRs in the transition domain and elongation zone of the cortex (fig. S12B), with 297 common targets of both BES1 or BZR1 and GTL1 being induced in the elongating cortex by BL treatment. Given the strong overlap between BES1 and GTL1 targets, we hypothesized that these TFs physically interact to regulate a common set of genes. Co-immunoprecipitation showed that GTL1-FLAG pulled down BES1-GFP (fig. S12C). These results suggest that BRs induce GTL1 and subsequently BES1 and GTL1 interact to control a common set of target genes. This type of feed-forward loop could provide a mechanism to amplify the BR signal and/or to direct BES1 to drive tissue-specific gene expression by interacting with other more specifically expressed TFs. scRNA-seq reveals cell-type-specific expression underlying gtl1 df1 phenotypes Our results indicate that gtl1 df1 mutants have reduced cortex cell elongation. On the other hand, gtl1 df1 mutants have longer trichoblasts (49). A downstream regulatory network that enables GTL1-mediated growth inhibition has been dissected in trichoblasts (49). To identify the cell-type-specific changes in gene expression underlying gtl1 df1 cortex phenotypes we performed scRNA-seq on gtl1 and df1 single mutants, and on the gtl1 df1 double mutant. Using pseudobulk  differential expression analysis, we detected relatively subtle changes in gtl1 or df1 single mutants compared to wild type (Fig. 6, A and B). In contrast, over 8,000 genes were differentially expressed in gtl1 df1 double mutants versus wildtype (Fig. 6C).
Over 1,000 genes were up-regulated across all developmental stages of the cortex of gtl1 df1, and an approximately equal number of genes were down-regulated. The majority of cortex DEGs were affected in the elongation zone (Fig. 6, D and E, fig. S12D and Data S3). Of the down-regulated genes in the cortex of the double mutant, 226 genes were also up-regulated by BL treatment. Furthermore, 31.3% of the core BR DEGs were down-regulated in the cortex of gtl1 df1, whereas only 6.8% were up-regulated (Data S3). These results suggest that GTL1 and DF1 promote the expression of a subset of BR-induced genes in the cortex.
Plotting gtl1 df1 DEGs along cortex pseudotime illustrated the down-regulation of several genes involved in cell elongation including CESA5 and AHA2 (Fig. 6E). These genes were significantly enriched for the GO term "cell wall organization or biogenesis" (fig. S12E). We next examined C/VIF2, because it is induced by BL in the cortex, but its expression decreased in cortex cells of gtl1 df1 (Fig. 6, F and G). A pC/VIF2-H2B-Venus reporter showed expression of C/VIF2 in the transition and elongation zone of the wild-type cortex, whereas its expression was reduced in the cortex of gtl1 df1 mutants ( Fig. 6H and Movie S1). The reduced expression of cell wall-related genes in gtl1 df1 mutants validates our celltype-specific BR GRNs and identifies a function of GTL1 in promoting cortex cell elongation in response to BRs.

Discussion
Understanding how hormone-mediated GRNs are controlled in space and time has the potential to enable the engineering of specific downstream responses to optimize plant growth under a changing environment (8,50). Plant hormones including BRs, auxin, gibberellins, and abscisic acid have been shown to exhibit tissue-specific responses (51)(52)(53)(54), but how the associated GRNs are modulated in different cell types at particular developmental stages is largely enigmatic. In this study, we profiled BR-responses across cell types, developmental stages and time points of treatment using scRNA-seq, providing a high-resolution map of signaling outputs. These data are publicly available as an interactive browser (https://shiny.mdc-berlin.de/ARVEX/). We identified the elongating cortex as a spatiotemporal context for BR signaling, where BRs activate cell wall-related genes and promote elongation. We further showed that HAT7 and GTL1 are BR-induced regulators along cortex trajectories that control cell elongation. These findings highlight the ability of single-cell genomics to identify context-specific TFs, a capability that could be leveraged to precisely engineer plant growth, development, and/or responses to stress. Our results reveal spatiotemporal BR responses and the underlying GRNs at unprecedented resolution.

Transgenic reporters
To generate new reporters for BR-responsive genes, we first added the FASTRED seed coat selection cassette (57, 58) and a MoClo (59) Level 1 acceptor site to the binary vector pICH86966 (Addgene plasmid #48075). pHAT7-HAT7-mCitrine and pGTL1-GTL1-mCitrine were assembled into this FASTRED destination vector using Level 1 BsaI golden gate assembly. To facilitate one-step promoterreporter construction, we assembled an AarI flanked RFP dropout using the overhangs described in the Mobius (60) upstream of Venus-H2B followed by the Ubiquitin10 terminator (tUBQ10), a plasma membrane marker (pUBQ10-mScarlet-LTI6-tNos), and a constitutive histone maker (pUBQ10-H2B-CFP-t19s). Promoters containing up to~3kb of sequence upstream of the ATG start codon for the gene of interest were PCR amilfied with AarI containing primers and used to replace the AarI-RFP module in golden gate reactions to generate Promoter-Venus-H2B constructs. Our Venus-H2B reporter included a Ubiquitin tag to decrease reporter perdurance. Although the plasma membrane marker and histone marker were included as positive controls in the constructs, they were not further analyzed in this study. Assemblies were confirmed by restriction digestion and sequencing, transformed into Agrobacterium, and used to transform Arabidopsis via floral dip (61). FASTRED positive T1 seeds were selected under a fluorescent dissecting scope and only lines with 3:1 segregation of seed coat fluorescence in the T2 generation were used. T2 lines with bright seed fluorescence were typically homozygous in our conditions. Therefore, we used bright T2 seeds or homozygous T3 seeds for experiments. We ensured that reporter signals were consistent across at least three independent transgenic lines.

Generation of mutant lines using multiplex CRISPR
We produced hat7 single mutants and hat7 hb13 hb20 hb23 quadruple mutants using FASTRED multiplex CRISPR constructs containing an intronized version of Cas9 (58,62). Two gRNAs were designed per gene using CHOP-CHOP (63). gRNA containing oligos were hybridized and cloned into pDGE sgRNA shuttle vectors using BpiI (Data S11). Each of the gRNA containing shuttle vectors were then assembled into pDGE666 (Addgene plasmid # 153231) using BsaI golden gate assembly, sequence verified, and transformed into wild-type Arabidopsis as described above. We selected FASTRED positive T1 seeds and subsequently screened FASTRED negative (putatively Cas9-free) T2 seeds for frameshift mutations using Sanger sequencing coupled with ICE analysis of CRISPR edits (64). The edits were similarly confirmed in the T3 generation and at least two homozygous alleles from independent lines were used for experiments.

Confocal microscopy
Confocal imaging for the majority of experiments was performed using a Zeiss 880 equipped with a 40X objective. Excitation and detection were set as follows: Venus and mCitrine, excitation at 488 nm and detection at 499-571 nm; GFP, excitation at 488 nm and detection at 493-558 nm; PI staining, excitation at 561 nm and detection at 605-695 nm. Confocal images were processed using the Fiji package of ImageJ (70). Tile scans were stitched and representative me-dian longitudinal sections for each image are shown. Identical settings were used for images that were directly compared. For BRI1 TSKO confocal, roots were imaged between a block of agar and cover glass in imaging chambers. Image acquisition was performed with a FluoView1000 inverted confocal microscope (Olympus) equipped with a dry 20X objective (NA 0.75) using 514 nm laser excitation and a spectral detection bandwidth of 500-530 nm for mCitrine and 535 nm laser excitation together with a spectral detection bandwidth of 570-670 nm for PI. For the confocal time-lapse video, 7-day-old wild type and gtl1 df1 seedlings expressing pC/VIF2-H2B-VENUS were placed on ½ MS agar blocks containing PI, placed side-by-side in chambered coverglass (Nunc Lab-Tek, ThermoFisher) and imaged under a vertical ZEISS LSM900 microscope equipped with a Plan-Apochromat M27 20x/ 0.8 n.a. objective. The root tip was imaged every 12 min and automatically tracked with the TipTracker software (71).
scRNA-seq profiling of Arabidopsis root protoplasts using the 10X Genomics chromium system scRNA-seq experiments were performed as previously described (21)  , centrifuged again at 500g for 3 minutes, and the pellet resuspended in washing solution at a concentration of~2000 cells/uL. We loaded 16,000 cells, with the aim to capture 10,000 cells per sample with the 10X Genomics Chromium 3' Gene expression v3 or v3.1 kits. Cell barcoding and library construction were performed following the manufacturer's instructions. cDNA and final library quality were verified using a Bioanalyzer High Sensitivity DNA Chip (Agilent) and sequenced on an Illumina NextSeq 500 or NovaSeq 6000 instrument to produce 100bp paired-end reads.
For BL scRNA-seq, we first grew plants on 1 µM BRZ to deplete endogenous BRs, then transferred plants to either a fresh BRZ plate or 100nM BL. We monitored the efficacy of these treatments using a constitutively expressed 35S-BES1-GFP line. In agreement with previous reports (12,74,75), BES1-GFP was predominantly present in the cytoplasm under low BR conditions resulting from BRZ treatment but accumulated in the nucleus following BL treatment (Fig. S1). We performed two separate BL scRNA-seq treatment experiments. The first consisted of a BRZ and 2 hour BL treatment. The second experiment included two additional replicates of BRZ and BL 2 hours along with the other time points in our time course (BL 0.5, 1, 4, and 8 hour treatments). Each of the BL treatments was staggered so that all samples were collected simultaneously. A total of 70,223 cells were recovered from the BL treatment scRNA-seq experiments. Wild type Col-0, bri1-T, and pGL2-BRI1-GFP/bri1-T were similarly profiled in a side-by-side scRNA-seq experiment under control conditions with two replicates per genotype, resulting in 34,861 cells. Lastly, scRNA-seq was performed on Wild type Col-0, gtl1, df1, and gtl1 df1 in duplicate under control conditions, resulting in 74,810 scRNA-seq expression profiles.
COPILOT uses a non-arbitrary scheme to remove empty droplets and dying or low-quality cells. We used one iteration of COPILOT filtering, which adequately separated high-quality cells from the background in our samples based on an examination of barcode rank plots. To address issues with doublets and outliers, the resulting high-quality cells were further filtered to remove the top 1% of cells in terms of UMI counts, and putative doublets were removed with DoubletFinder(80) using the estimated doublet rate from the 10X Genomics Chromium Single Cell 3' Reagent Kit user guide.
Normalization, annotation, and integration of scRNA-seq datasets Downstream analysis were carried out using Seurat version 3.1.5. Samples were first individually processed and examined. Data were normalized using SCTransform (81) and all detected genes were subsequently retained for analysis, except those from mitochondria, chloroplasts or those affected by protoplasting (absolute log2 fold-change >= 2) (21, 23). Principal component analysis (PCA) was performed by calculating 50 principal components using the RunPCA function (with approx=FALSE). UMAP non-linear dimensionality reduction was next calculated via the RunUMAP function using all 50 principal components with parameters n_neighbors = 30, min_dist = 0.3, umap.method = ' 'umap-learn' ', metric = ' 'correlation' '. These processing steps have been previously described (21) and are documented in jupyter notebooks as part of the COPILOT workflow.
To follow the developmental progression from the meristem to the elongation zone more closely, we updated the root atlas (21) developmental annotation to subdivide the meristem into the proliferation domain and transition domain as previously defined (22). The previous meristem annotation of the root atlas was based on correlation annotation by comparing each cell from scRNA-seq to bulk data from morphologically defined sections (82). On the other hand, HIGH PLOIDY2 was used to mark the meristem in a second bulk expression profile (83), which corresponds to the proliferation domain defined by Ivanov and Dubrovsky. Therefore  (85,86). We used the receiver operating characteristic (ROC) test implemented in Seurat FindMarkers to identify genes enriched in each developmental zone. A largely distinct set of markers was enriched in the transition domain ( fig. S2C and Data S2). These included genes involved in vesicle-mediated transport ( fig. S21D), in line with the observation that vesicle recycling activity is highest in this region (87). We transferred the cell type and developmental stage labels from the wild-type atlas (21) to each sample using the Seurat label transfer workflow (88,89). To align corresponding cell types and developmental stages, we integrated samples from each experiment using the Seurat reference-based integration pipeline (88,89). A sample from the atlas with the highest genes detected (sc_12) was used as a reference (21) and two previously described samples (dc_1 and dc_2) (23) were also included in each integration. PCA and UMAP were subsequently calculated for each integration object using the batch-corrected "integrated" assay as described above. Although sc_12, dc_1 and dc_2 were not used in any subsequent analysis, their inclusion at the integration step helped to generate a comparable UMAP projection among different integration objects that facilitates interpretability.

Plotting gene expression values on the UMAP projection
We subsequently examined changes in cell state caused by the BL treatments or in the mutants profiled by plotting the log-normalized, 'corrected' counts produced by the SCTransform function (Hafemeister and Satija, 2019) rather than the batch-corrected "integrated" values when visualizing changes in expression.
Pseudotime estimation and heatmaps of gene expression trends Cortex cells were extracted from the integrated Seurat objects (BR time course, bri1-T vs wild type and gtl1 df1 vs wild type). Pseudotime was then inferred on the SCT assay of the extracted cortex cells using CytoTRACE v0.

Pseudobulk differential expression analysis
Recent benchmarks point towards pseudobulk methods, which aggregate cell level counts for subpopulations of interest on a per-sample basis, as top performers for crosscondition comparisons in scRNA-seq (92,93). Therefore, we employed a pseudobulk approach implemented in muscat (Multi-sample multi-group scRNA-seq analysis tools) (92) to examine changes in each combination of cell type and developmental stage. Pseudobulk expression profiles were aggregated for each of these subpopulations by summing the raw counts using the aggregateData function. We then performed differential expression testing using the edgeR method (94) incorporated in the pbDS function. A term for the experimental batch and/or replicate was included in the contrast to adjust for potential batch effects. A gene was considered differentially expressed in a given subpopulation if the false discovery-rate adjusted p-value was <=0.05, absolute fold change was >=1.5 and detection frequency was >=10% in one of the conditions. Gene ontology enrichment analysis was conducted on the differentially expressed genes using the R package "gprofiler2" (95). Comparisons between DEG lists were performed using the GeneOverlap package (version 1.12.0; http://shenlab-sinai.github.io/shenlabsinai/). p-values for intersections between gene lists were computed using Fisher's exact test. Visualizations were generated using Seurat (88), ComplexHeatmap (96), and ggplot2 (97).
WOT differential expression along cortex cell wall + trajectories WOT constructs trajectories of cells from a reference time point by minimizing the difference over all genes (31). The algorithm requires as input the expression profiles of cells as well as an estimation of their proliferation rate. We estimated proliferation rates using imaging data (98), as previously described (99). As only the bottom 0.5 cm of each root was observed at each time, we expect some cells to exit the observed section due to proliferation. We estimated the number of cells that should exit the observed section based on the growth rate with the assumption that the section of root stays in equilibrium and assigned the calculated number of cells with the highest pseudotime a growth rate of zero so that they have no descendents in the observed root section at the next time point. We constructed trajectories using full gene expression profiles and evaluated the quality of the trajectories by checking the proportion of cells whose highest fate probability matched the annotation. We found that for 90% of cells the largest fate assigned by WOT matched the annotation, rising to 97% in the maturation zone where we have the greatest confidence in the annotation. The cell wall signature was calculated for each cell by taking the sum of Z-scores for each of the 107 BR-induced cell wall-related genes in the signature (GO:0071554), truncated to [-5,5]. We defined the cortex cell wall+ subset as cortex cells with a cell wall score greater or equal to 1. This threshold was chosen as it selected less than 5% of cells from other lineages while still retaining >20% of cortex cells at the 2 hour time point. Any cortex cell that did not belong to the cortex cell wall+ group was labelled as "cortex cell wall-". We performed differential expression on the cortex cell wall+ and cell wall-subsets at 2 hours, using WOT lineages to also perform differential expression on their putative ancestors and descendants. Statistical significance was evaluated using Welch's t-Test with adjusted p-values for multiple tests, requiring t FDR < 0.01. Results were ranked by the adjusted expression ratio p max p min + ε with ε = 0.1, where larger ε puts more emphasis on genes with non-zero expression in both groups. To generalize our WOT analysis we also constructed trajectories for each combination of cell type and developmental stage and performed differential expression analysis between each time point using the same process (Data S5).

Gene regulatory networks
In order to construct GRNs, we used CellOracle (v0.7.0) for single-cell GRN inference (48). In the first stage of the Cel-lOracle pipeline, a base GRN is defined, representing a global set of biologically plausible Transcription factor-Target interactions. We used publicly available scATAC-seq data from Arabidopsis roots (100) GSE155304:GSM4698760; (100) to determine regions of open chromatin. Cell Ranger ATAC (v1.2.0) was used to process raw scATAC seq data to call a peak-by-cell matrix. Cicero (v1.11.1) (101) was implemented to infer a co-accessibility map of chromatin regions. Transcription start sites were then annotated based on the Arabidopsis TAIR10 genome assembly. Finally, peaks with weak co-accessibility scores were filtered following instructions of CellOracle manual (https://morrislab.github.io/CellOracle.documentation/tutorials/base_grn.html).
To expand the number of TFs present in the base GRN we also included TF-Target interactions from DNA affinity purification sequencing (DAP-seq) (102) and a previously constructed integrative gene regulatory network (iGRN) (103). Our resulting base GRN contained 11.7 million interactions between 1,601 transcription factors and 31,019 target genes. In the second step of the CellOracle pipeline, a regularized machine learning approach is used to define active edges and their regulatory strength in clusters or subpopulations of scRNA-seq data. In this process, the expression of target genes is predicted based on regulatory transcription factor levels from the base GRN. Inactive edges with low predictive ability are pruned from the base GRN, revealing contextspecific GRN configurations (48). To test CellOracle on Arabidopsis root data, we first inferred GRN configurations for each of the 36 cell type and developmental stage combinations in our WT atlas (21) using the SCT normalized counts. We limited the base GRN to genes dynamically expressed along pseudotime for each cell type plus associated transcription factors (104). Each cell type GRN was then constructed with default parameters following the CellOracle manual. To filter network edges with the "filter_links" function, we retained the top 20,000 edges (pvalue <=0.01) for each subnetwork. This recovered known developmental regulators (Data S6 and Data S7) including MYB36 in the endodermis (105) and BRN1/BRN2 in the root cap (106), confirming that CellOracle analysis of Arabidopsis root scRNA-seq data can infer GRNs configurations for particular cell identities and states. We implemented similar procedures to infer context-specific GRN configurations for each cell type, developmental stage and time point of the BR time course samples (sc_43-50). We used transcription factors plus DEGs from BL 2 hour vs. BRZ pseudobulk analysis of each cell type/developmental zone combination. The resulting set of 201 GRN configurations spanned 767,970 edges between 1,164 transcription factors and 7,135 targets (Data S8). Network centrality measures were calculated using the built-in functions of the CellOracle pipeline (Data S9). The data needed to reproduce our results and jupyter notebooks demonstrating the processes are available on ARVEX (https://shiny.mdcberlin.de/ARVEX/).

Data and code availability
Single-cell RNA-seq data have been deposited at GEO: GSE212230.