Mycoplasma pneumoniae, a threatening pathogen with a minimal genome, is a model organism for bacterial systems biology for which substantial experimental information is available. With the goal of understanding the complex interactions underlying its metabolism, we analyzed and characterized the metabolic network of M. pneumoniae in great detail, integrating data from different omics analyses under a range of conditions into a constraint‐based model backbone. Iterating model predictions, hypothesis generation, experimental testing, and model refinement, we accurately curated the network and quantitatively explored the energy metabolism. In contrast to other bacteria, M. pneumoniae uses most of its energy for maintenance tasks instead of growth. We show that in highly linear networks the prediction of flux distributions for different growth times allows analysis of time‐dependent changes, albeit using a static model. By performing an in silico knock‐out study as well as analyzing flux distributions in single and double mutant phenotypes, we demonstrated that the model accurately represents the metabolism of M. pneumoniae. The experimentally validated model provides a solid basis for understanding its metabolic regulatory mechanisms.
A new genome‐scale metabolic reconstruction of M. pneumonia is used in combination with external metabolite measurement and protein abundance measurements to quantitatively explore the energy metabolism of this genome‐reduce human pathogen.
We established a detailed biomass composition for M. pneumoniae, thus allowing for growth simulations.
Using our metabolic model, we corrected the metabolic network topology and the functional annotation of key metabolic enzymes.
M. pneumoniae, unlike other laboratory‐grown bacteria, uses a high fraction of energy (up to 89%) for cellular maintenance and not for growth.
Simulating different growth conditions as well as single and double mutant phenotypes, we analyzed pathway connectivity and the impact of gene deletions on the growth performance of M. pneumoniae, highlighting the limited adaptive capabilities of this minimal model organism.
Representing cellular networks by mathematical models and gaining novel biological insights from subsequent in silico analysis and iterative experimental validation are major hallmarks of systems biology. The integration of experimental data from different sources into suitable mathematical models poses a formidable challenge, but allows to place metabolites and enzymes into their network context and to extract biologically relevant information for the examined system (Kitano, 2002). For this purpose, constraint‐based modeling approaches that allow the determination of possible network flux distributions provide a useful framework (Fell and Small, 1986; Savinell and Palsson, 1992a, 1992b; Feist et al, 2009; Oberhardt et al, 2009).
Flux balance analysis (FBA) is a mathematical method to determine metabolic fluxes within a constraint‐based model, that is, fulfilling the steady‐state condition. Thereby, the flux distribution is optimized toward so‐called objective functions, commonly energy production or growth in form of biomass production, for a given set of available nutrients (Varma and Palsson, 1994a; Kauffman et al, 2003; Reed and Palsson, 2003). FBA has been successfully used to, for example, predict the effect of mutations on Escherichia coli growth rates (Edwards and Palsson, 2000), to predict active pathways under different growth conditions (Covert et al, 2001), to improve biotechnology applications (Puchalka et al, 2008), and to understand infection processes (Oberhardt et al, 2008). More recently, in vivo reaction directionalities were assigned to the E. coli metabolic network, considering quantitative metabolite measurements and associated fluxes (Fleming et al, 2009). In the past year, several tools for constraint‐based modeling, including tools for data integration and extraction of relevant information from the modelled system, have been developed (Fleming and Thiele, 2011; Schellenberger et al, 2011; Thorleifsson and Thiele, 2011). In addition, the first reconstruction of the human metabolism has been published (Duarte et al, 2007; Rolfsson et al, 2011) being the largest metabolic reconstruction so far.
However, the sheer complexity of living systems and the scarce availability of unbiased, large‐scale quantitative data often prevent comprehensive organism‐wide studies. Thus, knowledge gain from modeling approaches is limited due to the difficulties in setting up validation experiments for specific model predictions. In recent years, the bacterium Mycoplasma pneumoniae, a human pathogen preferentially colonizing the pulmonary epithelium (Waites and Talkington, 2004), has been established as model organism for systems biology, providing the research community with detailed quantitative information on its genome, transcriptome, proteome, and metabolism (Güell et al, 2009; Kühner et al, 2009; Yus et al, 2009; Maier et al, 2011). M. pneumoniae combines several favorable properties for organism‐wide analyses. It evolved by massive genome reduction, resulting in a single, small chromosome of only 816 394 base pairs, encoding 689 proteins (Himmelreich et al, 1996; Dandekar et al, 2000). Its expressed proteome is of low complexity, spanning only three orders of magnitude in abundance (Maier et al, 2011). As a consequence of reduced genome and parasitic life, M. pneumoniae lacks many metabolic pathways, forcing it to acquire the necessary cell building blocks, among them amino acids, nucleobases, and fatty acids, from the environment (Yus et al, 2009). For ATP generation, M. pneumoniae relies on simple organic acid fermentation, due to the absence of TCA cycle and a functional respiratory chain (Pollack et al, 1997; Yus et al, 2009). The lack of most anabolic processes and rescue pathways known from more complex organisms and, consequently, the expected high linearity of its metabolic network render M. pneumoniae an ideal organism to study basic metabolic functions and to dissect energetic expenses. In addition, it can be grown autonomously in a laboratory environment and basic genetic tools, such as transposon mutagenesis to study gene essentiality, are available.
We present here a genome‐scale constraint‐based model of the M. pneumoniae metabolism (iJW145—in silico model including 145 genes, designed by JW). We established the biomass composition of an average M. pneumoniae cell based on quantitative experimental data. Performing FBA for growth simulations under a wide variety of conditions, we explored and characterized the metabolic network by an iterative cycle between model predictions and their experimental validation. We corrected the metabolic network annotation and the functional annotation of three metabolic enzymes. In addition, we quantitatively dissected the M. pneumoniae energy metabolism. We show that M. pneumoniae dedicates most of its energy to cellular homeostasis and not growth, as suggested for other bacteria (Schuetz et al, 2012), possibly as a consequence of its small size and its parasitic life style. We analyzed carbon fluxes in vivo and predicted metabolic fluxes in silico. Finally, we proved the predictive capacity of the model by predicting gene essentiality with high accuracy (95%) and specificity (98%). The analysis of the derived single and double mutant phenotypes provided insight into the adaptive capabilities of the metabolic network.
The presented model and the biological findings herein reported can be used for the design of future experiments and the development of novel engineering tools. Furthermore, this model will provide the basis to design dynamic models for metabolic sub‐networks and to relate the gene‐regulatory and the protein interaction networks to metabolism in M. pneumoniae.
A constraint‐based metabolic reconstruction is a union of (i) a stoichiometric metabolic model, (ii) a set of constraints for metabolic fluxes. and (iii) a list of genes responsible for the catalysis of reactions included in the model. Based on a curated wiring diagram for the metabolism of the bacterium M. pneumoniae (Yus et al, 2009), we built a genome‐scale constraint‐based metabolic reconstruction (iJW145_reconstruct—initial reconstruction of the network, Figure 1). The main subsystems of the metabolic network are energy producing pathways, amino acid, nucleotide, lipid, and cofactor metabolism as well as transport reactions. In addition, we considered RNA, DNA, and protein biosynthesis as part of the metabolic network (Supplementary Figure S1; Supplementary Table S1).
Missing transport reactions were added and reaction reversibilities for the metabolic network were defined (Supplementary Table S2; Supplementary information). We assessed the correct connectivity of the reconstructed network by individually testing the feasibility of the production of all metabolic intermediates in silico. Contrary to more complex organisms, the M. pneumoniae metabolic network is predominantly composed of linear pathway modules and with the exception of ubiquitous cofactors, such as AMP, ADP, ATP, H+, H2O, NAD+, NADH, Pi, PPi, only few metabolites interconnect the core metabolic routes (Supplementary Figure S1; Supplementary Table S3). For our study, this low connectivity of the M. pneumoniae metabolism facilitates the analysis of inter‐pathway crosstalk and enables the prediction of metabolic changes in response to environmental perturbations. Particularly, the absence of most catabolic and also a high number of anabolic routes in M. pneumoniae allows the direct relation of external metabolite measurements to intracellular metabolic fluxes and corresponding catabolic activity.
A possible connection of sugar metabolism to aspartic acid via its direct precursor oxaloacetate, which can be synthesized from pyruvate or malate, has been suggested (Manolukas et al, 1988). However, in an organism‐wide approach to quantify cellular proteins of M. pneumoniae, no proof for the presence of the proteins catalyzing the corresponding reactions has been found (Maier et al, 2011). To assure the completeness of the reconstructed network, we validated the proposed connection in vivo by determining the carbon flux from glycolysis into aspartate. To this end, we monitored the incorporation of 13C‐labelled carbon into aspartate. In agreement with previous results (Yus et al, 2009), even after 96 h of growth in medium containing 13C6‐labelled glucose as sole carbon source, no increase in heavy isotope‐labelled reporter ions for aspartate was observed (Supplementary Figure S2A). Thus, for M. pneumoniae cells grown in rich medium we could discard a link between glycolysis and amino‐acid metabolism involving aspartate as connecting metabolite.
To be able to predict metabolic flux distributions and to simulate cell growth with FBA, an accurate, quantitative representation of the biomass composition of an average M. pneumoniae cell is required. Therefore, experimental data on DNA content, RNA composition, and protein abundances were considered (Table I; Supplementary information) (Güell et al, 2009, 2011; Maier et al, 2011). For all three classes of macromolecules, we defined an artificial molecule reflecting their average cellular composition (Supplementary information). To accurately represent the composition of the cell membrane in our model, we experimentally analyzed the fatty acid profile of M. pneumoniae and quantified the fatty acid composition of the cytosol, the surrounding growth medium, as well as the cytoplasmic membrane directly (Supplementary Figure S3; Supplementary information). We found predominantly fatty acid chains of 16 and 18 carbons length in both the membrane and the cytoplasm (Supplementary Figure S3A). Based on our results, we designed an artificial molecule representing the average fatty acid composition to describe lipids in the M. pneumoniae biomass equation (Supplementary information). To account for the experimentally shown essentiality of vitamins and other cofactors, such as FAD, NAD+/NADH or folate (Yus et al, 2009), the end products of the secondary metabolism pathways were included into the biomass qualitatively, that is, in small arbitrary quantities, since a sensitivity analysis showed that a 10‐fold change in cofactor amounts does not significantly influence on the growth rate (Supplementary Figure S4; Supplementary Table S4; Supplementary information). Other cofactors, for example organic phosphate, were included based on the literature information (Table I; Supplementary information).
The biomass equation defining the macromolecular composition of an average M. pneumoniae cell in a general form reads:
Equation 1: Biomass (eq. 1)
DNA+RNA+proteins+lipids+bases+amino acids+fatty acids+cofactors→Biomass
The assembly of the stoichiometric network (iJW145_reconstruct) and the assignment of reaction reversibilities together with the definition of the biomass composition provided us with a version of our M. pneumoniae model that was able to simulate growth: iJW145_growth (Figure 1).
Model assessment and refinement
To validate the network structure of iJW145_growth and to avoid errors in the reconstruction (Reed et al, 2006), we determined the network behavior under different nutrition conditions. To this end, we defined several flux constraint sets based on own experimental data (Supplementary Figure S5) and the literature information (Supplementary Table S5; Supplementary information). As a general strategy, we minimized the possible number of flux constraints to not restrict the solution space of the model but to keep a high predictive capacity (Edwards et al, 2002; Covert and Palsson, 2003; Price et al, 2004). Maximum constraints were set to limit the uptake of all sugar sources and arginine (reflecting the available nutrients for the growth conditions and preventing unlimited ATP synthesis from arginine, Supplementary information), while minimum constraints account for mRNA and protein degradation rates, detoxification, and maintenance energy costs (Supplementary Table S5; Supplementary information).
Initial simulations for different growth conditions verified the network annotation (Supplementary Figure S1; Supplementary information). In addition, the observed in silico fluxes confirmed experimental data on mRNA expression of metabolic proteins (Güell et al, 2009) and for the design of a defined medium (Yus et al, 2009) (Supplementary information). However, we also identified several conflicts between model predictions and available experimental results (Hames et al, 2009; Yus et al, 2009). These conflicts were resolved by an iterative process, involving protein sequence comparison, additional experiments, literature mining, and repeated simulations using adjusted constraints (Figure 2; Supplementary information). In summary, the comparison of model results and experimental data referring to cellular redox states and to nucleotide metabolism led to the functional re‐annotation of three metabolic enzymes: MPN051—glycerol phosphate oxidase (Figure 2A), MPN394—NADH oxidase (Figure 2B), MPN256—no CTP synthase activity (Figure 2C), and guided the correction of the wiring diagram (correction of one reaction, deletion of two reactions, Supplementary information).
Simulation of biomass production
Batch culture growth experiments on glucose lead to a measurable decrease in medium pH (pH 7.8–5.5 during a 4‐day growth course) due to secretion of lactic and acetic acid (Yus et al, 2009). A shift from energetically favorable acetic acid production during early growth (4 ATP/glucose) toward predominantly lactic acid production during later growth stages (two ATP/glucose) has been observed (Supplementary Figure S5C–E). Concurrently, cellular lactate dehydrogenase levels (MPN674, LDH) increased five‐fold from 203 copies to above 1000 copies per cell during 4 days of growth (Maier et al, 2011). To correctly represent this metabolic shift in our model and to account for further carbon incorporated into glycolysis from other sources than glucose, we directly constrained the energetically favored acetic acid production. To this end, we fitted the maximum constraint for acetate production to the experimentally determined ratio of secreted lactic and acetic acid (Supplementary Figure S5E).
In laboratory experiments, M. pneumoniae can be grown in rich and defined medium (Chanock et al, 1962; Yus et al, 2009). In both conditions, it metabolizes glucose as major carbon and energy source. Alternatively, a variety of additional reduced carbon compounds, such as fructose, mannose, ribose, ascorbic acid, glycerol, glycerol 3‐phosphate (G3P) (Yus et al, 2009), and glycerophosphocholine (G3PC) can be metabolized (Schmidl et al, 2011). In addition, arginine can be used to produce ATP but we did not consider it as alternative carbon source since (i) its contribution to energy production is negligible, (ii) in vivo it did not permit growth (Yus et al, 2009), and (iii) only one enzyme involved in arginine fermentation has been detected (Maier et al, 2011). For the other carbon compounds, we adjusted the flux constraints of the model according to rich and defined growth medium compositions as described in Yus et al (2009) and applied FBA. We verified in silico the experimentally tested metabolic capabilities of M. pneumoniae (Table II; Yus et al, 2009; Schmidl et al, 2011). Under defined medium conditions (Supplementary Table S5) the model predicted growth on glucose, mannose, mannitol, ribose, and ascorbate. However, these strict constraints did not allow growth on fructose, glycerol, G3P, and G3PC (Table II). We identified the inability to provide the pentose phosphate pathway precursor fructose 6‐phosphate (F6P) for de novo nucleotide synthesis as cause for the observed in silico growth limitations (Supplementary Figure S1).
While in silico there is no difference in energy yield for M. pneumoniae grown on different carbon sources (assuming the same amount of carbon taken up for each sugar source, Supplementary information), in vivo the doubling times differed significantly and only glucose and mannose allowed robust growth (Yus et al, 2009). This apparent discrepancy can be explained by the abundances of the respective uptake and processing systems. While the glucose‐specific uptake protein (MPN207) has high copy numbers (∼385/cell), the known transport proteins for other sugars (fructose, ribose, ascorbate, mannitol, mannose, glycerol, and G3P) are 8–20 times less abundant (Maier et al, 2011). In fact, to allow relevant growth of M. pneumoniae on fructose in vivo the cells had to be adapted over several serial passages and showed significant overexpression of the proteins involved in fructose import and metabolism (Yus et al, 2009).
Cellular energy balance
Using constraints derived from non‐linear fittings to experimentally obtained metabolite quantification data (glucose consumption as well as acetic and lactic acid production for defining the carbon uptake and the acetic acid production, and total protein increase to restrict the growth rate) from batch culture growth (Supplementary Figure S5B–E), we applied FBA to our refined model to determine cellular ATP production at different growth stages (24, 36, 48, and 60 h). We found that during mid exponential growth (36 h after inoculation) ∼60 000 ATP molecules per cell and second are synthesized in silico (Supplementary Table S6; Supplementary Materials and methods).
Next, we quantitatively assessed the contribution of the available cellular energy to biomass production (eq. 1, Biomass) and other functions specified in the model, for example, protein turnover (Supplementary information). For this purpose, we compared in silico and in vivo doubling times (Figure 3A) at different growth stages. We found that simulation‐based results (between 2.3 and 3.8 h) differed significantly from experimentally measured values during the exponential phase in batch culture growth (between 19.7 and 59.7 h; Figure 3A), as well as from previous reports on microscope studies of single cell division (∼8 h; Seybert et al, 2006). These results suggest the existence of yet undefined ATP consuming reactions, contributing considerably to cellular energy homeostasis in batch culture growth.
To quantify the contribution of those additional cellular ATP sinks to the energy metabolism of M. pneumoniae, we defined a single unspecific energy consuming reaction. The minimum constraint of this reaction was fitted manually for each simulated time point in order to allow reproduction of in vivo doubling times. Strikingly, 71–88% of the total ATP available in silico is used for cellular tasks not directly involved in biomass production (Figure 3B; Supplementary Table S6). Consequently, depending on the growth time, only between 12 and 29% of the total available energy is used for cell growth. More precisely, at 36 h of batch culture growth M. pneumoniae uses 9.8% of its total energy for protein production and degradation (protein half‐life of 23 h; Maier et al, 2011) while 8.4% of the total ATP is dedicated to RNA production (mRNA half‐life of 1 min; Maier et al, 2011) and not even 0.1% to DNA synthesis (Figure 3C; Supplementary Table S6; Supplementary information). Lipid production consumes 0.5% of the available ATP while 5.9% are used for metabolic precursor uptake and the subsequent synthesis of secondary metabolites, such as vitamins, FAD, NAD(H), Pi, folate, and other defined functions such as detoxification (Figure 3C; Supplementary Table S6; Supplementary information).
To characterize the yet undefined energy sinks, we further either classified them as growth‐associated maintenance (GAM) and non‐growth‐associated maintenance (NGAM) (Pirt, 1965; Varma and Palsson, 1994b). To estimate the contribution of GAM to the total energy cost, we calculated upper boundaries for the ATP expenses related to post‐translational modifications (0.01%), DNA repair (0.01%), and chaperone‐assisted protein folding (0.25%) based on the available literature and experimental data (Figure 3C, Supplementary Table S6; Supplementary information; Drake et al, 1998; Naylor and Hartl, 2001; Maier et al, 2011; van Noort et al, 2012). Additionally, considering protein turnover costs and other defined expenses mainly related to the uptake and processing of cellular building blocks and reaction cofactors (as defined in the model), total expenses on GAM account for a maximum of 7.6% of the total available cellular energy.
Systematic literature screening for further cellular energy sinks identified proton translocation by the cellular ATPase as most significant quantifiable NGAM task (Kobayashi, 1985, Hutkins and Nannen, 1992). In fact, for lactic acid bacteria, the proton ATPase is the major contributor to maintain pH homeostasis (for a review, see Hutkins and Nannen, 1992). In these bacteria, the ATPase is predominantly involved in creating an optimal proton gradient across the cytoplasmic membrane to allow nutrient import and to maintain the intracellular pH in an acidic environment (Kobayashi, 1985; Hutkins and Nannen, 1992; Moreno et al, 1998). Considering the experimentally determined ATPase abundance (between 99 and 150 ATPase complexes per cell, Supplementary information; Maier et al, 2011) and assuming constant activity at the maximum catalytic rate reported (130 r.p.s.; Wu et al, 2010), we estimated the global ATP hydrolysis rate of the M. pneumoniae ATPase to be maximally ∼38 500 ATP per cell and second, accounting for up to 57% of the total available ATP at simulated 36 h of growth (Figure 3C; Supplementary Table S6; Supplementary information).
To further analyse the contribution of maintenance tasks to cellular energy homeostasis, we challenged growing M. pneumoniae cultures by incubating them for 36 h in medium with pre‐adjusted pH, ranging from pH 5.5 to 8.0. After harvesting the cultures, we measured growth medium glucose levels and protein content, reflecting cellular energy production and bacterial growth, respectively (Figure 4A; Supplementary information). Comparing protein production with glucose consumption for those cultures revealed that at physiological (high) pH, generated ATP is utilized for growth to a much larger extent than at low pH (Figure 4B; Supplementary information). These findings support our in silico energy analysis, particularly the predicted high cellular maintenance costs at later growth stages in batch culture (Figure 3B). Intriguingly, metabolically inactive cells in medium with low pH values resume activity when shifted back to physiological pH (Figure 4C). This indicates that growth arrest in acidic medium is not predominantly caused by cell death but rather by reversible metabolic stalling as a consequence of an unfavorable ion balance for ATP generation. Thus, in M. pneumoniae energy metabolism and biomass production are decoupled depending on growth medium acidity. Below a critical pH, cells enter into a reversible dormant metabolic state.
Taken together, we used our constraint‐based model and experimental results to quantitatively analyze the global energy balance of M. pneumoniae. Considering all quantifiable ATP consuming processes we are able to explain about 75–100% of the total growth stage‐dependent cellular energy expenses (Supplementary Table S6). Biomass production itself does account for about 11–22%, GAM for about 2–7%, and NGAM for 57–80% of the total ATP generated. This high fraction of energy dedicated to NGAM on first sight contradicts findings in E. coli where GAM expenses are far exceeding NGAM costs ∼2.5‐fold to ∼7‐fold (Varma and Palsson, 1994a; Feist et al, 2007). However, artificially adjusting E. coli doubling times to values measured for M. pneumoniae (20 h; Figure 3A) resulted in much higher NGAM expenses in E. coli, and a GAM/NGAM ratio approaching the value calculated for M. pneumoniae (Supplementary Table S7). Furthermore, our results coincide with a recently published whole‐cell model for M. genitalium (Karr et al, 2012) identifying protein and RNA production as the major ATP sinks in biomass production. Movement, which has been shown to be an energy consuming process in M. mobile (Jaffe et al, 2004), could not be estimated so far due to lacking information about the exact function of the gliding machinery and associated ATP costs.
It is important to note that in M. pneumoniae in silico even during the exponential growth phase 78–89% of the available energy is not directly used for the production of cell building blocks but for cellular homeostasis. This finding could also explain the large discrepancy of 44% between energy production and consumption observed in the whole‐cell model for M. genitalium (Karr et al, 2012).
Prediction of flux distributions for varying conditions in silico and in vivo analysis of central carbon metabolism
Introducing the maintenance energy sink reaction into our model and setting its constraint as described above completed the construction process of our metabolic model. The final model, iJW145, contains 306 reactions connecting 216 metabolites and 145 enzymes (Supplementary Figure S1; Supplementary Table S1). Its capability to simulate growth of M. pneumoniae in silico allowed us to reproduce experimentally determined doubling times.
Constraint‐based models describe cellular metabolism under steady‐state assumptions. Integrating experimental data on metabolic parameters and our information about energy homeostasis, we analyzed changes in metabolism during the exponential growth phase in silico. For this purpose, in addition to fitting sigmoidal functions to experimentally measured glucose consumption, protein production, as well as to lactate and acetate secretion during a 96‐ h time course (Supplementary Figure S5), we fitted a polynomial function to the in silico maintenance costs (Figure 3B; Supplementary Table S8). The fitted functions allow the calculation of metabolic constraints for any given time of a 4‐day batch culture growth experiment (Supplementary information).
Simulating growth under rich medium conditions with respectively determined sets of constraints (Supplementary Table S5; Supplementary information) permitted the prediction of metabolic flux distributions for different times of the exponential growth phase, thus providing information about the flux changes in batch culture growth (Figure 5A; Supplementary Table S9). To this end, we used representative flux distributions instead of considering multiple optima since a flux variability analysis (FVA) showed only negligible variability in the majority of the reactions (Supplementary Table S10, see below). We determined those changes between in silico flux distributions for 24, 36, 48, and 60 h of batch culture growth and found that the majority of the reactions (51.6%) show the same trend as biomass synthesis, that is, from 24 to 36 h the flux increases and from 36 to 48 h as well as from 48 to 60 h the flux decreases (Figure 5B). Another 2.6% of the fluxes change contrary to biomass synthesis, that is, first decrease (24–36 h) and then increase (36–48–60 h). Those findings support the notion that microorganisms generally optimize for growth (Neidhardt, 1996) and that the reductive genome evolution of M. pneumoniae eliminated most cellular functions not associated with growth and survival. For 11.4% of all reactions a constant flux increase is observed while 5.9% of the fluxes constantly decrease over time. The constantly changing reactions either belong to glycolysis, pyruvate metabolism, the energy producing arginine fermentation (amino‐acid metabolism) or are associated transport reactions (Figure 5C). This can be explained by the increase in maintenance costs during batch culture growth (Figure 3B) and the subsequent adaptation of the catabolic pathways as well as by the imposed acetate production constraints. Our results suggest the central carbon metabolism to be influenced rather by the external conditions challenging cellular maintenance functions than by the biomass production rate, which further supports our energetics analysis.
Around 1.6% of all reactions have the same flux at all times and two reactions show other changes not being used at all simulated time points. Both of them belong to the nucleotide metabolism. Our FVA revealed that around 14% of all reactions, mostly associated either directly or indirectly (by involving nucleotides as reaction cofactors) to nucleotide metabolism, can be used with differing fluxes to produce optimal solutions for FBA problems (Supplementary Table S10). For example, the two routes to produce deoxy‐CDP and deoxy‐GDP, respectively, are energetically equal and therefore do not alter growth rate and energy homeostasis, which explains the changes observed during the flux analysis. In general, the FVA revealed that the flux distributions predicted have only very low variability with most of the reactions either allowing no or only very little changes (Supplementary Table S10). A total fraction of 26.1% of the reactions is not used under the simulated rich medium conditions. A higher fraction of unused reactions was observed in amino‐acid metabolism, nucleotide metabolism, and pentose phosphate pathway (between 25 and 35%). Hence, these pathways could serve as possible rescue routes for the adaptation to diverse stress conditions. Lipid metabolism also contains 30% non‐active reactions, which in addition to being a direct effect of not including cardiolipin in the biomass (its catalyzing reactions are not used) is in agreement with detected protein quantification data of the respective pathway (Maier et al, 2011). Due to the simulated conditions, reactions involved in the processing of alternative sugars are not active which is in agreement with the low abundance of involved processing enzymes (Maier et al, 2011). In summary, the majority of the reactions of the metabolic network of M. pneumoniae are active even under favorable growth conditions reflecting the reduced genome and the simple linear network structure.
Usually, metabolic fluxes are studied in connection with transcriptomics data to reveal regulatory mechanisms, but in M. pneumoniae mRNA and protein levels correlate only moderately (Pearson's correlation coefficient of 0.41–0.51; Maier et al, 2011). Therefore, we integrated in silico flux changes along the exponential growth phase directly with protein abundance changes (Supplementary information). Aligning qualitative changes for protein abundances and fluxes along the exponential growth phase and in addition considering a 25% error rate for protein abundances (a two‐fold error was reported by Maier et al, 2011), we found that about 86% of all protein abundance changes agree with the in silico flux changes (Supplementary information). These findings are in agreement with a recent study analyzing the dynamic adaptation of Bacillus subtilis to nutritional shifts (Nicolas et al, 2012). The integration of information about post‐translational modifications (van Noort et al, 2012), however, did not lead to further conclusions about the influence of protein levels on metabolic regulation (Supplementary information).
To further validate the flux distributions predicted by the model, M. pneumoniae cells were pulse‐fed with heavy isotope‐labelled glucose (13C6H12O6) and the propagation of the labelled carbon atoms through glycolysis was observed by GC–MS (Supplementary information). To this end, we monitored reporter compounds for the pentose phosphate pathway (ribose 5‐phosphate (R5P)) and for lipid metabolism (G3P and glucose 1‐phosphate (G1P)) (Supplementary Figure S2A). Under rich medium conditions, we found a large excess of heavy isotope reporter ions for all detected intermediates of glycolysis already 15 s (earliest possible time point for reproducible measurements) after supplying 13C‐labelled glucose (Supplementary Figure S2B). This finding confirms the high fluxes predicted for glycolysis reactions in silico.
The model predicts slow influx from G3P and the pentose phosphate pathway into glycolysis and even slower outflux from glycolysis into lipids via G1P when compared with the overall speed of glycolysis in silico. We used R5P, a key intermediate of the pentose phosphate pathway, to check for outflux from glycolysis into the pentose phosphate pathway in vivo. Contrasting the high fluxes observed for glycolysis with 13C‐carbon saturation after 15 s, we did not observe a significant accumulation of labelled reporter ions for R5P on comparable timescales (Supplementary Figure S2A). However, after about 24 h supply of 13C6‐labelled glucose, intracellular R5P pools were fully labelled, indicating a slow, albeit steady flux from glycolysis into the pentose phosphate pathway in vivo.
Two independent branches connect glycolysis to lipid biosynthesis. The first branch involves G1P, a precursor providing hexose sugars for glycolipid synthesis. G1P gets synthesized from glucose 6‐phosphate (G6P) in a reaction catalyzed by the phosphoglucomutase (MPN066, 80 copies per cell; Maier et al, 2011). G1P is rapidly detected during flux experiments. However, the conversion of G6P into G1P is reversible and it cannot be excluded that G1P is only used as an extension of the G6P pool. For the second connection, G3P provides the polar head group to which fatty acids are covalently attached by ester bonds during phospho‐ and glycolipid synthesis (Supplementary Figure S1). When not imported from the growth medium or obtained from its precursor glycerol, G3P can be generated from dihydroxyacetonephosphate (DHAP) in the corrected GPO catalyzed reaction (compare Figure 2A; MPN051: 62 protein copies per cell; Maier et al, 2011). Analysis of the incorporation of labelled reporter ions in G3P revealed a comparatively slow conversion of DHAP into G3P, reaching saturation after 24 h of incubation (Supplementary Figure S2A) suggesting low flux between glycolysis and lipid biosynthesis as compared with flux velocity in glycolysis.
Summing up, our experimental analysis, in consistency with previous suggestions (Yus et al, 2009) and predictions from our own model as well as from the whole‐cell model for M. genitalium (Karr et al, 2012), showed that the flux through glycolysis far exceeds those shuffled to adjacent pathways such as the pentose phosphate pathway or lipid biosynthesis (Figure 6). Comparing in vivo results to flux distributions predicted in silico, we found that the model predicts well the major carbon flux through glycolysis. However, the flux directions between glycolysis and other metabolic pathways (namely between DHAP and G3P or between FBP+GAP and the pentose phosphate pathway) are inverted (Figure 6). This suggests a slight overestimation of the import of ribose and glycerol/G3P, which however should not significantly affect growth rates and overall energy metabolism due to the low fluxes in the branching reactions when compared with the speed of glycolysis (Figure 6).
Prediction of gene essentiality
To test the accuracy of the refined metabolic network on a global scale, we performed an in silico knock‐out study for genes encoding 131 metabolic proteins (Supplementary information; Supplementary Table S11). We systematically silenced, that is, limited to zero, all reactions catalyzed by the individual gene products and recorded the resulting growth ability for rich medium conditions (Supplementary Table S11; Supplementary information). Seventy‐three genes (56% of metabolic enzymes included in the prediction, Supplementary information) were identified as essential, because the respective in silico knock‐out either led to growth arrest (22%) or rendered the FBA infeasible (34%), that is, at least one minimum constraint could not be satisfied (Figure 7A; Supplementary Table SII). Conversely, 58 artificial gene knock‐outs (44%) resulted in objective values larger than zero, either showing no change as compared with the wild type (19%) or with a lower objective value as in the wild type representing a reduced fitness phenotype (25%), thereby predicting these genes to be not essential for growth and survival in M. pneumoniae (Figure 7A; Supplementary Table SII).
We evaluated the predictions by comparing them with a genome‐wide transposon mutagenesis‐based knock‐out screen in the closely related bacterium Mycoplasma genitalium (Glass et al, 2006) (Table III; Supplementary information). This comparison was based on the assignment of functional orthologs in the two mycoplasmas (Supplementary Table S12; Supplementary information). In a first, unbiased analysis of the in silico knock‐out results using gene essentiality in M. genitalium as sole criterion, we achieved 86% accuracy (correct predicted/total predicted) and 98% specificity (true negatives/(true negatives+false positives)) of our in silico essentiality prediction (Table III; Supplementary information). In case of a contradiction between prediction and this transposon study, we screened an M. pneumoniae transposon library (Halbedel and Stülke, 2007) for respective mutants (Supplementary Figure S7). Thus, we could confirm five model predictions that were not supported by the M. genitalium study alone, raising the prediction accuracy to 92%. Considering additional parameters, namely the literature data, and the simulation conditions, allowed the classification of another eight false‐negative hits (Supplementary information). Thus, the model predicts gene essentiality with a final accuracy of 95% and a specificity of 98% (Table III; Supplementary information). We conclude that our refined model of the M. pneumoniae metabolism (iJW145) possesses high predictive power regarding metabolic phenotypes.
To gain a more detailed quantitative understanding of the resulting phenotypic changes, we assayed the relative in silico flux changes of representative flux distributions for the individual reduced fitness knock‐out strains. As expected, most individual reaction fluxes (i.e., fluxes of all reactions under all different conditions) are downregulated (54%) or do not change (34%) in response to the different gene deletions in silico (Figure 7B). However, we also observed several (highly) upregulated reactions (10.4% of all fluxes) and a few changes in flux direction for reversible reactions (0.4%). Only 88 new fluxes in 32 reduced fitness mutant strains (1.1% of all fluxes) have been observed and they are all found in only 15 different reactions, highlighting—in agreement with our FVA—the limited ability of M. pneumoniae to dynamically adapt to perturbation conditions by using alternative metabolic routes for the synthesis of biomass components and ATP. We analyzed highly upregulated reactions as well as those with new or reversed fluxes in at least one in silico knock‐out strain further (Figure 7C). As expected, lactic acid production (M011) gets upregulated in all strains with gene knock‐outs in the acetate branch of pyruvate metabolism. Interestingly, the other 28 examined reactions belong either to the nucleotide metabolism (and reactions involving nucleotides as cofactors) or to the pentose phosphate pathway. This suggests that the non‐essential genes in those pathways allow M. pneumoniae to adapt to environmental changes and hence lead to a more robust metabolic network.
We used our final model (iJW145) to perform an in silico genetic interaction screen by predicting metabolic phenotypes for double knock‐outs of the 58 non‐essential metabolic genes (Supplementary Table S13). Genetic interaction screens are performed to assay network connectivity and to link functionally related genes of different metabolic pathways (Tong et al, 2001, 2004; Szappanos et al, 2011). Especially, the analysis of synthetic lethal and sick interactions, that is, those double mutants that cause cell death or reduced fitness, allows the identification of gene products that impinge on the same biological process (Hartman et al, 2001). The analysis of the synthetic lethal and sick interactions predicted in silico shows that genes involved in pyruvate metabolism have a global effect on the metabolic network behavior (Figure 8, for a list of gene names see Supplementary Table S14). Particularly, mpn674 encoding the lactate dehydrogenase has a strong influence on the phenotype fitness due to the limited acetic acid production. Pyruvate metabolism is central to cellular energy metabolism and mutations in the involved genes limit ATP availability for all cell functions. As expected, the double knock‐outs of the lactate dehydrogenase together with a gene involved in acetic acid synthesis are lethal for M. pneumoniae. Additionally, the genes coding for proteins involved in sugar uptake and processing can limit the synthesis of ATP and are also found enriched among synthetic lethal interactions (P‐value 6.4E−07). Furthermore, pentose phosphate pathway and folate metabolism genes, respectively, are found to be enriched among the synthetic lethal interactions (P‐values 0.00192 and 0.009981), confirming the existence of rescue routes for parts of these pathways suggested by the single knock‐out analysis.
In summary, flux activity analysis at different time points of the exponential growth phase and in in silico knock‐outs coincide suggesting the same metabolic pathways to be responsible for adaptation to perturbations in M. pneumoniae, namely nucleotide metabolism and pentose phosphate pathway. This can be explained by the relatively unchanging environmental conditions M. pneumoniae encounters in its natural habitat leading to the elimination of pathways not required for life in the lung. In agreement with a recently published study relating metabolic and gene co‐expression networks to predict gene essentiality in M. pneumoniae (Güell et al, 2012), our metabolic gene essentiality analysis identified the genes associated with the main catabolic pathways (glycolysis+pentose phosphate pathway) among the non‐essential genes as most important for the metabolic performance further highlighting the simplicity of the metabolic network of M. pneumoniae.
We developed a comprehensive metabolic model, iJW145, for M. pneumoniae. First, we reconstructed the metabolic network based on a manually curated metabolic map (Yus et al, 2009) giving rise to iJW145_recocnstruct. In the second step, we assigned reaction reversibilities and semi‐quantitatively determined the biomass composition of M. pneumoniae resulting in iJW145_growth. Third, the model was improved in an iterative process of in silico growth simulations and their evaluation based on a variety of experimental data and literature information. We were able to correct the structure of the metabolic network allowing the final model, iJW145, to reproduce experimental findings along the exponential growth phase (Yus et al, 2009).
The number of constraint‐based metabolic models of different organisms has been constantly increasing over the past years (for a list of validated models, see Feist et al, 2009 and Supplementary Table 2). For mycoplasmas, an automatic reconstruction for M. genitalium has been published (Suthers et al, 2009). However, since it has been shown that automatic reconstructions are often highly error prone (Reed et al, 2006; Henry et al, 2010), we integrated different experimental data directly during the manual reconstruction process. This enabled us not only to obtain an accurate metabolic reconstruction but also to revise the wiring diagram and the functional annotation of key enzymes.
After establishing the biomass composition of an average M. pneumoniae cell, we quantitatively analyzed its energy metabolism. Comparing our results with that of the recently published whole‐cell model for M. genitalium showed that both models identified protein and RNA production as well as their maintenance as the major energy sinks in biomass production (Karr et al, 2012). Most strikingly, we showed that M. pneumoniae—at least under laboratory conditions—dedicates only a small fraction of its ATP directly to biomass production. Alternative quantified ATP sinks include chaperone assisted protein folding, DNA maintenance, and post‐translational modifications. The combination of in silico calculations taking into account the ATPase protein copy number in M. pneumoniae and catalytic rates from other organisms reported in the literature, as well as the in vivo analysis of growth in medium at different pH suggest that the ATPase uses about 57–80% (growth stage depending) of the total generated energy to maintain a favorable proton gradient across the membrane and the intracellular pH. On the one hand, this surprising finding can be explained by the small size of M. pneumoniae. Membrane leaking and the transport of molecules across the membrane have a higher effect on cytoplasmic homeostasis in small organisms (M. pneumoniae has surface‐to‐volume ratio 2500 times higher than E. coli, Supplementary information; Supplementary Table S7). On the other hand, the continuous, growth‐associated secretion of acids poses an increasing pH maintenance burden on M. pneumoniae grown in batch culture. Therefore, in consistency with the abundance of the core components of the ATPase protein complex (Maier et al, 2011) and our experimental results for M. pneumoniae growth under pH stress (Figure 4A–B), we propose that at later growth stages most of the available energy is diverted toward cellular maintenance with the ATPase as a major energy sink of M. pneumoniae in batch culture growth, as also suggested for other lactic acid bacteria (Hutkins and Nannen, 1992).
In contrast to the high ATPase costs, the costs for protein folding and maintenance by molecular chaperones are surprisingly low taking into account their high cellular abundance (∼10% of the total cellular protein mass; Maier et al, 2011), as is the total amount of energy used for GAM (∼2–7% of the total generated ATP, ∼2.5–10% of NGAM). Despite it has been shown that GAM/NGAM estimations show considerable variance depending on the experimental data used (Varma and Palsson, 1994b; Feist et al, 2007; Orth et al, 2011), our finding is in agreement with recent results from the whole‐cell model in M. genitalium (Karr et al, 2012). Furthermore, we showed that the doubling time has major impact on the ratio between GAM and NGAM since assuming different doubling times while providing the same amount of nutrients leads to a significant alteration of this ratio (Supplementary Table S7). Summing up, we were able to explain 75–100% of the total ATP consumption during the exponential growth phase. Movement and attachment are the only known major energy consuming processes we could not assess in our analysis. M. mobile has been shown to use ATP for gliding (Jaffe et al, 2004) but no details on related energy consumption are available. The whole‐cell model for M. genitalium also showed that other cellular processes which are not considered in our metabolic model, such as chromosome condensation, RNA modification and processing, ribosome assembly, protein translocation, and the replication initiation only contribute to a minor extent to the cellular energy balance (<3%; Karr et al, 2012). The missing expenses can alternatively be attributed to (i) measurement errors in absolute cellular protein quantifications (a two‐fold error has been reported in Maier et al, 2011), (ii) the determined protein and mRNA turnover rates (we used average half‐life times Maier et al, 2011), or (iii) the estimation of cell doubling times (protein quantities in the beginning of a growth experiment, that is, up to 36 h after inoculation, are near the lower detection limit; Supplementary Figure S5A). When comparing our results with those of the whole‐cell model of M. genitalium (Karr et al, 2012), we find that both models agree in the general predictions on central carbon metabolism. Our findings on cellular energy homeostasis (Figures 3B and 4B) can also explain the discrepancy between produced and consumed energy found for M. genitalium in silico (44%; Karr et al, 2012).
The metabolism of most bacteria follows a single objective function: maximizing growth (Neidhardt, 1996; Buescher et al, 2012). A recently applied multidimensional optimality analysis revealed that metabolic flux states in bacteria additionally evolved to minimize the costs for adjusting to different environmental conditions (Schuetz et al, 2012). Our energy calculations and their integration with the model of M. pneumoniae represent a novel approach toward a more detailed understanding of cellular metabolism. They allow a quantitative dissection of cellular performance into energy gain, biomass production and other cellular, growth and non‐growth‐associated energy sinks. A comparison of maintenance energy costs between M. pneumoniae and E. coli revealed fundamental differences in their energy sink reactions, suggesting individual and characteristic energy expense profiles for different bacteria. We identified four parameters governing the composition of these energy expense profiles: the topology of the metabolic network, environmental conditions, growth rate, and cell size. In the case of M. pneumoniae, one can speculate if the large amount of energy diverted toward maintenance tasks and the associated slow growth in the laboratory additionally represent an adaption to its parasitic lifestyle on epithelial cells of the human lung.
Integrating flux distributions predicted in silico and an experimental analysis of the central carbon metabolism by 13C‐labelled glucose tracer experiments, we showed that glycolysis is directly connected to the pentose phosphate pathway and lipid biosynthesis but not to amino‐acid metabolism in M. pneumoniae in vivo. Model predictions and experimental results further agree that most of the carbon taken up is shuffled through glycolysis for the production of ATP during organic acid synthesis. However, in the flux directions of the minor inter‐pathway fluxes they disagree, indicating a slight overestimation of imported ribose and glycerol/G3P. For an exact in silico representation of also the minor fluxes, further experiments addressing nutrient uptake of M. pneumoniae from the environment are necessary but technically not feasible yet. Taking into account the comparatively small amounts of carbons shuffled via the routes interconnecting different metabolic pathways, we estimate the influence on the overall growth rate and energy balance to be very small.
One interesting general finding with respect to network dynamics is that under all simulated conditions (rich/defined medium, alternative sugars, knock‐outs) oxygen consumption is tightly coupled to acetic acid production. This prediction is in agreement with findings in other organisms, such as Lactococcus lactis, in which limited oxygen availability at later growth stages prevents the regulation of the cellular redox imbalance associated with acetic acid production, while the lactate dehydrogenase is released from its supposed oxygen‐dependent inhibition (Gottschalk, 1986; Neves et al, 2005). We propose that oxygen could also have a regulatory role in pyruvate processing in M. pneumoniae explaining the metabolic shift from mainly acetic to mainly lactic acid fermentation observed in vivo during a 4‐day batch culture growth experiment (Yus et al, 2009).
In silico knock‐out studies have been used in other organisms to predict gene essentiality (Reed and Palsson, 2003; Feist et al, 2007). With our prediction for M. pneumoniae we achieve slightly higher accuracy and specificity that has been obtained for E. coli so far (Feist et al, 2007). M. pneumoniae has an exceptionally high percentage of essential metabolic genes (56.6% versus. 19% in E. coli; Baba et al, 2006; Joyce et al, 2006) for which consequently no other gene can buffer the loss of function caused from gene deletion. We found that in M. pneumoniae nucleotide metabolism and the pentose phosphate pathway are the only pathways preserving rescue routes for gene deletion events. The most likely explanation for this lack of rescue routes in most pathways is the adaptation to parasitism in a specific relatively unchanging niche accompanied by the reductive genome evolution.
In vivo, synthetic genetic array analysis has been shown to allow the automating of the isolation and analysis of double mutants (Tong et al, 2001, 2004). We used in silico prediction of double mutant phenotypes to extract information about the specific effects caused by different gene deletions on the metabolic network, and the adaptive capabilities of M. pneumoniae. The results confirm the general findings of the single knock‐out simulations that have been evaluated based on independent experimental data. Especially in organisms for which no engineering tools for in vivo analysis of double mutant phenotypes are available, the applied in silico analysis provides a promising alternative to experimental approaches.
Materials and methods
The used model in sbml format including MIRIAM annotation is provided for download at http://nin.crg.es/serranolab/mycomap/ (user name: mycomap, password: bicha987) and can be accessed in the BioModels database with ID MODEL1301290000.
Constraint‐based modeling and flux balance analysis.
Constraint‐based modeling is an approach to analyze a (metabolic) system under steady‐state conditions. This means that the concentrations of the network components xi do not change over time, that is,
where N is the stoichiometric matrix of size mxn with m being the number of reactants and n the number of reactions and v an n‐dimensional vector containing the fluxes through the n reactions of the network. FBA is a mathematical method to determine a set of metabolic fluxes (the vector v) fulfilling the steady‐state condition and, at the same time, maximizing an objective function such as growth for a given set of available nutrients using linear programming (Varma and Palsson, 1994a; Kauffman et al, 2003). For our model, maximization of biomass production has been chosen as single objective function, since no other objectives have been revealed for M. pneumoniae so far.
We used the reconstruction and modeling platform ToBiN (Toolbox for Biochemical Networks; http://github.com/miguelgodinho/tobin). The initial reaction network was based on the list of reactions found in Yus et al (2009). Some changes had to be introduced to keep elements and charges balanced and to cope with reactions that can be represented in stoichiometric models only with considerable impact on the model complexity, such as RNA and DNA elongation reactions (see below). To simulate the exchange of compounds with the environment reactions producing or consuming given compounds have been defined (source and sink reactions). Appropriate reaction reversibilities and minimum–maximum flux constraints were imposed based on the experimental data and literature information (Supplementary Tables S2 and S5; Supplementary information). We used CellDesigner 4.1 (Kitano et al, 2005) to visualize the model (Supplementary Figure S1) and the abbreviations used can be found in Supplementary Table S15.
Definition of biomass equation.
Based on the general biomass equation (Results), the different components were identified and, if possible, quantified. According to Razin et al (1963), mycoplasma cells are composed of 54–62% protein, 12–20% lipids, 3–8% carbohydrates, 8–17% RNA, and 4–7% DNA. Assuming that one M. pneumoniae cell contains 10 fg of protein (Yus et al, 2009) and assuming proteins to compose 62% of the total cell mass, we were able to calculate the fractions of RNA and DNA. Based on sequence information and the biophysical properties of M. pneumoniae (Yus et al, 2009; Maier et al, 2011), we determined the DNA to account for 5% and the RNA to account for 6.5% of the total cell mass. Further, we assume that 20% lipids and 6.5% carbohydrates and other metabolites make up the missing 26.5% of the cell mass (Supplementary information). The lipid composition of mycoplasmas varies depending on the fatty acids provided with the medium (McElhaney and Tourtellotte, 1969; Pollack et al, 1970; Rottem, 1980). As these variations cannot be represented in a static model, an artificial molecule has been defined taking into account the average fatty acid composition determined (Supplementary Figure S3; Supplementary information, dataSheet). The synthesis of each macromolecule (RNA, DNA, and protein) is represented by a dedicated artificial reaction (Supplementary information). To consider mRNA and protein half‐lives, minimum constraints have been set on the respective degradation reactions. DNA repair as well as rRNA and tRNA degradation have been accounted for qualitatively by small overhang quantities of the respective components in the biomass equation (Supplementary information). Free bases and free amino acids (together accounting for ∼1.5% of the total biomass) were quantified by GC‐MS analysis (unpublished results). To account for the essentiality of vitamins and other cofactors, the end products of the secondary metabolism pathways are included into the biomass equation qualitatively (i.e., in small arbitrary quantities), whenever their precursors have been shown to be essential in the defined medium (Yus et al, 2009). This qualitative consideration is possible, as they are supposed to be low abundant and we proved in a sensitivity analysis that a 10‐fold change of their quantities does not change the general model behavior but only introduces a negligible change in the objective value and thus the in silico doubling yield. The missing proportion of biomass is accounted for by including a respective amount of G6P (representative for carbohydrates) into the biomass equation. To overcome missing experimental information and to restrict the complexity of the model, a number of approximations have been made which are fully listed in Supplementary information.
We varied different macromolecule fractions of the biomass based on the ranges reported by Razin et al (1963), namely protein from 54 to 62%, RNA from 6.5 to 10% (a further increase in the RNA fraction was not possible without altering other biomass fractions), and lipids from 12 to 20%. Furthermore, we tested the impact of a change in the amount of cofactors included in the biomass from 0.1‐ to 10‐fold as compared with the value given in the M. pneumoniae biomass function. In all cases, the total cell weight was maintained constant by accounting for the applied changes by a respective increase or decrease in G6P (representing carbohydrates). We tested the influence on the growth rate (keeping fixed all normal growth constraints, including the minimum constraint for maintenance costs) as well as on the maintenance energy (in this case, manually fitted the minimal constraint for maintenance energy to reproduce in vivo doubling times.
Growth simulations for the exponential growth phase (24–60 h after inoculation) have been accomplished using biomass production as objective function for the FBA. The resulting objective value ov gives information about the doubling time tdoub of an average M. pneumoniae cell. The relation of the tdoub and ov when assuming a constant cell number (1 g of cells is the default in the used modeling platform ToBiN) can be described by tdoub=1/ov. Thus, it is possible to distinguish between growth (ov larger than zero), catabolic activity (growth arrest) (ov equal to zero), and cell death (infeasibility of the FBA). The FBA solution is considered as infeasible if at least one of the minimum requirements specified cannot be satisfied under the given nutrient conditions.
Flux variability analysis.
FVA is a method to assess the potential alternative flux distributions supporting a given FBA objective. After performing FBA optimization, the objective flux (e.g., Biomass production) was fixed to the value 0.01% lower than the FBA result (this slight decrease in the objective flux aimed at avoiding model infeasibilities due to rounding errors) and for each model reaction two FBA runs were performed, one setting the flux maximization and one setting the flux minimization through the reaction as the objective. Reactions for which these two values differ span the space of alternative flux distributions. As these optimizations are not independent of each other, not every flux combination lying within the hypercube spanned by the minimal and maximal flux values supports the initial value of the objective (or even valid flux distribution).
All constraints and fluxes in ToBiN have the unit mmol × g(cells)−1 × h−1. To calculate the total ATP produced by M. pneumoniae in silico, we summed up the ATP production of the different catalytic pathways in M. pneumoniae using the following equation:
flux(prod(ATPtotal)) =flux(prod(LAC))+2 × flux(prod(ACE))+flux(prod(ornithine))
with prod:= production. To know the production per M. pneumoniae cell, we converted the obtained fluxes into molecules × cell−1 × s−1.
Comparison of qualitative changes in fluxes and protein abundances.
First, linear fittings to the in silico reaction fluxes obtained at t=24, 36, 48, and 60 h and to protein abundances measured at t=24, 36, 48, and 72 h during batch culture growth experiments in vivo (Maier et al, 2011, unpublished results). Second, we determined the qualitative overall change of fluxes and protein abundances during the exponential growth phase, considering proteins to change only if the measured abundance difference exceeds 25% of the abundance at t=24 h, thus accounting for the reported experimental error that would otherwise have a high impact especially on the changes of low‐abundant proteins (Maier et al, 2011). Finally, we aligned protein concentration changes with the change of the sum of fluxes of reactions catalyzed by the respective enzyme (Supplementary Figure S6).
Gene essentiality prediction.
For the gene essentiality prediction, the gene–protein relationship has been determined for all reactions for which the catalyzing enzyme is known. In each in silico knock‐out, all reactions catalyzed by the corresponding gene product have been limited to zero flux. A gene is considered as essential when its knock‐out leads to an objective value of zero (no growth but minimum constraints can be matched) or the infeasibility of the FBA (minimum constraints are not fulfilled). The genes encoding proteins catalyzing DNA degradation, protein folding and the ATPase reaction have been excluded from the essentiality prediction since their corresponding functions have not been modeled explicitly. All simulations of this section have been performed using rich medium conditions for 36 h growth time (Supplementary Table S5). Since the amount of energy used for maintenance tasks has been determined by manually fitting the minimum constraint of the respective energy consuming reaction to allow reproduction of the experimentally determined doubling time, those expenses have been neglected for the essentiality prediction. Otherwise, a knock‐out leading to significant slower energy production would result in infeasibility of the FBA. Subsequently, the obtained objective values give no information about the absolute doubling times but only about the relative changes in the growth rate between wild type and knock‐out simulation.
For the prediction of double mutant phenotypes, we applied the same strategy as for the single in silico gene knock‐outs, but simultaneously silenced the reactions catalyzed by two different non‐essential enzymes at a time. Double knock‐outs resulting in reduced fitness, that is, the objective value is smaller than for the two corresponding single knock‐outs alone, or in cell death, that is, the objective value equals zero or the FBA is infeasible, were considered for the analysis of synthetic lethal and sick interactions.
For the statistical analysis of accuracy and specificity of the gene essentiality prediction, we evaluated the prediction results based on a genome‐wide transposon study in M. genitalium (Glass et al, 2006), transposon screens in M. pneumoniae (this study), and the simulation conditions. Computationally and experimentally essential genes are considered as true positives, true negatives are computationally and experimentally not essential, computationally essential and experimentally non‐essential genes are defined as false positives and computationally non‐essential and experimentally essential genes accordingly false negative hits.
All sequence analyses have been performed using the Basic Local Alignment Search Tool (BLAST) for proteins (pBLAST) (Altschul et al, 1997). pBLAST (algorithm pblast) was used as M. pneumoniae uses the TGA codon to encode for tryptophan instead of indicating the end of a gene as in most other organisms. Protein sequences of related organisms (ordered for preference: other mycoplasmas, B. subtilis, L. lactis, E. coli) were obtained from KEGG (Kanehisa and Goto, 2000) or the National Center of Biotechnology Information (NCBI) (Tatusova et al, 1999) and used to perform pBLAST against the M. pneumoniae proteome. Alternatively, M. pneumoniae protein sequences were aligned to the nr‐DB to detect possible homologies. This has been done so (i) to identify the cofactors used by the GPO (MPN051), (ii) to shed light on the NOX isoform (MPN394), (iii) to confirm that a reaction converting UTP into CTP does not exist in M. pneumoniae, and (iv) to search for proteins possibly catalyzing phospholipid production. All pBLAST results are shown in Supplementary Information, pBLAST Results.
GC–MS analysis of fatty acids.
Fatty acids were targeted specifically by using tailored protocols. Depending on the case, growth medium, total cell content, cell pellet, or cytoplasm was analyzed as described in each protocol.
Fatty acid analysis was conducted on all three cellular preparations according to Ghanem et al (1991). Briefly, to a sample preparation (lyophilized in the case of cytoplasm and growth medium, humid in the case of cell pellet or total cell content) 1 ml of NaOH solution (3.75 M in a 1:1 (v/v) MeOH:H2O mixture) was added. The suspension was incubated for 5 min at 100°C, vortexed for 10 s and incubated for another 25 min at 100°C. Samples were cooled to room temperature and added with 1 ml MeOH:HCl solution (0.46:0.54 (v/v) of MeOH and 6 M HCl) followed by 1 ml MeOH:H2SO4 solution (0.46:0.54 (v/v) of MeOH and 50% H2SO4). This mixture was vortexed for 10 s, incubated at 80°C for 10±1 min, and cooled on ice immediately. Subsequently, 1.25 ml of a hexane:ether (1:1 v/v) mixture was added and mixed end‐over‐end for 10 min. The water layer was removed after centrifugation for 2 min at 3000, r.p.m. Then, to the organic layer 3 ml of a NaOH:NaCl mixture (0.3 M NaOH in 4.75 M NaCl) was added and the solution mixed end‐over‐end for 5 min. The aqueous phase was frozen and the organic phase transferred to a new tube. The content evaporated under nitrogen stream at 30°C, and reconstituted in 100 μl hexane/tert‐butyl methyl ether (1:1 v/v). GC‐MS was carried out on a 6890N gas chromatograph coupled with a 5973 MSD (Agilent Technologies, Palo Alto, CA, USA). The fatty acids methyl esters were separated on a Phenomenex Zebron ZB‐5 crosslinked 5% phenyl polymethyl siloxane column (15 m × 0.25 mm i.d., 0.25 μm film thickness). Helium was used as the carrier gas at a constant pressure of 5 p.s.i. A 1‐μl aliquot of the extract was injected into the system operated in split mode (split ratio 70:1). The GC temperature is ramped as follows: initial 150°C, increased to 240°C at 5°C/min, held for 3 min, thereafter increased to 310°C at 30°C/min, and held for 2 min. The injector and transfer line are kept at 280°C, the MS source at 230°C, and the quadrupole at 150°C. The mass range scanned is from 50 to 650 Da. Conventional fatty acids were identified by comparison with the analysis of 1 μl of a reference mixture (Supelco 37 Component FAME Mix) under the same conditions. Other fatty acid methyl esters were identified by monitoring the typical loss of 43 mass units (CH3‐CO−) from the parent ion.
GC–MS analysis of glycolytic intermediates and reporter compounds for adjacent pathways.
For the analysis of glycolysis intermediates (cytosolic extract), samples added with the internal standard (5 μl of methyltestosterone 10 ng/μl) were lyophilized and subsequently dried in a vacuum oven (500 mbar, 50°C) in the presence of diphosphorus pentoxide for at least 4 h. Then, aldehyde and ketone functional groups were converted into methyl oximes with 75 μl of methoxyamine hydrochloride in pyridine (2 g/L) at 40°C for 90 min with intermediate mixing. Subsequently, hydroxyl groups were converted into trimethylsilyl groups with 100 μl MSTFA (N‐Methyl‐N‐trifluoroacetamide) at 40°C for 50 min. Samples were transferred to glass inserts, spun for 5 min at 5000, r.p.m., and the supernatant transferred to a new vial for GC–MS analysis (Chan et al, 2011). In all, 4 μl of glycolytic products was analyzed (split ratio 1:10) using a HP‐Ultra1 crosslinked methyl‐silicone column, 16.5 m × 0.2 mm i.d., film thickness 0.11 μm (J&W Scientific, Folsom, CA, USA) in an Agilent 6890N gas chromatograph coupled to an Agilent 5973 mass selective detector. Helium was used as carrier gas at a constant pressure of 5 p.s.i. The GC temperature is ramped as follows: initial 70°C, held for 1 min, increased to 280°C at 6 °C/min, and held for 1 min at 280°C. For optimal sensitivity, the acquisition, performed in SIM mode, was split into four time segments with three characteristic ions per compound (plus the corresponding 13C isotopes using 15 ms dwell times) in each. Four different time segments ranged from 6 to 16 min (13 ions: m/z 103.1, 117.1, 211.1, 232.2, 236.2, 299.2, 315.2, 369.2, 371.3, 384.3, 386.3, 445.4, 448.4 for phospho‐groups, PEP, GAP, DHAP, and G3P), 16 to 19.5 min (15 ions: m/z 103.1, 104.1, 147.1, 160.1, 162.1, 191.1, 205.1, 207.1, 217.1, 220.1, 307.2, 310.2, 319.2, 323.2, 409.3, fructose, I.S., Glucose, and G1P), 19.5 to 28 min (11 ions: m/z 103.1, 104.1, 299.2, 315.2, 357.2, 359.2, 387.2, 459.3, 462.3, 471.3, 475.3: F6P, G6P, and R5P), and 28 to 32 min (9 ions: m/z 103.1, 104.1, 299.2, 315.2, 357.2, 359.2, 387.2, 459.3, 462.3: FBP) (bolded are the ions corresponding to the 13C‐labelled compounds). Identification of the sample constituents was based on the theoretical fragmentation pattern expected for each compound. Quantification was performed with respect to the corresponding calibration curve standardized against the internal standard (methyltestosterone‐MO‐TMS). Samples and the calibration curves were analyzed at least in triplicate.
The 64 pools of an ordered collection of M. pneumoniae transposon mutants generated by ‘haystack mutagenesis’ (Halbedel et al, 2006) were assorted into 10 groups. Then, genomic DNA extractions were performed using Illustrabacteria genomic KIT (GE). The disruptive insertions in genes mpn133, mpn321, mpn392, mpn533, and mpn595 were detected by PCR (Supplementary Figure S7). The fragments corresponding to junctions between the genes and the mini‐transposon were amplified using the primer 3JpMT85 and the primers 5MPN133, 5MPN321, 5MPN392, 5MPN533, and 5MPN595, respectively (Supplementary Table S16). The position of the transposon insertion in the different genes was determined by DNA sequencing.
pH Re‐buffering experiment.
To test for the influence of the medium pH on growth performance, M. pneumoniae cells were grown in batch culture in 75 cm2 culture flasks. Cells were grown in pre‐culture for 96 h in glucose containing medium, harvested by scraping and diluted into fresh growth medium. Medium pH was adjusted back to pH 7.7 after 4 days of growth by titration with sterile 1 M NaOH. Samples from growth medium supplemented with 1% glucose (55.5 mM) were taken at indicated time points (Figure 4C). Glucose and lactic acid concentrations were determined with enzymatic assays.
M. pneumoniae growth in medium at different pH.
M. pneumoniae cells were grown for 48 h in 75 cm2 culture flasks with Hayflick medium supplemented with 55 mM glucose. The medium of each growth flask was exchanged with Hayflick medium with a pre‐set pH (between pH 5.5 and 8) and grown for additional 36 h. Triplicates have been used for each pH. Metabolite and protein measurements were done as in Yus et al (2009).
We thank Jörg Stülke from the University Göttingen for allowing us to screen their transposon library. We thank past and present colleagues for helpful discussions; in particular, Javier Delgado for help on scripting; Konstantinos Michalodimitrakis for experimental supervision; and Hinnerk Eilers for information on protein activity. We thank Toni Hermoso for the conversion of Supplementary Figure S1 into a clickable metabolic map. This work was supported by the European Research council (ERC) advanced grant, the Fundacion Marcelino Botin and the Spanish Ministry of Research and Innovation to the ICREA researcher LS. JW is recipient of a laCaixa‐CRG fellowship.
Author contributions: JW developed and organized the project, designed the model, analyzed the results, prepared figures, and wrote the manuscript. JP designed the model and analyzed the results. ML‐S performed the mutant screens and prepared Supplementary Figure S7. JM conducted the mass spectrometry experiments. EY performed the pH stress experiments to validate the ATPase hypothesis. MG provided data extraction and analysis tools. RGG analyzed the mass spectrometry results. VMdS initiated the project and contributed to its development. LS and EK contributed to the project design and development and discussed results. TM participated in the project development, performed experiments, contributed to the figures, and wrote the manuscript.
Conflict of Interest
The authors declare that they have no conflict of interest.
Supplementary Figures S1–8, Supplementary Table Legends S1–16
17 Tables, the last table contains requested additional raw data
SBML Model with MIRIAM annotation. ID in the BioModels database: ID MODEL1301290000
Source data for Supplementary Figure S2
Source data for Supplementary Figure S3A
Source data for Supplementary Figure S3B
Source data for Supplementary Figure S3C
This is an open‐access article distributed under the terms of the Creative Commons Attribution License, which permits distribution, and reproduction in any medium, provided the original author and source are credited. This license does not permit commercial exploitation without specific permission.
- Copyright © 2013 EMBO and Macmillan Publishers Limited