- Open Access
Integration of enzymatic data in Bacillus subtilis genome-scale metabolic model improves phenotype predictions and enables in silico design of poly-γ-glutamic acid production strains
© The Author(s) 2019
- Received: 10 October 2018
- Accepted: 29 December 2018
- Published: 9 January 2019
Genome-scale metabolic models (GEMs) allow predicting metabolic phenotypes from limited data on uptake and secretion fluxes by defining the space of all the feasible solutions and excluding physio-chemically and biologically unfeasible behaviors. The integration of additional biological information in genome-scale models, e.g., transcriptomic or proteomic profiles, has the potential to improve phenotype prediction accuracy. This is particularly important for metabolic engineering applications where more accurate model predictions can translate to more reliable model-based strain design.
Here we present a GEM with Enzymatic Constraints using Kinetic and Omics data (GECKO) model of Bacillus subtilis, which uses publicly available proteomic data and enzyme kinetic parameters for central carbon (CC) metabolic reactions to constrain the flux solution space. This model allows more accurate prediction of the flux distribution and growth rate of wild-type and single-gene/operon deletion strains compared to a standard genome-scale metabolic model. The flux prediction error decreased by 43% and 36% for wild-type and mutants respectively. The model additionally increased the number of correctly predicted essential genes in CC pathways by 2.5-fold and significantly decreased flux variability in more than 80% of the reactions with variable flux. Finally, the model was used to find new gene deletion targets to optimize the flux toward the biosynthesis of poly-γ-glutamic acid (γ-PGA) polymer in engineered B. subtilis. We implemented the single-reaction deletion targets identified by the model experimentally and showed that the new strains have a twofold higher γ-PGA concentration and production rate compared to the ancestral strain.
This work confirms that integration of enzyme constraints is a powerful tool to improve existing genome-scale models, and demonstrates the successful use of enzyme-constrained models in B. subtilis metabolic engineering. We expect that the new model can be used to guide future metabolic engineering efforts in the important industrial production host B. subtilis.
- Genome-scale metabolic model
- Enzymatic data
- Bacillus subtilis
- Constraint-based methods
- Poly-γ-glutamic acid
Bacillus subtilis is the model organism for Gram-positive bacteria and one of the best-characterized bacteria. Its genome has been sequenced and the microbiological, molecular and genetic methodologies for its cultivation and manipulation are well established . Recent studies have generated comprehensive absolute proteome  and transcriptome datasets under multiple different growth conditions [3–5]. B. subtilis grows efficiently with low-cost carbon and nitrogen sources, it is Generally Recognized as Safe (GRAS) by the US Food and Drug Administration (FDA), and is an attractive chassis for synthetic biology applications thanks to its natural competence, easy chromosomal integration, and capability to efficiently secrete products such as enzymes and vitamins even without engineering [6, 7]. Moreover, B. subtilis has important applications in agriculture as plant growth promoting bacterium (PGPB) [8, 9]. Finally, this bacterium has the ability to produce a small, resistant and metabolically dormant spore, which allows easy storage and transportation of B. subtilis strains.
For further improvement of productivity of endogenous or exogenous target metabolites in B. subtilis, genetic modifications to redirect flux and fermentative process optimization are essential. A wide range of studies have been carried out aimed at optimizing B. subtilis as a cell factory for industrially relevant bioproducts such as riboflavin [10, 11], subtilisin [12, 13], synthetic xylanases , and poly-γ-glutamic acid (γ-PGA) [15, 16]. These studies have demonstrated industrially relevant improvements in production, but they have also shown that the rational identification of target mutations to optimize the metabolism of B. subtilis still represents a major challenge necessitating the use of trial-and-error approaches such as random mutagenesis and screening. In order to complement these labor and cost intensive approaches with rational design tools, multiple different versions of B. subtilis genome-scale metabolic models (GEMs) have been reconstructed during the past 10 years [17–20]. Some of these GEMs have been successfully used to drive metabolic engineering applications, such as the production of 3-hydroxypropanoic acid , riboflavin, cellulases (R,R)-2, -3-butanediol and isobutanol .
Although the available B. subtilis GEMs have generally been validated against experimentally determined growth rate data on different carbon substrates as well as by systematic gene essentiality analyses, all the existing models have a relatively low prediction accuracy for central carbon metabolic fluxes and key byproduct secretion rates. Since those models are already fairly complete in terms of representation of the majority of metabolic reactions that can take place in B. subtilis, it is likely that prediction accuracy could only be improved through the integration of additional information on bacterial metabolism that is not captured by the current GEMs.
Since the GEMs are generally analyzed through constraint-based methods, the predicted growth rate and the production of the target metabolite are mainly limited by the carbon source uptake rate, constrained based on experimental measurements. However, each metabolic flux is highly dependent on the concentration and kinetics of the enzymes catalyzing the reaction. For this reason, the predictions based only on uptake flux constraints may not agree with the experimental behavior. Different approaches for the integration of enzyme concentrations in GEMs have been developed to constrain the solution space and improve phenotypic predictions. FBAwMC  was one of the first approaches, based on the flux balance analysis (FBA) method, imposing concentration constraints for enzymes within the crowded cytoplasm to improve the prediction of growth rates of Escherichia coli under different growth media, without using the measurements of nutrient uptake rates. Other methods were proposed as an extension of FBAwMC, such as the metabolic modeling with enzyme kinetics (MOMENT) , which utilizes the kinetic parameters under the limitations of the total enzyme pool available. Similarly, Nilsson et al.  used an extension of FBA to predict the metabolic trade-offs in yeast, in which the sum of fluxes was constrained to the sum of the product of the maximum in vitro activity and the total enzyme mass. An alternative approach integrates quantitative measurements of protein and metabolite levels into GEM, by associating them with metabolic fluxes by using Michaelis Menten-like rate equations . The most recent and promising method, referred to as GECKO (Enzymatic Constraints using Kinetic and Omics data), uses enzymatic data, in the form of protein abundance and turnover number, as a new constraint for each metabolic flux, ensuring that fluxes do not exceed the maximum capacity in a given condition . GECKO was applied to Saccharomyces cerevisiae GEM and it was shown to allow more accurate predictions than the corresponding non-enzyme constrained GEM for growth rates under different carbon sources, for the simulation of the critical dilution rate for Crabtree effect, and for the metabolic behavior of a single-gene knockout strain. The GECKO approach is also relatively simple to implement computationally compared to methods that use more complex rate equations.
In this work, we integrate a set of available enzyme constraints (absolute protein levels and turnover numbers) for the reactions of central carbon metabolism into the iYO844 GEM of B. subtilis following the principles of GECKO. The aim was to improve the prediction accuracy of central carbon flux distributions and secretion rates of the main metabolic byproducts with the idea that the resulting model would also allow improved strain design. We use experimental secretion and intracellular flux data of wild-type and mutant strains to evaluate the accuracy of the enzyme-constrained model. We finally demonstrate the metabolic engineering potential of the model by designing knockout strategies for improving the production of γ-PGA, a biopolymer with promising features and applications as food, cosmetics and pharmaceutical additive. We show that the GECKO model predicts different optimal knockout strategies than the standard GEM, and that in vivo implementation of the suggested single-reaction knockouts results in a twofold improvement in γ-PGA production (both rate and titer) from glucose.
List of kcat values and protein quantifications integrated in iYO844 model
Organism of kcat data
Refs. for kcat data
g6p → f6p
1.55 × 10−5
dhap → g3p
1.28 × 10−5
g3p + nad +pi → 13dpg + h + nadh
5.77 × 10−5
13dpg + adp → 3pg + atp
3.60 × 10−5
3pg → 2pg
8.85 × 10−6
2pg → h2o + pep
3.17 × 10−5
g6p + nadp → 6pgl + h + nadph
8.05 × 10−6
accoa + h2o + oaa → cit + coa + h
2.51 × 10−5
icit + nadp → akg + co2 + nadph
1.10 × 10−4
fum + h2o → mal-l
7.29 × 10−6
mal-l + nad → h + nadh + oaa
1.06 × 10−4
accoa + pi → actp + coa
8.49 × 10−6
lac-l + nad → h + nadh + pyr
3.60 × 10−6
3pg + nad → 3php + h + nadh
1.90 × 10−5
h + oxa → co2 + for
6.21 × 10−7
micit → pyr + succ
6.80 × 10−8
akg +h → co2 + sucsal
6.80 × 10−8
As for the absolute protein quantifications, the data reported by Goelzer et al.  were used. These values are expressed in number of molecules per cell and are obtained from LC/MSE analysis . This dataset covers most of the cytosolic proteins in the B. subtilis 168 strain growing under aerobic batch conditions and in minimal media with different carbon sources. In particular, the measurements in minimal medium with glucose were retrieved from this study and converted into units of millimoles per gram of dry cell weight (mmol/gDW) by assuming 6.3 × 108 cells per ml per optical density at 600 nm (OD600)  and 0.48 gDW/L per OD600 . For each enzyme with known kcat, we used the upper limit of the 95% confidence interval of protein abundance value in order not to over-constrain the model predictions (Table 1). When the protein level quantification of a specific enzyme was not available, we assumed that the protein is present at levels under the detection limit, and the minimum value among all the protein level measurements in the same condition (i.e., 6.8 × 10−8 mmol/gDW) was used.
Integration of enzymatic data in the model
Since we considered reactions catalyzed by unique enzymes, the number of enzyme-constrained reactions (17) is equal to the number of enzymes.
In summary, each constrained metabolic reaction Rj includes a pseudo-metabolite representing enzyme usage, which is limited by protein abundance.
Following the GECKO approach  to implement the described method, as a first step, the iYO844 model was converted into an irreversible model and the constraint for uptake rate of glucose, the sole carbon source, was removed. Then, the stoichiometric matrix and the upper bound vector of the model were expanded by adding the kcat values and the known protein abundances (Table 1).
The kinetic data, in terms of turnover number, were integrated in the stoichiometric matrix (S) as in the standard approach. When using this approach, the total amount of cellular proteins (Ptotal) in the cell was assumed to be 0.55 g/gDW, corresponding to the value measured for E. coli , and the mass fraction of the accounted proteins (f) was computed equal to 0.0191, by summing the abundance, expressed as parts per million (ppm), of the 17 considered proteins, retrieved from PaxDB database . All the computations were implemented via MathWorks MATLAB R2012a and run with COBRA toolbox  using the Gurobi solver .
For the in silico simulation of each mutant strain, the reaction encoded by the knocked-out gene was imposed as inactive, namely with a flux equal to zero. All strains were simulated using both the previously published iYO844 model and the new enzyme-constrained model obtained in this work, hereafter referred to as ec_iYO844. The lower bound of glucose uptake rate was fixed to the specific experimental value for the iYO844 model (− 7.71 mmol/gDW/h for wild-type strain), based only on stoichiometric reactions and directionality, whereas it was set to an unlimited value (− 100 mmol/gDW/h) for ec_iYO844. The metabolic phenotype of wild-type B. subtilis strain was predicted by parsimonious flux balance analysis (pFBA) . In addition, mutant strains were simulated using the minimization of metabolic adjustment (MoMA) method , which minimizes the distance between the flux distributions in wild-type and mutant strains.
The pFBA, MoMA and FVA methods were run in Matlab using the available functions in the COBRA Toolbox .
Identification of deletion targets
The gene deletions that are required to optimize the production of γ-PGA were identified by using the iYO844 model and ec_iYO844. The MoMA method was used to find single or multiple gene deletions corresponding to the best trade-off between growth rate and secretion rate of the target metabolite as previously described .
Since the pgs operon, including the enzyme-encoding genes responsible for γ-PGA production, is not expressed in the laboratory strain modeled in iYO844, we added the production reaction reported in the GEM of Bacillus licheniformis  (0.77 glu-D + 0.23 glu-L → γ-pga). As an alternative approach to optimize γ-PGA production, we considered the flux maximization of three of its known precursors: 2-oxoglutarate (akg), d-glutamate (glu-d) and l-glutamate (glu-l). In this case, the secretion reactions of these precursors were added to the model.
Evaluation of prediction accuracy
Flux distribution and growth rate data of B. subtilis wild-type and single-gene/operon deletion strains, grown under M9 minimal medium with glucose, were retrieved from literature [43, 44]. In these datasets, internal fluxes, measured by 13C-labeling experiments, are reported for the main reactions of the central carbon pathway, together with growth, glucose uptake and acetate production rates. For each measured flux the coefficient of variation among the replicates, when available in literature, is typically lower than 10%, with only the reactions of pentose phosphate pathway and malate dehydrogenase (MDH) of the tricarboxcylic acid (TCA) cycle having slightly higher variability (< 38%).
Furthermore, 95 single-gene deletions, corresponding to genes included in the GEM and experimentally found to be lethal in genome-wide studies [44, 46], were simulated via MoMA with the iYO844 and ec_iYO844 models, and the percentage of correctly predicted essential genes (i.e., yielding a predicted growth rate lower than 0.05 h−1) was calculated and compared to experimental data.
Additional model analyses
In order to evaluate the impact of kcat variation on prediction accuracy of the enzyme-constrained model, a distribution of prediction errors was computed by simulating ec_iYO844 using 10,000 kcat sets, randomly extracted without replacement from the list reported in Table 1.
In addition, to investigate the minimum set of reactions, among the ones reported in Table 1, that must be enzyme-constrained to achieve the final accuracy of ec_iYO844, the model was studied by following a stepwise inclusion procedure for the constrained reactions. In particular, unless differently indicated, the prediction errors for wild-type and mutant strains, together with the accuracy of essential genes, were all taken into account as indices to evaluate the accuracy.
B. subtilis strains used in this study
trpC2 pheA1 swrA +
trpC2 pheA1 swrA+ degU32(Hy) (Smr)
trpC2 pheA1 swrA+ ∆odhAB (Catr)
trpC2 pheA1 swrA+ ∆sucCD (Tetr)
trpC2 pheA1 swrA+ degU32(Hy) (Smr) ∆odhAB (Catr)
trpC2 pheA1 swrA+ degU32(Hy) (Smr) ∆sucCD (Tetr)
Growth conditions and fermentation
The PB5383, PB5691 and PB5716 strains from a streaked LB agar plate (tryptone, 10 g/L; yeast extract, 5 g/L; NaCl, 10 g/L; agar, 15 g/L) were grown for 7 h in 2.5 mL of Penassay broth (Difco Antibiotic Medium 3) with 5 g/L glucose in an orbital shaking incubator at 37 °C, 250 rpm; cultures were diluted to an OD600 of 0.1 in 20 mL of E medium adapted from Leonard et al.  (citric acid, 12 g/L; glucose, 80 g/L; NH4Cl, 7 g/L; K2HPO4, 0.5 g/L; MgSO4·7H2O, 0.5 g/L, FeCl3·6H2O, 0.04 g/L; CaCl2·2H2O, 0.15 g/L; MnSO4·H2O, 0.104 g/L; tryptophan, 50 µg/mL; phenylalanine, 50 µg/mL) with 40 g/L of l-glutamic acid (EM+) or without l-glutamic acid (EM−) and incubated under the same conditions as above. In fermentation experiments, aliquots (500 µL) were collected for OD600 readings and γ-PGA extraction at the following time points: 0, 18, 25, 43, 49, 66 and 74 h. In growth experiments, pre-cultures in EM− were grown for 16 h and then diluted to OD600 0.1; finally, OD600 was monitored for 9 h during the exponential growth phase, thereby enabling the calculation of growth rate as previously described .
γ-PGA recovery and quantification
Each culture sample was centrifuged (16,000 rpm, 20 min, 4 °C) and the cell pellet was discarded. The supernatant was precipitated with three volumes of cold methanol and kept at − 20 °C for at least 12 h. For γ-PGA quantification, each sample was centrifuged (14,000 rpm, 15 min, 4 °C) and methanol was discarded; γ-PGA pellets were dried in a vacuum concentrator and then dissolved in deionized water. The concentration of γ-PGA was measured by spectrophotometer at 216 nm in a quartz cuvette, as previously reported . Absorbance measurements were carried out via the NanoPhotometer UV/Vis spectrophotometer (Implen, Schatzbogen, München, Germany), using deionized water as blank, and a standard calibration curve with purified γ-PGA (#G1049, Sigma Aldrich, dissolved in deionized water) to compute its concentration in g/L from absorbance values. Measurements were carried out for at least three independent fermentation experiments.
The same aliquots of γ-PGA produced by wild-type and mutant strains were also visually compared by agarose gel electrophoresis followed by methylene blue staining, as previously described .
Microsoft Excel and Matlab functions were used to carry out statistical tests. Analysis of variance (ANOVA) was used to compare γ-PGA concentration or production rate among the three producer strains (PB5383, PB5691 and PB5716). A priori pairwise orthogonal comparisons were carried out via individual T-tests between PB5383 vs knockout strains and between the two knockout strains.
Wilcoxon signed rank test was used to compare flux variability values for all the model reactions with non-zero FV values. The metabolic subsystems (retrieved from the iYO844 original model) over-represented in the list of reactions with relevant flux variability reduction were found by Fisher’s exact test.
A new enzyme-constrained GECKO model of B. subtilis (ec_iYO844) was constructed by integrating publicly available enzyme kinetic and proteomic data for a set of central carbon and related pathway reactions (Table 1). The prediction performance of ec_iYO844 was compared to that of the original model (iYO844), in terms of ability to capture flux distributions in wild-type and single-gene/operon deletion mutants, and essential genes. The two models were also compared in terms of flux variability analysis, to evaluate if enzyme constraints were able to decrease the variability of fluxes in one or more pathways. Finally, the predicted knockout strategies to improve γ-PGA production were compared between the models and the suggested single-reaction deletions were evaluated in vivo.
Evaluation of model prediction performance on wild-type strain
Evaluation of model prediction performance on mutant strains
As expected, with the model developed in this work, the fluxes through the CS reaction in the TCA cycle decreased in every mutant strain, resulting into higher consistency with experimental values, due to the flux constraint of 4.4 mmol/gDW/h imposed by enzymatic information (derived from Table 1).
Taken together, the results showed that new model is also superior in terms of mutant strain flux predictions, although with a percent improvement in prediction accuracy lower than what was observed for the wild-type strain. The strong assumption that proteomic profile is the same in wild-type (for which proteomic data was available) and single-gene/operon mutants may be partly responsible of the less accurate predictions obtained for mutant strains. For the enzyme-constrained model, we also used an alternative approach, similar to the MOMENT method , in which the kcat values for the 17 reactions were integrated and only the total amount of enzymes (g/gDW) was constrained (kc_iYO844 model). The results showed that this method has lower prediction error (0.55 median value) than iYO844 (0.67 median value), but, also in this context, the predictions obtained by ec_iYO844 are the most accurate (0.43 median value). The three models (iYO844, ec_iYO844 and kc_iYO844) were also solved via pFBA instead of MoMA, but conclusions do not change and the resulting prediction errors with pFBA are always higher than the respective ones with MoMA (data not shown).
Gene essentiality analysis
Predicted with iYO844
Predicted with ec_iYO844
Impact of the knowledge of enzymatic data
The impact of the specific types of enzymatic data used to develop the enzyme-constrained model on model performance was analyzed in two tests.
First, among the 17 enzyme-constrained reactions, we identified triose-phosphate isomerase (TPI), glyceraldehyde-3-phosphate dehydrogenase (GAPD_NAD), CS and PGCDr as the minimum set required to achieve the final prediction performance of the model developed, namely a prediction error equal to 0.27 for the wild-type strain and to 0.43 for the mutants, and 76% of accuracy for central carbon gene essentiality. This result, however, is expected to depend on the study on which the model is applied, with prediction accuracy increasing with the number of enzyme-constrained reactions. This outcome is also suggested by our gene essentiality analysis, in which a significant improvement was observed for central carbon pathway but not in other pathways that were not subjected to enzymatic constraints.
By considering the prediction error for wild-type and mutant strains as the only performance measure, the GAPD_NAD and CS constraints were sufficient to reach the maximum predictive power. The importance of these constraints could be due to the key role of glyceraldehyde 3-phosphate (converted into 3-phospho-d-glyceroyl-phosphate via the GAPD_NAD reaction) as a central node between glycolysis and pentose phosphate pathway, and the key role of citrate synthesis (from oxaloacetate and acetyl-CoA via the CS reaction) as central node between glycolysis and the entering of TCA cycle. To better elucidate the role of these two reactions, we performed three simulations in which the GAPD_NAD and CS reactions were individually or jointly constrained as the sole constraint(s) of the ec_iYO844 model. The model reached its lowest error only by constraining both reactions. The individual constraint of GAPD_NAD or CS was not able to fix the PPP flux distribution, while the improvement in the TCA cycle fluxes was observed by individually including the CS constraint.
Conversely, considering the accuracy of gene essentiality analysis, the TPI and PGCDr constraints were needed to correctly predict pfkA, eno and pgm as essential genes, which catalyze the PFK, ENO and PGM reactions, respectively (Fig. 6). To better investigate the role of these two constraints, we performed simulations in which TPI or PGCDr were constrained. Results showed that the two constraints individually contributed to the correct prediction of pfkA (TPI constraint), and of pgm and eno (PGCDr constraint).
Flux variability analysis
In silico metabolic engineering to increase γ-PGA flux
Poly-γ-glutamic acid is a biodegradable, water-soluble, non-toxic, and edible biopolymer that has a large number of biotechnological applications, ranging from biomedicine to bioremediation. For its unique proprieties, different studies for improving the microbial production of γ-PGA were carried out [52, 53]. The pgs operon, which includes the enzyme-coding genes responsible for the biosynthesis of this polymer, is not expressed in B. subtilis 168 laboratory strains. However, we previously demonstrated that pgs operon transcription could be activated in a B. subtilis 168 derivative by two mutations: swrA+ and degU32(Hy), thereby engineering γ-PGA production . In this work, our goal was to use genome-scale models to identify the genetic perturbations (gene deletions) that could increase the production of γ-PGA.
Even if the same single-reaction deletion targets were identified by both models, and similar in silico γ-PGA flux and growth rate upon SUCOAS deletion were predicted by both models, the effects of AKGD deletion had a significant difference between the two models (see Fig. 7). ec_iYO844 predicted a fivefold higher γ-PGA flux than iYO844, along with a significant growth rate reduction. Using iYO844, the γ-PGA flux upon AKGD deletion was even lower than the one predicted with SUCOAS deletion.
Regarding the predicted flux distribution of the AKGD deletion mutant, a major difference between iYO844 and ec_iYO844 was that the former predicted the activation of the 2-oxoglutarate decarboxylase (OXGDC) and the succinate-semialdehyde dehydrogenase (SSALy) bypass reactions. While the SUCOAS and AKGD reactions are downstream of the akg production in the TCA cycle, OXGDC and SSALy (akg + h → co2 + sucsal, and h2o + nadp + sucsal → (2) h + nadph + succ, respectively) are consecutive reactions forming a bypass for the production of succinate (succ) from akg (see Fig. 6). Conversely, the ec_iYO844 model predicted this bypass reactions as inactive due to an enzymatic constraint on OXGDC (see Table 1), for which the associated protein (MenD) was not experimentally detected in these conditions [33, 55]. For this reason, when AKGD is deleted, ec_iYO844 predicted a higher flux from akg towards glutamate (and subsequently γ-PGA) than iYO844, in which akg was predicted to be mainly converted into succ and to continue the TCA cycle, with a lower flux towards glutamate. We confirmed that the bypass from akg to succ was the source of this key diversity between the models by simply imposing a null value of OXGDC in iYO844, as it was essentially done in ec_iYO844 by constraining the flux to very low levels. As expected, results predicted an increase of γ-PGA flux in the AKGD deletion strain (see Fig. 7).
Using the same procedure, we searched for double-reaction deletions to further compare the two models. As expected from the analysis above, for iYO844 the best double deletion found was the AKGD and OXGDC (or SSALy) combination (see Fig. 7). Using ec_iYO844, the three best deletions were the AKGD combination with the phosphotransacetylase (PTAr), acetate kinase (ACKr) or pyruvate dehydrogenase (PDH) reactions (see Fig. 7). The additional deletion of PTAr, ACKr or PDH predicted an increase of γ-PGA flux by decreasing the flux towards acetate (see Fig. 6), also consistent with a previous metabolic engineering work . A similar conclusion could be drawn by considering the best double-reaction deletion target (AKGD + PDH) found using iYO844 with a null value for the OXGDC flux as additional constraint (data not shown). This target combination was immediately followed by combinations of AKGD with a set of PPP reactions, with lower γ-PGA flux (data not shown), not found using ec_iYO844, probably due to the major differences in the flux distribution of this pathway (as illustrated in Fig. 1). The results obtained by double-deletion analysis further confirmed that the flux through the OXGDC/SSALy bypass is a key difference between iYO844 and ec_iYO844, affecting the predicted deletion targets, but other enzymatic constraints can also play a role for specific strain design predictions.
In summary, the two models showed the same results in terms of recommended single-reaction deletion list (AKGD or SUCOAS), but the predicted fluxes showed quantitatively different values due to the constraints of the GECKO model. Such constraints also determine the difference between the two models in the suggested double-reaction deletion list, with iYO844 recommending the deletion of a reaction (OXGDC or SSALy) for which ec_iYO844 predicted a negligible flux.
In vivo metabolic engineering to improve γ-PGA production
The addition of 40 g/L of l-glutamate (EM+ medium), a direct precursor of γ-PGA, as in the original E medium composition, significantly improved biopolymer production compared to the EM− medium for all strains (Fig. 8d), while the growth profiles remained very similar both in presence and in the absence of glutamate (dashed lines Fig. 8a). The reference strain produced up to 32 g/L, while odhAB and sucCD mutants reached up to 43 and 38 g/L, respectively (not statistically significant difference compared to the reference strain; T-test). The per-cell production rate showed a significant difference between the reference strain (0.4 g/gDW/h) and the mutants, and between the two mutants (p-value < 0.05; T-test), with the odhAB deletion strain showing again higher production (0.7 g/gDW/h) than the sucCD mutant (0.43 g/gDW/h). The differences observed between fermentations in EM− and EM+ were expected, since l-glutamate is a direct precursor of γ-PGA; its presence results in increased γ-PGA flux, as previously measured . For this reason, a smaller percent difference between the knockout and parent strains was also expected in EM+ than in EM−.
Following the principles of the recently proposed GECKO method, we integrated enzymatic and proteomic data for 17 central carbon-related pathways in an existing genome-scale metabolic model of B. subtilis. We showed that the resulting model, ec_iYO844, provides more accurate prediction of metabolic phenotypes compared to the original non-enzyme constrained model. Both secretion and intracellular flux predictions are considerably better with ec_iYO844 than iYO844 for wild-type and, to a lesser extent, mutant strains. Notably, the large errors in wild-type flux predictions observed in the original model, i.e. for the pentose phosphate pathway (in which no flux was predicted), acetate secretion and TCA cycle (which were under- and over-estimated, respectively), are all fixed in ec_iYO844. Moreover, the new model does not require to fix glucose uptake rate to an experimentally determined value, as was the case with the original model, since the glucose uptake flux is determined by enzyme capacity constraints. This advantage has been previously highlighted for other models adding constraints to existing genome-scale model . A dataset of essential genes was also used to compare the accuracy of the initial and new model: the enzyme-constrained one provided a correct prediction (i.e., no growth) for more strains than the non-enzyme constrained model. The specific values retrieved for kcat of enzymes are essential to achieve the improved prediction performance, since a randomized set of kcat values could not provide a similar accuracy. However, not all the enzyme constraints for reactions were essential to reach the predictive accuracy of the ec_iYO844 model: constraints for only four (TPI, GAPD_NAD, CS and PGCDr) of the 17 reactions were sufficient under the conditions that we used in this study. It is important to highlight that this result is expected to be specific for the data used to evaluate predictions.
While the originally proposed GECKO method adopted a genome-scale constraining of reactions in yeast, in this study only 17 reactions involved in central carbon metabolism were enzyme-constrained based on manually curated enzyme kinetic and proteomic data. This small-scale constraining was sufficient to demonstrate the potential of this approach with the experimental data used, but larger-scale constraining would be expected to provide benefits for the prediction of reactions and metabolic engineering targets outside of central carbon metabolism. The limitations of small-scale constraining were apparent in the gene essentiality predictions for genes outside of central carbon metabolism, which were not affected by enzyme constraints for central carbon metabolic reactions. The primary challenge for developing constraints for reactions outside of central carbon metabolism is the limited availability of enzyme kinetic data for even relatively well-characterized organisms like B. subtilis. The poor availability of high quality enzyme kinetic data has been identified as a major impediment for developing improved computational models of cellular metabolism in general .
In contrast to essentiality predictions, enzyme constraints showed significant effect in reducing flux variability in reactions outside of central carbon metabolism. Out of the reactions with non-zero variability in the original iYO844 model, 81% of the reactions showed reduced flux variability. This demonstrates that enzymatic constraints can have an impact on different pathways apart from the one in which constrained are added. The result obtained for ec_iYO844 is consistent with the improvements observed in yeast with genome-scale enzyme constraining. The reduction in flux variability suggests that the ec_iYO844 model could be better for strain design than the original iYO844 model even for products that do not directly come from central carbon metabolism.
The evaluation of the models in a metabolic engineering study enabled the identification of two promising gene deletion targets in TCA cycle that led to the improvement in the production of γ-PGA, a biopolymer with a large number of applications. The model without and with the integration of enzymatic data predicted the same single-reaction deletion targets (AKGD and SUCOAS). However, differences were observed in the quantitative prediction of γ-PGA flux of mutant strains and in the list of recommended double-reaction targets. Among these, a remarkable difference was that the original model led to the selection of AKGD with OXGDC/SSALy as best deletion combination, while the enzyme-constrained model predicted that OXGDC/SSALy do not carry significant flux due to an undetectable protein level for one of the corresponding enzymes. Therefore, the strain design based on the original model was unnecessarily complex since the same result could be reached with a single-reaction knockout (AKGD).
Considering the single-reaction deletion targets, the inverse correlation of γ-PGA production with the level of the enzymes catalyzing SUCOAS and AKGD was previously observed by Yu et al.  in B. licheniformis, but, to our knowledge, no forward engineering studies have been done to confirm the effect of these knockouts. For this reason, the AKGD and SUCOAS reactions (encoded by odhAB and sucCD operons, respectively) were selected as deletion candidates in this study to test their effect on γ-PGA biosynthesis.
The two knockout strains (ΔodhAB and ΔsucCD) were constructed in vivo by implementing gene deletions in a previously engineered strain (PB5383), capable of γ-PGA production. Fermentation experiments showed that both mutants could reach a significantly higher (twofold) maximum γ-PGA concentration than PB5383 over a three-day growth in E medium without glutamate. The ΔodhAB strain also significantly outperforms the others in terms of per-cell γ-PGA production rate, with a two-fold improvement over PB5383. These results demonstrated that TCA cycle flux was successfully diverted towards glutamate, and subsequently toward γ-PGA production, by implementing the knockouts recommended by the model. As expected, the addition of l-glutamate in the medium increased γ-PGA production in all strains, since l- and d-glutamate are precursors of γ-PGA, but the ΔodhAB strain still showed improved production compared with the other strains.
In summary, we have developed a B. subtilis genome-scale metabolic model with significantly improved ability to predict metabolic flux distributions and phenotypes. The improved model relies only on a limited set of enzymatic constraints derived from measured protein expression profiles and kinetic parameters reported in the scientific literature. The application of this model to predicting metabolic phenotypes and metabolic engineering designs demonstrated the utility of the model and the value of relatively limited enzymatic constraints. This work paves the way for improving genome-scale metabolic models of other industrially attractive organisms by incorporating enzyme kinetic and proteomic data available in databases and the scientific literature.
MH conceived the study. IM performed the modeling and simulation work. IM, ER and MC performed the wetlab experiments. IM, LP, NS, PM, CC and MH designed the experiments, analyzed the data and wrote the manuscript. PM, CC and MH supervised the study. PM and CC provided reagents and materials. All authors read and approved the final manuscript.
The authors want to thank Jörg Stülke Lab for providing the genomic DNA of strains GP791 and GP1276, and Alessandra Albertini for her help in spectrophotometric assays.
The authors declare that they have no competing interests.
Availability of data and materials
The metabolic model and the experimental flux data used in this study are available in the GitHub repository (https://github.com/biosustain/ec_iYO844). The data collected in fermentation experiments (growth and γ-PGA concentration data) are available from the corresponding author on request.
Consent for publication
Ethics approval and consent to participate
This work was supported by Fondazione Cariplo through the Grant 2015–0397 “Conversion of industrial bio-waste into biofuels and bioproducts through synthetic biology”, Grant 2015–0400 “Rice straw valorization: recovery of inorganic and organic components” and by the Italian Ministry of Education, University and Research (MIUR): Dipartimenti di Eccellenza Program (2018–2022)—Dept. of Biology and Biotechnology "L. Spallanzani", University of Pavia. MJH and NS acknowledge funding from the European Union’s Horizon 2020 research and innovation program under Grant Agreement No. 686070 and the Novo Nordisk Foundation. The funders had no role in study design, data collection, analysis and interpretation, decision to publish, or preparation of the manuscript.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Kunst F, et al. The complete genome sequence of the Gram-positive bacterium Bacillus subtilis. Nature. 1997;390(6657):249–56.PubMedGoogle Scholar
- Muntel J, Fromion V, Goelzer A, Maabeta S, Mader U, Buttner K, Hecker M, Becher D. Comprehensive absolute quantification of the cytosolic proteome of Bacillus subtilis by data independent, parallel fragmentation in liquid chromatography/mass spectrometry (LC/MSE). Mol Cell Proteom. 2014;13:1008–199.Google Scholar
- Nicolas P, Mäder U, Dervyn E, Rochat T, Leduc A, Pigeonneau N, et al. Condition-dependent transcriptome reveals high-level regulatory architecture in Bacillus subtilis. Science. 2012;335(6072):1103–6.PubMedGoogle Scholar
- Hahne H, Mader U, Otto A, Bonn F, Steil L, Bremer E, Hecker M, Becher D. A comprehensive proteomics and transcriptomics analysis of Bacillus subtilis salt stress adaptation. J Bacteriol. 2010;192(3):870–82.PubMedGoogle Scholar
- Van Duy N, Mader U, Tran NP, Cavin JF, Tam LT, Albrecht D, Hecker M, Antelmann H. The proteome and transcriptome analysis of Bacillus subtilis in response to salicylic acid. Proteomics. 2007;7(5):698–710.PubMedGoogle Scholar
- Adams BL. The next generation of synthetic biology chassis: moving synthetic biology from the laboratory to the field. ACS Synth Biol. 2016;5:1328–30.PubMedGoogle Scholar
- Guiziou S, Sauveplane V, Chang HJ, Clerté C, Declerck N, Jules M, Bonnet J. A part toolbox to tune genetic expression in Bacillus subtilis. Nucleic Acids Res. 2016;44(15):7495–508.PubMedPubMed CentralGoogle Scholar
- Xiao-ying G, Chun-e H, Tao L, Zhu O. Effect of Bacillus subtilis and Pseudomonas fluorescens on growth of greenhouse tomato and rhizosphere microbial community. J Northeast Agric Univ (English Edition). 2015;22(3):32–42.Google Scholar
- Qiao J, Yu X, Liang X, Liu Y, Borriss R, Liu Y. Addition of plant-growth-promoting Bacillus subtilis PTS-394 on tomato rhizosphere has no durable impact on composition of root microbiome. BMC Microbiol. 2017;17(1):131.PubMedPubMed CentralGoogle Scholar
- Perkins JB, Pero JG, Sloma A. Riboflavin overproducing strains of bacteria. European Patent; 1991.Google Scholar
- Perkins JB, Sloma A, Hermann T, Theriault K, Zachgo E, Erdenberger T, Hannett N, Chatterjee NP, Williams V II, Rufo GA Jr, Hatch R, Pero J. Genetic engineering of Bacillus subtilis for the commercial production of riboflavin. J Ind Microbiol Biotechol. 1999;22:8–18.Google Scholar
- Wang JP, Yeh CM, Tsai YC. Improved subtilisin YaB production in Bacillus subtilis using engineered synthetic expression control sequences. J Agric Food Chem. 2006;54(25):9405–10.PubMedGoogle Scholar
- Navaneeth S, Bhuvanesh S, Bhaskar V, Vijay KP, Kandaswamy SKJ, Achary A. Optimization of medium for the production of subtilisin from Bacillus subtilis MTCC 441. Afr J Biotechnol. 2009;8:22.Google Scholar
- Gilbert C, Howarth M, Harwood CR, Ellis T. Extracellular self-assembly of functional and tunable protein conjugates from Bacillus subtilis. ACS Synth Biol. 2017;6:957–67.PubMedGoogle Scholar
- Osera C, Amati G, Calvio C, Galizzi A. SwrAA activates poly-gamma-glutamate synthesis in addition to swarming in Bacillus subtilis. Microbiology. 2009;155:2282–7.PubMedGoogle Scholar
- Scoffone V, Dondi D, Biino G, Borghese G, Pasini D, Galizzi A, Calvio C. Knockout of pgdS and ggt genes improves γ-PGA yield in B. subtilis. Biotechnol Bioeng. 2013;110(7):2006–122.PubMedGoogle Scholar
- Oh YK, Palsson BO, Park SM, Schilling CH, Mahadevan R. Genome-scale reconstruction of metabolic network in Bacillus subtilis based on high-throughput phenotyping and gene essentiality data. J Biol Chem. 2007;282:28791–9.PubMedGoogle Scholar
- Henry CS, Zinner JF, Cohoon MP, Stevens RL. iBsu1103: a new genome-scale metabolic model of Bacillus subtilis based on SEED annotations. Genome Biol. 2009;10(6):R69.PubMedPubMed CentralGoogle Scholar
- Tanaka K, Henry CS, Zinner JF, Jolivet E, Cohoon MP, Xia F, Bidnenko V, Ehrlich SD, Stevens RL, Noirot P. Building the repertoire of dispensable chromosome regions in Bacillus subtilis entails major refinement of cognate large-scale metabolic model. Nucleic Acids Res. 2012;41(1):687–99.PubMedPubMed CentralGoogle Scholar
- Hao T, Han B, Ma H, Fu J, Wang H, Wang Z, Tang B, Tao C, Zhao X. In silico metabolic engineering of Bacillus subtilis for improved production of riboflavin, Egl-237,(R,R)-2,3-butanediol and isobutanol. Mol BioSyst. 2013;9(8):2034–44.PubMedGoogle Scholar
- Kalantari A, Chen T, Ji B, Stancik IA, Ravikumar V, Franjevic D, Saulou-Bérion C, Goelzer A, Mijakovic I. Conversion of glycerol to 3-hydroxypropanoic acid by genetically engineered Bacillus subtilis. Front Microbiol. 2017;8:638.PubMedPubMed CentralGoogle Scholar
- Beg QK, Vazquez A, Ernst J, de Menezes MA, Bar-Joseph Z, Barabási AL, Oltvai ZN. Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity. Proc Natl Acad Sci. 2007;104(31):12663–8.PubMedGoogle Scholar
- Adadi R, Volkmer B, Milo R, Heinemann M, Shlomi T. Prediction of microbial growth rate versus biomass yield by a metabolic network with kinetic parameters. PLoS Comput Biol. 2012;8(7):e1002575.PubMedPubMed CentralGoogle Scholar
- Nilsson A, Nielsen J. Metabolic trade-offs in yeast are caused by F1F0-ATP synthase. Sci Rep. 2016;6:22264.PubMedPubMed CentralGoogle Scholar
- Yizhak K, Benyamini T, Liebermeister W, Ruppin E, Shlomi T. Integrating quantitative proteomics and metabolomics with a genome-scale metabolic network model. Bioinformatics. 2010;26(12):i255–i260260.PubMedPubMed CentralGoogle Scholar
- Sanchez BJ, Zhang C, Nilsson A, Lahtvee PJ, Kerkhoven EJ, Nielsen J. Improving the phenotype predictions of a yeast genome-scale metabolic model by incorporating enzymatic constriants. Mol Syst Biol. 2017;13(8):935.PubMedPubMed CentralGoogle Scholar
- Schomburg I, et al. BRENDA, the enzyme database: updates and major new developments. Nucleic Acids Res. 2004;3(suppl_1):D431–D433433.Google Scholar
- Wittig U, Kania R, Golebiewski M, Rey M, Shi L, Jong L, Algaa E, Weidemann A, Sauer-Danzwith H, Mir S, Krebs O, Bittkowski M, Wetsch E, Rojas I, Müller W. SABIO-RK-database for biochemical reaction kinetics. Nucleic Acids Res. 2012;40:D790–D796796.PubMedGoogle Scholar
- Davidi D, Milo R. Lessons on enzyme kinetics from quantitative proteomics. Curr Opin Biotechnol. 2017;46:81–9.PubMedGoogle Scholar
- Jin S, Sonenshein AL. Characterization of the major citrate synthase of Bacillus subtilis. J Bacteriol. 1996;178(12):3658–60.PubMedPubMed CentralGoogle Scholar
- Costa T, Steil L, Martins LO, Völker U, Henriques AO. Assembly of an oxalate decarboxylase produced under σK control into the Bacillus subtilis spore coat. J Bacteriol. 2004;186(5):1462–74.PubMedPubMed CentralGoogle Scholar
- Davidi D, Noor E, Liebermeister W, Bar-Even A, Flamholz A, Tummler K, Barenholz U, Goldenfeld M, Sholmi T, Milo R. Global characterization of in vivo enzyme catalytic rates and their correspondence to in vitro kcat measurements. Proc Natl Acad Sci. 2016;113(12):3401–6.PubMedGoogle Scholar
- Goelzer A, Muntel J, Chubukov V, Jules M, Prestel E, Nölker R, Mariasdassou M, Aymerich S, Hecker M, Noirot P, Becher D, Fromion V. Quantitative prediction of genome-wide resource allocation in bacteria. Metab Eng. 2015;32:232–43.PubMedGoogle Scholar
- Milo R. What is the total number of protein molecules per cell volume? A call to rethink some published values. BioEssays. 2013;35(12):1050–5.PubMedPubMed CentralGoogle Scholar
- Wang M, Weiss M, Simonovic M, Haertinger G, Schrimpf SP, Hengartner MO, von Mering C. PaxDb, a database of protein abundance averages across all three domains of life. Mol Cell Proteomics. 2012;11(8):492–500.PubMedPubMed CentralGoogle Scholar
- Schellenberger J, Que R, Fleming RM, Thiele I, Orth JD, Feist AM, Rahmanian S, Kang J, Hyduke DR, Palsson BØ. Quantitative prediction of cellular metabolism with constraint-based models: the COBRA Toolbox v2.0. Nat Protoc. 2011;6(9):1290–307.PubMedPubMed CentralGoogle Scholar
- Gurobi Optimization, LLC. Gurobi optimizer reference manual. 2018. https://www.gurobi.com.
- Lewis NE, Hixson KK, Conrad TM, Lerman JA, Charusanti P, Polpitiya AD, Adkins JN, Schramm G, Purvine SO, Lopez-Ferrer D, Weitz KK, Eils R, Konig R, Smith RD, Palsson BØ. Omic data from evolved E. coli are consistent with computed optimal growth from genome-scale models. Mol Syst Biol. 2010;6(1):390.PubMedPubMed CentralGoogle Scholar
- Segre D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci. 2002;99(23):15112–7.PubMedGoogle Scholar
- Mahadevan R, Schilling CH. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metab Eng. 2003;5(4):264–76.Google Scholar
- Park JH, Lee KH, Kim TY, Lee SY. Metabolic engineering of Escherichia coli for the production of L-valine based on transcriptome analysis and in silico gene knockout simulation. Proc Natl Acad Sci. 2007;104(19):7797–802.PubMedGoogle Scholar
- Guo J, Zhang H, Wang C, Chang JW, Chen LL. Construction and analysis of a genome-scale metabolic network for Bacillus licheniformis WX-02. Res Microbiol. 2016;167(4):282–9.PubMedGoogle Scholar
- Chubukov V, Uhr M, Le Chat L, Kleijn RJ, Jules M, Link H, Aymerich S, Stelling J, Sauer U. Transcriptional regulation is insufficient to explain substrate-induced flux changes in Bacillus subtilis. Mol Syst Biol. 2013;9(1):709.PubMedPubMed CentralGoogle Scholar
- Fischer E, Sauer U. Large-scale in vivo flux analysis shows rigidity and suboptimal performance of Bacillus subtilis metabolism. Nat Genet. 2005;37(6):636–40.PubMedGoogle Scholar
- Machado D, Herrgård MJ. Systematic evaluation of methods for integration of transcriptomic data into constraint-based models of metabolism. PLoS Comput Biol. 2014;10(4):e1003580.PubMedPubMed CentralGoogle Scholar
- Kobayashi K, Ehrlich SD, Albertini A, Amati G, Andersen KK, Arnaud M, Asai K, Ashikaga S, Aymerich S, Bessieres P, et al. Essential Bacillus subtilis genes. Proc Natl Acad Sci. 2003;100(8):4678–83.PubMedGoogle Scholar
- Smith JL, Goldberg JM, Grossman AD. Complete genome sequences of Bacillus subtilis subsp. subtilis laboratory strains JH642 (AG174) and AG1839. Genome Announc. 2014;2(4):e00663-14.PubMedPubMed CentralGoogle Scholar
- Leonard CG, Housewright RD, Thorne CB. Effects of some metallic ions on glutamyl polypeptide synthesis by Bacillus subtilis. J Bacteriol. 1958;76:499–503.PubMedPubMed CentralGoogle Scholar
- Zucca S, Pasotti L, Politi N, Casanova M, Mazzini G, De Angelis MGC, Magni P. Multi-faceted characterization of a novel LuxR-repressible promoter library for Escherichia coli. PLoS ONE. 2015;10(5):e0126264.PubMedPubMed CentralGoogle Scholar
- Zeng W, Chen G, Zhang Y, Wu K, Liang Z. Studies on the UV spectrum of poly (γ-glutamic acid) based on development of a simple quantitative method. Int J Biol Macromol. 2012;51(1–2):83–90.PubMedGoogle Scholar
- Mamberti S, Prati P, Cremaschi P, Seppi C, Morelli CF, Galizzi A, Fabbi M, Calvio C. γ-PGA hydrolases of phage origin in Bacillus subtilis and other microbial genomes. PLoS ONE. 2015;10(7):e0130810.PubMedPubMed CentralGoogle Scholar
- Luo Z, Guo Y, Liu J, Qiu H, Zhao M, Zou W, Li S. Microbial synthesis of poly-γ-glutamic acid: current progress, challenges, and future perspectives. Biotechnol Biofuels. 2016;9(1):134.PubMedPubMed CentralGoogle Scholar
- Cao M, Feng J, Sirisansaneeyakul S, Song C, Chisti Y. Genetic and metabolic engineering for microbial production of poly-γ-glutamic acid. Biotechnol Adv. 2018;36:1424–33.PubMedGoogle Scholar
- Yu W, Chen Z, Shen L, Wang Y, Li Q, Yan S, Zhong C-J, He N. Proteomic profiling of Bacillus licheniformis reveals a stress response mechanism in the synthesis of extracellular polymeric flocculants. Biotechnol Bioeng. 2016;113(4):797–806.PubMedGoogle Scholar
- Dawson A, Chen M, Fyfe PK, Guo Z, Hunter WN. Structure and reactivity of Bacillus subtilis MenD catalyzing the first committed step in menaquinone biosynthesis. J Mol Biol. 2010;401(2):253–64.PubMedPubMed CentralGoogle Scholar
- Feng J, Gu Y, Quan Y, Cao M, Gao W, Zhang W, Wang W, Song C. Improved poly-γ-glutamic acid production in Bacillus amyloliquefaciens by modular pathway engineering. Metab Eng. 2015;32:106–15.PubMedGoogle Scholar
- Konglom N, Chuensangjun C, Pechyen C, Sirisansaneeyakul S. Production of poly-γ-glutamic acid by Bacillus licheniformis: synthesis and characterization. J Metals Mater Miner. 2012;22:2.Google Scholar
- Nilsson A, Nielsen J, Palsson BO. Metabolic models of protein allocation call for the kinetome. Cell Syst. 2017;5(6):538–41.PubMedGoogle Scholar
- Gao H, Chen Y, Leary JA. Kinetic measurements of phosphoglucose isomerase and phosphomannose isomerase by direct analysis of phosphorylated aldose–ketose isomers using tandem mass spectrometry. Int J Mass Spectrom. 2005;240(3):291–9.Google Scholar
- Wierenga RK, Kapetaniou EG, Venkatesan R. Triosephosphate isomerase: a highly evolved biocatalyst. Cell Mol Life Sci. 2010;67(23):3961–82.PubMedGoogle Scholar
- Fillinger S, Boschi-Muller S, Azza S, Dervyn E, Branlant G, Aymerich S. Two glyceraldehyde-3-phosphate dehydrogenases with opposite physiological roles in a nonphotosynthetic bacterium. J Biol Chem. 2000;275(19):14031–7.PubMedGoogle Scholar
- D'Alessio G, Josse J. glyceraldehyde phosphate dehydrogenase, phosphoglycerate kinase, and phosphoglyceromutase of Escherichia coli simultaneous purification and physical properties. J Biol Chem. 1971;246(13):4319–25.PubMedGoogle Scholar
- Watabe K, Freese E. Purification and properties of the manganese-dependent phosphoglycerate mutase of Bacillus subtilis. J Bacteriol. 1979;137(2):773–8.PubMedPubMed CentralGoogle Scholar
- Brown CK, Kuhlman PL, Mattingly S, Slates K, Calie PJ, Farrar WW. A model of the quaternary structure of enolases, based on structural and evolutionary analysis of the octameric enolase from Bacillus subtilis. J Protein Chem. 1998;17(8):855–66.PubMedGoogle Scholar
- Olavarría K, Valdes D, Cabrera R. The cofactor preference of glucose-6-phosphate dehydrogenase from Escherichia coli-modeling the physiological production of reduced cofactors. FEBS J. 2012;279(13):2296–309.PubMedGoogle Scholar
- Singh SK, Miller SP, Dean A, Banaszak LJ, LaPorte DC. Bacillus subtilis isocitrate dehydrogenase—a substrate analogue for Escherichia coli isocitrate dehydrogenase kinase/phosphatase. J Biol Chem. 2002;277(9):7567–73.PubMedGoogle Scholar
- Ueda Y, Yumoto N, Tokushige M, Fukui K, Ohya-Nishiguchi H. Purification and characterization of two types of fumarase from Escherichia coli. J Biochem. 1991;109(5):728–33.PubMedGoogle Scholar
- Smith K, Sundaram TK, Kernick M, Wilkinson AE. Purification of bacterial malate dehydrogenases by selective elution from a triazinyl dye affinity column. Biochim Biophys Acta (BBA) Protein Struct Mol Enzymol. 1982;708(1):17–25.Google Scholar
- Shin BS, Choi SK, Park SH. Regulation of the Bacillus subtilis phosphotransacetylase gene. J Biochem. 1999;126(2):333–9.PubMedGoogle Scholar
- Garvie EI. Bacterial lactate dehydrogenases. Microbiol Rev. 1980;44(1):106.PubMedPubMed CentralGoogle Scholar
- Saski R, Pizer LI. Regulatory properties of purified 3-phosphoglycerate dehydrogenase from Bacillus subtilis. FEBS J. 1975;51(2):415–27.Google Scholar
- Svedruzic D, Liu Y, Reinhardt LA, Wroclawska E, Cleland WW, Richards NG. Investigating the roles of putative active site residues in the oxalate decarboxylase from Bacillus subtilis. Arch Biochem Biophys. 2007;464(1):36–47.PubMedPubMed CentralGoogle Scholar
- Liu S, Lu Z, Han Y, Melamud E, Dunaway-Mariano D, Herzberg O. Crystal structures of 2-methylisocitrate lyase in complex with product and with isocitrate inhibitor provide insight into lyase substrate specificity, catalysis and evolution. Biochemistry. 2005;44(8):2949–62.PubMedGoogle Scholar
- Calvio C, Celandroni F, Ghelardi E, Amati G, Salvetti S, Ceciliani F, Galizzi A, Senesi S. Swarming differentiation and swimming motility in Bacillus subtilis are controlled by swrA, a newly identified dicistronic operon. J Bacteriol. 2005;187(15):5356–66.PubMedPubMed CentralGoogle Scholar