Brassinosteroid gene regulatory networks at cellular resolution in the Arabidopsis root

1During development, cells pass through different states as they acquire identities and progress toward end-stage differentiation. Gene regulatory networks (GRNs) control this progression and must be tuned according to developmental stage, cell identity, and environmental conditions. 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. Technological advances in single-cell RNA-sequencing (scRNA-seq) and tissue-specific gene manipulations make it possible to address this challenge using the Arabidopsis root as a model system.

2Brassinosteroids are a group of plant steroid hormones that affect cell division and cell elongation during root growth. Brassinosteroids are sensed at the plasma membrane by BRASSINOSTEROID INSENSITIVE 1 (BRI1) family receptors, initiating signal transduction events that activate BES1 and BZR1 family transcription factors to control thousands of genes. The brassinosteroid GRN is typically represented singularly without consideration of cell specificity, even though brassinosteroids lead to different responses depending on the developmental context.

3By profiling brassinosteroid responses across the cell types and developmental stages of the root using scRNA-seq, we discovered that brassinosteroids affect gene expression in the elongating cortex. Reconstruction of cortex trajectories over a scRNA-seq time course showed that brassinosteroids trigger a shift from proliferation to elongation, which is associated with upregulation of cell wall-related genes. Loss of brassinosteroid signaling in the cortex using tissue-specific CRISPR reduced cell elongation. Our time-course data allowed us to propose brassinosteroid-responsive GRNs, which led to the identification of HAT7 and GTL1 as validated regulators of brassinosteroid response in the elongating cortex. These datasets represent 210,856 single-cell transcriptomes, providing a high-resolution view of brassinosteroid-mediated GRNs.

4To investigate spatiotemporal brassinosteroid responses in the Arabidopsis root, we used a sensitized system, which involved inhibiting brassinosteroid biosynthesis using Brassinazole (BRZ), then reactivating signaling with Brassinolide (BL), the most active brassinosteroid (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.5 cm root tips (containing meristem, elongation, and early differentiation zones) using the 10X Genomics Chromium system (Methods).

5To annotate cell types and developmental stages, we performed label transfer based on our single-cell atlas of the Arabidopsis root. We distinguished between two domains of the meristem: the proliferation domain, where cells have a high probability of dividing, and the transition domain, where cells divide less frequently but have not yet begun rapid expansion (fig. S2A–D and Data S2).

6After data integration, the 11 major cell types and eight developmental stages identified were logically arranged in 2D uniform manifold approximation and projection (UMAP) space as previously described for root datasets (fig. S3A–B). Marker genes characteristic of cell types and developmental stages remained enriched, suggesting that although brassinosteroids alter the expression of thousands of genes, cell identities can be successfully aligned through integration (fig. S3C).

7Previous studies have profiled brassinosteroid-responsive gene expression in bulk tissue or in a handful of cell types, conflating cell type and developmental stage. To obtain a 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 8,286 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 brassinosteroid-regulated genes (fig. S4C and Data S3).

8We found that 37% of DEGs were 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 brassinosteroids broadly influence gene expression, they modulate distinct sets of genes in different spatiotemporal contexts.

9Among the tissues with many DEGs was the epidermis, as previously described. Atrichoblasts, or non-hair cells in the epidermis were particularly affected, showing changes across both the meristem and elongation zone. Our data also indicated that brassinosteroids 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 and hydrotropism. However, it is unknown how brassinosteroids modulate gene expression in this tissue or what processes are affected. To address these questions, we focused on brassinosteroid-mediated gene expression in the elongating cortex.

10We found that brassinosteroid treatment led to 967 up-regulated genes and 1,156 down-regulated genes in the elongating cortex (Fig. 1B and fig. S4A). Gene ontology (GO) analysis indicated the up-regulated genes were enriched for genes related to “cell wall organization or biogenesis”, which is intriguing given the role of brassinosteroids in promoting cell elongation (fig. S4D). The cell wall-related DEGs included CELLULOSE SYNTHASES (CESAs), CELLULOSE SYNTHASE INTERACTIVE1 (CSI1), and cell-wall loosening enzymes such as EXPANSINS and XYLOGLUCAN ENDOTRANSGLUCOSYLASES. Cell-wall-related genes such as CESAs have been demonstrated to be direct targets of BES1 and BZR1, but their spatiotemporal regulation, especially in the cortex, has not been reported.

11To monitor their responsiveness to brassinosteroids, 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 FRUCTOSIDASE 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 brassinosteroid responses and raise the possibility that brassinosteroid induction of cell-wall-related genes is associated with cortex cell elongation.

12To better understand how brassinosteroids 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, and 8 hours (Fig. 1, D and E). These time points capture the rapid root elongation triggered by the re-addition of brassinosteroids. Waddington-Optimal Transport (WOT) is an analytical approach for analyzing expression trends over a scRNA-seq time course. WOT connects snapshots of gene expression to facilitate trajectory reconstruction, identifying putative ancestors for a given set of cells at earlier time points and descendants at later time points.

13To examine the trajectories leading to the activation of cell wall-related genes in the elongating cortex, we applied WOT and created a cell wall gene signature using 107 cell wall-related genes that were induced by BL in the elongating cortex. 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. 2A–B and fig. S6A–B), confirming that the cell wall score represents a brassinosteroid-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. In contrast, only 5% or fewer cells in other cell types exhibited scores this high (Fig. 2B). We therefore designated cells with a cell wall score of at least 1 as “responsive cortex cells” to indicate their exceptional brassinosteroid response (Fig. 2A).

14An 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, responsive cortex cells were sparse and predominately annotated as transition domain. Upon BL treatment, the annotation of responsive cortex cells shifted to the elongation zone (Fig. 2B). Using the cells at the 2 hour time point as a reference, we looked at the probability of cells being ancestors or descendants of responsive cortex cells (Fig. 2C). We also constructed a similar trajectory for the remaining cortex cells, which were designated “non-responsive cortex cells”. We visualized these probabilities on barycentric coordinates as triangle plots where vertices represent a 100% chance of becoming the labeled state, and interior positions indicate intermediate probabilities. Cells can be colored according to gene expression, cell types, or other quantities. The cortex non-responsive state was represented by the lower left vertex, the cortex responsive state on the lower right vertex, and the top vertex represented other states. This illustrated that cortex transition domain cells were predisposed toward the responsive state under BRZ conditions and subsequently shifted from transition domain to elongation zone annotation upon BL treatment (Fig. 2, D and E). These results suggest that brassinosteroids are involved in initiating the elongation of cortex cells via the activation of cell wall-related genes.

15To 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 responsive versus non-responsive states at each time point (see methods). Among the DEGs identified were known transcription factors in the brassinosteroid pathway, including BES1, BES1-INTERACTING MYC-LIKE1 (BIM1), and IBH1-LIKE 1 (IBL1); Fig. 2F). We also identified additional transcription factors including JACKDAW (JKD), which is involved in ground tissue specification, the class I HD-ZIP transcription factor HAT7 and GTL1. Because JKD was identified by WOT and a previous BL RNA-seq dataset but not pseudobulk differential expression analysis, we confirmed its brassinosteroid-responsiveness (fig. S6, C and D). These results indicate that WOT trajectories can identify brassinosteroid-responsive transcription factors that may be involved in regulating cell wall-related genes in the cortex.

16Since our results indicated that exogenous brassinosteroids lead to the activation of cell wall-related genes in the elongating cortex, we asked if this is also the case for endogenous brassinosteroids. A gradient of brassinosteroids is present along the longitudinal axis of the root, with low brassinosteroid levels in the proliferation domain. Brassinosteroid biosynthesis increases as cells enter the transition domain and peaks in the elongation zone, shootward of which is a brassinosteroid signaling maximum. Interpretation of this endogenous brassinosteroid gradient requires receptor BRI1 and its homologs BRL1 and BRL3.

17To identify differentially expressed genes, we performed two replicates of scRNA-seq on the brassinosteroid-blind bri1brl1brl3 triple mutant (bri1-T) along with paired wild-type controls (Fig. 3A–C). A previous study profiled single cells from bri1-T, suggesting potential brassinosteroid-responsiveness of the cortex. However, these data were from a single replicate, were compared to a wild-type sample from a different study, and did not resolve developmental stage-specificity. In contrast, our analysis across both cell types and developmental stages identified the elongating cortex as exhibiting differential gene expression (Fig. 3, 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” (fig. S7B). BL treatment increased the proportion of elongating cortex cells (fig. S7C), but fewer elongating cortex cells were observed in bri1-T compared to WT (fig. S7D).

18To further verify these changes in gene expression, we performed an additional side-by-side scRNA-seq experiment in which we inhibited endogenous brassinosteroids via BRZ treatment. We profiled 16,642 WT control and 14,320 WT BRZ-treated cells and detected 6,928 DEGs in BRZ vs control using pseudobulk differential expression analysis of each cell type/developmental stage combination (fig. S8A–C). Consistent with the idea that brassinosteroids promote cell wall-related gene expression in the elongating cortex, BRZ down-regulated genes in the elongating cortex were enriched for genes associated with “cell wall organization or biogenesis” (fig. S8D). Together our analysis of bri1-T and BRZ treatment indicates that endogenous brassinosteroid signaling promotes the expression of cell wall-related genes in the cortex associated with the onset of elongation.

19The epidermis is widely described as the major site for brassinosteroid-promoted gene expression in the root. Previous studies showed that epidermal expression of BRI1 was sufficient to rescue morphological phenotypes, including meristem size of bri1-T. To determine the extent to which brassinosteroid-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. We identified 8,188 DEGs in comparison with wild type (Fig. 3F) and 8,046 DEGs in comparison with bri1-T (Fig. 3G and fig. S9A–D), indicating that gene expression remains perturbed (fig. S10A–D) and that this is far from a complete rescue of the bri1-T phenotype.

20Characterization of cell-type-specific brassinosteroid signaling has relied on tissue-specific complementation lines, which led to conflicting results and overlooked the role of brassinosteroid signaling in the cortex. To selectively block brassinosteroid signaling in cell types of interest, we performed tissue-specific CRISPR of BRI1. We used a bri1 mutant complemented with pBRI1-BRI1-mCitrine (Fig. 4A) 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 (Fig. 4A–C and fig. S11A–C).

21Since our scRNA-seq data indicated that brassinosteroids promote the expression of cell wall-related genes in the elongating cortex, we hypothesized that loss of brassinosteroid 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. 4D–E).

22In contrast, epidermal knockout of BRI1 in pWER-BRI1-CRISPR lines resulted in both reduced meristem cell size and reduced mature cortex cell length (Fig. 4D–E), which is consistent with the reported role of epidermal brassinosteroid signaling and brassinosteroid-responsiveness across developmental zones of the epidermis in our scRNA-seq data.

23We measured the cell width of the tissue-specific CRISPR lines and found that pWER-BRI1-CRISPR lines have significantly wider cortex cells. However, despite their reduced cortex cell length, pCO2-BRI1-CRISPR lines did not have increased cell width (fig. S11D). This prompted us to perform 3D cell segmentations using MorphoGraphX. We found that the cortex cell volume was reduced in pCO2-BRI1-CRISPR roots but increased in pWER-BRI1-CRISPR roots (fig. S12, A and B). Brassinosteroids primarily control cell elongation rather than cell volume when signaling is affected across the entire root. On the other hand, the changes in volume observed in the tissue-specific CRISPR lines could be due to constraints between tissues or tissue-specific functions of brassinosteroids in influencing cell volume. These results indicate that in addition to the epidermis, brassinosteroid 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. This may explain the widening of cortex cells in pWER-BRI1-CRISPR lines.

24BRI1 was also reported to rescue bri1-T morphology when expressed in the developing phloem using the CVP2 promoter. However, 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 gene expression patterns 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. BRI1 driven by its native promoter was still present in the stele of our tissue-specific CRISPR lines when we observed phenotypic defects, suggesting that, unlike pCVP2-BRI1, native expression of BRI1 in the stele is not sufficient for brassinosteroid-induced cell elongation and root growth. These results confirm the role of the epidermis in brassinosteroid-regulated root growth and reveal the function of the cortex in brassinosteroid-mediated cell elongation, demonstrating how scRNA-seq can identify a spatiotemporal context for hormone signaling.

25To define a core set of genes associated with brassinosteroid 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 with development in wild-type cortex trajectories. The intersection of these three lists identified a core set of 163 brassinosteroid-responsive DEGs (Fig. 5A and Data S3). Consistent with regulation by brassinosteroids, 69% of the core DEGs are BES1 and BZR1 direct targets from ChIP experiments. Expression along cortex pseudotime illustrates induction by BL treatment and down-regulation in bri1-T (Fig. 5B). HAT7 and GTL1 were induced along these trajectories, suggesting a potential role for these transcription factors in controlling brassinosteroid-regulated gene expression in the cortex (Fig. 5C–E).

26To 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. 5F). We also observed HAT7 signals in the epidermis and endodermis, in line with expression patterns in our wild-type scRNA-seq atlas. HAT7 expression was decreased when brassinosteroid biosynthesis was inhibited with BRZ and restored upon BL treatment (fig. S13A).

27pGTL1-GTL1-mCitrine was more broadly expressed, with increasing levels in the cortex as cells progress from the transition domain to the elongation zone (Fig. 5G). GTL1-mCitrine expression was reduced by BRZ and increased by BL treatment (fig. S13B). These results confirm that brassinosteroids 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, suggesting that they may be part of the brassinosteroid-directed GRN activated as cells progress from proliferation to elongation.

28Previous studies inferred global or temporally resolved GRNs for brassinosteroid response, but they lacked cell type and developmental-stage specificity. To infer GRN configurations across our brassinosteroid time series, we used CellOracle (see methods), focusing on BL DEGs and associated transcription factors.

29Analysis of network importance scores, such as centrality measures, is a powerful approach to prioritize candidate regulators among DEGs. Since the cell wall signature peaked at 2 hours after BL treatment, we prioritized transcription factors with high network centrality scores in the elongating cortex at this time point. HAT7 was the top-ranked transcription factor in terms of degree centrality. HB13, HB20, and HB23 were also among the top 10 transcription factors (Fig. 6A–B and Data S4).

30We used CRISPR to generate hat7 loss-of-function mutants but did not observe phenotypes concerning cortex cell elongation (Fig. S14A–C). Together HAT7, HB13, HB20 and HB23 make up the alpha clade HD-ZIP I transcription factors. Since HB13, HB20 and HB23 are induced by brassinosteroids and predicted to regulate cell wall-related genes in our GRNs (Fig. 6B, fig. S15A–B and Data S5), 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. 6C–D and fig. S16A–C), providing evidence that HAT7 and its homologs are required for cell elongation. hat7 hb13 hb20 hb23 quadruple mutants also had wider mature cortex cells (fig. S16C) with increased volume compared to wild-type (fig. S17, A and B). Despite the decrease in final cell length, the root length of the quadruple mutant was not reduced (fig. S14A), suggesting that the decrease in cell length is at least partially compensated for by increased cell production. We found that hat7 hb13 hb20 hb23 quadruple mutants have an average of 15 more meristematic cortex cells compared to wild-type (fig. S16B), which is consistent with a compensatory increase in proliferation.

31We next investigated GTL1, the 5th highest-ranked transcription factor in the BL 2-hour elongating cortex GRN (Fig. 6A–B). Given that GTL1 was shown to function redundantly with DF1 in terminating root hair growth, we examined gtl1 df1 double mutants, finding shorter mature cortex cell lengths and shorter roots compared to wild-type (Fig. 6C–D and fig. S14A–C). DF1 was challenging to detect in scRNA-seq (fig. S15A) due to its low expression level. However, we observed increasing trends of DF1 expression along WOT trajectories in the BL time course, especially in cell wall responsive cortex cells (fig. S15B), which was verified using a pDF1-DF1-GFP reporter (fig. S15C).

32We examined the overlap between our CellOracle BL 2-hour elongating cortex GRN and a GRN reported by Clark et al. from bulk RNA-seq data in response to brassinosteroids. Only 43/31,330 edges are shared between the two networks. Moreover, HAT7, HB13, HB20, HB23, and GTL1 were among the top 10 transcription factors in terms of out-degree in the CellOracle elongating cortex GRN, which we validated through mutant analysis. On the other hand, none of these were among the top 100 transcription factors in the Clark et al. GRN (fig. S18, A and B), suggesting that they would have been difficult to prioritize from the bulk data alone. Together, our genetic analysis of HAT7 and GTL1 family transcription factors illustrates the power of GRN-mediated discovery of regulatory factors in spatiotemporal brassinosteroid response.

33Since BES1 is known to interface with other transcription factors in controlling brassinosteroid-regulated gene expression, we compared target genes for BES1 and BZR1 to ChIP targets of GTL1 and DF1. BES1 and BZR1 share 3,020 common targets with GTL1 and 2,490 common targets with DF1 (fig. S19A). When compared to brassinosteroid-regulated genes from scRNA-seq, BES1 and GTL1 targets showed the strongest enrichment in genes up-regulated by brassinosteroids in the transition domain and elongation zone of the cortex (fig. S19B), with 297 common targets of both BES1 or BZR1 and GTL1 being induced in the elongating cortex by BL treatment.

34Given the overlap between BES1 and GTL1 targets, we hypothesized that these transcription factors physically interact to regulate a common set of genes. Co-immunoprecipitation showed that GTL1-FLAG pulled down BES1-GFP (fig. S19C). These results suggest that brassinosteroids 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 brassinosteroid signal and/or to direct BES1 to drive tissue-specific gene expression by interacting with other more specifically expressed transcription factors.

35In addition to GTL1, we found that previously described BES1 interacting transcription factors including BRASSINOSTEROIDS AT VASCULAR AND ORGANIZING CENTER (BRAVO), BIM1, and MYB30, were differentially expressed in our scRNA-seq datasets in response to BL (Data S3). Among these, MYB30 was enriched in atrichoblasts of the epidermis and induced by BL treatment (fig. S20, A and B). CellOracle correctly predicted that MYB30 regulates LTPG2, a known target of MYB30 (Mabuchi et al., 2018), which also showed enrichment in atrichoblasts and temporal induction in response to brassinosteroids following MYB30 activation (fig. S20B). This supports the idea that cell type-specific patterns of brassinosteroid response may be at least partially explained by interactions between BES1 or BZR1 and other transcription factors.

36Our results indicate that gtl1 df1 mutants have reduced cortex cell elongation. On the other hand, gtl1 df1 mutants have longer trichoblasts. A downstream regulatory network that enables GTL1-mediated growth inhibition has been dissected in trichoblasts. 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 the gtl1 df1 double mutant. Using pseudobulk differential expression analysis, we detected relatively subtle changes in gtl1 or df1 single mutants compared to the wild type (fig. S21A–C). In contrast, 8,391 genes were differentially expressed in gtl1 df1 double mutants versus wild type (Fig. 7A).

371,077 genes were up-regulated across all developmental stages of the cortex of gtl1 df1, and 947 genes were down-regulated. The majority of cortex DEGs were affected in the elongation zone (Fig. 7, A and B, fig. S21C 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 brassinosteroid 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 brassinosteroid-induced genes in the cortex.

38Plotting gtl1 df1 DEGs along cortex pseudotime illustrated the down-regulation of several genes involved in cell elongation including CESA5 and AHA2 (Fig. 7C). These genes were enriched for the GO term “cell wall organization or biogenesis” (fig. S21D). 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. 7, D and E). 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. 7F and Movie S1). The reduced expression of cell wall-related genes in gtl1 df1 mutants validates our cell-type-specific brassinosteroid GRNs and identifies a function of GTL1 in promoting cortex cell elongation in response to brassinosteroids.

39Understanding 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. Plant hormones including brassinosteroids, auxin, gibberellins, and abscisic acid have been shown to exhibit tissue-specific responses, but how the associated GRNs are modulated in different cell types at particular developmental stages is enigmatic. In this study, we profiled brassinosteroid 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 brassinosteroid signaling, where brassinosteroids activate cell wall-related genes and promote elongation. We further showed that HAT7 and GTL1 are brassinosteroid-induced regulators along cortex trajectories that control cell elongation. These findings highlight the ability of single-cell genomics to identify context-specific transcription factors, a capability that could be leveraged to precisely engineer plant growth, development, and responses to stress. Our results reveal spatiotemporal brassinosteroid responses and the underlying GRNs.

40Arabidopsis accession Columbia-0 (Col-0) was used as a wild type. The following lines were previously described: bri1 GABI_134E10; bri1-116brl1brl3 triple mutant (bri1-T); pGL2-BRI1-GFP/bri1-T; gtl1-1 (WiscDsLox413-416C9), df1-1 (SALK_106258), and gtl1-1 df1-1; JKD-Ypet recombineering line. We produced hat7 single mutants and hat7 hb13 hb20 hb23 quadruple mutants using FASTRED multiplex CRISPR constructs containing an intronized version of Cas9.

41scRNA-seq experiments were performed as previously described. For each sample, ~0.5cm root tips from 7-day-old plants were harvested, and digested with protoplasting solution, and 16,000 cells were loaded on a 10X Genomics Chromium instrument, with the aim to capture 10,000 cells per sample with the 10X Genomics Chromium 3` Gene expression v3 or v3.1 kits.

42For BL scRNA-seq, we first grew plants on 1 μM BRZ to deplete endogenous brassinosteroids, then transferred plants to either a fresh BRZ plate or 100 nM BL. We performed two separate BL scRNA-seq treatment experiments. The first pilot experiment consisted of a BRZ and 2 hour BL treatment. The second experiment included two additional BRZ and BL 2 hours replicates and a single replicate of 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, 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. To test the effect of inhibiting endogenous brassinosteroid biosynthesis, we grew wild type on 1 μM BRZ or a mock DMSO control and performed two replicates of scRNA-seq spanning 30,962 cells for BRZ vs control analysis. Lastly, scRNA-seq was performed on wild type, gtl1, df1, and gtl1 df1 in duplicate under control conditions resulting in 74,810 scRNA-seq expression profiles.

43scRNA-seq data were aligned against the Arabidopsis TAIR10 reference genome. Quality filtering of cells was performed using the R package COPILOT (Cell preprOcessing PIpeline kaLlistO busTools). Downstream analysis were carried out using Seurat version 3.1.5, Waddington-Optimal Transport, muscat, tradeSeq, and GRNs inferred using CellOracle.

44To selectively block brassinosteroid signaling in cell types of interest, we performed tissue-specific CRISPR of BRI1 in a bri1 mutant complemented with pBRI1-BRI1-mCitrine. Two gRNAs targeting BRI1 were simultaneously expressed along with tissue-specific Cas9 to knock out BRI1 either in the epidermis and lateral root cap (pWER-BRI1-CRISPR) or in the cortex (pCO2-BRI1-CRISPR). For each root used for quantitative analysis, BRI1-mCitrine signal was acquired to confirm the efficiency of the tissue-specific knockout system. The supplementary materials provide all the methods details, including those summarized above.