Methane utilization in Methylomicrobium alcaliphilum 20ZR: a systems approach

Biological methane utilization, one of the main sinks of the greenhouse gas in nature, represents an attractive platform for production of fuels and value-added chemicals. Despite the progress made in our understanding of the individual parts of methane utilization, our knowledge of how the whole-cell metabolic network is organized and coordinated is limited. Attractive growth and methane-conversion rates, a complete and expert-annotated genome sequence, as well as large enzymatic, 13C-labeling, and transcriptomic datasets make Methylomicrobium alcaliphilum 20ZR an exceptional model system for investigating methane utilization networks. Here we present a comprehensive metabolic framework of methane and methanol utilization in M. alcaliphilum 20ZR. A set of novel metabolic reactions governing carbon distribution across central pathways in methanotrophic bacteria was predicted by in-silico simulations and confirmed by global non-targeted metabolomics and enzymatic evidences. Our data highlight the importance of substitution of ATP-linked steps with PPi-dependent reactions and support the presence of a carbon shunt from acetyl-CoA to the pentose-phosphate pathway and highly branched TCA cycle. The diverged TCA reactions promote balance between anabolic reactions and redox demands. The computational framework of C1-metabolism in methanotrophic bacteria can represent an efficient tool for metabolic engineering or ecosystem modeling.

Flux balance model and growth parameters. The model consists of three compartments (extracellular space, periplasm and cytoplasm) and 396 internal enzymatic reactions (represented by 407 genes), 16 transport reactions between cell compartments, 422 metabolites and 20 exchange rates for nutrients and excreted compounds (see Supplemental Material, Table S1).
Since the energetics of methane oxidation in 20Z R can be linked to aerobic or anaerobic respiration and fermentation, electron transfers were further refined by incorporating aerobic respiration pathways (from NADH/ H + , FADH 2 , and cytochrome C L ), anaerobic respiratory pathways, and fermentation pathways 16 . As the source of electrons is still unknown, we incorporated a pMMO-specific electron carrier (e-donor), which could be reduced via three pathways: 1) a ubiquinone (UQH 2 )/e-carrier to represent UQ-linked methane oxidation; 2) a cytochrome c L (cytC L )/e-carrier to represent direct coupling between pMMO and methanol dehydrogenase; and 3) cytc L and cytbc1/2e-carriers to represent previously proposed reverse (also known as uphill) electron transfer 18,24 . The electron transfer reactions were connected to the production of proton/ion motive force (PMF), with the following ratios: electron transfer from NADH/H + to ubiquinone linked to the transfer of 4 H + (or 0.4 PMF); electron transfer from ubiquinol to cytochrome bc1 (cytBc1) and to cytochrome c (cytC) linked to the transfer of 4 H + (or 0.4 PMF); electron transfer from cytochrome c (cytC) or from cytochrome c L (cytc L ) to cytochrome oxidase linked to the transfer of 2 H + (or 0.2 PMF). The electron transfer from ubiquinol to ubiquinol oxidases was also linked to the translocation of 4 H + (or 0.4 PMF). ATP production via oxidative phosphorylation was set as the following: 1 ATP produced per 3 H + translocated. The initial value of non-growth and growth-associated ATP maintenance were adapted from Methylomicrobium buryatense 5 G(B1) 24 .  Additional constraints were added as follows: 1) only water and CO 2 are included as the expected excretion products; excreted formate, acetate, lactate, and exopolysaccharides were incorporated as part of the biomass equation and included in the biomass flux; 2) the methane consumption rate was set at 11.7 ± 0.14 mmol g −1 DCW h −1 , (DCW, dry cell weight) based on measurements taken during continuous culture growth (at observed μ = 0.12 ± 0.01 h −1 ); and 3) the oxygen consumption rate was set at >11.7 and <20.6 mmol gDCW −1 h −1 . The lowest parameter of oxygen consumption is based on the stoichiometry of methane oxidation, as 1 mol of oxygen is needed for oxidation of 1 mol methane. The upper limit represents the highest oxygen consumption rate observed during methane-limited growth (unpublished data). The oxygen consumption rate during exponential growth in continuous culture was 14.7 ± 0.9 mmol g −1 DCW h −1 .
The initial metabolic model was further refined using enzymatic data (compiled in Table 3), and published 13 C-carbon-labeling analyses, and transcriptomic data 16 . A number of central metabolic pathway reactions, including a set of reactions for the anaplerotic carbon fixation and PPi-linked reactions 29 were validated and the following modifications were made: 1) both oxaloacetate and malate were set as key entry points for CO 2  (3) all variants of the RuMP and the serine cycle were incorporated into the model; 24,35 (4) all possible variants of the TCA cycle were included 15,36 .
An overview of the central metabolic pathways is shown in Fig. 1 (a detailed figure is provided as supplemental Figure S1). We ran a set of in silico flux balance analysis (FBA) experiments using the COBRA toolbox to validate the stoichiometric completeness of the proposed metabolic network at all steps of the model development.
Simulations for unconstrained network supported the direct coupling (electrons for methane oxidation comes from cytC L ) as the most optimal solution for methane oxidation, however the predicted O 2 /CH 4 consumption was low ( Table 4). As could be predicted, incorporation of PPi-reactions slightly improved biomass yield. An additional constraint on intracellular source of PPi, with biosynthesis reactions as the only source, lead to the reduction of biomass and an increase in O 2 consumption. However, even after the modification, the predictions differ from experimental measurements (Table 4). Specifically, the predicted growth rate was higher (0.145 vs 0.12), while predicted O 2 :CH 4 consumption rate was low (1.21 vs 1.26). Two possible explanations were explored further: 1) ATP maintenance in methanotrophic bacteria is higher than assumed 38 ; and 2) ubiquinol (UQH 2 ) serve as supplemental sources of electrons for methane oxidation. The ATP requirements in methane grown cells of M. alcaliphilum 20Z R were measured, and the contribution of UQH2 was investigated by further constraining metabolic networks to observed oxygen consumption rate. ATP maintenance requirement in M. alcaliphilum 20Z R . Methane-limitation experiments were carried out, in order to estimate non-growth-associated ATP maintenance 39 This study  various concentrations of methane from 0 to 2 mmol and the growth rate and methane consumption rate were measured. We found that 20Z R cells require 3.6 ± 0.3 mmol CH 4 gDCW −1 h −1 to sustain metabolic activity but not grow. Methane concentrations below this parameter led to decreases in cell density, while higher methane concentrations supported cell growth. Using the assumed ATP yield Y ATP, max = 6 mol ATP per one mol CH 4 consumed, the non-growth-associated ATP requirement was calculated as 21.6 mmol ATP g DCW −1 hr −1 . Growth-associated ATP maintenance for 20Z R was found by fitting the reconstructed model to the experimentally determined biomass yield. We plotted the CH 4 -uptake rate against the growth rate from continuous culture experiments (Table 1). Model simulations were carried out for various values of growth-associated ATP maintenance and the range of 90-110 mmol ATP g DCW −1 was found to match experimental data ( Fig. 2 67 3-Hexulose phosphate synthase* Ru5P 172 ± 9 0.13 68 Phosphoglucose isomerase 32 ± 2 -"-  While the updated maintenance coefficients improved the correlation between simulations and observed experimental data, the oxygen consumption rates were still below expected 1.26.

Model predictions depending on simultaneous constraints of O 2 and CH 4 consumption. A set
of additional in-silico simulations with constrained methane and oxygen consumption rates was carried out in order to investigate correlations between biomass production and O 2 /CH 4 ratios (Fig. 3). The optimal solution for biomass production directly correlated with oxygen consumption, reaching the optimal output at 1.    4 consumption ratio of 1.5. The switch to the redox arm leads to a decrease in biomass production. As the majority of the experimental data lay in the range of 1.25-1.4, optimal biomass production could not be modeled by any single mode, whether it was UQ-mediated or direct coupling. In summary, the most realistic behavior of the system could be modeled by dual simultaneous constraints on oxygen and methane consumption.
To further validate this prediction, we estimated the parameters for growth at low oxygen tension. It has been shown that upon oxygen limitation the strain switches to a fermentation mode 16 . We ran an initial unconstrained model with only a constraint on oxygen uptake (O 2 consumption rate was set to 12 mmol gDCW −1 h −1 , or O 2 / CH 4 = 1.02). This modification resulted in significantly reduced bacterial growth (0.06 h −1 ) and an extremely elevated excretion of acetate (3.83 mmol gDCW −1 h −1 ). According to the in-silico simulations, the phosphoketolase pathway is the most compelling mode of methane fermentation, with the majority of cellular ATP produced from an acetate kinase reaction. Both the growth rate and the increase in acetate excretion have been previously observed as main outcomes of O 2 -limiting growth 16,40 . Furthermore, the kinetic properties of purified acetate kinase from M. alcaliphilum 20Z support the in-silico prediction (Table 3 31 ).
Model validation: metabolomics data. The validity of the metabolic model was further assessed by metabolomics. Metabolomic profiles were obtained for cells grown in continuous culture on methane, as well as for exponentially grown cells from batch cultures on methane and methanol. Two hundred and fifty-three compounds were identified in bacterial cell pellets. The complete list of the non-targeted metabolites is presented in Table S2 (Supplemental materials).
The magnitude of changes between the two groups was compared to flux changes between the same two conditions predicted by the model. A positive and statistically significant correlation between predicted fluxes and experimental metabolomics data was observed. For example, the model predicts an increased carbon flux into the pool of amino acids under methanol growth, and a slight drop of the RuMP cycle. These predictions correspond to the directions of changes for metabolites in these pathways measured by biochemical profiling. Calculations of Spearman's index, which is proper to use in the presence of a small number of observations, between these data sets have confirmed the statistically significant correlations between model predictions and experimental measurements (Fig. 1C,D).
The comparison also highlighted a set of metabolites whose changes contradict the model predictions, including intermediates of the TCA cycle, pentose phosphate pathway, and pyruvate (Fig. 1C). Additional modifications of the model and in-silico simulations were carried out. Changing the phosphoketolase reaction directionality (from fructose 6-phosphate → acetylphosphate + erythrose 4-phosphate to acetylphosphate + glyceraldehyde phosphate → xylulose 5-phosphate) improved model predictions for intermediates of the pentose phosphate pathway. We also considered a complete TCA cycle for methanol and methane growth, which reduced the Spearman's index between metabolomics profiling and flux ratios to 0.9 (p-value = 2.7E-24). Simulations with a complete TCA cycle for methanol growth and an incomplete (knockout of 2-oxogluterate dehydrogenase) TCA cycle for methane growth led to the Spearman's index of 0.94 (p-value = 1.1E-30). However, changes in fumarate and succinate levels did not correlate with predicted flux changes, so the metabolic reactions producing or consuming those compounds were investigated further. Simulations with an unusual TCA cycle 41 -in which oxoglutarate is converted to glutamate and then to succinate via succinate-semialdehyde dehydrogenase and glutamate carboxylase-and, additionally, oxaloacetate converted to aspartate and then to fumarate via aminotransferase and aspartate lyase resulted in the Spearman's index of 0.99 (p-value = 2E-55).
Finally, a set of enzymatic analyses were set up to validate the updated metabolic network. The activity of phosphoketolase was measured in actively grown 20Z R cells demonstrating that the pathway is active upon growth on methane and methanol (Table 3); however, the highest activity of phosphoketolase is found in methane-grown cells. We also detected relatively high activity of succinate-semialdehyde dehydrogenase, supporting the hypothesis of a branched TCA cycle at the oxoglutarate node.

Concluding Remarks
Here we demonstrated the applicability of the systems approach to analyze and improve our understanding of the methane utilization network in M. alcaliphilum 20Z R . A genome-scale metabolic model of M. alcaliphilum 20Z R was constructed and refined using continuous cell culture parameters.
The flux balance analysis of the metabolic model quantified initial steps of methane oxidation, RuMP cycle, EMP/EDD pathways, TCA cycle; major routes for nucleotide, amino acid, and lipid biosynthesis, and fatty acid metabolism; as well as biomass growth under different conditions of the strain cultivation.
The model analysis has also identified the crucial role of an additional constraint on the O 2 consumption rate. Due to the high plasticity of the redox reactions in the strain, the additional constraint on O 2 -input helped the model to navigate across the metabolic network and correctly reproduce experimentally observed parameters (growth rate and corresponding yields), and thus make accurate predictions.
Furthermore, we show that the integration of non-targeted metabolomics data with metabolic flux predictions can help to refine the central metabolic pathways, by highlighting nodes poorly resolved by in-silico reconstruction. The model simulations suggest that PPi-dependent reactions in the central metabolic pathways might contribute to energy preservation upon oxygen limitation, as well as improve carbon conversion efficiency. Global, non-targeted, metabolomic profiling combined with in-silico simulations uncovered a set of highly interconnected nodes at the level of acetyl-CoA, fumarate, and glutamate. Further enzymatic studies confirmed the carbon flux from acetyl-CoA to xulylose-5-phosphate via phosphoketolase, from glutamate to succinate via succinate-semialdehyde dehydrogenase and glutamate carboxylase, and from aspartate to fumarate via aspartate lyase. The model will serve to inform hypothesis-driven strain engineering strategies, targeting enhanced methane conversion rate and efficiency.

Methods
Metabolic modeling. To build the metabolic framework for 20Z R we used the whole genome sequence of Methylomicrobium alcaliphilum 20Z (NC_016112.1) 15 and the previously published metabolic reconstruction for the closely related species, Methylomicrobium buryatense 5 G(B1) 24 . Each reaction in these metabolic pathways was checked for mass and charge balances except reactions coupled with electron transfer system. The model was developed in a format compatible for FBA [42][43][44] . Moreover, to standardize reaction and metabolite IDs in our model according to BiGG Models ID Specification and Guidelines 45 , we used iAF692 model data 46 for almost all of entities and KEGG abbreviations for enzymatic reactions, for which we could not find corresponding BiGG reactions. To link genes with reactions, the NCBI Reference Sequence annotation, NC_016112.1, was applied. To mathematically model the methane utilization network and to solve FBA optimization problems we employed GNU Linear Programming Kit (GLPK) (http://www.gnu.org/software/glpk/) solver in MATLAB using the wellknown COBRA toolbox designed for flux balance model reconstruction and analyses 47 . To visualize metabolic networks and obtained fluxes we employed the Escher web-tool ( Figure S1  20 ml/L of phosphate solution (5.44 g KH 2 PO 4 ; 5.68 g Na 2 HPO 4 ) and 20 ml/L of 1 M carbonate buffer. The following methane mixtures were used for bioreactor studies: (i) 5% CH4: 3.5% O 2 to represent oxygen-limiting conditions; (ii) 5% CH 4 : 5%O 2 , to represent optimal growth; and (iii) 2.5%CH 4 :10% O 2 to represent methane-limiting conditions. Methane mixtures were directly purged into bioreactor culture at 0.1 sL h −1 rate. In batch cultures, methane (99.9%, Airgas) was injected into vials to represent 10-20% of the headspace. For methanol growth cultures were supplemented with 0.5% methanol.
Cultivation. Culturing was carried out in either closed vials (50 ml cultures in 250 ml vials, or 200 ml cultures in 1 L bottles) with shaking at 200 r.p.m. or bioreactor cultures (fed-batch or turbidostat). Two types of bioreactors were used: 1) a DASbox mini bioreactor (0.5 L working volume; 200 ml culture) with two individual bioreactor units, each having automatic temperature, pH, and DO controls, a sample port for measuring OD, and a coupling to a BlueSens sensor system for simultaneous measuring off-gases (CH 4 , O 2 , and CO 2 ); or 2) a 2.7 L bench top BioFlo 110 modular bioreactor (New Brunswick Scientific, Edison, NJ, USA).
Batch cultures were grown in 125 ml, 250 ml, or 1.2 L bottles. The headspace:medium was set at 4:1 ratio. CH 4 , O 2 , and CO 2 composition of the headspace were determined using an SRI GC system. Control set of samples include only ininoculated media. The methane and oxygen consumption and CO 2 production rates were calculated by estimating decline(or increase) of the corresponding compounds over time. The data were analyzed to assess yield (Y), growth rate, and O 2 /substrate ratios. In addition, samples (15 ml) were taken for metabolomics.

ATP maintenance experiments.
To estimate growth-dependent and non-growth associated ATP maintenance cells of 20Z R were pre-grown to OD = 0.36 ± 0.05 and transferred to new vials supplemented with various amounts of methane added to headspace (0, 0.02, 0.04, 0.08, 0.1, 0.2, 0.4, 0.8, 1.6 mmol). Each experiment was represented by 3-6 biological replicates. Vials with cultures were placed back onto shaking platform. OD and CH 4 , O 2 and CO 2 composition were measured every 30 min -1 h over 12 hours. The growth rates of methanotrophic bacteria correlated with the initial concentrations of substrate added. No growth was observed in cultures supplemented with methane below 3.6 mmol g −1 CDW h −1 .
Methane utilization for maintenance and growth is mathematically expressed as: 4 4 , where m CH 4 is the methane requirement for non-growth-associated maintenance activities, Y max is the maximum biomass yield including growth-associated maintenance, and μ is the growth rate 39 . These ATP requirements were determined from a plot of methane uptake versus optimal computed growth rate, given the experimentally determined values for m CH 4 and Y max . m CH 4 has been observed to have value equals to 3.6 mmol CH 4 g DCW −1 h −1 , while Y max has a value 0.6 gDCW per g of CH 4 . Using the assumed ATP yield Y ATP, max = 6 mol ATP per one mol CH 4 consumed and the y-intercept (3.6 mmol CH 4 g DCW −1 hr −1 ), the non-growth-associated ATP requirement was calculated as 21.6 mmol ATP g DCW −1 hr −1 . To estimate growth-associated ATP maintenance the model was simulated by varying only the CH 4 consumption rate upon different values of growth-associated ATP maintenance as a stoichiometric coefficient in the biomass equation ('BOF' in Supplemental Material, Table S1). Figure 2B demonstrates agreement between the equation and the flux balance approach and the corresponding growth-associated ATP requirement is the range of 90-110 mmol ATP g DCW −1 .
Dry cell weight measurement. Cultures (150 ml) from bioreactors were centrifuged to collect the biomass. After careful removal of the liquid phase, tubes of known weight with biomass were weighed (to obtain wet-cell biomass weight), lyophilized overnight using a Labconco freeze-dry system and weighed again. The observed DCW parameters were as follows: 1 L of cell culture with OD = 1 corresponds to 0.336 ± 0.025 g DCW).
Biomass composition. Samples of lyophylized biomass were submitted to AminoAcids (aminoacids.com) for complete amino acid analys, and to Matrix Genetics (Seattle, WA, http://matrixgenetics.com) for FAME derivatization and GC-MS analysis. The intracellular concentration of glycogen was determined using anthron reagent 49 ; ectoine concentration was determined using an HPLC assay 50 . For glycogen measuremnts 20z R biomass (dry) was lysed by incubation of with potassium hydroxide (10% KOH) at 95 °C for 1 h. After neutralization and centrifugation (14000 rpm @ 15 min, 4 °C), ethanol was added to the supernatants (two volumes) and samples were refrigerated overnight. The precipitated glycogen was resuspended in 1 ml of DI H 2 O and centrifuged again (14000 rpm, 15 min, 4 °C). 1 ml of 0.2% (w/v) anthrone reagent (95% sulfuric acid) was added to 0.2 ml resuspended glycogen and incubated for 15 minutes, after which the sample was measured using a Janeway 632OD spectrophotometer and 620 nm. The concentrations of glycogen were estimated by comparing to glycogen standards prepared in the same manner as cell samples. Extracellular polysaccharides were measured in samples of broth using a sulfuric acid assay 51 . Extracellular concentrations of formate, acetate, and lactate were measured as previously described 16 . Non-targeted metabolite profiling. Metabolomic analysis of cells and spent supernatant from cultures of the M. alcaliphilum 20Z R was performed using Metabolon's untargeted global biochemical profiling platform [http://www.metabolon.com/technology.aspx] and according to the published protocol 52 . The data represent accurate relative quantification only for individual compounds across all the samples which contain it. The original scale data represent raw ion counts for the integrated peaks specifically representing each compound. The scaled data is generated by simply dividing each value for a compound by the median (among detected samples), then imputing any null values with the minimum detected value. Thus, the scaled data are centered around 1.0 (median) for all compounds, but the variation among the samples is not affected 53,54 . Statistical analyses were performed on natural log-transformed data.
SCIeNtIfIC RepoRTs | (2018) 8:2512 | DOI:10.1038/s41598-018-20574-z Statistical comparison between predicted flux changes and biochemical profiling. Statistical comparison between methane-grown cells (methane batch) and methanol culture condition was initially conducted to reveal key differences in metabolic signatures. The magnitudes of changes between two groups were compared to predicted flux changes between the same two conditions. An enzymatic reaction with the highest value of flux for biosynthesis of corresponding metabolite was used to calculate flux ratio between two growth conditions. In order to estimate Spearman's correlation index, predicted flux ratios and metabolic signatures were ranked according to the method 55 . 3D plots reflected metabolomics and model ranks taking into account metabolic pathways for measured substances were generated with assistance of Plotly tool (Plotly Technologies Inc. Collaborative data science. Montréal, QC, 2015; https://plot.ly).
Isolation of soluble protein fraction. Cell lysates were generated from biomass cultivated under above described conditions. Cells were pelleted via centrifugation for 1 min at 13,000 g and culture supernatants were discarded. Cell pellets were quenched in liquid nitrogen, thawed and solubilized on ice in 1 mL of lysis buffer (50 mM MOPS, pH 7.4, 50 mM NaCl, supplemented with 1X cOmplete Protease Inhibitor Cocktail solution (Roche Diagnostics Corporation, Indianapolis, IN)). The cells were then sonicated on ice at 4 °C, at 90% power setting for 30 seconds × 3 cycles, with a one minute cool-down period between sonication cycles using a Braun-Sonic-L ultrasonicator. Lysates were cleared via centrifugation at 16,000 × g at 4 °C for 5 minutes, and the supernatants were isolated for use in subsequent enzymatic analysis.
Enzyme activity determination. Cell lysates were generated from biomass cultivated under above described conditions. Cells were pelleted via centrifugation for 1 min at 13,000 g and culture supernatants were discarded. Cell pellets were quenched in liquid nitrogen, thawed and solubilized on ice in 1 mL of lysis buffer (50 mM MOPS, pH 7.4, 50 mM NaCl, supplemented with 1X cOmplete Protease Inhibitor Cocktail solution (Roche Diagnostics Corporation, Indianapolis, IN)). The cells were then sonicated on ice at 4 °C, at 90% power setting for 30 seconds × 3 cycles, with a one minute cool-down period between sonication cycles using a Braun-Sonic-L ultrasonicator. Lysates were cleared via centrifugation at 16,000 g at 4 °C for 5 minutes, and the supernatants were isolated for use in subsequent enzymatic analysis.
Phosphoketolase activity was measured using a hydroxymate assay that detects acetyl phosphate 35,57,58 . Briefly, 100 μL reactions in reaction buffer [30 mM potassium phosphate (pH 6.5), L-cysteine hydrochloride (2 mM), sodium fluoride (20 mM), sodium iodoacetate (10 mM), thiamine pyrophosphate (1 mM), and fructose 6-phosphate (30 mM) (Sigma)] were initiated by the addition of 25 μg total protein. Reactions were incubated at 37 °C for 1 h. After incubation, 50 μL of 2 M hydroxylamine hydro-chloride (pH 7.0) was added to the reaction and incubated at room temperature for 10 min. Reactions were then terminated by adding 300 μL of a 1:1 mixture of 2.5% ferric chloride in 2 M hydrochloric acid and 10% trichloroacetic acid. A color change due to ferric-hydroxamate formation was measured at A505nm, and the concentration of acetyl-P formed during the initial reaction was calculated by regression analysis compared to known standards.

Availability of data and material.
Developed genome-scale model (xls) is available on the web-site: http:// sci.sdsu.edu/kalyuzhlab/.