Characterization of the Genome-scale Metabolic Reconstruction of P. freudenreichii
Structural characteristics of the computational model (iMR558)
The in-silico metabolic network of P. freudenreichii includes 558 genes representing nearly 23% of the strain’s open reading frames (ORFs) and two percent of its genome (see Table 1). The metabolic network also accounts for 648 unique metabolites participating in 715 elementary and charge balanced reactions (see Table 2). While 606 reactions represent intracellular bioconversions, 109 are transporters responsible for translation of the metabolites across the cell membrane. Of the total number of reactions, 673 (94%) received gene-protein-reaction (GPR) associations either based on the annotation data or in the process of model curation and refinement. An additional 41 reactions without GPR associations were also incorporated based on experimental evidence on their presence and/or their essentiality for supporting growth. Metabolic reactions are distributed between two cellular compartments, including the cytoplasm and the extracellular space. Two additional demand reactions were added into the network to account for intracellular accumulation of cobalamin and riboflavin. The computational version of the metabolic network also includes 95 exchange reactions (one per extracellular metabolite) to define the boundaries of the network with its in silico environment.
Table 1
Main features of the P. freudenreichii DSM20271 genome sequence by Koskinen et al.15 Other Collection Numbers: Niel 1928T, ATCC 6207, CCM 1857, NCDO 564, NCIB 5959, NCTC 10470, P104 (Iowa State University Collection);
Sequence name
|
RM25
|
Accession number
|
GCA_000940845.1
|
Genome size (bp)
|
2,649,166
|
Predicted genes
|
2353
|
Protein coding genes
|
2320
|
GC%
|
67.34
|
Table 2
Characteristics of the computational metabolic model of P. freudenreichii (iMR558);
Total Genes
|
558
|
Total Pathways
|
63
|
Total Metabolic Reactions
|
715
|
Intracellular
|
606
|
Transport
|
109
|
Exchange
|
95
|
Spontaneous
|
4
|
Total Metabolites
|
745
|
Intracellular
|
648
|
Extracellular
|
98
|
Dead-end metabolites
|
225
|
iMR558 contains 184 reactions that were added during the gap-filling process in the metabolic network. Indeed, these reactions are necessary to consume carbon sources and unblock the production of biomass components and other known byproducts of metabolism, including vitamins. During the gap-filling process, 55 cytoplasmic reactions (associated with 21 genes) were added to the model and annotated through homology investigations for the responsible coding sequence. These reactions were mostly associated with lipid metabolism and biosynthesis of vitamin B12 precursors. Eighty-eight transporter reactions responsible for protein-facilitated transport of the metabolites across the cell membrane were annotated through extensive search in biochemical databases and homology predictions. An additional 11 reactions were incorporated into the model to account for spontaneous bioconversions, diffusion of metabolites through the membrane, carbon dioxide equilibrium, nongrowth-associated ATP and biomass equations. The final version of the metabolic network still contains 31 orphan reactions, representing the current knowledge gap regarding the metabolism of P. freudenreichii at the genetic level. Among orphan gap-filler reactions, 16 were transporters, two were required for ubiquinone biosynthesis, five for cell wall biosynthesis, six for lipid biosynthesis and two for the biosynthesis of vitamin precursors.
iMR558 includes seven high-level metabolic subsystems (see Fig. 1). The largest group is amino acid metabolism, accounting for 24% of total reactions, followed by carbohydrate metabolism with 18% and coenzyme metabolism with 17% of all reactions. The relatively high percentage of reactions in coenzyme groups correlates with the known capability of the strain to produce various vitamin compounds. The complete list of reactions, metabolites, and the GPR associations included in the reconstruction together with the SBML version of the model is provided in supplementary materials 1 and 2. Details about manually annotated genes included in the model are provided as supplementary material 3.
Biomass composition
A well-characterized biomass composition of B. subtilis was adopted as a template for the development of the biomass composition of P. freudenreichii DSM20271. The biomass equation was adjusted to account for some experimentally identified specifications of P. freudenreichii, including DNA and fatty acid compositions. DNA composition was obtained from DSM20271 genome sequence data. The compositions of RNA, protein, lipids and the cell wall were assumed to be similar to those of B. subtilis.
The fatty acid composition was adapted to experimentally reported data39–41. P. freudenreichii accumulates various types of fatty acids for the assembly of the cytoplasmic membrane. Total lipid analysis of DSM20271 has revealed that almost 97% of this subspecies is composed of saturated branched and unbranched fatty acids C15:0 and C17:041. While anteiso-C15:0 fatty acid is the most abundant fatty acid in this strain, anteiso-C17:0 and n-C17:0 also have major contributions to the biomass lipid composition39–41. The molecular composition of the lipid species in P. freudenreichii biomass, including glycerophospholipids, including cardiolipin, phosphatidylethanolamine, phosphatidylglycerol, and phosphatidylinositol, were hence adjusted to represent the experimentally reported data.
A previous study reported that glutamic acid, alanine, meso-isomers of diaminopimelic acid, muramic acid, and glucosamine constitute the peptidoglycan of P. freudenreichii42. By inspecting the genome, we were able to identify genes putatively responsible for the biosynthesis of peptidoglycan, comprised of N-acetylglucosamine and N-acetylmuramic acid disaccharides, and the corresponding reactions were included. While the overall composition of the cell wall was adopted from B. subtillis, the peptidoglycans are hence specific to P. freudenreichii.
Lipoteichoic acid was excluded from the biomass equation because of a lack of genetic and biochemical evidence in P. freudenreichii43.
P. freudenreichii uptake and secretory system
Detailed characterization of the uptake and secretory system of an organism is key to developing an accurate constrained-based metabolic reconstruction of its metabolism. Despite the experimental data on the capability of P. freudenreichii to catabolize various types of sugars, the membrane transporters of several saccharides and alcohols were not explicitly identified from the DSM20271 genome. Based on homology modeling, Loux et al explored the uptake mechanism of a number of carbohydrates in several Propionibacterium subspecies10. While their study did not include strain DSM20271, it revealed that the subspecies of the P. freudenreichii genera have diverse stain-specific sugar utilization capabilities10. To identify DSM20271-specific transporters, databases including TCDB, TransportDB, KEGG and UNIPROTKB were extensively searched followed by rigorous BLAST analysis. As a result, a total of 88 proteins were reannotated as putative transporters of P. freudenreichii DSM20271.
The secretion of propionate and acetate transporters is a characteristic of P. freudenreichii. However, there is a lack of knowledge about the mechanism of propionate transport across the cell membrane. The previous genome-scale metabolic network of Bacillus subtilis includes a proton symport reaction for these organic acids. As we could not find a hypothetical protein for their transport based on BLAST, we used the same mechanism as in B. subtilis.
Curation of P. freudenreichii carbohydrate metabolism
A majority of pathways relevant to catabolism of carbohydrates can be annotated from the strain genome. However, some carbohydrates, such as galactose, adonitol and glycerol, could not be catabolized by the initially developed metabolic network due to crucial gaps in the corresponding transport/bioconversion routes. Using homology analysis together with literature data, the gaps were eliminated, and the final model gained the ability to grow on these substrates. Nonetheless, despite our efforts, the gaps in catabolism of erythritol, salicin and esculin could not be filled due to the limited knowledge of their degradation pathways and unfruitful cross-species sequence similarity search.
Curation and refinement of central and energy metabolism
The process of high-quality metabolic reconstruction of an organism entails detailed characterization of the (central) metabolic and bioenergetic pathways and integration of the relevant data from genomic and biochemical references into the reconstruction. Many enzymes in the central metabolism of P. freudenreichii have been previously characterized44,45. In line with experimental evidence46–48, the genomic data confirm the activities of the Embden-Meyerhof pathway (EMP), both oxidative and nonoxidative branches of the pentose phosphate pathway (PPP), and the fermentative pathways leading to the synthesis of propionate and acetate.
Consistent with its facultative anaerobic nature, P. freudenreichii has a structurally complete TCA cycle45,48−50, which may justify the observed oxidative phenotypes with free oxygen and the ability of the strain to grow under aerobic conditions51–54. Genome annotation data suggest that strain DSM20271 contains all the enzymes of the electron transport chain essential for both anaerobic and aerobic metabolism, including NADH dehydrogenase/complex I (NDH-1) (NADH: ubiquinone reductase), membrane-anchored fumarate reductase/succinate dehydrogenase (SdhA, SdhB, SdhC), cytochrome bd complex (cydA and cydB), cydC subunit, which is required for the synthesis of the functional form of cytochrome bd, ATPase (atpA, atpB, atpC, atpD, atpE, atpF, atpG, atpH), and all genes indispensable for biosynthesis of heme, an essential cofactor for the cytochrome bd complex. In addition, transcriptomic and proteomic analyses have revealed that all components of the electron transport chain, including NADH dehydrogenase/complex I (NADH: quinone reductase) and fumarate reductase and ATPase47,54, are expressed in P. freudenreichii under both anaerobic and aerobic conditions.
The Wood-Werkman cycle (WWC) is a key part of central and energy metabolism in P. freudenreichii. All the essential enzymes for carbon flow through WWC were identified from genome annotation data with the exception of propionyl-CoA:succinate CoA transferase (2.8.3.-). This enzyme is key to the production of propionic acid through this pathway, and biochemical evidence supports the presence of this activity during P. freudenreichii growth on glucose55. By homology query, we were able to identify a putative corresponding gene with 46% sequence similarity with scpC of Escherichia coli (strain K12) and embedded it in the model.
P. freudenreichii has shown the ability to grow on minimal amino acid and peptide-free media46,56. Hence, the strain must be capable of de novo synthesis of all essential amino acids. This capability is represented in the strain’s genome, where most genes required for biosynthesis of amino acids from the central metabolism are present. However, our draft network, which was reconstructed based on genome-annotation data, was not able to produce some amino acids. The gaps were resolved by identifying and adding some missing spontaneous reactions in amino acid metabolism, restoring the model capacity in de novo amino acid biosynthesis.
Secondary metabolites biosynthesis
Many of the value-added bioproducts are secondary metabolites. Reconstruction of the genome-scale metabolic network allows for a detailed understanding of the secondary metabolic pathways and how they stem from and are connected to the primary metabolism, which in turn inspires bioproduction improvement strategies. Corresponding to the documented capability of P. freudenreichii in producing various types of vitamins, the complete biosynthetic pathways for several vitamins, including vitamin B2 (riboflavin), vitamin B3 (niacin) and vitamin B12 (cobalamin), which are known to be produced by P. freudenreichii, were extensively curated and incorporated into iMR558. In addition, while the ability of P. freudenreichii to produce vitamin B6 (pyridoxine) and vitamin B9 (folic acid) is also evident, this capability is not reflected in the genomic data due to the limited understanding of these pathways at the genetic level, and our extensive BLAST study did not improve the relevant genetic resolution. Furthermore, although the vast genetic gaps in the biosynthesis pathways of vitamin B5 (pantothenate) and vitamin B7 (biotin) vitamin B1 (thiamine) comply with the essentiality of their presence in the medium for growth 56,57 we conducted a BLAST search to complete their biosynthesis pathway as much as possible.
The P. freudenreichii pathway for the biosynthesis of cobalamin has been the subject of several studies58–60. However, the biosynthetic pathway of cobalamin in current biochemical databases includes many gaps and topological inconsistencies that prevent it from being functional in a computational model. On the other hand, the literature on the characterization of the enzymes involved in cobalamin biosynthesis and details about the cofactor specificity of these reactions are scarce. Therefore, an extensive effort was made to curate the vitamin B12 biosynthesis pathway by collecting together pieces of information from various resources.
The cobalamin coenzyme (alternatively named vitamin B12 coenzyme) comprises two covalently linked structures called upper and lower ligands. The upper ligand, adenosine-GDP-cobinamide, is composed of a cobalt ion-containing corrin ring whose de novo biosynthesis pathway in P. freudenreichii is initiated by the conversion of L-glutamate to aminolevulinate through the C5 pathway61. The lower ligand is N1-(alpha-D-ribosyl)-5,6-dimethylbenzimidazole, which is biosynthesized from reduced flavin mononucleotide (FMN) by the Blub/CobT2 gene62,63. The bioavailable form of cobalamin coenzyme, which is essential to the bioactivity of many human cells, is the result of covalent binding between dimethylbenzimidazole and adenosyl cobinamide64. There are also other analogs of cobalamin --pseudovitamin-B12–, sharing the same upper ligand bound with various lower ligands. These metabolites are used in some bacterial species as cofactors for the activation of cobalamin-dependent enzymes65,66, including methyl-malonyl-CoA mutase67, an enzyme in the propionate biosynthesis pathway. (See Fig. 2.)
As seen from Fig. 2, the biosynthesis of the upper ligand (adenosine-GDP-cobinamide) is initiated from L-glutamate followed by the assembly of a corrin ring and the subsequent insertion of a cobalt ion into the structure (though an anaerobic pathway) and ultimately joining of an aminopropanol arm to the assembly59. Most of the enzymatic steps in this pathway are predicted by the strain’s genome annotation with the exception of CbiE (EC: 2.1.1.289) and CbiT (EC: 2.1.1.196). Our BLAST search identified two ORFs putatively responsible for these steps, including RM25_0741 and RM25_0737, both sharing 34% sequence identity with the bifunctional fused protein CbiET from Bacillus megaterium.
After adding these GPRs, the pathway of adenosine-GDP-cobinamide biosynthesis seemed to be complete, and it was not able to carry any flux during the model simulation. Our gap analysis revealed two structural deficiencies leading to the blockage of flux through this pathway: 1) inability of the network to metabolize or exert two byproducts of the enzymatic reactions, including S-adenosylhomocysteine and acetaldehyde, and 2) a gap in the enzymatic steps for biosynthesis of the aminopropanol arm, which is attached to adenosyl cobyrinate hexamide to constitute adenosyl-cobinamide phosphate. Both of these network deficiencies lead to violation of the pseudo steady-state assumption governing the standard FBA and thereby nonfunctionality of the pathway in silico.
Acetaldehyde dehydrogenase, which catalyzes the conversion of acetaldehyde to acetyl-CoA, has already been reported to be active in P. freudenreichii68. Our similarity search identified both RM25_1259 and RM25_0854 as homologs of acetaldehyde dehydrogenase of Escherichia coli (strain K12) with 35% amino acid sequence identity (e-value 9.52E-62), leading to elimination of the corresponding gap from the network.
Several enzymatic steps in the biosynthesis of adenosine-GDP-cobinamide account for the conversion of S-adenosyl-L-methionine (SAM) to S-adenosylhomocysteine (SAH). While SAM can be supplied by L-methionine at the expense of one ATP (EC: 2.5.1.6), the fate of the produced S-adenosylhomocysteine is not represented in the genome annotation data. Our similarity search identified a weak homology between RM25_1277 and ahcY (EC: 3.3.1.1) of Thermoplasma volcanium (sequence identity = 34%), which, if valid, can restore the required SAM/SAH balance for the biosynthesis of adenosine-GDP-cobinamide. Due to the lack of better evidence, we filled the network gap accordingly.
As mentioned above, the biosynthesis of adenosine-GDP-cobinamide requires attachment of an aminopropanol arm to adenosyl cobyrinate hexamide to enable the formation of adenosyl-cobinamide phosphate. There are two forms of aminopropanol arm, i.e., (R)-1-Aminopropan-2-ol and (R)-1-aminopropan-2-yl-phosphate, which are synthesized from glycine and threonine, respectively. There is biochemical evidence that (R)-1-aminopropan-2-ol cannot be metabolized in P. freudenreichii DSM2027169, which is consistent with our unfruitful BLAST search. On the other hand, one of the two enzymatic steps for the synthesis of (R)-1-aminopropan-2-yl-phosphate (L-threonine 3-O-phosphate decarboxylase; EC: 4.1.1.81) has already been predicted by genome annotation data. We were also able to identify RM25_1014 as the putative encoder of L-threonine kinase (EC: 2.7.1.177) with 33.5% similarity with the bifunctional L-threonine kinase/L-threonine-O-3-phosphate decarboxylase (MM_2060) of Methanosarcina mazei Gö170. Filling this gap accordingly enabled the vitamin B12 pathway to carry flux during FBA simulation. Guided by the COBRA method, therefore, we were able to develop the first evidence-based functional pathway of vitamin B12 biosynthesis.
P. freudenreichiimetabolic model simulation and validation
Carbon source utilization
P. freudenreichii has demonstrated the capability to ferment a wide range of sugars. Loux et al10 examined the ability of 21 subspecies of P. freudenreichii to utilize 50 different sugars. While all the strains were capable of degrading glucose, glycerol, mannose, galactose, inositol, erythritol, adonitol and esculin, the ability to degrade lactose, xylitol, melibiose, gluconate, ribose, D-fructose, arbutine, L-arabitol and L-arabinose was variable among the strains, indicating the strain-specific profile of carbohydrate fermentation in the P. freudenreichii genus. As such, and because of the absence of P. freudenreichii DSM20271 from the study of Loux et al10, we gathered the carbon source utilization data of the strain from the scarce literature and applied them in our model curation process.
Simulation of carbon source utilization was carried out using flux balance analysis (FBA). Table 3 compares the carbon source utilization profile predicted by iMR558 with the results of the wet-lab experiments. With the exception of erythritol, the in silico predictions were found to match the experimental data.
Table 3
Carbohydrate fermentation ability of P. freudenreichii DSM20271
Metabolite
|
In silico
|
In vivo
|
Reference
|
D-Glucose
|
+
|
+
|
15,99
|
D-Fructose
|
+
|
+
|
15,99
|
Trehalose
|
-
|
-
|
99
|
D-Xylose
|
-
|
-
|
99
|
D-Ribose
|
+
|
+
|
99
|
Cellobiose
|
-
|
-
|
100
|
Adonitol
|
+
|
+
|
15,99
|
L-Arabinose
|
+
|
+
|
101,102
|
L-Rhammmose
|
-
|
-
|
99,100
|
D-Mannose
|
+
|
+
|
15
|
Salicin
|
-
|
-
|
99
|
Sucrose
|
-
|
-
|
100
|
Starch
|
-
|
-
|
100
|
Glycerol
|
+
|
+
|
15,71,99
|
D-Mannitol
|
-
|
-
|
99
|
L-Sorbitol
|
-
|
-
|
100
|
Lactose
|
-
|
-
|
15,100
|
Melibiose
|
-
|
-
|
99
|
Galactose
|
+
|
+
|
15,99
|
Sucrose
|
-
|
-
|
100
|
L-Lactate
|
+
|
+
|
56
|
Maltose
|
-
|
-
|
99
|
Inositol
|
+
|
+
|
15,99
|
Erythritol
|
-
|
+
|
15,99
|
Our in silico strain was not able to produce biomass by using glycerol as the sole carbon source. While corresponding phenotypes were not reported from wet-lab experiments, the same behavior was observed from P. acidipropionic glycerol-limited culture35. In agreement with experimental data, however, our model was able to utilize glycerol for growth on complex culture medium containing amino acids and/or additional carbon sources such as glucose6,71. P. freudenreichii was also able to grow in silico in the presence of adenosine as described by Glatz and Anderson56 and in a medium containing nicotinic acid, cysteine and tryptophan defined by Crow72.
Can P. freudenreichii utilize amino acids as nitrogen or sulfur sources? The literature does not present any direct evidence of such capability. However, P. freudenreichii is able to grow on YEL medium without additional nitrogen or sulfur sources47,73, implying that the strain should be able to derive its required nitrogen and sulfur from the amino acids present in yeast extract. This hypothesis is corroborated by our simulations indicating that some amino acids, including asparagine, alanine, aspartate, glutamate, serine, threonine, cysteine and glutamine, may be metabolized as nitrogen source alternatives to ammonium. Our simulations also identified cysteine as the only amino acid-based source of sulfur for P. freudenreichii.
Growth and byproduct secretion
Well-curated strain-specific GEMs can be used to explore the metabolic capabilities of relevant organisms. Insights from such studies can uncover the potential of a strain to be applied in value-added bioprocesses. Hence, we used iMR558 as an in silico representative of P. freudenreichii to evaluate its metabolic capabilities, including biomass and organic compound production. Flux balance analysis (FBA) and flux variability analysis (FVA) were used to simulate microbial growth and estimate flux distribution over the metabolic network. To enable validation of the simulation results, the in silico culture conditions matched those reported from the experimental growth cultures74. In doing so, the model was constrained to the culture conditions applied in the study of de Vries et al74. Hence, the glucose uptake rate (GUR) was set to 1 mmol/gDW/hr, and tryptophan and other amino acids present in casein hydrolysate, including glutamine, glutamate, proline, lysine, leucine, and aspartate75–77, were allowed to be untaken up to 0.05 mmol/gDW/hr.
The goal of the simulation is to predict the growth rate and byproduct secretion in P. freudenreichii under different culture conditions compared with the experimental data. However, microbial growth and metabolic phenotypes are highly influenced by the ATP required for maintenance functions (ATPM). Equally important, the demand for ATPM is not constant and can vary depending on the limiting substrate and other nutrients present in the medium. Given that the ATPM data for P. freudenreichii are not available for any growth culture, we calibrated the model by setting the ATPM at the values that corresponded to the experimentally observed growth yield values74. The calibrated model was then used to predict the yields of acetate and propionate as the major byproducts of P. freudenreichii. The simulation results in comparison to the experimental data are presented in Table 4.
Table 4
Predicted molar growth yields and byproduct yields of P. freudenreichii in synthetic culture with glucose and glycerol (optimal, anaerobic cultivation)
Substrate
|
Growth (gDW/mol substrate)
|
Acetate yield (mol/mol substrate)
|
Propionate yield (mol/mol substrate)
|
ATP maintenance (in silico)
|
Reference
|
In vivo
|
In silico
|
In vivo
|
Most Probable value (in silico): min and max FVA
|
In vivo
|
Most Probable value (in silico): min and max FVA
|
Glucose
|
65
|
65
|
0.52
|
0.53
|
1.06
|
1.03
|
0.23
|
74
|
Glycerol
|
24.1
|
24.1
|
0.00
|
0.07–0.07
|
1.07
|
0.93
|
0.76
|
74
|
ATPM values of 0.23 and 0.76 were obtained for the strain’s growth on glucose- and glycerol-limited cultures, respectively. Although no relevant experimental data for P. freudenreichii are available, the corresponding experimental values for a strain of the same family, i.e., Propionibacterium acidipropionici have been reported to be 0.76 and 0.36, respectively35. Given that ATPM can vary by several orders of magnitude in various organisms (E. coli (7.6 mmol/gdw/hr)78, B. subtilis (10 mmol/gdw/hr)32 and P. putida (12 mmol/gdw/hr)79), the identity of the order of magnitude of the predicted ATPM with the corresponding experimental values from close strain points to the consistency of the model prediction.
The calibrated model also predicts the yields of acetate and propionate on both glucose- and glycerol-limited cultures very accurately. The ratio of propionate to acetate in glucose culture was also calculated to be 1.94, falling within the range of experimentally reported values (1.5–3.6)50,80,81.
To examine the robustness of the model predictions, we also calculated the variability range of acetate and propionate production using FVA under optimal growth conditions. Interestingly, both the minimum and maximum fluxes through the organic acids coincided with the FBA results. The identity of both results suggests that not only the strain’s central metabolism operates under optimal growth conditions but also the production of acetate and propionate in P. freudenreichii is growth-associated.
Batch culture simulation
Whereas the FBA calculates the metabolic fluxes under steady state assumption, the industrial strains are mostly cultivated under batch/semibatch conditions during which the concentrations of the medium constituents constantly change. Dynamic flux balance (dFBA) is a modified version of FBA developed to enable simulation of batch culture using constraint-based metabolic models82–85. Herein, we used dFBA to predict the kinetics of growth11 and flux distribution in P. freudenreichii using iMR55846. To enable comparison with experimental data, our in silico culture conditions were adapted to the study of Meurice et al., who conducted anaerobic batch growth experiments to evaluate P. freudenreichii propionate and acetate yields on glucose as the sole carbon source46. Previous studies have shown that the decreased culture pH due to the secretion of organic acids can influence P. freudenreichii growth86–88. In the Meurice et al. experiment, however, pH remained relatively constant during the first 100 hours of fermentation, and thus, we conducted the simulation within the same time frame to enable in silico-in vivo comparison, independent of the pH effect.
These initial parameters for dFBA simulation were calculated from the experimental dataas follows: inoculum density (ICD) of 0.0044 gr and initial glucose concentration (IGC) of 28.668 mmol. Glucose uptake flux was calculated by equations (1) and (2)89 and specifically calibrated to enter the system with the corresponding value of -1 mmol/gDW/hr. rglucose, S, µ, X and Δt correspond to the glucose uptake flux, glucose molar concentration, specific growth rate, biomass concentration and fermentation time interval, respectively.
![](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAToAAAA9CAYAAAA9MKVQAAASJElEQVR4Ae2ajZXUOgxGaYEaaIEeKIEaaIEO6IAOqIAKaIAG6IAe9p27jwtarexkfrKTycrnhExsWZY/SZ+VLG8eujUCjUAjcHAE3hx8f729RqARaAQemug6CBqBRuDwCDTRHd7FvcFGoBFoousYaAQagcMj0ER3eBdftsGfP38+fPz48eHz58+XKerZjcAVEfj69evDhw8fHn78+LFKaxPdKphep9CXL18e3r59+/Dt27cSgN+/fz+OQYQEHdenT58e+xjrdh0EOGw4aMSYO76h/zU3SO7du3ePMbeEQxPdEkKvdJzEguRGyUQ/QcYFERJ03Jnz5s2b4bxXCufZ28YP4MkBAsZcHCz0cX/tjQP1/fv3i2TXRPfaI6XYP8lEIo0qOaZAcJBartx4paC/2+UIgD9+oHrLDYzButvD46EKTjM8mug6Up4hwKsRRDZqVHMEFlVGbhDfqArMsv08R8DK7devX88EwTgfMs+EXlEHsVgdvELQRCcSfX9EQBKbnY5WfP3qtG3QcOBwoPTBsYzz9+/fH7EaxW0T3TKGdy3B6w8Jk5OF1yH6c7XgNyHIbNSYQwJyzV5vR/PpR0f8wD4K0JmOrcb2YhdVChjzDeo1V29r/QFWxHTVmugqVA7UN3r94dWUwMjNKiL352eIkvlcBOIpiei3J8iNk3hv35s4FCAX9nZLAgZTsMEO/DU7fLJ/jvS81h/6rNr780ivpLrvbhEgUfL3NhKI5KlOP/qZs6ZBWCbi2qqDSo458fsegZyJEjnsg0Rv0VgXLHIl/NK2gIuHD/YsVdAcbJVfz7H7Uh/41nANDNf4Q5ywO7cmuozIgZ4JMJIjkgrb83tG9dc85E9JFNbwJF0zz9expeC36qtsfAkXmTRLa2Ensqdc55B3rKBHlR0Jfqr/Zvub+UASHNmC3tFbw2zN0dgafyhT2dREN0L2AP0Gan79MmkgvNzOSZT4irVEYFaAed38TCIRsLnSy3JbPa/FQTuxde21hNFoT/pt9EcgsMIGbLpGc2+VD4gpMJrtBVtm46fYuMYfTXSnIHog2VH1ZAVWBfCagKogci2Ce9QYO1f/SOcW/dp5TuW1hT3qhHj2gB9xc81qzf2N7mv90UQ3QvDg/VUwcsKSLJBd1RjjqhpzR8lvtVERHYlhEKKbqo5nrnzix1fBvBZjVDNcVi3cIVl0YcOoIYc+16Uioa+qjqxWqop3pP+a/diZcVE/+LGH2JB1X3kMf4iPvsEfYIUsY/nAm/kAvKzKox/jutGW0T628Afrgg+6c6sjOkv1890hQPBWSWEwEOBVc7wakwCqMYiTwM9Jgyx9Jhw2oYdnriwv+SAXiQ57efZ1nETlN/3oIQGZQ19u9GEbNkJeyFvV5j/UMBed6KoSJuve4hlbK9IefVsFQ+SxOR5g7Bt/sl9JCeIBK3Q5J8fCyAfsFV3GgdjTx2XTTuzJ/kVmK3/oU+2I9ya6iMaBfhtsMfAJTIPbRCJoY4NMCNAYuI5LAHmOgV+RjHO5S0ZLBKKN7MFmwnDHPmzhsrlf92U/iY08CR+b8lGH45AfxHCLBjbYiw3uGTvox5dcsV8b8Rfz4uEQ5ZgHBvjAfu4VNuisfOBaszFl0Btjz/6t/DHbC2s30emBg90JeJLVgCPoSG4Cgn4u++LWq4RxXBJEJ4lI4qCH35GUlM93ZXN/fpZQTcg4LkFlIhglHzZibyZX9WTSXkqYaMsWvyEC9obNXNiPn/gNLhUm2MEhg0zlB/fEOPptM1/PfCCmI1sk60i6runca/tj5E/XbaITiYPdCShObwKKgI5BWfXF7ZNoXKPmfPTGxBnJ088cEg2bltpsfQktV49VtWj1UK0paWf7TRjWuXUDX6/ov8ouiQmcc0OHRBnHZpX4zAdLB5akm320pT/wMXaNcGqii54/0G8C+9xkNTFytXMJPAb/kk0SIolbNSuCHNBV8kmK1ZpUSczJzTlgcE9tRkzuKVd71eHAnmc+kKxG/mH+iHS149r+MF4rvfqwiU4kDnS/RlXiq2+ueM6FyQpqiUAkxFwNuC4EDknFZqDn5DOx8prKV5XeiEjjenv7PSMmbB3tCSyryn3mA6vA2SE4It0t/MGBx3qVL6OfmugiGgf5bUDNTrg1W2U+Vc8sqNfoQcbvTEvyVgMVwUpQ+duP+83kaH8kOhLDxK/wGSX/kt23HJeYRn465XBgHzMfVAdWXFfSrYjn2v7Ar8Qn9uYKP/ujiS4jcoBnnE4QLDl/zVbRQ9BmclkzVxnsqJLN8XiXhOgjaSIZmSgzQmMtKzsSkHXVwRj6JQYqX/oY505DHhkae89rPQ7s7B9x0edxP/SxJzHRdPvFBv+CN23mA0nQ12DwiRV21osPPLSu6Q90YSfrrWlNdGtQapmLEDD4c7JVSk1akodANkmQhXBJWhPS+b6q8wrDPJMQ8rKSRJf60Ike5UkaGzKSHetJgI7v8S6+7sf9Y6vkkgkbDKmGuMBIwmPOzAesFedBkBEjfsfxqPeW/mii22PkHswmEycSymyLEFEkOGXpG53gszHmZH0k+kgX/ZlMtWGv99F+7I9k5B5mY+CVMXMeumYYOV6tiY5b+OMv0WEUJ65sbBBwEnBScMpV791ufu3dExN9l7wOrV2v5W6PgD4fJc7tLWwLjo7AX6KDxCh5YXlIiFMY4oPoIEEIcO2JvASaryC5nF6ad8RxSQDM117xdWDvmBA77IvXo26NwK0Q+Et0GsCpS2BCbFsRkR80IdVux0bA19b43ejYO+7d7RGBZ0QHuUF013hNHW2Y0x0i7bYNAmsrw5ZbX0U3VveHVcyuZ0RntbXV9xRfZXhl63ZcBGYfuo+7697ZXhF4RnT+OX4rg/kjB6fjPX1n2goL9B79G92W2LXuRmAtAk+I7txqC/LyDwyQGL95BeavteiMzW82/lWXMapH5zMnf7vzuyHVZtVcn9dhXzGqilE5ZVhr9B2SfkkfeXSjM9qtLdiHba7PfaTXOX1vBBqBl0PgCdHxwZikPqXaIqGZQ6JLUPy3Eckhb0VCkwC50wdZqCuv7396rP7qK3H6HxfRB4FlovOVHHlkuOzLH8q1w/WibCY6ZdkD+3d99t+tEWgE9oHAk2yUoHLij0wl6Ulokjy2UT8yVDsQUdUkNMgjNskRMoxNYoawYoPk4v/Rk4xiH/IQU2U/8yuiog8is7l+3r/zJX7l+34cBPbsW+JxlGN4wLi9J29civcTojNBYzLPwPDVLhshkLkyQ64iFtcYERpOgyBzs3/JXuQySaGLefSz79is9DLhRhl+qzfuH530c+21cRDFK+IX+/nd7SkCYIJv8+H2VOo6TxzMoxikPx/crGqOEde5MNAq4j3HvGO3vo9shmuw+dyYfEJ0p2wSgyqSQAcOYCwbZWXlK2Fer6r2dFwOLMk0V3NZJzZgS56PnGPZ6ayJLc6LRKZ+11cvMuwPhzB35DDn3/IukbO/bKv7ZmzPZH0L/Ihr8Mlxfaotxk6MEWKQ2ImNcdbLZMdz9pvzsJECg/FMhOjDp/g2Xpfux7UvuWObZDbSA24VHiP52H820fmamas2lAtmXIjfIwJkjI0CfiauETmyLvI5CPKaylV2OlYRL1WOhFCB61wDhj0TrOiKFVK2Zy/P+ojgiU28Cbp72Ee0fcvf+HtELqeuK2HGeRUxMZ7JbkZyyKOHQ5fYxce5mWdLeZPnbfVMjMkL5FIuOvK659p/NtGZ6Nxj87TKpxMy9LGZqkmc2QEmZD51XD/3Z90zOe2pKjb1oJ/g4YptpjfK7fW3wZWJzsqiSe6f54gP4jbH+j+J034RdzGhTd7sC7VKdpAXcchz1cgdCwXiFpuzLDL0z2K+0r1FHzaS39hsPEZcRmviB3A4ZQ8164xWCP0menQ+yYGhACngYcpjf0WAyKgvEhc6JKOoZySfZaJcJlAJOdqJ/diem3uK/ZW9jIvB3onC5Ir7Z0+zRIr7f02/wQhcruVTYizmjYf8LHGNtzgv+4BcifkjiUQ59lJVelHmpX6TkxKxxLyG6MzTGLtLNj/P6qUZf8ZNFAxjYQxlYQmkeh3EwRIdDosyOpI+9MHwACHJ0MeYmxMYvpExRsMmxmMwEDwEKesaSOi1z7nMVydr208fstkByrKesuydIMqk+gey3d2wlb3R9JuBtztjNzSIuMCfXnEpfDs6uKNc/B31xX5+GzfcbcTWjHzMNWINf1U+Ql/WQRzrX9dCxhyybw93ccl5NrKNvMcv5t5Izv6ziQ4FlpuAacJLWJUzHAPsSHLoIjjolwwNBEiDPtZgftwYY85BBpAqkkGXlaH6IxkKBnIAGHXy270p511iQycXc7VbmT3f8QF2gwX4Vtjt2f5r2ca+2T9Y5ERjjP4cr9Xaxhm60MM84i7GrDlgH3OQI3aqRh5F32hrzi9iNNuIDLqJUxvPVew7fqu7OGT8R/aIY9zbSJb+s4gOUEdJgWNxTLf9I8DhQuBzcWitaVYrJuqaOTMZAjwn7Uw+jmED87Hp0mZFm3HwMGedWZMQ46Ho3JiMsRJBJ4nNhSw4RLLiOZKc62eyY//IVT7xoGau/mYtZLF1aV+uufUdO4jDtUSn/FrSPpnorAKq4BJIAOy2fwQ88Qmwyp/VDgjEU+QrHfYZrKN4YXxml6f66NB1nTV34zqSEvPW7Nd95CTVvqgTQvLtAnkIR0KMJMna9I/2Rj/jNG3HL9UVCw/WQAYCjHY9Khr8gw/Y4ynXqYfXCMOBSY+2sI9RJZznnUx0Oj4r4tnTahac1bzue3kESDASTn/GSmJmDfJrg2umhzGSFX1VwtFHIJMAo0aiM/8a8SYBZF3iM7KBfmWireIbKy0PllNJYLb2S4zpJ/a59pKE19p3LtFhz5p2MtGxaQIwnj4Eh4HCeLf9I0CASFhWGXuy2niCMF6iUeHEysc1wYl4HzVin3HmkqwQNEmOPq5IaubOSNdr7t8d0eEMHGb5rZNJmniivWan7X3vJqIkYiUekzLvAd96OU8ZnhnL/fRVOiEHdeVxnn3lkzyUdT1k7MsVmDLeldM25P2tDM/EscRvP3eJLs9RBv3MhdQ8PLC/qlJ9xXRu3/8hII5guKadKj8+qtas1jJ3h4BVRSQY+2avG/Fgi3OZwxikRLLTCEIC1jlZL+shX5FLnCd50BcTAH3M5Rq9ctOvTczlN1Uid+yLzaSBoHLzEMhzlJvNHZGjc/v+DwFxjH7+N/r8l/I5tp5L/t/TRDdC5oD9EBTkkBPaigYSmDXmRhn0STSSD5UMhGJDnnm5WbU5P45LvNWYciQEeiPpOqYtcZ8mxswWZHLTzqpCQ1a9cS36wRSir3TmNfr5H45riW5NjERcn0dgHO3fh0EAQpB0cvJJdJDAjFwYr17vAIl+9OdxqrJIjgIqgVRE5fe5asz5Vmk+e5d4sh2MY3+VSJJmVYGpb1Y5WNGSfMiDIXvmudsyAmCmD8CNQ6XyRdS0JkaifBNdROOgv/kuRSB55eoDUnCMe/Xdy4TPc4WMAIVI4lwJFJ250Yd81XzlrcbsWyKtaAdzIE3mVISF7b52qz/eGZuNs09wEUPWyOtHff37fwTwiZhV98pXYrfkE+W815HmaN8bgT8IWIFBeLmR1JBIrqI4memvyHFEVDNydN0R6WoHRJkbVRZr5ldQCTC+bue5s71n2X7eHgHj6pSKuYlue78cYoVZBTb6XuLrRSbHEVEBlEFckaNAjohHvRVpaX9+JRrZ7lreIc+KQB3v+8sggP+o5k71RRPdy/jn7lcZVWBsTELL39QISOblNiIq5ByLlVcmJ0kr65XoMklKnlVyaHsm46ybvfGKW5Folu3nbRAgDsAfP+SYWFrxeRQuzejxV4eABOI3kxxkI0KL5Mgc50lUPnsHWF5/mSfxMIb+2BiPpOX8bCdzICiTQ5JSnnFtUX8cs887r8asy6V9jvV9WwQ4rMAdf818NLKiiW6ETPf/RcDvWJCQAWf1RtBFQvs76c9fOSEp53CnSWY8Q54SEGN+S6OPcQI7f4vxjweQTRzHFsa4mMs8E0Mb6Y9ENbMl7iX+Rq+kH/v793YIEA+XHC5NdNv55lCaeR2ENEjweKJCeJFs4qYlGoI0/hWSORCMxBnnoNu18jzlJEDGJVzH1I1N8b/KqBP719qizr7fPwJNdPfvw95BI9AILCDQRLcAUA83Ao3A/SPwH3BxWAebF97KAAAAAElFTkSuQmCC)
Eq. (1)
Eq. (2)
The assumption of biomass optimality was maintained during simulation. Nongrowth-associated maintenance ATP was calibrated to account for 1.4 mmol/gDW/hr.
As seen from Fig. 3, while the glucose level dropped to 23 mmol during the first 100 h of batch cultivation, the propionic acid and acetate concentrations increased up to 5.8 and 2.4 mmol, respectively. The ratio of propionic acid to acetate is hence calculated to be 2.4, which compares well with the experimentally reported value of 2.2 mmol46. The growth curve also compares satisfactorily with the corresponding experimental data46. The validation of our model thus implies that iMR558 can accurately capture the key metabolic behaviors of P. freudenreichii under batch cultivation.
Flux variability analysis: the flexibility of flux distribution over the metabolic network
FBA calculates an optimal solution for the growth rate and flux distribution within the feasible flux space. However, due to the underdetermined nature of linear programming, the calculated optimal flux distribution is not unique, and thus, there are infinitive mathematical solutions for the FBA problem12,18,90. Nevertheless, the entire infinitive solutions are constrained between the boundaries of flux space, which is represented by the intervals according to which fluxes through all reactions can vary without negatively affecting the optimal growth. Although the flux distribution obtained from FBA is valuable to predict the metabolic phenotypes of an organism, more in-depth insight into the capacities of the metabolic network can be gained if the span of flux through each reaction that can support optimal growth is explored. In the COBRA methodology, this calculation is called flux variability analysis18 (FVA; see Methods section). The results of FVA hence can be interpreted in terms of the flexibility of the metabolic network to maintain growth optimality by adapting the internal fluxes. The result of FVA can also delineate the extent of uncertainty to the predictions of a typical constraint-based model that does not include regulatory and kinetic constraints12,18,90.
Overall, 262 (37%) reactions had nonzero lower/upper bounds, indicating that a large number of reactions can carry flux during strain anaerobic growth on minimal glucose-limited media. In addition, 247 (35%) reactions were found to have the same lower and upper found values, indicating nonflexibility of flux through them to satisfy the optimal anaerobic growth. The highest number of flexible reactions belonged to the pyruvate metabolism pathway.
A number of flux intervals (178) were also found to include zero. The corresponding reactions hence may be inactive when the strain grows anaerobically. The highest number of reactions with zero-including intervals was related to the vitamin B12 biosynthesis pathway. A total of 305 (43%) reactions were found to have zero lower and upper fluxes, implying their inactivity during anaerobic growth on glucose. These reactions are mostly involved in catabolism of alternative carbon sources, transport processes and fatty acid metabolism. Therefore, while they may not be required for growth on glucose-limited minima medium, they may be potentially essential for growth under alternative conditions. However, many of them stand for blocked reactions that cannot carry flux under any conditions tested (the statistics of dead-end metabolites and their reaction should be reported.)
Suboptimality analysis of P. freudenreichii metabolism
There is evidence suggesting that bacterial metabolism tends to operate virtually optimally when cultivated on minimal medium37. The optimality of bacterial metabolism in limited cultures is largely captured by the FBA37 of well-curated GEMs. FBA is a linear-programming-based algorithm predicting the optimal value of a metabolic objective function (OF) and a single flux distribution (among infinitive) supporting the optimal OF. Given the more restricted nature of and energy resource limitation in anaerobic metabolism, anaerobic strains have been shown to be more adherent to optimal utilization of resources to maximize their access to energy. Our FBA of iMR558 with growth rate as an objective function showed the consistency of P. freudenreichii metabolic function with optimality theory (see above). However, bacterial strains tend to also operate suboptimally, particularly when accessing excess carbon/energy resources. The suboptimal function of metabolism leaves a large room for the flexibility of carbon flux distribution over the metabolic network, which in turn can result in the activity of otherwise inactive pathways and thereby alternative metabolic phenotypes absent from optimal growth conditions. The analysis of suboptimality is hence essential to explore the full metabolic capacity of a strain. This is of particular importance in the context of secondary metabolism, as the secretion of most secondary metabolites is triggered only under suboptimal conditions, widening the boundaries of secondary metabolite production. In this study, to gain insight into the suboptimal behavior of P. freudenreichii metabolism, we constrained the growth rate to 95% of its FBA-predicted maximum and then performed FVA for all metabolic reactions. Interestingly, no change in the activity of silent reactions was observed, indicating the robust adaptation of the metabolic structure to optimal function. Nonetheless, significant metabolic flexibility, represented by an expanded allowable flux interval in many reactions, was observed. Of particular, acetate and propionate secretion, once exhibiting a rigid flux under optimal growth conditions, gained a significant variation interval. Suboptimality also allowed a positive succinate secretion flux, justifying the detection of this compound in some P. freudenreichii cultues52,55,71. In addition, the upper limit of the vitamin B12 production rate increased under suboptimal growth, providing an additional explanation for the broad range of reported experimental values. Our inspection hence demonstrates the utility of suboptimality analysis of metabolic networks in exploring a full range of metabolic phenotypes that may or may not be reflected in experimental literature and understanding the systems mechanisms behind alternative metabolic behaviors.
Functional analysis of P. freudenreichii central metabolism
In the previous sections, we characterized the capability of iMR558 in predicting the global metabolic features of P. freudenreichii, including growth and byproduct secretion. At the same time, knowledge of flux distribution over the key metabolic pathways is highly informative for deciphering the reported complexities associated with strain behaviors as well as successful strain engineering. Therefore, herein, we zoomed into the metabolic network at the pathway and reaction levels to explore the distribution behavior of flux across the central metabolism. In doing so, in addition to the FBA, we also use FVA to gain insight into the potential flexibility of flux behavior. However, while a reaction may take alternative flux values within an interval determined by FVA, the probability of the occurrence of arbitrary fluxes is not constant within this range. Therefore, there is a need to describe the occurrence of random variable flux values in the context of probability theory and obtain the probability distribution of flux, which is rigorously controlled by the topology of the solution space. The probability distribution of the fluxes within the network is obtained by uniform random sampling (RS)17,19, which enables us to determine the theoretically most likely state of the flux distribution. Comparison with the experimental data when available will reflect the validity of the model as a reprehensive of P. freudenreichii metabolism.
Figure 4 illustrates the distribution of carbon flux over the P. freudenreichii central metabolic network. For each reaction in this figure, the point flux value predicted by FBA, the allowable flux interval under optimal growth conditions predicted by FVA and the most likely flux value predicted by RS are reported.
Metabolic state under anaerobic growth
Experimental studies have yielded a consistent picture of the P. freudenreichii central metabolic activity under anaerobic growth conditions, that is, the strong activity of EM, WW, and acetic acid biosynthesis pathways to support cell bioenergetics requirements along with very limited activity of TCA and PP pathways. This figure is well represented by our COBRA-based flux distribution simulation, as illustrated in Fig. 4. Our GEM, however, can be used to gain more detailed insight into the anaerobic metabolic function of strains, including alternative flux distribution patterns and the flexibility of metabolic behavior. Our simulation shows that all glycolytic reactions carry a rigid flux under optimal growth except for phosphofructokinase (PFK) and enolase (ENO). PFK catalyzes the conversion of glucose-6-phosphate to fructose-diphosphate. The flexibility of flux through PFK stems from the presence of both ATP-dependent and PPI-dependent PFKs, resulting in a potential flux trade-off between these two reactions. FVA calculates a zero minimum flux for the ATP-dependent and a positive minimum flux for the PPI-dependent PFK, indicating the nonessential activity of the former and essential activity of the latter for optimal growth. This was confirmed by our knockout experiment, which indicated no impact on growth for the removal of ATP-dependent ATP, where an 8% growth rate reduction due to PPI-dependent knockout was observed. This figure suggests the dependence of the strain’s metabolism on the PPI as an energy carrier. The extent of this dependency can be quantified by the observed twofold magnitude of flux through PPI-PFK compared with ATP-PFK, as obtained from random sampling of the flux space. The preference of a PPI-dependent flux through glycolysis entails the availability of an adequate PPI pool. As seen from the flux distribution map, the majority of PPIs are compensated by two reactions downstream of glycolysis both around the phosphoenolpyruvate (PEP) node, including phosphoenolpyruvate carboxykinase (PEPCTrP) (EC: 4.1.1.38) and pyruvate-phosphate dikinase (PPDK) (EC: 2.7.9.1). Based on flux sampling, these two reactions recover over 84% of PPI consumed in glycolysis. Nonetheless, within the interval of 0.228–0.961 mmol/gDW.h, the flux through PPI-PFK has a tradeoff of as much as 0.733 mmol/gDW.h with ATP-PFK, inducing a relatively large flux variability throughout the network, especially under suboptimal conditions.
Downstream of the glycolysis pathway, pyruvate functions as a hub node in the distribution of flux across the central metabolic network. Pyruvate supplies a carbon/energy resource for the potential function of three energy generating pathways, i.e., the TCA cycle, WWC, and acetic acid biosynthesis.
The rigidity of flux through WWCs confirms the bioenergetic advantage of WWCs over potential alternative routes. We were interested in exploring the extent of this advantage by the COBRA analysis tools. To this end, we introduced specific gene knockouts leading to complete removal of WWCs. Interestingly, not only did the strain maintain its growth capability, but the electron transport chain remained functional in the absence of WWC. Furthermore, the growth was dropped no more than 10%. Hence, according to our results, the anaerobic respiration capability of the strain is independent of the WWC activity, which provides fumarate as a terminal electron acceptor. This observation is consistent with genomic and biochemical evidence indicating the diversity of electron acceptors in P. freudenreichii. Our results also suggest that Propionibacteria may have developed a robust metabolic faculty (perhaps prior to the evolution of WWCs) to secure their bioenergetic requirements and that the energetic advantage associated with WWCs may have been overestimated.
The activity of the TCA cycle in P. freudenreichii and its role in strain bioenergetics have been the focus of several studies. The mode of activity of the TCA cycle can modulate the bioenergetics of an organism and the flux distribution over the metabolic network. P. freudenreichii has the entire set of enzymes for full activity of the TCA cycle and the electron transport chain. While the activity of the TCA cycle under aerobic growth conditions is expected, the potential activity of this pathway and its contribution to energy generation under anaerobic conditions are not well understood. Our flux distribution result shows that under the anerobic growth regime, all reactions in the oxidative branch of the TCA cycle carry no or marginal flux. This observation agrees with experimental studies that found no evidence of TCA cycle activity under anaerobic conditions in glucose-limied medium46,91.
Substrate-level/Oxidative Phosphorylation
Both the genome and the structure of the P. freudenreichii metabolic network show that the strain has the entire set of enzymes for the activity of the election transport chain. This is compatible with the strain’s respiration capability under aerobic conditions. However, both experimental and our in silico analysis show that the electron transport chain significantly contributes to cell bioenergetics under an anaerobic regime. Table 5 shows the contribution of the reactions in the network to energy production in the anaerobic state. As predicted by our model, the mechanism of the strain’s anaerobic respiration involves transport of electron menaquinol to fumarate by the activity of fumarate reductase generating menaquinone. Menaquinone then accepts electrons from NADH by proton-translocating NADH dehydrogenase. The produced proton gradient generated by the latter is then consumed by F-type H+-transporting ATPase to generate ATP. While parallel experimental data are not reported for P. freudenreichii, the same mechanism has been observed in P. acidipropionici and P. acne67,92,93. In addition, P. freudenreichii may also utilize sulfate or ferrous iron as electron acceptors94, and given the denitrification capability of P. freudenreichii DSM20217, nitrate may also play the role of an electron acceptor. The activity of the electron transport chain is limited in the aerobic regime due to the diminished synthesis of cytochrome units when exposed to oxygen.
Table 5
Metabolic enzymes involved in ATP production, with their corresponding reactions, functional roles, and flux span under anaerobic conditions, minimal media with glucose
Enzyme
|
EC. No.
|
Reaction ID.
|
Function
|
Gene Association
|
Mean value (mmol/gDW/hr)
|
minFLux (mmol/gDW/hr)
|
maxFLux (mmol/gDW/hr)
|
AckA
|
2.7.2.1
|
ACKr
|
Acetate kinase
|
RM25_0078
|
-0.229 ± 0.128
|
-0.415
|
0.643
|
PGK1
|
2.7.2.3
|
PGK
|
Phosphoglycerate kinase
|
RM25_1432
|
1.888 ± 0.000
|
1.883
|
1.887
|
ATPase
|
7.1.2.2
|
ATPS4r
|
ATP Synthase
|
RM25_1128 and RM25_1130 and RM25_1127 and RM25_1129 and RM25_1131 and RM25_1132 and RM25_1133 and RM25_1134
|
0.958 ± 0.000
|
0.955
|
0.960
|
Pyk
|
2.7.1.40
|
PYK
|
Pyruvate kinase
|
RM25_1458
|
0.615 ± 0.153
|
0
|
1.426
|
AckA
|
2.7.2.1
|
PPAKr
|
Propionate kinase
|
RM25_0078
|
0.179 ± 0.128
|
0
|
1.048
|
Metabolic states during the transition from an anaerobic to an aerobic growth regime
Experimental studies have demonstrated the capacity of P. freudenreichii metabolism to support growth in the presence of oxygen. This capability is now well understood at the genetic level, as the genomic data show that the strain possesses all enzymes for full aerobic growth. The strain’s growth and metabolic flux distribution response to the level of oxygen exposure are controversial among different investigations. Our GEM proves a useful tool to shed light on the strain’s metabolic states under different physiological and environmental conditions. During our simulations, we found a relatively complex flux distribution behavior, depending on the extent of the oxygen exposure.
The onset and augmentation of oxygen exposure is accompanied by a constant decline in WWC activity and thereby propionate biosynthesis and a gradual shift from anaerobic to aerobic respiration. By further increasing oxygen exposure, a threshold is reached at which both the WWC and TCA cycles are inactive and the respiration mode is almost aerobic. An additional oxygen supply corresponds to the start of TCA cycle activity, which is constantly increased with the level of oxygen availability. Our simulation of glucose-limited growth over a wide range of oxygen supplies hence shows a switch-like relative activity of the WWC and TCA cycles during the transition from anaerobic to aerobic growth conditions, with no metabolic state in which both cycles are coactivated. In accordance with the experimental data, in silico simulation shows that P. freudenreichii metabolism can operate in a full aerobic fashion. However, due to the toxicity of oxygen for the strain, the maximum oxygen uptake and the associated growth rate observed in silico may not be realized in vivo.
A specific feature of acetic acid metabolism in propionibacteria is that it is linked to both WWCs and the TCA cycle by the function of succinyl-CoA:acetate CoA-transferase (SCACT) (EC 2.8.3.18). This reversible enzyme plays a key role in metabolic flexibility and regulates cell bioenergetics under both anaerobic and aerobic growth conditions. Figure 5 illustrates the profile of acetate production vs. OUR. Under low OUR (micro- or semiaeration regime), acetate production increases with oxygen availability. This interval is identical to the interval within which the WWC and propionate fluxes as well as the forward activity of SCACT descend from maximum to zero. The maximum acid production rate (50% of the flux through glycolysis) coincides with inactivation of both WWCs and TCAs. In this state, ATP is synthesized merely by the activity of phosphoglycerate kinase, pyruvate kinase and acetate kinase. With the increase in OUR, the acetate production pattern changes suddenly from an increasing trend to a sharply diminished trend. This change corresponds to the aforementioned switch-like transition of activity from the WWC to TCA cycle and from a forward SCACT to an inverse one. The inverse activity of SCACT leads to recycling of carbon to both oxidative and reductive branches of the TCA cycle. Given that the biosynthesis of every mole of acetate has already generated one mole of ATP, its recycling seems to realize significant energy efficiency, especially by supporting the TCA-mediated increase in flux through the electron transport chain. Confirming our results, a higher biomass and acetate yield in aerobic cultures (compared to anaerobic cultures) is accompanied by a lower yield of propionate, as previously observed in experimental studies51–53, 95. To inspect the extent of this potential efficiency, we knocked out SCACT and simulated growth at maximum OUR. Surprisingly, the growth dropped no more than 10%. This observation once again shows the bioenergetic robustness of P. freudenreichii and the importance of genome-scale analysis.
Although the decreasing trend of the acetate production rate may not be observed in experimental studies due to the inactivity of the electron transport chain at high oxygen pressures, the results of a study by Quesada-Chanto et al.96 confirms the acetate production trend predicted by our model. Their results suggest that during a gradual increase in exposure to oxygen after a certain point, a slightly lower acetate yield is obtained.
In addition to the electron transport chain, the expression of a pyruvate oxidase in cultures of P. freudenreichii has recently become a controversial route of oxygen consumption. Two research works postulated that the activity of pyruvate oxidase is responsible for consuming oxygen at low pressures54,97. However, our simulations suggest that oxygen consumption relies solely on cytochrome bd oxidase at both low and high oxygen uptakes.
Biosynthesis of other organic acids
Optimal FVA anticipated a zero variability range for the secretion rates of malate and citrate. However, secretion of trace amounts of lactate, pyruvate and succinic acid was predicted to be possible. However, propionic acid and acetate secretion rates exhibit flux rigidity illustrated by the tight flux intervals. These FVA predictions are in agreement with the results of a previous experimental report on the detection of acetate and propionate as the only major fermentation end products in minimal cultures of P. freudenreichii46.
The CO2 fixation capacity
As seen from Fig. 4, PEP has the potential to transform alternatively into pyruvate via pyruvate kinase (EC: 2.7.1.40) or into oxaloacetate via diphosphate-phosphoenolpyruvate carboxykinase (EC: 4.1.1.38). A large flux interval is predicted for both reactions, which is consistent with their competitive action for mediating phosphoenolpyruvate concentration. Although the flux through pyruvate kinase is bounded by zero, sampling of flux space shows that these reactions may be active simultaneously, although flux through the direct route is more favorable. The activity of the diphosphate-phosphoenolpyruvate carboxykinase toward the formation of oxaloacetate is accompanied by CO2 fixation. However, according to our results, CO2 fixation is not marginal in contrast with previous reports48, and the flux through this reaction is no less than 20% of the total flux through glycolysis. Furthermore, the activity of diphosphate-phosphoenolpyruvate carboxykinase is essential to achieve the maximum biomass yield.
Cobalamin coenzyme production
A major biotechnological importance of P. freudenreichii is relevant to its potential for industrial production of vitamin B12. In silico characterization of the strain capacity for cobalamin coenzyme biosynthesis can inform the development of economically viable cobalamin coenzyme bioprocesses. The COBRA methodology is efficiently used to elucidate how and under which constraints microbial metabolism can direct resources toward secondary metabolites. Hence, we used constraint-based analysis of iMR558 to characterize cobalamin coenzyme biosynthesis in P. freudenreichii compared with experimental data and to explore the strain capacity for cobalamin coenzyme production.
While the upper ligand of cobalamin coenzyme can be synthesized under both anaerobic and (micro)aerobic culture regimes, the biosynthesis of DMBI requires (mico)aerobic growth conditions63. Consistent with the experimental data, in silico simulation of iMR558 showed zero flux through DMBI synthesis unless small amounts of oxygen uptake were allowed, whereas the upper ligand was found to be produced regardless of aeration regime.
We next performed FBA, FVA and RS to evaluate the capacity of P. freudenreichii in cobalamin coenzyme production. As seen from Table 6, FBA predicts no positive flux through adenosyl cobinamide synthesis. However, a positive upper limit of flux was predicted by FVA, indicating the capability of the strain to produce adenosyl cobinamide under optimal growth conditions. The RS also predicted the most likely value of flux through adenosyl cobinamide to be 0.004 µmol/gDW/hr per GUR of 1 mmol/gDW/hr. Therefore, the yield of cobalamin coenzyme production is calculated to be 180 µg B12/gDW, falling well within the range of experimentally observed values (Table 6). Higher values are also reported in the literature. The higher production yields can be attributed to the complex medium used for cultivating P. freudenreichii containing corn-steep liquor or yeast extract, whereas we used a glucose-limited medium for our in silico experiments52,98. As a secondary metabolite, cobalamin coenzyme competes with growth for energy, and hence, the coupling of its production with growth is less probable.
Table 6
Cobalamin production on glucose during anaerobic cultivation (99% optimal: captured by FVA and sampling)
Growth (gDW/mol substrate)
|
Minimum flux (mmol/gDW/h)
|
Maximum flux (mmol/gDW/h)
|
Most probable (mmol/gDW/h)
|
Yield on Biomass (in silico) (µg B12/gDW) (Most probable)
|
Yield on Biomass (in vivo) (µg B12/gDW)
|
Reference
|
0.0392
|
0
|
8.69 E-05
|
4.23 E-06 ± 4.04 E-06
|
180 ± 172
|
127,180,201
|
98
|
We also used robustness analysis to evaluate the potential effect of proton exchange on cobalamin production. Our simulation showed an inverse relationship between proton export flux and the amount of cobalamin produced. Consistently, experimental data have shown that vitamin B12 production can be reduced in acidic environments, where the highest production level was found in neutral environments (pH = 7)53.
Gene and Reaction Essentiality Analysis of P. freudenreichii Metabolism
Rational metabolic engineering of biotechnologically important strains for improved value-added metabolite production often relies on knockout of particular genes from the biocatalyst genome. Considering the time and cost associated with gene knockout experiments, prior knowledge of gene essentiality can accelerate strain improvement efforts by helping eliminate lethal knockout scenarios. Accurate and well-curated GEMs are invaluable tools to predict gene knockouts that potentially inhibit strain growth.
We conducted a gene/reaction essentiality study of P. freudenreichii by sequentially silencing genes/reaction in iMR558 and simulating anaerobic growth (GUR = 1 mmol/gDW/hr). Gene deletion analysis identified a core set of 154 genes (28% of the total) associated with 308 reactions as essential. Additionally, the knockout of a separate set of 34 genes associated with 39 reactions was found to reduce the growth rate by 61–99% compared with the wild-type strain.
In our reaction essentiality analysis, we also identified 220 (31% of the total) reactions as essential.
Of the total essential reactions, 39 are associated with isozymes, requiring double/multiple gene knockout for their deletion, and 22 have no gene association. These data imply that nearly 80% of reactions identified as essential at the metabolic level are also essential at the genetic level. Furthermore, 15 nonlethal reaction removals were identified to reduce the growth rate below 95% of that in the wild-type strain.
Overall, the results of iMR558 essentiality analysis imply a relatively low robustness of the P. freudenreichii metabolic network against gene/reaction perturbation, with nearly one-third of all genes/reactions being absolutely essential for vital metabolic functions (comparison to other strains). The ratio of the essential genes to the essential reactions is 0.7, which is expected to increase when genes responsible for orphan essential reactions will be identified. Although comparison of in silico gene/reaction deletion results with in vivo knockout studies can be used to validate the structure of the metabolic models, the lack of corresponding data for P. freudenreichii limits our ability for such validation. However, our gene knockout simulation may provide insight into the development of more efficient metabolic engineering scenarios.
Study strengths and limitations
The metabolic network of P. freudenreichii was reconstructed under very limited availability of biochemical and physiological data of this strain. Therefore, a major part of our metabolic reconstruction relies on genome annotation data. However, extensive literature-based manual development and curation of the network and careful reannotation of many ORFs render the model sufficiently accurate to represent a wide range of strain phenotypes reported in the literature. With rigorous reannotation of P. freudenreichii ORFs, we tried to extend the strain’s metabolic network toward an ideal genome-scale network as much as possible. Nonetheless, the annotation of the strain’s whole genome is still incomplete, where almost 40% of the ORFs do not have an annotation or assigned putative function. Therefore, the strain may feature other metabolic capabilities that are not captured by our reconstruction. Future availability of more extended high-throughput and omics data will provide opportunity for further improvement of the reconstruction in terms of genome coverage, phenotype prediction, and quantitative accuracy.