Genome-scale metabolic reconstruction and metabolic versatility of an obligate methanotroph Methylococcus capsulatus str. Bath

Model reconstruction

The metabolic reconstruction of Methylococcus capsulatus Bath (Mcap) was carried out using its annotated genome sequence obtained from NCBI. The complete reaction network was built using primary and review literature, biochemical databases such as KEGG (http://www.genome.jp/kegg/) and MetaCyc (http://metacyc.org/), and automated GSMM pipelines such as ModelSEED (Henry et al., 2010) and BioModels (Chelliah, Laibe & Le Novere, 2013). The model was manually constructed on a pathway-by-pathway basis, beginning from the core set of reactions (central metabolic pathways) and expanding it manually to include the annotated reactions and pathways associated with Mcap. Annotations for reaction IDs and compound IDs used in the constructed model were obtained from the KEGG database. The neutral formula for the compounds was obtained from KEGG database or EMBL-ChEBI database. The charge on the metabolites and charged formula were calculated at pH 7.2. All of the reactions included in the network were both elementally and charged balanced and were categorized into either reversible or irreversible. If available, the directionality for each of the reactions was determined from primary literature data, or from the KEGG database. The other reactions were assumed to be reversible. However, some reactions were made irreversible to resolve futile cycles in the model. The gene-protein-reaction assignments were made from the genome annotation and the KEGG database. Spontaneous reactions were included in the reconstruction if evidence suggested their presence, such as the presence of at least the substrate or product in the reconstruction, or on the basis of literature evidence. The final model included a list of metabolic reactions, required enzymes and their respective genes, gene protein reaction association, the major subsystem of the reaction, metabolite name, and compartment (p- periplasm or c- cytoplasm, e- extracellular). The model was converted into an SBML format recognizable by COBRA Toolbox for FBA.

Biomass composition

Biomass components, their respective percentages in dry cell weight (DCW), references, and corresponding coefficients in biomass equation are summarized in Table S1. The relative percentages of deoxynucleotides were calculated using the GC content (63.58%) reported for Methylococcus capsulatus Bath (Ward et al., 2004). The amino acid composition in proteins, relative percentages of ribonucleotides, and ash were assumed to be similar to that in Methylomicrobium buryatense 5GB1 (De la Torre et al., 2015; Gilman et al., 2015). As literature references on the composition of carbohydrates, cell wall components (LPS and peptidoglycan), vitamins, and cofactors in DCW for Methylococcus capsulatus Bath are unavailable, these components were assumed to be in the same proportion as reported for Methylomicrobium buryatense 5GB1 (De la Torre et al., 2015; Gilman et al., 2015; Puri et al., 2015), Methylomicrobium alcaliphilum 20Z (But et al., 2013; Ivanova et al., 2006; Kalyuzhnaya et al., 2013; Khmelenina, Kalyuzhnaya & Trotsenko, 1997), Methylomonas methanica (Khmelenina et al., 1994; Trotsenko & Shishkina, 1990), Methylobacterium extorquens AM1 (Peyraud et al., 2011) or Escherichia coli (Wientjes, Woldringh & Nanninga, 1991), in proportions as shown in the Table S1. The lipid content was assumed to be the same to that in Methylomicrobium buryatense 5GB1 (De la Torre et al., 2015) (Table S1). The fatty acid composition of the phospholipids was compiled from primary literature (Makula, 1978; Trotsenko & Shishkina, 1990). For phosphatidylcholine and cardiolipin, the building blocks (fatty acids) were assumed to be the same as phosphatidylglycerol and ethanolamine and each species was supposed to be in equal proportion. The fatty acid composition for other lipid species was considered to be similar to that of phospholipids. Calculations related to the biomass equation are provided in Table S1.

Model simulation and analysis

Energy parameters and main energetic assumptions

Bacteria require energy in the form of ATP for growth (e.g., macromolecular syntheses such as DNA, RNA, and Protein) and cellular maintenance. Depending on the growth rate of cells, growth associated ATP maintenance could vary, hence in this study, two different conditions were tested (i) ATP maintenance value of 23 mmol/gDCW (same as the low value used for E. coli (De la Torre et al., 2015; Varma, Boesch & Palsson, 1993)), and was considered as low ATP maintenance, (ii) ATP maintenance value of 59.81 mmol/gDCW, and was considered as wild type network (same as high value used for E. coli (Feist et al., 2007)). To account for non-growth-associated maintenance requirement, an ATP hydrolysis reaction was added to the model and the flux of this reaction was fixed to 8.39 mmol/gDCW/h (Feist et al., 2007). The ATP is produced by the ATP Synthase complex which involves the translocation of protons between periplasm of cell membrane and cytosol.

Evaluation on nitrogen and carbon sources

To evaluate the metabolic versatility of Mcap, different nitrogen and carbon sources were tested and the growth rate was observed. For testing different carbon sources, the methane uptake was set to zero and the flux for the respective carbon source was set to 28.10 mmol/gDCW/h (the methane uptake at the optimal growth of 0.37 h⁻¹) unless mentioned otherwise in the result section. Similarly, in the case of different nitrogen sources, the flux value for the uptake of the default nitrogen source, that is, nitrate, was set to zero and the lower bound of the reaction flux for the respective nitrogen source was set equal to the uptake rate of nitrate in optimal growth conditions (0.37 h⁻¹).

Identification of Co-sets

Correlated sets (Co-sets) are the sets of reactions whose fluxes remains functionally correlated, that is, remain correlated in all possible states of the network under a condition. The Co-sets for Mcap model were calculated by the COBRA Toolbox function “identifyCorrelSets.” We used uniform random sampling using the Artificial Centering Hit-and-Run sampler while keeping the maximum correlation (R²) to 1.0 for the sampling to be used for “identifyCorrelSets.”

Characteristics of the metabolic reconstruction

The metabolic reconstruction for Methylococcus capsulatus Bath (Mcap) “iMC535” was carried out based on its annotated genome sequence, available literature, information from biochemical databases such as KEGG and BRENDA, and automated GSMM pipelines such as ModelSEED and BioModels (see Methods section for more details). The final constructed model contained 899 reactions corresponding to 535 genes with 865 metabolites distributed over two different cellular compartments: cytoplasm and periplasm (Table 1). The extracellular space was used as a pseudo-compartment. Depending on the cellular localization, each metabolite was placed in one or more of the three compartments, and the flux across the outer and inner membranes was enabled by transport reactions. Out of the 899 reactions, 766 were associated with metabolic conversions, while 133 were required for the transport or exchange of nutrients. A total of 46 non-annotated metabolic reactions were included in the model because either they were known to occur in Mcap based on the literature reports, or were required to fill the gaps in the network (Table S2).

All reactions constituting the model were distributed into 14 systems on the basis of major metabolic properties of Mcap (Fig. 1). The majority of the reactions were associated with amino acid, exchange, carbohydrate, lipid, and nucleotide metabolism. The information for each of the pathways identified downstream of core metabolism was verified manually. Central metabolic pathways, including carbohydrate metabolism, energy metabolism, and biosynthetic pathways for fatty acids, amino acids, nucleic acids, lipids, cofactors, and vitamins along with other correlated pathways, were manually reconstructed. Similarly, the pathways related to nitrogen uptake, utilization of sulfate, ammonia, nitrate, urea, and phosphate were also manually reconstructed. Each reaction was tested for mass balance and irreversibility constraints. The complete functional model “iMC535” is provided in Table S3 in a Spreadsheet format, and in File S1 in SBML format.

Metabolic versatility and energy source determination

Flux balance analysis is widely used to analyze the reconstructed metabolic networks. The model performance was evaluated by comparing in silico performance with previously reported experimental observations for growth and response to different carbon, nitrogen, and energy sources.

Methane utilization

A summary of methane metabolism and central metabolic pathways is shown in Fig. 2. In Mcap, methane can be oxidized by both sMMO and pMMO depending on the copper concentration in the environment. It has been reported that the growth rates for Mcap lie between 0.25 h⁻¹ and 0.37 h⁻¹ (Joergensen & Degn, 1987). However, the experimental values of methane and oxygen consumption rates at the maximum-reported growth rate are not available. When we used the methane consumption rate of 18.5 mmol/gDCW/h, which has been reported for a related bacterium (Methylomicrobium buryatense 5GB1) (De la Torre et al., 2015), the growth rate was simulated to be 0.23 h⁻¹ (similar to Mb5GB1), and the oxygen uptake at this growth was 23.26 mmol/gDCW/h. The methane uptake rate of 28.10 mmol/gDCW/h was predicted at the optimal growth rate of 0.37 h⁻¹, and was used for all the downstream analysis. Therefore, the methane intake rates of Mcap are expected to be higher than that for Methylomicrobium buryatense 5GB1. As the FBA maximized the biomass yield on methane, in reality the biomass yield is expected to be lower. That is, the methane intake rate is very likely higher than that suggested by FBA. The model-calculated ratio of uptake values of oxygen and methane corroborates with the available experimental values which suggests that almost an equal or a slightly higher amount of oxygen is required as compared to methane for the optimum growth of Mcap (Islam et al., 2015).

The growth rates of 0.28 h⁻¹ and 0.37 h⁻¹ were observed in the case of pMMO knockout (low copper) and sMMO knockout (high copper) simulated conditions, respectively (Table 2). These values corroborate well with the experimentally observed growth range of 0.25–0.37 h⁻¹ (Joergensen & Degn, 1987) in these conditions. In terms of carbon conversion efficiency (CCE), a 34.14% increase is observed from low-copper to high-copper simulated conditions, which is very close to the previously reported experimental value of 38% increase (Leak & Dalton, 1986). Since, NADH is a limiting factor for the growth efficiency in Mcap (Leak & Dalton, 1986), the lower growth observed in low-copper conditions when sMMO are active can be attributed to the utilization of NADH by sMMO for the oxidation of methane. On the contrary, in high-copper conditions, the higher growth efficiency is attributed to a reduced NADH requirement for methane oxidation, since pMMO uses quinones for the oxidation of methane (Larsen & Karlsen, 2016; Shiemke et al., 1995; Zahn & DiSpirito, 1996). Additionally, iMC535 also predicts the release of 0.37–0.39 moles of CO₂ per mole of methane consumed at the optimum growth range of 0.25–0.37 h⁻¹, which is in agreement with the previously-reported experimental value of more than 0.3 moles of CO₂ per mole of methane consumed (Whittenbury, Phillips & Wilkinson, 1970). These results further validate the reconstructed model.

Table 2:

Computational predictions for different knockout simulations and growth conditions.

	Biomass flux (h−1)	O2 consumption (mmol/gDCW/h)	O2:C ratio	CO2 production	Nitrate consumption (mmol/gDCW/h)	CCE (%)
Wild type network	0.37	34.18	1.22	12.54	3.41	55.36
Low ATP maintenance	0.44	29.85	1.06	9.49	4.07	66.23
Knockout simulations
sMMO knockout	0.37	34.18	1.22	12.54	3.41	55.36
pMMO knockout	0.28	39.77	1.42	16.49	2.54	41.30
EMP knockout	0.36	34.49	1.23	12.77	3.36	54.57
Carbon source
Acetate	0.59	21.12	0.38	31.42	5.42	44.09
Methanol	0.43	16.40	0.58	9.91	3.98	64.74
Formaldehyde	0.28	11.67	0.42	16.49	2.54	41.30

DOI: 10.7717/peerj.6685/table-2

Note:

The intake flux for each of the carbon source was set as 28.10 mmol/gDCW/h.

After the oxidation of methane into methanol, the next step of formaldehyde formation is catalyzed by ferricytochrome c-dependent methanol dehydrogenase, which can be utilized by one of the following pathways: Tetrahydromethanopterin pathway (THMPT), Tetrahydrofolate (THF)-linked carbon transfer pathways, RuMP cycle, or Pentose phosphate pathway (PPP). In Mcap, formaldehyde is oxidized to CO₂ via either THMPT or THF pathway and is assimilated in biomass through the RuMP. The flux distribution shows that 56.2% of the formaldehyde goes to the RuMP pathway and 43.2% goes to the THMPT cycle. Thus, the RuMP pathway is 1.2 times more active than THMPT pathway (Fig. 3). Mcap has three variants for pathways downstream of the RuMP cycle: Entner–Doudoroff (ED) pathway, EMP pathway, and Bifidobacterial shunt (BS) (Kao et al., 2004). All these downstream pathways are interconnected through the transaldolase reactions of the PPP and result in ribulose-5-phosphate regeneration or pyruvate synthesis. We observed that the flux through PPP and ED pathway is almost three times more than that of glycolysis (Fig. 3). ED pathway is important as it supports the only known NADP-linked reaction in C-1 oxidation in Mcap via 6-phosphogluconate dehydrogenase, and thus could serve as the NADPH source for biosynthesis (Kao et al., 2004). The model did not show any growth in ED pathway knockout simulations, which is mainly because 6-phospho-D-gluconate (an intermediate of this pathway) is required for the growth of this bacterium, and can only be synthesized by the enzyme phosphogluconolactonase present in the ED pathway. Additionally, the gene for gluconokinase, which provides an alternate source for 6-phospho-D-gluconate, is also absent in the Mcap genome and hence, in the absence of any alternate source, the ED becomes essential for Mcap. On shutting down the EMP pathway, the predicted growth rate was only slightly reduced to 98.56% of the optimal growth (Table 2).

The ratio of the activity of the enzymes involved in the ED (6-Phosphogluconate dehydrase/phospho-2-keto-3-deoxygluconate aldolase) and EMP pathway (Phosphofructokinase and Fructose diphosphate aldolase) has been reported as 3.3:1 (Kao et al., 2004). With methane intake of 28.10 mmol/gDCW/h (at a growth rate of 0.37 h⁻¹), our model showed flux through ED and EMP pathways at a similar ratio of 3.4:1 without any external constraints, and hence corroborates very well with the available literature evidence. Additionally, in the robustness analysis, it was observed that the ED:EMP ratio can go up to a maximum of 10.9:1 for the optimal growth of 0.37 h⁻¹. All these observations are consistent with the biochemical studies of gammaproteobacteria, which point toward the role of ED pathway as a major route for C-1 assimilation even though the energy requirement of this pathway is greater (Anthony, 1982; Boden et al., 2011; Strom, Ferenci & Quayle, 1974). The flux through BS under normal simulation conditions was zero, suggesting that this pathway was non-essential in this organism.

An intriguing property which makes this bacterium unique is the presence of genes of both RuMP cycle and serine cycle. While the former plays a major role in biomass production, the role of the serine cycle is still unclear. It has also been hypothesized in earlier studies that serine cycle cannot operate as the major assimilatory pathway in Mcap unless glyoxylate, which is depleted during the cycle, is regenerated. However, no genes were identified in the Mcap genome which can convert acetyl-CoA to glyoxylate (Chistoserdova, Vorholt & Lidstrom, 2005). In our unconstrained network, we observed that all the reactions involved in the serine cycle had non-zero flux values and hence, were active. However, it is also important to note that all these reactions were also involved in other major pathways, and hence this observation points toward a secondary role of reactions of this cycle. In the case of serine cycle knockout simulations in iMC535, it was observed that only one reaction of this cycle was essential for the growth of Mcap, and had a role in glycerate synthesis and purine metabolism. These observations of flux distribution are consistent with the fact that Methylococcus capsulatus bath is type I methanotroph and a member of the gammaproteobacteria family which utilizes RuMP pathway instead of the serine pathway for formaldehyde assimilation.

Analysis of essential genes and flux variability in Mcap

iMC535 was used as a framework to identify the candidate essential genes in Mcap. To simulate deletion of the reactions associated with each gene, the flux through the reaction was set to zero and the biomass reaction was maximized by FBA. Out of the 535 metabolic genes, the deletion of 161 (30.1%) genes was predicted to be lethal (biomass flux less than 10⁻⁶) (Table S4; Fig. 4). As expected, gene deletions in central metabolic processes had a highly negative impact. However, most of the single gene deletions for methane oxidation were not lethal, since multiple copies of these genes were present in the genome. Interestingly, nitrite reductase gene was found essential for Mcap when it was grown on nitrogen sources such as nitrite and nitrate, whereas this gene was not essential when the model was grown on amino acids as the nitrogen source. The majority of the essential genes belonged to pathways associated with amino acid metabolism, carbohydrate metabolism, nucleotide metabolism, and metabolism of vitamins and cofactors. In addition to single-gene deletion study, double gene deletion study was also carried out. A total of 103 unique gene pairs showed multiple lethal/sick interactions, that is, growth rate less than 0.01 h⁻¹ with respect to either of the single gene deletion condition, out of which 76 interactions were lethal (Table S5). The genes showing lethal interactions were mostly involved in amino acid metabolism, exchange, carbohydrate metabolism, lipid metabolism, nucleotide metabolism, and metabolism of vitamins (Table S5).

In a metabolic network, optimal flux through objective function can be achieved through a number of alternative flux distributions if the alternate routes for the same function are available in the bacterium. Such variability in the solution space can be explored using FVA, which computes the range of allowed flux values for each reaction while keeping the optimum value of growth rate. This analysis revealed that 99 out of 383 active enzymatic reactions could show more than 0.1 mmol/gDCW/h flux change without affecting the growth rate. Other 15 enzymatic reactions showed variability between 10⁻⁴ and 10⁻¹ mmol/gDCW/h. A total of 158 enzymatic reactions with zero flux could potentially carry flux in an alternate optimal solution space. On the other hand, 269 reactions had fixed flux (less than 10⁻⁴ mmol/gDCW/h variation) at optimal growth value and represented 70.23% of the active network. High variability was observed for reactions related to glucose, fructose, respective phosphates and bisphosphates because alpha-D, beta-D, and D compounds of above-mentioned species were considered as different entities. Glycolysis, EMP, EDD, serine cycle, fatty acid synthesis, etc., showed high variability due to alternate reactions or presence of the above entities.

Network topology and reaction subsets

In the metabolic network, only a few metabolites were found to be highly connected, whereas the majority of the metabolites showed limited connections. The highly connected metabolites were referred to as hubs, which are crucial to the network since they connect the other metabolites (with few connections) to the central metabolism. Top 35 highly connected metabolites out of 827 unique metabolites in the iMC535 model are shown in (Fig. S1). Ranking of most of the highly connected metabolites in Mcap is similar to that observed in other models (De la Torre et al., 2015). As apparent, Mcap model also had the common metabolic hubs such as water, energy carrying molecules (ATP, NADH, NADPH, etc.), CoA, pyruvate, etc., which were expected due to their common functions in the metabolic network. Interestingly, amino acids, such as glutamate, glutamine, aspartate, alanine, and serine also showed high connectivity (Patel & Hoare, 1971).

In a metabolic model, reactions with linearly correlated fluxes form the correlated reaction set (Co-set) (Marashi & Hosseini, 2015; Price, Schellenberger & Palsson, 2004). Earlier studies on reaction correlation showed that uniformly sampled fluxes can be used to explore the regulatory properties of reactions in metabolic models (Barrett, Price & Palsson, 2006; Barrett, Herrgard & Palsson, 2009). There are 88 Co-sets present in the Mcap model, out of which the largest Co-set had 146 correlated reactions. Further analysis of this Co-set revealed that these reactions were involved in the synthesis of biomass precursors, and hence showed a correlation in all the conditions considered for this analysis. There were 17 reactions present in Co-set 2 participating in fatty acid biosynthesis, 12 reactions in Co-set 3 participating in lipopolysaccharide biosynthesis, six reactions in Co-set 4 participating in methane metabolism and five reactions in Co-set 5 participating in sulfur metabolism. From Co-set 6 onward, there were 7, 14, and 62 Co-sets having four, three and two reactions, respectively. The summary and detailed information of Co-sets and its participating reactions is given in Table S6.