Universal Plant Phosphoproteomics Workflow and Its Application to Tomato Signaling in Response to Cold Stress*

Phosphorylation-mediated signaling transduction plays a crucial role in the regulation of plant defense mechanisms against environmental stresses. To address the high complexity and dynamic range of plant proteomes and phosphoproteomes, we present a universal sample preparation procedure that facilitates plant phosphoproteomic profiling. This advanced workflow significantly improves phosphopeptide identifications, enabling deep insight into plant phosphoproteomes. We then applied the workflow to study the phosphorylation events involved in tomato cold tolerance mechanisms. Phosphoproteomic changes of two tomato species (N135 Green Gage and Atacames) with distinct cold tolerance phenotypes were profiled under cold stress. In total, we identified more than 30,000 unique phosphopeptides from tomato leaves, representing about 5500 phosphoproteins, thereby creating the largest tomato phosphoproteomic resource to date. The data, along with the validation through in vitro kinase reactions, allowed us to identify kinases involved in cold tolerant signaling and discover distinctive kinase-substrate events in two tomato species in response to a cold environment. The activation of SnRK2s and their direct substrates may assist N135 Green Gage tomatoes in surviving long-term cold stress. Taken together, the streamlined approach and the resulting deep phosphoproteomic analyses revealed a global view of tomato cold-induced signaling mechanisms.

Phosphorylation-mediated signaling transduction plays a crucial role in the regulation of plant defense mechanisms against environmental stresses. To address the high complexity and dynamic range of plant proteomes and phosphoproteomes, we present a universal sample preparation procedure that facilitates plant phosphoproteomic profiling. This advanced workflow significantly improves phosphopeptide identifications, enabling deep insight into plant phosphoproteomes. We then applied the workflow to study the phosphorylation events involved in tomato cold tolerance mechanisms. Phosphoproteomic changes of two tomato species (N135 Green Gage and Atacames) with distinct cold tolerance phenotypes were profiled under cold stress. In total, we identified more than 30,000 unique phosphopeptides from tomato leaves, representing about 5500 phosphoproteins, thereby creating the largest tomato phosphoproteomic resource to date. The data, along with the validation through in vitro kinase reactions, allowed us to identify kinases involved in cold tolerant signaling and discover distinctive kinase-substrate events in two tomato species in response to a cold environment. The activation of SnRK2s and their direct substrates may assist N135 Green Gage tomatoes in surviving long-term cold stress. Taken  Protein phosphorylation is crucial for plant cells in perceiving and responding to environmental stimuli through transduction of signals from receptor kinases to targets (1,2). Compared with those of humans, plant kinases are double in number and diversity, highlighting the importance of the plant kinome and phosphoproteome in regulating responses to both abiotic and biotic stresses (3,4). Therefore, profiling the phosphoproteomic changes in response to environmental stresses is an efficient way to understand and to delineate a global view of plant defense mechanisms. Cold stress is a major environmental factor that affects the growth, distribution, and yield of many important crops growing in tropical or subtropical areas (5)(6)(7). Under prolonged cold temperatures, plant cells alter the expression of thousands of genes to reach a cold acclimation status. Many important transcription factors are expressed under cold stress to regulate plant cold acclimation (5,8). For example, the increase of transcription factor ICE1 expression plays a role in the modulation of coldresponsive genes (CORs) such as the expression of three CBF genes in Arabidopsis (9). Under cold stress, ice1 mutant plants had reduced expression of CBF genes and displayed the cold sensitive phenotype (9). Besides the important role of transcriptional factors in cold acclimation, many plant kinases are activated and positively regulate plant freezing tolerance at the post-translational level. One of the canonical events is the activation of MEKK1-MKK2-MAPK4/6 cascades in Arabidopsis under short periods of cold treatment, which has been linked to enhanced freezing tolerance (10,11). For example, a phosphoproteomics study revealed that the phosphorylation level of Thr31 on MKK2 was highly increased in the coldtolerant banana (Musa app. Dajiao) but not in the cold-sensitive Cavendish banana (12). Another example is the serine/ threonine kinase OST1, one of the core components of the Abscisic acid (ABA) 1 pathway, which modulates ICE1 protein turnover through phosphorylation at Ser278 (13). Phosphorylation of this site prevents the degradation of ICE1 protein under cold stress, which promotes cold tolerance in Arabidopsis. These examples suggest that the plant phosphoproteome and kinome are involved in the regulation of molecular events that trigger cold acclimation.
As the tomato is one of the most important horticultural crops in the world, we utilized our proteomic approach to investigate the underlying mechanisms of their cold tolerance. Considering the limited knowledge of the molecular mechanisms that regulate tomato cold tolerance, different tomato subtypes with distinct cold tolerances are suitable as a model system to study these mechanisms. Previous reports have systematically compared the cold tolerances of tomatoes in the view of the transcriptome (14,15). Distinct gene expression patterns were observed from different cold tolerant tomato variants, indicating the complexity of cold-induced molecular mechanisms in tomatoes. However, there has been no large-scale study to characterize the roles of the tomato kinome and phosphoproteome in the regulation of cold tolerance. We recently observed that a novel cultivated tomato, N135 Green Gage (cultivar, Solanum lycopersicum), is tolerant of prolonged cold exposure (4°C), whereas a wild tomato, Atacames (PIM, Solanum. pimpinellifolium), displays a coldsensitive phenotype under cold conditions. Thus, these tomatoes are ideal materials to study the important cold tolerance signaling pathways in tomato.
Mass spectrometry (MS) has emerged as a powerful technology for identifying thousands of phosphorylation sites in a single shot from mammalian cells and extracellular vesicles (16 -19). However, significant analytical difficulty is encountered in plant phosphoproteomics because of the high dynamic range of the plant proteome, the rigidity of plant cell walls, and the interference from chlorophyll and secondary metabolites (20 -22). These challenges hamper the sensitivity and efficiency of detecting low abundance phosphorylation events in plants through MS. In addition, resources concerning the tomato phosphoproteome are still limited by the lack of a full genomic sequence of the tomato. Stulemeijer et al. only reported 50 phosphoproteins from tomato seedlings using TiO 2 enrichment and LC-MS/MS analysis, which was partially attributed to poor spectrum and protein sequence matching (23). To address the limits of global plant phosphoproteomics, we have carefully evaluated the performance of several techniques that have been previously employed in each proteomic sample preparation step. We introduce here a universal sample preparation protocol, which significantly increased the coverage and depth of the plant phosphoproteome. This protocol was then applied to study the phosphoproteomic perturbation of two tomato varieties under prolonged cold treatment. This in-depth phosphoproteomic resource reveals the phosphorylation sites implicated in kinase activation and cold-responsive gene expression. Upon coupling this data with in vitro kinase screening, we discovered a connection between SnRK2s activation and cold tolerance through phosphorylation of their downstream kinases, which sheds light upon which tomato phosphoproteins are critical for conferring cold tolerance.

EXPERIMENTAL PROCEDURES
Experimental Design and Statistical Rationale-We designed and compared six sample preparation workflows detailed in Table I. Three technical replicates were performed for each protocol, and one LC-MS/MS injection was performed for each replicate resulting in a total of 18 analyses. The Tris-HCl protocol was used as a standard control to compare the number of identified phosphopeptides with the other strategies. The number of unique phosphopeptides identified from each replicate of the strategies was calculated using Microsoft Excel after removal of the redundant phosphorylated sequences in the MaxQuant output evidence file.
To access the tomato phosphoproteome and to quantify the phosphoproteomic changes in terms of long-term cold stress (4°C, 5 days), three biological replicates were collected from the 30-day-old leaves of N135 Green Gage and Atacames tomatoes treated maintained at either room temperature or 4°C (cold stress) for 5 days, resulting in a total of 12 samples. Stable isotope dimethyl labeling was performed to quantify the phosphorylation perturbation upon cold stress. The leaves from plants treated at room temperature were set as controls and the corresponding peptides were labeled using light formaldehyde (CH 2 O), whereas the samples exposed to cold stress were set as treatments and labeled using heavy formaldehyde (C 13 D 2 O). The same amount of a control and a treatment replicate were pooled together, and the labeled phosphopeptides were purified and subsequently separated into six fractions. Each fraction was analyzed by one LC-MS/MS injection, resulting in a total of 36 raw files. The number of identified phosphopeptides was calculated as described above. The statistical analysis is described in the Data Analysis subsection.
Plant Materials and Cold Treatment-Tomatoes were grown in a greenhouse at Purdue University. Tomato seeds were first germinated in soil, and the seedlings with identical growth were transferred to separate pots. Four-week-old tomato plants (about 30 cm height) were subjected to mock or cold treatment (4°C) for 5 days. Equal amounts of leaf tissue from cold tolerant and sensitive tomato varieties were collected for further electrolytic leakage and phosphoproteomic analyses.
Electrolytic Leakage Measurement-The electrolyte leakage (EL) was measured as follows: at least three fully expanded leaves from each tomato variety were detached and immersed in 50 ml tubes with 25 ml distilled water. After gentle shaking overnight, the initial electrolyte conductivity (E1) was measured with a conductivity meter. Next, the samples were autoclaved. After cooling to room temperature, the conductivity (E2) was measured again. The relative electrolyte leakage was calculated as: E1/E2 ϫ 100.
Real-time PCR-RNA was extracted from 1-month-old tomato plants grown in soil under mock treatment or 4°C for the indicated time. 1 g RNA was used for reverse transcription using M-MLV reverse transcriptase (Promega). Real-time PCR was performed using iQ SYBR Green Supermix (Bio-Rad, Hercules, CA) on a CFX96 realtime PCR detection system (Bio-Rad). Tomato Actin 7 (ACT7) was used as the internal reference for all reactions.
Filter Aided Sample Preparation-A 200 g aliquot of lysate was loaded into a Microcon 10k Da (MilliporeSigma, Burlington, MA) filter tube using the FASP protocol (26). The GdnHCl buffer was replaced by 100 l of 8 M urea buffer; The urea buffer was replaced by 100 l 50 mM ammonium bicarbonate (ABC). Lys-C was added onto the filter at a 1:100 (w/w) ratio for 3 h at 37°C, and trypsin was subsequently added to a final 1:100 (w/w) ratio overnight. The digests were collected with two washes of 100 l 50 mM ABC buffer, and then the peptides were acidified with 10% trifluoroacetic acid (TFA) to pH ϳ3 and desalted using a 50 mg Sep-Pak C18 column.
Methanol-Chloroform Precipitation-The protein precipitation was performed as previously described with some modifications (27). First, 200 g of lysate was added to four volumes of methanol, followed by an equal volume of chloroform with mixing. Three volumes of ddH 2 O were added to the tube with mixing. The solution was centrifuged at 16,000 ϫ g for 3 min. The upper aqueous layer was removed, the protein pellet was washed with four volumes of methanol, and the tube was centrifuged again. The supernatant was discarded, and the precipitated protein pellet was air dried.
Protein Digestion-Protein extracts (200 g) were 5-fold diluted for the SDC-SLS protocol (28), 8-fold diluted for the urea protocol, or 10-fold diluted for the GdnHCl protocol using 50 mM TEAB. Proteins were then digested with Lys-C (FUJIFILM Wako Chemicals, Richmond, VA) in a 1:100 (w/w) enzyme-to-protein ratio for 3 h at 37°C, and trypsin (Sigma-Aldrich) was added to a final 1:100 (w/w) enzymeto-protein ratio overnight. The SDC and SLS were separated from digested peptides by acidifying the solution using 10% TFA, adding 100%ethyl acetate (50% final volume), and removing the organic solution. The aqueous layer was dried in a SpeedVac, and digests were then suspended with 0.1% TFA and desalted using a 50 mg Sep-Pak C18 column (Waters, Milford, MA).
Dimethyl Labeling-Dimethyl-labeling was performed as previously described (29). The tryptic peptides were dissolved in 100 l of 100 mM TEAB and mixed with 4 l of 4% 13 CD 2 O or 12 CH 2 O, and then 4 l of freshly prepared 0.6 M sodium cyanoborohydride was immediately added. The mixture was agitated for 60 min at room temperature. The reaction was stopped by adding 16 l of 1% ammonium hydroxide on ice and agitating the mixture for 1 min. The dimethyl labeled peptides were then mixed with 20 l of 10% formic acid (FA), and the two labeled samples were combined and desalted using a Sep-Pak C18 column.
PolyMAC Enrichment-Phosphopeptides were enriched using a PolyMAC-Ti kit (Tymora Analytical, West Lafayette, IN) (30). Briefly, the digested peptides were resuspended with 20 l of ultrapure water, 200 l of Loading buffer (50 mM glycolic acid in 0.5% TFA and 95% ACN) was added, and the sample was vortexed. Next, 50 l of the PolyMAC/Magnetic Capture beads were added to the sample and incubated for 20 min. The solvent was removed using a magnetic separator rack, and the beads were washed with 200 l of Washing buffer 1 (1% TFA in 80% ACN) for 5 min. The beads were then incubated with 200 l of Washing buffer 2 (80% ACN in ddH 2 O) for 5 min, and the solvent was removed as before. Finally, the beads were incubated twice with 100 l of Elution buffer (400 mM NH 4 OH in 50% ACN) for 5 min each. The eluates were combined and completely dried in a SpeedVac.
Basic pH Reverse-phase Fractionation-Basic pH reverse-phase fractionation was performed with some modifications as previously described (31). First, 2 mg of Magic C18-AQ beads (5 m particle size) were suspended in 100 l of methanol and loaded into a 200 l StageTip with a 20 m polypropylene frit. The C18 StageTips were activated with 100 l of 40 mM NH 4 HCO 2 , pH 10, in 80% ACN and equilibrated with 100 l of 200 mM NH 4 HCO 2 , pH 10. The isolated phosphopeptides were suspended with 200 mM NH 4 HCO 2 , and the C18 StageTips were washed with 100 l of 200 mM NH 4 HCO 2 , pH 10. The bound phosphopeptides were fractionated from the StageTip with 50 l of 8 buffers containing different ACN concentrations: 5%, 9%, 13%, 17%, 21%, and 80% of ACN in 200 mM NH 4 HCO 2 , pH 10. The eluted phosphopeptides were dried and stored at Ϫ20°C.
Tomato SnRK2E Protein Expression-Tomato SnRK2E coding sequence was amplified from tomato cDNA with forward primer 5Ј-CGGAATTCATGGATCGGACGGCAGTGACA 3Ј and reverse primer 5Ј-GCGTCGACTTACATTGCATAGACAATCTC. The PCR product was digested with EcoRI and SalI and inserted into a pGEX4T-1 vector. The construct was then introduced into BL21 cells and expressed in E. coli. Recombinant GST-SnRK2E protein was purified using glutathione-agarose beads (GST) (Sigma-Aldrich) according to the manufacturer's instructions.
In vitro Kinase Reaction-The in vitro kinase reaction was performed based on the siKALIP approach with some modifications (32). First, 200 g of Lys-C digested peptides were dephosphorylated using TSAP (Roche, Basel, Switzerland) in a 1:100 (w/w) enzyme-topeptide ratio at 37°C overnight. The reaction was terminated by heating the sample at 75°C for 10 min. The dephosphorylated peptides were desalted using a Sep-Pak C18 column and then suspended in kinase reaction buffer (50 mM Tris-HCl, 10 mM MgCl 2 , 1 mM DTT, and 1 mM ␥-[ 18 O 4 ]-ATP, pH 7.5). The SnRK2E kinase (1 g) was incubated with the desalted peptides at 25°C overnight. The kinase reaction was quenched by acidifying with 10% TFA to a final concentration of 1%, and the peptides were again desalted with a Sep-Pak C18 column. The heavy 18 O-phosphopeptides were further digested by trypsin at 37°C for 6 h and enriched by PolyMAC-Ti as described above. The eluates were dried in a SpeedVac prior to LC-MS/MS analysis.
LC-MS/MS Analysis-The peptides were dissolved in 5 l of 0.3% FA with 3% ACN and injected into an Easy-nLC 1000 (Thermo Fisher Scientific). Peptides were separated on a 45 cm in-house packed column (360 m OD ϫ 75 m ID) containing C18 resin (2.2 m, 100Å, (Bischoff Chromatography, Leonberg, Germany) with a column heater (Analytical Sales and Services, Flanders, New Jersey) set at 50°C. The mobile phase buffer consisted of 0.1% FA in ultra-pure water (buffer A) with an eluting buffer of 0.1% FA in 80% ACN (buffer B) run over a linear 60 min (method comparisons) or 90 min (large-scale phosphoproteomics) gradient of 5-30% buffer B at a flow rate of 250 nL/min. The Easy-nLC 1000 was coupled online with a LTQ-Orbitrap Velos Pro mass spectrometer (Thermo Fisher Scientific). The mass spectrometer was operated in the data-dependent mode in which a full MS scan (from m/z 350 -1500 with the resolution of 30,000 at m/z 400) was followed by the 10 most intense ions being subjected to collision-induced dissociation (CID) fragmentation. CID fragmentation was performed and acquired in the linear ion trap (normalized collision energy (NCE) 30%, AGC 3e4, max injection time 100 ms, isolation window 3 m/z, and dynamic exclusion 60 s).
Data Processing-The raw files were searched directly against a Solanum lycopersicum database from the Sol Genomics Network (http://solgenomics.net) version 3.0 (33) (30,868 entries) with no redundant entries using MaxQuant software (version 1.5.5.1) (34) with a 1% FDR cutoff at the protein, peptide, and modification levels. The first peptide precursor mass tolerance was set at 20 ppm, and MS/MS tolerance was set at 0.6 Da. Search criteria included a static carbamidomethylation of cysteines and variable modifications of (1) oxidation on methionine residues, (2) acetylation at the N terminus of proteins, and (3) phosphorylation on serine, threonine or tyrosine residues. For 18 O-phosphopeptides, heavy phosphorylation (ϩ85.979 Da) was also set as a variable modification. Searches were performed with full tryptic digestion and allowed a maximum of two missed cleavages on the peptides analyzed from the sequence database. Dimethyl-labeling quantitation was performed using MaxQuant by setting the multiplicity as 2, and the match between runs function was enabled with a match time window of 1.0 min. The false discovery rates of proteins, peptides and phosphosites were set at 1% FDR. The minimum peptide length was six amino acids, and a minimum Andromeda score was set at 40 for modified peptides. A site localization probability of 0.75 was used as the cut-off for localization of phosphorylation sites.
Data Analysis-All data were analyzed using the Perseus software (version 1.5.6.0) (35). The intensities of both light and heavy dimethyllabeled phosphorylation sites were extracted from MaxQuant output file. The unlocalized phosphorylation sites (localization probability Ͻ 0.75) were filtered out through Perseus. The intensities of class I phosphorylation sites (localization probability Ͼ 0.75) were log2 transformed, and the quantifiable phosphorylation sites were selected from the identification of all triplicate replicates in at least one sample group. The missing intensity values were replaced from a normal distribution (width ϭ 0.2 and down shift ϭ 1.8) for statistical analysis. The significantly enriched phosphorylation sites were selected by the ANOVA test with a permutation-based FDR cut-off of 0.05. Normalization was carried out by subtracting the median of the log2 transformed phosphopeptide intensities. Principal component analysis (PCA) of the tomato samples was performed with a Benjamin-Hochberg FDR cutoff of 0.05. For hierarchical clustering, the intensities of the ANOVA significant phosphorylation sites (FDR Ͻ 0.05) were first z-scored and clustered using Euclidean as a distance measure for column and row clustering. The number of clusters was set at 300, with a maximum of 10 iterations and 1 restart.
All the localized and significantly changed phosphorylation sites were submitted to Motif-X (36) to determine the enriched phosphorylation motifs with the ITAG 3.0 database as background. The significance was set at 0.000001, the width was set at 13, and the minimum number of occurrences was set at 20. The phosphorylation clusters in Fig. 4C were visualized using GProX version 1.1.13 (37) using fuzzification value of 2 and 100 iterations.

A Universal Sample Preparation Workflow for Plant Phos-
phoproteomics-An optimized sample preparation protocol for plant samples is critical to enable the identification of low abundant phosphoproteins in plant cells and to deepen the coverage of the plant phosphoproteome. To improve the depth of the plant phosphoproteome, it is necessary to evaluate the key factors during sample preparation that affect phosphopeptide identification, such as the efficiency of plant protein extraction, removal of interferences, and the efficiency of enzymatic digestion (Fig. 1A). One of the major challenges is that the toughness of the plant cell wall reduces the efficiency of plant protein extraction compared with mammalian cells. To address this problem, various groups have reported the use of Tris-HCl with salts (24), Tris-HCl with 30% sucrose (38), or SDS coupled with Filter Aided Sample Prep (FASP) (39) for extraction of proteins from plant tissues. However, these results showed varying identification efficiencies and yields, meaning that the lysis conditions for mature plant tissues are still not optimized. Because guanidine hydrochloride (GdnHCl) has been used as a protein denaturing reagent in plant RNA extraction protocols with great success (40), we hypothesized that GdnHCl is an efficient lysis reagent that can be used to break down the cell wall and to extract the proteins in the cells. To compare the efficiency of lysis strategies for mature plant tissues, we tested three lysis protocols: (1) Tris-HCl protocol, (2) SDC-SLS protocol, and (3) GdnHCl protocol (see Table I Table S1-S3). By combining triplicate one-shot MS analyses, the GdnHCl protocol permitted the identification of Ͼ4000 phosphopeptides more than the other protocols in quick LC-MS analyses (Fig. 1C). We observed that GdnHCl lysis buffer generated the deepest green color of lysate compared with the other protocols, suggesting that the GdnHCl buffer may provide the best lysis efficiency for breaking down the cell wall and organelle membranes (supplemental Fig. S1A). Therefore, we evaluated the protein yields using the three lysis buffers and further performed gene ontology (GO) analysis to determine the significantly enriched cellular component categories from the three protocols using Fisher's exact test (FDR Ͻ 0.05). The results show the GdnHCl buffer extracted the highest protein amount of the three protocols (supplemental Fig. S1B). Moreover, the cytosolic and organelle phosphoproteins are highly enriched in the GdnHCl protocol compared with those using the Tris-HCl and SDC-SLS protocols, indicating that using the GdnHCl protocol provides the best phosphoprotein extraction yield of the three protocols (supplemental Fig. S1C).
Another crucial factor that affects phosphopeptide enrichment is presence of nucleic acids, phospholipids, photosynthetic pigments, and secondary metabolites in plant tissues. These interferences increase the sample complexity, resulting in low phosphopeptide purification efficiency and specificity. It is difficult to remove these interferences from digested peptides using C18 cartridges because most of them have hydrophobicity values like those of digested peptides, suggesting that it is necessary to eliminate them prior to enzymatic digestion. Nucleic acids may also interfere with phosphopeptide enrichment and identification (41). Thus, we compared the performance of FASP and methanol-chloroform precipitation to remove the aforementioned unwanted molecules before enzymatic digestion (Table I). We found that  using 10k Da MWCO filters to remove the interferences does not improve the identification numbers ( Fig. 1D and supplemental Table S4), which may be attributed to high sample loss caused by the multiple washing steps in the FASP protocol. On the other hand, performing methanol-chloroform precipitation further increases the number of identified phosphopeptides by 20% (5141 phosphopeptides) compared with the GdnHCl direct digestion protocol without the precipitation step (4242 phosphopeptides) (Fig. 1D and supplemental Table S5), demonstrating that it is beneficial to eliminate the interferences in plant tissues prior to the phosphopeptide enrichment.
The next step that influences phosphopeptide identification is the efficiency of enzymatic digestion. Masuda et al. reported that the phase-transfer surfactant (PTS) protocol significantly improves the solubility of hydrophobic proteins and the activities of proteases compared with a 8 M urea protocol using a small amount of cell lysate (25,28). Thus, we compared two digestion protocols: first, the urea digestion protocol and second, the PTS digestion protocol ( Table I). The lysate solution was divided into two parts for protein precipitation. Each precipitated protein pellet was suspended in one of the two buffers with several pulses of sonication and diluted for subsequent enzymatic digestions. The identification results show that using SDC and SLS as surfactants improves the number of identified phosphopeptides (5446 phosphopeptides) compared with urea digestion protocol (5141 phosphopeptides) (Fig. 1E). In addition, the PTS protocol produces a much lower percentage of missed cleavages (17% versus 27%), meaning that the protease enzymatic activities are enhanced in the SDC-SLS buffer compared with the urea digestion buffer (Fig. 1F).
Overall, the new universal protocol with the optimized procedures in protein extraction, precipitation, and tryptic digestion (supplemental Fig. S2A) allowed us to enlarge the coverage of tomato phosphoproteome more than 2-fold compared with the standard protocol, with the identification of ϳ8100 unique phosphopeptides from triplicate single-shot 1-hr LC-MS analyses of mature tomato leaves (supplemental Fig.  S2B and supplemental Table S6).
Two Cold-responsive Phenotypes of Tomatoes-We applied the universal sample preparation procedure to study the tomato global phosphoproteomic changes in response to cold stress and the underrepresented signaling mechanisms involved in tomato cold tolerance. We compared two tomato varieties: Atacames (PIM, Solanum pimpinellifolium), a wild species that displays a cold-sensitive phenotype, and N135 Green Gage (Solanum lycopersicum), a cultivated tomato with a cold-tolerant phenotype. To compare their cold tolerance under prolonged cold exposure, the two tomato varieties were germinated in soil for about one month and then were transferred to 4°C or room temperature for 5 days. The photographs illustrate that both two tomato varieties grew well under room temperature (supplemental Fig. S3A). How-ever, after 5 days of cold treatment (4°C), only a few leaves of N135 Green Gage (S. lycopersicum) exhibited slight wilting, whereas an aggravated wilting of leaves was observed in Atacames (S. pimpinellifolium) (supplemental Fig. S3A). We further measured the degree of electrolytic leakage (EL) from the leaves of two tomato varieties, which represents the damage to the plasma membrane and the sensitivity to cold stress. At room temperature, the two tomato varieties had similar electrolytic leakage values (ϳ 6.6%). However, the EL value of Atacames was significantly increased to ϳ60% after 5 days of 4°C treatment, which is about 3-fold higher than that of N135 Green Gage (ϳ16%) (supplemental Fig. S3B). The above results indicate that N135 Green Gage had a higher tolerance to long-term cold stress. We further evaluated the expression level of cold-responsive gene CBF3, which is an important component of cold stress resistance in plants. Consistent with the observed cold phenotypes, the tomato CBF3 gene is expressed at higher levels in N135 Green Gage than in Atacames after cold treatments (supplemental Fig. S3C). These findings demonstrate that Atacames (S. pimpinellifolium) was more sensitive to cold stress when compared with N135 Green Gage (S. lycopersicum). These data prompted us to compare the phosphoproteomic changes between the two tomato varieties to uncover the signal transductions critical for surviving prolonged cold stress.
Profiling of Tomato Phosphoproteomic Changes under Cold Stress-To understand the underlying signals involved in cold stress tolerance, we performed MS analyses to measure the phosphoproteomic changes of the leaves of N135 Green Gage and Atacames treated at 4°C for 5 days versus control plants grown at room temperature (25°C) (Fig. 2A). Each experiment was performed in biological triplicates. We used the optimized protocol for protein extraction, precipitation, and enzymatic digestion (supplemental Fig. S2A). The digested peptides (from 1 mg protein) were dimethyl-labeled, and labeled peptides were mixed. The pooled phosphopeptides were enriched using PolyMAC-Ti, and the purified phosphopeptides were fractionated into six fractions by basic pH reverse-phase StageTips and subsequently analyzed by a 2 h LC-MS/MS method on a LTQ Orbitrap Velos Pro mass spectrometer. We identified ϳ5500 phosphoproteins and Ͼ30,000 unique phosphopeptides, corresponding to ϳ23,500 unique phosphorylation sites (ϳ14,400 class I sites) (Fig. 2B, supplemental Fig. S4A, and supplemental Table S7-S8). Singly phosphorylated peptides are dominant (87%) in the identified phosphopeptides (Fig. 2C). Among all the class I phosphorylation sites, ϳ85% are localized at serine residues, and ϳ1.5% are at tyrosine residues (Fig. 2C). To our knowledge, this is the first large-scale tomato phosphoproteomic study with such extensive coverage.
To evaluate the quality of the biological replicates of each sample, we performed a principal component analysis (PCA), which shows that all biological replicates are tightly together and each sample is well-separated; this means that the sam-ples cluster by the nature of sample but not by the batch (Fig.  3A). The PCA analysis also indicated that the two tomato varieties have distinct phosphoproteomic profiles (50%, PC1), and cold stress rewired the signal transduction of the two tomato varieties (ϳ30%, PC2). To gain insight into the phosphoproteomic changes in response to cold stress in the two tomato varieties, we performed an ANOVA test controlled with FDR at 0.05 to locate the phosphorylation sites with significant change upon cold stimulation (supplemental Table S9). Unsupervised hierarchical clustering further confirmed the PCA result that two major clusters are composed of biological replicates of the two tomato species (supplemental Fig. S4B). The heat map reveals six main clusters of the most changed phosphorylation sites (Fig. 3B). Clusters 1 and 2 contain the phosphorylation sites which are relatively activated in Atacames (S. pimpinellifolium) and in Atacames under cold stress, respectively. The phosphorylation motif analysis (p Ͻ 0.000001) shows that the [-pS-P-] motif is dominant in these clusters (Fig. 3C), indicating that cold-induced kinases which recognize this motif are crucial in the signaling of Atacames. Cluster 3 represents the phosphorylation events induced un-der cold stress in both tomato varieties, and the [-R-X-X-pS-], [-pS-P-], [-L-X-X-X-X-pS], and [-G-pS-] motifs are enriched in this cluster (Fig. 3C, supplemental Fig. S4C). Clusters 4 and 5 are the phosphorylation sites up-regulated in N135 Green Gage under cold conditions. We found one acidic ([-pS-D-]) and one basic phosphorylation motif ([-R-X-X-pS-]) enriched in clusters 4 and 5 but not identified in clusters 1 and 2, indicating that these proteins might be uniquely involved in cold-tolerant mechanisms in the N135 Green Gage tomato (Fig. 3C). The sixth cluster shows the repressed phosphorylation sites in the two tomatoes in a cold environment. Many phosphorylation sites were significantly up-or downregulated in response to cold stress, implying that the activities of kinases and phosphatases were significantly perturbed.
Kinase-Mediated Molecular Events under Cold Stress-We next evaluated the perturbation of kinase activity in the two tomato species in response to cold stress to help us understand the cascade of signaling pathways that regulate cold tolerance. Considering that phosphorylation on kinases plays an important role in mediating kinase activities, we compared FIG. 2. Profiling the phosphoproteome of tomato leaves in response to cold stress. A, The tomatoes leaves were grown in three biological replicates for 5 days at room or cold temperature. Proteins were extracted, precipitated, and digested using the optimized protocol. The digested phosphopeptides were dimethyl-labeled, combined, and enriched using PolyMAC. To enlarge the phosphoproteome coverage, the isolated phosphopeptides were fractionated using basic pH reverse-phase StageTips prior to LC-MS/MS analysis. B, The total number of identified phosphoproteins, phosphopeptides, and phosphorylation sites in the phosphoproteome of tomato leaves. C, Distribution of the number of phosphates and amino acid residues for all detected phosphopeptides. the phosphorylation changes of 669 phosphorylation sites on 292 identified kinases in two cold-treated tomato varieties (Fig. 3D). Interestingly, we found that the two tomato varieties have unique phosphorylation expression patterns in terms of cold stress perturbation. For example, in Atacames (S. pimpinellifolium), many phosphorylation sites in the activation loops of mitogen-activated protein kinases (MAPK) are highly activated. This includes MAPK 1/2 (pT219 and pY221), which is the homolog of MAPK3/6 in Arabidopsis (pT221 and pY223); MAPK 9 (pT191 and pY193), which is the homolog of AtMAPK1 (pT191 and pY193); MKK 1 (pT31), which is the homolog of AtMKK2 (pT31) (supplemental table S5A). The activation of MAPKs in cold-treated Atacames tomato could explain why the [-pS-P-] motif is more enriched in the cold- induced phosphorylation sites. It is therefore possible that Atacames has the same signaling mechanisms as Arabidopsis in which cold stress is known to induce phosphorylation of the activation loops of AtMKK4/5-MAPK3/6 and AtMKK1/2-MAPK4/6 to mediate cold tolerance (10). In our data, we also identified two sites (pS436 and pS440) with the [-pS-P-] motif on bHLH transcription factor 87 (Solyc06g068870.3.1) that were highly phosphorylated after 5 days of cold treatment (supplemental Fig. S5B). bHLH transcription factor 87 is the homolog protein of AtICE1 in Arabidopsis, and this transcription factor is one of the substrates of MAPK3/6 in cold-treated Arabidopsis (42), indicating the importance of MAPKs signaling in Atacames under cold stress.
To identify the kinases responsible for the enriched acidic and basic phosphorylation motifs in cold-treated N135 Green Gage (clusters 3 and 4), we searched the up-regulated phosphorylation sites (Ͼ 2-fold increase) on the kinases in N135 Green Gage in response to cold treatment ( Fig. 3D and supplemental Table S5A). We found that several phosphorylation sites on the casein kinases such as pS355 of Solyc07g040720.3.1, pS310 of Solyc09g066440.3.1, and pS347 of Solyc12g007280.2.1-the homologs of AtCKL2, AtCKL6, and AtCKL4, respectively-were significantly increased in cold-treated N135 Green Gage leaves (Fig. 3D).
Interestingly, we observed that cold stress highly induced two phosphorylation sites on two kinases (pS154 of Solyc01g103940.3.1 and pS173 of Solyc01g108280.3.1) in N135 Green Gage that are the homologs of SnRK2A and SnRK2E in Arabidopsis, respectively (Fig. 3D). SnRK2s are critical components of ABA signaling, and their activities are strictly regulated in plant cells to switch between growth and stress response statues (43). It has been reported that SnRK2s autophosphorylate in the activation loop after ABA blocks the activity of phosphatase 2C; the autophosphorylation activates their kinase function, resulting in the phosphorylation of many of their downstream substrates in response to ABA exposure and osmotic stresses in Arabidopsis (24). Furthermore, we also observed that two phosphorylation sites on previously reported SnRK2s direct substrates in Arabidopsis were only identified in N135 Green Gage. The first site is pS58 of abscisic acid responsive elements binding factor 2 (ABF2) (Solyc01g108087.1.1), which is an important transcription factor involved in activation of ABA responsive gene expression in Arabidopsis (44). The second site is pS878 of enhancer of mRNA decapping protein 4 (VCS) (Solyc07g065950.3.1), which has the [-R-X-X-pS-] phosphorylation motif; as suggested by its name, VCS is one of the components of the mRNA decapping complex (45). It has been reported that vcs mutant plants display defects in shoot apical meristem formation and a severe temperature-dependent phenotype (46). VCS is a substrate of SnRK2s in Arabidopsis, and many stress-responsive genes were misregulated in the VCS-knockdown plants (45). These findings are consistent with our motif analysis results that the [-R-X-X-pS-] sequence-the AtSnRK recognition motif-was enriched in clusters 4 and 5. In summary, SnRK2s were highly activated in N135 Green Gage under a cold environment, implying that cold stress induces the activation of SnRKs and the activated SnRK2s might be crucial in regulating cold-responsive signaling in N135 Green Gage.
Besides the high levels of phosphorylation on CKLs and SnRK2s in N135 Green Gage, we also observed that several phosphorylation sites (pS303, pS602, pS609, pS613, pS1387, and pS1390) on the serine/threonine-protein kinase ATM (Solyc04g049930.3.1) were significantly up-regulated in cold-treated N135 Green Gage compared with Atacames. Activation of ATM through phosphorylation is a key step in the cellular response to DNA damage under oxidative stresses (47). Plant cells tend to accumulate reactive oxygen species (ROS) in response to cold stress, and ROS act as signal transducers to activate downstream gene expression and kinase activity to achieve cold acclimation (48). It has been reported that the phosphorylated histone variant H2AX, one of the substrates of ATM, is a component of repairing complexes at DNA damage sites (49). We observed two phosphorylation sites (pS121 and pS124) in cluster 5 at the C-terminal of histone H2A (Solyc05g014835.1.1) were only phosphorylated in N135 Green Gage, implying that cold stress may induce DNA repair signaling in N135 Green Gage. It is possible that ATM kinases play a role in sustaining genome integrity during cold acclimation in N135 Green Gage.
Identification of SRK2E Substrates in Tomato by Stable Isotope Labeled Kinase Assay Linked Phosphoproteomics (siKALIP) Strategy-These findings prompted us to systematically identify the direct substrates of SnRK2E in N135 Green Gage to delineate the downstream signal transductions of SnRK2E under cold stress. To this end, we utilized a recently introduced siKALIP approach which couples in vitro kinase assay screening with in vivo kinase-dependent phosphoproteomics to identify direct substrates of SnRK2E (32) (Fig. 4A). For in vitro SnRK2E screening, the leaves of N135 Green Gage were lysed and digested according to the optimized protocol. Then, 200 g of digested proteins were dephosphorylated using temperature-sensitive alkaline phosphatase (TSAP), and the reaction was terminated by heating the samples at 75°C for 5 min. The dephosphorylated peptides were incubated with recombinant SnRK2E (Solyc01g108280.3.1) and ␥-18 O 4 -ATP. Finally, the heavy 18 O-phosphopeptides phosphorylated by SnRK2E were enriched through Poly-MAC-Ti and analyzed by LC-MS/MS. A total of 1,770 class I heavy 18 O-phosphorylation sites were identified from the in vitro kinase screen (32) (supplemental Table S10). We first analyzed the enriched phosphorylation motifs of SnRK2E using Motif-x algorithm. Several classic SnRK2s phosphorylation motifs including [-L-X-R-X-X-pS-] and [-R-X-X-pS-] are overrepresented in the SnRK2E phosphorylated sites (Fig.  4B), suggesting the kinase screening has high specificity. Next, we overlapped the 1398 heavy phosphorylation sites that are identified in at least two of the three replicates with the in vivo phosphorylation sites up-regulated (Ͼ 2-fold) after cold treatment in N135 Green Gage. A total of 75 phosphorylation sites corresponding to 72 phosphoproteins were selected as candidate SnRK2E substrates (Fig. 4C, supplemental Table S11). Among them, the homologs of the known SnRK2s substrates in Arabidopsis, ABF2 and VSC, were both phosphorylated by SnRK2E, validating the sensitivity and specificity of the siKALIP approach. Interestingly, we also identified multiple kinases as potential direct substrates of SnRK2E in N135 Green Gage that may be involved in downstream cold acclimation signaling. For example, we observed that SnRK2E directly phosphorylated two previously mentioned kinases, pS355 of CKL2 and pS609 of ATM ( 4D). These findings suggest that cold stress may activate cold tolerance responses through kinases such as the SnRK-CKL2 and SnRK-ATM axes in N135 Green Gage. DISCUSSION The high complexity and dynamic range of the plant proteome is a major challenge in studying plant phosphoproteomics. We achieved efficient coverage of the plant phosphoproteome by enhancing the lysis efficiency, removing the interfering species, and optimizing the digestion protocols, which permitted identification of many low abundant phosphorylation sites from both cytosol and organelles without depletion of high abundance proteins such as RuBisCo prior to PolyMAC enrichment. This optimized protocol is applicable to studies of global phosphoproteomic remodeling of any given plant species after environmental or genetic perturbation. We have demonstrated that more than 14,000 class I phosphorylation sites corresponding to 5000 phosphoproteins can be identified in two novel tomato species that display distinct tolerance phenotypes after prolonged cold stress, which provides an informative resource to the community to study cold acclimation mechanisms in tomatoes.
More importantly, this work not only provides the deepest quantitative analysis of the tomato cold-induced phosphoproteome to date but also partially reveals the underlying mechanisms conferring the cold tolerance in tomatoes at the posttranslational level. Through phosphorylation motif analysis, we identified known molecular mechanisms in cold-treated Atacames (S. pimpinellifolium) and further discovered new kinase-substrate relationships that may regulate cold tolerance in N135 Green Gage (S. lycopersicum). This analysis shows that the two tomato varieties have distinct kinome remodeling during prolonged cold stress. For example, several MAPKs are activated in Atacames in response to a cold environment, whereas in N135 Green Gage two SnRK2 kinases are activated. The data is supported by a previous report that the concentration of ABA is significantly increased in S. lycopersicum after 12 h and continues rising until 72 h under exposure to cold temperature (10°C) (50). Because ABA is the upstream activator of the kinase SnRK2E, we propose that the induction of SnRK2E activity is attributed to the high production of ABA in N135 Green Gage under prolonged cold stress.
The overrepresented acidic phosphorylation motif in coldtreated N135 Green Gage suggests that kinases recognizing the acidic motif maybe involved in cold acclimation signaling in N135 Green Gage. It has been reported that the expression of CKL2 is induced during drought stress and ABA treatment, and CKL2 participates in stomatal closure to prevent water loss by phosphorylating and physically interacting with actin depolymerizing factor 4 (ADF4) to prevent its degradation (51). Considering that phosphorylation levels on the casein kinases increase after cold stimulation in N135 Green Gage, it is possible that the casein kinases are activated in this species to alleviate water loss caused by prolonged cold exposure. Furthermore, casein kinases are known to recognize acidic phosphorylation motifs, and the motif [-pS-D-] is significantly enriched in the clusters of phosphorylation sites associated with cold-treated N135 Green Gage (clusters 4 and 5). The combination of our phosphoproteomic data with in vitro SnRK2E kinase-substrate screening enabled us to link the cold-induced acidic ([-pS-D-]) motif with the activation of SnRK2E in N135 Green Gage. Our results suggest that SnRK2E phosphorylates CKL2 under cold stress in N135 Green Gage, and CKL2 further phosphorylates many downstream substrates with an acidic motif to regulate cold tolerance. We also observed that SnRK2E directly phosphorylates several sites of ATM and histone H2A in N135 Green Gage, suggesting a potentially important role for the SRK-ATM-H2A axis in N135 Green Gage cold stress.
In summary, we have developed a novel phosphoproteomic pipeline for profiling the tomato phosphoproteome and have illustrated its application to the study of signal transduction in tomatoes under cold exposure. More importantly, we characterized the induction of different kinases in two tomato varieties which show distinct cold phenotypes. Our analyses identified novel tomato signaling cascades that regulate cold acclimation. Together, we expect that this large-scale quantitative phosphoproteomic data will serve as a useful database to aid in identification of key enzymes that trigger cold tolerance in tomatoes.