Latent Pathway Activation and Increased Pathway Capacity Enable Escherichia coli Adaptation to Loss of Key Metabolic Enzymes*

The ability of biological systems to adapt to genetic and environmental perturbations is a fundamental but poorly understood process at the molecular level. By quantifying metabolic fluxes and global mRNA abundance, we investigated the genetic and metabolic mechanisms that underlie adaptive evolution of four metabolic gene deletion mutants of Escherichia coli (Δpgi, Δppc, Δpta, and Δtpi) in parallel evolution experiments of each mutant. The initial response to the gene deletions was flux rerouting through local bypass reactions or normally latent pathways. The principal effect of evolution was improved capacity of already active pathways, whereas new flux distributions were not observed. Combinatorial changes in capacity and pathway activation, however, led to different intracellular flux states that enabled evolution in three of the four parallel cases tested. The molecular bases of the evolved phenotypes were then elucidated by global mRNA transcript analyses. Activation of latent pathways and flux changes in the tricarboxylic acid cycle were found to correlate well with molecular changes at the transcriptional level. Flux alterations in other central metabolic pathways, in contrast, were apparently not connected to changes in the transcriptional network. These results give new insight into the dynamics of the evolutionary process by demonstrating the flexibility of the metabolic network of E. coli to compensate for genetic perturbations and the utility of combining multiple high throughput data sets to differentiate between causal and noncausal mechanistic changes.

All biological systems are capable of short term response to environmental changes and, on longer time scales, to evolutionary adaptation. Due to the high growth rates and large numbers of individuals, adaptive evolution of microbes under laboratory conditions rapidly leads to improved growth phenotypes. This principle is exploited for evolutionary engineering (1,2) and experimental testing of general evolutionary principles (3)(4)(5). Combining experimental evolution with targeted genetic perturbations then allows one to pose specific questions on robustness and adaptability of biological networks. In response to such externally introduced deletions, microorganisms can invoke a number of different strategies to adjust their functionality. For metabolic net-works, these strategies can involve a local bypass of the deleted reaction, complete redirection of flux, reassignment of enzymes to catalyze the deleted reaction, activation of silent genes, or activation of otherwise down-regulated pathways. Whereas all of these mechanisms could potentially cope with genetic perturbations, the exact mechanisms utilized during evolution are poorly understood. Thus, a central question in adaptation and evolution is to determine whether evolving cells refine their existing pathway usage or whether they invoke major metabolic changes such as the activation of latent pathways.
The molecular basis of such evolutionary processes is now experimentally traceable with the ability to rapidly improve microbial phenotypes using laboratory evolution (from weeks to a few months) (6 -8) coupled with the principal accessibility of the underlying causes through various "omics" methods or genome resequencing (9). As a particularly popular tool, simultaneous transcript level monitoring of all genes within the genome by DNA microarray technology was used to identify altered gene expression in evolved Escherichia coli and Saccharomyces cerevisiae strains (4,10,11). Altered expression levels, however, do not distinguish between cause and effect and thus cannot directly reveal mechanistic links between altered expression and phenotype. In particular, when considering evolution of metabolic functions, more direct information on intracellular flux rerouting would be necessary to reveal the molecular mechanisms that cause a given improved phenotype. Such in vivo reaction rates are accessible through methods of 13 Cbased metabolic flux analysis (12), which have been used successfully to identify functional flux states in various microbes (13)(14)(15)(16)(17). Potentially, flux data can fill the gap between the intrinsically noisy and indirect transcriptome, proteome, or metabolome data and the actual phenotype (18,19). Thus, the combination of transcript profiles and quantitative intracellular flux data can provide greater insight into biological processes at the molecular level by implicating gene expression changes to altered phenotypes through association with measured flux data.
Beyond maximizing growth rates of wild-type strains on "exotic" substrates through adaptive evolution (6,7,20,21), rapid recovery of high growth rates was demonstrated for metabolic gene deletion mutants of E. coli (3). Replicates of evolved mutants exhibited phenotypic characteristics that suggested the selection of different biochemical mechanisms during parallel evolution under identical conditions. The molecular bases, however, remained unknown because multiple flux scenarios could explain the improved growth phenotypes. Here we evolved four E. coli knock-out mutants affected in metabolic key branch points, phosphoglucose isomerase (pgi), phosphoenolpyruvate carboxylase (ppc), phosphate transacetylase (pta), and triose-phosphate isomerase (tpi), for several hundred generations under exponential growth conditions on glucose. Since these lesions might be bypassed by at least two different routes, we used metabolic flux and global gene expression analysis to identify the metabolic mechanisms responsible for the improved phenotypes.

EXPERIMENTAL PROCEDURES
Strain Construction-The ppc and tpi mutants of MG1655 were described previously (3), and the pgi and pta mutants were constructed by in-frame gene deletions via homologous recombination facilitated by the red recombinase system (22) starting with the E. coli wild-type MG1655 (ATCC, Manassas, VA). The plasmids pKD46, pKD13, and pCP20 were used to introduce the recombinase gene, homologously recombine with the target gene, and remove antibiotic resistance markers, respectively. Each knock-out was confirmed by PCR with genomic DNA.
Adaptive Evolution-Evolution of constructed deletion mutants of E. coli was conducted in 250 ml of M9 minimal medium supplemented with 2 g/liter of glucose in 500-ml Erlenmeyer flasks using magnetic stir bars for aeration at 37°C. M9 medium contained (per liter of deionized water) 0.8 g of NH 4 Cl, 0.5 g of NaCl, 7.5 g of Na 2 HPO 4 ⅐2H 2 O, and 3.0 g of KH 2 PO 4 . The following components were sterilized separately and then added (per liter final volume of medium): 2 ml of 1 M MgSO 4 , 1 ml of 0.1 M CaCl 2 , 0.3 ml of 1 mM filter-sterilized thiamine HCl, and 10 ml of a trace element solution containing (per liter) 1 g of FeCl 3 ⅐6H 2 O, 0.18 g of ZnSO 4 ⅐7H 2 O, 0.12 g of CuCl 2 ⅐2H 2 O, 0.12 g of MnSO 4 ⅐H 2 O, and 0.18 g of CoCl 2 ⅐6H 2 O. At the start of evolution, initial precultures of each mutant were grown overnight in LB medium before being transferred to minimal medium for adaptive evolution. Duplicate evolution experiments were started from the same parental deletion mutant. In evolution cultures, cells were grown overnight and allowed to reach midexponential growth with an optical density at 600 nm (A 600 ) below 0.5 before being diluted by passage into fresh medium. The dilution factor at each passage was adjusted daily to account for changes in growth rate. The optical density was typically at an A 600 Յ 2.4 ϫ 10 Ϫ6 . This process of batch growth and serial passage was conducted for 30 days for the pta mutants (ϳ800 generations), 45 days for the ppc mutants (ϳ750 generations), and 50 days for the pgi (ϳ800 generations) and tpi (ϳ600 generations) mutants, where the ppc and tpi evolution experiments were reported previously (3). This process of evolution resulted in eight evolved mutants with end points designated as ptaE1, ptaE2, ppcE1, ppcE2, pgiE1, pgiE2, tpiE1, and tpiE2. The number of generations was estimated on a daily basis by calculating the starting optical density of each batch culture and determining how many doublings occurred during batch growth until being passed into fresh medium. Cultures were evolved until a stable growth rate was achieved for more than 5 days. This process of serial passage maintained a state of prolonged exponential growth so that no culture entered stationary phase. Duplicate cultures were evolved concurrently under identical conditions. 13 C-Labeling Experiments-Frozen glycerol stock cultures were used to inoculate LB complex medium. After 8 h of incubation at 37°C and constant shaking, LB precultures were used to inoculate M9 medium precultures that were grown overnight for inoculation of cultures for physiological or 13 C-labeling experiments. Aerobic batch cultures containing 30 ml of M9 medium were inoculated (1:100 -1:200) in 500-ml baffled shake flasks and incubated on a gyratory shaker at 250 rpm and 37°C. For 13 C-labeling experiments, glucose was added either entirely as the 1-13 C-labeled isotope isomer (Ͼ99%; Euriso-top, GIF-sur-Yvette, France) or as a mixture of 20% (w/w) U-13 C (Ͼ98%; Isotech, Miamisburg, OH) and 80% (w/w) natural glucose.
Cell growth was monitored by following the A 600 . Glucose and acetate concentrations were determined enzymatically using commercial kits (Beckman-Coulter (Zurich, Switzerland) or Dispolab (Dielsdorf, Switzerland)). Other organic acids in culture supernatants were detected by high pressure liquid chromatography analysis (PerkinElmer Life Sciences) at a wavelength of 210 nm, using a Supelcogel C8 column (4.6 ϫ 250 mm) at 30°C and a mobile phase of 2% (v/v) sulfuric acid at a flow rate of 0.3 ml/min.
The following physiological parameters were determined during the exponential growth phase as described previously (23): maximum growth rate, biomass yield on glucose, specific glucose consumption rate, and specific byproduct production rates, using a predetermined correlation factor of 0.44 g of cellular dry weight per liter and A 600 unit.
Metabolic Flux Ratio (METAFoR) 3 Analysis by Gas Chromatography-Mass Spectrometry-Samples for gas chromatography-mass spectrometry analysis were prepared as described previously (24). Briefly, aliquots of 13 C-labeled batch cultures were withdrawn during the midexponential growth phase (A 600 ϭ 0.8 -1.2). Cell pellets were hydrolyzed in 6 M HCl at 105°C for 24 h in sealed microtubes. The hydrolysates were dried under a stream of air at around 60°C and then derivatized at 85°C in 30 l of dimethylformamide (Fluka, Buchs, Switzerland) and 30 l of N-(tert-butyldimethylsilyl)-N-methyl-trifluoroacetamide with 1% (v/v) tert-butyldimethylchlorosilane (Fluka) for 60 min (25). Derivatized amino acids were analyzed on a series 8000 GC, combined with an MD 800 mass spectrometer (Fisons Instruments, Beverly, MA). The gas chromatography-mass spectrometry-derived mass isotope distributions of proteinogenic amino acids were then corrected for naturally occurring isotopes (24). The corrected mass distributions were related to the in vivo metabolic activities with previously described algebraic equations and statistical data treatment, which quantified 12 largely independent ratios of fluxes through converging reactions and pathways to the synthesis of seven intracellular metabolites (24).
In [1-13 C]glucose experiments with tpi mutants, 13 C label occurred in position 3 of pyruvate, but no position 3 13 C label occurred in the upstream metabolites phosphoglycerate and phosphoenolpyruvate (PEP). The normal flux ratio definition of METAFoR analysis (24) would relate this label to the Entner-Doudoroff (ED) pathway. This pathway, however, would introduce 13 C label at the position 1 of pyruvate. The normally inactive methylglyoxal bypass, which channels dihydoxyacetone phosphate (DHAP) molecules to pyruvate, in contrast, would precisely introduce such 13 C label at the position 3 of pyruvate (26). To quantify the amount of pyruvate originating from DHAP through the methylglyoxal bypass (pyruvate through methylglyoxal bypass), we determined the mass isotopomer distribution vector (MDV) of glycerol, which is identical to the MDV of DHAP 1-3 (1-3 indicates that carbon atoms 1-3 of DHAP are considered). The base fragment m 0 ϭ 377, which corresponds to a derivatized glycerol molecule (two tert-butyldimethylsilyl and one dimethylsilyl derivatization chain) was used. To assess the relative contribution of methylglyoxal bypass to pyruvate synthesis (f), the MDV of DHAP 1-2 , serine 2-3 , and pyruvate 2-3 are used as follows.
The MDV of DHAP 1-2 was not measured; however, the labeling enrichment for this fragment will be similar to the one of DHAP 1-3 , since position 1 of DHAP will be bearing the 13 C-labeled atom. Hence, the MDV of DHAP 1-2 can be calculated, and Equation 1 can be used to determine f. This approach neglects, however, the contribution of unla-beled two carbon molecules coming through the ED pathway. To take the ED pathway contribution into account, we applied an iterative process. The fraction of pyruvate molecules synthesized through the ED pathway was calculated from the previously determined methylglyoxylate flux ratio and the MDV of pyruvate 1-3 , DHAP 1-3 , serine 1-3 , and the MDV of a three-carbon unit with one labeled 13 C atom as follows.
where f 1 represents the fraction of pyruvate molecules originating from DHAP, f 2 represents one-half of the fraction derived through the ED pathway, and (1 Ϫ f 1 Ϫ f 2 ) represents the fraction derived through serine. Activity of the ED pathway introduces supplementary unlabeled pyruvate 2-3 molecules. To account for these, the MDV of serine 2-3 in Equation 1 was replaced by the MDV of a combination of unlabeled C 2 molecules and serine 2-3 . The fraction of pyruvate 2-3 molecules originating from DHAP 1-2 was recalculated, and the whole process was repeated iteratively until stable values for both of the fractions of pyruvate 2-3 molecules originating from DHAP 1-2 and of pyruvate molecules derived through the ED pathway were obtained (around five iterations). The S.E. values for these two determined flux ratios were evaluated numerically. Indeed, since all individual components in the mass distribution vectors have an S.D. value (24), normally distributed random values for these individual components were chosen using the MATLAB function normrnd (The Mathworks), with the constraint that the sum of the elements of a mass distribution is equal to 1. These new mass distribution vectors were used to determine the two ratios by repeating the process 1000 times in a MATLAB-based program. The mean values and the S.D. values for the two ratios were determined from these 1000 estimations. The mean values were within 1% of the calculated ratios, and the S.D. value was used as the error for the ratio.
Moreover, in the absence of triose phosphate isomerase in the tpi mutants, the hypothesis used to calculate the flux ratios PEP through the pentose phosphate (PP) pathway and serine through Embden-Meyerhof-Parnas (EMP) pathway is not valid anymore. Therefore, a [1-13 C]glucose experiment was used to determine the relative contribution of the EMP pathway to DHAP synthesis. In such a setup, all DHAP molecules derived through the EMP pathway contain a 13 C atom at position 1, whereas the PP pathway will generate unlabeled DHAP molecules. Therefore, the relative contribution of the EMP pathway to DHAP can be determined as follows.
where f represents the contribution of the EMP pathway to DHAP synthesis, and (1 Ϫ f) is the contribution of the oxidative branch of the PP pathway to DHAP synthesis. The error on this ratio was determined using error propagation (24). 13 C-Constrained Net Flux Analysis-Intracellular net fluxes were estimated with a stoichiometric model that contained all major pathways of central carbon metabolism (25). For the tpi mutants, the previously described stoichiometric model was augmented with the methylglyoxal bypass based on the METAFoR results. The network considered was similar to the one depicted in Fig. 1. For all mutant analyses, the deleted reactions were kept in the network to obtain independent evi-dence for their in vivo absence (or evidence for the takeover by another gene). Only for the ppc mutants was the reaction from 2-oxoglutarate to fumarate removed from the network, since METAFoR analysis demonstrated a complete absence of cyclic tricarboxylic acid cycle operation (see Supplemental Table 1). The reaction matrix consisted, for the different strains, of 25-29 unknown fluxes and 21-24 metabolite balances (including the three experimentally determined rates of glucose uptake, acetate, and biomass production).
To solve this underdetermined system of equations with 4 -5 degrees of freedom, the following seven calculated flux ratios were used as additional constraints, as was described previously (25): serine derived through the EMP pathway, pyruvate derived through the ED pathway, oxaloacetate originating from PEP, PEP originating from oxaloacetate, pyruvate originating from malate (upper and lower boundaries), and PEP derived through the PP pathway (upper boundary). The first four ratios were used as equality constraints, whereas the others were used only as boundary constraints. When active, based on the METAFoR data, the glyoxylate shunt was also considered in the network, and the ratio oxaloacetate originating from glyoxylate was implemented as an upper bound. 4 The ratios DHAP derived through the EMP pathway and pyruvate through the methylglyoxal bypass were used as equality constraints for the tpi mutants, using the following equations: the fraction of pyruvate derived through the methylglyoxal bypass.
Fluxes into biomass were calculated from the known metabolite requirements for macromolecular compounds (28) and the growth ratedependent RNA and protein content (29). The sum of the weighed square residuals of the constraints from both metabolite balances and flux ratios was minimized using the MATLAB function fmincon. The residuals were weighed by dividing through the experimental error (25). The computation was repeated at least five times with randomly chosen initial flux distributions to ensure identification of the global minimum, and the system always converged to the same solution. An extended version of the software FiatFlux was used to calculate all metabolic flux ratios and net fluxes (30). mRNA Transcriptional Profiling-Affymetrix (Santa Clara, CA) E. coli antisense genome arrays were used for all transcriptional analyses. Each experimental condition was tested in triplicate using independent cultures and processed following the manufacturer's recommended protocols. Six replicates of the wild-type strain grown on glucose were used for the reference point. Briefly, cultures were grown to midexponential growth phase (A 600 Ϸ 0.5). 3 ml of culture was added to 6 ml of RNAprotect (Qiagen, Valencia, CA), and RNA was isolated using RNeasy kits (Qiagen, Valencia, CA) following the manufacturer's instructions. Total RNA yields were measured using a spectrophotometer (A 260 ), and quality was checked by visualization on agarose gels and by measuring the sample A 260 /A 280 ratio. cDNA synthesis, fragmentation, and terminal labeling were conducted as recommended by Affymetrix. Raw .CEL files were analyzed using robust multiarray average (31) for normalization and calculation of probe intensities.
Expression values were then assessed for statistically significant differential expression using t tests. After conducting pairwise t test comparisons between evolved mutants and wild type, those genes meeting a 5% false discovery rate-adjusted p value cut-off were chosen as having statistically significant changes in gene expression. This resulted in selection of subsets of differentially expressed genes for each tested evolution mutant.
These subsets of differentially expressed genes were then organized into known regulon structures (32) for further analysis. The probability (p value) of the observed regulon enrichment of differentially expressed genes was calculated using the hypergeometric distribution (33), where N (equal to 4345) represents the total number of E. coli genes listed on the Affymetrix GeneChip, r is the total number of genes that are a part of the regulon, n is the number of differentially expressed genes, and y is the number of genes that are differentially expressed and a member of the regulon.
In Vitro Enzyme Activity-Crude cell extracts for in vitro enzyme assays were prepared from pellets of 45-ml culture aliquots. Pellets were resuspended in 4 ml of 0.9% (w/v) NaCl and 10 mM MgSO 4 , passed three times through a French pressure cell (1.2 ϫ 10 8 pascals), and centrifuged at 11,000 ϫ g. In vitro activity of the methylglyoxal bypass was determined by enzyme assays for methylglyoxal reductase and two coupled enzyme assays for methylglyoxal synthase-methylgloxal reductase and methylglyoxal synthase-glyoxalase. For methylglyoxal reductase, the assay of Saikusa et al. (34) was adapted: 50 mM Tris-HCl (pH 7.5), 0.125 mM NADPH, 10 mM methylglyoxal, and 100 l of crude cell extract per ml of assay. The assay for methylglyoxal synthase-methylglyoxal reductase was similar to the previous one, with the sole difference being the use of 0.75 mM DHAP instead of methylglyoxal. For methylglyoxal synthase-glyoxalase, the assay of Hoper and Cooper (35) was adapted: 40 mM imidazole (pH 7.0), 1.65 mM glutathione (pH 7.0), 0.75 mM DHAP, and 100 l of crude cell extract per ml of assay. For all three assays, the reaction was initiated by the addition of either methylglyoxal or DHAP and consumption of NADPH (at 340 nm), and formation of S-lactoylglutathione (240 nm) was detected in the two first and second assays, respectively. The extinction coefficients for NADPH and S-lactoylglutathione were 6.2 and 3.4 mM Ϫ1 cm Ϫ1 , respectively. The protein content in crude cell extracts was determined with the biuret reaction.

Physiological Characterization of Evolved Knock-out Mutants-To
identify the metabolic routes chosen by adaptive evolution to cope with a gene knock-out, we selected four gene knock-outs that severely affect the flux distribution at key branch points of glucose catabolism: pgi, ppc, pta, and tpi (Fig. 1). Two parallel cultures of each mutant were evolved for 30 -50 days (600 -800 generations), under conditions of exponential batch growth in M9 minimal medium with 2 g/liter glucose. The end points of evolution were defined as having unaltered physiology for Ͼ100 generations, and these mutant populations were designated as ptaE1, ptaE2, ppcE1, ppcE2, pgiE1, pgiE2, tpiE1, and tpiE2. Akin to the previously reported ppc and tpi evolution experiments (3), the pgi and pta mutants rapidly evolved improved phenotypes (Table 1). Most knock-out mutations severely reduced the specific growth and glucose consumption rates, but all evolved mutants largely recovered the wild type-like rates at the end point of evolution. Whereas all mutants evolved to improved phenotypes, some parallel evolved cultures displayed different phenotypes; e.g. the evolved pgi mutants pgiE1 and pgiE2 exhibited very different overflow metabolism and biomass yield, and the ptaE1 and ptaE2   mutants differed significantly in their glucose uptake rates (Table 1). Thus, the evolutionary end points are convergent and reproducible with respected to the main selective pressure for high growth rate, but the underlying physiological states are not necessarily the same. Metabolic Flux Analysis-Which are the actual metabolic mechanisms chosen by adaptive evolution, and do all mutants evolved in parallel rely on the same mechanisms to cope with a given lesion? To answer these questions, intracellular fluxes were quantified from 13 Clabeling experiments (25).
As expected (24,36,37), glucose catabolism in the unevolved pgi mutant is rerouted from the EMP to the PP and ED pathways to bypass the lesion (Fig. 2A). Furthermore, operation of the otherwise inactive glyoxylate shunt in the unevolved pgi mutant was consistent with earlier reports (37,38). Whereas the severely altered flux distribution was maintained, evolution more than doubled the absolute flux level to about 65% of the wild-type glucose uptake rate (Table 1). In contrast to the unevolved pgi mutant, initial glucose catabolism proceeded almost exclusively through the PP pathway in both evolved mutants ( Fig. 2A).  Table 1 (25). Physiological data were allowed to vary within the confidence intervals given in Table 1. Reactions encoded by deleted genes are highlighted in gray and were not removed from the flux model. Reactions not considered in the network because of independent 13 C flux ratio evidence for their absence are represented by a minus sign. Relative to the glucose uptake rate, confidence intervals were between 10 and 40% for the transhydrogenase flux, below 20% for the tricarboxylic acid cycle, and less than 10% for all other fluxes. A complete list of determined absolute fluxes is shown in Supplemental Table 3.
In the lower part of metabolism, however, the two evolved network topologies differed significantly. Whereas the pgiE2 mutant flux distribution was similar to the unevolved mutant with an active glyoxylate shunt and no acetate secretion, the pgiE1 mutant was more similar to the wild type with acetate secretion and full tricarboxylic acid cycle flux instead of the glyoxylate shunt flux. Thus, the pgi mutants evolved to improved phenotypes through rather different intracellular flux scenarios. Unexpectedly, these flux adaptations occurred far away from the deleted gene, possibly suggesting the need for downstream metabolic adjustment in relation to the initially implemented flux distribution. Knock-out of phosphoglucose isomerase forced glucose catabolism primarily through the NADPH-producing PP pathway, resulting in NADPH production encompassing the biosynthetic needs of the cells (36). Hence, in contrast to the wild type, where the membrane-bound transhydrogenase PntAB converts NADH to NADPH to meet the NADPH requirements of the cells, in the pgi mutants, the soluble transhydrogenase converts the excess NADPH to NADH ( Fig. 2A).
Generally, PEP carboxylase inactivation precludes growth on glucose as the sole carbon source, because this anaplerotic reaction replenishes tricarboxylic acid cycle intermediates that are withdrawn for biosynthesis (39). Physiological suppressor mutants occur rapidly, however, and the slow glucose growth phenotype of the unevolved, suppressed ppc mutant was almost fully recovered through adaptive evolution with indistinguishable physiology in the two evolved mutants (Table 1). In all ppc mutants, the normally glucose-repressed glyoxylate shunt replaced the anaplerotic function of PEP carboxylase (Fig. 2B), as was described previously (40). The flux through the shunt, however, was in 40% excess of the anabolic demand for the biomass precursors oxaloacetate and 2-oxoglutarate. In the unevolved and the ppcE1 mutants, this excess flux was catabolized through the PEP-glyoxylate cycle with the glyoxylate shunt and PEP carboxykinase as key reactions (38). Additionally, malic enzyme contributed to this cycle in all three ppc mutants by converting malate to pyruvate. Thus, the latent glyoxylate shunt functionally replaced the ppc mutation and was immediately invoked upon deletion of ppc. Evolution of the ppc mutants led to subtle differences in the metabolic network topology by using either PEP carboxykinase or malic enzyme as was necessary to balance the excess precursors being generated through the active glyoxylate shunt.
Blocking the main acetate secretion route in pta mutants was counteracted by secreting pyruvate instead of acetate ( Table 1). As had been observed for different pta mutants, small amounts of acetate were still produced (41). Since the phenotype of the unevolved mutant was otherwise similar to the wild type, growth physiology and network topology were largely unaltered in the evolved mutants ( Fig. 2C and Table 1). Nevertheless, the end point flux states were detectably different in the two evolved pta mutants with significantly lower absolute fluxes and lower pyruvate secretion in the ptaE2 mutant. The increased tricarboxylic acid cycle and EMP pathway expressions, which had been observed for a double pta-ackA knock-out mutant (42), were not reflected in the flux states of the pta mutants.
Knock-out of the triose-phosphate isomerase in the tpi mutant affects a stoichiometrically equal splitting of the glycolytic flux into glyeraldehyde phosphate and DHAP. To prevent internal accumulation of DHAP, tpi mutants convert DHAP to pyruvate through the normally inactive methylglyoxal bypass (26,43) (Fig. 2D). Since the present 13 C data provide only indirect evidence for methylglyoxal bypass fluxes, its activation was confirmed through in vitro enzyme data. Compared with the wild type, the unevolved tpi mutant exhibits about 2.5-fold higher in vitro activities in the glyoxalase I branch of the methylglyoxal bypass ( Table 2). Evolution more than doubled the overall flux level to about 80% of the wild-type glucose uptake rate (Table 1). This is probably a direct consequence of improved methylglyoxal bypass fluxes through the glyoxalase I branch with about 4-fold increased in vitro activities ( Table 2), indicating that the glutathione intermediate branch is more suitable for higher fluxes through the methylglyoxal bypass than the two consecutive oxidation-reduction reactions in the methylglyoxal reductase branch. Thus, the normally latent methylglyoxal bypass functionally replaced the introduced tpi mutation. Whereas the unevolved strain exhibited slow growth, evolution led to greatly improved growth through implementation of metabolic adjustments (methylglyoxal bypass) downstream of the introduced lesion needed to accommodate different metabolite pools present due to the absence of triose-phosphate isomerase, which resulted in almost tripled absolute fluxes.
mRNA Transcriptional Profiling-Whereas the flux data identified metabolic mechanisms that cope with lesions and the changes brought about by evolution, they cannot reveal the actual genetic basis. In an attempt to elucidate the molecular mechanisms responsible for the observed changes in network topology during evolution, we generated differential genome-wide transcription profiles of the evolved pgi, ppc, and tpi mutants that exhibited major flux differences. All results and analyses were restricted to statistically significant expression changes (as determined by p value cut-offs selected for a false discovery rate of 5%). Comparing the wild type with the evolved populations, 1205 differentially expressed genes were identified. All significant changes in gene expression were qualitatively correlated with the flux changes through the encoded reaction by relative comparison of evolved expression levels and flux changes against the wild-type base-line measurements.
In the pgiE1 and pgiE2 mutants, 19 fluxes representing 35 genes and 26 fluxes representing 53 genes, respectively, changed significantly when compared with the wild type. A qualitative agreement between flux and expression changes occurred in 7 (of 35) and 35 (of 53) of these genes for mutants pgiE1 and pgiE2, respectively ( Fig. 3A and Supplemental Table 2). Besides the expected decrease in expression of pgi, expression of the glycolytic genes pfkA and gapA was consistently reduced and correlated with the lower glycolytic flux in both mutants. Consistent with the differences in phenotype and flux profiles, a large number of expression differences were found between the two mutants in the lower portion of central metabolism. In particular, decreased expression of the tricarboxylic acid cycle genes icdA, sucABCD, and Adaptive Evolution of E. coli MARCH 24, 2006 • VOLUME 281 • NUMBER 12 sdhABCD and, less pronounced, of all glycolytic enzymes in pgiE2 was consistent with the flux data. Moreover, increased expression of the glyoxylate shunt gene glcB correlated with the activated shunt in this pgiE2 mutant, and no significant difference in expression was seen in the pgiE1 mutant with an inactive shunt. Thus, there is a genetic basis for the much lower tricarboxylic acid cycle and much higher glyoxylate shunt in the pgiE2 mutant and the absence of significant changes in the pgiE1 mutant. E. coli has two transhydrogenases, a soluble and a membrane-bound, which are used to balance the NADPH and NADH pools (44 -46). The soluble transhydrogenase, encoded by udhA, converts NADPH to NADH, whereas the membrane-bound transhydrogenase, encoded by pntAB, converts NADH to NADPH and accounts, in batch cultures of the wild type, for around 40% of the required NADPH production (36). Expression changes for the two transhydrogenases mostly correlated with the flux distribution for both evolved mutants. Indeed, the expression of the membrane-bound transhydrogenase showed a statistically significant decrease in expression, since NADPH was produced in excess to the biosynthetic requirements, whereas the expression of the soluble transhydrogenase increased, but surprisingly was significant only for pgiE2 that had the lower NADPH to NADH conversion ( Fig. 2A).
In both evolved ppc mutants, altered expression of 14 genes was qualitatively correlated with flux changes (Fig. 3B and Supplemental Table  2). Common to both mutants was altered expression in the lower part of central metabolism with decreased expression of ackA and the tricarboxylic acid cycle genes sucABC and sdhABCD and increased expression of fumB and the glyoxylate shunt gene aceA. Thus, decreased ace-tate secretion, activation of the glyoxylate shunt, and decreased tricarboxylic acid cycle fluxes were probably directly determined through altered expression of key genes in these pathways.
In the tpiE1 and tpiE2 mutants, 15 fluxes representing 42 genes and 16 fluxes representing 43 genes, respectively, changed significantly when compared with the wild-type. Qualitative correlation between expression and flux changes, however, was observed for only two genes (one of which was decreased expression of tpi) in tpiE1 and four genes in tpiE2 ( Fig. 3C and Supplemental Table 2). Notably, the major flux change in the mutants compared with the wild-type, activation of the normally latent methylglyoxyal bypass, was probably genetically determined, because both mutants increased expression of gloA (Ͼ2-fold). This view is further supported by the about 10-fold higher in vitro activity of the gloA-encoded glyoxylase I ( Table 2). The decreased glycolytic and increased tricarboxylic acid fluxes were not reflected by changes in gene expression.
Since higher glyoxylate shunt fluxes appeared to have a genetic basis, we were interested to see whether this was due to a particular mechanistic change affecting only the shunt genes or if other expression changes were correlated in individual genes or in a regulatory cascade. Hence, we searched for genes with statistically significant expression changes in the same direction in the pgiE2, ppcE1, and ppcE2 mutants with an active glyoxylate shunt. Expression of 38 genes was increased, and expression of 132 genes was decreased exclusively in the pgiE2, ppcE1, and ppcE2 mutants but not in the other three mutants investigated. Mostly, these genes were involved in metabolic pathways that branch off from the tricarboxylic acid cycle, including increased expres- sion of tRNAs associated with glutamine (glnX), asparagine (asnTUVW), and methionine (metVWYZ) and decreased expression in some of the biosynthetic pathways for aspartate (ansA) and methionine (metABCE). In addition, redox metabolism was affected in the activated glyoxylate shunt mutants with increased expression of the hyaA, torA, torC, torY, and torZ genes that are involved in quinone biosynthesis from the tricarboxylic acid cycle intermediate ␣-ketoglutarate. This coordinated pattern of expression changes may indicate one or more mutations in the transcriptional network that controls expression of the glyoxylate shunt and many other genes. Whereas the active glyoxylate shunt was clearly relevant for the observed phenotypes, it remains unclear whether the related changes contribute to the mutant phenotypes.
Changes in transcriptional regulatory units were further investigated by using significant expression changes in each strain to calculate the probability that a statistical change has occurred in a regulon. p values for each evolved mutant were calculated for 124 regulons using the hypergeometric distribution at a p value cut-off of 5%. Screening for regulatory changes that could correspond to physiological changes, we found that both evolved pgi mutants exhibited consistent down-regulation of the LeuO regulon (associated with leucine biosynthesis) and that all three mutants utilizing the glyoxylate shunt (pgiE2, ppcE1, and ppcE2) down-regulated the MetJ regulon (associated with methionine biosynthesis) ( Table 3). Whereas production of the amino acids leucine and methionine is essential, it appears that a change in the regulatory network was induced to balance the amino acid with other biological demands. In the case of the pgi mutants, lower availability of pyruvate (a precursor to leucine) could force a reduction in leucine production to allow pyruvate to fulfill other metabolic needs. In the case of mutants utilizing the glyoxylate shunt, activation of the PEP(pyruvate)-glyoxylate cycle (38) was found, which involves oxaloacetate (a precursor to methionine) and thus could limit the production of methionine. Thus, in two cases, adaptive mechanisms occurred in the transcriptional regulatory network that are closely connected to the observed changes in the metabolic network.

DISCUSSION
Several strategies can be invoked in response to externally introduced genetic perturbations that confer genetic robustness to the network. The initial response to the four investigated deletions was a local rerouting of fluxes around the lesion, which involved activation of latent pathways in several mutants. Although cultures were followed for several hundred generations, there was not a single case where evolution invented a new solution (e.g. activation of a silent gene or reassignment of enzyme function). Instead, the primary effect of evolution over the observed time frame was an increase in the capacity of already active pathways. As a consequence, the overall fluxes increased in those mutants that were severely affected by the mutation (pgi, ppc, and tpi mutants). In several cases, the first metabolic response in the unevolved mutants, which involved the activation of pathways that are usually inactive or almost inactive, such as the glyoxylate shunt or the ED pathway, was even lost during adaptive evolution. Examples are the loss of glyoxylate shunt activity in pgiE1 and the two evolved tpi mutants and the loss of local flux rerouting through the ED pathway in the evolved pgi mutants. Hence, the immediate strategy invoked by E. coli in response to an externally introduced genetic perturbation is not always optimal.

TABLE 3 Regulon changes in the MetJ and LeuO regulons
p values calculated using the hypergeometric distribution indicate significance of expression changes within each regulon. The mean fluorescence intensity of the wild-type expression level is given as the Log2 value. Expression changes for each gene within the regulon are given as Log2 ratios. Shaded boxes showed statistical change (p value Ͻ 0.05) in the regulon.
However, no single optimum exists, as demonstrated by the parallel evolved pgi mutants that differ in anaplerotic reactions and tricarboxylic acid cycle fluxes.
Whereas there are a number of ways for an organism to compensate for a gene deletion, one hypothesis is that an introduced genetic perturbation should force a redistribution of fluxes through the network as an immediate rescue solution, and evolution would then entail a process of refining this newly established initial state (47). Generally, our results fully support this hypothesis, but we were surprised to find relatively large variations in the flux states of parallel evolved mutants. This appears to reflect the robustness of the metabolic network and its ability to utilize different means to not only survive but to improve growth characteristics. These results provide experimental support to computational results stating that bacterial metabolic networks have many different means of achieving similar, equivalent functionality (alternate optimal solutions) (48). It is likely that biological noise and stochastic variations play some role in the development of these diverse flux states, but it is presently unclear how these properties are integrated into evolutionary dynamics, and selection is not clearly defined.
The distribution of molecular fluxes is regulated by multiple mechanisms at several levels that include gene expression, posttranscriptional control, enzyme kinetics, and allosteric control. To determine whether transcriptional modifications were the cause of the altered flux distributions, gene expression levels were systematically compared with flux levels. A particularly high correlation was found for the complete activation or inactivation of latent pathways, such as the glyoxylate shunt and the methylglyoxal bypass. In those cases, expression of at least one gene correlated qualitatively with the flux changes through these pathways. Moreover, flux changes in the tricarboxylic acid cycle also correlated mostly with altered gene expression, e.g. the reduced fluxes in parts of the tricarboxylic acid cycle for pgiE2 and both evolved ppc mutants. The lack of correlation in the tpi mutants might be related to the less pronounced changes in absolute tricarboxylic acid cycle fluxes than for the other three mutants, whose succinyl-CoA synthetase and succinate dehydrogenase fluxes were essentially zero ( Fig. 2 Table 3). No flux-expression correlation was found for the three pathways of initial glucose catabolism, not even in the pgi mutants with greatly increased PP pathway fluxes. Solely pfkA and gapA expression correlated with the reduced EMP pathway flux in both pgi evolved mutants. Hence, glyoxylate shunt, methylglyoxal bypass, and tricarboxylic acid cycle flux changes appear to be controlled, at least in part, at the transcriptional level. Flux through the PP and EMP pathways and PEP carboxylase, in contrast, were not controlled at the transcription level.

and Supplemental
Similarly, a strong qualitative correspondence between gene expression and metabolic fluxes for the glyoxylate shunt had previously been observed in S. cerevisiae, whereas tricarboxylic acid cycle and PP pathway fluxes only partially correlated (49). Comparison of flux and gene expression for E. coli grown under anaerobic conditions on xylose and glucose showed a similar absence of correlation for the PP pathway; however, in contrast to our results, the EMP pathway showed a strong flux-gene expression correlation (27). The absence of correlation in our data might be related to the relatively small changes in absolute EMP pathway fluxes. Hence, gene expression changes are not always manifested in the expressed phenotype, since attributes such as translational efficiency, allosteric control, or changes in enzyme kinetics will not be reflected in individual mRNA transcript levels. Thus, the combined analysis of gene expression data with flux data is one means of pinpointing mRNA transcript and regulatory changes that may be causal to observed phenotypes. In the case of this study, interpreting the gene expression data in the context of flux measurements allowed us to focus on a small number of meaningful expression changes out of the thousands of observed expression changes.
Generally, E. coli is able to rapidly recover to nearly wild-type growth from a severely crippled phenotype resulting from initial genetic perturbations by utilization of existing (though sometimes dormant) pathways parallel to the lesion. This conclusion can probably be extrapolated to mutations in the central metabolism of many organisms that, like E. coli, feature a highly interconnected network of core reactions but may differ for organisms with simpler metabolism and for mutations in less redundant parts of the network (e.g. biosynthetic routes to essential compounds). Metabolic capacity of the activated parallel pathways appears to be analogous to the pathways used in the wild type; however, downstream metabolic adjustments are often needed to refine usage of the new pathways. These adjustments could be implemented to alleviate intracellular metabolite pools resulting from new pathway usage (pgi, ppc, and tpi mutants). It was even observed that strains evolved in parallel frequently utilized different means of achieving this refinement, sometimes leading to surprisingly large differences in flux states (pgiE1 and pgiE2). Overall, we have found that the metabolic network of E. coli is robust in response to a genetic perturbation, with metabolic adjustments occurring in two phases: an initial flux rerouting to compensate for the lesion and a subsequent downstream adjustment to optimize flux rerouting. Combined analysis of multiple data types was pivotal to unravel mechanistic details of these downstream adjustments.