Global Profiling of Phosphorylation Reveals the Barley Roots Response to Phosphorus Starvation and Resupply

Phosphorus (P) deficiency is a major threat to the crop production, and for understanding the response mechanism of plant roots, P stress may facilitate the development of crops with increased tolerance. Phosphorylation plays a critical role in the regulation of proteins for plant responses to biotic and abiotic stress; however, its functions in P starvation/resupply are largely unknown for barley (Hordeum vulgare) growth. Here, we performed a global review of phosphorylation in barley roots treated by P starvation/resupply. We identified 7,710 phosphorylation sites on 3,373 proteins, of which 76 types of conserved motifs were extracted from 10,428 phosphorylated peptides. Most phosphorylated proteins were located in the nucleus (36%) and chloroplast (32%). Compared with the control, 186 and 131 phosphorylated proteins under P starvation condition and 156 and 111 phosphorylated proteins under P resupply condition showed significant differences at 6 and 48 h, respectively. These proteins mainly participated in carbohydrate metabolism, phytohormones, signal transduction, cell wall stress, and oxidases stress. Moreover, the pathways of the ribosome, RNA binding, protein transport, and metal binding were significantly enriched under P starvation, and only two pathways of ribosome and RNA binding were greatly enriched under Pi resupply according to the protein–protein interaction analysis. The results suggested that the phosphorylation proteins might play important roles in the metabolic processes of barley roots in response to Pi deficiency/resupply. The data not only provide unique access to phosphorylation reprogramming of plant roots under deficiency/resupply but also demonstrate the close cooperation between these phosphorylation proteins and key metabolic functions.


INTRODUCTION
Phosphorus (P) is an essential macronutrient for plant growth, development, metabolism, and regulatory processes (Lopez-Arredondo et al., 2014;Cong et al., 2020). However, inorganic phosphate (Pi), the only form of P that can be acquired by plant roots, is not readily taken up by roots because of its low diffusion rate and instability of forms (Wu et al., 2013). Globally, almost 30% of the cultivated land experiences Pi deficiency (Macdonald et al., 2011). Therefore, a lot of Pi fertilizer is applied to low-Pi soils to increase crop yield; however, crops can acquire, at most 30% of applied Pi, which has led to the excess application of Pi fertilizers with the phenomenon of eutrophication of surface waters, and soil hardening (Correl, 1998;Macdonald et al., 2011;Lopez-Arredondo et al., 2014). In addition, increasing consumption of Pi fertilizer exacerbates the depletion of non-renewable rock phosphate (Johnston et al., 2014). Thus, it is critical to understand the mechanism of Pi homeostasis and to screen Pi-efficient crop genotypes for promoting Pi-use efficiency.
In order to cope with low-Pi conditions, plants have evolved a range of elaborate strategies at the morphological, physiological, and molecular levels to enhance their acquisition efficiency and utilization efficiency of Pi . Especially, these strategies are closely related to root traits including architectural, morphological, physiological, and symbiotic traits (Niu et al., 2012;Sawers et al., 2017;Kafle et al., 2019). In addition, highand low-affinity Pi transporters are involved in Pi uptake and transport activities (Parra-Almuna et al., 2020). For example, the best understood, phosphate transporter 1 (PHT1) is mainly expressed in roots and is involved in the absorption, distribution, and remobilization of phosphate (Zhang et al., 2019a;Parra-Almuna et al., 2020). In addition, changes of root system architecture involved in plant response to P deficiency are regulated by various factors, including phytohormones, sucrose, and transcriptional modification (Chiou and Lin, 2010;Dai et al., 2016). Transcriptomic, proteomic, and metabolomic platforms have been applied to characterize the Pi stress response (Sha et al., 2016;Luo et al., 2018;Ren et al., 2018). A proteomic study showed that many differentially accumulated proteins took part in carbon and energy metabolism, signal transduction, secondary metabolism, and stress defense associated with low-Pi tolerance in barley (Nadira et al., 2016). Recently, posttranscriptional modification was shown to play a key role in the plant response to P deficiency (Park et al., 2014;Yang et al., 2019Yang et al., , 2020Wang et al., 2021). Arabidopsis thaliana PHOSPHATE2 (AtPHO2) and NITROGEN LIMITATION ADAPTATION (AtNLA) as E2 conjugase and E3 ubiquitin ligase enzymes, respectively, were found to degrade several PHT1 members Lin et al., 2013;Park et al., 2014). In rice, decreased phosphorylation levels of OsMAPK6, OsCK2, and five calcium-dependent protein kinases (CDPKs) were found in response to phosphate starvation . The OsPP95, a protein phosphatase type 2C (PP2C), regulates dephosphorylation and transportation of phosphate transporters to maintain P homeostasis in rice .
Protein phosphorylation is widely used to study posttranslational modifications (PTMs) and mainly effects the hydroxyl groups of threonine, serine, and tyrosine in eukaryotic cells (Reinders and Sickmann, 2005). In order to understand the regulatory mechanisms of the Pi tolerance, numerous phosphorylation proteins have been identified in many plants, including A. thaliana, rice, and wheat Yang et al., 2019). Barley (Hordeum vulgare L.) is one of the oldest cereal crops and is ranked fourth in the world cereal production. It is widely used as a major food source, livestock feed, and raw material for malting and brewing (FAO., 2018) and displays wide adaptation to extreme environments. After many years of high-nutrient breeding, crop varieties with good tolerance to limited nutrients are often eliminated, and most research on P is derived from model species, not crops (Oldroyd and Leyser, 2020). Our previous efforts to identify two barley lines with different P efficiency were screened from cultivated barley (Ren et al., 2016), but there is little underlying knowledge of phosphorylation modification under Pi deficiency/resupply in barley.
In this study, a high P-efficiency barley genotype GN121 was used for a comprehensive phosphorylation analysis under Pi deficiency/resupply. We adopted high-resolution liquid chromatography tandem mass spectrometry (LC-MS/MS) linked to highly sensitive immune-affinity antibody analysis and powerful bioinformatics, which were used to identify the substrates and describe the functional characteristics to explore the interplay of phosphorylation in barley roots. To the best of our knowledge, this is the first comprehensive analysis of phosphorylation in barley roots. Our research not only greatly extends the list of phosphorylated proteins in barley but also widens the knowledge of P stress response and defense mechanisms in crops.

Barley Materials and Growth Conditions
The low-Pi-tolerant barley GN121, which reportedly has strong resistance to low Pi stress, was used for phosphoproteomics work (Ren et al., 2018). The nutrient solution and growth conditions of the hydroponic culture are described by Ren et al. (2018). Briefly, surface-sterilized seeds were germinated in a Petri dish on a double-layered filter paper and grown in a greenhouse for 10 days, and then transferred to modified Hoagland hydroponic nutrient solution containing 0.397 mM KH 2 PO 4 (high Pi; +Pi) or 0.0397 mM KH 2 PO 4 (low Pi; -Pi) Pi concentrations. Seedlings (10 days old) were subjected to low Pi for 48 h and then resupplied with high Pi for 48 h. The roots were sampled after Pi starvation for 6 h (P6) and 48 h (P48), and 6 h (R6) and 48 h (R48) of Pi resupply conditions, the roots before treatment were used as control (CK). All roots were directly harvested into liquid nitrogen and stored at −80 • for subsequent protein extraction.

Protein Extraction
The samples were fully ground to powder in liquid nitrogen before being transferred to a 5-ml centrifuge tube. Then, 2 ml lysis buffer [8 M urea, 1% Triton-100, 10 mM dithiothreitol (DTT), and 1% Protease Inhibitor Cocktail (Calbiochem, Darmstadt, Germany)] were added to the powder, and the mixtures were sonicated three times on ice using a high-intensity ultrasonic processor. The remaining debris was removed by centrifugation at 20,000×g for 10 min at 4 • C. The proteins were precipitated with cold 20% trichloroacetic acid for 2 h at −20 • C. After centrifugation at 12,000×g for 10 min at 4 • C, the supernatant was discarded. The remaining precipitate was washed three times with cold acetone. The protein was redissolved in 8 M urea, and the protein concentration was determined with a Pierce BCA Protein Assay Kit (Pierce, Thermo Fisher Scientific, Bonn, Germany) kit according to the instructions of the manufacturer.

Trypsin Digestion
The procedures were as described by Zhong et al. (2017). Dithiothreitol was added to the protein solution to a final concentration of 5 mM and reduced at 56 • C for 30 min. Then, this mixture was alkylated with iodoacetamide to 11 mM for 15 min at room temperature in darkness. The treated protein samples were then diluted by adding 100 mM triethylammonium bicarbonate (TEAB, Sigma, Germany) to urea at <2 M concentration. Finally, trypsin was added at 1:50 trypsinto-protein mass ratio for the first digestion overnight and 1:100 trypsin-to-protein mass ratio for a second 4-h digestion.

Enrichment for Phosphorylated Peptides
To enrich modified peptides, tryptic peptides dissolved in the NETN buffer (100 mM NaCl, 1 mM EDTA, 50 mM Tris-HCl, 0.5% NP-40 and pH 8.0) were incubated with prewashed antibody beads (Lot number 001, PTM BIOLABS, Chicago, IL, USA) at 4 • C overnight with gentle shaking. The beads were washed four times with the NETN buffer and twice with ddH 2 O. The bound peptides were eluted from the beads with 0.1% trifluoroacetic acid. Finally, the eluted fractions were combined and vacuum-dried. The resulting peptides were desalted with C18 Zip Tips (Millipore, Darmstadt, Germany) according to the instructions of the manufacturer for the LC-MS/MS analysis.

The LC-MS/MS Analysis
The tryptic peptides were dissolved in 0.1% formic acid (solvent A), directly loaded onto a home-made reversed-phase analytical column (15 cm length, 75 µm i.d.). The gradient comprised an increase from 6 to 23% solvent B (0.1% formic acid in 98% acetonitrile) over 26 min, 23-35% in 8 min, climbing to 80% in 3 min, and then holding at 80% for the last 3 min, all at a constant flow rate of 400 nl/min on an EASY-nLC 1000 UPLC System (Bruker Daltonics).
The peptides were subjected to nanospray ionization (NSI) source followed by tandem mass spectrometry (MS/MS) in Q Exactive TM Plus (Thermo Fisher Scientific, Bremen, Germany) coupled online to the UPLC system. The electrospray voltage applied was 2 kV. The m/z scan range was 350-1,800 for a full scan, and intact peptides were detected in the Orbitrap at a resolution of 70,000. Peptides were then selected for MS/MS using normalized collision energy (NCE) setting as of 28, and the fragments were detected in the Orbitrap at a resolution of 17,500. A data-dependent procedure that alternated between one MS scan followed by 20 MS/MS scans with 15.0 s dynamic exclusion. Automatic gain control (AGC) was set at 5E4. Fixed first mass was set as 100 m/z.

Database Search
The resulting MS/MS data were processed using the Maxquant search engine (v.1.6.6.0; Max Plank Institute of Biochemistry, Germany). The MS/MS spectra were searched against the human UniProt database (http://www.ebi.ac.uk/GOA/) concatenated with reverse decoy database. Trypsin/P was specified as a cleavage enzyme allowing up to four missing cleavages. The mass tolerance for precursor ions was set as 20 ppm in the First search and 5 ppm in the Main search, and the mass tolerance for fragment ions was set as 0.02 Da. Carbamidom108ethyl on Cys was specified as fixed modification, and acetylation modification and oxidation on Met were specified as variable modifications. The false discover rate was adjusted to <1%, and the minimum score for modified peptides was set >40. A label-free quantification (LFQ) method was based on extracted ion currents (XICs) and was used to quantify phosphorylated protein abundance.

Bioinformatics Annotation Analysis
Gene Ontology (GO) annotation proteome was based on three categories (biological process, cellular component, and molecular function) and derived from the UniProt-GOA database (http:// www.ebi.ac.uk/GOA/). Identified protein domain functional description was annotated by InterProScan (http://www.ebi.ac. uk/interpro/) based on the protein sequence alignment method. Eukaryotes database of Wolfpsort (http://www.genscript.com/ psort/wolf_psort.html) was used to identify the subcellular localization of proteins. Soft MoMo (http://meme-suite.org/ tools/momo), a motif-x algorithm, was used to analyze the model of sequences constituted with amino acids in specific positions of modify-13-mers (six amino acids upstream and downstream of the site) in all protein sequences. STRING database version 10.1 (https://www.string-db.org/) was deployed to analyze the protein-protein interaction (PPI) of the phosphoproteins identified in the current study, applying a confidence score of 0.7 (high confidence).

Identification of Phosphorylated Sites and Proteins in Barley Roots
To investigate the phosphorylation profiling of barley (GN121) roots under Pi deficiency/resupply, a combination of iTRAQbased quantitative proteomic and LC-MS/MS method was used to identify phosphorylated proteins and phosphorylation sites (Figure 1). The mass spectrometry data have been deposited at the ProteomeXchange (http://proteomecentral. proteomexchange.org/cgi/GetDataset) with dataset identifier PXD022053 and PXD022077. A total of 7,710 phosphorylation sites were identified, associated with 3,070 proteins that were quantified (Figure 2A and Supplementary Table 1). The results indicated that 50.1% contained a single phosphorylated site, 22.7% contained two phosphorylated sites, 10.6% contained three phosphorylated sites, and 6.1% contained four phosphorylated sites. Additionally, the phosphorylated sites of 50 (1.48%) proteins were above 10, and there were at least 20 phosphorylated sites in four proteins ( Figure 2B). In addition, we counted the phosphorylation sites on serine, threonine, and tyrosine. A total of 11,538 phosphorylation sites were detected, with 9,715 on serine (84.2%), 1,712 on threonine (14.8%), and 111 on tyrosine (0.9%) sites ( Figure 2C). The distribution patterns of phosphorylation types were consistent with other previous reports on wheat and rice (Zhang et al., 2014a;Zhong et al., 2017).

Functional Characterization and Subcellular Localization of Phosphorylated Proteins in Barley Roots
To better understand the potential roles of phosphorylation in barley roots, the GO analysis was used to assign them to biological process, molecular function, and cellular component (Figure 4 and Supplementary Table 3). In the biological process, most phosphorylated proteins were involved in metabolic processes (33%), response to stimulus (18%), and cellular processes (16%) ( Figure 4A). The most common molecular functions were binding (61%) and catalytic activity (30%) ( Figure 4B). For the cellular component, a majority of phosphorylated proteins were related to organelle (39%), cell (38%), and membrane (17%) ( Figure 4C).

Design of Venn Diagram and Functional Analysis of Phosphorylated Proteins
The change of phosphorylation pattern can be used to explore the potential function of phosphorylation in biological processes (Hou et al., 2017). We identified differentially phosphorylated proteins (DPPs) and sites by at least 1.5 quantification ratio change or more (p < 0.05) within each treatment and constructed a Venn diagram based on DPPs (Figure 5 and Supplementary Table 4).
To further understand the characteristics of DPPs, analyses of GO (biological process, molecular function, and cellular component), Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and protein domains were performed (Figures 6, 7 and Supplementary Tables 5-7). The results showed that the GO terms of DPPs were significantly enriched at P6/CK, including regulation of translational initiation, cellular response to iron ion, PEPC activity, actin cytoskeleton, and nuclear lamina ( Figure 6A). The GO terms of DPPs were significantly enriched at P48/CK, such as ribosomal large subunit assembly, cellular response to aluminum ion, rRNA binding, phosphoenolpyruvate carboxy kinase activity, and PEPC activity ( Figure 6B). Consistently, for the P resupply treatment at R6/CK, the GO terms were mainly enriched in cellular protein complex disassembly, cellular response to aluminum ion, structural constituent of the nuclear pore, nuclear pore outer ring, and RNA binding ( Figure 6C). For R48/CK, the GO terms mainly enriched were 5S rRNA binding, cellular response to interferongamma, kinase binding, rRNA binding, and perinuclear region of cytoplasm ( Figure 6D).
The KEGG pathway analysis showed that the protein synthesis and degradation (RNA transport, mRNA surveillance pathway, and ABC transporters) were significantly enriched for P6/CK, P48/CK, and R6/CK. Furthermore, the metabolic process (nitrogen and purine metabolism) and biosynthesis (monoterpenoid biosynthesis) were significantly enriched for P48/CK and R48/CK (Figures 7A-D). Domain enrichment analysis showed that the most common forms of phosphorylated proteins were the domains of adenosine/AMP deaminase, UBA/TS-N, PB1, Villin headpiece, and La (Figures 7E,F). Above all, these results proved that phosphorylated proteins widely participate in diverse pathways in response to the Pi change.

Kinases and Phosphatase Proteins
Protein kinases and protein phosphatases are two vital antagonistic components involved in protein phosphorylation (Uhrig et al., 2013). In the present study, a total of 40 and 10 DPPs were identified as protein kinases and phosphatases under Pi starvation/resupply, respectively. In detail, under the Pi starvation treatment, 21 kinases and three phosphatases were identified at P6/CK, and 15 kinases and four phosphatases were identified at P48/CK. Consistently, under the Pi resupply treatment condition, 18 kinases and three phosphatases were identified at R6/CK, and 19 kinases and two phosphatases were identified at R48/CK. Interestingly, these protein kinases and protein phosphatases had characteristics of poly-phosphorylation sites; For example, a protein kinase superfamily protein (HORVU2Hr1G110900.9) showed significantly upregulated phosphorylation at a total of six sites under four different treatments, suggesting that it was a potential component of the   Table 8). Some of these enzymes have been shown to be extensively involved in the response to low Pi stress Hu, 2016;Jiang et al., 2018).

Characteristics of the Different Phosphorylation States of One Protein in Response to Different Treatments
In this study, seven DPPs showed specific phosphorylation characteristics in response to both Pi starvation and resupply (Supplementary Table 9

Protein-Protein Interaction Network of Phosphorylated Proteins
To better understand the interactions between phosphorylated proteins, we created PPI for DPPs under Pi starvation and resupply stages, respectively (Figure 8 and Supplementary Table 10). To improve the reliability of PPI analysis, we retained all interactions that had a confidence score ≥0.7 (i.e., high confidence). The results showed that 94 DPPs were mapped to the protein interaction database under Pi starvation ( Figure 8A). Of these, 56 were upregulated and 38 were downregulated, and seven DPPs were identified with node degree ≥20, including 60S ribosomal protein L5-1 (HORVU5Hr1G092630.1), nuclear cap-binding protein subunit 2(HORVU6Hr1G056260.5), and splicing factor 3B subunit 1 (HORVU6Hr1G020850.3). In general, the proteins were enriched for the pathway terms: RNA binding, ribosome, protein   transport, and metal binding. For the Pi resupply stage, there were 84 DPPs as nodes, of which 59 were upregulated and 25 were downregulated, with the proteins mainly enriched in RNA binding and ribosome terms ( Figure 8B). The phosphorylated proteins clearly had close interactions in the RNA-binding and ribosome terms under Pi starvation/resupply, and there were more phosphorylated proteins and metabolic processes involved in interactions for the Pi starvation treatment.

DISCUSSION
Reversible protein phosphorylation plays key roles in response to stress signals and in intracellular coordinating responses (Hu et al., 2015). In plants, the root is the main organ absorbing Pi and perceiving external nutrient stimuli (Marchsner, 2012). Recently, evidence suggests that protein phosphorylation is extensively involved in P absorption and transport responses to Pi starvation/resupply (Mehta et al., 2020;Wang et al., 2020a,b). The present study identified 7,710 phosphorylation sites on 3,373 proteins in barley roots, of which 327 DPPs were identified that were involved in a variety of metabolic processes and biological processes, such as signaling transduction pathway, carbon fixation, and carbohydrate metabolism pathway. The PPI network analysis also indicated that a wide range of interactions was regulated by phosphorylation.

Phosphoproteins Involved in Phytohormone Signaling Pathway
Many phytohormones are involved in P signaling (Rubio et al., 2009;Wendrich et al., 2020). ABA plays a pivotal regulatory role in plant abiotic stress responses and can trigger major changes in plant physiology Raghavendra et al., 2010;de Zélicourt et al., 2016), mediated by the RCAR-PP2C-SnRK2 regulatory modules that are involved in reversible protein phosphorylation regulating processes . The PP2Cs are key negative regulators of ABA-signaling pathways, which can inhibit the activity of SnRKs to cause the signaling cascade (Zhang et al., 2014a). In our study, two phosphatase 2C family proteins were identified, of which, HORVU1Hr1G093230.1 was significantly upregulated for 48 h of Pi starvation, indicating that the protein might actively participate in Pi tolerance by increasing phosphorylation. The HORVU7Hr1G049260.4 was significantly upregulated for 6 h under Pi starvation and resupply conditions, indicating that this PP2C protein had short-term activity for the P response. In addition, the KEGG pathway enrichment analysis showed that these PP2Cs were enriched in phytohormone signal transduction. The ABA receptors PYR/PYL/RCAR inhibited the activity of PP2Cs in an ABA-dependent manner (Rubio et al., 2009). The structure of the PYR1 monomer resembles a folded hand and completely wraps ABA in its central hydrophobic pocket by the interactions of polar and hydrophobic character to induce the PP2Cs to combine on the surface composed of gate and latch loops (Santiago et al., 2009). In the present study, we identified a PYR1(HORVU3Hr1G031380.1) protein that was downregulated at Ser219 under Pi starvation, but which was not identified under Pi resupply condition. We used the SWISS-MODEL server (https://swissmodel.expasy.org/) to predict the 3D protein structure models of PYR1, which showed that the phosphorylated site was present in the C-terminal part of helix α5 (Supplementary Figure 2); previous research showed that the Ser152 mutation caused a reduction in the stability of the ternary ABA-receptor-PP2C complex. Notably, Ser152 is located in the β7-α5 loop and does not contact the ABA molecule ). Thus, we speculate that the Ser219 phosphorylated in helix α5 of PYR1 had a positive function in the ABA signaling pathway.
In addition to ABA, ethylene took a key role in mediating adaptive responses to Pi stress that involved in the morphological changes of the root systems (Iqbal et al., 2013;Dubois et al., 2018). Roots treated with P-deficient stress generally produced more ethylene than P sufficient ones (Gilbert et al., 2000). EIN2 is the central positive regulator in the ethylene signaling pathway, which indicates CTR1 kinase-dependent dephosphorylation at several sites that were triggered by ethylene. Meanwhile, the EIN3/EIL1-dependent ethylene responses were stimulated by endoplasmic reticulum-nucleus translocation (Qiao et al., 2012). Most recently, the target of rapamycin (TOR) kinase has been reported to accurately regulate the phosphorylation of EIN2, and EIN2 protein plays an important role in regulating nutrient metabolism in plants under different nutrient levels and stress (Fu et al., 2021). In our research, the EIN2 (HORVU5Hr1G050330.2) was downregulated in the P6/CK at Ser269, but was upregulated in the R6/CK at Ser462. Moreover, the KEGG pathway enrichment analysis showed that EIN2 was enriched in the MAPK signaling pathway and plant hormone signal transduction. Hence, we speculate that the ethylene response to Pi starvation was triggered by the EIN2 dephosphorylation, but the ethylene reaction was prevented under the Pi resupply condition. Interestingly, phosphorylated EIN2 showed short-term activity at different phosphorylated sites in Pi starvation and resupply.

Phosphoproteins Involved in the Carbohydrate Metabolism
In general, Pi that is absorbed by plant roots is esterified into hexose phosphate and ATP, which can provide energy to drive protein phosphorylation modification (Plaxton and Tran, 2011;Yang et al., 2019). However, P deficiency results in the promotion of glucose metabolism in plants, whereas the content of ATP decreases (Raghothama, 1999;Cai et al., 2012). Sucrose, a major regulator of plants, is transported from shoots to roots via the phloem in response to P starvation and is involved in morphological changes in roots, including lateral root and root hair development (Liu and Vance, 2010;Lei et al., 2011). However, sucrose synthesis is mainly catalyzed by SPS (Wu et al., 2014;Nemati et al., 2018), a plant-specific enzyme, with its activity controlled by a complex feedback mechanism, including phosphorylation/dephosphorylation (Wu et al., 2014). In this study, a SPS-1F (HORVU6Hr1G028330.1) was only identified in the Pi starvation treatment and showed a significantly downregulated phosphorylation level at Ser216, but with no significant difference after Pi resupply. Previous research showed that the phosphorylation sites of spinach and maize SPS were regulated by light/dark cycle and occurred at Ser158 and Ser162, respectively (McMichael et al., 1993;Takahashi et al., 2000). In addition, osmotic pressure also leads to significant changes of phosphorylation level in plant SPS (Toroser and Huber, 1997). Overall, the SPS was activated by dephosphorylation and inactivated by phosphorylation. Therefore, we suggest that SPS dephosphorylation may participate in the regulation of sucrose in barley under Pi starvation. Another example consistent with this phenomenon is TPS, a key enzyme in the trehalose biosynthetic process (Tischler et al., 2013). We identified four TPSs with significant changes in phosphorylation levels, among which HORVU3Hr1G071570.3 and HORVU6Hr1G084190.1 showed specific phosphorylation under Pi starvation condition, whereas HORVU4Hr1G066510.7 showed significantly upregulated specific phosphorylation under Pi resupply. In particular, the TPS (HORVU1Hr1G076480.3) showed phosphorylation at four sites (Ser5, Ser49, Ser70, and Ser71) under Pi starvation and showed a significant upregulation at Ser70 under Pi resupply. We speculate that these proteins were beneficial to the response of roots to P starvation by enhancing phosphorylation at more sites.
In the glycolysis pathway, PEPC is crucial, contributing to the plant tolerance of Pi starvation (Jiang et al., 2018). It catalyzes the metabolism of organic acids in roots to facilitate soil Pi solubilization in the rhizosphere, providing a glycolytic bypass (together with malate dehydrogenase and NAD-malic enzyme) to the ADP-limited cytosolic pyruvate kinase and Pi circulation of PEPC byproduct (Vance et al., 2003;Takahashi-Terada et al., 2005;Plaxton and Podestá, 2006;Chen et al., 2007). Previous studies showed that the serine residue of PEPC was specifically phosphorylated by PEPC kinase, for which the activity was allosterically controlled by second messengers, such as Ca 2+ (Nimmo, 2003;Izui et al., 2004;Boxall et al., 2017). In this study, PEPC1(HORVU5Hr1G055350.3) showed significantly upregulated phosphorylation at Ser10 and Ser18 in the Pi starvation treatment. In support of our findings, PEPCs were researched to be phosphorylated in A. thaliana and Lupinus albus during Pi starvation Chen et al., 2007), and there was no detectable PEPC phosphorylation after Pi resupply (Gregory et al., 2009). In addition, we also identified a PEPC3 (HORVU7Hr1G080510.2) that showed downregulated phosphorylation at Ser13 for both Pi starvation and resupply. Previous studies found that plant PEPCs were activated by phosphorylation of the highly conserved Ser11 residue located near the N termini of PEPC1 and PEPC2, and the process reduced the sensitivity of PEPC to allosteric inhibitors, and increased its affinity to PEP (Li et al., 2020). Comparing PEPC sequences of barley with that in A. thaliana showed that two barley phosphorylated sites of Ser18 in PEPC1 and Ser13 in PEPC3 corresponded to the conserved site Ser11 in A. thaliana (Supplementary Figure 3). Notably, Ser10, a serine phosphorylation site, was identified for 6 h after Pi resupply in our research. There were three serine sites in PEPC1 and four serine sites in PEPC3 of barley that were close to the Ser11 site in A. thaliana. Thus, we speculate that these potential phosphorylation modification sites may play an active role at different stages of P starvation.
Fructose-bisphosphate aldolase is a key enzyme in glycolysis and gluconeogenesis. A previous study in drought-tolerant tomato showed that the gene encoding fructose-bisphosphate aldolase was downregulated by drought stress (Gong et al., 2010). In maize, drought, heat, and combined drought and heat stress resulted in increased phosphorylation levels of fructosebisphosphate aldolase, but decreased protein expression (Hu et al., 2015). A proteomic study in rice roots showed that fructosebisphosphate aldolase, phosphoglycerate mutase, and other related proteins related to glucose metabolism were upregulated by low Pi stress (Hu, 2016). In our study, we found that a fructose-bisphosphate aldolase 2 (HORVU3Hr1G088540.1) had downregulated phosphorylation for 6 h under Pi starvation. We speculated that reversible phosphorylation of fructosebisphosphate aldolase might play an important role in plant adaptation to a wide range of abiotic stress, including Pi starvation.

Phosphoproteins Involved in Signal Transduction
In numerous plant responses to abiotic stresses, Ca 2+ is regarded as the second messenger in intracellular signaling (Davies, 2014). In P-deficient rice, Ca 2+ concentration decreased in roots and increased in shoots and was possibly involved in regulating the Pi starvation response ; in Chlamydomonas reinhardtii, two CDPKs (CDPK1 and CDPK3) were regulated by P and nitrogen starvation (Motiwalla et al., 2014). As Ca 2+ sensors, CDPKs, calmodulin, and calceneurin B-like proteins play important and extensive roles in plant responses to several abiotic stresses (Luan et al., 2002;Hamel et al., 2014). Especially, CDPKs participate in various signal transduction pathways, such as ABA, drought, high salt, oxidative bursts, and osmotic stress (Motiwalla et al., 2014). In the current study, we identified two CDPKs that had significantly changed phosphorylation levels. In detail, CDPK33 (HORVU2Hr1G046660.6) showed upregulated phosphorylation in two sites (Ser13 and Ser32) for 6 h under the Pi starvation condition; CDPK19 (HORVU5Hr1G110900.3) showed increasing phosphorylation trends for 48 h after Pi starvation and resupply, indicating temporal specificity of phosphorylation. The CDPK family protein HORVU6Hr1G066470.4 showed upregulated phosphorylation level at Ser39 for 6 h under the Pi starvation condition, indicating that upregulated phosphorylation of some CDPKs might increase their activity in response to Pi starvation. However, calcium-dependent lipid-binding (CaLB) family protein HORVU2Hr1G036050.1 had downregulated phosphorylation at Ser474 under Pi starvation. In addition, we found that the phosphorylation level of calmodulin-binding transcription activator 1 (HORVU2Hr1G029780.1) and calmodulin-binding transcription activator 2 (HORVU4Hr1G078620.5) showed significantly upregulated phosphorylation at different serine phosphorylation sites under Pi starvation/resupply. These results indicate potential differences between periods and changes in P content in the phosphorylation modification of Pi homeostasis. As CaLB, annexins (such as AtANN1, OsANN1, OsANN3, and OsANN10) may be involved in Ca 2+ channel regulation, peroxidase, and ATPase/GTPase activities in response to abiotic stress (Richards et al., 2013;Qiao et al., 2015;Li et al., 2019;Gao et al., 2020). In this study, annexin 4 (HORVU1Hr1G057170.11) had significantly increased phosphorylation levels under the Pi resupply condition. This indicates that annexin 4 may have a positive function in sensing changing signal responses to Pi resupply. The ubiquitin proteasome system regulates PTM in plants and mainly participates in plant hormone responses and abiotic stress (Liu et al., 2020). However, phosphorylation usually functions as a tuning switch to govern ubiquitination pathways because of the lower energetic cost of a single kinase and ATP relative to ubiquitination, which has a higher energetic cost with its cascade of E1, E2, and E3 enzymes plus ubiquitin (Filipčík et al., 2017). In this study, we identified six ubiquitin-related proteins that had significantly changed phosphorylation level under Pi starvation/resupply, such as an ubiquitin family protein (HORVU1Hr1G052260.3), an ubiquitin-specific protease family C19-related protein (HORVU2Hr1G100660.2), ubiquitin-conjugating enzyme 23 (HORVU3Hr1G004520.12), ubiquitin carboxyl-terminal hydrolase 13 (HORVU3Hr1G075220.5), ubiquitin carboxylterminal hydrolase 26 (HORVU4Hr1G071690.11), and ubiquitin-protein ligase 2 (HORVU5Hr1G033930.8). In particular, ubiquitin-specific protease family C19-related protein and ubiquitin-conjugating enzyme 23 had significantly downregulated and upregulated phosphorylation under Pi resupply, respectively; and ubiquitin family protein (HORVU1Hr1G052260.3) phosphorylation was significantly upregulated at Ser31 at P48/CK, but downregulated at R48/CK. Previous studies have found that ubiquitin phosphorylation occurred on most eligible serine, threonine, and tyrosine residues (Swaney et al., 2013); the phosphorylation of Ser65 of ubiquitin, the phospho-ubiquitin event driving Parkinson's disease key enzymes of PINK1 or Parkin-Ubl is crucial to activate their activities and respond to mitochondrial damage (Ordureau et al., 2014;Pickrell and Youle, 2015). Thus, we speculate that these ubiquitin proteins may trigger a more sophisticated response to P starvation by regulating the phosphorylation level.
Many MAPK cascades are involved in plant responses to P (such as OsMAPK4, OsMAPK6, TaMPK6, and TaMPK1) (Wen et al., 2014;Yang et al., 2019), and the phosphorylation/dephosphorylation of the MAPK cascade is crucial for activation of enzymes that participate in transmitting and amplifying signals from specific environmental stimuli (Jonak et al., 1994). We found that a MAPK1 (HORVU7Hr1G097740.1) had downregulated phosphorylation at Ser193 under Pi resupply.

Phosphoproteins Involved in Cell-Wall Stress
Primary cell walls consist of cellulose cross-linked by hemicelluloses or pectins (Tenhaken, 2015), and previous research has reported the involvement of the cell wall in plant responses to Pi starvation (Zhu et al., 2012). In this study, we found that three cellulose synthase-related proteins [cellulose synthase 1 (CESA1, HORVU0Hr1G038120.1) and cellulose synthase family proteins (HORVU3Hr1G071770.1 and HORVU6Hr1G050750.17)] had upregulated phosphorylation in response to Pi starvation/resupply. Compared with CK, the phosphorylation level of CESA1 showed at least nine quantification ratio changes at Ser226 for P6, P48, and R6. However, we find no significant phosphorylation for 48 h after Pi resupply, suggesting that phosphorylated CESA1 might play important role in the regulation of cellulose synthase in response to Pi starvation. In addition, the phosphorylation of CESA1 plays a critical role in regulating cell elongation, especially for rapidly expanding tissues (Chen et al., 2010). A cell-wall protein AWA1-like (HORVU1Hr1G089280.1) was also identified in P48/CK, with increased phosphorylation at Ser520. These results suggested that phosphorylation participated in cell-wall stress responses in barley roots under Pi starvation.
Many receptor-like kinase members have been found to perceive abiotic stress signals in the cell wall, and a large number of signaling proteins, including ABA, second messengers, reactive oxygen species (ROS), and phosphorylating transcription factors, are involved in the process of transmitting stress signals into the cell (Lindner et al., 2012;Wolf et al., 2012). We found that nine receptor kinases [receptor kinase 2 (HORVU4Hr1G015770.3, HORVU6Hr1G069980.1

Phosphoproteins Involved in Oxidase Stress
The Plant NADPH oxidases participate in responses to various biotic and abiotic stresses by facilitating the production of ROS . In tomato under Pi starvation, the malondialdehyde and ROS levels increased to induce oxidative stress (Zhang et al., 2019b). Recently, a study demonstrated the close relationship between oxidative stress and root growth in barley (Wang et al., 2019). Maintaining root growth is a crucial strategy for plants to adapt to Pi starvation, which requires differential cell-wall synthesis and remodeling, including the activity of antioxidative enzymes (Tenhaken, 2015), such as peroxidase and glutathione S-transferase. We found that a peroxidase superfamily protein (HORVU5Hr1G037160.3) and a glutathione Stransferase family protein (HORVU1Hr1G017050.3) had downregulated phosphorylation under Pi resupply. This result indicated that the change in the phosphorylation level of the antioxidative enzymes may regulate their activity to maintain the ROS balance and protect plants from damage.
Respiratory burst oxidase homolog (rboh) genes can encode NADPH oxidase and have been identified in diverse plant species. In A. thaliana, a total of 10 rboh genes (AtRbohA-J) have been identified: AtrbohB plays a key role in the seed after-ripening; AtRbohC, AtRbohH, and AtRbohJ are essential for polar root hair growth (Mangano et al., 2017); and AtrbohD and AtrbohF mainly function in stress responses and phytohormone signaling (Torres et al., 2002;Kwak et al., 2003;Torres and Dangl, 2005;Desikan et al., 2006). It was reported that the phosphorylation of RbohF induced its activity and triggered the Ca 2+ -ROS signaling network (Kimura et al., 2012). In the present study, a respiratory burst oxidase protein F (HORVU3Hr1G087210.2) showed upregulated phosphorylation at Thr627 under Pi starvation, suggesting that phosphorylation could participate in ROS metabolism under Pi starvation.

Model of the Phosphorylation Mechanisms Underlying Pi Starvation/Resupply in Barley Roots
To deeply research the biological processes of P response regulated by phosphorylation, we mainly analyzed the phosphorylated proteins involved in root growth and development under Pi starvation/resupply (Figure 9). We found that these phosphorylated proteins participated in diverse regulation of the root P responses, including phytohormone signaling, carbohydrate metabolism, signaling translation, cell-wall stress, and oxidases stress. For example, in phytohormone pathway, two phosphoproteins (PP2Cs) were upregulated and two phosphoproteins (PYR1 and EIN2) downregulated under Pi starvation condition. For signaling translation, most proteins were activated by phosphorylation under Pi starvation conditions, and Ca 2+ had one of the most important roles. In particular, phosphoproteins related to cell-wall stress took an active role in the root Pi starvation responses. Additionally, these biological processes involved in the root P responses are linked by complex modes of action to perform their functions, and these need further exploration.

CONCLUSION
This is the first extensive data on phosphorylation in barley roots under Pi starvation/resupply. A total of 7,710 phosphorylation sites corresponding to 3,373 proteins were identified. The phosphoproteins played key roles in diverse biological processes and metabolic pathways, especially in phytohormone signaling, carbohydrate metabolism, and Ca 2+ signaling pathways, and the phosphoproteins involved in signaling pathways differed between Pi starvation and resupply. Our results not only broaden the range of metabolic processes known to be regulated by phosphorylation under Pi starvation/resupply conditions, but also provide new information for examining the functions of phosphorylation in barley roots. However, further studies are needed to reveal the effect of phosphorylation under Pi stress and to explain the potential mechanisms of response to Pi nutrition that are regulated by protein phosphorylation.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. The mass spectrometry data have been deposited at the ProteomeXchange with dataset identifier PXD022053 and PXD022077.

AUTHOR CONTRIBUTIONS
ZM, JW, CL, and PR carried out the proteomic analysis and drafted the manuscript. LY, BL, YM, and XM participated in material culture and performed the statistical analysis. XS and HW conceived of the study and participated in its design. HW, ES, and KY helped to draft the manuscript. All authors read and approved the final manuscript.