The systemin signaling cascade as derived from phosphorylation time courses under stimulation by systemin and its inactive Thr17Ala (A17) analog

Systemin is a small peptide with important functions in plant wound response signaling. Although transcriptional responses of systemin action are well described, the precise signaling cascades involved in its perception and signal transduction are poorly understood at the protein level. Here we use a phosphoproteomic profiling study involving stimulation time courses with systemin and its inactive analogon A17 to reconstruct a systemin-specific kinase/phosphatase signaling network. The time course analysis of systemin-induced phosphorylation patterns revealed early events at the plasma membrane, such dephosphorylation of H+-ATPase, rapid phosphorylation of NADPH-oxidase and Ca2+-ATPase. Later responses involved transient phosphorylation of small GTPases and vesicle trafficking proteins, as well as transcription factors. Based on a correlation analysis of systemin-specific phosphorylation profiles, we predict substrate candidates for 56 systemin specific kinases and 18 phosphatases. Among the kinases are several systemin-specific receptor kinases as well as kinases with downstream signaling functions, such as MAP-kinases. A regulatory circuit for plasma membrane H+-ATPase was predicted and confirmed by in-vitro activity assays. In this regulatory model we propose that upon systemin treatment, H+-ATPase LHA1 is rapidly de-phosphorylated at its C-terminal regulatory residue T955 by phosphatase PLL5, resulting in the alkalization of the growth medium within 2 minutes of systemin treatment. We further propose that the H+-ATPase LHA1 is re-activated by MAP-Kinase MPK2 later in the systemin response. MPK2 was identified with increased phosphorylation at its activating TEY-motif at 15 minutes of treatment and the predicted interaction with LHA1 was confirmed by in-vitro kinase assays. Our data set provides a valuable resource of proteomic events involved in the systemin signaling cascade with a focus on predictions of substrates to systemin-specific kinases and phosphatases.


Introduction
Almost 30 years ago, the quest for signaling molecules mediating systemic defense responses after local injury by insect herbivores culminated in the discovery of systemin as the first peptide with signaling function in plants (1,2). The 18-amino-acid oligopeptide was established as an essential component of the wound signaling pathway that is responsible for the systemic regulation of defense gene expression (3,4). Systemin was initially described as the long-sought hormonal signal that is released at the site of wounding, that travels through the vasculature and induces the defense response in distal, unwounded tissues (2,5). However, this model had to be modified when it was shown that systemin rather acts locally at the site of wounding, where it induces and amplifies the production of jasmonates as systemic signals for defense gene activation in distal tissues (6,7). Considering its paracrine immune-modulatory activity, systemin is thus better described as a plant cytokine, a phytocytokine (8) rather than a wound hormone. Systemin is synthesized as a larger precursor protein from which it is proteolytically released (4,9). Whether or not precursor processing and systemin secretion are regulated processes still awaits an answer. On the other hand, if systemin is released passively simply as a result of tissue disruption, it could also be addressed as a damageassociated molecular pattern, a DAMP (8).
Local action of systemin is initiated by its interaction with a saturatable binding site at the cell surface (10,11). Purification and characterization of the binding protein tentatively placed the receptor into the family of leucin-rich repeat receptor-like kinases (LRR-RLKs), but the protein initially declared as the systemin receptor later turned out to be the tomato ortholog of BRI1, the brassinosteroid receptor in Arabidopsis (12)(13)(14)(15). Recent work identified the systemin receptor as SR470, a LRR-RLK closely related to known pattern recognition receptors (16). The early cellular responses to systemin do in fact resemble those typically triggered by microbe-associated molecular patterns (MAMPs) including an increase in cytosolic calcium, extracellular alkalization, plasma membrane depolarization, an oxidative burst, and ethylene production (17)(18)(19). By still unknown mechanisms, these early signaling events translate into the activation of the octadecanoid pathway for jasmonate production (20)(21)(22). The locally produced jasmonate signal is then perceived in distal leaves resulting in the systemic activation of defense responses (23). In addition to the slow-moving jasmonate signal, rapidly propagating electrical and calcium signals have been linked to the activation of jasmonate signaling and activation of wound response gene expression in systemic tissues (24)(25)(26).
Despite almost 30 years of research, it is still unclear how early responses at the plasma membrane are activated by systemin, and how these early signaling events including the influx of calcium, plasma membrane depolarization, and the production of reactive oxygen species are linked to downsteam events, like jasmonate production and defense gene activation. Since systemin is perceived by a receptor kinase, it is reasonable to assume that subsequent events are regulated by phosphorylation and can thus be captured by phosphoproteomics. Large scale phosphoproteomics was shown to be a powerful tool to identify global patterns and novel players in plant signaling networks. For example, time-courses studies of phosphorylation patterns in response to sucrose (27) has led to the identification of a novel kinase, SIRK1, regulating aquaporins (28). The cellular response to nitrate was well studied by large-scale phosphoproteomics to conclude about the cellular kinase network involved (29). Also for plant hormones, such as brassinolide, large-scale phosphoproteomics time-course studies revealed novel insights into the intracellular signaling network (30). Phosphoproteomics resulted in the identification of FERONIA as the receptor of rapid alkalization factor, and the plasma membrane H + -ATPase as a downstream target (31). In a similar approach, we used time-resolved phosphoproteomics in a systemin-responsive cell culture system to get further insight into the systemin signaling cascade and the sequence of early signaling events.

Experimental Design and statistical rationale
Cell suspension cultures of Solanum peruvianum (200 ml) were subjected to stimulation with 10 nM systemin (sys), 10 nM of the inactive Thr17Ala analog (A17), or injected with equal volume (1 ml) of water (w). For each treatment, a batch of cells was harvested after 2, 5, 15 and 45 minutes ( Figure  1A). Untreated cells were harvested as controls at time point 0. All treatments and time points were sampled from three independent replicates using independent batches of cells in each experiment. Data analysis was carried out jointly for all samples by averaging quantitative information for each time point and treatment across the three replicates. The same cell culture materials were used for phosphoproteome profiling as well as for activity assay. For phosphopeptide profiling, microsomal membranes were isolated from each sample, digested with trypsin, enriched for phosphopeptides and subjected to untargeted data-dependent acquisition by LC-MS/MS.

Cell suspension cultures and alkalization assay
The Solanum peruvianum cell suspension culture was kindly provided by Georg Felix (University of Tübingen). Cells were subcultured weekly by diluting 8 ml of the suspension in 70 ml fresh Murashige-Skoog medium containing Nitsch vitamins (32), 0.34 g/l KH2PO4, 5 mg/l 1naphthaleneacetic acid and 2 mg/l benzylaminopurine and grown on a rotary shaker at 26 °C with 120 rpm. For the alkalization assay continuous measurements of extracellular pH was performed in 10 ml of cultured cells 6-8 days after subculture. Synthetic systemin and systemin-A17 peptides (Pepmic; Suzhou, China) were added from a 200-fold concentrated stock solution in water.

Microsomal fraction extraction
The microsomal fraction was extracted according to (33) with minor modifications. A total of 1 g of harvested cells was homogenized in 10 ml extraction buffer (330 mM sucrose, 100 mM KCl, 1 mM EDTA, 50 mM Tris-MES pH 7.5 and fresh 5 mM DTT) in the presence of 1 % v/v proteinase inhibitor mixture (Serva) and 0.5 % of each phosphatase inhibitor cocktails 2 & 3 (Sigma-Aldrich, Germany) in a Dounce Homogenizer. At least 200 strokes were performed. The homogenate was filtered through four layers of miracloth and centrifuged for 15 min at 7500 × g at 4 °C. The supernatant was collected and centrifuged for 75 min at 48,000 × g at 4 °C. The pellet was washed with 5 ml of 100 mM Na2CO3 then again was centrifuged for 75 min at 48,000 x g at 4 °C. The microsomal membrane pellets used for LC-MS/MS were re-suspended in 100 µl UTU (6 M urea, 2 M thiourea and 10 mM Tris-HCl pH 8). Microsomal membrane pellets used for ATPase assay were re-suspended in 100 µl of re-suspension buffer (330 mM sucrose, 25 mM Tris-MES pH 7.5, 0.5 mM DTT). Protein concentrations were determined using Bradford assay (Roth, Germany) with BSA as protein standard. Samples were stored at −80 °C until further use.

Protein preparation for LC-MS/MS
Microsomal proteins were predigested for three hours with endoproteinase Lys-C (0.5 µg/µl; Thermo Scientific) at 37 °C. After 4-fold dilution with 10 mM Tris-HCl (pH 8), samples were digested with 4 µL Sequencing Grade Modified trypsin (0.5 µg/µl; GE Life Sciences) overnight at 37 °C. The reaction was stopped by addition of trifluoroacetic acid (TFA) to reach pH ≤ 3 before further processing.

Phosphopeptide enrichment
Phosphopeptides were enriched over titanium dioxide (TiO2, GL Sciences) with some modifications as described (34). TiO2 beads were equilibrated with 100 µl of 1 M glycolic acid, 80% acetonitrile and 5 % TFA. The digested peptides were mixed with the same volume of the equilibration solution and incubated for 20 minutes with 2 mg TiO2 beads with continuous shaking. After brief centrifugation, the supernatant was discarded and the beads were washed 3 times with 80 % acetonitrile and 1% TFA. Then the phosphopeptides were eluted with 240 µl of 1% ammonium hydroxide. Eluates were immediately acidified by adding 50 µl 10% formic acid to reach pH<3. Enriched phosphopeptides were desalted over C18-Stage tip (35) prior to mass spectrometric analysis.

LC-MS/MS analysis of peptides and phosphopeptides
Tryptic peptide mixtures were analyzed by LC/MS/MS using nanoflow Easy-nLC1000 (Thermo Scientific) as an HPLC-system and a Quadrupole-Orbitrap hybrid mass spectrometer (Q-Exactive Plus, Thermo Scientific) as a mass analyzer. Peptides were eluted from a 75 μm x 50 cm C18 analytical column (PepMan, Thermo Scientific) on a linear gradient running from 4 to 64% acetonitrile in 240 min and sprayed directly into the Q-Exactive mass spectrometer. Proteins were identified by MS/MS using information-dependent acquisition of fragmentation spectra of multiple charged peptides. Up to twelve data-dependent MS/MS spectra were acquired for each full-scan spectrum recorded at 70,000 full-width half-maximum resolution. Fragment spectra were acquired at a resolution of 35,000. Overall cycle time was approximately one second.
Protein identification and ion intensity quantitation was carried out by MaxQuant version 1.5.3.8 (36). Spectra were matched against the Tomato proteome (Solanum lycopersicum, ITAG2.4, 34725 entries) using Andromeda (37). Thereby, carbamidomethylation of cysteine was set as a fixed modification; oxidation of methionine as well as phosphorylation of serine, threonine and tyrosine was set as variable modifications. Mass tolerance for the database search was set to 20 ppm on full scans and 0.5 Da for fragment ions. Multiplicity was set to 1. For label-free quantitation, retention time matching between runs was chosen within a time window of two minutes. Peptide false discovery rate (FDR) and protein FDR were set to 0.01, while site FDR was set to 0.05. Hits to contaminants (e.g. keratins) and reverse hits identified by MaxQuant were excluded from further analysis. Raw files and identified spectra were submitted to ProteomeXChange and are available to the public under the accession number PXD010819.

Mass spectrometric data analysis and statistics
Reported ion intensity values were used for quantitative data analysis. cRacker (38) was used for label-free data analysis of phosphopeptide ion intensities based on the MaxQuant output (evidence.txt). All phosphopeptides and proteotypic non-phosphopeptides were used for quantitation. Within each sample, ion intensities of each peptide ions species (each m/z) were normalized against the total ion intensities in that sample (peptide ion intensity/total sum of ion intensities). Subsequently, each peptide ion species (i.e. each m/z value) was scaled against the average normalized intensities of that ion across all treatments. For each peptide, values from three biological replicates then were averaged after normalization and scaling. In case of nonphosphopeptides, protein ion intensity sums were calculated from normalized scan scaled ion intensities of all proteotypic peptides. A list of quantified phosphopeptides including normalized ion intensities, standard deviation and number of spectra is available as Supplementary Table 1.

H + -ATPase activity assay
Plasma membrane (PM) H + -ATPase activity was determined by resuspending 5 µg of microsomal proteins in 10µl resuspension buffer and 40 µl of reaction buffer containing 25 mM Tris-HCl pH 6.3, 50 mM KCl, 5 mM MgCl2, 5 mM CaCl2, 0.1 mM Ouabain (Thermo Fisher Scientific, USA), 1 mM NaN3, 100 nM Concanamycin A (Sigma-Aldrich, Germany), 0.02 % Brij58 (Sigma-Aldrich, Germany), 0.5 µg/µl BSA and 10 mM DTT. The reaction was initiated by the addition of 1 mM ATP, proceeded for 2 hours at room temperature, and stopped with 150 µl of stopping reagent (8% ascorbic acid, 2.1% HCl and 0.7% ammonium molybdate). After 5 min, 150 µl of 2% tri-sodium citrate and 2% acetic acid solution was added and incubated for 15 min at room temperature. The absorbance as a measure of inorganic phosphate was measured at 640 nm using plate reader (Tecan Spark). ATPase activity was measured with and without 400 µM orthovanadate (Na3VO4), and the difference between these two activities was attributed to the PM H + -ATPase.

Recombinant protein expression
The whole coding sequence of MPK2 (Solyc08g014420.2.1) was amplified by PCR using the primers ATTACCATGGCAGATGGTTCAGCTCCG (forward) and ATATCTCGAGCATGTGCTGGTATTCGGGATT (reverse) with S. lycopersicum leaf cDNA as template. The PCR product was cloned into pCR2.1-TOPO vector (Thermo Scientific) and amplified in E. coli. The MPK2 insert was excised with NcoI and XhoI (underlined in the primer sequence) and ligated to pET21d (+) (Merck; Darmstadt, Germany) to produce 14420HispET21d, which was transformed into E.coli strain BL21 RIL (Agilent Technologies; Waldbronn, Germany) to be expressed as C-terminally 6xHis-tagged recombinant protein. MPK2 was expressed in unmodified form and as constitutively active mutant (CA) with two amino acid substitutions D216G/E220A (39).
The PCR products were cloned into pET21d (+) as described above to produce His76100pET21d, which was transformed into E.coli strain BL21 RIL to be expressed as N-terminally 6xHis-tagged recombinant protein.
Purification of His-tagged proteins E.coli BL21 RIL cells were grown at 37 °C , 200 rpm to an A600 value of 0.8 and then induced with 1 mM isopropyl-β -d-thiogalactopyranoside (IPTG) for 4 h at 30 °C. Cells were pelleted by centrifugation and lysed in Bugbuster™ protein extraction reagent (Merck; Darmstadt, Germany) by continuous agitation for 20 min at room temperature. After centrifugation at 12,000 g for 20 min at 4 °C, the supernatants were collected. The supernatants were subjected to affinity purification on Ni-NTA agarose following the manufacturer's instructions (Qiagen; Hilden, Germany).

In-vitro kinase activity assay
Recombinant MPK2 (2 µg) was incubated with peptide substrate (GLDIETIQQSYTV; amount as indicated for each experiment) overnight in 50 µl of reaction buffer (20 mM HEPES pH 7.0, 10 mM MgCl2, 2 mM DTT and 0.1 mg/ml BSA) with 1 µM ATP at 30 °C. Kinase activity was measured as ATP consumption by addition of 50 µl of Kinase-Glo® Plus Reagents (Promega; Fitchburg, WI), 10 min incubation at 30°C, and subsequent recording of the decrease in luminescence (Tecan Spark luminometer). To detect the phosphorylated peptides by mass spectrometry, the reaction was performed as described with 20 µg of peptide substrate. The reaction was stopped by addition of TFA until pH ≤ 3 and reaction products were purified over C18-Stage tips (35) prior to mass spectrometric analysis.

In-vitro phosphatase activity assay
Recombinant PLL5 phosphatase (1 µg) was incubated with the phosphorylated peptide substrate (GLDIETIQQSYT(ph)V; amount as indicated for each experiment) in 50 µl of reaction buffer (20 mM Tris-HCl pH 7.5, 5 mM MgCl2, 1mM EGTA, 0.02 % β-mercaptoethanol and 0.1 mg/ml BSA) at room temperature for 30 minutes. Released phosphate was detected by addition of 150 µl of 8 % ascorbic acid, 2.1 % HCl and 0.7 % ammonium molybdate. After 5 min at room temperature the reaction was stopped by addition of 150 µl 2 % tri-sodium citrate and 2 % acetic acid (40). After incubation for 15 min at room temperature the reduced phoshomolybdenum complex was quantified at 640 nm using a plate reader (Tecan Spark) based on a standard curve of 0 -0.1µmol K2HPO4. To detect the unphosphorylated peptides by mass spectrometry, the reaction was performed as above with 200 µg of phosphorylated peptide substrate. The reaction was stopped by addition of TFA to reach pH ≤ 3 and reaction products were purified over C18-Stage tips (4) prior to mass spectrometric analysis.

The inactive A17 systemin analogon allows definition of systemin specificity
Early cellular responses to systemin include a rapid influx of calcium and depolarization of the plasma membrane which are necessary and sufficient for systemin signaling and induction of the wound response (17-19, 41, 42). Systemin-induced changes in plasma membrane ion permeability can be assayed conveniently in tomato (Solanum peruvianum) cell cultures as an alkalization of the culture medium (18) (Figure 2A). This cell culture system was used here to analyze dynamic changes of protein phosphorylation in the systemin signaling network.
Systemin responses, including extracellular alkalization, the induction of ethylene biosynthesis, and expression of wound response genes are inhibited by A17, a systemin derivative in which the penultimate threonine residue is substituted by alanine (18,43,44). In our cell culture system, stimulation with 10 nM systemin resulted in a rapid alkalization of the medium from pH 5.5 to pH 6.4 within about 8 minutes (Figure 2A). A17 was completely inactive at this concentration. When A17 was added to the systemin treatment, the alkalization response was reduced with increasing concentrations of A17 ( Figure 2B). To explain the antagonistic activity of A17, it was suggested that the peptide may interact with the binding site of the systemin perception system (43). Consistent with this proposition, binding to the putative receptor was shown to rely on the N-terminus of systemin, while the C-terminus, including Thr17, is required for activation (10). The receptor was recently identified as SYR1, and A17 was confirmed as a competitive systemin antagonist (16). By comparing phosphoproteomic responses triggered by A17 and by systemin, we will be able to filter the data set for systemin-specific phosphorylation changes. We therefore included the inactive A17 analog as a control in our experiments.

Description of the phosphoproteomics data set
A total of 3312 phosphopeptide species were identified. Most of the identified phosphopeptides were singly phosphorylated and 89.5% of identified phosphorylation sites were phosphoserines (Figure 1 B, C). Quantitative information under at least one treatment condition was obtained for 2960 phosphopeptide species matching 1309 protein models in at least one biological replicate. Thereby, the majority of phosphopeptides were identified at all time points ( Figure 3A). The high overlap of identified phosphopeptides across the time course indicates high analytical reproducibility and high potential for good quantitative coverage. At each sampled time point, the majority of phosphopeptides was identified under all treatment conditions ( Figure 3B). However, at 5 minutes, 15 minutes and 45 minutes, a large number of phosphopeptides was identified only after systemin stimulation. Far fewer phosphopeptides were found to be specific for the A17 or water control treatments.

Phosphorylation time profiles
Treatment of suspension cell cultures with systemin resulted in characteristic changes of phosphorylation patterns. K-means clustering was used to group these phosphorylation responses ( Figure 4A; Supplementary Table 2). A high number of 827 phosphopeptides was rapidly dephosphorylated upon systemin treatment (cluster A). Three clusters were characterized by a transient increase in phosphorylation, with 465 phosphopeptides peaking at 2 minutes (cluster B), 501 at 5 minutes (cluster C), and 705 at 15 minutes after systemin treatment (cluster D). Further 1001 phosphopeptides were classified as not responsive (cluster F). The stimulation of suspension cells with A17 also resulted in characteristic changes in phosphorylation ( Figure 4B). Interestingly, upon A17 stimulation no phosphopeptides displayed a transient increase in phosphorylation at 5 minutes, i.e. cluster C is missing. Instead, 1042 phosphopeptides showed an increase in phosphorylation at 45 minutes (cluster E). This response was not observed upon systemin treatment. Stimulation of the suspension cells with water resulted in response groups very similar to the clusters observed upon A17 treatment ( Figure 4C). Particularly, cluster C with a phosphorylation peak at 5 minutes also was missing, and cluster E with maximum phosphorylation at 45 minutes was observed. The missing cluster C and emergence of a "later" maximum in cluster E in both A17 treatment and water stimulation suggests that the control treatments A17 and water do trigger some cellular response, but weaker and delayed compared to systemin treatment.
To substantiate the observed phosphorylation profiles as physiologically relevant, we analyzed the activity of the plasma membrane H + -ATPase ( Figure 5) that has previously been implicated in the systemin-induced alkalization response (19). Two phosphopeptides of the C-terminal regulatory domain of plasma membrane ATPases Solyc07g017780 (AHA1) and Solyc03g113400 (LHA1) were identified with different time profiles under all three treatments. The peptides encompass the penultimate threonine residue, phosphorylation of which facilitates 14-3-3 protein binding and proton pump activation (45,46). In agreement with the observed alkalization of the growth medium upon systemin treatment, the peptides were slightly dephosphorylated at 2 minutes of treatment ( Figure 5 A, B). Dephosphorylation of the regulatory threonine residue was associated indeed with a measureable decrease in H + -ATPase activity at 2 to 5 minutes of systemin treatment ( Figure 5C). In contrast, A17 and water treatment induced a very transient increase in C-terminal phosphorylation and H + -ATPase activity at 2 minutes. We conclude that the time profiles presented in Figure 4 present precise quantitative measurements of phosphorylation status with corresponding effects on protein activity.

Processes involved in the systemin response
To better understand the processes at each time point, we carried out an over-representation analysis (Fisher Exact Test) of the phosphoproteins present in each of the clusters A to F under systemin stimulation ( Figure 6A). Proteins with functions in signaling (MapMan bin 30; p-value < 3.7E-19), transport (bin 34; p-value <2.3E-24), cellular processes (bin 31; p-value <3.0E-16), RNA related processes (bin 27; p-value <8.4E-6), and cell wall proteins (bin 10; p-value <3.9E-5) were highly over-represented in several clusters (Supplementary Table 3

Systemin-specific responses
To define precisely which of the phosphorylation events are systemin-specific, we compared for each phosphopeptide the phosphorylation time profile induced by systemin with those induced by A17 or water ( Figure 6B). For example, 649 out of 727 phosphopeptides in systemin-induced cluster A (dephosphorylated after 2 min), were classified in a "later" cluster (B to F) upon A17 treatment, or were not identified at all (7 phosphopeptides). Thus, these proteins were considered systemin-specific. Only 71 phosphopeptides showed the same rapid de-phosphorylation kinetics (i.e. also classified as cluster A) when treated with A17. In contrast, 665 of the 1001 phosphopeptides that were clustered as non-responsive to systemin treatment (cluster F), upon A17 treatment showed a peak phosphorylation in an earlier cluster. A total of 161 proteins were found with equal classification upon A17 and systemin treatment, and 83 phosphopeptides were not identified in the A17 treatment. In total, 1515 phosphopeptides (45% of all identified phosphopeptides) were classified as systemin-specific.  Table 4). Among the systemin-specific phosphopeptides, 105 and 18 phosphopeptides matched 78 protein kinases and 18 phosphatases, respectively (Supplementary Table 5).
Not all phosphopeptides from the 78 systemin-specific kinases met the requirement of at least four quantitative values across the time course and were therefore excluded from further analysis. The phosphorylation profiles of the peptides from the remaining 56 kinases and the 18 phosphatasederived peptides were then used to construct a correlation network of systemin-specific responses. Each phosphopeptide of a kinase or phosphatase was correlated against all other phosphopeptides from the functional categories over-represented in the systemin-induced response clusters ( Figure  6A). A stringent cutoff-value of r=0.85 was used to display individual proteins connected by the systemin-specific kinases or phosphatases ( Figure 7A, B). For each kinase/phosphatase in the networks, node size reflects its degree, i.e. the number of interacting target proteins. Specifically, for the proteins in cluster A with fast dephosphorylation response, phosphopeptides from a MAP-triple kinase (S(ph)SNATVEASGLSR, Solyc11g033270), a protein phosphatase 1 regulatory subunit (LSS(ph)PAAQS(ph)PSVSTK, Solyc02g070260) and a protein phosphatase C (LTNS(ph)PDVEIK, Solyc10g005640) had the highest degree in the correlation network (210, 134, 118 respectively). These proteins had highest numbers of systemin-specific correlation partners. Other kinases with high degree (number of interactions) were an uncharacterized kinase family protein (AIS(ph)LPS(ph)SPR, Solyc02g076780, degree 118), a calcium-dependent protein kinase (CDPK)related kinase (TES(ph)GIFR, Solyc10g078390, degree 115), and a receptor like kinase (S(ph)FENEIR, Solyc09g072810, degree 110). Among proteins displaying transient phosphorylation at 2 minutes (cluster B) we found CBL-interacting protein kinase (T(ph)TCGTPNYVAPEVLSHK, Solyc04g076810, degree 112), serine/threonine protein kinase (NQGS(ph)PSDTCSESDHK, Solyc01g108920, degree 110) and cell division protein kinase (DRLDS(ph)IDGK, Solyc07g063130, degree 102) with highest degrees. For proteins with transient phosphorylation at 5 minutes (cluster C), the protein kinases with highest degree were the ethylene receptor (GGS(ph)QTDSSISTSHFGGK, Solyc05g055070, degree 103), and two receptor kinases (SIDLFTDVSEEGLS(ph)PR, Solyc09g083210, degree 49; ENPGDS(ph)GSLVQPGHDIEK, Solyc12g010740, degree 44). Kinases with highest degree among the proteins with transient phosphorylation at 15 minutes (cluster D) were a serine/threonine protein kinase (SEEVDGS(ph)SSIR, Solyc11g033270, degree 210), the activating motif of a mitogen activated protein kinase (VTSETDFM(ox)T(ph)EY(ph)VVTR, Solyc08g014420, degree 128), and a uncharacterized protein kinase (S(ph)AAGTPEWM(ox)APEVLR, Solyc09g009090, degree 122).
Interestingly, these protein kinases with high degree often were identified with several phosphorylation sites which showed different response profiles during the systemin treatment. For example, the MAP-triple kinase (MAPKKK) Solyc11g033270 (degree 210) was identified with two phosphorylation sites outside the kinase domain: S(ph)SNATVEASGLSR (Ser859) showed rapid dephosphorylation (cluster A) and SEEVDGS(ph)SSIR (Ser345) showed transient phosphorylation at 15 minutes (cluster D).
Time profile of the rapidly dephosphorylated MAPKKK peptide S(ph)SNATVEASGLSR showed highest correlation with dephosphorylated peptides of an ubiquitin ligase (LEEGS(ph)SPEQR, Solyc05g007820), a YAK1 homologous protein kinase (TVYSY(ph)IQSR, Solyc03g097350) and a RNA recognition motif containing protein (SSDS(ph)QELTTTELK, Solyc02g062290), and a protein phosphatase 2C family protein (LTSS(ph)PEAEIK, Solyc06g065580). The finding that kinases and phosphatases correlate with the same phosphopeptide suggest kinase-phosphatase activation/inactivation loops upon systemin stimulation. At 15 minutes of systemin treatment, the MAPKKK phosphopeptide SEEVDGS(ph)SSIR of protein kinase Solyc11g033270 displayed a transient increase in normalized intensity. This correlated highest with a transient phosphorylation of a calmodulin domain containing protein (DYS(ph)ASGYSSR, solyc06g053530) and a transcriptional regulator (S(ph)PIPLEAVK, Solyc01g097890). However, while our data set reveals correlative relationships between the phosphorylation sites that may be indicative of enzyme-substrate interaction, it must be kept in mind that the data does not directly allow to derive the activity change of the respective proteins without further experiments.
Within this network, we then displayed the degree distribution of phosphopeptides from systeminspecific kinases/phosphatases, other kinases/phosphatases and the putative targets for systemininduced clusters A to F ( Figure 7C, D). It becomes apparent that systemin-specific kinases in clusters B, C, and D (which showed a transient phosphorylation increase) showed a higher degree than kinases without systemin specificity, but kinases in generally had higher degree than their target proteins ( Figure 7C). In contrast, for phosphatases, higher than average degree was observed for nine systemin-specific phosphatases in cluster A (rapid dephosphorylation response). In cluster F (no systemin induced phosphorylation response), higher than average degree was observed for nine other phosphatases (those not linked to the systemin response). This observation supports the conclusion that the kinases identified as systemin-specific, are indeed responsible for the phosphorylation events that peak at 2 minutes (cluster B), 5 minutes (cluster C) and 15 minutes (cluster D) of stimulation, while systemin-specific phosphatases mediate rapid dephosphorylation responses (cluster A). A full list of kinase-target relationships is available in Supplementary Table 6. Selected candidate targets for systemin-specific kinases and phosphatases are addressed below.
The interaction of the C-terminal peptide of H + -ATPase LHA1 with candidate kinase MPK2 and candidate phosphatase PLL5 was investigated by in vitro activity assays (Figure 8 B-E). Recombinant MPK2 was expressed either in wild-type form (WT) or as a constitutively active mutant (CA) with D216G and E220A amino acid substitutions (39). Kinase activity, assayed as ATP consumption, was observed only in presence of GLDIETIQQSYTV peptide substrate, and was higher for CA compared to WT ( Figure 8B). Increasing the substrate peptide concentration from 10 µM to 200 µM resulted in increasing kinase activity ( Figure 8C). Phosphorylation of the C-terminal LHA1 peptide substrate at T955 by recombinant MPK2 was confirmed by mass spectrometric analysis (supplementary figure  1A,C). Recombinant PLL5 was able to dephosphorylate synthetic phosphopeptide GLDIETIQQSY(pT)V in vitro ( Figure 8D). PLL5 activity, assayed as release of inorganic phosphate, was observed only in presence of the phosphopeptide substrate, and increased with increasing substrate concentration ( Figure 8E). De-phosphorylation of the LHA1-derived phosphopeptide substrate by PLL5 phosphatase was confirmed by mass spectrometric analysis (supplementary figure 1B,C). Interestingly PLL5 phosphatase was identified in systemin-induced cluster A (rapid dephosphorylation), while MPK2 was identified in systemin-induced cluster D with maximum phosphorylation of peptide VTSETDFM(ox)T(ph)EY(ph)VVTR at 15 minutes. The identified phosphopeptide contained the typical MAP-Kinase activating motif (pT)E(pY). Thus, the phosphatase may well be involved in dephosphorylation of LHA1 at early systemin stimulation, while MPK2 is a candidate kinase for rephosphorylation ( Figure 5A) with concomintant increase in H+-ATPase activity ( Figure 5C) at later time points (15 min) of systemin treatment.

Respiratory burst oxidases
Another protein with highly systemin-specific phosphorylation response was the respiratory burst oxidase. Phosphopeptides were identified for two isoforms of respiratory burst oxidases (DVFSEPS(ph)QTGR, Solyc03g117980; NLSQM(ox)LS(ph)QK, Solyc06g068680). Potential interacting kinases include receptor kinase (Solyc09g064270 and Solyc02g07000) and a CBL-interacting protein kinase (Solyc06g068450, r=0.90). Phosphorylation of the respiratory burst oxidase occurred at early time points (2 minutes) with a strong transient increase of phosphorylation for both peptides. No such transient increase in phosphorylation of the respiratory burst oxidase was observed under A17 treatment (Supplementary Figure 2A).

Processes at the cell wall
Among the target proteins with most significant systemin-induced phosphorylation responses were cellulose synthase proteins. Four phosphopeptides matching cellulose synthase like proteins (SSS(ph)RLNLSTR Solyc04g077470, SHS(ph)GLMR Solyc08g076320, SSS(ph)ESGLAELNK Solyc09g057640, GLIDSQSLSSS(ph)PVK Solyc09g075550) were identified as potential targets of systemin-specific kinases and phosphatases (Supplementary Figure 2B). Systemin-induced phosphorylation changes of cellulose synthase-like proteins were either rapid dephosphorylation (cluster A) or transient phosphorylation at 2 minutes (cluster B). The transient phosphorylation of an uncharacterized receptor kinase at 2 minutes (NNTS(ph)SVSPDSVTAK, solyc12g036330) showed strongest correlation with the transient phosphorylation of peptide S(ph)SEGDLTLLVDGKPK of cellulose synthase-like protein solyc04g077470, while dephosphorylation of peptide GQLPS(ph)GQVVAVK matching an uncharacterized receptor-like cytoplasmic kinase (Solyc05g024290) showed strongest correlation with dephosphorylation of peptide GLIDSQSLSSS(ph)PVK of cellulosesynthase-like protein (Solyc09g075550). These two peptides also showed a high correlation to a POLTERGEIST-LIKE 1 phosphatase (Solyc08g077150, peptide FVS(ph)PSQSLR; r=0.92).

Vesicle trafficking and G-Protein signaling
In total, 52 kinases and 8 phosphatases were found with high correlations to phosphopeptides of vesicle trafficking proteins. While phosphopeptides of adaptins and syntaxins were rapidly dephosphorylated (cluster A), most phosphopeptides of vesicle-associated membrane proteins (VAMPs) were transiently phosphorylated at 2 minutes of systemin stimulation (cluster B). For example, phosphopeptide NS(ph)DDGYSSDSILR in cluster A matching a KEU-homolog (Solyc10g008490) showed strong correlation (r=0.99) to phosphopeptide VNS(ph)LVQLPR matching the catalytic subunit of a protein phosphatase (Solyc10g008490). Phosphopeptide VVYIS(ph)PHS(ph)SPGHSEDAFK of VAMP Solyc04g015130 with transient phosphorylation at 2 minutes showed high correlation with a protein phosphatase 2C Solyc07g053760 and protein kinase Solyc06g068920.
Several phosphopeptides of small GTPase proteins were also identified and these proteins may also have functions in vesicle budding and fusions. Most prominent members of this group included the ARF GTPase activator Solyc05g023750 and RAB GTPase activator Solyc12g009610. The phosphopeptide NS(ph)SELGSGLLSR of the ARF-GTPase activator was identified as being rapidly dephosphorylated (cluster A), while phosphopeptide TLS(ph)SLELPR of the RAB-GTPase activator Solyc12g009610 was transiently phosphorylated at 2 minutes of systemin treatment. Phosphopeptide AELS(ph)VNR of a RAB-GTPase Solyc09g098170 was transiently phosphorylated at 5 minutes. Phosphopeptide NS(ph)SELGSGLLSR correlated highest with phosphopeptide ALIGEGS(ph)YGR of receptor kinase Solyc12g098980, phosphopeptide TLS(ph)SLELPR showed highest correlation with receptor kinase Solyc07g056270. Phosphorylation of the RAB-GTPase Solyc09g098170 correlated highest with receptor kinase Solyc07g055130 and the protein phosphatase 2C Solyc05g055790.

Transcription factors
Phosphopeptides of several transcription and splicing factor families (bin 27) were over-represented with similar abundance in all systemin-response clusters (supplementary Table 2). The highest number of identified phosphopeptides and highest degree for this functional protein group in the network (37) was observed for the RNA-processing factor SERRATE (Solyc01g009090) with matching phosphopeptides S(ph)PPFPPYK and NRS(ph)PPHHPPGR. Phosphopeptide S(ph)PPFPPYK was transiently phosphorylated at 2 minutes and showed best correlation to phosphorylation sites of kinases CTR1 (Solyc09g009090), a MAP-triple kinase (Solyc12g088940), and a CBL-interacting kinase (Solyc04g076810). Phosphopeptide NRS(ph)PPHHPPGR was transiently phosphorylated at 5 minutes of systemin treatment. Phosphorylation sites of receptor kinases Solyc02g023950 and Solyc09g072810, a CDK-related kinase Solyc10g078390 and a MAP-triple kinase Solyc11g033270 showed high correlation with that second SERRATE phosphorylation site.

Discussion
The aim of our work was to (i) identify early events induced by systemin and (ii) identify kinases and phosphatases involved in these signaling events. The maturation of the signaling peptide systemin is beginning to be understood (9), and the primary receptor of systemin was recently identified (16). While individual early responses to systemin perception, such as an increase in cytosolic calcium (17), ROS production (16,48), the modulation of the H + -ATPase activity (19), or the activation of MAPkinase signaling (49) have been identified, the proteome-wide effects of systemin treatment were not yet studied.

The use of inactive analogon A17 to define systemin-specific responses
Antagonistic and inactive peptides-such as mutant peptide variants, chemically modified peptides or peptide mimetics, are important tools for the characterization of peptide-receptor interaction and to discover endogenous peptide function (50)(51)(52). For example, the dominant-negative effect produced by antagonistic peptides can be used for the characterization of peptide function when suitable mutants are not available, and in the case of genetic redundancy (53,54). The use of inactive analoga in studying receptor signaling pathways provides an elegant way to define specific responses and rule-out common unspecific treatment responses. Comparing the responses to systemin and the antagonistic A17 analog, approximately half of the identified phosphopeptides could be identified as not specific to the systemin treatment and were excluded from the reconstruction of the systemin signaling network. Including the A17 analog thus allowed us to gain insight into systemin-specific protein signaling processes.

Systemin-induced modulation of the plasma membrane H+-ATPase
Systemin-induced alkalization of the growth medium has been well described (18,19,51). Alkalization of the growth medium is a typical process connected with inhibition of the plasma membrane H + -ATPase, as it also occurs upon stimulation by flagellin and other microbe-associated molecular patterns (31,55,56). By pairwise correlation of phosphopeptide abundance profiles we identified the H + -ATPase LHA1 as putative target of 12 kinases and two phosphatases, while receptor like kinase Solyc11g011020 was the only kinase putatively involved in the phosphorylation of LHA1. These kinase-substrate predictions were confirmed by in-vitro kinase and phosphatase assays. MPK2 and PLL5 phosphatase were found to indeed act on the respective peptide of LHA1. While such verification can only be done for a limited number of kinases or phosphatases, the results obtained for MPK2 and PLL5 phosphatase confirm the general validity of our correlation analysis. We thus consider the correlation analysis for prediction of kinase/phosphatase-substrate interactions as a powerful means to predict new regulatory circuits in the systemin signaling network. Thereby we identified a novel connection of the MPK pathway with the plasma membrane H + -ATPase LHA1. MPK2 is activated by phosphorylation of the typical TEY motif at 15 minutes of systemin treatment and, therefore, is a primary candidate for re-phosphorylation of the H + -ATPase at later times of systemin stimulation. Phosphorylation by MPK2 occurred at the penultimate regulatory threonine residue resulting in activation of the proton pump. MPK2-mediated activation of LHA1 is thus likely to be responsible for re-acidification of the extracellular space dampening the initial alkalization response.
While knowledge of specific interactions of phosphatases with their substrates is still limited in plants, several phosphatases have been implicated in the regulation of the plasma membrane proton pump. This includes a type 2A phosphatase activity partially purified from maize roots (57), an unidentified type 2C phosphatase in Arabidopsis seedlings and Vicia faba guard cells (58), and PP2C clade D phosphatases that are involved in proton pump regulation during auxin-induced acid growth of Arabidopsis hypocotyls (59). PP2C-D3 (alias PP2C38; At4g28400) also was shown to dephosphorylate BIK1 attenuating immune signaling through the FLS2 and EFR pattern recognition receptors (60). Likewise, EFR-mediated signaling also is downregulated by Arabidopsis PLL4 and PLL5 phosphatases, which belong to PP2C clade C (61). We show here that the tomato PLL5 homolog, is (i) involved in systemin signaling and (ii) dephosphorylates LHA1. The precise mechanisms of PP2C phosphatase regulation and signaling are still unclear, especially if and how they are activated or inactivated by reversible phosphorylation (62). Under the assumption that similar to PP2A proteins (62), PLL5 phosphatase may be activated by de-phosphorylation (consistent with its identification in systemin-induced cluster A), it is a primary candidate for proteins involved in the downregulation of H + -ATPase activity and the resulting alkalization of the growth medium.

Reconstruction of early vs. late responses from plasma membrane to the nucleus
The recently identified plasma membrane located LRR-RLK receptors for the systemin peptide were not identified among the 3312 phosphopeptides from this work, while several other receptor kinases, among them also the phytosulfokine receptor PSKR2 (Solyc07g063000) or FERONIA (Solyc09g015830) were quantified in our experiments. Possibly, some important receptor kinases were not enriched in the membrane purification protocol due to their strong embedding within the cell wall. Nonetheless, our work identified 57 systemin specific kinases and 18 phosphatases to be involved in systemin signaling, likely downstream of the primary systemin receptor SYR1 (16). Substrates for these kinases and phosphatases were identified by pairwise correlations. Correlation analysis was previously used to reconstruct topology of phosphorylation networks induced by external nutrient supply (63). Thereby, plasma membrane proteins were identified as signal induction layers within the network, while cytosolic and nuclear proteins were considered as an effector layer.
Considering the absence of phosphopeptides from recently identified systemin receptors, the earliest harvesting time point of 2 minutes may have been already too late to capture the immediate first systemin responses. Ligand-induced phosphorylation of ligand-perceiving receptor kinases is known to occur within seconds (64). The identified phosphorylation patterns at 2 minutes may thus represent early responses of the systemin signaling cascade. This point of view is supported by the observation of rapid (de)phosphorylation changes for several cytosolic proteins such as a CBLinteracting kinase, MAP-triple kinase, and CDPKs already at 2 minutes (clusters A and B), while a large number of receptor kinases with systemin specific phosphorylation patterns peak only after 5 minutes. Also observed very rapidly after systemin stimulation are plasma membrane-located proteins that are known to belong to the systemin signaling cascade, including the proton pump that is rapidly de-phosphorylated after systemin treatment, and two RBOHD isoforms showing a strong increase of phosphorylation after 2 minutes, compared to G-protein signaling or transcription factor phosphorylation at later times.

Conclusions
Our study presents 56 systemin specific kinases and 18 systemin specific phosphatases. Putative substrates were identified by pairwise correlations of systemin induced phosphorylation time profiles. Thereby we identified novel candidates for the systemin signaling pathway. Two of these candidates, MPK2 and PLL5 were confirmed as kinase and phosphatase for reversible phosphorylation of the C-terminal regulatory site of the plasma membrane H + -ATPase. We highlighted several proteomic processes upon systemin stimulation ranging from plasma membrane processes to cytosolic and nuclear signaling, which will provide a valuable resource for further indepth functional studies in future elucidation of systemin signaling cascades.