Activating Silent Glycolysis Bypasses in Escherichia coli

All living organisms share similar reactions within their central metabolism to provide precursors for all essential building blocks and reducing power. To identify whether alternative metabolic routes of glycolysis can operate in E. coli, we complementarily employed in silico design, rational engineering, and adaptive laboratory evolution. First, we used a genome-scale model and identified two potential pathways within the metabolic network of this organism replacing canonical Embden-Meyerhof-Parnas (EMP) glycolysis to convert phosphosugars into organic acids. One of these glycolytic routes proceeds via methylglyoxal and the other via serine biosynthesis and degradation. Then, we implemented both pathways in E. coli strains harboring defective EMP glycolysis. Surprisingly, the pathway via methylglyoxal seemed to immediately operate in a triosephosphate isomerase deletion strain cultivated on glycerol. By contrast, in a phosphoglycerate kinase deletion strain, the overexpression of methylglyoxal synthase was necessary to restore growth of the strain. Furthermore, we engineered the “serine shunt” which converts 3-phosphoglycerate via serine biosynthesis and degradation to pyruvate, bypassing an enolase deletion. Finally, to explore which of these alternatives would emerge by natural selection, we performed an adaptive laboratory evolution study using an enolase deletion strain. Our experiments suggest that the evolved mutants use the serine shunt. Our study reveals the flexible repurposing of metabolic pathways to create new metabolite links and rewire central metabolism.


Introduction
Excluding some marked exceptions (e.g., [1][2][3]), the central carbon metabolism of all living organisms can be divided into few interacting pathways, i.e., glycolysis, the pentose phosphate pathway, and the tricarboxylic acid (TCA) cycle, each representing highly conserved reaction patterns. The emergence of these conserved patterns can be explained in two complementary ways. First, the structure of central metabolism mostly reflects a metabolic network that could have emerged under primordial conditions where primitive sys-tems prevailed solely by relying on abiotic catalysis [4,5]. This primordial network could have served as the backbone supporting the emergence of the first cells and, as such, was "frozen" at the origin of life, offering only little flexibility for adaptation throughout the specific evolutionary trajectory of current organisms. Alternatively, and complementarily, central metabolism might constitute an optimal solution for interconverting essential cellular metabolites under a given set of biochemical constraints, including pathway length, favorable thermodynamics, proper kinetics, chemical properties of metabolic intermediates, avoidance of radical enzymes, and toxic intermediates [6][7][8][9][10]. In this scenario, it might be that many pathway solutions have existed which were outcompeted by the ones providing optimal ATP yield (Embden-Meyerhof-Parnas (EMP)) or highest rates (Entner-Doudoroff (ED)).
The universal structure of central metabolism facilitates the study of cellular physiology by enabling to apply knowledge gained from one organism to the understanding of another. On the other hand, this conserved metabolism also severely restricts the metabolic and chemical space we can easily explore and impacts bioproduction by limiting the set of key starting metabolites. Therefore, there is a growing interest in restructuring central metabolism for exploring alternative metabolic networks and pave the way to new bioproduction capabilities, e.g., by replacing EMP glycolysis with the stoichiometrically favorable nonoxidative glycolysis [11][12][13][14]. Such novel routes can recruit enzymes from heterologous sources or even integrate new-to-nature metabolic conversions [15]. A possible strategy of this sort is to exploit the native enzymes and pathways of an organism for implementing new pathways that can replace segments of central metabolism. There are two main advantages for this approach. From a practical point of view, it relies only on enzymes that were optimized during evolution to operate within the desired cellular environment. From a scientific perspective, the recruitment of enzymes belonging to different cellular processes to establish novel metabolic networks is akin to the emergence of pathways during evolution and thus provides a platform to explore such key evolutionary events [16].
Here, we followed a systems biology approach of exploring if Escherichia coli's metabolic network contains alternative routes bypassing EMP glycolysis. Engineering these pathways can aid in defining the barriers in which metabolism can operate, providing a better understanding of why glycolysis is shaped in the way it is. Furthermore, we hoped to extend the available design space synthetic metabolism can operate in to tailor metabolism for biotechnological needs. Thus, in this study, we aimed to fully replace the endogenous EMP glycolysis by relying on native enzymes only. An in silico analysis identified multiple possible pathways allowing the required conversion of phosphosugars into pyruvate. We engineered two of these routes in E. coli strains lacking key EMP glycolysis enzymes. First, we show that a methylglyoxal-dependent route [17][18][19] can carry all glycolytic flux to support a high growth rate. We also implemented the "serine shunt," which bypasses glycolysis via serine biosynthesis and deamination to pyruvate. We demonstrate that the operation of this synthetic route relies on a delicate balance between the rate of serine production and consumption, avoiding the inhibitory accumulation of this amino acid. In parallel, we submitted an E. coli strain lacking EMP glycolysis to adaptive evolution, which selected for the emergence of the serine shunt after~140 days of cultivation.

Computational Analysis to Identify Glycolytic
Bypasses in E. coli. In order to identify possible pathways that can bypass EMP glycolysis resorting to native E. coli enzymes only, we used a constraint-based approach based on the genomescale metabolic model of this bacterium (iML1515)) except for a few alterations described in supplementary methods (section titled "Alterations to the iML1515 model").
This approach is very similar to the one used recently for finding latent carbon fixation cycles in E. coli [20]. Namely, we set up a mixed integer linear problem-(MILP-) based optimization problem which simultaneously looks for solutions that balance an objective reaction (here, glycerol to pyruvate) while minimizing the total number of reactions.
The full MILP is formulated as follows: Minimize Such that v obj = −1, where the vector v contains the relative reaction rates, z are the Boolean reaction indicators, x are the log-scaled metabolite concentrations, and g 0 is a vector of all the reactions' standard Gibbs free energies in units of RT, i.e., g i 0 = Δ r G ′ 0ðiÞ/RT (R is the gas constant and T is the temperature in Kelvin). v obj is a new reaction that we add to the iML1515 model with the following chemical formula: glycerol + n ADP + n P i ⟶ pyruvate + n ATP (where n is the required ATP yield). In the first constraint, the rate of the objective reaction (v obj ) is set to be exactly -1, ensuring that any pathway solution would exactly balance it, and therefore, the overall reaction in the pathway would be glycerol to pyruvate (and the required number of ATPs). All the other constraints are based on the standard thermodynamic fluxbalance analysis algorithm [21]. Together, these set of constraints allow only solutions that satisfy mass balance and the second law of thermodynamics, within a given range of flux and metabolite concentration bounds. The parameter β (which we set arbitrarily to 10) sets the upper bound for the rate of each single reaction in the pathway relative to the objective reaction. The lower and upper bounds on the concentrations of most metabolites are set to 1 μM and 10 mM, with the exception of 15 central metabolites and cofactors which are confined to more specific ranges based on physiological data (see Supplementary Table 4).
The objective function (not to be confused with the objective reaction) of the MILP is the sum of all z i , i.e., the total number of active reactions in the solution. Note that as a preprocessing step, we split all reactions to a forward 2 BioDesign Research and backward reaction and therefore all (unidirectional) rates must be positive. The constraint v -βz ≤ 0 ensures that each reaction indicator (z i ) can be equal to 0, only if the rate is 0. We do not need to care about it being equal to 1 even if a reaction is not active, since the optimization goal (which minimizes the sum of all indicators) will prevent that. Running the MILP solver once will yield the optimal pathway in terms of overall length (number of activate reactions). In order to explore the space of possible (suboptimal) glycolysis bypasses, we then iterate the solution space using integer cuts to eliminate all the potential solutions within a radius of 2 around each one of the previously found solutions (using the same method as in [20]). We stopped the search after 10 solutions were found (for each possible ATP yield).  [23], which was used as wildtype reference ( Table 1). The deletions were carried out by λ-Red recombineering using kanamycin resistance cassettes generated via PCR using the FRT-PGK-gb2-neo-FRT (Km) cassette (Gene Bridges, Germany) for deletion of sdaA, sdaB, tdcB, tdcG, and mgsA. For the deletion of dld, gldA, gloA, aldA, hchA, and ppsA, pKD3 (chloramphenicol) and pKD4 (kanamycin) were used as a template for amplification of deletion cassettes [24]. Primer pairs used are indicated in Supplementary Table S3. Cell preparation and transformation for gene deletion was carried out as described [23,25]. The coding sequences of the WT sequences of serA and the mutated genes were amplified by PCR using the primer pairs serA-pet-F and serA-pet-R (Supplementary Table S3). The amplified fragments were cloned into a modified pET16b expression vector (Table 1) by using In-Fusion cloning kit (Takara, Shiga, Japan). The sequence of the inserts of the resulting plasmids was verified by Sanger sequencing. For exchanging the native promoter of sdaA with a constitutive strong promoter, a chloramphenicol resistance cassette was amplified from pKD3 by using the primer pair SdaA-ProEx-F and CAP-sdaA-R. In a second PCR, the promotor sequence (5 ′ -ACCTATTGACAATTAAAGGCTAAAAT GCTATAATTCCAC-3′ [25]) was amplified from pZ-ASS using the primer pair pS-bridge and SdaA-ProEx-R. The purified PCR products were used in a fusion PCR together with the primer pair SdaA-ProEx-F and SdaA-ProEx-R. This resulted in a promoter exchange cassette containing 50 bp flanks to integrate into the intergenic region between nudL and sdaA. To integrate the feedback-resistant version of 3-phosphoglycerate dehydrogenase (SerA * , H344A N346A N364A) and simultaneously replace the native promotor of the gene with a strong constitutive one, a chloramphenicol resistance cassette was amplified from pKD3 using primer pair serA * -ProEx-F and CAP-SerA * -R. A PCR product containing a constitutive strong promotor (5′-ACCTATTGACAATTAAAGGCTAAAAT GCTATAATTCCAC-3′ [25]) and the ser A * gene was amplified using primer pair pS-bridge and serA * -ProEx-R. The purified PCR products were used in a fusion PCR together with the primer pair serA * -ProEx-F and serA * -ProEx-R, resulting in a chloramphenicol cassette containing ser A * behind a strong promotor and a 50 bp flank upstream of the CAT cassette to integrate into the intergenic region behind rpiA.

Evolution in GM3-Driven Long-Term Continuous
Culture. For evolution experiments, precultures of Δeno strain were obtained in permissive minimal MA medium supplemented with 20 mM glycerol and 10 mM succinate. The preculture was used to inoculate the growth chambers (16 mL per chamber) of two parallel independent GM3 devices [26]. A continuous gas flow of sterile air through the culture vessel ensured constant aeration and growth in suspension by counteracting cell sedimentation. The cultures were grown in the corresponding medium under turbidostat mode (dilution threshold set to 80% transmittance (OD ≈ 0:4, 37°C) until stable growth of the bacterial population. The cultures were then submitted to a conditional medium swap regime. This regime enabled gradual adaptation of the bacterial populations to grow in a nonpermissive medium which contained 20 mM glycerol only. Dilutions of the growing cultures were triggered every 10 minutes with a fixed volume of medium calculated to impose a generation time of 3.5 hours on the cell population, if not otherwise stated. The growing cultures were fed by permissive or nonpermissive medium depending on the turbidity of the culture with respect to a set OD threshold (OD 600 value of 3 BioDesign Research 0.4). When the OD exceeded the threshold, a pulse of nonpermissive medium was injected; otherwise, a pulse of permissive medium. When the cultures grew on glycerol only (100% nonpermissive medium), the feeding mode was set to turbidostat to increase growth rates. Three isolates were obtained on agar-solidified medium containing glycerol as sole carbon source from both evolution experiments and further analyzed.
2.6. Genomic Analysis of Evolved Strains. Pair-end libraries (2 × 150 bp) were prepared with 1 μg genomic DNA from the evolved strains as well as from the ancestor Δeno strain and sequenced using a MiSeq sequencer (Illumina). The PALOMA pipeline, integrated in the platform Microscope (https://mage.genoscope.cns.fr/microscope/home/), was used to map the reads against E. coli K12 wildtype strain MG1655 reference sequence (NC_000913.3) for detecting . breseq pipeline [28] was applied to map the reads against the reference for identification of genomic variants, including SNPs and insertion-deletion polymorphisms (INDELs).

Growth
Experiments. 4 mL M9 medium containing 10 mM glycerol and 40 mM succinate (permissive growth condition) was used as precultures for growth experiments. Strains were harvested (6,000 * g, 3 min, RT) and washed three times in M9 medium without carbon source. Cultures were inoculated into the M9 media to an OD 600 of 0.01 in a 96-well microtiter plate (Nunclon Delta Surface, Thermo Scientific). Each well contained 150 μL of culture and 50 μL mineral oil (Sigma-Aldrich) to avoid evaporation. Growth monitoring and incubation at 37°C was carried out in a microplate reader (EPOCH 2, Bio-Tek). In the program, 4 shaking phases of 60 seconds were repeated three times (linear shaking 567 cpm (3 mm), orbital shaking 282 cpm (3 mm), linear shaking 731 cpm (2 mm), and orbital shaking 365 cpm (2 mm)). After the shaking cycles, absorbance at 600 nm was measured. Raw data were converted to 1 cm wide standard cuvette OD values according to OD cuvette = OD plate /0:23. Matlab was used to calculate growth parameters. All experiments were carried out in at least three replicates. Average values were used to generate the growth curves. Variability between triplicate measurements was less than 5% in all cases displayed.
2.8. 13 C Isotopic Labelling Experiments. 13 C-isotope tracing was performed to deduce carbon flux. As described previously [25], 2 mL of an early stationary culture, grown in M9 minimal medium containing 1,6-13 C 2 -glucose (Sigma-Aldrich, Taufkirchen, Germany) as sole carbon sources, was pelleted and washed in ddH 2 O. The cell pellet was hydrolyzed in 1 mL 6 N HCl at 95°C for 24 h. Subsequently, the HCl was removed by evaporation at 95°C under an airstream. The hydrolyzed biomass was resuspended in 1 mL ddH 2 O. Amino acid masses were analyzed after separation by ultraperformance liquid chromatography (Acquity, Waters, Milford, MA, USA) using a C18-reversed-phase column (Waters, Eschborn, Germany) as described previously [25,29]. Mass spectra were acquired by an Exactive mass spectrometer (Thermo Scientific, Dreieich, Germany). Data was analyzed using Xcalibur (Thermo Scientific, Dreieich, Germany). Amino acid standards (Merck, Darmstadt, Germany) were used to determine specific retention times.
2.9. RNA Extraction and cDNA Synthesis. Biological triplicates were cultured on M9 minimal medium with either 20 mM glycerol (WT and Δtpi, NGR2 iso1) or 4 mM glycerol and 40 mM succinate (WT, Δtpi, Δpgk) as carbon source. In exponential phase,~2:5 × 10 8 cells were harvested, mixed with two volumes of RNAprotect™ Bacteria Reagent (Qiagen, Hilden, Germany), pelleted down, and stored at -20°C. Total RNA was extracted using the RNeasy Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions. In brief, enzymatic lysis was followed by oncolumn removal of genomic DNA with RNase-free DNase™ (Qiagen, Hilden, Germany) and spin-column-based purification. Integrity and concentration of the isolated RNAs were confirmed by gel electrophoresis and NanoDrop™ One (Supplementary Table S5 A), respectively. Reverse transcription was performed on 500 ng RNA with the qScript™ cDNA Synthesis Kit (QuantaBio, Beverly, MA USA). A RT-negative control has been produced for each sample as reported in the Supplementary Table S5, B. The resulting nucleic acid stock concentration was 25 ng/μL and it was stored at -20°C.

Expression Analysis by Reverse Transcriptase
Quantitative PCR (RT-qPCR). Primers for qPCR have been designed with Eurofins Genomics' qPCR Assay Design Tool (https://www.eurofinsgenomics.eu/en/ecom/tools/qpcr-assaydesign/), based on Prime+ of the GCG Wisconsin Package enhanced with additional parameters for perfect probe design. Primers used in this experiment are shown in Supplementary  Table S3. The gene encoding 16S rRNA has been widely used as a bacterial reference gene [30]. It is expressed from rrsA, rrsB, rrsC, rrsD, rrsE, rrsG, and rrsH operons [31], and we chose it as well-established transcript for expression normalization of cells in exponential growth. Primer specificity was assessed by conventional PCR on template cDNA using DreamTaq Green™ polymerase (Thermo Fisher Scientific, Dreieich, Germany) ( Figure S6). 125 pg of nucleic acid (0.5 μL 1 : 100 stock solution), 1 μM primers (2μL 2,5 μM primer mix), and 2.5 μL Maxima SYBR™ Green/ ROX qPCR Master Mix (Thermo Scientific, Dreieich, Germany) were used per reaction. Three technical replicates were analyzed for each biological replicate. The experiment was carried out with the Applied Biosystems 7900HT Fast Real-Time PCR System (Applied Biosystems) in clear MicroAmp™ Optical 384-well reaction plates with Barcode (Applied Biosystems). Each plate included the respective RNA (RT-) and no-template (NTC) negative controls. Differences in expression levels were calculated according to the 2 -ΔΔCT method [32][33][34]. Reported data represents the 2 -ΔΔCT value that was calculated for each sample individually relative to the average of all biological WT replicate ΔCt (Ct(GOI)-Ct(rrsA)) values. RT-controls showed that residual gDNA contributed to an average of less than 3% of the signal for mgsA expression analysis, with a maximum of 7.3%. RT-controls for serA, serB, serC, and sdaA expression analysis indicate a wider gDNA contribution to the signal across samples, with an average 8.5%, with one maximum of 30%. Data shown are not corrected for this gDNA contamination, but gDNA signal is listed in Supplementary Material qPCRs Analysis + RT-controls. Primer efficiency 5 BioDesign Research was assessed by qPCR of 6 consecutive tenfold dilutions of 25 ng/μL using wildtype genomic DNA. Efficiency was calculated as E = ð10 ð−1/slopeÞ − 1Þ × 100. The linear range of primers was globally between Cq 9.04 and 27.72. Linear regression showed that primers have approximate efficiencies between 89 and 103.5%, with R 2 between 0.9988 and 1. The efficiencies of all primers were within the 90-105% range, with the exception of serB (89%).
2.11. Protein Expression and Purification. The His-taged WT and mutated SerA proteins were expressed in E. coli BL21 (DE3) Codon+ (Invitrogen). The cells were grown in 400 mL terrific broth containing 100 μg/mL carbenicillin at 37°C until they reached an OD 600nm = 2 upon which expression for 16 h at 20°C was induced by addition 500 μM IPTG. The cells were harvested by centrifugation for 30 min at 10000g at 4°C. Cell pellets were frozen at -80°C for one night. Thawed cells were then suspended in 32 mL of Buffer A (50 mM phosphate (Na/K), 500 mM NaCl, 30 mM imidazole, 15% glycerol, and pH 8.0) and lysed for 30 min at room temperature after addition of 3.6 mL of Bug Buster (Novagen), 32 μL DTT (dithiothreitol) 1 M, 320 μL Pefabloc 0.1 M (Millipore), and 23 μL Lysonase (Novagen). Lysate was clarified at 9000g for 45 min at 4°C and then loaded onto a 5 mL HisTrap FF column preequilibrated in buffer A. The protein was eluted in buffer B (50 mM phosphate (Na/K), 500 mM NaCl, 250 mM imidazole, 1 mM DTT, 15% glycerol, at pH 8.0) and desalted on a gelfiltration column Hi Load 16/60 Superdex 200 pg in buffer C (50 mM Tris, 50 mM NaCl, glycerol 15%, 1 mM DTT, at pH 8.0). The protein was frozen and stored at -80°C if not immediately used for assays.

Characterization of SerA Kinetic Parameters.
Assays were performed using a Safas UV mc2 double beam spectrophotometer at room temperature using quartz cuvettes (0.6 cm path length). Assays of SerA-catalyzed reduction of 2-ketoglutarate were conducted in 40 mM potassium phosphate, 1 mM DTT, at pH 7.5. Kinetic parameters for 2-ketoglutarate were determined by varying its concentration (from 0.005 to 1 mM) in the presence of a saturating concentration of NADH (250 μM). The reactions were monitored by recording the disappearance of NADH at 340 nm (molar extinction coefficient = 6220 M −1 :cm −1 ). Kinetic constants were determined by nonlinear analysis of initial rates from duplicate experiments using SigmaPlot 9.0 (Systat Software, Inc.).

Characterization of Serine Inhibition of SerA 2-
Ketoglutarate Reductase Activity. Standard enzyme assays were conducted with both substrates at saturating concentrations (250 μM NADH and 600 μM 2-ketoglutarate) in 120 μL final volume. Reaction mixes contained 5 μg of purified enzyme. Serine was added at concentrations ranging from 0.5 to 10 μM.

In Silico Analysis of Potential Pathways Bypassing
Glycolysis in E. coli. We aimed to uncover latent metabolic routes that can potentially replace the canonical EMP glycolysis and offer new connections between phosphosugar and organic acid metabolism. We used the latest metabolic model of E. coli from the BiGG database [35] and systematically searched for all thermodynamically feasible combinations of native reactions that can convert the feedstock glycerol into pyruvate, an organic acid from which all metabolites of the TCA cycle can be derived (Methods). We chose glycerol as the feedstock as it can be directly converted into the simplest phosphosugars-the triose phosphates dihydroxyacetone phosphate (DHAP) and glyceraldehyde 3-phosphate (GAP)-such that its conversion to pyruvate is expected to follow a 1 : 1 stoichiometry. Using this approach, we were able to identify more than 100 thermodynamically feasible routes that can potentially convert glycerol to pyruvate (https://gitlab .com/elad.noor/glycolysis-bypass/-/blob/master/results/ figureS2.pdf). These routes and their variants, e.g., routes that share the same intermediates and general metabolic conversions but use different cofactors (e.g., quinones instead of NAD + ), represent combinations of four broad strategies to convert glycerol to pyruvate ( Figure 1): (i) oxidation of triose phosphates to phosphoglycerate, followed by conversion to phosphoenolpyruvate and pyruvate, as within canonical EMP glycolysis; (ii) generation of pyruvate from phosphosugars via the Entner-Doudoroff (ED) pathway (blue arrows in Figure 1); (iii) conversion of triose phosphates to methylglyoxal, which is then oxidized to pyruvate (magenta arrows in Figure 1); and (iv) conversion of 3-phosphoglycerate (3PG) towards serine biosynthesis, followed by serine deamination to pyruvate (green arrows in Figure 1). Since the ED pathway has been analyzed and compared to EMP glycolysis in multiple previous studies [7,[36][37][38], we decided to focus on the methylglyoxal-dependent and serine-dependent routes and test their potential of bypassing the EMP in E. coli. Notably, our approach differs from a previous study which used a similar computational method but included all KEGG database reported reactions (not only from E. coli), limiting their search to pathways yielding at least 1 ATP [8]. This method would not find the routes, identified by our approach, via serine or methylglyoxal which are ATP-neutral or even waste ATP, respectively.

A Glycolytic Pathway via
Methylglyoxal. E. coli is known to channel flux towards methylglyoxal upon accumulation of triose phosphates which results from phosphate depletion limiting the activity of GAP dehydrogenase or from excessive carbon intake [17,39]. The operation of the methylglyoxal bypass thus enables the cells to adapt to imbalanced central metabolism [40]. Considering the high toxicity of methylglyoxal [41,42], it was not clear whether the methylglyoxal bypass could indeed support all glycolytic flux without resulting in severe growth defects.
To assess the capability of methylglyoxal metabolism to replace glycolysis, we first constructed a Δtpi strain. This strain, when fed with glycerol as sole carbon source, requires the activity of the methylglyoxal shunt to support almost the entire carbon assimilation flux (different to when using glucose as feedstock). To our surprise, when cultivated on glycerol as sole carbon source, the Δtpi strain immediately grew 6 BioDesign Research without the need for any dedicated overexpression of methylglyoxal synthase (MgsA) or evolution (yellow line in Figure 2(a)). As previously described, a Δtpi strain utilizes glucose-derived GAP via EMP and DHAP via the methylglyoxal route [43]. On glycerol, the tpi deletion can theoretically be bypassed by the ED pathway. Here, a GAP is "borrowed" to condense with a DHAP to make fructose 1,6-bisphosphate via fructose 1,6-bisphosphate aldolase (Fba) and then use the ED pathway to generate GAP and pyruvate, which both can recover the GAP condensed with DHAP by Fba. In order to exclude this theoretical option, we deleted the genes zwf and eda individually in the Δtpi strain to block the ED pathway. The resulting strains Δtpi Δzwf and Δtpi Δeda grew similar to the Δtpi strain, ensuring no contribution of the ED pathway to the metabolic bypass in the tpi deletion strain. An even stronger confirmation of the activity of the methylglyoxal pathway is provided by the fact that the deletion of mgsA in the Δtpi strain abolished growth (Figure 2(b)). However, one alternative explanation is that the deletion of mgsA causes a DHAP accumulation which inhibits growth. We aimed at excluding that this secondary effect, rather than the absence of the methylglyoxal bypass, is responsible for the inability of the Δtpi ΔmgsA strain to grow on glycerol. Therefore, we analyzed growth of the Δtpi ΔmgsA strain on medium additionally containing succinate. Here, it grew similar to the Δtpi strain; thus, toxic effects caused by the absence of mgsA, e.g., by DHAP accumulation, can be excluded. We then constructed a strain lacking EMP glycolysis by deleting phosphoglycerate kinase (Δpgk). As expected, the Δpgk strain was only able to grow if two carbon sources were provided, feeding both parts of the metabolism (divided by the pgk deletion), e.g., glycerol and succinate. As the ED pathway could theoretically bypass the pgk deletion, we additionally blocked it by deleting the 2-keto-3-deoxygluconate 6phosphate aldolase gene (eda), resulting in strain Δpgk Δeda. In contrast to the Δtpi strain, the Δpgk Δeda strain required overexpression of the gene coding for methylglyoxal synthase (mgsA) from a plasmid to grow on glycerol as sole carbon source (green vs. blue line in Figure 2(a)).
Our results indicate that the methylglyoxal bypass can indeed completely replace glycolysis, despite the high reactivity of its namesake intermediate. Notably, when glycerol is the carbon source, DHAP is expected to accumulate to a high concentration in the Δtpi strain, where its metabolism is completely blocked, while in the Δpgk Δeda strain such accumulation can be prevented by DHAP conversion to other phosphosugars ( Figure 1). As it has been reported that accumulation of DHAP leads to the expression of mgsA and activation of the methylglyoxal bypass in E. coli   Figure S2). These results suggest that growth of the Δtpi strain fed on glycerol is enabled by an elevated cellular DHAP concentration which enforces a high carbon flux through the methylglyoxal pathway without requiring any change of expression of mgsA.
In the case of the Δtpi strain fed on glycerol, pyruvate as the final product of the methylglyoxal pathway is needed to supply all parts of central carbon metabolism through gluconeogenesis, TCA cycle, and anaplerotic reactions. Contrarily, in the case of the Δpgk Δeda strain overexpressing mgsA, pyruvate should be essential for the operation of TCA and anaplerotic reactions and for providing 3-phosphoglycerate as precursor of serine and glycine biosynthesis. Thus, in both strains, a deletion of phosphoenolpyruvate synthase (ppsA), which is essential to provide PEP from pyruvate, e.g., anaplerotic reactions, should impede methylglyoxal-dependent growth ( Figure 1). Indeed, the Δtpi ΔppsA strain and the Δpgk Δeda ΔppsA (overexpressing mgsA) strain lost the ability to grow on glycerol (Figures 2(b) and 2(c)), thus confirming the pivotal activity of the methylglyoxal bypass.
Methylglyoxal can be converted to pyruvate via three routes (Figure 1, magenta arrows): (i) methylglyoxal attachment to glutathione, followed by hydrolysis to D-lactate and oxidation to pyruvate; (ii) direct electron rearrangement of methylglyoxal to give D-lactate, which is oxidized to pyruvate; (iii) reduction of methylglyoxal to lactaldehyde, followed by oxidation to D-lactate and pyruvate. Previously, the glutathione-dependent pathway was suggested to be the dominant route [19,40,44,45]. However, an unequivocal proof for this claim is missing. To finally settle this question, we used the Δtpi and Δpgk Δeda strains to generate multiple deletion strains, each carrying the deletion of a different enzyme of the methylglyoxal catabolism. We found that the absence of HchA, GldA, or AldA did not impair growth on glycerol via the methylglyoxal bypass (Figure 2(b)), indicating that the two glutathione-independent routes described above contribute very little, if at all, to methylglyoxal metabolism. On the other hand, the absence of GloA or Dld completely abolished growth on glycerol (Figure 2(b)). These findings confirm the glutathione-dependent route to be indispensable and sufficient for methylglyoxal metabolism. In many cases, 13 C-labeling experiments are valuable to verify a pathway's activity, while in the case of the methylglyoxal route, no difference in the labeling pattern compared to EMP glycolysis was expected (Supplementary Fig. S3a); a 13 C-labeling experiment in the Δtpi was carried out using 1,6-13 C 2 -glucose. The Δtpi strain was used because other strains were not able to grow on glucose. 1,6-13 C 2 -glucose is converted via upper glycolysis to 3-13 C GAP and 3-13 C DHAP. In the Δtpi strain, GAP is converted to pyruvate via EMP glycolysis and DHAP is catabolized via the methylglyoxal route. As can be seen in Supplementary Fig. S3b, the labeling pattern was similar to the WT and matches the predictions, showing that the 13 C stays in the C-3 position; thus, no other pathway which involves a carbon rearrangement is active converting the isolated DHAP to pyruvate.

The Serine Shunt
Can Replace EMP Glycolysis. The second glycolytic bypass identified by our in silico analysis is the serine shunt which combines biosynthesis and degradation
To assess the feasibility of this pathway, we generated an enolase deletion strain (Δeno), which cannot operate EMP glycolysis while still being able to generate 3phosphoglycerate from which the serine biosynthesis pathway begins (Figure 1). For growth on a minimal medium, the Δeno strain requires two carbon sources for feeding both upper and lower parts of central carbon metabolism (e.g., glycerol and succinate or xylose and succinate). We assumed that overexpression of the four enzymes involved in serine biosynthesis from 3-phosphoglycerate and degradation to pyruvate would enable the activity of the serine shunt and support growth of the Δeno strain on glycerol only (Figure 1, green arrows). Therefore, we constructed a plasmid constitutively expressing a synthetic operon containing the four following genes: (i) a variant of the native gene coding for 3-phosphoglycerate dehydrogenase which was engineered to remove its native allosteric inhibition by serine (ser A * , coding for SerA H344A N346A N364A, catalyzing the first step in serine biosynthesis) [46]; (ii) serC, coding for phosphoserine aminotransferase; (iii) serB, coding for phosphoserine phosphatase; and (iv) sdaA, coding for the major serine deaminase isozyme in E. coli, converting serine into pyruvate [47]. However, the resulting plasmid (p-ser A * -serB-serC-sdaA) failed to support growth of the Δeno strain on glycerol. Moreover, PCR analysis of the transformants revealed that the strains did not carry the whole operon but only fractions of it. We assume that this occurred due to toxic effects of the overexpression of the serine biosynthesis genes and hence resulted in gene inactivation/removal from the plasmid.
Serine's toxicity is well known and can be attributed, at least partially, to its deamination to the highly reactive and toxic compound hydroxypyruvate [48]. Indeed, we found that addition of even small amounts of serine (1-8 mM) substantially delayed the growth of the WT (Figure 3(a)) and the Δeno strain (Figure 3(b)) with glycerol and succinate as co-carbon sources. The Δeno strain seems to be more sensitive to serine than the WT strain, as its growth was completely inhibited at serine concentrations of~8 mM. Further deletion of the genes coding for serine deaminases (ΔsdaA ΔsdaB ΔtdcB ΔtdcG [49]), thus removing serine sinks, increased the sensitivity of both the WT and the Δeno strain further: even small amounts of serine (0.5 mM) completely inhibited growth on glycerol and succinate (Figures 3(c) and 3(d)). Interestingly, in one of the Δeno strain cultures incubated in the presence of 48.8 mM serine, the cells began to grow after >100 hours (Figure 3(b)). This suggests the emergence of mutations enabling the cells to tolerate this high level of serine. We isolated three individual serine-tolerant strains (G3 mutants) and sequenced their genomes. In the resulting reads, we found sequence differences in pcnB encoding poly(A) polymerase I and leuB encoding 3-isopropylmalate dehydrogenase of leucine biosynthesis, both not directly connectable to the serine insensitivity of the mutants. But, in all three isolated strains, a somewhat different genomic region was multiplied several fold (as indicated by a 4-8 fold increased sequencing cover-age, Supplementary Table S1 and Supplementary Fig. S4) (G3  #2). Notably, these amplified genomic regions harbor the genes sufA and sufB, coding for subunits of the iron-sulfur cluster scaffold complex [50]. Among its other functions, this cluster serves as an essential component of the primary serine deaminase enzyme SdaA [51]. Hence, it seems likely that increasing the availability of iron-sulfur clusters in the cell contributed to a faster and more efficient degradation of serine, making the cells more tolerant to this amino acid. To test if the increased tolerance to serine would enable the mutated strain to grow on glycerol via the serine shunt, we overexpressed the four genes described above. Indeed, we found that transformation of the G3 mutant strains with p-ser A * -serB-serC-sdaA enabled them to grow on glycerol as sole carbon source (Figure 4(a)). It therefore seems that the above-mentioned changes supported a strong metabolic sink for toxic serine to enable high flux via serine biosynthesis and degradation with minimal adverse effects.
In order to decrease and stabilize the rate of serine biosynthesis to avoid its accumulation, we pursued a strategy of gene expression from the chromosome rather than from a plasmid in the naïve Δeno strain. The genes ser A * and sdaA were overexpressed from the chromosome (see method section), while the native expression was preserved for the genes serB and serC. We found that overexpression of ser A * from the chromosome in the Δeno strain is possible only after the chromosomal overexpression of sdaA is established, again pointing to the necessity of a balanced serine synthesis and degradation activity for growth via the serine shunt (data not shown). Finally, to avoid any carbon flux through the methylglyoxal bypass and the ED pathway, we deleted mgsA and the gene encoding glucose 6-phosphate dehydrogenase (Δzwf). The Δeno ΔmgsA Δzwf strain overexpressing ser A * and sdaA from the chromosome was indeed able to grow on glycerol with no need for prior adaptation (Figure 4(b)). As expected, the deletion of ppsA, which blocks the conversion of pyruvate to phosphoenolpyruvate, abolished growth (Figure 4(b)), thus confirming that this essential metabolite cannot be produced by some variant of EMP glycolysis that bypasses or replaces the enolase reaction. Overall, our results indicate that the serine shunt can be implemented in the cells to replace the canonical EMP glycolysis using a rational engineering approach, where the rate of serine biosynthesis is balanced with that of serine degradation, consequently avoiding deleterious accumulation of this amino acid.

Continuous Culture Evolution
Results in the Emergence of the Serine Shunt. As an alternative option to the rational engineering of a synthetic glycolysis bypass, we resorted to continuous culture protocols to investigate which of the three bypass routes-ED pathway, methylglyoxal bypass, or serine shunt-would emerge naturally during long-term cultivation of the Δeno strain under selective conditions (we chose to start with the Δeno strain as it has the potential to activate all three mentioned bypass routes). We applied a medium-swap regime 9 BioDesign Research using GM3 continuous cultivation devices [26,52] to a growing population of the Δeno strain, fed alternatively with a permissive medium containing both glycerol and succinate, and a nonpermissive medium containing only glycerol (Methods). When the turbidity of the culture, measured in real time, was below a predefined value, the culture was diluted with the permissive medium; otherwise, the nonpermissive medium was used to dilute the culture. Such a medium-swap protocol enables gradual adaptation of the bacterial population to conditions it initially cannot grow in [26,52], in our case proliferation without succinate as substrate of the lower part of the carbon metabolism.
We conducted the adaptive evolution experiment in two parallel cultures, both of which adapted to grow on glycerol as sole carbon source after 130-140 days amounting to more than 900 generations (Figure 5(a)). We then  10 BioDesign Research applied a turbidostat mode-cultivating the cells solely on glycerol and diluting the culture every time a predefined turbidity is reached-to evolve the culture towards a higher growth rate ( Figure 5(a)). We isolated two strains from each of the evolved cultures and sequenced their genomes. Strains isolated from each of the evolved cultures showed highly similar mutation profiles (Supplementary  Table S2). All sequenced strains harbored a nonsynonymous mutation in serA, either replacing T372 with asparagine or replacing L370 with methionine. No mutation was observed in genes coding for enzymes participating in the methylglyoxal shunt or the ED pathway. This suggested that the adaptive evolution of the Δeno strain led to the emergence of the serine shunt, rather than to any other of the possible glycolytic bypass routes. However, neither known bypasses (ED and MG) nor other latent bypasses can be excluded at this point. We found that the isolates from only one of the evolved cultures could stably grow on glycerol alone when cultivated within a 96-well plate ( Figure 5(b)). The strain iso1 was further characterized.
To exclude that an unknown pathway which involves carbon rearrangement is active in the iso1 strain, 13 C-labeling experiments were performed. The iso1 strain was grown in 1,6-13 C 2 -glucose which is converted via upper glycolysis to 2 molecules of 3-13 C-GAP which then is converted to 3-13 Cphosphoglycerate. Regardless which route is taken, EMP, ED, methylglyoxal, or serine shunt, the labeling pattern in amino acids derived from respective metabolites should be identical (Supplementary Fig. S3) unless a pathway involving

BioDesign Research
carbon rearrangement is active. Indeed, the labeling pattern of the iso1 strain was identical to the labeling of the WT control in all amino acids analyzed ( Supplementary Fig. S3). Thus, carbon rearrangement can be excluded.
The four reactions of the serine shunt usually carry much less flux than needed for the serine shunt to operate as sole glycolytic option. To analyze if the expression of the corresponding genes is altered in the iso1 strain, we performed qPCR experiments after growth on glycerol. The transcript level of serC and serB genes was increased in the iso1 strain compared to the WT reference, showing a 11-and 2.5-fold increased RNA level of serC and serB, respectively. Surprisingly, serA and sdaA expressions did not significantly differ from the WT control (Supplementary Fig. S5). To ensure that no glycerol is converted via the alternative bypasses, we deleted their key reactions in the iso1 strain. Neither deletion of mgsA nor zwf abolished growth (Figures 6(a)-6(c)), indicating that the glycolysis bypass in the iso1 strain does not rely on methylglyoxal bypass or the ED pathway. Despite our continuous efforts, we were unable to delete sdaA, which suggests that serine deamination is essential for growth or survival also in rich medium (e.g., Medium X, see Methods), pointing again to the serine shunt as the bypassing route.
To further analyze the operation of the serine shunt in the evolved strain, we restored the WT version of the gene serA in strain iso1 and performed growth tests on glycerol. The recombinant SerA WT derivatives of iso1 had lost the ability to grow on glycerol, but the growth phenotype appeared very leaky. Interestingly, we observed contrasting growth behaviors between various experimental replicates ( Figure 6(d)): some were unable to use glycerol as sole carbon source, and others started to grow slowly after a short lag-phase, while a few grew faster after a similar delay. These observations strongly indicated the emergence of mutations. Indeed, Sanger sequencing of PCR-products of the serA locus of every isolate from a growing culture revealed that SerA was mutated. Each of the tested strains contained one of the following mutations: A367T, H342Y, L332P, or R338C, while none of them harbored the SerA variants L370M or T372N identified in isolates from the long-term evolution experiments. Since these mutants arose rather quickly, within hours or days, we concluded that they already appeared in the precultures, which contained glycerol and succinate. Replacing glycerol with xylose for alimenting the upper part of central metabolism in the preculture prevented the emergence of mutants growing on glycerol as sole carbon source within the duration of the experiment (>5 days) (Figure 6(e)).
Finally, we deleted the serA gene in the iso1 strain. A serA deletion blocks serine biosynthesis and, hence, turns the iso1 strain into a serine auxotroph. This strain allowed us to test the serine demand of the strain and compare it to a WT deleted in serA (ΔserA). If a different pathway than the route via serine biosynthesis would be used, the serine demand of both strains, iso1 ΔserA and ΔserA, would be similar. But, if growth of the glycolytic bypass active in the iso1 strain is dependent on the serine shunt, the supplemented serine would account for all biomass produced from lower metabolism, compared to serine derived biomass building blocks in the ΔserA strain only. As expected, both strains, iso1 ΔserA and ΔserA, did not grow without serine supplementation, and, when a gradient of serine was supplemented to the medium containing 20 mM glycerol, the iso1 ΔserA strain required much higher serine concentrations to reach similar ODs compared to the ΔserA strain ( Figure 6(f), Supplementary Fig. 7). These results indicate that the majority of the serine provided to the iso1 ΔserA strain is used for providing carbon to feed lower metabolism, e.g., TCA cycle intermediates.

12
BioDesign Research biosynthesis and degradation pathway as glycolytic bypass: e.g., inability to delete the serine deaminase gene sdaA in the evolved strains, mutations fixed in serA during the evolution experiment, and spontaneous mutations arising in serA in the evolved genetic context of iso1 strain where the WT allele of serA had been restored. Interestingly, all the mutations identified in SerA in the course of our study are located in the C-terminal allosteric regulatory domain of the enzyme and thus could possibly modulate the negative allosteric feedback of L-serine on SerA activity. To test whether this hypothesis was correct, we characterized both the kinetic parameters and the regulatory effect of L-serine on enzyme activity of the purified SerA WT, the SerA variants which emerged during the evolution experiments (SerA T372N, culture NGR1) and (SerA L370M, culture NGR2), as well as the triple mutation SerA variant (H344A N346A N364A), which was used for the rational engineering of the serine shunt and was previously reported to be feedback-resistant [46]. The promiscuous 2-ketoglutarate reductase activity of SerA was monitored, which, in contrast to the oxidation of 3-phosphoglycerate, is a thermodynamically favorable reaction and is known to be also regulated by L-serine [53]. All the variants tested exhibited 2-ketoglutarate reductase activities comparable with that of the WT enzyme (Table 2), which is in accordance with the fact that none of the mutations were close to the catalytic site of the enzyme [46]. As expected, the activity of WT SerA was strongly decreased in the presence of micromolar concentrations of L-serine (0.5-10 μM) ( Table 2). This concentration range corresponds to the previously published IC 50 for L-serine, which was determined to be between 2 and 10 μM [53]. By contrast, the SerA variants showed a significantly lower sensitivity to the presence of L-serine, as demonstrated by the minimal decrease of the activity of L370M, T372N, and H344A N346A N364A SerA variants with increasing L-serine concentration (Table 2). Our findings in both, the isolated mutants as well as in the rationally engineered strains, support that the operation of the serine shunt is highly dependent on the presence of a feedback-resistant SerA variant. Kinetic properties of the downstream serine deaminase SdaA (K M of~2.7 mM for L-serine) support this conclusion. The K M is three orders of magnitude higher than the IC 50 of SerA for L-serine [51]; thus, a feedbackinhibited SerA variant might not allow sufficient flux via the pathway to achieve L-serine concentrations high enough for SdaA activity.

Discussion
Our computational analysis revealed the presence of multiple feasible glycolytic bypasses in E. coli's native metabolic network. All of these were combinations of the canonical catabolic EMP and ED pathways and two cryptic routes, which so far had not been described as a significant bypass to EMP glycolysis in any organism (in the case of the serine shunt) or were only described as a sink for excess carbon (i.e., the methylglyoxal route) [18,19,43]. Here, we were able to engineer these two glycolytic alternatives in E. coli, relying exclusively on native enzymes. Previous studies observed the emergence of the methylglyoxal pathway in a Δtpi strain during growth on glucose, thus splitting carbohydrate metabolism between EMP glycolysis from GAP and methylglyoxal metabolism from DHAP [18], and the serine shunt was proposed to improve acetyl-CoA precursor supply for the bioproduction of poly(3-hydroxybutyrate) [54]. But as of yet, to the best of our knowledge, both pathways have not been described as a complete glycolytic bypass. In addition to the rational engineering of the pathways, a directed evolution experiment evolving an enolase deletion strain was conducted, leaving all bypass options open for natural selection. This experiment resulted in the emergence of the serine shunt. Intuitively, a bypass via the ED pathway, which operates in E. coli when growing on gluconate [55], might seem more likely to occur. But, evolving the ED pathway to bypass the glycolytic blockade presumably requires the accumulation of several adaptive mutations susceptible to increase the gluconeogenic flux, to inactivate/downregulate the 6-phosphogluconate dehydrogenase and to enhance phosphogluconate dehydratase (edd) and 2-keto-3-deoxygluconate 6-phosphate aldolase (eda) activity. In contrast, and in line with the outcome of the evolution experiment, engineering the functional serine shunt required only the balanced overexpression of a feedback-resistant SerA variant and of SdaA. Curiously, we showed that the methylglyoxal pathway needed only the overexpression of MgsA in the Δpgk strain; therefore, it was expected that the methylglyoxal pathway could be established most easily as a glycolytic bypass by evolution. One reason for the selection of the serine shunt instead of the methylglyoxal pathway during evolution of the enolase deletion strain might be the higher intracellular concentration of 3-phosphoglycerate (~4 mM) compared to DHAP (~0.5 mM) during growth on glycerol [56]. The 3-phosphoglycerate concentration, determined for WT E. coli, is likely even higher in the Δeno strain; this together with the lower ATP cost of the serine shunt might have triggered the fixation of mutations establishing the serine shunt as Δeno bypass. Noteworthy, in a previous study, a pyruvate auxotrophic strain was evolved for the emergence of pyruvate generating reactions; here, the underlying pathway remained unclear but computational analysis revealed several pathway options, among them also the serine and the methylglyoxal option [57].
As expected, reverting the serA mutation which emerged in the evolved Δeno strain abolished growth on glycerol as sole carbon source. However, after a short period of incubation, growth associated with the appearance of various point mutations in the serA gene was restored. This indicates that the strain underwent adaptation processes during the evolution, possibly priming its metabolism for the use of the serine shunt. These can also be seen in the qPCR results of the involved genes. However, these adaptations are not directly reflected in any of the mutations identified in the genome sequencing. A further indication of the activity of the serine shunt was provided by the strain's high demand of serine supplementation to the medium, when serA was deleted in the evolved strain (iso1 ΔserA). When compared to a WT ΔserA reference, significantly higher serine concentrations were needed by iso1 ΔserA in order to reach similar OD (600 nm). All together our experiments provide strong evidence for the activity of the serine shunt in iso1.
Both of the glycolytic bypasses contain toxic intermediates. In our experiments, only the toxic effects of L-serine caused a stress response in the strain transformed with the plasmid overexpressing serine biosynthesis and degradation genes. On the other hand, no additional overexpression of methylglyoxal degradation enzymes was necessary to implement the methylglyoxal bypass, indicating the sufficiently fast native glutathione-dependent conversion of the highly reactive methylglyoxal to D-lactate. Besides involving a toxic intermediate, the serine shunt presents a thermodynamic barrier at the level of 3-phosphoglycerate dehydrogenase, which additionally is highly inhibited by low concentrations of L-serine.
While many reactions of central metabolism can be catalyzed using different electron acceptors, e.g., PQQdependent glucose oxidation [58], or cosubstrates, e.g., PPi instead of ATP for phosphorylations by phosphofructokinase or pyruvate-phosphate dikinase [59], the general structure of central metabolism remains unchanged. Why some glycolytic pathways and their variants emerged and manifested throughout all organisms to convert sugars into biomass building blocks, while other biochemical possibilities either never appeared in nature or were discarded during evolution, is not well understood. Some insights came from previous studies of the connection between sugar catabolism and the organism's environment (i.e., anaerobic or aerobic). There is evidence that ATP yield outweighs protein cost of a pathway when operating under anaerobic conditions; hence, the EMP pathway is dominantly used under these conditions. Contrarily, obligate aerobic organisms prefer the ED pathway, which ensures a lower protein cost and a higher thermodynamic driving force. Accordingly, facultative anaerobic organisms like E. coli possess both pathway options [7]. The methylglyoxal pathway and serine shunt, besides their structural differences, significantly differ in their ATP balance, i.e., consumption of ATP (methylglyoxal pathway) versus no consumption/formation of ATP (serine shunt) per pyruvate generated. These differences in ATP yield render both pathways infeasible in the absence of terminal electron acceptors, e.g., under fermentative conditions. This might provide an explanation for the absence of these pathways in nature. However, the nonphosphorylative ED pathway operates as the route for glucose to pyruvate conversion in some archaea and hence similar to the serine shunt produces no net ATP yield [60].
In the past decade, important breakthroughs were achieved in the field of synthetic metabolism. The implementation of C1 assimilatory routes for the fixation of CO 2 [61,62], formate [63], or formaldehyde [64], enabling biomass production from renewable resources, are examples of audacious attempts to rewire central carbon metabolism. While primarily motivated by the need to provide sustainable options for the industrial production of fuels and chemicals, these studies accelerate the understanding of the underlying principles of metabolism. Other examples like the implementation of phosphoketolase-dependent nonoxidative glycolysis [12,13] illustrate the pace synthetic biology has adopted to define the limits and barriers which might have shaped the structures of central metabolism at the onset of cellular life. Our approach, combining rational design of nonnatural pathways and experimental evolution, is in line with this effort trying to "make sense of the biochemical logic of metabolism" [6] and at the same time opens opportunities for new bioproduction routes of potential industrial relevance, e.g., for optimizing production of serine, cysteine, or 1,2-propanediol.

Data Availability
The data supporting the presented findings are available within the paper or in its Supplementary files. Strains and plasmids used here are available on request from the corresponding author. Data from the following public repositories were used in this study: BiGG (http://bigg.ucsd.edu/), KEGG (https://www.kegg.jp/), and eQuilibrator (http:// equilibrator.weizmann.ac.il/).

Additional Points
Code Availability. The code used in this study can be found at Zenodo (https://zenodo.org/record/6509926).

Conflicts of Interest
The authors declare no competing interests.