Phosphoproteomic analysis of the response of maize leaves to drought, heat and their combination stress

Drought and heat stress, especially their combination, greatly affect crop production. Many studies have described transcriptome, proteome and phosphoproteome changes in response of plants to drought or heat stress. However, the study about the phosphoproteomic changes in response of crops to the combination stress is scare. To understand the mechanism of maize responses to the drought and heat combination stress, phosphoproteomic analysis was performed on maize leaves by using multiplex iTRAQ-based quantitative proteomic and LC-MS/MS methods. Five-leaf-stage maize was subjected to drought, heat or their combination, and the leaves were collected. Globally, heat, drought and the combined stress significantly changed the phosphorylation levels of 172, 149, and 144 phosphopeptides, respectively. These phosphopeptides corresponded to 282 proteins. Among them, 23 only responded to the combined stress and could not be predicted from their responses to single stressors; 30 and 75 only responded to drought and heat, respectively. Notably, 19 proteins were phosphorylated on different sites in response to the single and combination stresses. Of the seven significantly enriched phosphorylation motifs identified, two were common for all stresses, two were common for heat and the combined stress, and one was specific to the combined stress. The signaling pathways in which the phosphoproteins were involved clearly differed among the three stresses. Functional characterization of the phosphoproteins and the pathways identified here could lead to new targets for the enhancement of crop stress tolerance, which will be particularly important in the face of climate change and the increasing prevalence of abiotic stressors.


Introduction
In the past decades, progress on increasing crop yields using semi-empiric breeding and genetics has reached a plateau, which is largely linked to increasingly adverse environmental conditions, especially drought and heat stress (Parent and Tardieu, 2014). In the field, it is often the simultaneous occurrence of several abiotic stressors, which are most lethal to crops. Heat, drought and their combination are the severer stressors for crops and are responsible for most of production losses (Lobell et al., 2011;Suzuki et al., 2014).
Moreover, global climate change may increase the occurrence and distribution of these stressors, causing a further reduction of productivity (Rasul et al., 2011).
Recent studies have demonstrated that plant responses to the combinations of two or more stressors are unique and cannot be directly extrapolated from the response to a single stress. In Arabidopsis response to combination stress, some 61% of the transcriptome changes are not predictable from the response to single stress, and plants prioritized between potentially antagonistic responses for only 5-10% of the responding transcripts (Rasmussen et al., 2013). In wheat response to combination of drought and heat stress, proteomic analysis indicates that few common proteins are observed responding to single and multiple high-temperature events (Yang et al., 2011); we obtained similar results in maize (Hu et al., 2010). So, the simultaneous occurrence of several stress results in highly complex responses of plants; extrapolated the response to combined stresses is largely controlled by different, and sometimes opposing, signaling pathways that may interact with and inhibit each other (Vile et al., 2012;Suzuki et al., 2014).
Plants respond to stress with a wide range of modifications that cause to changes at the morphological, cellular, physiological, biochemical, and molecular levels (Lopes and Reynolds, 2011;Aprile et al., 2013). Overall, protein phosphorylation plays a critical role in regulating many biological functions, including stress responses by signal transduction. Phosphorylation and dephosphorylation can switch many regulatory proteins and enzymes on and off, thus control a wide range of cellular processes and signal relays . In recent years, global quantitative analysis of protein expression and phosphorylation has been performed using iTRAQ-based quantitative proteomic and LC-MS/MS methods (Alvarez et al., 2014;Han et al., 2014), which facilitate the simultaneous detection of changes in protein expression and phosphorylation levels under control and stressed conditions. Large-scale phosphoproteomic analyses have been conducted on crops, especially crops response to stress. In wheat seedlings leaves response to drought stress, some phosphoproteins related to drought tolerance and osmotic regulation exhibit significant phosphorylation changes; there are commonalities and differences of phosphoproteins in different cultivars of bread wheat under drought stress Zhang et al., 2014). In maize leaves response to drought stress, a total of 138 phosphopeptides display highly significant changes and most phosphorylation changes do not reflect protein abundance variation. These proteins influenced epigenetic control, gene expression, cell cycle-dependent processes and phytohormonemediated responses. Bonhomme et al. (2012). These results provide a series of phosphoproteins and phosphorylation sites and a potential network of phosphorylation signaling cascades in wheat seedling leaves. However, little information is available regarding the changes of phosphoproteins and the phosphorylation sites of many stress-responsive protein kinases under combined drought and heat stress. Therefore, characterizing posttranslational modifications of these proteins in crops response to combined stress are very important for understanding the associated signaling pathways and mechanisms for tolerance to stress.
In many regions of the world, maize is an important cereal crop grown mainly in semi-arid environments that are characterized by water scarcity and high temperature, two conditions that usually occur simultaneously in the field. Exposing plants to a combined stress may lead to agonistic or antagonistic responses or cause to some responses that are potentially unrelated to those response to the corresponding single stress. So, to analyze such responses, we performed iTRAQ-based phosphoproteomic analyses in maize exposed to heat, drought and their combination. Furthermore, bioinformatics analyses were conducted to confirm the constitutive nature of the different stress-responsive phosphorylated proteins. This work would provide a basis for further elucidating maize endurance to drought, heat and their combination stress.

Plant Material and Stress Treatments
Maize seeds (Zhengdan 958) were used in the experiments. Zhengdan 958 is a high-yield maize hybrid that is grown in China. The seeds were surface-sterilized for 10 min in 2% hypochlorite, washed in distilled water and germinated on moistened filter paper. Maize seedlings were grown in Hoagland's nutrient solution in a light chamber under 400 µmol m −2 s −1 photosynthetically active radiation, a 14/10-h day/night cycle, a day/night temperature of 28/22 • C, and a relative humidity of 75%. When the fifth leaves were fully expanded, the seedlings were subjected to various treatments.
Drought stress was imposed by placing the seedlings in PEG solution (−0.7 MPa) for 8 h at 28 • C and 40% relative humidity. Heat stress was applied by raising the temperature from 28 to 42 • C at an interval of 2 • C/h and then kept at 42 • C for 1 h, for a total of 8 h. Therefore, each stress treatment lasted 8 h. The combined stress consisted of simultaneous treatment with PEG and heat stress. The control seedlings were kept at 28 • C and 75% relative humidity. Then, the expanding leaves (the fifth from the bottom) of the treated and untreated seedlings were sampled, immediately frozen in liquid N 2 , and stored at −80 • C until analysis. Three biological replicates were performed for each treatment.

Protein Extraction
Total proteins from the fifth newly-expanded leaves of the maize seedlings were extracted according to the method reported by Wang et al. (2013) and Zhang et al. (2014). Approximately 0.5 g of fresh leaves from each biological replicate were ground into a fine power in liquid N 2 in a mortal and further ground in a 4 ml SDS buffer (30% sucrose, 2% SDS, 100 mM Tris-HCl, pH 8.0, 50 mM EDTA-Na 2 , 20 mM DTT) and 4 ml phenol (Tris-buffered, pH 8.0) in a 10 ml tube, followed by the addition of 1 mM phenylmethanesulfonyl fluoride (PMSF) and PhosSTOP Phosphatase Inhibitor Cocktail (one tablet/10 ml; Roche, Basel, Switzerland) to inhibit protease and phosphatase activity. The mixture was thoroughly vortexed for 30 s and the phenol phase was separated by centrifugation at 14,000 × g and 4 • C for 15 min. The upper phenol phase was pipetted into fresh 10 ml tubes and 4-fold volumes of cold methanol plus 100 mM ammonium acetate were added. After centrifugation at 14,000 × g and 4 • C for 15 min, the supernatant was carefully discarded and the precipitated proteins were washed twice with cold acetone. Finally, the protein mixtures were harvested by centrifugation. Using a 2-D Quant Kit (Amersham Bioscience, America) containing bovine serum albumin (BSA) (2 mg/ml) as the standard, we carried out the measurement of protein content. To enhance the quantitative accuracy, extracted proteins from every biological replicate were adjusted to the same concentration for the subsequent analysis.

Protein Digestion and iTRAQ Labeling
Protein digestion was performed according to the FASP procedure described by Wisniewki et al. (2009) andLv et al. (2014), and the resulting peptide mixture was labeled using the 4-plex iTRAQ reagent according to the manufacturer's instructions (Applied Biosystems). Briefly, 200 µg of protein from each sample was mixed with 30 µl of STD buffer (4% SDS, 100 mM DTT, 150 mM Tris-HCl pH 8.0). The detergent, DTT and the other low-molecular-weight components were removed using UA buffer (8 M urea, 150 mM Tris-HCl pH 8.0) with repeated ultrafiltration (Microcon units, 30 kD). Then, 100 µl of 0.05 M iodoacetamide in UA buffer was added to block reduced cysteine residues, and the samples were incubated for 20 min in darkness. The filters were washed three times with 100 µl of UA buffer, then twice with 100 µl of DS buffer (50 mM triethylammoniumbicarbonate at pH 8.5). Finally, the protein suspensions were digested with 2 µg of trypsin (Promega) in 40 µl of DS buffer overnight at 37 • C, and the resulting peptides were collected as a filtrate. The peptide content was estimated via UV absorption at 280 nm using an extinction coefficient of 1.1 per 0.1% (g/l) solution, which was calculated based on the proportion of tryptophan and tyrosine residues in vertebrate proteins.
For labeling, each iTRAQ reagent was dissolved in 70 µl of ethanol and added to the respective peptide mixture. The samples were called control (under no stress), Drought, Heat and DH (combined drought and heat stress) and were labeled with reagent and vacuum dried.

Enrichment of Phosphorylated Peptides using TiO 2 Beads
The labeled peptides were mixed, concentrated using a vacuum concentrator and resuspended in 500 µl of loading buffer (2% glutamic acid/65% ACN/2% TFA). Then, TiO 2 beads were added and then the sample was agitated for 40 min. The sample was centrifuged for 1 min at 5000 g, yielding the first set of beads. The supernatant from the first centrifugation was mixed with more TiO 2 beads, which were treated as before, yielding the second set of beads. Both sets of beads were combined and washed three times with 50 µl of wash buffer I (30% ACN/3%TFA) and then three times with 50 µl of wash buffer II (80% ACN/0.3% TFA) to remove the remaining non-adsorbed material. Finally, the phosphopeptides were eluted with 50 µl of elution buffer (40% ACN/15% NH 4 OH), lyophilized and subjected to MS analysis.

Mass Spectrometry
For nanoLC-MS/MS analysis, 5 µl of the phosphopeptide solution was mixed with 15 µl of 0.1% (v/v) trifluoroacetic acid, then 10 µl of the mixture was injected into a Q Exactive MS (Thermo Scientific) equipped with an Easy-nLC (Proxeon Biosystems, now Thermo Scientific). The peptide mixture was loaded onto a C18 reversed phase column (15 cm long, 75 µm inner diameter, RP-C18 3 µm, packed in-house) in buffer A (0.1% formic acid) and separated using a linear gradient of buffer B (80% acetonitrile and 0.1% formic acid) over 240 min at a flow rate of 250 nl/min, which was controlled by IntelliFlow technology. The peptides were eluted with a gradient of 0-60% buffer B from 0 to 200 min, 60 to 100% buffer B from 200 to 216 min, and 100% buffer B from 216 to 240 min.
For MS analysis, the peptides were analyzed in positive ion mode. MS spectra were acquired using a data-dependent, top-10 method by dynamically choosing the most abundant precursor ions from the survey scan (300-1800 m/z) for HCD fragmentation. Determination of the target value was based on predictive automatic gain control (pAGC). The duration of dynamic exclusion was 40.0 s. Survey scans were acquired at a resolution of 70,000 at m/z 200 and the resolution for HCD spectra was set to 17,500 at m/z 200. The normalized collision energy was 27 eV, and the under fill ratio, which specifies the minimum percentage of the target value that is likely to be reached at the maximum fill time, was defined as 0.1%. The instrument was run with peptide recognition mode enabled.

Data Analysis
MS/MS spectra were searched against the Uniprot_Zea_mays database (62977 sequences, downloaded June 14th, 2013) and a decoy database using Mascot 2.2 (Matrix Science), which was embedded in Proteome Discoverer 1.4. For protein identification, the following options were used: peptide mass tolerance = 20 ppm; MS/MS tolerance = 0.1 Da; enzyme = trypsin; missed cleavage = 2; fixed modification: carbamidomethyl (C); iTRAQ4plex (K); iTRAQ4plex (N-term); variable modification: oxidation (M), phosphorylation (S/T/Y). The score threshold for peptide identification was set at a 5% false discovery rate (FDR). The PhosphoRS site probability was estimated on the probability (0-100%) of each phosphorylation site being truly phosphorylated. The PhosphoRS site probabilities above 75% indicate that a site is truly phosphorylated.

Bioinformatics
The molecular functions of the identified proteins were classified according to their gene ontology annotations and their biological functions. The subcellular localization of the unique proteins identified in this study was predicted using the publicly available program WolfPsort (http://wolfpsort. org). Protein-protein interaction networks were analyzed using the publicly available program STRING (http://string-db.org/). STRING is a database of known and predicted protein-protein interactions. The interactions include direct (physical) and indirect (functional) associations, and they are derived from four sources: the genomic context, high-throughput experiments, coexpression and previous knowledge. STRING quantitatively integrates the interaction data from these sources for a large number of organisms and, where applicable, transfers information between these organisms.
Motif-X online software (http://motif-x.med.harvard.edu/ motif-x.html) was used to find phosphorylation site motifs in the identified maize proteins and to predict the specificity of these motifs based on the identified phosphopeptide sequences. The parameters were set to peptide length = 21, occurrence = 5, and statistical significance for p-values of less than 0.000001.

Statistical Analysis
The phosphoproteins assays were the mean of three replicates. Means were compared by one-way analysis of variance and Duncan's multiple range test at 5% level of significance.

Phosphopeptide Identification Under drought, Heat and Combined Stress
Maize plants at the five-leaf stage were subjected to drought, heat and their combination stress. Multiplex iTRAQ-based quantitative proteomic and LC-MS/MS methods were performed on the total proteins extracted from the fifth newly-expanded leaves, resulting in the identification of 1367 unique phosphopeptides (Figure 1). These phosphopeptides corresponded to 1039 proteins and contained 2171 nonredundant phosphorylation sites, of which 1313 (60.48%) were serine (S) residues, 649 (29.89%) were threonine (T) residues and 209 (9.63%) were tyrosine (Y) residues at an estimated false discovery rate of 5%. In detail, based on a significant linear regression (p < 0.05) and a threshold of ≥1.5-fold or ≤0.66-fold change ratio of stress-induced phosphorylation levels compared with that of the control, the phosphorylation level of 172 phosphopeptides had a significant change under heat, of which 77 were up-regulated and 95 were down-regulated; the phosphorylation level of 149 phosphopeptides had a significant change under drought, of which 69 were up-regulated and 80 were down-regulated; and the phosphorylation level of 144 phosphopeptides had a significant change under the combination stress, of which 70 were up-regulated and 74 were down-regulated. The phosphopeptides with significant changes in phosphorylation levels corresponded to 282 proteins, of which 46 proteins were common under three stress conditions ( Table 1), 69 proteins were common under heat stress and the combined stress (Table S1), 24 proteins were common under drought stress and the combined stress (Table S2), 15 proteins were common under drought and heat stress (Table S3), 75 proteins were only identified under drought stress (Table S4), 30 proteins were only identified under heat stress (Table S5), and 23 proteins were only identified under the combined stress (Table S6).

Heat Shock Proteins
In this study, the phosphorylation levels of seven heat shock proteins (HSPs), including five small HSPs (sHSPs) and two HSP70s, changed significantly under drought, heat or combined both stresses. The phosphorylation levels of three sHSPs (B4G250, BF976 B6T649) were greatly upregulated under heat and the combined stress. However, under drought, the phosphorylation levels of the three sHSPs did not change significantly. It was worth mentioning that the two phosphopeptides of B4G250 differed significantly in phosphorylation level. In contrast to the three sHSPs, the phosphorylation level of the sHSP C0P2N6 was down-regulated by heat and combined stresses but was not significantly affected by drought stress; the phosphorylation level of sHSP B4FR07 was down-regulated by drought stress. Regarding the two HSP70s, the phosphorylation level of B6SZ69 was significantly upregulated by drought stress and down-regulated by heat stress but had no significant change under the combined stress. The phosphorylation level of K7WBH2 was significantly up-regulated by drought stress, down-regulated by heat stress and had no significant change under the combined stress.

Responses of Kinases and Phosphatases to the Three Stresses
It was also worth mentioning the responses of enzymes, including kinases and phosphatases, to stresses. Under heat ( Table 2), 31 phosphopeptides, including 24 different enzymes, were identified, of which five were kinases and three were phosphatases. Under drought stress ( Table 2), 29 phosphopeptides, including 22 unique enzymes, were identified, of which seven were kinases and five were phosphatases. Under the combined stress (Table 2), 29 phosphopeptides corresponding to 21 unique proteins, including five kinases   Table 2).
To elucidate the interactions of the protein kinases/phosphatases and HSPs with other proteins under the three stress treatments, protein-protein interaction analysis of significantly changed phosphoproteins was conducted using STRING software (Figures 2-4). Under drought and the combined stress, AGC_PKAPKG_like.1-ACG kinases (4335426, including homologs of PKA, PKG and PKC) had interactions with transcription factor HY5 (4327123). Protein phosphatase FIGURE 2 | The protein-protein interaction network analysis among significantly changed proteins in maize seedlings under drought using String software. 2C (4341433, 4332374) showed an interactions with initiation factor 2 subunit family domain containing protein (4348536). Under drought, AGC_PKAPKG_like.1-ACG kinases also had interactions with heat shock protein DnaJ (OsJ_01978), initiation factor 2 subunit family domain containing protein (4348536), CBS domain-containing protein (43327739), phospholipase C (4352524) and transcription initiation factor IIF (4348228); protein phosphatase 2C also showed an interaction with heat shock protein DnaJ (OsJ_01978). Under drought and heat stress, a dnaK family protein (HSP70, 4332413) exhibited an interaction with oxygen-evolving enhancer protein 1 (LOC_Os01g31690.1). Under three stress treatments, phospholipase C (4352524) exhibited an interaction with an ethylene-responsive element-binding protein (4349277), and chaperone protein dnaJ 10 (sHSP40, 4346080) interacted with a MYB family transcription factor (4335542). These findings indicated that the AGC_PKA/PKG_like.1-ACG kinases, protein phosphatase 2C and phospholipase C may play an important role in protein substrate phosphorylation/dephosphorylation and that HSPs may play an important protective role as chaperone proteins under drought stress, heat stress and combined both stresses (maize protein query sequences corresponded to rice protein query sequences; see Supplementary Table S7 for drought   stress, Table S8 for heat stress, and Table S9 for the combined stress).

The Characteristics of the Different Phosphorylation States of One Protein in Response to Three Stresses
In this study, 19 phosphoproteins ( Table 3) were found to have different phosphorylated states. Importantly, these peptides had specific phosphorylation characteristics in response to drought, heat and the combination of both stresses. In particular, phosphoenolpyruvate carboxykinase (C0P3W9) had 11 different phosphorylation sites, and the phosphorylation of sAPStPkR, sAPSTPk, sAPStPkRsAPTtPIk, and gEAAAQGAPstPR changed significantly under the three stresses; sAPttPIk, sAPTTPIk, gEAAAQGAPStPR, and eVDYADNsVTENTR changed significantly under heat and the combined stress conditions; rsAPTtPIk and sAPStPkR changed significantly under drought and the combined drought and heat stress conditions; and gGAHsPFAVAISEEER changed significantly only in response to drought stress. The phosphorylation sites of the other 18 proteins were similar to these. This result shows the diversity of the phosphorylation sites and their specificity in response to different stress treatments.   Y stands for significant changes; N stands for no significant changes. "-"indicates that the phosphorylation level was not detected.

Identification of Phosphorylation Motifs in the Phosphopeptides
In the present study, 160, 159, and 161 proteins were identified as having significantly changed phosphorylation levels under drought, heat and their combination, respectively.
To determine whether the phosphorylated versions of these proteins shared common phosphorylation site motifs, Motif-X online software (http://motif-x.med.harvard.edu/motif-x.html) was used to predict the motif specificity of these proteins based on the identified phosphorylation sites. An analysis of serine (S), threonine (T) and tyrosine (Y) residues (Table 4) showed that in response to all the three stress treatments, some peptides had the motifs RRxS and xSPx in common, of which the residues adjacent to the phosphorylated S were enriched for arginine (R) and proline (P); In response to heat stress and the combined stress, the motifs PxTP and RT were common, of which the residues adjacent to the phosphorylated T were enriched for arginine (P) and proline (R). These results are the first to demonstrate a high sensitivity and specificity for threonine sites under heat and combined drought and heat stress.

The Signaling Pathways Associated with Phosphorylated Proteins Under Various Stress Conditions
All identified phosphoproteins were classified using gene ontology (GO) annotation software and were further categorized into three functional groups: molecular function, biological process and cellular component. The results of the GO analyses for drought, heat and the combined stresses are shown in Figures 5-7, respectively. The most common molecular functions were binding activity and catalytic activity, and the most common biological processes were cellular process and metabolic process. Moreover, the most common biological processes and molecular functions that the proteins performed were predicted to occur in the organelles, the cytoplasm and the nucleus. According to biological process analysis using the BLAST2GO program, among these phosphoproteins with significantly changed phosphorylation state, for drought stress, 16 were categorized as "response to stimulus, " 12 were involved in transport and 11 were identified as DNA binding proteins that function as transcription factors ( Figure 5B, Table 5); for heat stress, 12 were categorized as "response to stimulus, " nine were transporter proteins, and 14 were DNA binding proteins that function as transcription factors ( Figure 6B, Table 5); for the combined stress, nine were categorized as "response to stimulus, " 10 were transporter proteins, and eight were DNA binding proteins that function as transcription factors ( Figure 7B, Table 5). In addition, the phosphorylation states of a cell division cycle 5-like, phospholipase C, MYB-CC-type transcription factor, BEL1-type homeodomain protein, bZIP transcription factor superfamily protein, and glycine-rich protein 2 were common for the three stresses; a dnaJ 15-like chaperone protein, disease resistance protein rpp13, and rac GTPase-activating protein 1 were specific to drought; the phosphorylation states of a probable receptor-like protein kinase At5g56460like, translocase of chloroplast chloroplastic-like isoform x1, uncharacterized LOC100501590, carbohydrate transporter sugar porter transporter, probable metal-nicotianamine transporter ysl12-like, uncharacterized LOC100501590, probable receptor-like protein kinase At5g56460-like, protein stichel-like, hypothetical protein ZEAMMB73_938746 and an expressed methyl binding domain-containing protein were specific to heat stress. 11 proteins were specific to the combined stress ( Table 5).
According to KEGG analysis, for drought stress, these phosphoproteins with significantly changed phosphorylation state involved mainly in the protein processing, photosynthesis and carbon metabolism pathways ( Figure 5D). There were 11 proteins (group accessions: Q8W149, C0P9H7, B6SL90,  C0HIN5, B6SZ69, B4FNM4, B6T1H0, B6UIM2, B4FQU2, K7V1I2, B4G250) involved in protein processing, and seven proteins (group accessions: B4FAW3, B4FRF2, P27789, C0P3W9, E9NQE1, C0PD30, and B4FJG1) involved in photosynthesis and carbon metabolism. For heat stress, these phosphoproteins with significantly changed phosphorylation state involved mainly in signaling pathways related to protein processing, inositol phosphate-glycerophospholipid metabolism and photosynthesis ( Figure 6D). Specifically, there were 10 proteins (group accessions: B6SZ69, BF4F976, B6T649, B4G250, BF4F7Z5, Q8W149, B4FQK5, C4J9U8, B6T1H0, and B4FQU2) involved in protein processing, three proteins (group accessions: B6U8P0, B8A1A6, and B4FY17) involved in inositol phosphateglycerophospholipid metabolism, and five proteins (group accessions: B6SQV5, B4FRF2, C0P3W9, C0PD30, B4FKM0) involved in photosynthesis and carbon metabolism. For drought and heat combined stress, these phosphoproteins with significantly changed phosphorylation state involved mainly in signaling pathways related to photosynthesis and carbon metabolism and protein processing ( Figure 7D). In detail, there were eight proteins (group accessions: C0P3W9, B7ZYP6, C0PD30, E9NQE1, B6SKI1, B6SQV5, B7ZYP6, and B4FKM0) involved in photosynthesis and carbon metabolism, and there were eight proteins (group accessions: B6T346, K7TLV1, K7TUM2, K7V1I2, B4FQK5 B6T1H0, B6TG30, K7V792) involved in protein processing. These results indicated that the signaling pathways related to protein processing and to photosynthesis and carbon metabolism played an important role under all three stress conditions. In addition, the differentially phosphorylated proteins that were related to phosphateglycerophospholipid metabolism were primarily involved in the signaling pathways induced by heat stress.

Changes in Receptor Protein Phosphorylation
Receptors are proteins that are either embedded in the plasma membrane or localized to the cytoplasm or nucleus of a cell.
Receptors enable the body to detect changes in the internal or external environment. In this study, the phosphorylation levels of seven receptor proteins significantly changed under the three stresses. The phosphorylation level of a vacuolarsorting receptor 3-like protein (B6U4K3) was reduced by the three stress treatments; that of a probable receptor-like protein kinase At5g56460-like (C0PHB9) was reduced by heat and the combined drought and heat stress; that of the TPA: leucine-rich repeat receptor-like protein kinase family protein (B7ZYR5) was reduced by the individual drought and heat stress treatments, that of the receptor-like protein kinase HERK1-like (C0PND4) was increased by all three stress treatments; and that of the calciumsensing receptor (B6TVL4) and the gibberellin receptor GID1l2 (B6TY90) were increased by the combined stress.

Discussion
Reversible protein phosphorylation is a ubiquitous regulatory mechanism that plays critical roles in transducing stress signals to bring about coordinated intracellular responses. In this study, we report the comprehensive analysis of the phosphorylation changes in maize leaves response to drought, heat and their combination using iTRAQ-based quantitative proteomic and LC-MS/MS methods.

Phosphorylation Regulatory Network in Maize leaves Response to Three Stress Treatments
The interplay between phosphatases and kinases strictly controls biological processes such as metabolism, transcription, cell cycle progression, differentiation, cytoskeletal arrangement and cell movement, apoptosis, intercellular communication, and immunological functions (Johnson, 2009;Pjechová et al., 2014). In this study, five kinases and three phosphatases were identified under heat stress, three kinases and two phosphatases were under drought stress, and three kinases and three phosphatases were under combined heat and drought stress. Nevertheless, according to the analysis of protein-protein interactions among significantly changed phosphoproteins, only AGC_PKAPKG_like.1-ACG kinases, phospholipase C and protein phosphatase 2C were found to have interactions with some phosphoproteins, which involved in gene expression, protein synthesis, and stress tolerance. Our analysis suggested potential multiple phosphorylation regulatory mechanisms of these phosphoproteins for further experimental validation. Furthermore, phospholipase C and protein phosphatase 2C were showed to have a significant phosphoryltion changes under three stress conditions. Previous study also show that the phosphorylation levels of phospholipase C and protein phosphatase 2C is obviously affected Hwang et al., 2014;Wei et al., 2014;Zhang et al., 2014), which indicated that the signal pathway related to phospholipase C and protein phosphatase 2C played an important role in plants response to stress. Besides, we found phosphorylation events in other important kinase. For example, serine threonine-protein kinase wnk4-like (K7TZQ1) changed significantly under three stress conditions. Ser/Thr phosphorylation plays key roles in the regulation of plant growth and development. In wheat response to drought stress, phosphoproteome analysis also reveals two serine threonineprotein kinases . Calcium-dependent protein kinase (C4J038, CDPK) changed under drought stress. In plants, calcium is a ubiquitous second messenger in signal transduction cascades. Most of the known Arabidopsis calcium-stimulated protein kinase activities are related to CDPKs (Cheng et al., 2002). In wheat, among nine phosphorylated CDPKs identified, CDPK7 was phosphorylated in response to drought stress .

Enzymes Involved in the Gluconeogenesis Pathway
Phosphoenolpyruvate carboxykinase (PEPCK, C0P3W9) plays an important role in organic acid metabolism. In the cytoplasm, PEPCK and malate dehydrogenase can synthesize malate from glycolytically derived phosphoenolpyruvate (Fortes et al., 2011). In this study, PEPCK had 7 significantly phosphorylated peptides under heat stress, 5 under drought stress and 8 under the combined stress. Regardless of drought, heat and their combination, the phosphorylation level of these peptides increased or decreased under the same stress. Importantly, although the phosphorylation level changed across the different treatments, the protein expression level did not. Similarly, in drought-stressed Pinus halepensis, PEPCK activity increased without an increase in its transcription or translation (Fontaine et al., 2003;Hýsková et al., 2014). These results demonstrated that PEPC phosphorylation plays an important role in plants' responses to abiotic stress.
Fructose-bisphosphate aldolase is a key enzyme in the pathways of gluconeogenesis. In drought-tolerant tomato response to drought, the gene encoding this enzyme is down-regulated by drought stress (Gong et al., 2010). In the present study, the three stresses all increased the phosphorylation level of this enzyme, and heat and combined drought and heat stress decreased the protein expression. Gluconeogenesis consumes plenty of energy thus down-regulation of the fructosebisphosphate aldolase could inhibit gluconeogenesis for keeping energy in stressed plants.

Protein Degradation by Ubiquitination
Protein degradation via the ubiquitin/26S proteasome system is the main protein degradation pathway; it plays a crucial role in removing misfolded or damaged proteins and in controlling the abundance of certain regulatory proteins during abiotic stress (Vierstra, 2003;Stone, 2014). Ubiquitination is a multistep process involving the sequential action of three enzymes: E1 (ubiquitin activating enzyme), E2 (ubiquitin conjugating enzyme), and E3 (ubiquitin ligase). It has been shown that in plants, certain E3 ubiquitin ligases are involved in transcriptiondependent resistance to high temperature and drought stress (Kim and Kim, 2013;Liu et al., 2014). Saccharomyces bearing a mutation in RSP5, which encodes an essential E3 ubiquitin ligase, are hypersensitive to heat stress (Uesugi et al., 2014). In humans, DNA damage induces the phosphorylation-dependent degradation of the E3 ubiquitin ligase TRIM24 in the nucleus, which disrupts the interaction between TRIM24 and p53 and activates p53 (Jain et al., 2014). However, in the present study, the E3 ubiquitin ligases rglg2-like isoform x1 and upl4-like were dephosphorylated under heat stress and the combined drought and heat stress, but the expression level of rglg2-like changed only minimally (upl4-like was not identified in the protein expression analysis). These results indicate that decreasing the phosphorylation level of the 26S proteasome protein complex can promote the eventual degradation of regulatory proteins that are involved in heat stress.

Phosphoproteins Involved in Water, Sugar and H + Transport
Transporter proteins are important for turgor pressure maintenance and water potential regulation, which are crucial for the growth and survival of plants under biotic and abiotic stress. For example, plasma membrane intrinsic proteins (PIPs) are primary channels that mediate water uptake in plant cells, and they are regulated by phosphorylation. In this study, two aquaporin isoforms were differentially phosphorylated under heat and the combined stress: phosphorylation of the aquaporin PIP2-7 and the integral membrane protein nod26-like were reduced and increased, respectively. Other results also demonstrated that the aquaporin CsPIP2-1, which is involved in phosphorylation-dependent salt-and drought-stress responses in Camelina (Jang et al., 2014), and PIPs from maize shoots are phosphorylated on serine residues by a calcium-dependent kinase both in vitro and in vivo (Van Wilder et al., 2008). Together, these results indicate that phosphorylation plays an important role in the activity of PIPs and, consequently, of water channels.
H + -ATPase, which belongs to the cation transport ATPase (P-type) family, is related to H + electrochemical gradient across the plasma membrane, which regulates cell growth and response to stress (Schaller and Oecking, 1999;Bobik et al., 2010). In Nicotiana tabacum, the phosphorylation of the H + -ATPases leads to the increase of this enzyme activity (Bobik et al., 2010). In wheat, the H + -ATPase identified in two cultivars showed up-regulated phosphorylation levels .

The Sites of Phosphorylation Differ Across the Three Stress Treatments
Analysis of phosphorylation site patterns showed a striking distinction among the three stress treatments, indicating that drought, heat, and the combined drought and heat stress resulted in different sites being phosphorylated ( Table 1, Tables S1-S6). Under the three stress treatments, the phosphopeptides of each protein had different phosphorylation levels and regulatory patterns under the three stress treatments (Table 3).
Further, analysis of these motifs for significant phosphorylation phosphoproteins showed that [SP] and [RxxS] shared in common under three stress treatments, which were showed.
[PxTP] and [RT] shared in common under heat stress and combined double stresses.
[AxS] were specific to combined double stresses. These results provide substantial novel insight into phosphorylation and dephosphorylation and signal transduction; this is the first study to dissect the differences between drought, heat and combined drought and heat stress. Importantly, these results reveal a previously unappreciated quantity and diversity of phosphorylated protein isoforms, and uncover novel signaling pathways.

Phosphorylation Changes of Heat Shock Proteins
Heat shock proteins (HSPs) have "chaperone-like" activities and accumulate in response to various stresses (Eisenhardt, 2013). At least five families of chaperones are found in higher plants: the HSP100 (Clp) family, the HSP90 family, the HSP70 (DnaK) family, chaperonins (GroEL and HSP60), and the small (20-40 kD) HSP (sHSP) family. sHSPs play important and extensive roles in plant defenses against abiotic stress (Mu et al., 2013). The phosphorylation of sHSPs has been demonstrated in maize mitochondria (Lund et al., 2001) and in barley (Hordeum vulgare) seed endosperm (Slocombe et al., 2004). In this study, the phosphorylation level of seven HSPs, including five sHSPs (B4G250, B4F976, B6T649, C0P2N6, and B4FR07) and two HSP70s, changed significantly under drought, heat or combined drought and heat stress. Notably, under heat stress and combined drought and heat stress, the phosphorylation level of the identified sHSPs (except B4FR07) was up-regulated. Taken together, these results indicate that the phosphorylation of sHSPs appears to be important for the regulation of sHSP function in plant responses to heat and to combined drought and heat stress.
These results are the first to demonstrate the phosphorylation of plant sHSPs under the combination of drought and heat stress.

Receptor Proteins
Most proteins that are synthesized on the rough endoplasmic reticulum are delivered to various cellular destinations, including the vacuoles and lysosomes. Such sorting involves the recognition of targeting signals on the proteins by receptors. The vacuolarsorting receptor (VSR) is involved in sorting clathrin-coated vesicles from the Golgi apparatus to the vacuoles. In pumpkin, a putative vacuolar sorting receptor, PV72, undergoes a Ca 2+dependent conformational change (Watanabe et al., 2002). Here, the phosphorylation level of vacuolar sorting receptor 3-like (B6U4K3) was significantly reduced under the three stress treatments, suggesting that stress affects vacuolar sorting in higher plants. The calcium-sensing receptor modulates the cytoplasmic Ca 2+ concentration and is crucial for proper stomatal regulation in response to elevated levels of external Ca 2+ . In this study, a calcium-sensing receptor (B6TVL4) located in the chloroplast thylakoid membrane had significantly increased phosphorylation levels under the combined drought and heat stress, but little is known about the role of calciumsensing receptors in plant responses to stress. Gibberellin (GA) perception is mediated by GID1 (GA-INSENSITIVE DWARF1), a receptor that shows similarity to hormone-sensitive lipases. The discovery of GID1 has yielded new insight into how GA is perceived (Hirano et al., 2008). In this study, the phosphorylation level of the gibberellin receptor GID1L2 (B6TY90) clearly increased under combined drought and heat stress. However, at present, the molecular mechanism of GID1L2's function in GA signaling, and especially in GA-dependent regulation of plant responses to stress, is unclear. The probable receptorlike, At5g56460-like protein kinase; a TPA: leucine-rich repeat receptor-like protein kinase family protein; and the receptor-like protein kinase HERK1-like were also identified, although their mechanisms of action require further study.

Conclusions
Compared to single drought or heat stress, their combination had a different effect on phosphorylation sites and level of phosphoproteins. Of the 282 phosphoproteins identified in the present study, only 46 were found to show similar changes under all stress conditions. Of the 12 phosphorylation motifs identified, only two were common under the three stresses. In particular, these phosphoproteins involved in signaling pathways were different among drought, heat and their combination. Therefore, our results could not only provide new information for understanding the mechanisms of crop tolerance to the combined stress, but also contribute to the identification of crop cultivars with increased tolerance to increasing climate variability.

Author Contributions
XH and WW conceived and designed the research. FZ and NL performed the experiments, DZ, GZ and XH analyzed the data. LW and CL contributed reagents/materials/analysis tools. XH and DZ wrote the paper. All authors read and approved the manuscript.