Untargeted metabolomics unravels functionalities of phosphorylation sites in Saccharomyces cerevisiae

Abstract

Background

Coordinated through a complex network of kinases and phosphatases, protein phosphorylation regulates essentially all cellular processes in eukaryotes. Recent advances in proteomics enable detection of thousands of phosphorylation sites (phosphosites) in single experiments. However, functionality of the vast majority of these sites remains unclear and we lack suitable approaches to evaluate functional relevance at a pace that matches their detection.

Results

Here, we assess functionality of 26 phosphosites by introducing phosphodeletion and phosphomimic mutations in 25 metabolic enzymes and regulators from the TOR and HOG signaling pathway in Saccharomyces cerevisiae by phenotypic analysis and untargeted metabolomics. We show that metabolomics largely outperforms growth analysis and recovers 10 out of the 13 previously characterized phosphosites and suggests functionality for several novel sites, including S79 on the TOR regulatory protein Tip41. We analyze metabolic profiles to identify consequences underlying regulatory phosphorylation events and detecting glycerol metabolism to have a so far unknown influence on arginine metabolism via phosphoregulation of the glycerol dehydrogenases. Further, we also find S508 in the MAPKK Pbs2 as a potential link for cross-talking between HOG signaling and the cell wall integrity pathway.

Conclusions

We demonstrate that metabolic profiles can be exploited for gaining insight into regulatory consequences and biological roles of phosphosites. Altogether, untargeted metabolomics is a fast, sensitive and informative approach appropriate for future large-scale functional analyses of phosphosites.

Background

Adaptation to changes in the environment is crucial for the survival of any organism and adjustment of protein levels and activities is pivotal for this process. Long-term adaptations via (post-)transcriptional regulation or irreversible protein modification typically affect protein abundance. Fast adaptations are mediated by allosteric regulators and a chemically diverse group of covalent post-translational modifications (PTMs). In eukaryotes, phosphorylation stands out among PTMs for the broad scope of regulation and the complexity of its network with over 160 known kinases and phosphatases in Saccharomyces cerevisiae [1]. Recent advances in proteomics allow for deep interrogation of phosphoproteomes, leading to large catalogues of phosphorylation sites (phosphosites) in different species [2–4], suggesting that as many as 45% of all proteins in eukaryotes may be phosphorylated [5]. While the number of mapped phosphosites is ever increasing, we still rarely know whether phosphorylation events have a functional role. A substantial fraction of all phosphosites was suggested to be non-functional, resulting from stochastic off-target kinase activity [6, 7]. Typically, only a handful of all observed phosphosites in any organism has been functionally validated and characterized [5, 8] because the success in phosphoproteomics in detecting phosphosites has largely outpaced the performance and throughput of the available approaches for analysis of phosphosites functionality.

Present workflows for functional analysis of phosphosites typically start with phosphoproteomics to generate catalogues of modified amino acids. Computational approaches then distinguish between relevant and non-functional phosphosites on the basis of different criteria including kinase recognition motives [9, 10], conservation patterns of phosphosites across species [6, 7], correlations of changes in phosphorylation state with metabolic flux [11], or temporal patterns of phosphosite appearances upon disturbance [12]. After prioritization of interesting phosphosites, their functionality remains to be determined in laborious follow-up experiments. Typically, these experiments take advantage of site specific phosphosite mutants and include a combination of phosphorylation assays, growth and viability assays, localization studies, expression studies and enzymatic in vitro assays [13–20]. Targeted measurements of intracellular metabolites have been used in some cases to investigate phosphoregulation of metabolism [11, 16, 20]. All these methods have been powerful in elucidating the importance of single phosphosites. However, they are not suitable for high throughput analysis as the readout of most approaches is specific to the investigated protein and would imply tedious laboratory work to make them applicable in larger scale. Although growth and viability assays allow high throughput, their sensitivity for assessing functionality is unclear. Given that phosphorylation often only fine-tunes the activity of a protein, interfering with phosphorylation might not necessarily penetrate to the phenotypic level.

The recent developments in genome editing render the generation of large numbers of mutants with modified phosphosites for functional analysis feasible [21–24]. Hence, the limitation shifts to a suitable high throughput approach for subsequent mutant analysis, in particular when phenotypes are more subtle, as can be expected for most point mutations. In principle the various omics technologies would be applicable, but in practice they are typically too laborious or expensive. A notable exception is mass spectrometry-based metabolomics. In particular, it offers the potential to identify small differences in metabolite profiles well before a change of the phenotype is observed [25, 26]. Since many of the phosphoproteins listed in phosphosite databases are part of the metabolic network, metabolomics is a promising approach for probing their functionality.

Here we investigate the potential of untargeted metabolomics for higher throughput identification and characterization of phosphosite functionality. We introduced pointmutations into 26 phosphosites from 25 proteins in S. cerevisiae, several of which were known to be functionally regulated by phosphorylation. The proteins included metabolic enzymes and regulatory proteins implicated in the TOR or HOG signaling, two pathways closely related to metabolism [27–29]. The phosphoresidues were mutated to amino acids that mimic a constitutively dephosphorylated and, for a selected subset, phosphorylated state, resulting in a total of 32 phosphomutants. Applying untargeted metabolomics, we retrieved most previously characterized phosphoproteins as functionally relevant and suggest functionality for additional phosphosites. We also showed that untargeted metabolomics aids in unveiling the biological roles of functional phosphosites in regulating metabolic pathways. Our results demonstrate that untargeted metabolomics is a highly sensitive method for high throughput detection of functionalities of phosphosites.

Results

Phenotypic analysis reveals functionalities of strong-impact phosphosites

For the purpose of defining a set of phosphomutants for phenotypic and metabolic analysis, we selected single phosphosites or multiple phosphosites in close proximity (hereafter jointly termed as phosphosites) on metabolic enzymes or proteins involved in TOR or HOG signaling. The 26 chosen phosphosites were located on 25 proteins that represent different functional classes, such as metabolic enzymes, kinases or transcription factors (Table 1 and Additional file 1). The two signaling pathways were selected for their metabolic function, rendering them amenable for a metabolomics approach. TOR signaling controls the cellular response to nitrogen availability, having a substantial impact on the state of metabolism when metabolizing nitrogen sources of different quality [27, 29]. The main role of the HOG signaling pathway is sensing of extracellular hyperosmolarity and eliciting a stress response, which in S. cerevisiae primarily consists of a metabolic adaptation to increase the concentration of the osmolyte glycerol together with cell cycle and transcriptional modulation [28]. The total set of the 26 phosphosites contained 13 functionally characterized and 13 uncharacterized phosphosites. We classified phosphosites as characterized when a phosphosite mutant has been reported to have an effect in any type of functional analysis. Most uncharacterized phosphosites were selected based on differential phosphorylation data from phosphoproteomics studies in either different growth conditions or upon kinase deletion [4, 11, 30] (Additional file 1). The chosen phosphosites were mutated to alanine whenever serine or threonine was the phosphoresidue, or phenylalanine if tyrosine was the phosphoresidue, thereby abolishing phosphorylation of the target amino acid (Table 1). For six out of the 26 phosphosites the residues were additionally mutated to glutamic acid to mimic a constitutively phosphorylated state of the protein, resulting in a total of 32 mutants. Throughout this study, mutants with abolished phosphorylation are referred to as “OUT” mutants, while phosphomimic mutants are referred to as “IN” mutants. Up to three independent mutant replicates per phosphosite were generated (Additional file 1). Phosphomutant replicates with broad metabolic effects exhibited highly similar metabolic profiles (Additional file 2: Figure S1), suggesting that genetic manipulations were precise.

Table 1 Phosphoprotein mutant set

To identify condition-specific regulation, the 32 phosphomutants were grown on four combinations of carbon and nitrogen sources that require different activities of the mutated proteins, such as glycolytic or gluconeogenic fluxes, nitrogen sources of different quality, and low osmotic stress on 0.22 M sodium pyruvate (Table 2). We quantified growth rates and maximum cell density and qualitatively assessed lag times and other abnormalities in the growth curves of microtiter plate cultures (Fig. 1a and Additional file 3). Six mutants in four different phosphosites exhibited an aberrant phenotype under at least one condition (Fig. 1b). Three of these phosphosites were already previously characterized, but functionality of the Tip41 phosphosite was hitherto unknown. The phosphomutants of Pda1 and Pbs2 showed varying effects under most conditions, whereby the IN mutant was more strongly affected in both proteins. The IN mutant of the kinase Pbs2 showed a very similar behavior as its target Hog1 OUT, supporting the hypothesis that phosphorylation of Pbs2 S508 inactivates the protein, while in Hog1, phosphorylation of S174 and T176 in the activation loop activates this kinase [15, 31, 32]. Overall, we retrieved three out of the 13 known phosphoregulated proteins under the tested conditions. Thus, growth phenotypes provide strong evidence for phosphosite functionality but would require testing of many conditions.

Table 2 Carbon and nitrogen sources of media for growth and metabolomics experiments

Fig. 1

Growth phenotype of S. cerevisiae phosphomutants. a Example curves for analyzed growth traits of wild-type (green) and mutant (red). b Growth phenotype in different growth conditions. Results that differ significantly from the wild-type are highlighted. Red and green boxes indicate higher or lower results for the respective growth trait of the mutant compared to the wild-type. Black boxes indicate a general difference in shape. Mutants with a growth phenotype are highlighted in blue. Abbreviations: μ: growth rate. CD: maximum cell density in one growth phase. L: lag time. S: abnormalities in shape of growth curve

Untargeted metabolomics reveals changes in phenotypically silent phosphomutants

Next, we tested whether metabolic functionality of phosphosites whose genetic perturbation does not propagate to physiological phenotypes could be identified by untargeted metabolomics. Intracellular metabolites were extracted with hot ethanol from exponentially growing cultures in the four before used conditions using 96-well cultivation. The extracts were injected into a Time Of Flight mass spectrometer (TOF-MS) and metabolic profiles in the m/z range of 50 to 1000 Dalton were recorded [33]. Across all conditions, 122–259 ions were annotated to 183–342 metabolites using a genome-wide reconstruction model of S. cerevisiae [34] (Additional file 2: Table S1). The data was processed to remove intensity drifts during measurements and an OD-specific fold change was calculated for every ion of each metabolic profile. The fold changes of the replicates of one mutant were compared to the fold changes of all other samples using a 2-sample t-test and the median log2 fold change over all replicates was determined. Ions with a corrected p-value < 10⁻³ and a |log2 fold change| > 0.3785 were considered as changing significantly. The selected log2 fold change cutoff corresponds to a 30% change and represents the 97.5% quantile across the fold changes of all ions from the wild-type dataset (Additional file 2: Figure S2A).

Typically, metabolite concentrations changed less than two-fold in the phosphomutants (Additional file 2: Figure S2B) and the majority of changes were in arginine, proline, pyrimidine, purine or lysine metabolism (Additional file 2: Figure S4). While 31% of the mutants exhibited no significant (p-value < 10⁻³, |log2 fold change| > 0.3785) metabolic response, most exhibited changes in some metabolites, and five of the six mutants with physiological phenotypes (Hog1 OUT, Pbs2 IN, Pda1 OUT, Pda1 IN and Tip41 OUT) showed changes in more than 15 metabolites in at least one condition (Fig. 2a). In many cases these broadly changing metabolic profiles were presumably the consequence of altered growth rates (Fig. 2b). Generally, phosphomutants with previously characterized phosphosites featured more changes than the unknown subset, suggesting that the latter contains more nonfunctional sites (Fig. 2c). Metabolic profiles of phosphosites with strong metabolic impact often showed a tendency towards condition-specific effects, demonstrating the regulatory role of phosphosites under particular conditions. Hog1 OUT for example showed broad effects only when grown on pyruvate as the carbon source, presumably because of the elevated osmotic stress in this condition, and Tip41 showed almost exclusively effects when growing on proline as sole nitrogen source (Fig. 2a). Overall, our data suggest that metabolomics can detect subtle consequences of deregulated phosphorylation that is without detectable growth phenotypes.

Fig. 2

Overview on changing ions in metabolic profiles. a Number of changing ions of each mutant and condition. Only mutants with significantly changing ions are depicted (|log2 fold change| > 0.3785 and p-value < 10⁻³). b Relation between number of changing ions and mutant growth rate. c Distribution of changing ions in different mutant subclasses with no (none), little (1–4), intermediate (5–19) or many (>19) changing ions

Metabolic profile analysis indicates phosphosite functionality

While we did not detect changes in direct reactants of the 15 phospho-mutated enzymes after applying a strict cutoff, we found metabolite changes for 12 mutants to be enriched in specific pathways, suggesting phosphorylation to be functionally important for these pathways (Additional file 4). For example, the OUT mutant of Yap4, a transcription factor implicated in stress responses, showed changes enriched in purine metabolism when growing on glucose/glutamine (p-value 0.0007). Metabolic effects of the Pda1 phosphomutants predominantly affected arginine and proline metabolism (p-value 0.005 for Pda1 IN and p-value 0.0004 for Pda1 OUT respectively in glucose/NH₄ ⁺, p-value 0.0003 for Pda1 IN in glucose/glutamine). A particularly clear case was Lys20, one of two homocitrate synthase isoenzymes that catalyze the conversion of oxoglutarate to homocitrate, the first step in lysine biosynthesis. Lysine metabolism was found to be exclusively enriched (p-value 0.006) in Lys20 OUT when growing on glucose/proline. While lysine biosynthesis is known to be subject to transcriptional and allosteric control, it was recently shown that Lys20 activity is also regulated by phosphorylation [26, 35]. Metabolite responses in Lys20 OUT were observed only during growth on the poor nitrogen source proline and significant (p-value < 10⁻³, |log2 fold change| > 0.3785) changes were restricted to the final pathway product lysine and its precursor saccharopine (Fig. 3). Consistent with previous results [26], increased lysine levels in Lys20 OUT suggest phosphorylation at position S395/T396 to deactivate the enzyme. Since lysine biosynthesis flux is rather low during slow growth on proline, our data indicates that the pathway capacity is down-regulated by phosphorylation in this condition.

Fig. 3

Metabolite level changes in lysine biosynthesis of the Lys20 OUT mutant growing on glucose/proline. Coloring indicates the strength and direction of the fold changes. Significant fold changes are highlighted with a bold frame (p-value < 10⁻³, absolute log2 fold change > 0.3785). Metabolites in grey were not detected

To identify more global impact of phosphosites on metabolism, we used two correlation approaches to analyze the metabolic profiles. In the first, we correlated log2 fold changes between all 496 possible mutant-pairs for each condition. The correlation coefficients of all mutant-pairs showed a normal distribution around 0, with a mean of 0.014 (Additional file 2: Figure S5A). The four correlation coefficients from each mutant-pair were compared against the distribution of all derived coefficients to identify significantly correlation pairs. Five mutant-pairs correlated significantly across the four conditions (p-value < 0.01), suggesting these phosphosites to exhibit similar effects on metabolism (Additional file 5). Two of these five pairs were mutant variants of the same protein (Gln3 site1 OUT – Gln3 site2 OUT and Rim15 OUT – Rim15 IN), emphasizing similar functionality of the introduced mutations. Rim15 OUT, implicated in cell cycle regulation in response to nutrient, and Pda1 IN were negatively correlated in all conditions, whereby the strongest opposing effects were largely in amino acid metabolism. The two other cases were Pfk2 OUT associating with Ste50 OUT as well as Yap4 OUT.

In the second correlation, we quantified the similarity of each mutant with itself across all four conditions to identify phosphosite functionalities based on subtle but consistent changes in metabolites across conditions. The correlation coefficients of all condition-pairs showed an average above zero, emphasizing that our phosphomutant set contains mutants with condition-independent metabolic effects (Additional file 2: Figure S5B). Out of the eight phosphomutants that were recovered in this analysis, we already described strong metabolic changes predominantly under one condition for three mutants (Pda1 OUT, Hog1 OUT and Pbs2 IN, see Additional file 6). The calculated correlation coefficients were particularly weak for Hog1 OUT and Pbs2 IN. We suspect that despite the predominant effects in specific conditions, this weak similarity is driven by underlying subtle, but condition-independent, metabolic changes. One such example in Pbs2 IN is the slight depletion of glycerol across conditions, emphasizing the role of Pbs2 in controlling glycerol accumulation via the HOG pathway even under non-stress conditions [36]. Additionally, the Gpd2 OUT, Gut1 OUT, Sko1 OUT and Lys20 OUT mutants associated with themselves across conditions. These mutants are cases without strong effects under any tested condition, but with subtle condition-independent changes, which suggest a possible fine-tuning function of the phosphosites. The respective metabolites driving condition-independent correlations were diverse, e.g. affecting nucleotide metabolism in Gut1 OUT and lysine in Lys20 OUT. As expected, we also recovered Pda1 IN that has effects throughout all conditions.

To systematically interpret the evidence for functionality of a given mutated phosphosite, we ranked the functionality evidence based on a score that reflects the number of changing ions, enriched pathways, correlation with other mutants, and correlation with itself for each phosphosite (Fig. 4). The number of changing ions was scored between 0 and 3. The changes were classified as either little (1–4 ions), affecting equal or less ions than in the wild-type, where 4.75 ions were changing on average, intermediate (5–19 ions), or many (>19 ions) changes, affecting a broad part of metabolism with changes in > 10% of all averagely detected ions. The growth condition with the most changes was considered for scoring. The three other aspects were scored between 0 and 2 according to p-values. A score was attributed based on the most significant p-value obtained in any growth condition. A p-value of at least 10⁻² was scored with a 1, while a p-value below 10⁻³ was scored with a 2. Supporting our approach, most of the already characterized phosphosites scored very high (Fig. 4). The three characterized cases with low scores can be explained by incomplete mutation of all relevant phosphosites (e.g. Ste50) or lack of protein relevance under the tested condition (e.g. Crz1 and Put3) [18, 37, 38]. Our data suggests that the reported phosphoregulation of Pfk2 at position S163 [11] resulted, at least in part, from a non-endogenous increase of Pfk2 levels from the employed plasmid construct (Additional file 2: Figure S6). Since we cannot exclude that the reported phenotype of this Pfk2 variant originated from increased protein abundance, we did not consider Pfk2 S163 as a known functional phosphosite and found only weak functionality evidence in the Pfk2 IN mutant. Overall, the metabolomics data identified most of the already characterized phosphosites as functional with a total score of 3 or higher, one novel regulatory phosphosite with high confidence (Tip41 S79), and several other novel functional phosphosites with lower confidence. The much higher number of phosphosite functionality identification by metabolomics compared to growth phenotype analysis underlines the value of metabolomics in assessing at least metabolic functionality. In the next sections, we exploit the metabolomics data to hypothesize specific functions for phosphosites in the glycerol dehydrogenases Gpd1 and Gpd2, and the MAK kinase Pbs2.

Fig. 4

Phosphomutant functionality scoring based on metabolic profiles. The number of changing ions was scored between 0 and 3 according to the following scheme: 0 for no change; 1 for 1–4 changes; 2 for 5–19 changes; 3 for > 19 changes. The remaining three criteria were scored between 0 and 2 based to on the strongest p-value: 0 for p-value > 10⁻²; 1 for p-value < 10⁻²; 2 for p-value < 10⁻³. Mutants exhibiting a growth phenotype in the growth analysis are highlighted in blue

Phosphorylation of Pbs2 at S508 negatively regulates Hog1 activity

Phosphorylation at position S514 and T518 activates the MAPK kinase Pbs2, which in turn activates the Hog1 kinase in response to osmotic stress [28, 39] (Additional file 2: Figure S7). Besides these two well-studies phosphosites, the here investigated S508 phosphosite was suggested to be important for the inactivation of Pbs2 after adaptation to osmotic stress [15]. Our approach scored both Pbs2 IN and its target Hog1 OUT as regulated by phosphorylation (Fig. 4). The strongest metabolic responses of both mutants occurred during growth on pyruvate, presumably due to the osmotic stress of increased sodium concentration in this medium. The striking similarity of the Hog1 OUT and Pbs2 IN profiles and the very different Pbs2 OUT profiles strongly support the hypothesis that phosphorylation at position S508 deactivates Pbs2, which causes reduced Hog1 activity (Fig. 5a and b).

Fig. 5

Metabolic profiles, HOG signaling activation and phenotype of Pbs2 phosphomutants. a Log2 fold changes of all ions corresponding to deprotonated metabolites of Hog OUT compared to Pbs2 IN in pyruvate/NH₄ ⁺ media. b Log2 fold changes of all ions corresponding to deprotonated metabolites of Hog OUT compared to Pbs2 OUT in pyruvate/NH₄ ⁺ media. c Time course Western blot of Hog1 activity in wild-type and Pbs2 phosphomutants, measured as phosphorylated Hog1 (P-Hog1) upon osmotic stress. Total Hog1 (t-Hog1) was measured as control. d Phenotype analysis upon osmotic stress and exposure to different temperatures. Spotting assay of Pbs2 phosphomutants and Hog1 OUT in synthetic complete (SC)-glucose, SC-glucose supplemented with 0.5 M NaCl, SC-glucose at 42 °C and SC-glucose at 16 °C. The assay was done three times, the results of one exemplary experiment are shown

We further investigated the presumably inactivating role of S508 by exposing the Pbs2 mutants to elevated concentrations of NaCl and monitored HOG pathway activity by following Hog1 phosphorylation over time (Fig. 5c). Consistent with the above reasoning, Pbs2 IN failed to phosphorylate and activate Hog1 upon osmotic stress. However, no delay in inactivation of Hog1 was observed for Pbs2 OUT compared to wild-type dynamics, suggesting that deactivation of Pbs2 was not impaired in the Pbs2 OUT mutant. Since the HOG signaling pathway is known to cross-talk with other signaling pathways, we wondered whether phosphorylation of S508 had other physiological roles besides the deactivation of Pbs2 after adaptation to osmotic stress. One such interconnected pathway is the cell wall integrity pathway [40, 41]. To test a link of Pbs2 S508 phosphorylation to cell wall integrity, we grew the mutants at different temperatures to impose different cell wall stress levels. We found that growth of Pbs2 OUT was exclusively impaired at elevated temperatures, while Pbs2 IN and Hog1 OUT showed growth defect most strongly when growing at 16 °C (Fig. 5d). A change in cell wall stress leading to conditional phenotypes of the Pbs2 phosphomutants suggests that phosphorylation of S508 is important for proper activation of HOG signaling under specific cell wall stress conditions and plays a critical role in the cross-talk between the two signaling pathways.

Phosphorylation of Gpd1 and Gpd2 influences arginine metabolism

The glycerol dehydrogenases Gpd1 and Gpd2 catalyze the conversion of dihydroxyacetone phosphate to glycerol-3-phosphate, the first step of glycerol synthesis. Although Gpd1 and Gpd2 have distinct roles, they can at least partially compensate for each other [42, 43]. In addition to transcriptional regulation, both enzymes are subject to phosphoregulation [11, 20]. In our screen, both glycerol dehydrogenases were detected as regulated phosphoproteins, with Gpd2 OUT ranking slightly higher in the metabolic profile analysis than Gpd1 OUT (Fig. 4). Although we did not detect changes in the direct reactants of either of the mutated isoenzymes, both OUT mutants exhibited specific metabolic responses throughout most conditions with an enrichment in arginine and proline metabolism (p-value 0.002), in particular when growing on proline as nitrogen source (Additional file 4). Unexpectedly, we detected the strongest changes in citrulline and ornithine, along with some less pronounced changes in arginine that we confirmed by targeted LC-MS/MS measurements (Fig. 6a, b). The striking correlation of Gpd2 OUT with itself across all conditions demonstrates that this phosphoregulation is not specific to the proline condition (Additional file 6).

Fig. 6

Metabolite responses in arginine biosynthesis of the Gpd1 OUT and Gpd2 OUT mutants. a Metabolite response in arginine biosynthesis under all tested conditions. Coloring indicates strength and direction of fold changes. Metabolites in grey were not detected in the respective conditions. Metabolites without color boxes were not detected. b Quantification of arginine and citrulline levels by LC-MS/MS for Gpd1 OUT and wild-type in glucose/NH₄ ⁺ medium. c Possible metabolic routes for G3P

How glycerol metabolism influences arginine metabolism remains unclear. The OUT mutants have more active glycerol dehydrogenases and hence generate more G3P than the wild-type [11, 20], but neither of the known routes from G3P connects to arginine metabolism (Fig. 6c). Although the observed metabolite changes suggest a possible blockage at Arg1 or Arg4 in the glycerol dehydrogenases phosphomutants (Fig. 6a), the influence is more likely to occur indirectly via a regulatory interaction. Previously, transcriptional and translational regulation of arginine metabolism has been suggested to be implicated in adaptation to osmotic stress [44]. In our study, Pbs2 IN and Hog1 OUT, which have impaired HOG signaling and cannot activate Gpd1 by dephosphorylation, showed strong changes mostly in arginine but also citrulline upon osmotic stress (log2 fold change arginine of 1.845 and 2.463 for Pbs2 IN and Hog1 OUT respectively, see Additional files 7 and 8). These metabolites change in opposite directions compared to the Gdp1 OUT and Gpd2 OUT mutants where arginine is decreased while citrulline is increased, supporting a phosphorylation-dependent influence of the glycerol dehydrogenases on arginine metabolism.

Discussion

The large networks of kinases and phosphatases represent possibly the most complex post-translational regulation system and play key roles in essentially all cellular processes. Given the rapidly growing collection of reported phosphosites, a major limitation is assessing phosphosite functionality and generating hypotheses on their precise roles, at a scale and speed that matches the high-throughput of modern mapping methods [2, 3, 5, 8, 12, 45]. Here we demonstrate that untargeted metabolomics of point mutants with phospho IN or OUT modifications has the potential to address this functionality problem. From phosphosite mutants of 26 proteins involved in metabolism, TOR and HOG signaling as a proof-of-concept, we were able to recover 10 of the 13 positive controls by metabolomics, compared to only three when relying on quantitative growth assays. Importantly, metabolomics was able to detect also more subtle metabolic effects that were phenotypically silent. The other three known phosphoproteins were presumably not detected because we did not mutate all relevant phosphosites (Ste50 OUT) or did not chose the appropriate conditions that required activity of these proteins (Crz1 OUT and Put3 OUT) [18, 37, 38].

Of the remaining 13 so far uncharacterized phosphosites, metabolomics provides evidence that those on Cys3, Gln3, Gut1, Nbp2, Pfk2, Tco89 and Tip41 are generally functional. At least Gln3 and Tip41 had been suspected to be regulated by phosphorylation, although no phosphosites were identified so far [46–48]. The strongest evidence in our data was for S79 on Tip41, an essential regulatory protein in TOR signaling whose role upon activation of TOR is not yet clearly established [29, 49]. The exclusively strong effect of Tip41 OUT during growth on a poor nitrogen source suggests that S79 is phosphorylated when the TOR complex is inactive and becomes dephosphorylated upon activation of TOR. Since our recovery of positive controls was fairly good, we believe that the newly tested phosphosites for which we could not obtain functionality evidence are indeed mostly without functional relevance. The employed, metabolomics-based scoring system allowed classification of functionality, but does not directly identify the regulatory consequences of phosphorylation. Nevertheless, in several cases the metabolite profiles provided starting hypotheses for the specific biological function without prior knowledge. For instance accumulation of pathway end products suggested that phosphorylation of Lys20, the first enzyme in lysine biosynthesis, specifically regulates lysine pathway flux. The correlation of Pbs2 IN and Hog1 OUT metabolite profiles strongly suggested that the two phosphorylation events have a similar function in mediating HOG signal flow. However, changes in metabolism were not always intuitive and possibly not direct consequences of the mutated phosphosites, as observed in the Gpd1 and Gpd2 mutants. The main benefit of the metabolomics data in the process of understanding phosphosite functionality is to point out relevant cases efficiently, which can then be investigated in more depth by other methods or in combination with proteomics and transcriptomics data.

Even though about half of the tested phosphosites appeared to be functional, most phosphomutants had no detectable growth phenotype and exhibited relatively modest metabolic changes, both in number and magnitude. This apparent resilience to point mutations does not come entirely unexpected given that the much more drastic intervention of 118 kinase and phosphatase deletions in S. cerevisiae resulted in few detectable growth phenotypes where only 26% of the mutants exhibited a significant transcriptional response and about 40% had no detectable metabolite response [26, 50]. Several reasons are likely to contribute to this robustness, including the large overlap and redundancy both in kinase targets and phosphosites on a given protein [2, 51–53]. Additionally, adaptation to a new steady state in our experiments will mask many direct regulatory effects and our population-based measurements will have missed, for example, cell cycle specific regulation events.

Conclusions

In this study, we investigated the potential of untargeted metabolomics along with growth analysis to identify functionality of phosphosites in a higher throughput. While growth analysis was suitable to detect functionality of a small number of phosphosites, metabolomics allowed recovering most of the previously characterized sites, along with some novel functional phosphosites and aided in hypothesis generation on biological roles of functionally relevant sites. Our results demonstrate the potential of untargeted metabolomics to experimentally identify functionality and conditional relevance of post-translational modifications at high throughput to complement current computational approaches [2, 52–55].

Methods

Strains used in this study

All phosphosite mutants were generated in S. cerevisiae wild-type FY4 using the delitto perfetto method [21], the oligonucleotides used in this study are listed in Additional file 9. Briefly, in a first step, a cassette containing a positive and a negative selection marker for growth on geneticin or galactose containing media was amplified from pCORE-Kp53 using primers with flanking region homologous to the genomic region of the phosphosite of interest. The cassette was inserted via homologous recombination into the region of interest into wild-type FY4. Next, the entire cassette was replaced via homologous recombination with oligonucleotides containing the desired point mutation, resulting in genomic point mutants of candidate phosphosites free of selection markers. Mutants were selected after confirmation of proper integration of the desired mutation by colony PCR and sequencing of the manipulated genomic region +/− 200 base pairs. Up to three independent mutant replicates per phosphosite were generated. Since no differences in the metabolic profiles were observed between replicates of phosphomutants with strong metabolic effects, the genetic approach was considered reliable and data from untargeted metabolomics and growth phenotype analysis of the replicates was merged for analysis.

Growth phenotype analysis

All strains were grown in phthalate-buffered Verduyn minimal media supplemented with different carbon and nitrogen sources [56]: glucose/NH₄ ⁺ media contained 20 g l⁻¹ glucose and 5 g l⁻¹ (NH₄)₂SO₄, pyruvate/NH₄ ⁺ media contained 24.42 g l⁻¹ sodium pyruvate and 5 g l⁻¹ (NH₄)₂SO₄, glucose/glutamine media contained 20 g l⁻¹ glucose, 5 g l⁻¹ K₂SO₄ and 1 g l⁻¹ glutamine and glucose/proline media contained 20 g l⁻¹ glucose, 5 g l⁻¹ K₂SO₄ and 1 g l⁻¹ proline.

Precultures with biological triplicates of each strain were grown at 30 °C, 300 rpm, in deep 96-well plates in 1 ml of media, using one glass bead (Ø 4 mm) per well to improve mixing efficiency. The preculture was used for inoculation of microtiter plates at a starting optical density of 595 nm (OD₅₉₅) of 0.03. Growth was monitored via light scattering using the Biolector® plate reader at 30 °C, 1000 rpm. All experiments were run in duplicates. The growth rate for all replicates was determined by a linear fit to the log-transformed OD₅₉₅ growth curve. Maximum cell densities were determined as the end point of growth at the entry of the stationary phase. Clear outliers such as non-growing replicates and single replicates with different growth signal e.g. due to technical reasons were removed manually. A two-sided 2-sample t-test assuming unequal variance was performed to compare the significance of difference in growth rate and maximally reached cell density between mutants and wild-type. The p-values were corrected for multiple testing as described by Storey [57]. Growth rates or maximally reached cell densities of at least 10% difference to the wild-type and a p-value < 10⁻³ were considered significant. Abnormal lag time and growth curve shape were qualitatively assigned.

Metabolite extraction, sample measurement and data processing for untargeted metabolomics

Metabolites were extracted from deep 96-well cultures when the average OD₅₉₅ of the plate was around 1 and all cultures had doubled at least twice. The entire plates was centrifuged for 1 min at −9 °C, 4000 rpm. The supernatant was discarded, 150 ul of 80 °C 75% (v/v) ethanol in 10 mM ammonium acetate at pH 7.5 was added to each well. For extraction, plates were incubated for 5 min at 80 °C with three vortexing steps. The plates were centrifuged again for 5 min and the supernatant containing the metabolites was transferred to fresh plates for analysis. The metabolite extraction and measurements of extracts were conducted in two batches of mutants (Additional file 1). Mutants were grown in triplicates on the 96-well plate while wild-type was inoculated in four to five replicates per plate. Each culture plate was grown and extracted twice on different days yielding to day replicates for every mutant.

Metabolites were measured by direct flow double injection of extracts on an Agilent 6550 series quadrupole TOF MS with the help of a GERSTEL MPS2. Ions within a mass/charge ratio range of 50–1000 were measured [33]. The ions were annotated either using the genome-wide reconstruction model of S. cerevisiae metabolism by Dobson et al. (2010) or, for pathway enrichment analysis, the S. cerevisiae reactants defined in the KEGG database [34, 58]. Only ions corresponding to deprotonated metabolites were considered for further analysis. The entire annotation is provided in Additional file 10. All data analysis was conducted using Matlab (The Mathworks, Natick).

Measured ion intensities of samples taken at the exponential phase were corrected for potential signal drifts during the measurement run by applying locally weighted scatter plot smoothing (LOWESS). For further processing and retrieval of changes in metabolites only samples in a restricted OD₅₉₅ range between 0.6 - 1.7, with marginal variation between conditions and sampling batches, were considered. This OD range contained the majority of all exponentially growing samples. A fitted intensity curve was used as OD-wise reference curve to compare each ion of every sample to the average intensity at its specific sampling OD. This OD-wise comparison allows to account for potential impacts of the extraction OD on ion intensity and comparison of samples with different extraction ODs. Fitting was done by smoothing the measured intensities of all samples vs. the OD₅₉₅ of the samples using LOWESS. Including all samples for the calculation of the reference curve enables a thorough covering of the entire OD range and renders the reference curve robust. Further, LOWESS was performed in robust mode to prevent individual strong metabolic changes from having a dominating impact on the reference curve. An adapted procedure for the calculation of the OD-wise reference curve was chosen for batch 2, as the number of Pda1 OUT and Pda1 IN mutant strains strongly dominated this batch, being the only mutants generated in three replicates in batch 2. Therefore, three Pda1 OUT and three Pda1 IN replicates were chosen randomly from each of the two deep 96-well plate replicate for the calculation of the reference curve, thereby adapting the number of replicates for this mutant to the remaining phosphomutants in batch 2. Random sample choice was done 1000 times, and the final reference curve for batch 2 was calculated as the median of the references curves of all permutations.

For every ion we determined the fold change compared to the intensity of the ion-specific reference curve at the corresponding sample-OD. Fold changes for all samples of each mutant were pooled and ion-wise outliers that deviate more than two standard deviations from the mean fold change were excluded. The median log2 fold change was determined and a p-value was assigned by comparing the fold changes of one mutant to all other samples using a two-sided 2-sample t-test assuming unequal variance (Additional files 7 and 8). The p-values were corrected for multiple testing as described before [57]. A cutoff on the fold change was applied for removal of very small changes. For this purpose, the absolute log2 fold changes of all wild-type ions were pooled, and the 97.5% quantile was determined to be at 0.381. Accordingly, corresponding to a 30% change, the log2 fold change cutoff was defined at 0.3785. Ions with a log2 fold change stronger than +/−0.3785 and a corrected p-value < 10^-3 were considered significantly changing.

Analysis of untargeted metabolomics data

Pathway enrichment analysis was performed based on the approach detailed in Kühne et al. (2015) [59]. Briefly, metabolites passing a relaxed cutoff (p-value 0.1 and absolute log2 fold change of 0.1375, corresponding to a 10% change) were sorted according to their absolute change, and a hypergeometric test was performed to calculate the enrichment of each metabolic pathway. For the enrichment calculation, the hit subset was defined starting with only the top hit of the generated list of sorted metabolites. For each subsequent enrichment analysis, the hit subset was extended by adding the next metabolite from the sorted list. Eventually the most significant result of all generated enrichments based on the different hit subsets was considered. Pathway definition was retrieved from KEGG [58]. Ions that were annotated to several metabolites within a given metabolic pathway were lumped to one entry.

To retrieve the mutant-pair correlation over all conditions, a pair-wise Pearson correlation coefficient was calculated for all mutant-pairs for each condition. Calculation was based on log2 fold changes of ions corresponding to deprotonated metabolites. The correlation coefficients of all mutant-pairs showed a normal distribution that centered around 0 (Additional file 2: Figure S5A). The four correlation coefficients from each mutant-pair were compared using a 2-sample t-test against this distribution of correlation coefficients of all remaining mutant-pairs in order to find meaningful mutant-pair correlations. Mutant-pairs with a p-value below 0.01 were considered as correlating significantly with each other throughout all conditions. For analyzing self-correlation of mutants across conditions, the pair-wise Pearson correlation coefficient was calculated for each mutant for every condition-pair combination, resulting in six correlation coefficients per phosphomutant. Calculation was based on log2 fold changes of ions corresponding to deprotonated metabolites. The six correlation coefficients of each mutant were compared against an average of zero performing a t-test. Mutants with a p-value below 0.01 were considered as correlating significantly with themselves across conditions.

Targeted LC-MS/MS analysis

For targeted LC-MS/MS analysis, 500 ml shake flasks containing 50 ml glucose/NH₄ ⁺ medium were inoculated with overnight precultures grown in the same medium. At OD₅₉₅ 1 +/− 0.15, 1 ml of culture was quenched by mixing with 4 ml 60% (v/v) methanol precooled at -40 °C. 100 ul uniformly-labeled ¹³C yeast extract was added as internal standard and metabolites from cell pellets were extracted by incubation with 1 ml 75% (v/v) 80 °C ethanol for 3 min, with three vortexing steps. Samples were dried and resuspended in MilliQ water and analyzed for quantitative targeted analysis of specific metabolites as described before with the exception that citrulline was normalized to the ¹³C-arginine signal for lack of ¹³C-citrulline signal in the standard [60]. Mean and standard deviation of three wild-type and four Gpd1 OUT replicates was determined.

Spotting assay

Yeast cells from SC-galactose plates were dissolved in 0.9% (w/v) NaCl and OD₅₉₅ was adjusted to 1. Serial 1:10 dilutions were prepared and 3.5 ul were spotted onto SC-glucose plates with or without 0.5 M NaCl. Plates without NaCl were incubated at 16 °C, 30 °C or 42 °C, while plates with NaCl were incubated at 30 °C for 2–4 days.

Hog1 phosphorylation analysis in the Pbs2 mutants

Yeast strains were grown in YPD at 30 °C to an OD₆₆₀ of 0.6. Cells were stressed with NaCl (final concentration 0.4 M) and samples were taken at the indicated time points. Proteins were extracted by glass bead lysis, separated on a 10% (w/v) SDS page and transferred to a PVDF membrane. The membrane was probed with anti-ph-p38 (Cell Signaling 9215 L) and anti-total Hog1 (Santa Cruz Biotechnology sc6815).