A Comprehensive Proteomic Survey of ABA-Induced Protein Phosphorylation in Rice (Oryza sativa L.)

abscisic acid (ABA) is a key phytohormone regulating plant development and stress response. The signal transduction of ABA largely relies on protein phosphorylation. However; little is known about the phosphorylation events occurring during ABA signaling in rice thus far. By employing a label-free; MS (Mass Spectrometry)-based phosphoproteomic approach; we identified 2271 phosphosites of young rice seedlings and their intensity dynamics in response to ABA; during which 1060 proteins were found to be differentially phosphorylated. Western-blot analysis verified the differential phosphorylation pattern of D1, SMG1 and SAPK9 as indicated by the MS result; suggesting the high reliability of our phosphoproteomic data. The DP (differentially phosphorylated) proteins are extensively involved in ABA as well as other hormone signaling pathways. It is suggested that ABA antagonistically regulates brassinosteroid (BR) signaling via inhibiting BR receptor activity. The result of this study not only expanded our knowledge of rice phosphoproteome, but also shed more light on the pattern of protein phosphorylation in ABA signaling.


Introduction
Abscisic acid (ABA) is an isoprenoid hormone playing critical roles in various plant biological processes. Limited cellular water availability is a major abiotic factor leading to ABA accumulation and signal transduction [1]. Under dehydration or other types of osmotic stress, ABA promotes the stomata closure in guard cells to maintain water homeostasis in plants. The desiccation during seed maturation also induces high levels of ABA accumulation, which leads to seed dormancy and germination inhibition [2]. ABA signaling is a broad topic consisting of ABA metabolism and transport, ABA perception and signal transduction, as well as the ABA signal response and modulation [1]. Based on knowledge from Arabidopsis, a new "PYR/PYL/RCAR-PP2C-SnRK2" cascade model for ABA signaling has been proposed and validated [3,4]. The three core components of this model play vital but different roles in ABA signaling. PYR (PYrabactin Resistance 1)/PYL (PYR1-like)/RCAR (Regulatory Components of ABA Receptor) are water soluble ligand-binding proteins belonging to START (star-related lipid-transfer)-domain superfamily, and have been shown to be ABA receptors by two independent research groups [5,6]. PYR/PYL/RCAR harbors binding sites for PP2Cs (2C-type phosphorylated after treatment with five different phytohormones including ABA [27]. In an attempt to identify the SnRK2 protein kinase substrates in Arabidopsis, two independent labs conducted phosphoproteomic experiments in snrk2.2/2.3/2.6 triple mutants as well as the wild type [28,29]. Umezawa et al. [28] (2013) identified 35 peptides, including AtMPK1, AtMPK2, AREB1 and SNS1, which were specifically phosphorylated in Wild-type but not in the triple mutant plants treated with ABA. Similarly, Wang et al. [29] (2013) found that 58 proteins were not sensitive or less sensitive to ABA treatment in the triple mutants. Further in vitro kinase assays confirmed that most of the 58 proteins were SnRK2 substrates involved in various biological processes. In addition to the SnRK2 deficient plants, a paper was published very recently reporting the rapid phosphoproteomic effects of ABA on wild type and ABA receptor quadruple mutants in Arabidopsis [30].
Rice is a major crop feeding over half of the global population. It has been fully recognized that ABA plays critical roles in rice seed germination, grain filling, biotic and abiotic stress responses [31,32]. Therefore, clarifying the molecular mechanism of ABA signaling in rice would be of great biological as well as economic significance for rice genetic improvement. Unfortunately, the protein phosphorylation events, in particular the phosphosites and intensity dynamics, in ABA signaling have been barely studied. In this study, we focused on the profiling of young rice seedling phosphosites, phosphopeptides and phosphorylation intensity dynamics in response to ABA by employing a label-free, MS-based phosphoproteomic approach, and aimed to reveal the biological significance of protein phosphorylation in ABA signaling from the phosphoproteomic data.

Profile of the Abscisic Acid (ABA)-Induced Phosphoproteome in Rice
Profiling the phosphorylation proteins, peptides and sites is fundamental and crucial to understand the roles of protein phosphorylation in ABA signaling in rice. In the current study, a MS-based, quantitative, label-free phosphoproteomic approach was employed to identify the phosphopeptides and phosphosites from 14-day-old young seedlings of Nipponbare (Oryza sativa, ssp. japonica) upon ABA treatment. Because exogenous ABA could be efficiently up taken by intact rice roots [33], the ABA treatment was conducted by immersing young seedlings in water containing 100 µM ABA. Previous study indicated that robust rice proteomic response to exogenous ABA started at 3 h and lasted over 6 h [34], we therefore chose 3 and 12 h to represent the ABA early and late signaling stages respectively, while seedlings without ABA treatment was used as a control (CK). Samples in each time point were conducted with three biological replicates. Prior to the protein extraction for MS analysis, we tested the transcriptional level of Rab16a (LOC_Os11g26790), a marker gene of ABA induction in rice, by quantitative RT-PCR (Real Time Polymerase Chain Reaction) [35]. As shown in Figure 1A, the level of Rab16a was significantly increased in response to ABA, indicating a valid ABA treatment in the experiment. In this study, a total of 1774, 1344 and 1653 phosphosites from 1704, 1279 and 1577 phosphopeptdies were identified in CK, 3 and 12 h, respectively ( Figure 1B). After putative redundancy was filtered and removed, collectively 2271 phosphosites on 2162 phosphopeptides were profiled on 1086 phosphoproteins. Over 90% of the phosphopeptides in all the three samples carried only one phosphorylation modification, close to 10% of the phosphopeptides were phosphorylated twice, whereas triple phosphorylation was very rarely detected ( Figure 1C). Additionally, we counted ratios phosphorylation which occurred on serine, threonine and tyrosine respectively. Phosphoserine was the predominant phosphorylations type and accounted for over 94% of the phosphosites, around 5% of the phosphosites were on threonine, and only 0.2%-0.3% of the phosphosites were on tyrosine ( Figure 1D).

Differentially Phosphorylated Proteins in Response to ABA Treatment
A differential phosphorylation pattern usually indicates the important regulatory roles of the phosphoprotein in the corresponding biological process. Upon ABA treatment, a total of 1549 peptides were found to be differentially phosphorylated (Fold change > 2, p < 0.05). There were 366, 143 and 185 peptides which were specifically phosphorylated in CK, 3 and 12 h respectively. Meanwhile, 130, 332 and 76 peptides were found to be specifically not phosphorylated in CK, 3 and 12 h (Figure 2A). At the protein level, we identified 1060 DP proteins, including 268, 92, 123 proteins which were specifically phosphorylated and 91, 230, 55 proteins which were specifically not phosphorylated, in CK, 3 and 12 h respectively ( Figure 2B). Among the DP proteins, we found 96 kinases and 19 phosphatases whose function are specifically for protein phosphorylation or dephosphorylation. In addition, the DP proteins also covered 65 transcription factors, including 9 bZIP family members (Table S1). It is known that bZIP is a major group of ABA responsive factors that could be phosphorylated by SnRK2 kinases in ABA signaling. The ABA-induced phosphorylation pattern strongly implied that these bZIPs are functionally involved in ABA signaling. By blasting against the potential SnRK2 target proteins obtained from Arabidopsis [28,29], we found that 154 of the DP proteins are close orthologs of these SnRK2 targets (E value < 1 × 10 −10 , score ≥ 80 and identity ≥ 80), suggesting their involvement in the ABA core signaling pathway (Table S2). Based on the phosphorylation intensity pattern at the three time points, a hierarchical clustering analysis divided the DP proteins into seven clades, indicating their divergent roles in ABA signaling ( Figure 2C; Table S3). For example, proteins in clade V had a low level of phosphorylation at CK, while the phosphorylation intensity increased at 3 h and reached the maximum number at 12 h, indicating their roles in ABA response in the late stage. Meanwhile, proteins in clade VII are proposed to be involved in the early ABA signaling because their phosphorylation intensity significantly elevated at 3 h, but soon dropped back to a low level at 12 h.

Differentially Phosphorylated Proteins in Response to ABA Treatment
A differential phosphorylation pattern usually indicates the important regulatory roles of the phosphoprotein in the corresponding biological process. Upon ABA treatment, a total of 1549 peptides were found to be differentially phosphorylated (Fold change > 2, p < 0.05). There were 366, 143 and 185 peptides which were specifically phosphorylated in CK, 3 and 12 h respectively. Meanwhile, 130, 332 and 76 peptides were found to be specifically not phosphorylated in CK, 3 and 12 h (Figure 2A). At the protein level, we identified 1060 DP proteins, including 268, 92, 123 proteins which were specifically phosphorylated and 91, 230, 55 proteins which were specifically not phosphorylated, in CK, 3 and 12 h respectively ( Figure 2B). Among the DP proteins, we found 96 kinases and 19 phosphatases whose function are specifically for protein phosphorylation or dephosphorylation. In addition, the DP proteins also covered 65 transcription factors, including 9 bZIP family members (Table S1). It is known that bZIP is a major group of ABA responsive factors that could be phosphorylated by SnRK2 kinases in ABA signaling. The ABA-induced phosphorylation pattern strongly implied that these bZIPs are functionally involved in ABA signaling. By blasting against the potential SnRK2 target proteins obtained from Arabidopsis [28,29], we found that 154 of the DP proteins are close orthologs of these SnRK2 targets (E value < 1 × 10 −10 , score ≥ 80 and identity ≥ 80), suggesting their involvement in the ABA core signaling pathway (Table S2). Based on the phosphorylation intensity pattern at the three time points, a hierarchical clustering analysis divided the DP proteins into seven clades, indicating their divergent roles in ABA signaling ( Figure 2C; Table S3). For example, proteins in clade V had a low level of phosphorylation at CK, while the phosphorylation intensity increased at 3 h and reached the maximum number at 12 h, indicating their roles in ABA response in the late stage. Meanwhile, proteins in clade VII are proposed to be involved in the early ABA signaling because their phosphorylation intensity significantly elevated at 3 h, but soon dropped back to a low level at 12 h.
In addition, the co-phosphorylation pattern of the proteins in each clade hinted that the members may work in the same or similar pathways of ABA signaling.
In addition, the co-phosphorylation pattern of the proteins in each clade hinted that the members may work in the same or similar pathways of ABA signaling.

Motif-X and GO (Gene Ontology) Analysis of Differentially Phosphorylated Proteins
Substrate specificity of kinases is largely determined by the motifs harboring the phosphosite. Therefore, a view of the conserved motifs around phosphosites of the DP proteins may provide useful clues to find the kinases catalyzing phosphorylation on them. With the aid of a written program, we extracted a 13 amino acid-long sequence with the phosphosite in the center for all the 1636 unique phosphosites on DP proteins, including 1484 phosphoserines, 148 phosphothreonines and 4 phosphotyrosines. Subsequently, 1589 AA sequences were extracted and inputted into an online tool Motif-X (http://motif-x.med.harvard.edu/) [36] for conserved motif detection and resulted in the identification of 4 phosphoserine motifs (sequences number n > 100) and 1 phosphothreonine motif (n > 100) ( Figure 3A-E). As found in many other cases, proline-directed motifs (sP) and (tP) are most commonly present around the phosphosites of our DP proteins [37][38][39][40]. Proteins with a (sP) motif have been reported to be substrates of MAPKs, SnRK2s (sucrose non-fermenting1-related protein kinase 2), RLKs (receptor-like kinases) and many other kinases [41]. We also identified two basophilic motifs (Rxxs) and (LxRxxs) which could be recognized by MAPKK, CaMK (calmodulin-dependent protein kinase)-II and protein kinase A [41]. Interestingly, (sF) has been rarely reported by previous reports in many species including rice [37,38,41]; this study together with two of our previous cases we found it as an over-represented motif in rice leaf, seed and seedlings [39,40]. This observation might be ascribed to the difference in methods of protein extraction and MS identification.

Motif-X and GO (Gene Ontology) Analysis of Differentially Phosphorylated Proteins
Substrate specificity of kinases is largely determined by the motifs harboring the phosphosite. Therefore, a view of the conserved motifs around phosphosites of the DP proteins may provide useful clues to find the kinases catalyzing phosphorylation on them. With the aid of a written program, we extracted a 13 amino acid-long sequence with the phosphosite in the center for all the 1636 unique phosphosites on DP proteins, including 1484 phosphoserines, 148 phosphothreonines and 4 phosphotyrosines. Subsequently, 1589 AA sequences were extracted and inputted into an online tool Motif-X (http://motif-x.med.harvard.edu/) [36] for conserved motif detection and resulted in the identification of 4 phosphoserine motifs (sequences number n > 100) and 1 phosphothreonine motif (n > 100) ( Figure 3A-E). As found in many other cases, proline-directed motifs (sP) and (tP) are most commonly present around the phosphosites of our DP proteins [37][38][39][40]. Proteins with a (sP) motif have been reported to be substrates of MAPKs, SnRK2s (sucrose non-fermenting1-related protein kinase 2), RLKs (receptor-like kinases) and many other kinases [41]. We also identified two basophilic motifs (Rxxs) and (LxRxxs) which could be recognized by MAPKK, CaMK (calmodulin-dependent protein kinase)-II and protein kinase A [41]. Interestingly, (sF) has been rarely reported by previous reports in many species including rice [37,38,41]; this study together with two of our previous cases we found it as an over-represented motif in rice leaf, seed and seedlings [39,40]. This observation might be ascribed to the difference in methods of protein extraction and MS identification. By using the protein subcellular localization prediction tool WoLF PSORT (http://www. genscript.com/wolf-psort.html) [42], we analyzed protein distribution in various cellular compartments. The results showed that nuclear, chloroplast and cytoplasmic are the three most enriched cellular compartments which accounted for 45.4%, 22.5% and 16.9% of the DP protein destinations. It was also found that 9.2% and 3.3% of the DP proteins are localized in the plasma membrane and mitochondria respectively, while less than 3% of the DP proteins are distributed in the remaining six cellular compartments analyzed ( Figure 3F). It has been known that ABA signal activates a lot of ion channels and transporters, which are usually lipid-soluble proteins located in plasma membrane. The detection of relatively less number of plasma membrane proteins may attribute to the phenol method used for protein preparation, which majorly extracted the water-soluble proteins but excluded some of the lipid-soluble proteins. Therefore, phosphoproteomic analysis of proteins isolated from various cellular fractions may be necessary to obtain a more comprehensive result.
To assign the biological relevance of the DP proteins, GO (Gene Ontology) analysis of them was performed in terms of "cellular component", "molecular function" and "biological process" and visualized by WEGO (Web Gene Ontology Annotation Plot) [43]. In the perspective of "cellular component", categories of "cell" and "external encapsulating structure" were enriched, while "cell part", "intracellular" and "membrane" were under-represented when compared with the distribution ratio of the rice whole genome encoding proteins ( Figure 4). Interestingly, we found that the ABA DP proteins showed over-represented molecular functions in "catalytic". Given that proteins in "catalytic" are mostly enzymes like kinases and phosphatases, the enriched distribution suggested the extensive involvement of kinases and phosphatases in the ABA response. On the contrary, the "binding" was less favored for the ABA-induced DP proteins ( Figure 4). In terms of "biological process", DP proteins were enriched in "metabolic process" and "response to stimuli", suggesting that ABA predominantly affects the signaling and metabolism via protein phosphorylation. Meanwhile, we found that DP proteins were less preferentially distributed in "cellular process" and "cellular component organization". By using the protein subcellular localization prediction tool WoLF PSORT (http://www.genscript. com/wolf-psort.html) [42], we analyzed protein distribution in various cellular compartments. The results showed that nuclear, chloroplast and cytoplasmic are the three most enriched cellular compartments which accounted for 45.4%, 22.5% and 16.9% of the DP protein destinations. It was also found that 9.2% and 3.3% of the DP proteins are localized in the plasma membrane and mitochondria respectively, while less than 3% of the DP proteins are distributed in the remaining six cellular compartments analyzed ( Figure 3F). It has been known that ABA signal activates a lot of ion channels and transporters, which are usually lipid-soluble proteins located in plasma membrane. The detection of relatively less number of plasma membrane proteins may attribute to the phenol method used for protein preparation, which majorly extracted the water-soluble proteins but excluded some of the lipid-soluble proteins. Therefore, phosphoproteomic analysis of proteins isolated from various cellular fractions may be necessary to obtain a more comprehensive result.
To assign the biological relevance of the DP proteins, GO (Gene Ontology) analysis of them was performed in terms of "cellular component", "molecular function" and "biological process" and visualized by WEGO (Web Gene Ontology Annotation Plot) [43]. In the perspective of "cellular component", categories of "cell" and "external encapsulating structure" were enriched, while "cell part", "intracellular" and "membrane" were under-represented when compared with the distribution ratio of the rice whole genome encoding proteins (Figure 4). Interestingly, we found that the ABA DP proteins showed over-represented molecular functions in "catalytic". Given that proteins in "catalytic" are mostly enzymes like kinases and phosphatases, the enriched distribution suggested the extensive involvement of kinases and phosphatases in the ABA response. On the contrary, the "binding" was less favored for the ABA-induced DP proteins (Figure 4). In terms of "biological process", DP proteins were enriched in "metabolic process" and "response to stimuli", suggesting that ABA predominantly affects the signaling and metabolism via protein phosphorylation. Meanwhile, we found that DP proteins were less preferentially distributed in "cellular process" and "cellular component organization".

Validation of the Phosphorylation Patterns by Western-Blot
To validate the MS identified phosphorylation sites of the DP proteins, we manually checked the spectrum map of D1 (LOC_Os05g26890), SMG1 (LOC_Os02g54600) and SAPK9 (LOC_Os12g39630) ( Figure S1) and synthesized their antibodies for Western-blot analysis. Because phosphorylation could increase the molecular weight and lag the target protein band shift in electrophoresis, it is expected to have a bigger band on the top of the original protein band if protein phosphorylation has occurred. Samples were also treated with CIAP (Calf Intestine Alkaline Phosphatase) to confirm that the lagged band is the phosphorylated target protein, instead of a non-specific background band in the Western-blot. β-Tubulin was used as an internal control for the semi-quantification of signals in each sample. The western-blot results showed that the relative phosphorylation band signals of SMG1 at 3 and 12 h were decreased to 49% and 45% respectively when CK was used as the control, which is in consensus to the 43% and 66% in our phosphoproteomic data ( Figure 5, Table 1). The result of D1 was similar to SMG1. The MS identification revealed that SAPK9 was lowly phosphorylated at CK, but slowly increased phosphorylation at 3 and 12 h. Accordingly, we found that the phosphorylation band signals of SAPK9 at CK, 3 and 12 h were 1, 1.21 and 2.69 respectively ( Figure 5, Table 1). 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.

Validation of the Phosphorylation Patterns by Western-Blot
To validate the MS identified phosphorylation sites of the DP proteins, we manually checked the spectrum map of D1 (LOC_Os05g26890), SMG1 (LOC_Os02g54600) and SAPK9 (LOC_Os12g39630) ( Figure S1) and synthesized their antibodies for Western-blot analysis. Because phosphorylation could increase the molecular weight and lag the target protein band shift in electrophoresis, it is expected to have a bigger band on the top of the original protein band if protein phosphorylation has occurred. Samples were also treated with CIAP (Calf Intestine Alkaline Phosphatase) to confirm that the lagged band is the phosphorylated target protein, instead of a non-specific background band in the Western-blot. β-Tubulin was used as an internal control for the semi-quantification of signals in each sample. The western-blot results showed that the relative phosphorylation band signals of SMG1 at 3 and 12 h were decreased to 49% and 45% respectively when CK was used as the control, which is in consensus to the 43% and 66% in our phosphoproteomic data ( Figure 5, Table 1). The result of D1 was similar to SMG1. The MS identification revealed that SAPK9 was lowly phosphorylated at CK, but slowly increased phosphorylation at 3 and 12 h. Accordingly, we found that the phosphorylation band signals of SAPK9 at CK, 3 and 12 h were 1, 1.21 and 2.69 respectively ( Figure 5, Table 1). 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.

Discussion
In plants, it has been widely accepted that ABA signals are transmitted majorly through the "PYR/PYL/RCAR-PP2C-SnRK2" cascade, in which protein phosphorylation/dephosphorylation plays essential roles [3]. As the first step toward understanding ABA signaling, profiling the ABA-induced phosphoproteins, phosphosites and phosphorylation dynamics is fundamental. So far, most of our knowledge on this issue was obtained from fragmented kinase function characterization researches by using traditional biochemical assays. Some primary phosphoproteomic methods such as 2D-DIGE (two dimension difference gel electrophoresis) were applied for the investigation of phosphorylation patterns in rice upon ABA treatment, but yielded limitedphosphoproteins data [34,44]. Compared with 2D-DIGE, high resolution LC-MS (liquid chromatograph-mass spectrometer)/MS has been proven to be more sensitive, efficient and accurate in phosphoprotein identification. In the current study, we demonstrated a high through-put, accurate, phosphoproteomic identification of the ABA-induced protein phosphorylation events in rice by utilizing a label-free, MS-based approach. A total of 1086, 844 and 1019 phosphoproteins were profiled from young rice seedlings upon ABA treatment at three time points, representing the largest dataset of ABA-induced phosphoproteins in plants. In addition, our MS data provided highly reliable phosphorylation intensity quantification, which has been verified by the Western-blot analysis of SMG1, D1 and SAPK9. With the quantification data, 1060 phosphoproteins including kinases, phosphotases, transcription factors and so on, were found to be differentially phosphorylated in response to ABA application, suggesting that ABA signaling is a highly complex process with a wide range of proteins involved. Meanwhile, the DP pattern also revealed the potential roles of these proteins in ABA response.
According to the "PYR/PYL/RCAR-PP2C-SnRK2" model, protein phosphorylation starts from the SnRK2s' autophosphorylation, while no phosphorylation would occur on ABA receptors and PP2Cs (more specifically, the clade A PP2Cs) as they only competitively bind with ABA or each other to decide the phosphorylation status of SnRK2s [3]. The rice genome contains 12 ABA receptors, 10 clade A PP2Cs and 10 SnRK2 genes [23,45,46]. Consistent with previous reports, we did not find any ABA receptors or clade A PP2Cs in our DP protein list. On the contrary, the phosphorylation of 2 SnRK2 including SAPK6 and 9 were significantly induced by ABA treatment.

Discussion
In plants, it has been widely accepted that ABA signals are transmitted majorly through the "PYR/PYL/RCAR-PP2C-SnRK2" cascade, in which protein phosphorylation/dephosphorylation plays essential roles [3]. As the first step toward understanding ABA signaling, profiling the ABA-induced phosphoproteins, phosphosites and phosphorylation dynamics is fundamental. So far, most of our knowledge on this issue was obtained from fragmented kinase function characterization researches by using traditional biochemical assays. Some primary phosphoproteomic methods such as 2D-DIGE (two dimension difference gel electrophoresis) were applied for the investigation of phosphorylation patterns in rice upon ABA treatment, but yielded limitedphosphoproteins data [34,44]. Compared with 2D-DIGE, high resolution LC-MS (liquid chromatograph-mass spectrometer)/MS has been proven to be more sensitive, efficient and accurate in phosphoprotein identification. In the current study, we demonstrated a high through-put, accurate, phosphoproteomic identification of the ABA-induced protein phosphorylation events in rice by utilizing a label-free, MS-based approach. A total of 1086, 844 and 1019 phosphoproteins were profiled from young rice seedlings upon ABA treatment at three time points, representing the largest dataset of ABA-induced phosphoproteins in plants.
In addition, our MS data provided highly reliable phosphorylation intensity quantification, which has been verified by the Western-blot analysis of SMG1, D1 and SAPK9. With the quantification data, 1060 phosphoproteins including kinases, phosphotases, transcription factors and so on, were found to be differentially phosphorylated in response to ABA application, suggesting that ABA signaling is a highly complex process with a wide range of proteins involved. Meanwhile, the DP pattern also revealed the potential roles of these proteins in ABA response.
According to the "PYR/PYL/RCAR-PP2C-SnRK2" model, protein phosphorylation starts from the SnRK2s' autophosphorylation, while no phosphorylation would occur on ABA receptors and PP2Cs (more specifically, the clade A PP2Cs) as they only competitively bind with ABA or each other to decide the phosphorylation status of SnRK2s [3]. The rice genome contains 12 ABA receptors, 10 clade A PP2Cs and 10 SnRK2 genes [23,45,46]. Consistent with previous reports, we did not find any ABA receptors or clade A PP2Cs in our DP protein list. On the contrary, the phosphorylation of 2 SnRK2 including SAPK6 and 9 were significantly induced by ABA treatment. Kobayashi et al. [23] (2004) demonstrated that SAPK8, 9 and 10 protein accumulations is induced by ABA treatment. There is also literatures suggesting that ABA up-regulate the transcription levels of SAPK4 and 6, but the final effects in the protein level remain unclear [24,47]. Assuming that an induced expression pattern usually indicates the function in the process, it is rationale for us to propose that SAPK6, 8, 9 and 10 are functionally involved in the ABA signaling. Nevertheless, our study did not detect any obvious phosphorylation change in SAPK8 and 10. To explain this discrepancy, we checked the mRNA expression level of SAPK8 and 10 in different tissues, and found that the expression of both genes in young seedlings is in extremely low in abundance, which may result in the poor or omitted MS identification of them. Other possible reasons could be that SAPK8 and 10 respond to ABA in an earlier or later stage which is beyond the range checked in this study, or the growth conditions, such as high humidity (90%) during hydroponic culture down-regulated the ABA signaling pathway prior to the ABA treatment. Therefore, developing more sensitive MS identification equipments and methods, and also extending the time range of ABA treatment will be helpful for us to thoroughly profile the ABA induced phosphorylation events.
Phytohormone regulation highly complex process, in which the reception of one hormone usually triggers the synergistic or agonistic interplay with other hormones. It has been known that brassinosteroids (BRs) act antagonistically with ABA in promoting seed germination and post-germinative growth. The BR signaling repressor BIN2 (Brassinosteroid INsensitive 2) can stabilize and phosphorylate the positive ABA regulator ABI5 (Abscisic Acid Insensitive 5) in presence of ABA, thus positively regulate ABA response. The antagonism from BR to ABA is probably mediated by the inactivation of the BIN2-ABI5 cascade when BRs are applied to the plants [48]. However, a view of the ABA-BR antagonism from the ABA side has not yet been achieved. In the current study, we found that SERK2/OsBAK1 (LOC_Os04g38480), an ortholog of Arabidopsis BR co-receptor BAK1 (BRI1-Associated receptor Kinase 1) was down-phosphorylated by ABA treatment. In Arabidopsis, the BAK1-BRI1 complex binds with BR to sense the signal [49]. Suppressing of SERK2 in rice leads to severe BR insensitivity and other BR-deficient phenotypes such as altered plant architecture and immune response, indicating SERK2 plays key roles in BR signaling, possibly as a receptor [50,51]. The ABA-induced dephosphorylation of SERK2 provided strong evidence that ABA antagonistically regulates BR signaling via inhibiting the BR receptor activity. In addition to BR components, several proteins involved in the signaling of other phytohormones were differentially phosphorylated as well. For example, D1/RGA1 (LOC_Os05g26890) is a key regulator of the G protein-dependent GA signal transduction pathway. d1 plants show pleiotropic phenotypes including dwarf, erected panicles, small and roundish seed shape due to the reduced cell proliferation [52]. Our phosphoproteomic data revealed that the phosphorylation intensities of D1/RGA1 at 3 and 12 h were decreased to 30% and 48% respectively when compared with that of the CK, which is in agreement with the ABA-GA antagonistic effects. Moreover, MHZ7 (LOC_Os07g06130) involved in ethylene signaling, OsDOS (LOC_Os01g09620) and C3H12 (LOC_Os01g68860) participating in JA (jasmonic acid) signaling were all differentially phosphorylated during the ABA treatment, suggesting the extensive cross-talk between ABA and these phytohormones [53,54].
Downstream of the constructed "PYR/PYL/RCAR-PP2C-SnRK2" ABA signaling pathway, the ABA signal is finally transmitted by prtoein phosphorylation from SnRK2 to ABF transcription factors, of which many have been reported to be bZIPs [25,[55][56][57]. In the current study, we detected 10 differentially phosphorylated bZIPs, 6 of which showed an up-phosphorylation pattern upon ABA induction. bZIP72 is a member of the group A OsbZIP family with ABRE trans activity in yeast [55]. qRT-PCR analysis showed that bZIP72 was highly expressed in most of the rice tissues and can be induced by ABA, ethylene, MeJA, PEG, salinity and cold. More interestingly, bZIP72 over-expression lines became hypersensitive to ABA, and showed retarded seedling growth, delayed seed germination and enhanced drought tolerance when compared with the Wild-type. ABF1 (ABA responsive element binding factor 1) is another reported ABA signaling-related bZIP gene. The expression of ABF1 is induced by ABA and various abiotic stresses. In the allelic mutants Osabf1-1 and Osabf1-2, some ABA responsive genes were suppressed, and the plants became more sensitive to drought and salinity treatments than wild type plants [57]. Though it has been known that ABA induces the quantity of bZIP72 and ABF1 in the mRNA level, how these two bZIPs are involved in ABA signaling remains unclear. Their DP pattern as revealed in this study may suggest that bZIP72 and ABF1 participate in the ABA signaling via the phosphorylation activation from the upstream components in the pathway, possibly SnRK2s. More importantly, the revealed phosphosites of these transcription factors provided valuable information for studies to elucidate the mechanism of these proteins in ABA signaling. For example, bZIP72 harbors over 40 potential phosphosites as predicted by Rice Phospho 1.0 tool [58]. However, only the Ser71 and Ser116 of bZIP72 were detected and up-phosphrylated under ABA treatment in our phosphoproteomic data (Table 1). Thus, Ser71 and Ser116 are very likely to be the target sites of bZIP72 in ABA signaling, which will be further investigated in our lab in the future.

ABA Treatment
Fourteen-day-old Nipponbare young seedlings were obtained from pure water culture in a growth chambers (28 • C, 90% relative humidity, 14 h daytime and 10 h night time). For the ABA treatment, healthy, intact young seedlings were kept in a beaker with water containing 100 µM ABA (Sigma, St. Louis, MO, USA). The whole seedling samples were collected at CK, 3 and 12 h respectively, and immediately used for RNA isolation and protein extraction.

RNA Isolation and qRT-PCR (Quantitative Real Time-Polymerase Chain Reaction)
The total RNA of the whole young seedlings was extracted by Trizol (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions. DNaseI (Takara, Dalian, China) was added in the total RNA for eliminating DNA contamination. The DNA free total RNA was quantified by a Nanodrop spectrometer (Thermo, Waltham, MA, USA) and 2 µg RNA was used for each reverse transcription by using M-MLV reverse transcriptase (Takara, Dalian, China) with Oligo(dT) 20 (Invitrogen, Carlsbad, CA, USA) as primer. Quantitative real-time RT-PCR analysis was conducted in a total volume of 10 µL containing 2 µL reverse-transcribed product above, 0.2 µM each primer and 5 µL THUNDERBIRD SYBR ® qPCR Mix (Toyobo, Shanghai, China) on a CFX96 touch realtime PCR detection system (Bio-Rad, Hercules, CA, USA). The ubiquitin (GenBank accession No. AF184280) was used as an internal control. The relative expression level of tested genes was calculated based on C t values by following the 2 −∆∆Ct method [59]. Primers sequences are as follows: Ubif: 5 -GCTCCGTGGCGGTATCAT-3 ; Ubir: 5 -CGGCAGTTGACAGCCCTAG-3 ; Rab16af: 5 -GAGGGAGGAGCACAAGACC-3 and Rab16ar: 5 -ATTCCATCATCCTCAGAC GAG-3 .

LC-MS/MS (Liquid Chromatograph-Mass Spectrometer/Mass Spectrometer) and Data Analysis
This was done by following the methods of Qiu et al. [40] (2016) in Beijing Proteome Research Center. LC-MS/MS analyses were performed on an Easy-nLC 1000 liquid chromatography system (Thermo, Waltham, WA, USA) coupled to a Q-Exactive Plus via a nano-electrospray ion source (Thermo). 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 (C18, 1.9 µm, 120 Å, Maisch GmbH) with a linear 5%-35% acetonitrile gradient at 500 nL/min. Survey scan were acquired after accumulation of 8.06 × 10 2 ions in Orbitrap for m/z 300-1400 using a resolution of 70,000 at m/z 400. The top 20 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. The dynamic exclusion of previously acquired precursor ions was enabled at 18 s. Raw spectral data were processed for phosphopeptide identification and phosphosite quantification in Proteome Discoverer 1.4.1.14 suites with Mascot search engine against the rice genome annotation project database: (ftp://ftp.plantbiology.msu.edu/pub/data/Eukaryotic_ Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/). 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 at high identification confidence (False discovery rates < 1%) and E Value < 0.05 by the Percolator tool of the Protein Discoverer package. Low-scoring peptides (Mascot score ≤ 20) were excluded from the analysis when they were not further supported by additional high-scoring identifications in other replicates or experiments. All phosphopeptide hits were automatically re-analyzed by the phosphoRS software within the Protein Discoverer software suite. PhosphoRS probability higher than 90% was required for a phosphorylation site to be considered 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. Extracted ion chromatogram (sum of areas under the curves of MS1 chromatograms of each peptide identified) was used for normalization across all runs. Within-group means were calculated to determine fold changes. The differentially phosphorylated protein was defined to have over two fold changes in the average intensity with credible student's t-test (p < 0.05). The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium [60] via the PRIDE partner repository with the dataset identifier PXD004939.

Conclusions
Here, we report the profiling of 2271 phosphosites on 2162 phosphopeptides from young rice seedlings and the quantification of their intensity dynamics in response to ABA. 1060 proteins were found to be differentially phosphorylated during the ABA treatment. The DP proteins are extensively involved in ABA as well as other hormone signaling pathways. It is suggested that ABA antagonistically regulates BR signaling via inhibiting the BR receptor activity. Western-blot analysis verified the differential phosphorylation pattern of D1, SMG1 and SAPK9 as indicated by the MS result. Our result of this study expands our knowledge of rice phosphoproteome, but also provides insight into the mechanism of protein phosphorylation in ABA signaling.