Genome Scale Constraint-Based Metabolic Reconstruction of Propionibacterium Freudenreichii DSM 20271


 Propionibacterium is an anaerobic bacterium with a history of use in the production of Swiss cheese and, more recently, several industrial bioproducts. While the use of this strain for the production of organic acids and secondary metabolites has gained growing interest, the industrial application of the strain requires further improvement in the yield and productivity of the target products. Systems modeling and analysis of metabolic networks are widely leveraged to gain holistic insights into the metabolic features of biotechnologically important strains and to devise metabolic engineering and culture optimization strategies for economically viable bioprocess development. In the present study, a high-quality genome-scale metabolic model of P. freudenreichii ssp. freudenreichii strain DSM 20271 was developed based on the strain’s genome annotation and biochemical and physiological data. The model covers the functions of 23% of the strain’s ORFs and accounts for 711 metabolic reactions and 647 unique metabolites. Literature-based reconstruction of the central metabolism and rigorous refinement of annotation data for establishing gene-protein-reaction associations renders the model a curated omic-scale knowledge base of the organism. Validation of the model against experimental data indicates that the reconstruction can capture the key structural and functional features of P. freudenreichii metabolism, including the growth rate, the pattern of flux distribution, the strain’s aerotolerance behavior, and the change in the mode of metabolic activity during the transition from an anaerobic to an aerobic growth regime. The model also includes an accurately curated pathway of cobalamin biosynthesis, which was used to examine the capacity of the strain to produce vitamin B12 precursors. Constraint-based reconstruction and analysis of the P. freudenreichii metabolic network also provided novel insights into the complexity and robustness of P. freudenreichii energy metabolism. The developed reconstruction, hence, may be used as a platform for the development of P. freudenreichii-based microbial cell factories and bioprocesses.


Introduction
Propionibacterium freudenreichii is a gram-positive microaerophilic bacterium with a long history of being used in the manufacturing of Swiss cheese 1 . Being studied as a model organism to understand the physiology of Propionibacteria, the strain has recently received renewed attention for its potential biotechnology application. Propionibacterium is capable of producing several industrial bioproducts, including cobalamin, folic acid, trehalose and propionic acid [2][3][4][5][6] . Recent studies have demonstrated the health-supporting use of P. freudenreichii in edible products, relevant to anticancer and antimutagenic effects 7,8 . While this bacterium is capable of growing in both fermentation and anaerobic respiration modes, its aerotolerance depends on the level of oxygen exposure. P. freudenreichii is classically divided into two subspecies based on two features: 1) fermenting lactose and 2) reducing nitrates. Data in previous studies shows that while ssp. shermanii can metabolize lactose as a carbon source but not reduce nitrate as a nitrogen source, and the inverse is apparently true of ssp. freudenreichii 9 . However, recent studies have shown that some strains of these subspecies have both or neither of these capabilities 10 , and thus, the abovementioned classi cation might be inaccurate.
The central metabolism of P. freudenreichii has been explored by several research groups. Interest in the study of Propionibacteria metabolism may partially be relevant to their metabolic peculiarities, particularly the presence of the Wood-Werkman Cycle, which quality constraint-based models have demonstrated the ability to replicate the metabolic behaviors of the relevant organisms and generate testable hypotheses to decipher complexities associated with their physiology. Systems analysis of metabolic ux distribution and bioenergetic pro ling of industrially relevant organisms can facilitate the generation of novel metabolic engineering strategies to enhance the yield and productivity of bioprocesses.
In this study, a constraint-based metabolic model of Propionibacterium freudenreichii DSM20271 (iMR558) was reconstructed based on the strain's published genome sequence 15 and physiological data and according to the high-quality metabolic reconstruction criteria 16 . Herein, rst, the detailed steps for in silico reconstruction and analysis of P. freudenreichii metabolism, including crurtion and re nement of the central and secondary metabolic pathways, are presented. A set of constraint-based analysis methods 17 , including ux balance analysis (FBA) 11,16 , ux variability analysis (FVA) 18 and Monte Carlo ux sampling 19, are used for systems-level characterization of P. freudenreichii metabolism. The behavior of the in silico strain is validated against the experimental data. Constraint-based analysis tools are then applied to estimate the metabolic capacities of the strain, including cobalamin biosynthesis, to gauge the exibility of metabolic ux distribution under optimal and suboptimal growth conditions and to characterize the metabolic states of the strains during the transition from an anaerobic to an aerobic growth regime.
While addressing several inconsistencies reported from various experimental studies of P. freudenreichii, our model re ects the overall experimentally observed behaviors of the strain under various conditions, the usefulness of COBRA in exploring the capacities of the strain, and the need for much further high-throughput multiomic studies to understand the complexities associated with the strain's metabolism.

Metabolic Network Reconstruction
To the extent that experimental and bioinformatics data permit, our metabolic network reconstruction complies with the protocol documented by Thiele and Palsson for reconstruction of high-quality genome-scale metabolic netwroks 16 .
Reconstruction of the metabolic network of Propionibacterium freudenreichii ssp. freudenreichii DSM 20271 was initiated by retrieving annotation data from KEGG [20][21][22] and PATRIC RASTtk-enabled Genome Annotation Service 23 . Based on annotation data, a catalog of metabolic reactions associated with their encoding genes and catalyzing enzymes was developed, resulting in a draft metabolic network. The draft network was then used as a template to develop a re ned and functional metabolic reconstruction.
Due to the inaccuracy and/or insu ciency of some predictions by automated bioinformatic tools, the draft reconstruction assembled upon annotation data often involves several de ciencies, impeding it to function as an in silico operative representation of the metabolism. Some of these de ciencies include discrepancies in functional assignment to the genes, ambiguities in gene regulatory relationships, missing links in crucial pathways, undetermined cofactor speci cities of the enzymes and unclear directionality of the reactions. Identifying and resolving these de ciencies requires manual curation of the draft reconstruction based on available data in the literature and biochemical resources as well as assumptions of relevant hypotheses. To this end, comprehensive knowledge bases including Metacyc 24 , BRENDA 25 , UNIPROT 26 and eQuilibliator 27 together with literature were extensively consulted to address the missing link, determine reaction directionalities, and remove inconsistencies in the network. Our network curation process involved nding and lling the gaps in the network, correcting the functional assignments of misannotated genes, re nement of gene-protein-reaction (GPR) associations, determining the enzymes' substrate and cofactor speci cities, determining the likely direction of reactions, and eliminating futile cycles and arti cial reactive loops. Gap lling and network re nement were accompanied by establishing the corresponding GPR association whenever possible. The process of genetically evident gap lling was carried out by cross-species homology analysis using a BLAST search (https://blast.ncbi.nlm.nih.gov/Blast.cgi) 28 . Some metabolic reactions were also included in the reconstruction based on literature evidence, although the corresponding coding genes were not identi ed.
Because the uptake and secretion system of P. freudenreichii is poorly characterized, there was a need to explore relevant resources, including the TCDB database 29,30 , ABCdb repertory 31 , and UNIPROT 26 , to nd potential transporters of the strain. For many transported metabolites, these resources did not introduce responsible genes. Therefore, transport mechanisms were mostly identi ed by BLAST-based homology studies in TCDB. The uptake/secretion mechanisms of experimentally evident transported metabolites that were not detected by homology study were assumed to be similar to those in B. subtilis 32 .
Typical genome-scale metabolic network reconstructions (GENREs) include a biomass equation representing virtually all constituents of the cell and are used as an objective function for the simulation of growth. A comprehensive biomass composition of P. freudenreichii has not been published. However, previous analyses have shown that the growth and ux distribution phenotypes are not signi cantly sensitive to biomass composition 33 . Therefore, the biomass equation of P. freudenreichii was formulated by modifying the experimentally determined biomass composition of the gram-positive bacterium B. subtilis 32 . Growthassociated ATP maintenance (GAM) in the biomass equation was assumed to be 35 based on the P. freudenreichii theoretical calculated ATP yield (28.8) by Papoutsakis et al 34 .
A reaction was also added to the metabolic network to account for nongrowth-associated ATP maintenance (NGAM). The initial value of ux through the ATP maintenance (ATPM) reaction was assumed to be 0.76 based on a previously determined ATPM for P. acidipropionici 35 . The ATPM ux, however, was adjusted to calibrate the model by replicating the experimentally observed growth yield.
All of the reactions included in the metabolic network were mass balanced by correcting stoichiometric coe cients and were charge-balanced by adding proton or water molecules.

Computational modeling of the metabolic network
The developed metabolic network was converted to a computational constraint-based model (CBM) according to the COBRA methodology. The topology of the metabolic network was mathematically represented as a stoichiometric matrix of S m,n, where m denotes the number of metabolites aligned in the rows and n denotes the number of reactions aligned in the columns. In the stoichiometric matrix, each element of S ij is the coe cient of the ith metabolite in the jth reaction. When data are available, each reaction is associated with its catalyzing enzyme, which in turn is linked to its encoding gene(s). The relationship between genes, proteins, and reactions is mathematically represented by Boolean logical regulatory rules.
For each extracellular metabolite, an exchange reaction was included in the computational model to facilitate de ning the environment of the in silico strain. The ux space of the metabolic reactions was constrained by limiting to zero the lower bounds of all irreversible reactions and all exchange reactions but those relevant to the medium constituents. An additional constraint was imposed by setting the lower bound of the exchange ux of the limiting substrate (glucose) to the experimental value. An SBML version of the computational model is given in supplementary 1.

Model Curation and Re nement
Network gap lling and model re nement were carried out by an iterative approach. The initial model was rst examined for its capability in simulating the growth and byproduct secretion phenotypes observed from P. freudenreichii experimental cultures. As mentioned above, the initial draft model was not able to replicate the metabolic phenotypes due to the gaps in biosynthesis of the precursors of the biomass and byproducts. Therefore, an iterative process of gap nding-gap lling-simulation was run until the network was rendered functionally operative. Further model prediction improvement was also carried out by inspecting and re ning the metabolic network as described earlier. We also revised the reactions of the model for cofactor consumption, as it has been reported that this feature affects the metabolic behavior of the genome-scale model 36 .

Model Simulation and Analysis
The model was simulated and analyzed for the purposes of 1) curating the metabolic network, 2) validating the reconstruction, and 3) exploring the metabolic capabilities of P. freudenreichii. Model analysis was carried out by using a series of constraint-based methods, including ux balance analysis (FBA), ux variability analysis (FVA), random sampling (RS) and robustness analysis (RA), as described below.

Flux Balance Analysis:
In metabolism, most metabolites are involved in more than one reaction. This implies that the stoichiometric constraint governing the function of the metabolic network can be satis ed by an in nite number of ux distribution states. FBA is a linear-programming optimization technique that calculates one steady-state ux solution corresponding to the optimal value of a special ux considered the metabolic objective function (OF). Given that bacterial strains are evidently assumed to evolve to maximize growth, biomass is often considered the OF of metabolism, enabling FBA to simulate strain growth and the associated metabolic ux distribution in silico. Nevertheless, any reaction other than the biomass equation can also be the subject of optimization by FBA, allowing estimation of the metabolic network capacity for its biosynthesis.
A typical FBA problem can be formulated as follows: maxc T ν, Subject to S ij .ν j = 0, Flux Variability Analysis: FBA can only calculate one among in nitive ux states corresponding to the optimal values of an OF. However, these in nitive ux states are all contained in a close solution space. This implies that each reaction ux can only vary within a limited interval to avoid violating the optimality of the OF. FVA is an FBA-based technique to calculate the interval within which each reaction ux can change while still supporting the FBA-determined OF value. For an FBA model with biomass as the objective function, the FVA problem is formulated as follows: max/min ν j , Subject to Sij .ν j = 0, c T ν = Z biomass While the assumption of optimal bacterial growth on substrate-limited bacteria has experimental support 37 , the strains are also found to grow suboptimally, especially under different batch and/or rich cultures. FVA can be used to quantify the suboptimal activity of metabolism by constraining the upper bound of the biomass production rate to values below the maximum 18 . In our study, we performed all FVA-based suboptimality analyses by restricting the growth rate to 95% of its FBA-determined maximum value.

Sampling of Solution Space
Although all reaction ux values within the FVA-determined interval will yield an identical OF value, the probability distribution of ux within the feasible range is not uniform. Random sampling of the ux space can elucidate the structure of the ux distribution, which in turn can be analyzed to obtain the most likely ux values within the allowable space 19 . Whenever informative, we sampled the ux space at 1500 random points (more than 2 × number of extended model reactions) using the Monte Carlo algorithm 19,38 . The ux distribution of the reactions of interest was then inspected for skewness, and the mean, standard deviation (SD), and median of each ux was calculated from the normalized ux distribution.

Computational environment and tools
All simulations were carried out in the MATLAB environment by COBRA Toolbox 17 . The linear programming solver 'CPLEX' was used as the optimization solver.

Gene and Reaction Essentiality Analysis
The essential gene/reaction for growth was explored by sequentially removing each gene/reaction and performing FBA.
Gene/reaction knockouts leading to a zero growth rate were considered to be essential, and those that only reduced growth were considered to be growth supporting. Knockouts with no effect on growth were considered to be nonessential.

Results And Discussion
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 re nement. 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 ribo avin. The computational version of the metabolic network also includes 95 exchange reactions (one per extracellular metabolite) to de ne the boundaries of the network with its in silico environment.  Dead-end metabolites 225 iMR558 contains 184 reactions that were added during the gap-lling 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-lling 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 nal 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-ller reactions, 16 were transporters, two were required for ubiquinone biosynthesis, ve 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 identi ed speci cations 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 data [39][40][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:0 41 . 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 composition [39][40][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. freudenreichii 42 . 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 speci c to P. freudenreichii.
Lipoteichoic acid was excluded from the biomass equation because of a lack of genetic and biochemical evidence in P.

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 identi ed from the DSM20271 genome. Based on homology modeling, Loux et al explored the uptake mechanism of a number of carbohydrates in several Propionibacterium subspecies 10 . While their study did not include strain DSM20271, it revealed that the subspecies of the P.
freudenreichii genera have diverse stain-speci c sugar utilization capabilities 10 . To identify DSM20271-speci c 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 nd 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 nal 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 lled due to the limited knowledge of their degradation pathways and unfruitful cross-species sequence similarity search.
Curation and re nement 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 characterized 44,45 . In line with experimental evidence 46-48, the genomic data con rm 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 cycle 45,48−50, which may justify the observed oxidative phenotypes with free oxygen and the ability of the strain to grow under aerobic conditions [51][52][53][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), membraneanchored 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 ATPase 47,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 ow through WWC were identi ed 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 glucose 55 . 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 media 46,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 (ribo avin), 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 re ected 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 studies 58-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 speci city 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 pathway 61 . The lower ligand is N1-(alpha-D-ribosyl)-5,6-dimethylbenzimidazole, which is biosynthesized from reduced avin 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 assembly 59 . 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 identi ed 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 ux during the model simulation. Our gap analysis revealed two structural de ciencies leading to the blockage of ux through this pathway: 1) inability of the network to metabolize or exert two byproducts of the enzymatic reactions, including Sadenosylhomocysteine 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 de ciencies 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. freudenreichii 68 . Our similarity search identi ed 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 identi ed 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 lled the network gap accordingly.
There is biochemical evidence that (R)-1-aminopropan-2-ol cannot be metabolized in P. freudenreichii DSM20271 69, which is consistent with our unfruitful BLAST search. On the other hand, one of the two enzymatic steps for the synthesis of (R)-1aminopropan-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 Simulation of carbon source utilization was carried out using ux balance analysis (FBA). Table 3 compares the carbon source utilization pro le 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. 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 culture 35 . 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 glucose 6,71 . P. freudenreichii was also able to grow in silico in the presence of adenosine as described by Glatz and Anderson 56 and in a medium containing nicotinic acid, cysteine and tryptophan de ned by Crow 72 .
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 sources 47,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 identi ed cysteine as the only amino acid-based source of sulfur for P. freudenreichii.

Growth and byproduct secretion
Well-curated strain-speci c 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 ux variability analysis (FVA) were used to simulate microbial growth and estimate ux distribution over the metabolic network. To enable validation of the simulation results, the in silico culture conditions matched those reported from the experimental growth cultures 74 Table 4. 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 uxes 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. 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 ux was calculated by equations (1) and (2) 89 and speci cally calibrated to enter the system with the corresponding value of -1 mmol/gDW/hr. rglucose, S, µ, X and Δt correspond to the glucose uptake ux, glucose molar concentration, speci c growth rate, biomass concentration and fermentation time interval, respectively.
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. Flux variability analysis: the exibility of ux distribution over the metabolic network FBA calculates an optimal solution for the growth rate and ux distribution within the feasible ux space. However, due to the underdetermined nature of linear programming, the calculated optimal ux distribution is not unique, and thus, there are in nitive mathematical solutions for the FBA problem 12,18,90 . Nevertheless, the entire in nitive solutions are constrained between the boundaries of ux space, which is represented by the intervals according to which uxes through all reactions can vary without negatively affecting the optimal growth. Although the ux 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 ux through each reaction that can support optimal growth is explored. In the COBRA methodology, this calculation is called ux variability analysis 18 (FVA; see Methods section). The results of FVA hence can be interpreted in terms of the exibility of the metabolic network to maintain growth optimality by adapting the internal uxes. 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 constraints 12,18,90 .
Overall, 262 (37%) reactions had nonzero lower/upper bounds, indicating that a large number of reactions can carry ux 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 non exibility of ux through them to satisfy the optimal anaerobic growth. The highest number of exible reactions belonged to the pyruvate metabolism pathway.
A number of ux 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 uxes, 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 ux 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 medium 37 .
The optimality of bacterial metabolism in limited cultures is largely captured by the FBA 37 of well-curated GEMs. FBA is a linearprogramming-based algorithm predicting the optimal value of a metabolic objective function (OF) and a single ux distribution (among in nitive) 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 exibility of carbon ux 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, signi cant metabolic exibility, represented by an expanded allowable ux interval in many reactions, was observed.
Of particular, acetate and propionate secretion, once exhibiting a rigid ux under optimal growth conditions, gained a signi cant variation interval. Suboptimality also allowed a positive succinate secretion ux, justifying the detection of this compound in some P. freudenreichii cultues 52,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 re ected 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 ux 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 ux across the central metabolism. In doing so, in addition to the FBA, we also use FVA to gain insight into the potential exibility of ux behavior. However, while a reaction may take alternative ux values within an interval determined by FVA, the probability of the occurrence of arbitrary uxes is not constant within this range. Therefore, there is a need to describe the occurrence of random variable ux values in the context of probability theory and obtain the probability distribution of ux, which is rigorously controlled by the topology of the solution space. The probability distribution of the uxes within the network is obtained by uniform random sampling (RS) 17,19, which enables us to determine the theoretically most likely state of the ux distribution.
Comparison with the experimental data when available will re ect the validity of the model as a reprehensive of P. freudenreichii metabolism. Figure 4 illustrates the distribution of carbon ux over the P. freudenreichii central metabolic network. For each reaction in this gure, the point ux value predicted by FBA, the allowable ux interval under optimal growth conditions predicted by FVA and the most likely ux 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 gure is well represented by our COBRA-based ux 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 ux distribution patterns and the exibility of metabolic behavior. Our simulation shows that all glycolytic reactions carry a rigid ux under optimal growth except for phosphofructokinase (PFK) and enolase (ENO). PFK catalyzes the conversion of glucose-6-phosphate to fructose-diphosphate. The exibility of ux through PFK stems from the presence of both ATP-dependent and PPI-dependent PFKs, resulting in a potential ux trade-off between these two reactions. FVA calculates a zero minimum ux for the ATP-dependent and a positive minimum ux for the PPI-dependent PFK, indicating the nonessential activity of the former and essential activity of the latter for optimal growth. This was con rmed 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 gure suggests the dependence of the strain's metabolism on the PPI as an energy The rigidity of ux through WWCs con rms 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 speci c 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 ux 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 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 signi cantly 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 Ftype 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. acne 67,92,93 . In addition, P. freudenreichii may also utilize sulfate or ferrous iron as electron acceptors 94, and given the denitri cation 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. 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 speci c 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 exibility and regulates cell bioenergetics under both anaerobic and aerobic growth conditions. Figure 5 illustrates the pro le 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 uxes as well as the forward activity of SCACT descend from maximum to zero. The maximum acid production rate (50% of the ux 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 signi cant energy e ciency, especially by supporting the TCA-mediated increase in ux through the electron transport chain. Con rming 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 studies [51][52][53]95 . To inspect the extent of this potential e ciency, 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 con rms 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 pressures 54,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 ux rigidity illustrated by the tight ux 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.
The CO2 xation capacity As seen from Fig. 4 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 e ciently 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 conditions 63 . Consistent with the experimental data, in silico simulation of iMR558 showed zero ux 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 ux through adenosyl cobinamide synthesis. However, a positive upper limit of ux 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 ux 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 experiments 52,98 . As a secondary metabolite, cobalamin coenzyme competes with growth for energy, and hence, the coupling of its production with growth is less probable. 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 ux 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 identi ed 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 identi ed 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 identi ed as essential at the metabolic level are also essential at the genetic level. Furthermore, 15 nonlethal reaction removals were identi ed 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 identi ed. 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 e cient 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 su ciently 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.

Conclusions
Our work presents a high-quality full functional metabolic reconstruction and constraint-based model of P. freudenreichii ssp. freudenreichii DSM 20271. A comprehensive literature-based curation strategy was adopted for network development and re nement, leading to a computational model capable of replicating several behavioral characteristics of P. freudenreichii metabolism, long being discussed in the literature. Our model predicted the strain's growth and distribution of carbon ux over the metabolic network, consistent with the relevant available common knowledge. The shift in the metabolic behavior of P. freudenreichii during the transition from an anaerobic to an aerobic growth regime was not only predicted in accord with experimental observations but also new insight into the response of the strain to the level of oxygen exposure and deep robustness of the energy generation from the electron transport chain irrespective of the activity of WWC, and some other pathways were revealed in silico. The reconstruction also includes a well-curated vitamin B12 biosynthesis pathway compatible with FBA requirements, which can be used in other metabolic reconstructions to screen and engineer high-yielding hosts for biosynthesis.
The developed reconstruction hence can serve as a platform to gain a more in-depth understanding of the physiology of P. freudenreichii and to develop the strains to an industrial producer of value-added bioproducts.