The Nipponbare (Oryza sativa L. ssp japonica) plants used in this study were grown by hydroponics method in growth chambers (90% relative humidity, 30/28°C, 14 h light/10 h dark cycle). Two-week old young seedlings with hydroponics were soaked into water containing 10 μM epibrassinolide (Cat No. E1641, Sigma, St. Louis, MN, USA) for 24 hours, during which samples were collected at different time points (0 h, 3 h, 6 h, 12 h, and 24 h) and immediately stored in liquid nitrogen until use. Three biological replicates were performed for each treatment and the control.

Total RNA of the samples after BR treatment was isolated using Trizol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s manual. Four micrograms of total RNA was performed for reverse transcription using first strand cDNA synthesis Kit (Toyobo, Shanghai, China). The primer pairs were set to detect the expression level of ILI1 (LOC_Os04g54900) (forward: 5′ ATGTCGAGCAGCCGGAGGTC 3′, reverse: 5′ CGTCTCGCTGAGGTTGTCC 3′), BUI1 (LOC_Os06g12210) (forward: 5′ CGACGACGAAGCTGCTGA 3′, reverse: 5′ CGCCTGGGCTGTTGTGAT 3′), IBH1 (LOC_Os04g56500) (forward 5′ CCGCCGAACCCTAACCCTAG, reverse 5′ CAGGAAGTGGAAGGCCAGCAT) and ubiquitin gene was used as an internal control (primer forward 5′ CACCCTGGCTGACTACAACA 3′, reverse 5′ TTCTTCTTGCGGCAGTTGAC 3′). Real-time quantitative RT-PCR was performed in a total reaction volume of 20 ml (10 μL THUNDERBIRD SYBR^® qPCR Mix (Toyobo, Shanghai, China), 1 μL cDNA, 1 μL primers, and 8 μL water) on the Bio-Rad CFX96 real-time PCR detection system (Bio-Rad, Hercules, CA, USA). The relative expression level was calculated by the 2^-ΔΔCT method. The experiment was performed in three biological replicates.

The total protein extraction and digestion were performed by exactly following the methods described by (Hou et al., 2015). For phosphopeptides enrichment, 1 mg digested peptides were dissolved with binding buffer (80% ACN, 5% TFA, 1 M lactic acid), then incubated with 4 mg TiO₂ beads (GL sciences, Torrance, CA, USA) for three times, each time for 30 minutes at room temperature. The totally 12 mg beads were then washed with 500 μL binding buffer for twice and 500 μL washing buffer (80% CAN, 0.5% TFA) for twice. All TiO₂ beads were transferred into a 200 μL homemade StageTip containing two layers of C18 solid phase extraction disk (3M, St. Paul, MN, USA). The StageTip was centrifuged at 300 g for 10 min to discard the flow through. The enriched phosphopeptides were incubated with 150 μL elution buffer (40% ACN, 15% NH₃H₂O) for four times. The combined 600 μL eluates were subsequently dried to ∼5 μL in a SpeedVac and reconstituted with 5% MeOH in 1% TFA solution for LC-MS/MS analysis.

This was done in Beijing Proteome Research Center. LC-MS/MS analyses were performed on an Easy-nLC 1000 liquid chromatography system (Thermo, Waltham, MA, USA) coupled to a Q-Exactive Plus via a nano-electrospray ion source (Thermo, Waltham, MA, USA). The peptide mixture was eluted from a 360-μm ID × 2 cm, C18 trap column and separated on a homemade 100 μm ID × 10 cm column with a linear 5–35% acetonitrile gradient at 500 nl/min. Survey scan were acquired after accumulation of 3e⁶ ions in Orbitrap for m/z 300–1400 using a resolution of 70,000 at m/z 400. The top 20 most intense precursor ions were selected for fragmentation in the HCD cell at normalized collision energy of 27%, and then fragment ions were transferred into the Orbitrap analyzer operating at a resolution of 17,500 at m/z 400. For the phosphopeptide identification and phosphosite quantification, raw spectral data were processed in Proteome Discoverer 1.4.1.14 suites with Mascot search engine against the rice genome annotation project database¹. The mass tolerance was set at 20 ppm for precursor, and 50 mmu for the tolerance of product ions. Oxidation (M), Acetyl (Protein-N term), and Phospho(S/T/Y) were set as variable modifications, and Carbamidomethyl (C) as static modification in the Mascot searches for phosphopeptides. Two missed cleavage on trypsin was allowed. The results were filtered for peptide with False discovery rates <1% by the Percolator tool of the Protein Discoverer package. PhosphoRS software within the Protein Discoverer software suite were used for an automatically re-analyzing of all the phosphopeptide hits, and those phosphorylation sites with a PhosphoRS probability higher than 90% were accepted as localized. Only those peptides which were phosphorylated in at least two of the three biological replicates were considered as truly phosphorylated. The phosphorylation is quantified based on the peak area under the ion intensity by using precursor ions area detector in PD1.4.1.14. Within-group means were calculated to determine fold changes. The differentially phosphorylated (DP) protein was defined to have over twofold changes in the average intensity with credible student’s t-test (P < 0.05).

CIAP treatment was performed by adding 1 μL of calf intestinal alkaline phosphatase (Takara, Dalian, China) into 20 μg protein of each sample for 30 min at 37°C. Then, Western blot was conducted following (Hou et al., 2015). The prepared protein samples were resolved in 10% SDS-polyacrylamide gel, and transferred onto a 0.45 μm polyvinylidene fluoride fluoropolymer (PVDF) membrane (Millipore, Darmstadt, Germany) by using an electrophoretic blotting system (Bio-Rad, Hercules, CA, USA). The immune-blot were detected with corresponding primary antibodies (1:1000 dilution), secondary antibody IgG conjugated with HRP (1:20,000 dilution), and visualized by using the enhanced chemiluminescence (Pierce, Waltham, MA, USA). β-tubulin protein was used as the internal control. The band intensities were quantified using the ImageJ software according to the instructions². All the sample intensities were first normalized to the control β-tubulin, and then calculated based on the ratio to set the relative level of 0 h into 1. The primary antibody against OsBZR1, D1, GID2, SAPK9, and SMG1 were commercially synthesized by immunizing rabbits and affinity purified by GenScript Company (Nanjing, China). The antigenic determinant peptide sequences are D1/RGA1: CSRSHSLSEAETTK; SMG1: MRPGGPPSLRAGLQC; SAPK9: MERAAAGPLGMEMPC; GID2: MSQPAELSREENVYC, and BZR1: CRPPKIRKPDWDVDP. Anti-β-tubulin (Cat No. M20005) were purchased from Abmart Company (Shanghai, China).

Phytohormone signaling is majorly transferred via protein phosphorylation cascades (Walton et al., 2015). In an effort to explore the roles of protein phosphorylation in BR signaling, we profiled the phosphoproteome of Nipponbare (Oryza sativa L. ssp japonica) 14 DAG (Day After Germination) seedlings by using a quantitative, label-free phosphoproteomic approach. Given that obvious genomic response to BR started at 3 h after treatment and robust physiological responses were observed after 12 h (Goda et al., 2002; Deng et al., 2007), three time points 0 hour (0 h), 3 hour (3 h), and 12 hour (12 h) after BR treatment, which represent the control, early BR response and late BR response, respectively, were selected for the proteome assay. Prior to the proteomic assay, we checked the mRNA expression level of three BR responsive genes, including two BR-induced genes ILI1 (Increased Leaf Inclination 1) (LOC_Os04g54900) and BU1 (Brassinosteroid Upregulated 1) (LOC_Os06g12210) and one BR-repressed gene IBH1 (ILI1 binding bHLH)(LOC_Os04g56500), at five time points (Tanaka et al., 2009; Zhang et al., 2009). As indicated in the qRT-PCR result, the expression levels of the three genes were vigorously altered, confirming a valid BR treatment on the plants in this experiment (Supplementary Figures S1A–C). The LC-MS/MS assay identified 3,412 phosphosites on 3,179 phosphopeptides at 0 h, 2,980 phosphosites on 2,780 phosphopeptides at 3 h and 2,507 phosphosites on 2,347 phosphopeptides at 12 h (Figure 1A). Of the phosphosites detected at 0 h, 89.7% were phosphoserine, 9.9% were phosphothreonine, and 0.4% were phosphotyrosines (Figure 1B), and the proportions were similar at the other two sampling times. Emerging evidences have shown that the distribution ratios of phosphorylation amino acids are much conserved in plants, despite the great variety of species, tissue and treatment applied (Han et al., 2014; Lv et al., 2014; Wang K. et al., 2014; Zhang M. et al., 2014; Hou et al., 2015; Qiu et al., 2016). In the three samples, 92.9–93.3% of the peptides carried only one phosphorylation group, 6.5–6.8% of the peptides carried two phosphorylation modifications, whereas only less than 0.3% of peptides had more than two phosphorylation modifications (Figure 1C). The 0, 3, and 12 h phosphopeptides were corresponded to 1668, 1500, and 1354 phosphoproteins, respectively, which represented the largest scale identification of BR responsive phosphoproteins in plants thus far.

As revealed by previous mega data analysis, featured motifs around the phosphosites provided valuable information of the kinase-substrate specificities. In this study, Motif-X³ (Chou and Schwartz, 2011) was employed to search the enriched motifs around phosphosites identified. The full set of 3,659 distinct phosphosites (3,283 phosphoserines, 361 phosphothreonines, and 15 phosphotyrosines) revealed at all three sampling times was submitted for analysis by motif-X software. As shown in Figure 2, at least five phosphoserine motifs (n > 200) and one phosphothreonine motif (n > 200) were identified (Figures 2A–E). As reported by many studies, [sP], [sxS], and [Rxxs] were frequently recurring (respectively, 1,498, 955, and 830 times) motifs, while the motif [Kxxs] (a basic S type) was recorded 266 times. [sF] was annotated to be a low frequency motif in plants by a previous meta-analysis (van Wijk et al., 2014). However, it seemed to be not the case, at least in rice, because two recent studies both found that [sF] was over-represented in rice leaves and developing seeds (Hou et al., 2015; Qiu et al., 2016). In the current study, we detected 252 hits of [sF] in rice young seedlings, which again supported that [sF] is a highly conserved motif for phosphorylation. [tP] was the only motif associated with phosphothreonine here (Figure 2F), while we did not found any conserved motifs for phosphotyrosine possibly due to the limited number of input sequences.

NetsurfP⁴ has been routinely employed to predict the surface accessibility and potential conserved secondary structures based on primary amino acid sequences (Petersen et al., 2009). By using the Netsurfp online tool, we found that 3093 of the 3283 (94.2%) phosphoserines and 254 of the 361 (70.4%) phosphothreonines lay on an exposed surface of the proteins, which is in agreement with the assumption that a surface location facilitates the access for the phosphorylation catalyzing enzyme(s) (Figure 2G). Nevertheless, all of the 15 phosphotyrosines were likely buried within the protein, displaying a divergent pattern from the other two phosphosite types. For the protein secondary structures, 327 (10.0%), 78 (2.4%), and 2797 (77.6%) of the phosphoserines were located in α-helix, β-strand and coil structures, respectively (Figure 2H). A similar distribution obtained for the phosphothreonines, implying a preference for coils with respect to the phosphorylation of both serine and threonine. The pattern for the phosphotyrosines was rather different: 73.3 and 26.7% of the locations were present in, respectively, an α-helix and a β-strand, with none present in a coil structure. The divergent results in surface accessibility and secondary structures implied different mechanisms of phosphorylation between serine/threonine and tyrosine.

Phosphorylation is a reversible, highly dynamic post-translational modification. Therefore, an altered phosphorylation pattern usually indicates the potential function of phosphorylation in the corresponding biological process. In this study, we collectively identified 3434 phosphopeptides which showed a DP pattern among the three time points of BR treatment (|log2(fold-change)|≥ 1, P < 0.05) (Table 1 and Supplementary Table S1). Of these, 598 were phosphorylated at 0 h, 214 at 3 h and 219 at 12 h, while the equivalent numbers of non-phosphorylated peptides were 136, 151, and 589 (Figure 3A); the remaining 1,527 peptides were phosphorylated at the three time points, but to a significantly varied degree (P < 0.05). A search of the set of known rice proteins allowed the DP peptides to be mapped onto 1,821 proteins, among which, respectively, 260, 92, and 100 were phosphorylated at the three time points, and respectively, 70, 81, and 252 were non-phosphorylated (Figure 3B). Of potential interest was that about 10% (191/1,821) of the proteins were either kinases or phosphatases, enzymes which are involved in protein phosphorylation/dephosphorylation, suggesting that these proteins are potential components of the BR signaling cascade. And indeed, several known BR signaling-related kinases, including OsGSK1 and OsSERK1 were detected to be DP (Koh et al., 2007; Li et al., 2009; Park et al., 2011). Rice genome contains at least 555 epi-genetic controlling factors (Gendler et al., 2008). Hou et al. (2015) reported that HDT701 and 27 other epigenetic controlling factors were DP in response to the Xoo infection, from which the proposal was that phosphorylation switch overriding the epi-genetic regulation may be a very universal model in the plant disease resistance pathway. In consistence with the report, the BR treatment also significantly altered the phosphorylation of 54 epigenetic controlling factors (involved in DNA methylation, histone methylation, histone acetylation and chromatin remodeling) and 118 transcription factors belonging to the families such as bHLHs, bZIPs, C3Hs, and Mybs. An analysis based on the CELLO algorithm showed that 52, 15, 14, and 10% of the DP proteins were located in the nucleus, chloroplast, cytoplasm, and plasma membrane, respectively. In contrast, proteins in the remaining 5 compartments such as mitochondrial, ER, and golgi accounted only less than 10% in total (Figure 3C).

Based on the phosphorylation intensity, a hierarchical clustering analysis divided the DP proteins into seven groups (Figure 3D). The various phosphorylation dynamic tendencies suggested different roles of the DP proteins in each group. For example, the members of group V were most strongly phosphorylated at 3 h, and much less so at 0 and 12 h, indicating their involvement in the early response to BR. Those in group VII, in contrast, may be involved in the later response, because their level of phosphorylation remained low throughout the sampling period. To our surprise, the majority of DP proteins belonged to group I, in which the level of phosphorylation declined over time.

By using the String 10.0⁵ (Search Tool for the Retrieval of Interacting Genes/Proteins) (Szklarczyk et al., 2015) and Cytoscape visualization (Shannon et al., 2003), we constructed PPI networks for DP kinases/phosphatases and all the BR induced DP proteins, respectively. Figure 4A depicted the sub-networks of DP kinases/phosphatases with high confidence score of 0.9. Eight PP2A phosphatases formed 11 edges (interaction relationships) sub-network, of which PP2A-2 (LOC_Os03g59060) and PP2Ac-4 (LOC_Os10g27050) were in the center. In Arabidopsis BR signaling, PP2As are critical phosphatases to dephosphorylate BZR1 and BZR2, which helps to retain the functional BZR1 and BZR2 in the nucleus and confer plants with BR response. The DP pattern suggested critical roles of these PP2As in BR response, whereas the predicted interaction relationships implied that PP2As may work in forms of protein complexes in BR signaling. A CDKG-2-CDKF-1-R2-SNT7 pathway was revealed in our PPI analysis. CDKG-2 (LOC_Os04g41100), CDKF-1 (LOC_Os06g22820), and R2 (LOC_Os05g32600) are cycling-dependent kinases functioning in cell proliferations, indicating this to be a key pathway in BR-regulated cell proliferation. In addition, we also detected a putative OsMEK1-MA3K.14-SMG1 pathway which is related to rice cell proliferation as well as stress resistance (Figure 4A). The full set of 1,821 DP proteins produced an interactome map composed of 459 nodes (proteins) and 946 edges. The various PPI sub-networks are associated with cellulose biosynthesis (Figure 4B), oxidation-reduction reactions (Figure 4C), phytohormone signaling (Figure 4D), and vascular ATP synthesis/transportation (Figure 4E).

To validate the MS identified phosphorylation status, we performed Western-blot analysis for seven selected DP protein in the time course of BR treatment, including BZR1 (LOC_Os07g39220), D1 (LOC_Os05g26890), GID2 (LOC_Os02g36974), SAPK9 (LOC_Os12g39630), and SMG1 (LOC_Os02g54600). Due to the phosphates attachment, phosphorylated protein bands migrate more slowly through the gel than the unmodified proteins. By using β-tubulin as an internal control, the intensity of each phosphorylated band was semi-quantified. Fit well with our expectations, all the five DP proteins showed similar phosphorylation tendencies as the phosphoproteomic data indicated, though the extent of phosphorylation may vary from each other (Figure 5). The reduced phosphorylation status of BZR1, D1, GID2, and SMG1 in response to the BR treatment was confirmed by the Western blot result, as was the induced increase in the level of phosphorylation for SAPK9. Moreover, the absence of the slower-migrating bands following treatment with CIAP (Calf Intestine Alkaline Phosphatase) confirmed their identity as phosphorylated proteins. The results above strongly suggested that our phosphoproteomic data is highly reliable and such a MS-based quantification strategy could be applied for the phosphorylation dynamic detections in other biological processes.

According to the constructed BR signaling pathway in A. thaliana, the initial sensing of BR involves the cell membrane receptor BRI1, from which the signal is passed through the phosphorylation-mediated cascade BSK1-BSU1-BIN2/GSK1, finally reaching the transcription factor BZR1 which triggers the BR response (Wang et al., 2012; Wang W. et al., 2014). The present phosphoproteomic analysis revealed that both of the BR signaling core components OsGSK1 (Glycogen Synthesase Kinase3-like 1) (LOC_Os01g10840) and OsBZR1 (LOC_Os07g39220) were DP. OsGSK1 belongs to the plant GSK3/SHAGGY-like protein kinase family, and is a close ortholog of AtBIN2. Meanwhile, OsBZR1, encoding a transcription factor with a DUF822 domain, shared extensive sequence homology with AtBZR1. The ectopical expression of OsGSK1 led to a stunted plant growth, which mimicked the typical phenotype of the BR-deficient mutants. Therefore, OsGSK1 might serve as a negative regulator of BR signaling, just as its ortholog BIN2 does in A. thaliana (Koh et al., 2007). The function of BZR1 also appears to be conserved between A. thaliana and rice. OsBZR1 RNAi suppressing lines showed dwarfism, erect leaves and reduced BR sensitivity, which is almost identical to the bzr1-D in A. thaliana (Wang et al., 2002; He et al., 2005).

A very interesting cytoplasm-nucleus shuttling model has been proposed for the BR signaling in the step from BIN2 to BZR1. In absence or low level of BR, upstream components of the cascade phosphorylate BIN2 to activate its kinase activity, which in turn catalyzes the phosphorylation on BZR1. The phosphorylation on BZR1 promotes its binding with a 14-3-3 protein, which allows the transport of BZR1 from the nucleus to the cytoplasm, where OsBZR1 remains non-functional. However, when the supply of BR is adequate, BIN2 is dephosphorylated, thereby losing its kinase activity; meanwhile, a PP2A de-phosphorylates BZR1, releasing it from its complex with the 14-3-3 protein, and allowing it to remain in the nucleus (Gampala et al., 2007). Bai et al. (2007) validated the 14-3-3 protein-mediated cytoplasm-nucleus shuttling mechanism of OsBZR1 in rice. Despite that they emphasized the binding of 14-3-3 to OsBZR1 retained its subcellular localization in cytoplasm under low BR level, the mechanistic basis of this binding has not be elucidated. The present experiment has shown that the extent of the phosphorylation affecting OsGSK1 and OsBZR1 was gradually reduced as a response to the BR treatment, in accordance with the behavior of BIN2 and BZR1 in A. thaliana. Perhaps more importantly, the differential phosphorylation of OsBZR1 implies that phosphorylation/dephosphorylation is likely to represent the means to determine whether it binds with 14-3-3 (as is the case in A. thaliana), which would mean that the BR signal transduction pathway is highly conserved between mono- and dicotyledonous plants. It is noteworthy that we failed to detect any phosphorylation on most of the reported BR signaling core components such as BRI1, BAK1, BSK1 and BSU1. In fact, a phosphoproteomic work of the BR signaling in A. thaliana encountered the same dilemma (Lin et al., 2015). The neglect of these proteins in the phosphoproteomic identification may be an artifact of low protein abundance and/or inadequate detection sensitivity. To overcome this, phosphorpoteomic works on more time points during the BR induction or different rice tissues might be necessary for wider proteomic coverage.

Plant growth and development is fine-tuned by profound cross-talk between the various phytohormones. The interaction between BR and other hormones has been extensively studied (Nemhauser et al., 2006; Vanstraelen and Benkova, 2012). Though GA and BR both are growth-promoting hormones, research in A. thaliana revealed an antagonistic relationship between them (Bouquin et al., 2001). In the rice root, exogenously supplied BR simultaneously represses certain GA synthesis genes, while promoting GA homeostasis genes, thereby finally reducing the level of bioactive GA present (De Vleesschauwer et al., 2012). In addition, the interaction of the GA and BR signaling pathway may also contribute to the complexity of their cross-talk. Some BR-insensitive mutants were found to be compromised in sensing GA signal, while some GA-hypersensitive mutants even showed enhanced sensitivities to BR (Gallego-Bartolome et al., 2012). Here, the two GA signaling proteins GID2 (LOC_Os02g36974) and D1/RGA1 (LOC_Os05g26890) were both down-phosphorylated in response to the BR treatment. gid2 and d1/rga1 exhibited GA-insensitivity and typical GA-defective phenotype such as dwarf, erected panicles, indicating that both of them are positive regulators of GA signaling (Ueguchi-Tanaka et al., 2000; Gomi et al., 2004). Through a SCF^GID2-proteasome pathway, GID2 mediates the degradation of GA signaling repressor SLR1 to activate the GA signal transduction (Gomi et al., 2004). D1/RGA1 is involved in GA and BR signaling, and it participates in the GA signaling majorly via a G protein-dependent pathway (Wang et al., 2006). SLR and GID2 act epistatic over D1 (Ueguchi-Tanaka et al., 2000). As we mentioned above, previous studies focused to explain the BR-GA antagonism from the view of quantity changes of GA synthesis and/or signaling genes in the transcriptional/translational levels. Nevertheless, the BR-induced down-phosphorylation of GID2 and D1/RGA1 is suggestive of a novel mechanism whereby BR inhibits the GA response by shutting down GA signaling. In addition to the GA signaling components, we also found that OSRK1/SAPK6 and SAPK4 involving in the ABA signaling (Chae et al., 2007; Diedhiou et al., 2008), CPT1 participating in auxin response (Haga et al., 2005) and OsEIN2 related to ethylene signaling (Jun et al., 2004) were DP in this study, suggesting a critical role of protein phosphorylation in the interplay between BR and other phytohormones.

One of the most direct effects of BR on plant is the altering of plant architecture, including plant height, leaf angle and tiller numbers (Zhang C. et al., 2014). Though BR defective mutants usually display unfavorable agronomic traits like dwarf, decreased leaf angle and more tillers, proper manipulation of the BR related genes was regarded as an effective way to improve crop architecture and eventually increase the yield (Sakamoto et al., 2006). Therefore, sorting out the protein connecting BR and plant architecture establishment will help us to construct the regulatory network and provide more gene resources for genetic improvement. As we expected, several rice architectures controlling proteins were found in the DP protein list. Cellulose synthesis for cell wall is a major determinant for the building of plant architecture and mechanic strength. Cellulose is synthesized at the plasma membrane by a complex containing multiple CESAs (cellulose synthase catalytic subunits) in plants. In rice, at least 10 CESAs (OsCESA1-10) have been identified, the contribution of three of which (OsCesA4, A7 and A9) has been defined by analysis of loss-of-function mutants, which all express a brittle leaf phenotype (Tanaka et al., 2003). Our PPI analysis of the DP proteins revealed a cellulose synthesis network containing five CESAs, which included OsCESA1, 2, 7, 8, and 9. The suggestion is therefore that BR influences cellulose synthesis and thereby has a profound effect on the architecture of the plant. Consistent with this notion, it was observed that several other cellulose synthesis-related proteins were DP. For example, BC3 (Brittle Culm3), a dynamin-related protein, is responsible for the membrane trafficking between the plasma membrane and intracellular compartments, which is an important process that regulates the deposition and metabolism of cellulose of the second cell wall. Mutation of BC3 led to inferior mechanical properties in rice plants (Xiong et al., 2010). The BR-induced DP pattern of multiple CESAs and other cellulose-related proteins in the current study strongly implied that BR regulates rice architectures by controlling cellulose synthesis via a phosphorylation-dependent mechanism.

Reactive oxygen species (ROS) production, including superoxide radical, hydroxyl radical and hydrogen peroxide (H₂O₂) is implicated as important regulatory and signaling elements for plants to adapt to unfavorable environments. It has been reported that application of exogenous BRs modified a wide range of antioxidant enzymes to counter the oxidative stress led by various environmental stimulus (Bajguz and Hayat, 2009). For example, application of BR promoted the activities of superoxide dismutase (SOD), catalase (CAT), ascorbate peroxidase (APX) in maize seedlings under water stress as well as in rice exposed to saline stress (Li et al., 1998; Núñez et al., 2003). Nevertheless, it is not clear yet that how BRs promote the antioxidant enzyme activities in response to stresses. Interestingly, we found that several DP antioxidant enzyme proteins such as CATA (LOC_Os02g02400), SOD (LOC_Os03g22810), and PEX5 (LOC_Os08g39080), which hinted that BR-induced protein phosphorylation may be involved in the adjustment of the activities of these antioxidant enzymes. Peroxisomes are single membrane-bound organelles where generation or degradation of ROS occurs. Peroxisomal proteins are originally translated in the cytoplasm and transported into peroxisomes. PEX5 is a PTS (Peroxisomal targeting signal) receptor protein, which recognizes the PTS and is responsible for the translocation of the peroxisomal proteins. In A. thaniana, Atpex5 mutant loses germinability in the absence of sucrose, but could be rescued by PEX5. PEX5 was also able to translocate PTS-containing proteins into the peroxisome by interacting with OsPEX7p (Lee et al., 2006). SOD is a type of key enzymes in antioxidant defense, which catalyze the partitioning of the superoxide radical into oxygen or less harmful H₂O₂, thus to alleviate the cell damages caused by oxidative stresses. Catalases act as scavengers of H₂O₂, and could further decompose it into oxygen and water. In rice, it has been known that CATB functions in an ABA-dependent manner to prevent the excessive accumulation of H₂O₂ under water stress (Ye et al., 2011). CATs are also targets of excess copper toxicity, which leads to retarded seed germination in rice (Ye et al., 2014). A very recent study suggested that the physical association-dissociation of GLO (Glycolate Oxidase; H₂O₂ producer) and CAT serve as a specific machinery to modulate H₂O₂ levels in rice (Zhang et al., 2016). Moreover, a potential PEX5-CATA-SOD interaction relationship was suggested by our PPI analysis (Figure 4C). Hence, we proposed that exogenous BR alters the activity of antioxidant enzymes via changing the phosphorylation status and/or intensity on antioxidant enzymes, and finally confers plants with resistance to environmental stresses.