Integrated omics networks reveal the temporal signaling events of brassinosteroid response in Arabidopsis

Brassinosteroids (BRs) are plant steroid hormones that regulate cell division and stress response. Here we use a systems biology approach to integrate multi-omic datasets and unravel the molecular signaling events of BR response in Arabidopsis. We profile the levels of 26,669 transcripts, 9,533 protein groups, and 26,617 phosphorylation sites from Arabidopsis seedlings treated with brassinolide (BL) for six different lengths of time. We then construct a network inference pipeline called Spatiotemporal Clustering and Inference of Omics Networks (SC-ION) to integrate these data. We use our network predictions to identify putative phosphorylation sites on BES1 and experimentally validate their importance. Additionally, we identify BRONTOSAURUS (BRON) as a transcription factor that regulates cell division, and we show that BRON expression is modulated by BR-responsive kinases and transcription factors. This work demonstrates the power of integrative network analysis applied to multi-omic data and provides fundamental insights into the molecular signaling events occurring during BR response.

When BRs are present, the BRI1/BAK1 complex activates a kinase signaling cascade resulting in the inactivation of BIN2 and the dephosphorylation of BES1/BZR1. This allows BES1/BZR1 to regulate target gene expression in the nucleus (15,40,49,57,(62)(63)(64). While there are many conserved GSK3 phosphorylation sites in BES1/BZR1 proteins, the exact sites that are phosphorylated by BIN2 and the sites responsible for negative regulation of BES1 activity are not well defined.
The role of BR in modulating cell division is welldocumented, particularly in the Arabidopsis root. It has been shown that BL has a dose-dependent effect on cell division in the root meristem: roots treated with higher levels of BL have more meristematic cells. Additionally, the bes1-D gain-of-function and bri1-116 loss-of-function mutants have altered cell cycle progression, which implicates that plant cyclins may be involved. Accordingly, it has been shown that CYCLIN D3;1 (CYCD3;1) is induced by BL and can rescue cell division defects in the bri1-116 mutant (14,19). BR has also been implicated in stem cell division and maintenance: roots treated with BL have excessive QC divisions and altered expression patterns of QC cell identity markers (14,30,53).
The BR signaling pathway in Arabidopsis is highly dependent on the levels of protein phosphorylation, modulation of protein levels, and downstream transcriptional regulation. Therefore, we determined the dynamic response to BR signaling in Arabidopsis by performing largescale transcriptome and (phospho)proteome profiling of seedlings treated with BL for different lengths of time. We then inferred a set of integrated, TF-centered Gene Regulatory Networks (GRNs) using our newly-developed Spatiotemporal Clustering and Inference of Omics Networks (SC-ION) pipeline. These networks illustrated how the phosphorylation state of TFs is important for predicting their downstream target genes. Additionally, SC-ION allowed us to infer one network per time point and visual-D R A F T ize the early and late (phospho)protein-transcript regulations in BR response. By combining these TF-centered networks with a correlation-based kinase signaling (i.e. kinase-centered) network, we illustrated the temporal cascade of BR response starting with kinase signaling and ending with differential transcript abundance. We were able to use this network to predict and experimentally validate previously uncharacterized BIN2 phosphorylation sites on BES1. Additionally, through network motif analysis, we identified a number of TFs putatively involved in BR response. In particular, we identified a C2H2-like TF whose mutants displayed hypersensitivity to BR, longer roots, and more cell divisions. The combination of this mutant's developmental phenotype as well as its putative role in BR response led us to name this TF BRON-TOSAURUS (BRON). Together our results provide a comprehensive guide to molecular signaling changes that occur in response to BR.

Results
Generating an integrative omics dataset of temporal BR response. To investigate the temporal response to BRs, we established a treatment system in which seedlings were sensitized to BRs by pre-treatment with 1 µM of the BR biosynthesis inhibitor brassinazole (BRZ) (3) for 7 days to reduce background BR signaling. The 7-day-old seedlings were then treated with a mock solvent or 1 µM BL for six different lengths of time (15 min, 30 min, 1 hr, 2 hr, 4 hr, 8 hr) (Fig 1A) . To confirm efficacy of the BL treatment, we assayed BES1 by western blot (Fig. S1). In BRZ treated seedlings, we found that BES1 predominantly exists in its phosphorylated form, while BL-treated seedlings showed an accumulation of dephosphorylated BES1 over time. Specifically, we observed an increase in the amount of dephosphorylated BES1 as early as 15 minutes after treatment, and the phosphorylated form of BES1 was undetectable by 1 hour after treatment (Fig.  S1). This demonstrates the expected BES1 response and thus the efficacy of BL treatment.
Based on the dynamic nature of BL response, we expected many transcripts, proteins, and phosphosites would have differential responses depending on the length of BL treatment. Thus, we performed multi-omics profiling of BL-treated and mock-treated seedlings at each of the six timepoints ( Fig 1A). We used 3' QuantSeq (35) to measure transcript levels and quantified protein abundance and phosphorylation level by performing two-dimensional liquid chromatography-tandem mass spectrometry (2D-LC-MS/MS) on Tandem Mass Tag (TMT) labeled peptides (18,32,(46)(47)(48).
To facilitate analysis of complex multi-run proteomics datasets, we constructed an analysis pipeline for quantitative proteomics data called TMT Normalization, Expression Analysis, and statistical Testing (TMT-NEAT). Our pipeline, which works on data generated from any organism, takes the TMT reporter ion intensity values (i.e. MaxQuant proteinGroups or PTM_Sites) and a metadata file containing sample information and TMT-labeling scheme as input. It first cleans the data by removing contaminants and appropriately labels the intensity data using the metadata file. Second, TMT-NEAT performs sample loading (within-run) and internal reference (between-run) normalization to eliminate batch effects (42). Third, it provides multiple qualitative plots such as hierarchical clustering and principal component analysis to visualize differences between biological groups. Finally, it performs differential expression analysis on the normalized values using a user-supplied p-or q-value threshold ( Fig 1B). TMT-NEAT is publicly available and can be run through an RShiny Graphical User Interface (GUI) (see Methods).
Using these methods, we identified 32,549 transcripts, 9,035 protein groups, and 26,950 phosphosites (arising from 5,648 phosphoproteins) across the six timepoints ( Fig 1C, Datasets S1-S3). We found that the number of differentially expressed (DE) transcripts, proteins, and phosphosites varied depending on the time point. When we examined the transcript data, we found that only 208 transcripts (17 up, 191 down) were DE in response to BL within the first 15 minutes, whereas 454 protein groups (214 up, 250 down) and 590 phosphosites (237 up, 353 down) were DE at this same time point. In addition, many more transcripts (2,653 total: 1,247 up, 1,406 down) were DE beginning at 30 minutes, which led us to speculate that the early BR response is predominantly post-transcriptional. This is supported by the role of BES1/BZR1, which must be dephosphorylated in order to enter the nucleus and transcriptionally regulate downstream genes in response to BR (16,62,63).
We next performed Gene Ontology (GO) analysis on the DE transcripts, proteins, and phosphosites at early (1 hour or earlier) and late (after 1 hour) timepoints (Fig S1, Dataset S4). We found that multiple BR response terms are enriched at different time points and for different DE gene-products. For example, BR mediated signaling pathway and cellular response to BR stimulus terms are enriched in the DE phosphosites at both early and late time points. Additionally, BR biosynthetic process and BR homeostasis are enriched in the DE transcripts at early timepoints, and the term for response to brassinosteroid is enriched in both DE transcripts and phosphosites. We also see enriched terms that have been linked to BR response such as cellular response to stress, response to auxin mediated signaling pathway, and defense response (38). Taken together, our omics profiling captures how different gene products temporally respond to BR treatment in Arabidopsis.
Predictive networks illustrate the temporal cascade of BR response. We have previously shown that integrating mRNA, protein, and phosphorylation data sets greatly  improves the predictive power of reconstructed gene regulatory networks (GRNs) (54). However, integrating these multi-omics data types remains challenging. Thus, we next developed a network inference pipeline, which we named Spatiotemporal Clustering and Inference of Omics Networks (SC-ION), to integrate any number of different types of expression profiles into one cohesive, predictive GRN ( Fig. 2A, see Methods). SC-ION builds on our MATLAB-based pipeline Regression Tree Pipeline for Spatial, Temporal, and Replicate data (RTP-STAR) (7,65), which is an adaptation of the GENIE3 (20) network inference method and functions only on transcriptomic datasets. In SC-ION, we further improve on RTP-STAR by: 1) incorporating Dynamic Time Warping (DTW) clustering for temporal data (13) and Independent Component Analysis (ICA) clustering for non-temporal data (37); 2) allowing the user to provide separate regulator and target matrices for integration of DE gene-products; 3) integrating any number of different types of expression profiles into one GRN; and 4) providing our pipeline as an RShiny GUI ( Fig. 2A).

D R A F T
Here, we used SC-ION to infer two separate transcription factor (TF)-centered networks of BR response. In these TF-centered GRNs, TFs serve as "regulators" used to infer their "target" genes. In the first network, which we call the abundance network (blue, solid edges, Fig 3), TF protein abundance (when quantified) or TF transcript abundance (when cognate protein was not quantified) was used as the "regulator" value to infer their "target" transcript abundance. In the second phosphosite network (green, dashed edges, Fig 3), we inferred the "target" transcript abundance using TF phosphosite intensities as the "regulator" value. For each of these networks, we took advantage of our temporal data by 1) clustering the gene-products using DTW to create temporally-informed regulatory modules and 2) inferring one network per time point to visualize the regulations unique to each time point (Fig 2B). When connecting the time-pointspecific networks, we found that there were distinct clusters of early and late predictions. In addition, there is cross-communication between the time points, where early regulators feed-forward into late regulators, and conversely where late regulators feed-back onto early regulators ( Fig 2B).
In addition to our SC-ION-generated TF-centered GRNs, we used our correlation-based approach to infer a kinase signaling network (purple, dotted edges, Fig 3) (55). In this network, we considered kinases with DE phosphosites in their p-loop (also termed activation loop) domain as potential regulators, as these phosphorylated kinases should be active (1) and thus useful for predicting kinase-signaling (4,44,55). In agreement with our previous work in maize (55), we observed that the correlation between kinase protein abundance and kinase p-loop phosphorylation intensity (i.e. activation state) greatly varies depending on the time point (Fig S3), motivating our use of p-loop site intensities rather than simply kinase abundance in our network. This "kinase-centered" network complements our TF-centered GRNs by predicting kinase-dependent signaling events.
When we merged our kinase signaling network with our TF-centered networks, we found that genes with known roles in the BR response pathway (39) were significantly enriched in the list of network regulators (Hypergeometric test, p<0.001). Additionally, we noticed that this merged network illustrates the temporal cascade of BR response across these different omics levels (Fig. 3, Dataset S4). Our inferred network places the kinase signaling interactions (purple) towards the top of the network (early in time).
Next are the TF phosphosite-level regulations (green), followed by the TF abundance-level regulations (blue). Thus,

D R A F T
our network predicts that BR response begins with kinase signaling, followed by transcriptional regulation via modulation of TF phosphorylation and/or abundance. This network prediction is in agreement with what we currently know about BR signaling, which begins with the kinases BRI1 and BAK1 initiating a series of (de)phosphorylation events leading to the regulation of downstream target genes by TFs including BES1 and BZR1. Our timepointspecific networks further illustrate which regulations are predicted to occur in early or late BR response ( Fig 2B).
Integration of kinase-signaling and TF-centered networks reveals the BIN2-BES1 signaling cascade in response to BR. We extracted the first-neighbor network of BES1 to illustrate how our multi-omics profiling and integrative network inference can elucidate the signaling events that occur in response to BR (Fig. 3). This network contains kinases that are predicted to be upstream of BES1 phosphorylation, TFs that are predicted to transcriptionally regulate BES1, and downstream targets of BES1. Our kinase-signaling network predicted that Serine 179 and Serine 180 (S179, S180) on BES1 are downstream of known BR-responsive kinases such as BIN2, BAK1, and BSK1 (16, 26, 50) ( Fig S3). In addition, while not present in our kinase signaling network, S171 differentially accumulates in our BR timecourse (Dataset S3), and all three sites are largely conserved in BES1 and its homologs ( Fig S3). Thus, we mutated these three sites from serine to alanine (phosphorylation null) and tested BES1 phosphorylation status in Nicotiana Benthamiana. While more than half of the wild-type BES1-FLAG is phosphorylated one day after transfection, BES1 3SA -FLAG exists mainly in dephosphorylated form ( Fig S3). When co-expressed with BIN2, the phosphorylated form of wild-type BES1-FLAG was increased, while phosphorylated form of BES1 3SA -FLAG was increased to a lesser extent, which supports the conclusion that S179, S180 and S171 are BIN2 phosphorylation sites ( Fig S3).
In order to elucidate the functional consequences of the phosphorylation on S179, S180 and S171, we transformed BES1-FLAG and BES1 3SA -FLAG to bri1-301, a weak allele of BR receptor BRI1 loss-of-function mutant, and assayed phosphorylation status in response to 100 nM BL. We found that BES1 3SA -FLAG has lower levels of dephosphorylated BES1, and the amount of dephosphorylated BES1 increased more in response to BL treatment in BES1 3SA -FLAG than in BES1-FLAG plants. Consistent with these observations, 18 out of 24 (75%) 3-week-old T1 plants overexpressing BES1 3SA -FLAG showed stronger gain-of-function BR mutant phenotype, with elongated leaf petioles and curly leaves, compared to the 18 bri1-301 plants overexpressing wild-type BES1-FLAG which did not show this phenotype ( Fig S3). Further, it has been shown that the mutation of S173 in BZR1 (equivalent to S171 in BES1) alone did not change the phosphorylation status of BZR1 (12). Thus, these results support that BIN2 phosphorylates BES1 at S179, S180 and S171 to inhibit its function in BR-regulated signaling.
Downstream of kinase signaling, BES1 is predicted to regulate different downstream targets depending on whether its phosphosite levels (green, dashed) or transcript abundance (blue, solid) is used in the SC-ION pipeline. We mined Chromatin ImmunoPrecipitation (ChIP)-chip and ChIP-Seq datasets on BES1 and its homolog, BZR1 (40,49,64), and found that 88 out of 144 (61%) predicted first-neighbor targets of BES1, in our network, are directly bound by BES1 or its homolog BZR1 (Fig  3, orange circles; Table S1). Some of these validated genes with known roles in BR response include IAA3/SHORT HYPOCOTYL 2 (SHY2) (29), XYLOGLU-CAN:XYLOGLUCOSYL TRANSFERASE 33 (XTH33) (49,67) and SMALL AUXIN UP RNA 26 (SAUR26) (40). Thus, we were able to use our integrative omics network to identify previously unknown BIN2 phosphorylation sites on BES1 as well as validate its putative direct downstream targets identified by ChIP.
Network motif analysis predicts TFs involved in BR response. We next leveraged the network prediction to identify candidate genes involved in mediating the response to BR. We used the Network Motif Score (NMS) (7,65) to classify genes in the TF-centered networks based on their presence in certain biological motifs, such as feedback and feed-forward loops. Genes with higher NMS scores have been shown to have a more important role in the biological process of interest (2,7,21,34,65). Accordingly, BES1 had the 25th highest NMS score (top 5%), illustrating that we could use the NMS to identify BR response regulators. We chose three TFs, ANTHOCYAN-LESS 2 (ANL2), TCX2, and BRON (AT1G75710), that had high (all in the top 35%) NMS scores in either the TF abundance or phosphosite GRNs (Dataset S5). We then examined subnetworks for each of these TFs, starting with kinase signaling and ending with transcriptional regulation.
SC-ION predicts that our first TF of interest, TCX2, and BES1/BZR1 HOMOLOG 2 (BEH2) regulate each other in a feedback loop (Fig S4A). Thus, we reasoned that TCX2 may regulate cell division in response to BR and treated the tcx2-2 and tcx2-3 mutants with 100 nM BL. In WT plants, the addition of BL causes a dramatic reduction in the root length. However, we did not find a significant difference in BL response in either of the tcx2 mutant alleles compared to WT (Fig S4, Table S2).
We then focused on the subnetwork for ANL2, which predicts that ANL2 and BES1 regulate each other in a feedback loop (Fig S4). To test if ANL2 may be involved in BR response, we treated WT, anl2-2, and anl2-3 plants with 100 nM BL. We found that anl2 mutant roots shorten more than WT when treated with BL, demonstrating that the anl2 mutants are hypersensitive to BL (Fig S4, Table S2). It has also been shown that BL can induce excessive QC divisions in the root (14,30,53). Thus, we checked the anl2 mutant roots for QC divisions as a secondary BR response phenotype: however, we found that anl2 mutant roots do not display excessive QC divisions (Fig S4).

D R A F T
The TF BRONTOSAURUS (BRON) regulates cell division in response to BR. Our last TF of interest, BRON, is predicted to be regulated by TFs whose phosphorylation is dependent on multiple BR-signaling kinases such as BAK1, BSKs, and SERKs in our network (Fig 4A). We obtained a weak (bron-1) and a strong (bron-2) allele for BRON (Fig S5) and examined their BR response as well as their root development (Fig 4B-E). We found that roots from both alleles shorten more in response to BL than wild type, demonstrating that bron mutant roots are hypersensitive to BL treatment ( Fig 4B, Table  S3). Additionally, we observed that bron mutants have longer roots with more meristematic cells (Fig S5), as well as significantly more divisions in the QC (Fig 4C-E). bron-2 also displays excessive columella divisions and a higher number of undivided (i.e. actively dividing) CEI, while bron-1 only shows excessive endodermis divisions, potentially due to its weaker effect on BRON expression (Fig 4E). This led us to hypothesize that BRON could regulate cell division in response to BR.

D R A F T
To determine the transcripts modulated by BRON, we performed RNA-seq on root tissue from the bron-2 mutant (Dataset S7). We found that a range of BR-responsive genes were enriched in this mutant, particularly those genes repressed by BL after 15 minutes and induced by BL after 1 or 2 hours (p < 0.001) (Dataset S8). Further, we found that most of the genes have lower expression in the bron-2 mutant, suggesting that BRON transcriptionally activates these genes (Fig S6).
Given the role of cyclins in regulating cell division (56), we specifically examined the expression of cyclins in the bron-2 mutant. We found that four cyclins are differentially expressed in the bron-2 mutant: three are repressed by BRON (CYCD3;1, CYCP4;1, CYCP4;2). Further, two of these cyclins, CYCD3;1, and CYCP4;1, are significantly induced by BL at 4 and 8 hours after treatment ( Fig 5A). Importantly, it has been shown that CYCD3;1 is induced by BL and contributes to BL-regulated cell division (14,19). We also examined the cell-type-specific expression of CYCD3;1, CYCP4;1, and BRON in the root stem cells and mature root cells (7,28), and we found key differences in where CYCD3;1 and CYCP4;1 are co-expressed with BRON. For example, CYCP4;1 is coexpressed with BRON only in the QC. In contrast, BRON and CYCD3;1 are co-expressed in the Epidermis/Lateral Root Cap (Epi/LRC) initials and the mature columella ( Fig  5B). Taken together, these results suggest that BR incudes the expression of CYCD3;1 and CYCP4;1, and therefore cell division, through the repression of BRON (Fig 5C).

Discussion
Here, we used a systems biology approach to unravel the temporal response to BR in Arabidopsis. By generating and integrating omics datasets, we were able to quantify how transcript, protein, and phosphorylation levels change in response to BR over time. We found that most of the (phospho)protein response occurs in the earlier timepoints (before 1 hour), whereas there are sets of early-and late-responsive transcripts (Fig 1). This suggested that BR first triggers the phosphorylation of TFs which then go on to regulate transcripts at later timepoints. Our predictive network reconstructed using SC-ION corroborated this hypothesis, illustrating how BR-responsive kinase signaling leads to phosphorylation of TFs and downstream transcriptional regulation ( Fig  3). Further, our time-point specific networks allow us to predict which regulations occur early and late in BR response (Fig 2).
One well-known example of BR response involves the TFs BES1/BZR1, which are de-phosphorylated in response to BR, allowing them to move to the nucleus and transcriptionally regulate downstream genes. It is well-established that the GSK3 kinase BIN2 phosphorylates BES1/BZR1 to negatively regulate their functions through the promotion of the cytoplasmic retention, inhibition of DNA binding and protein degradation of BES1/BZR1 (39). Although there are more than 20 predicted sites, the exact BIN2 phosphorylation sites on BES1 are not functionally characterized. Through our omics profiling and network analysis, we predicted that S171, S179 and S180 in BES1 are potentially phosphorylated by BIN2. Further mutational analysis indicated that the phosphorylation of the three sites contributed to BES1 phosphorylation by BIN2 and are important for BES1 activity. In contrast, the mutation of S173 in BZR1 (equivalent to S171 in BES1) affected BZR1 interaction with 14-3-3 and nuclear localization but did not change its phosphorylation status (12). Our study therefore identified three key BIN2 phosphorylation sites that negatively regulate BES1 activity and hence BR signaling (Fig S3).
We further used our network to predict novel roles for the TFs ANL2 and BRON in BR response (Figs 4 and 5, Fig  S4). We found that mutants of both TFs are hypersensitive to BL treatment, but only bron mutants display excessive cell divisions. Given anl2's hypersensitivity to BL, it would be interesting to investigate its role in BR response in the future, given that it likely has a different role than BRON. We additionally found a significant overlap between genes differentially expressed in the bron-2 mutant and genes induced or repressed by BR. Specifically, we observed a significant overlap in disease resistance, drought response, and stress response genes which were repressed by BL 15 minutes after treatment and repressed in the bron-2 mutant (Fig S6). Crosstalk between BR and drought response is well known and involves other TFs such as WRKY46/54/70, RD26 and TINY (6, 38, 59, 61), but the temporal aspects of this signaling have not been previously examined. Our observation that BRON is regulated by BLs as early as 15 minutes after treatment suggests that inhibition of stress responses quickly follows activation of the BR pathway.
We found that cyclin genes, specifically CYCD3;1 and CYCP4;1, are induced by BL and repressed by BRON, suggesting that BR induces these cyclins through its repression of BRON (Fig 5C). These results are supported by multiple studies which elucidate the role of CYCD3;1 in BR-induced cell division (14,19) and describe how BR specifically induces QC division in the root (14,30,53 Expression in root stem cells Expression in mature root cells

D R A F T
mature columella. We found that bron-2 mutants have excessive columella divisions (Fig 4E), suggesting that perhaps BRON represses cell division in the columella and lateral root cap through CYCD3;1. It would be of interest in the future to unravel BRON's cell-type-specific repression of cell division in response to BR.
Taken together, our temporal, integrated omics dataset, kinase-signaling and TF-centered networks can be used as a resource to identify additional genes like BRON which are implicated in BR response in Arabidopsis.

Plant Materials and Growth Conditions.
The Arabidopsis accession Columbia-0 was used as the wild type control in all experiments. The T-DNA insertion mutants tcx2-2 (SAIL_808_H08), tcx2-3 (SALK_021952), anl2-2 (SALK_000196C), and anl2-3 (SAIL_418_C10) were described previously (7). The bes1-D and bri1-301 mutants were also described previously (11,53,60,63 For bron-2 transcriptional profiling, total RNA was isolated from approximately 2 mm of five-day-old Col-0 and bron-2 root tips using the RNeasy Micro Kit. cDNA synthesis and amplification were performed using the NEBNext Ultra II RNA Library Prep Kit for Illumina. Libraries were sequenced on an Illumina HiSeq 2500 with 100 bp single-end reads. Reads were mapped to the TAIR10 genome using Cufflinks (51). Differential expression was performed using PoissonSeq with a p-value cutoff of 0.05. Raw sequencing data are deposited at the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE157000, accession token ivabyiuypzwxvil).
Protein Extraction and Digestion. The proteomics experiments were carried out based on established methods (46)(47)(48). Protein was extracted from aliquots of the tissue used for transcriptome profiling and digested into peptides with trypsin and Lys-C using the phenol-FASP method detailed in (46,48). Resulting peptides were desalted using 50 mg Sep-Pak C18 cartridges (Waters), dried using a vacuum centrifuge (Thermo), and resuspended in 0.1% formic acid. Peptide amount was quantified using the Pierce BCA Protein assay kit.
Tandem Mass Tag (TMT) Labeling. The TMT labeling strategy used in this experiment is provided in Table  S3. 45 µg of peptides were taken from each individual sample, pooled, and then split into two pooled references. TMT10plex TM label reagents (ThermoFisher, Lot UD280154) were used to label 200 µg of peptides, from each sample or pooled reference, at a TMT:peptide ratio of 0.2:1 as described in (48). After 2 hours incubation at room temperature the labeling reaction was quenched with hydroxylamine. Next, the ten samples were mixed together, an aliquot of 75µg of peptides was reserved for protein abundance profiling, and the remaining peptides were used for phosphopeptide enrichment and stored at -80ºC.

D R A F T
Phosphopeptide Enrichment. The TMT-labeled phosphopeptides were first enriched using the High-Select TiO 2 Phosphopeptide Enrichment Kit (Thermo) using the manufacturers protocol. The High-Select Fe-NTA Phosphopeptide Enrichment Kit (Thermo) was then used on the flowthrough from the TiO 2 enrichment to enrich additional phosphopeptides. The manufacturers protocol for the Fe-NTA kit was used except the final eluate was re-suspended with 50 µL 0.1% formic acid. The eluates from the TiO 2 and Fe-NTA enrichments were combined and stored at -80ºC until analysis by LC/MS-MS.

LC/MS-MS.
An Agilent 1260 quaternary HPLC was used to deliver a flow rate of 600 nL min-1 via a splitter. All columns were packed in house using a Next Advance pressure cell, and the nanospray tips were fabricated using a fused silica capillary that was pulled to a sharp tip using a laser puller (Sutter P-2000). 25 µg of TMT-labeled peptides (non-modified proteome), or 10 µg TiO 2 enriched peptides (phosphoproteome), were loaded onto 20 cm capillary columns packed with 5 µM Zorbax SB-C18 (Agilent), which was connected using a zero dead volume 1 µm filter (Upchurch, M548) to a 5 cm long strong cation exchange (SCX) column packed with 5 µm PolySulfoethyl (PolyLC). The SCX column was then connected to a 20 cm nanospray tip packed with 2.5 µM C18 (Waters). The 3 sections were joined and mounted on a Nanospray Flex ion source (Thermo) for on-line nested peptide elution. A new set of columns was used for every sample. Peptides were eluted from the loading column onto the SCX column using a 0 to 80% acetonitrile gradient over 60 minutes. Peptides were then fractionated from the SCX column using a series of 18 and 6 salt steps (ammonium acetate) for the non-modified proteome and phosphoproteome analysis, respectively. For these analyses, buffers A (99.9% H2O, 0.% formic acid), B (99.9% ACN, 0.1% formic acid), C (100 mM ammonium acetate, 2% formic acid), and D (2 M ammonium acetate, 2% formic acid) were utilized. For each salt step, a 150-minute gradient program comprised of a 0-5 minute increase to the specified ammonium acetate concentration, 5-10 minutes hold, 10-14 minutes at 100% buffer A, 15-120 minutes 10-35% buffer B, 120-140 minutes 35-80% buffer B, 140-145 minutes 80% buffer B, and 145-150 minutes buffer A was employed.
Eluted peptides were analyzed using a Thermo Scientific Q-Exactive Plus high-resolution quadrupole Orbitrap mass spectrometer, which was directly coupled to the HPLC. Data dependent acquisition was obtained using Xcalibur 4.0 software in positive ion mode with a spray voltage of 2.10 kV and a capillary temperature of 275°C and an RF of 60. MS1 spectra were measured at a resolution of 70,000, an automatic gain control (AGC) of 3e6 with a maximum ion time of 100 ms and a mass range of 400-2000 m/z. Up to 15 MS2 were triggered at a resolution of 35,000 with a fixed first mass of 120 m/z for phosphoproteome and 115 m/z for proteome. An AGC of 1e5 with a maximum ion time of 50 ms, an isolation window of 1.3 m/z, and a normalized collision energy of 33. Charge exclusion was set to unassigned, 1, 5-8, and >8. MS1 that triggered MS2 scans were dynamically excluded for 45 or 25 s for phospho-and non-modified proteomes, respectively.
Carbamidomethyl cysteine was set as a fixed modification while methionine oxidation and protein N-terminal acetylation were set as variable modifications. For the phosphoproteome "Phospho STY" was also set as a variable modification. The sample type was set to "Reporter Ion MS2" with "10plex TMT selected for both lysine and N-termini". TMT batch-specific correction factors were configured in the MaxQuant modifications tab (TMT Lot UD280154). Digestion parameters were set to "specific" and "Trypsin/P;LysC". Up to two missed cleavages were allowed. A false discovery rate, calculated in MaxQuant using a target-decoy strategy (Elias and Gygi, 2007), less than 0.01 at both the peptide spectral match and protein identification level was required. The 'second peptide' option identify co-fragmented peptides was not used. The match between runs feature of MaxQuant was not utilized. Raw proteomics data have been deposited on MassIVE and can be accessed at the link ftp://MSV000085606@massive.ucsd.edu with reviewer username "MSV000085606_reviewer" and password "Clark_BL".
Statistical analysis was performed using TMT-NEAT Analysis Pipeline version 1.3 (https://github.com/nmclark2/TMT-Analysis-Pipeline). This pipeline takes the "proteinGroups" (protein abundance) or "Phospho(STY)Sites" (phosphoproteome) tables output from MaxQuant as well as a metadata file detailing the TMT labeling scheme and sample information as input. Example input files are provided in the GitHub repository. First, the MaxQuant output table is trimmed to only include the needed information for statistical analysis, and the columns are re-labeled using the provided sample information. Contaminants are removed at this stage. Next, data are normalized using the sample loading normalization and internal reference normalization methods such that samples can be compared across runs (42). Quantitative plots such as boxplots, hierarchical clustering, and principal components analysis are provided for quality control. Finally, statistical analysis is performed using PoissonSeq (27) and histograms of p-or q-value dis-D R A F T tributions are generated. Proteins and phosphosites were categorized as differentially accumulating if they had a pvalue < 0.05 and fold-change > 1.1.
GO Analyses. GO analysis on the DE transcripts, protein, and phosphoproteins was performed using PAN-THER (33). Genes were separated depending on whether they were induced or repressed by BL in early (1 hr or prior) or late (after 1 hr) time points. Biological process GO terms were considered significantly enriched if they had a corrected p-value 0.05 (Dataset S4).
Gene Regulatory Network Inference. TF-centered gene regulatory networks were inferred using SC-ION version 2.0 (https://github.com/nmclark2/SCION). SC-ION builds on the RTP-STAR pipeline (7,65) by incorporating DTW and ICA clustering (13,37) and integration of different data types. SC-ION uses an adapted version of GENIE3 (20) which allows for separate regulator and target data matrices (54). This allows the user, for example, to use protein abundance data for regulators (TFs) and transcript data for targets (all genes). SC-ION takes regulator and target lists and regulator and target data matrices as input. In addition, SC-ION takes a clustering data matrix which can be different from the regulator and target data matrices. This allows the user to cluster genes based on different datatypes. A version of SC-ION without this clustering step is also available. SC-ION outputs a table of the predicted regulations as well as a weight for each edge, where a higher weight indicates higher confidence in that inferred edge (20). This table can be imported into software such as Cytoscape (45) for network visualization. Test input files for SC-ION are provided on GitHub.
Two TF-centered networks were inferred. In the first network, TF protein abundance (when quantified) or TF transcript abundance (when cognate protein was not quantified) was used as the "regulator" value to infer their "target" genes' transcript abundance. In the second phosphosite network, we inferred the transcript levels of "target" genes using TF phosphosite intensities as "regulators." In both networks, the clustering matrix was constructed by combining the "regulator" (TF abundance) and "target" (transcript levels) data so that genes were clustered based on the protein/phosphosite levels of the regulators and transcript levels of the targets. Clustering was then performed using the DTW method given the temporal nature of our data. Six sub-networks were inferred for each TF-centered network, where each sub-network represents the regulations predicted to happen using only one time point (15 min, 30 min, 1 hr, 2 hr, 4 hr, 8 hr). To achieve this, only the genes DE at each time point were used in these individual subnetworks. This allows us to denote which regulations are predicted to happen early or late in BR response. The six subnetworks were merged in Cytoscape to form the final abundance and phosphosite networks. Cytoscape was also used to create the merged abundance, phosphosite, and kinase signaling networks (Dataset S5).
We have previously used the Normalized Motif Score (NMS) to predict biologically important genes in GRNs from Arabidopsis (7,65). Four different motifs were used to calculate the NMS for the merged abundance and phosphosite networks separately: feed-forward loops, feedback loops, diamond motifs, and bi-fan motifs. First, the number of times a gene appeared in each motif was counted using the NetMatchStar app (43) in Cytoscape. Then, the counts were normalized to a scale from 0 to 1 and summed to calculate the NMS for each gene.
Activation loop domains, also called p-loop domains, in protein kinases were identified using a modified version of the pipeline described in (55). Briefly, all 35,386 protein sequences available in the TAIR10 annotation were searched for kinase domains using The National Center for Biotechnology Information batch conserved domain search tool (31). From this list of 1,522 proteins with identified kinase domains, 878 were also annotated with activation loop (p-loop) coordinates by the search tool. The kinase domains of proteins lacking the p-loop coordinates were aligned using MAFFT (22). The resulting alignment was manually searched for the well conserved p-loop beginning (DFG) and end (APE) motifs. An extra 482 p-loop coordinates were obtained, for a total of 1,360 protein kinases with p-loop coordinates.
The kinase signaling regulatory network was inferred using a previously described correlation-based approach (55). Kinases with phosphosites in the p-loop domain that were DE in response to BL were used as the potential regulators. All genes with phosphosites that were DE in response to BL were used as the potential targets. We used phosphosite intensity rather than abundance for the regulators as it has been shown that phosphosite intensity has greater predictive power (55) (Fig S2). Pearson and Spearman correlations were calculated for each regulatortarget pair, and edges were kept for those pairs with Pearson correlation 0.5 or Spearman correlation 0.6.

Mutant BES1 Cloning and Protein Level Detection.
The three Serine to Alanine mutations were introduced by two-step PCR using the primers provided in Table  S4. Two BES1 fragments were generated and combined to form the full-length mutant BES1. The full-length mutant BES1 was then cloned to Pro35S:FLAG vector to generate Pro35S:BES1 S3A -FLAG. Pro35S:BES1-FLAG was subcloned from Pro35S:BES1-GFP (63).
Agrobacterium containing Pro35S:YFP-BIN2 was used for co-expression. Leaf discs were collected 24 hours after infiltration and flash frozen in D R A F T liquid nitrogen. The samples were ground in 2xSDS buffer and resolved on SDS-PAGE. Anti-Flag antibody was used for western blotting.
To test BES1 phosphorylation status in response to BR, T1 seeds from transgenic plants overexpressing Pro35S:BES1-FLAG or Pro35S:BES1 S3A -FLAG were germinated on 1/2MS plates containing 1uM BRZ and 50mg/ml Kanamycin. Ten-day-old seedlings were collected and subjected to 100nM BL treatment. DMSO was used as control. Twelve seedlings were used for each treatment. The samples were ground in 2xSDS buffer and resolved on SDS-PAGE. Anti-Flag antibody was used for western blotting.
BR Phenotyping Methods. BL phenotyping were carried out as previously described (58). Seeds were sterilized for 4 hours in a Nalgene Acrylic Desiccator Cabinet (Fisher Scientific, 08-642-22) by mixing 200mL bleach (8.25% sodium hypochlorite) with 8mL concentrated hydrochloric acid to generate chlorine gas. Seeds were then resuspended using 0.1% agarose solution for plating. Control (BL0; DMSO solvent only) or BL100 treated (100nM Brassinolide; BL, Wako chemicals) were plated on ½ LS plates supplemented with 1% (w/v) sucrose. After seeds were plated, the plates were sealed with breathable tape (3M Micropore) and placed in the dark at 4°C for 5 days for stratification. Plants were grown for 7 days at 22°C under continuous light. Plates were imaged with an Epson Perfection V600 Flatbed Photo scanner at a resolution of 1200 DPI and root length was then measured in ImageJ.
bron Mutant Genotyping and Root Phenotyping. RT-qPCR was performed on bron-1 and bron-2 mutant alleles as described in (7) to measure BRON expression. Total RNA was isolated from approximately 2mm of five-day-old Col-0, bron-1 and bron-2 root tips using the RNeasy Micro Kit (Qiagen). qPCR was performed with SYBR green (Invitrogen) using a 7500 Fast Real-Time PCR system (Applied Biosystems) with 40 cycles. Data were analyzed using the Ct (cycle threshold) method and normalized to the expression of the reference gene UBIQ-UITIN10 (UBQ10). qPCR was performed on two technical replicates of three independent RNA samples (biological replicates). Primers used for qPCR are provided in Table  S4.
Confocal imaging was performed on a Zeiss LSM 710. Cell walls were counterstained using propidium iodide (PI). The number of meristematic cells and cell divisions were manually counted.

Statistics.
A generalized mixed linear model with penalized quasi-likelihood (glmmPQL in R) was used to determine the genotype x treatment interaction p-values for the BL phenotyping experiment. In this model, the genotype and treatment effects were considered fixed, and the experiment date and plate number were considered random effects with a Gaussian error distribution. Hypergeometric and/or Chi-squared tests were used for test for enrichment in the RNAseq DE gene lists. For bron mutant root phenotyping, a two-tailed Wilcoxon test was used for statistical significance as some of the data did not follow a normal distribution. To select p-and q-value cutoffs for the large-scale omics experiments, we used the distribution of p-and q-values generated from the statistical tests as described in (5,8).