Global Phosphoproteomic Analysis Reveals the Defense and Response Mechanisms of Jatropha Curcas Seedling under Chilling Stress

As a promising energy plant for biodiesel, Jatropha curcas is a tropical and subtropical shrub and its growth is affected by one of major abiotic stress, chilling. Therefore, we adopt the phosphoproteomic analysis, physiological measurement and ultrastructure observation to illustrate the responsive mechanism of J. curcas seedling under chilling (4 °C) stress. After chilling for 6 h, 308 significantly changed phosphoproteins were detected. Prolonged the chilling treatment for 24 h, obvious physiological injury can be observed and a total of 332 phosphoproteins were examined to be significantly changed. After recovery (28 °C) for 24 h, 291 phosphoproteins were varied at the phosphorylation level. GO analysis showed that significantly changed phosphoproteins were mainly responsible for cellular protein modification process, transport, cellular component organization and signal transduction at the chilling and recovery periods. On the basis of protein-protein interaction network analysis, phosphorylation of several protein kinases, such as SnRK2, MEKK1, EDR1, CDPK, EIN2, EIN4, PI4K and 14-3-3 were possibly responsible for cross-talk between ABA, Ca2+, ethylene and phosphoinositide mediated signaling pathways. We also highlighted the phosphorylation of HOS1, APX and PIP2 might be associated with response to chilling stress in J. curcas seedling. These results will be valuable for further study from the molecular breeding perspective.


Introduction
As one of the most critical limiting factors, low temperature affects the plant growth and development broadly as well as yield, quality, postharvest life and geographic distribution [1]. Chilling tolerance is the ability of a plant to tolerate low temperature (0-15 • C) without injury or damage [2]. Although it is possible to enhance the physical and physiochemical tolerance according to cold acclimation [3], plant species origin from tropical and subtropical areas, such as Oryza sativa, Zea mays and Lycopersicon esculentum, are sensitive to chilling stress and easily damaged by chilling temperature [4]. The cellular structure and physiological characterization of tropical plants were both changed in response to chilling stress, especially the photosynthetic organelle-chloroplast [5][6][7]. Compared to PSI, PSII was likely more sensitive to chilling temperature under moderate light in tropical trees. In order to protect PSI at such stress situation, the PSII photoinhibition appeared and PSII reaction centers was closed [6]. Subsequently, the recovery of PSII from low temperature depended on its ability to maintain PsaA, Cyt b6/f and D1 protein at photoinhibitory conditions [8]. Low temperature stress induced the changes of a variety of protein kinases and transcription factors in plants [5,9,10]. It was reported that increases in the cytosolic transient calcium flux played a vital

Ultrastructure Change of J. curcas Seedling under Chilling Stress
To assess the damage of J. curcas seedling under chilling treatment, the fourth-leaves at each treatment (C0 h, C6 h, C24 h and R24 h) were prepared for ultrastructural observation. In the leaf cell from C0 h, the morphological structures were intact chloroplasts with numerous embedded starch granules, which mainly distributed closely to the cell membrane. Besides, intact vacuole and mitochorina could also be observed clearly (Figure 2, C0 h). As a universal symptom, the obvious manifestations of chilling stress were chloroplast swelling, a distortion of vacuole, a reduction in the size of starch granules (Figure 2, C6 h), prolonged chilling treatment leaded to grana unstacked and starch granules continued to diminish with time till disappeared completely ( Figure 2, C24 h). However, the mitochorina from the J. curcas seedling showed relative stationary and was not as sensitive as chloroplast and vacuole in the whole chilling treatment process ( Figure 2). After recovery for 24 h, the chloroplasts of J. curcas seedling returned to normal morphological status ( Figure 2, R24 h).

Ultrastructure Change of J. curcas Seedling under Chilling Stress
To assess the damage of J. curcas seedling under chilling treatment, the fourth-leaves at each treatment (C0 h, C6 h, C24 h and R24 h) were prepared for ultrastructural observation. In the leaf cell from C0 h, the morphological structures were intact chloroplasts with numerous embedded starch granules, which mainly distributed closely to the cell membrane. Besides, intact vacuole and mitochorina could also be observed clearly (Figure 2, C0 h). As a universal symptom, the obvious manifestations of chilling stress were chloroplast swelling, a distortion of vacuole, a reduction in the size of starch granules (Figure 2, C6 h), prolonged chilling treatment leaded to grana unstacked and starch granules continued to diminish with time till disappeared completely ( Figure 2, C24 h). However, the mitochorina from the J. curcas seedling showed relative stationary and was not as sensitive as chloroplast and vacuole in the whole chilling treatment process ( Figure 2). After recovery for 24 h, the chloroplasts of J. curcas seedling returned to normal morphological status ( Figure 2, R24 h).

Phosphoprotein Identification and Phosphorylated Site Location
In total, 3101 phosphopeptides with 3101 phosphorylated sites corresponding to 1810 phosphoproteins were identified (Table S1) and the proportions of pS, pT and pY sites were calculated as 89.4%, 9.8% and 0.8%, respectively ( Figure S1). The mass spectrometry proteomics data have been deposited in the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the PRIDE partner repository [46] with the dataset identifier PXD011438. In order to evaluate analytical reproducibility, a range of quality control measures were taken for the three biological replicates of each condition before the comparison analysis of the phosphorylation levels between the chilling stress condition and the control. The result of Pearson correlation analysis showed that the four samples were individually clustered confidently with their replicates ( Figure S2). Only the phosphopeptides identified from all biological repeats were used for further analysis.

Screening Phosphoproteins with Phosphorylation Level Significantly Changed
The intensity of each phosphopeptide was normalized as the Zhang et al. described [28]. According to the ANOVA analysis, 996 phosphorylated sites corresponding to 805 phosphoproteins were screened out (Table S2). Compared to untreated J. curcas seedling (C0 h), 308, 332 and 291 phosphorylation sites corresponding to 279, 313 and 270 phosphoproteins were differentially changed after chilling for 6, 24 h (C6 h, C24 h) and recovery for 24 h (R24 h), respectively (Tables S3-S5 and Figure 3A). There were 44, 76, 76 and 112 phosphorylation sites corresponding to induced, up regulated, down regulated and depressed regulation after 6 h chilling treatment, prolonged the chilling treatment time to 24 h, the number of phosphorylation site for each regulation has changed, which are 38, 94, 67 and 133, respectively ( Figure 3A). To be mentioned, when recovery for 24 h, the number of induced phosphorylation sites increased obviously to 76 while the depressed ones decreased to 63 ( Figure 3A). Further venny charts indicated that a number of differential phosphorylation sites were overlapped between C6 h and C24 h when compared to C0 h, however, only few differential phosphorylation sites were overlapped no matter the comparison between C6 h and R24 h or C24h and R24 h ( Figure 3B).

Phosphoprotein Identification and Phosphorylated Site Location
In total, 3101 phosphopeptides with 3101 phosphorylated sites corresponding to 1810 phosphoproteins were identified (Table S1) and the proportions of pS, pT and pY sites were calculated as 89.4%, 9.8% and 0.8%, respectively ( Figure S1). The mass spectrometry proteomics data have been deposited in the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the PRIDE partner repository [46] with the dataset identifier PXD011438. In order to evaluate analytical reproducibility, a range of quality control measures were taken for the three biological replicates of each condition before the comparison analysis of the phosphorylation levels between the chilling stress condition and the control. The result of Pearson correlation analysis showed that the four samples were individually clustered confidently with their replicates ( Figure S2). Only the phosphopeptides identified from all biological repeats were used for further analysis.

Screening Phosphoproteins with Phosphorylation Level Significantly Changed
The intensity of each phosphopeptide was normalized as the Zhang et al. described [28]. According to the ANOVA analysis, 996 phosphorylated sites corresponding to 805 phosphoproteins were screened out (Table S2). Compared to untreated J. curcas seedling (C0 h), 308, 332 and 291 phosphorylation sites corresponding to 279, 313 and 270 phosphoproteins were differentially changed after chilling for 6, 24 h (C6 h, C24 h) and recovery for 24 h (R24 h), respectively (Tables S3-S5 and Figure 3A). There were 44, 76, 76 and 112 phosphorylation sites corresponding to induced, up regulated, down regulated and depressed regulation after 6 h chilling treatment, prolonged the chilling treatment time to 24 h, the number of phosphorylation site for each regulation has changed, which are 38, 94, 67 and 133, respectively ( Figure 3A). To be mentioned, when recovery for 24 h, the number of induced phosphorylation sites increased obviously to 76 while the depressed ones decreased to 63 ( Figure 3A). Further venny charts indicated that a number of differential phosphorylation sites were overlapped between C6 h and C24 h when compared to C0 h, however, only few differential phosphorylation sites were overlapped no matter the comparison between C6 h and R24 h or C24h and R24 h ( Figure 3B).

Phosphoprotein Identification and Phosphorylated Site Location
In total, 3101 phosphopeptides with 3101 phosphorylated sites corresponding to 1810 phosphoproteins were identified (Table S1) and the proportions of pS, pT and pY sites were calculated as 89.4%, 9.8% and 0.8%, respectively ( Figure S1). The mass spectrometry proteomics data have been deposited in the ProteomeXchange Consortium (http://proteomecentral.proteomexchange.org) via the PRIDE partner repository [46] with the dataset identifier PXD011438. In order to evaluate analytical reproducibility, a range of quality control measures were taken for the three biological replicates of each condition before the comparison analysis of the phosphorylation levels between the chilling stress condition and the control. The result of Pearson correlation analysis showed that the four samples were individually clustered confidently with their replicates ( Figure S2). Only the phosphopeptides identified from all biological repeats were used for further analysis.

GO Annotation of Significant Changed Phosphoproteins
All the identified significant changed phosphoproteins (805) were used for GO annotation. The distribution pie charts for biological process, cellular component and molecular function are shown in Figure 4. In the whole treatment process, the cellular protein modification process, transport, cellular component organization and signal transduction were significantly overrepresented from the biological process perspective. Cytoskeleton and ribosome were significantly overrepresented from the cellular component perspective. Nucleotide binding, protein binding and kinase activity were significantly overrepresented from the molecular function perspective (Figure 4).

GO Annotation of Significant Changed Phosphoproteins
All the identified significant changed phosphoproteins (805) were used for GO annotation. The distribution pie charts for biological process, cellular component and molecular function are shown in Figure 4. In the whole treatment process, the cellular protein modification process, transport, cellular component organization and signal transduction were significantly overrepresented from the biological process perspective. Cytoskeleton and ribosome were significantly overrepresented from the cellular component perspective. Nucleotide binding, protein binding and kinase activity were significantly overrepresented from the molecular function perspective (Figure 4).

GO Annotation of Significant Changed Phosphoproteins
All the identified significant changed phosphoproteins (805) were used for GO annotation. The distribution pie charts for biological process, cellular component and molecular function are shown in Figure 4. In the whole treatment process, the cellular protein modification process, transport, cellular component organization and signal transduction were significantly overrepresented from the biological process perspective. Cytoskeleton and ribosome were significantly overrepresented from the cellular component perspective. Nucleotide binding, protein binding and kinase activity were significantly overrepresented from the molecular function perspective (Figure 4).

Conservation Analysis of the Significant Changed Phosphoproteins
The sequences of the 805 significant changed phosphoproteins were used as queries to blast phosphoprotein databases that were constructed using data sets in the Plant Protein Phosphorylation DataBase (P3DB) [47] and PhoPhAt [48]. O. sativa and A. thaliana were compared against J. curcas to determine the degree of conservation of phosphorylation sites among different plant species. The thresholds were set as score ≥ 80, E-value < 1 × 10 −10 and identity ≥ 30%. In all, 561 (69.7%) of the 805 phosphoproteins had phosphorylated orthologs in the two species, 163 (20.2%) had phosphorylated orthologs in only one species ( Figure 5 and Table S6), 81 phosphoproteins had no phosphorylated orthologs in both of the two species (Table S6).

Conservation Analysis of the Significant Changed Phosphoproteins
The sequences of the 805 significant changed phosphoproteins were used as queries to blast phosphoprotein databases that were constructed using data sets in the Plant Protein Phosphorylation DataBase (P3DB) [47] and PhoPhAt [48]. O. sativa and A. thaliana were compared against J. curcas to determine the degree of conservation of phosphorylation sites among different plant species. The thresholds were set as score ≥ 80, E-value < 1 × 10 −10 and identity ≥ 30%. In all, 561 (69.7%) of the 805 phosphoproteins had phosphorylated orthologs in the two species, 163 (20.2%) had phosphorylated orthologs in only one species ( Figure 5 and Table S6), 81 phosphoproteins had no phosphorylated orthologs in both of the two species (Table S6).

Conservation Analysis of the Significant Changed Phosphoproteins
The sequences of the 805 significant changed phosphoproteins were used as queries to blast phosphoprotein databases that were constructed using data sets in the Plant Protein Phosphorylation DataBase (P3DB) [47] and PhoPhAt [48]. O. sativa and A. thaliana were compared against J. curcas to determine the degree of conservation of phosphorylation sites among different plant species. The thresholds were set as score ≥ 80, E-value < 1 × 10 −10 and identity ≥ 30%. In all, 561 (69.7%) of the 805 phosphoproteins had phosphorylated orthologs in the two species, 163 (20.2%) had phosphorylated orthologs in only one species ( Figure 5 and Table S6), 81 phosphoproteins had no phosphorylated orthologs in both of the two species (Table S6).

Analysis of Phosphorylation Motifs of Significant Changed Phosphopeptides
The kinase related phosphorylation motifs of the significantly changed phosphopeptides were identified by employing WebLogo and motif-X. Briefly, the significant changed phosphopeptides were centered at the phosphorylated amino acid residues of each experimental group (Table S7) and then were submitted for Weblogo analysis and phosphorylation motif extraction. Nine phosphorylation motifs were enriched from the four experimental groups ( Figure 6). Those phosphorylation motifs were then searched in relevant databases to find the specific protein kinases.
[sP] and [tP] motifs were the proline-directed motifs, which were potential substrates of MAPK [49].
[sF] contained the phenylalanine residue and was the minimal MAPK target motif [30]. The kinase related phosphorylation motifs of the significantly changed phosphopeptides were identified by employing WebLogo and motif-X. Briefly, the significant changed phosphopeptides were centered at the phosphorylated amino acid residues of each experimental group (Table S7) and then were submitted for Weblogo analysis and phosphorylation motif extraction. Nine phosphorylation motifs were enriched from the four experimental groups ( Figure 6). Those phosphorylation motifs were then searched in relevant databases to find the specific protein kinases.
[sP] and [tP] motifs were the proline-directed motifs, which were potential substrates of MAPK [49]. [LxRxxs] was basic motif representative of CDPK substrate and the motif [Rxxs] was a potential substrate for CDPK-II, which was also the 14-3-3 binding motif [52][53][54]. The two motifs [Rxxs] and [Kxxs] could be assigned to motif [-(K/R)-x-x-(pS/pT)-] and were identified as a phosphorylation motif of the SnRK2 or CDPK in plants by previous study [55]. [sF] contained the phenylalanine residue and was the minimal MAPK target motif [30].

Protein-Protein Interaction (PPI) Analysis of Significant Changed Phosphoproteins
The PPI network of the significant changed phosphoproteins identified in the current study were analyzed by STRING (http://string-db.org, version 9.1). A total of 514 KOGs representing 610 phosphoproteins (Table S8) were used to construct the PPI network. In order to improve the reliability of the PPI analysis, the confidence score was set at the highest level (≥0.900). Finally, a complex PPI network that contained 319 nodes and 1924 edges was displayed through Cytoscape ( Figure S3). With the aim to further extract the key potential interacting proteins from the whole PPI network, 85 KOGs representing 183 significantly changed phosphoproteins, which related to signal transduction, posttranslational modification and intracellular trafficking, transport (Tables S9 and  S10), were chosen and centered to construct the subnetwork (Figure 7). The result showed that KOG0583 (protein id: 897 and 2462), KOG0841 (protein id: 1724 and 2765) and KOG0070 (Protein id: 5) were the centered phosphoproteins in each functional group (Table S9 and Figure 7).

Protein-Protein Interaction (PPI) Analysis of Significant Changed Phosphoproteins
The PPI network of the significant changed phosphoproteins identified in the current study were analyzed by STRING (http://string-db.org,version9.1). A total of 514 KOGs representing 610 phosphoproteins (Table S8) were used to construct the PPI network. In order to improve the reliability of the PPI analysis, the confidence score was set at the highest level (≥0.900). Finally, a complex PPI network that contained 319 nodes and 1924 edges was displayed through Cytoscape ( Figure S3). With the aim to further extract the key potential interacting proteins from the whole PPI network, 85 KOGs representing 183 significantly changed phosphoproteins, which related to signal transduction, posttranslational modification and intracellular trafficking, transport (Tables S9 and S10), were chosen and centered to construct the subnetwork (Figure 7). The result showed that KOG0583 (protein id: 897 and 2462), KOG0841 (protein id: 1724 and 2765) and KOG0070 (Protein id: 5) were the centered phosphoproteins in each functional group (Table S9 and  . Protein-protein interaction (PPI) network of signal transduction, posttranslational modification and intracellular trafficking, transport related phosphoproteins by STRING. Nodes with orange, red and blue background color represent the KOGs of differential phosphoproteins related with signal transduction, posttranslational modification and intracellular trafficking, transport, respectively.

An Overview of Response and Defense Mechanisms of J. curcas Seedling under Chilling Stress
Based on the above results, 111 phosphoproteins (Table 1) with significant changes (with credible ANOVA analysis and fold change ≥ 2) were chosen to figure out a systematic chilling response and defense pathway in J. curcas seedling ( Figure 8). The ion stress signal was transferred by chilling sensors on plasma membrane into cells and led to increase of Ca 2+ . Subsequently, ABA, ethylene, MAPK, Phosphatidylinositol, CDPK and 14-3-3 signal pathways were activated by phosphorylation modification in J. curcas seedling under chilling stress. These signals were then transferred into nucleus and induced the expression changes of response and defense related genes. Under chilling stress, photosynthesis was depressed, which resulted in the proteins associated with PSI and PSII to be significantly changed at the phosphorylation level. Channels and transporters on the membrane which were associated with ion, auxin, H2O2 also regulated through phosphorylation or dephosphorylation. According to change at phosphorylation level, the misfolded proteins were possibly handled either by refolding or degradation. Therefore, an ubiquitination mediated degradation pathway centered on E3s was displayed in the schematic representation. Protein-protein interaction (PPI) network of signal transduction, posttranslational modification and intracellular trafficking, transport related phosphoproteins by STRING. Nodes with orange, red and blue background color represent the KOGs of differential phosphoproteins related with signal transduction, posttranslational modification and intracellular trafficking, transport, respectively.

An Overview of Response and Defense Mechanisms of J. curcas Seedling under Chilling Stress
Based on the above results, 111 phosphoproteins (Table 1) with significant changes (with credible ANOVA analysis and fold change ≥ 2) were chosen to figure out a systematic chilling response and defense pathway in J. curcas seedling ( Figure 8). The ion stress signal was transferred by chilling sensors on plasma membrane into cells and led to increase of Ca 2+ . Subsequently, ABA, ethylene, MAPK, Phosphatidylinositol, CDPK and 14-3-3 signal pathways were activated by phosphorylation modification in J. curcas seedling under chilling stress. These signals were then transferred into nucleus and induced the expression changes of response and defense related genes. Under chilling stress, photosynthesis was depressed, which resulted in the proteins associated with PSI and PSII to be significantly changed at the phosphorylation level. Channels and transporters on the membrane which were associated with ion, auxin, H 2 O 2 also regulated through phosphorylation or dephosphorylation. According to change at phosphorylation level, the misfolded proteins were possibly handled either by refolding or degradation. Therefore, an ubiquitination mediated degradation pathway centered on E3s was displayed in the schematic representation.

Discussion
We performed a comprehensive analysis of J. curcas seedling under chilling stress by employing the phosphoproteomic approach. Phosphorylation changes frequently occurred under chilling stress and the number of significantly changed phosphoproteins continually increased within the chilling duration ( Figure 3 and Tables S2-S5). As a low temperature-sensitive species, J. curcas appeared to respond to and defend against chilling stress mainly through reversible protein phosphorylation, which was stimulated by Ca 2+ , phytohormones and phosphoinositide. The whole regulated network illustrated in this study mainly included signal transduction, ion transport and protein ubiquitination.  Table 1. Black solid lines with arrows represent catalytic reaction, while dotted lines with arrows represent material transport. Blue solid lines with arrows represent direct relationship, while magenta lines with arrows represent ubiquination. Black solid lines with stubs represent negative regulation. Green dotted lines with arrows in chloroplast represent transport between PSII and PSI. The proteins labeled with red colors represent significantly changed phosphoproteins among C0 h, C6 h, C24 h and R24 h. The area outlined with dotted square represents ER associated degradation (ERAD).

Discussion
We performed a comprehensive analysis of J. curcas seedling under chilling stress by employing the phosphoproteomic approach. Phosphorylation changes frequently occurred under chilling stress and the number of significantly changed phosphoproteins continually increased within the chilling duration ( Figure 3 and Tables S2-S5). As a low temperature-sensitive species, J. curcas appeared to respond to and defend against chilling stress mainly through reversible protein phosphorylation, which was stimulated by Ca 2+ , phytohormones and phosphoinositide. The whole regulated network illustrated in this study mainly included signal transduction, ion transport and protein ubiquitination. At the early chilling treatment period, protein phosphorylation might involve in the regulation of ion transport for homeostasis. Prolonging the time of chilling treatment, the physiological injury appeared, which resulted in the proteins related chloroplast movement were phosphorylated against chilling ( Figure 2 and Table 1).

SnRKs Played a Central Role in the Chilling Responsive Signal Pathways
When plants were exposed to chilling stress, the osmotic stress mediated by ABA accompanied and the corresponding signal pathways were triggered [4]. According to the PPI analysis (Figure 7), KOG0583 represented SnRKs (Table S9, protein id: 897 and 2462) showed the central position, which involved in the ABA mediated signal pathway and seemed to act as a vital defense approach against chilling stress in J. curcas seedling. JcSRK2a (id 2462) were observed to be significantly up-regulated phosphorylated in the chilling treatment process and return to normal condition at the recovery process, which might indicate that the JcSRK2a phosphorylation was activated by osmotic stress caused by chilling treatment. Multiple sequence alignments between JcSRK2a and 10 subfamily of SnRK2 from A. thaliana showed the JcSRK2a shared the highest similarity (86.44%) with AtSnRK2.4 which could be activated by osmotic stress [56]. Furthermore, overexpression of TaSnRK2.4 in A. thaliana enhanced the plant tolerance to drought, salinity and low temperature [57]. The JcSRK2a (KDP38831.1) shared 79% homology with TaSnRK2.4 (ACU65228.1) and the phosphorylation site Ser348 of JcSRK2a was relative conserved comparing to the site of TaSnRK2.4 ( Figure S4). These results provided a possible clue for us to improve low temperature resistance of J. curcas. In other hand, as the member of the SnRK1 subfamily [58], KIN10 is one of central transcriptional integrators of stress and energy signaling [59,60]. We observed KIN10 (id 897) was also observed to be phosphorylated changed as the SnRK2a showed. In consideration of these results, we hypothesized the phosphorylation of SnRK2a and KIN10 mainly involved in defense against chilling stress and promoting catabolism respectively in J. curcas seedling.
ABI5 was highly inducible by either ABA treatment or stress [61]. In the presence of ABA, ABI5 accumulation was promoted by inducing KEG degradation, which might be accomplished by KEG self-phosphorylation [62]. In this study, the phosphorylation level of ABI5-2 (id 58) shared the similar change pattern as the SRK2a (id 2462) displayed ( Figure S5), which might suggest the phosphorylation of ABI5-2 and SRK2a were both triggered by ABA. Additionally, ABI5 also underwent sumoylation at K391 residue to protect it from degradation [19], thereby possibly accelerating ABA-mediated growth inhibition in J. curcas seedling under chilling stress. Coincidently, as a negative regulator, JcKEG (id 2703) showed down regulation trend at phosphorylation level, which possibly indicated that ABA had induced KEG degradation. This assumption was supported by previous studies [18,62]. In response to ABA, TRAB1 phosphorylation at Ser 102 was essential and the TRAB1 phosphorylation level increased in a very short period and declined thereafter in O. sativa [16,63]. The JcTRAB1 (id 2400) phosphorylation level changed in a similar trend with OsTRAB1 (Table 1). Therefore, we compared the phosphorylation site of JcTRAB1 (KDP38330.1, id 2400) with OsTRAB1 (BAA83740.1). The result showed the phosphorylation site of JcTRAB1 was different with OsTRAB1, however, these two phosphorylation sites were very-close and both of them were in the conserved region II (Figure 9). Moreover, the phosphorylation motifs of JcTRAB1 and OsTRAB1 were [RxxS] and [LxRxxs] respectively, which were potential substrates of SnRK2 or CDPK in plants [55]. It was assumed the JcTRAB1 phosphorylation at Ser in this conserved region was possibly related to the ABA response.
OsTRAB1 (BAA83740.1). The result showed the phosphorylation site of JcTRAB1 was different with OsTRAB1, however, these two phosphorylation sites were very-close and both of them were in the conserved region II (Figure 9). Moreover, the phosphorylation motifs of JcTRAB1 and OsTRAB1 were [RxxS] and [LxRxxs] respectively, which were potential substrates of SnRK2 or CDPK in plants [55]. It was assumed the JcTRAB1 phosphorylation at Ser in this conserved region was possibly related to the ABA response. The neighbors of SnRKs were chilling stress related protein kinases, which participated in multiple biological processes. Cytosolic Ca 2+ transiently increased as the result from chilling stress [64]. CDPK played a key role in abiotic stress and mediated in Ca 2+ signal transduction [65]. The phosphorylation of OsCDPK13 had been indicated as an important signal event in response to cold The neighbors of SnRKs were chilling stress related protein kinases, which participated in multiple biological processes. Cytosolic Ca 2+ transiently increased as the result from chilling stress [64]. CDPK played a key role in abiotic stress and mediated in Ca 2+ signal transduction [65]. The phosphorylation of OsCDPK13 had been indicated as an important signal event in response to cold stress [66]. In this study, JcCDPK8 (id 290) was inducible phosphorylated at Ser 43 only under chilling stress (Table 1), which might indicate that JcCDPK8 functioned in the phosphorylation status in response to chilling stress. In ethylene mediated signal pathway, ethylene responses were negatively regulated by receptors family [67]. Two proteins EIN4 (id 521) and EIN2 (id 1204) were phosphorylated changed. EIN4 and EIN2 were responsible for the negative and positive regulations respectively in ethylene mediated signal [68]. Although it was unclear whether the protein expression level of EIN4 and EIN2 had changed, their changes at phosphorylation level might suggest they were likely related to the ethylene signaling. Additionally, phosphorylation changes were also implicated in MAPK mediated signaling. Phosphorylated proteins involved in MAPK signaling including EDR (id 2304), MKP1 (id 2287) and HOS1 (id 2597) were all significantly changed at phosphorylation level. MAPKs were implicated in response to low temperature condition [69][70][71][72], as well as interacted with MKP1 [73]. In this study, we did not detect the significantly change of MAPK at phosphorylation level. However, MKP1 (id 2287) showed an elevated phosphorylation level under chilling treatment (id 2287), suggesting a possible link with MAPK function in J. curcas. MAPK subsequently interacted with ICE1, which was degraded by an E3 ligase HOS1 [21]. The phosphorylation level of HOS1 (id 2597) was increasing under chilling stress, possibly resulting in ICE1 degradation. It might be a possible reason to explain why no significantly changed ICE1 at the phosphorylation level was detected in this study.

Phosphorylation Participated in the Regulation of Phosphoinositide Metabolism
As the primary lipid-derived signals, phosphoinositides were involved in various plant responses to surrounding environment [74]. A total of 13 participants in phosphoinositide metabolism were differentially phosphorylated in this study (Table 1). They participated in reactions leading to the generation of different phosphoinositide species. Among them, PLC2 was the decisive enzyme in phosphoinositide metabolism [75], it catalyzed the hydrolysis of PI(4,5)P2 into DAG and IP3, which subsequently acted as two important second messenger molecules and resulted in release of Ca 2+ into the cytoplasm [76]. Recent studies showed the expressions of JcPLC2 and BpPLC2 were constitutive even under low temperature stress [27,77]. In J. curcas seedling, PLC2s (id 3110, 3180) were both differential phosphorylated at recovery stage instead of the chilling stage. Their phosphorylation sites were Ser312 and Ser278 respectively. The phosphorylation site Ser312 from PLC2 (KDP44118.1, id 3110) was located in the catalytic Y domain ( Figure S6), which have been reported in AtPLC2 mainly handled through degradation by ubiquitin-proteasome system at continuous chilling (24 h) and recovery period.

Photoinhibition and Chloroplast Movement
In response to the chilling stress, changes could be observed at the physiological and cellular level in J. curcas seedling within 24 h chilling treatment (Figure 1; Figure 2), which was consistent with the previous description [5]. Chloroplast was the first and most severely affected organelle when plants were exposed to freezing stresses. The photosyntic rate of J. curcas seedling was dramatically decreased ( Figure 1B) and the chloroplasts were swelling after 24 h chilling treatment (Figure 2), which indicated serious damage of the photosystem in J. curcas seedling. At chilling temperature, it was well known that photoinhibition of photosynthesis was enhanced [92] and the primary target of photoinhibition was PSII [93] and phosphorylation of PSII centers increased the stability of PSII complexes and concomitantly improved their survival under stress conditions [94]. In this study, as an important component of PSII, D1 (id 3541) was rapidly phosphorylated under chilling treatment for 6 h, which might result from an imbalance between energy supply and utilization in chloroplasts of chilling-sensitive J. curcas seedling. The rate of PSII repair was reduced by inhibiting D1 synthesis at the translation-elongation stage if the repair was not efficiently scavenged [95]. Coincidently, with the time of chilling treatment prolonged, the phosphorylation level of D1 dramatically decreased. It might suggest the phosphorylation of D1 was closely related to stability of PSII complexes. During recovery from chilling-induced photoinhibition in leaves of J. curcas seedling, we observed that the CP29 (id 3634) dephosphorylation matched very well with those of PSII recovery. These results might suggest that PSII reactivation from low temperature photoinhibition was closely related to CP29 dephosphorylation. Similar result had been illuminated at O. sativa during recovery from chilling induced photoinhibition [96].
Chloroplast movement was an efficient strategy to alleviate PSII damage and maintained maximal photosynthetic output [97]. In plants, temperature was an important factor in modulating chloroplast relocations and PHOTs were considered to mediate not only stomatal opening but also chloroplast movement [98]. The protein level of PHOT1 decreased while the PHOT2 slightly increased at low temperature in A. thaliana [99]. In this study, the phosphorylation level of PHOT1 (id 3118) was down-regulated while the PHOT2 (id 3490) was up-regulated after chilling treatment for 6 h and 24 h, which possibly indicated that phosphorylation related to the activity of PHOT1 and PHOT2. However, both of the PHOT1 and PHOT2 maintained the phosphorylation level previously even after the stimulus was removed (Table 1). Besides the PHOTs, the chloroplast movement was also associated with several proteins [100]. CHUP1 might function at the periphery of the chloroplast outer membrane and most likely represented an essential component for chloroplast movement. Similarly, the mutant pmi1 exhibited severely attenuated chloroplast movements under low and high light intensities, indicating that PMI1 was necessary for both chloroplast accumulation and avoidance movements [101]. Nevertheless, we observed little change of CHUP1 (id 1231 and id 1233), PMI1 (id 2153 and id 2154) at the phosphorylation level after chilling treatment in J. curcas seedling, subsequently, all of them inducible regulated at the phosphorylation level in the recovery stage. To date, there were relatively few physiological and molecular evidences for explaining chloroplast movement through phosphorylation. Therefore, we speculated that the activities of CHUP1 and PMI1 might be regulated by other post-translational modification, such as glycosylation or sumoylation.

Plant Material and Chilling Treatment
The uniform mature J. curcas seeds with same provenance were collected from Xishuangbanna Tropical Botanical Garden, Chinese Academy of Sciences, Yunnan Province. The mature J. curcas seeds were sowed and grown in a mixed soil of peat and vermiculite (1:1) in a growth chamber under 28/25 • C (day/night) temperature regime, a photon flux density of 150 µmol m −2 ·s −1 throughout a 14-h photoperiod and a relative humidity of 80%. Four-leaf stage seedlings were used in this experiment. Low temperature treatments were started at 11:00 AM on the first day by setting the temperature to 4 • C, which was reached about 50 min later. The temperature was returned to 28 • C after 24 h of chilling treatment. The fourth true leaves were detached from the top of J. curcas seedling at each time point (4 • C for 0, 6, 24 h and then returned to 28 • C for 24 h), which indicated as C0 h, C6 h, C24 h and R24 h respectively for each biological replicate. one part of the collected leaves were used for tissue preparation of transmission electron microscopy and other part of the collected leaves from the five plants were frozen in liquid nitrogen and stored at −80 • C until use. However, in consideration of the functional maturity, the first true leaves were chosen for measurements of photosynthesis, stomatal conductance and transpiration. The measurement was carried out by using a portable gas analysis system, LI-COR 6400 with a light-emitting diode light source (LICOR Inc., Lincoln, NE, USA). The measurement conditions were set as follows: block temperature 20 • C, photon flux density 1200 µmol m −2 ·s −1 , humidity 80% and the CO 2 concentration was stabilized by 5 L surge flask.

Tissue Preparation for Transmission Electron Microscopy and Protein Preparation
The fourth true leaves from C0 h, C6 h, C24 h and R24 h were fixed in 2.5% glutaraldehyde in 100 mM phosphate buffer (pH 7.0) for 4 h at room temperature. Then the samples were handled as the previous operation [38]. Leaf total protein was extracted according to the method of Shen et al. [102] with minor modifications. Approximately 500 mg leaves of each sample was homogenized in 2 mL of the homogenization buffer containing 20 mM Tris-HCl (pH 7.5), 250 mM sucrose, 10 mM ethylene glycol tetra acetic acid, 1mM phenylmethylsulfonyl fluoride, 1% dithiothreitol (DTT) and 1% Triton X-100 on the ice. The homogenate was collected into an Eppendorf tube and centrifuged at 10,000× g for 10 min at 4 • C. The supernatant was transferred to a fresh tube and precipitated by adding 10% cold trichloroacetic acid on ice for over 30 min. The mixture was centrifuged at 15,000 g for 10 min at 4 • C and the supernatant was discarded. After washed three times with acetone containing 1% DTT, the pellet was collected by centrifugation, air-dried and then suspended in SDT buffer containing 4% SDS, 1 mM DTT, 100 mM Tris-HCl (pH 7.5), sonicated for 2 min (with cooling on ice) and boiled for 15 min. The protein mixtures were harvested via centrifugation at 15,000× g for 20 min at 4 • C to remove insoluble material. Protein concentrations of experimental samples were quantified according to Bradford method [103]. Albumin (A5503, Sigma, St. Louis, MO, USA) was used as a standard for protein quantification. Three independent biological replicates were performed independently for each experimental sample.

Phosphopeptide Enrichment
The same amount of extracted protein mixture in each sample was directly reduced with DTT, alkylated with iodoacetamide and subsequently digested with endoproteinase Lys-C and trypsin as previously described [104]. The enrichment procedure for the phosphopeptides was performed as followed. Tryptic peptides (5 mg) were dissolved in 400 µL of loading buffer containing 65% acetonitrile (ACN)/2% trifluoroaceticacid (TFA) saturated with glutamic acid and incubated with an appropriate amount (tryptic peptide: TiO 2 = 1:1, w/w) of TiO 2 beads (GL Sciences, Tokyo, Japan) for 40 min. The mixture was centrifuged at 10, 000× g for 3 min at 4 • C and the supernatant was discarded. After washing with 800 µL wash buffer (65% ACN/0.1% TFA) for 40 min, then centrifugation as above, then the phosphopeptides were eluted twice with 800 µL elution buffer (500 mM NH 4 OH/60% ACN), centrifugation and the eluate was dried and reconstituted in 0.1% formic acid/H 2 O for MS analysis.

LC-MS/MS Analysis
The enriched phosphopeptides were separated on a self-packed C18 reverse-phase column (70 µm inner diameter, 150-mm length, Column Technology, Fremont, CA, USA) that directly connected to the nanospray ion source to an Q Exactive mass spectrometer (Thermo Fisher Scientific, San Jose, CA, USA) running in the positive ion mode. The pump flow was split to achieve a flow rate of 1 µL/min for sample loading and 300 nL/min for MS analysis. The mobile phases consisted of 0.1% formic acid (A) and 0.1% formic acid and 90% ACN (B). A four-step linear gradient of 2% to 5% B in 5 min, 5% to 22% B in 55 min, 22% to 90% B in 5 min and keeping 90% B for 10 min. The spray voltage was set at 2.0 kV and the temperature of the heated capillary was 270 • C. For data acquisition, each MS scan was acquired at a resolution of 60,000 (at 400 m/z) with the lock mass option enabled. A lock mass function was used to obtain high mass accuracy. The 12 most intense precursor ions were selected for collision-induced fragmentation in the linear ion trap at normalized collision energy of 37%. The threshold for precursor ion selection was 500 and the mass window for precursor ion selection was 2.0 Da. The dynamic exclusion duration was 120 s, the repeat count was 1 and the repeat duration was 30 s. Three biological replicates were performed independently from sample collection to the phosphopeptide identification using LC-MS/MS.

Protein Identification
The raw files were processed with MaxQuant (version 1.3.0.5) and searched against the NCBI Jatropha protein database (update to 20180509, 41287 entries) concatenated with a decoy consisting of reversed sequences. The following parameters were used for database searches: trypsin/P was chosen as enzyme specificity, Carbamidomethyl (C) was selected as a fixed modification, Oxidation (M), Acetyl (Protein N-term) and Phospho (STY) were selected as variable modifications. Up to two missing cleavage points were allowed. The precursor ion mass tolerances were 7 ppm and the fragment ion mass tolerances for the MS/MS spectra were 20 ppm. The false discovery rate (FDR) was set to ≤0.01 for peptide and protein. The minimum peptide length was set to 6. The localization of the phosphorylation site was based on the PTM scores that assigned probabilities for each of the possible sites according to their site-determining ions. The software MaxQuant was used to calculate the PTM scores and PTM localization probabilities. Only when the phosphorylation site of localization probabilities (p ≥ 0.75) and the PTM score (≥5) was considered as the reliable. An FDR of 0.01 was used for phosphorylation site identification.

Label Free Quantification and Screening of Phosphopeptides with Significant Changes at the Phosphorylation Level
We integrated the ion intensities over its chromatographic elute profile and employed MaxQuant software to calculate the quantification of phosphopeptides on the basis of a label free approach. For each phosphopeptide, its intensity was normalized to the mean of the intensities of all phosphopeptides within each biological replicate. The screening terms of phosphopeptides with significant changed at phosphorylation level were listed as below: (1) phosphopeptides detected in all three biological replicates; (2) phosphopeptides with credible ANOVA analysis or Student's-test analysis (FDR<0.05); (3) phosphorylation localization probability ≥ 0.75 and (4) phosphorylation site score difference ≥ 5.

Bioinformatics
Protein function was annotated through Blast2GO software (http://www.blast2go.com/ b2ghome), significantly enriched phosphorylation motifs were extracted from phosphopeptides with confidently identified phosphorylation residues using the Motif-X algorithm (http://motif-x.med. harvard.edu/). The phosphopeptides were centered at the phosphorylated amino acid residues and aligned and seven positions upstream and downstream of the phosphorylation sites were included. Because the upload restriction of Motif-X is 10 Mb, a FASTA format data set (nearly 10 Mb) containing the protein sequences from the J. curcas protein database was used as the background database to normalize the scores against the random distributions of amino acids. The occurrences threshold was set to 5% of the input data, set at a minimum of 20 peptides and the probability threshold was set to p < 10 −6 . The phosphoproteins blasted by the National Center for Biotechnology Information (NCBI) were used to obtain the KOG numbers of those proteins by eggNOG (http://eggnogdb.embl.de/). A data set containing all the KOG numbers was then used for PPI by using the Search Tool for the Retrieval of Interaction Genes/Proteins (STRING) database (http://string-db.org/) and the PPI network was displayed by Cytoscape software (Version 3.6.3).

Conclusions
In summary, the woody oil plant species, J. curcas seedling was sensitive to chilling stress. It was studied by employing a comparative phosphoproteomic analysis, physiological measurement, ultrastructure observation under chilling stress and recovery. For the 805 significantly changed phosphopeptides, 9 phosphorylation motifs were extracted, which were mainly regulated by CDK, SnRK2, MAPK and CDPK. A complex PPI subnetwork were constructed and indicated the crucial roles of SnRK, 14-3-3 and ADP-ribosylation factor in the responsive network. Our results showed the phosphorylation was an essential event and responsible for chilling response and defense in J. curcas seedling. Consequently, Ca 2+ , ABA, ethylene, phosphoinositide and 14-3-3 mediated signal pathways cross linked for chilling response. In response to chilling stress, the phosphorylation of SnRK2a, ABI5 and TRAB1 were possibly related to ABA mediated signal pathway. Additionally, the phosphorylation level of transport, photoinhibition and chloroplast movement related proteins were also significantly regulated. We also highlighted that the phosphorylation of JcHOS1, JcKEG, JcAPX, JcPIP2 and JcPI4K might activate the potential functions to defense against chilling stress. Finally, we depicted a schematic presentation including signal transduction, metabolism and ion transport on the basis of changes at phosphorylation level, which is valuable for us to understand the chilling response and defense network in J. curcas seedling.

Conflicts of Interest:
The authors declare no conflict of interest.

GO
Gene Ontology PS Photosystem ABA Abscisic acid SnRK2 Sucrose non-fermenting 1-related protein kinase 2 ABI5 ABA Insensitive 5 HOS1 High expression of osmotically responsive gene 1 CBF C-repeat binding factor ICE1 Inducer