Evolution of gene knockout strains of E. coli reveal regulatory architectures governed by metabolism

Abstract

Biological regulatory network architectures are multi-scale in their function and can adaptively acquire new functions. Gene knockout (KO) experiments provide an established experimental approach not just for studying gene function, but also for unraveling regulatory networks in which a gene and its gene product are involved. Here we study the regulatory architecture of Escherichia coli K-12 MG1655 by applying adaptive laboratory evolution (ALE) to metabolic gene KO strains. Multi-omic analysis reveal a common overall schema describing the process of adaptation whereby perturbations in metabolite concentrations lead regulatory networks to produce suboptimal states, whose function is subsequently altered and re-optimized through acquisition of mutations during ALE. These results indicate that metabolite levels, through metabolite-transcription factor interactions, have a dominant role in determining the function of a multi-scale regulatory architecture that has been molded by evolution.

Fig. 1

Evolution of knockout (KO) strains from a pre-evolved (i.e., optimized) wild-type strain. a Experimental design using adaptive laboratory evolution (ALE) and enzyme knockouts to investigate system re-optimization following major metabolic perturbations. b An isolated wild-type (wt) E. coli (MG1655 K-12) previously evolved on glucose minimal media at 37 °C was used as the starting strain for knockouts of key metabolic genes and subsequent re-evolution, or systems re-optimization. c Reactions disabled by the enzyme knockouts included the phosphotransferase sugar import system (ptsHIcrr), phosphoglucose isomerase (pgi), 6-phosphogluconate dehydrogenase (gnd), triophosphate isomerase (tpi), and succinate dehydrogenase complex (sdhCB). d Adaptive laboratory evolution trajectories of the initial reference knockout and evolved knockout lineages. e Counts of significantly different system components found for each evolved knockout relative to the unevolved knockout. Counts of metabolomic, transcriptomic, and fluxomic data are given as the average and standard deviation of the percent of significant features compared to all features measured for the lineage; counts for mutations are given as the average and standard deviation of the number of significant features (see Methods for criteria for significance)

Fig. 2

A multivariate analysis of biological network components as represented by different omics data types. a Partial least squares discriminatory analysis (PLS-DA) revealed a common trend in the two most dominant components: the primary component (PC1) most often corresponded to a movement away from (dashed line) and back to (solid line) evolved optimal fitness (i.e., optimal system configuration), while the secondary component (PC2) most often corresponded to a diversity among evolved optimal fitness states of different lineages (i.e., optimal system configurations). PLS-DA scores plots of the reference strain, initial knockout, and evolved endpoints for each lineage for metabolomics (b), transcriptomics (c), and fluxomics (d) data. The strain lineages denoted on the top of b also refer to the corresponding graphs below in c and d. All of the KO lineages matched the trend described above in the metabolomics data, one eKO did not match the trend in four of the five KO lineages in the expression data (i.e., all but eSdhCB), and one or more eKO did not match the trend in each of the KO lineages in the fluxomics data (see Methods for thresholds)

Fig. 3

Classification of changes in omics data between the reference strain (Ref), the unevolved knockout strains (uKOs) and evolved knockout strains (eKOs). a Individual components were mapped onto six profiles according to their abundance in the Ref, uKOs, and eKOs in both positive and negative directions. The six profile types are shown and include novel, overcompensation, partially restored, reinforced, restored, and unrestored. The novel profile sought to categorize system components that were changed only as a result of adaptation. The restored and partially restored profiles sought to categorize system components that were initially perturbed to suboptimal levels post-KO. The overcompensation profile sought to categorize system components that overshot a restored profile to levels even lower/higher than those in Ref. The reinforced profile sought to categorize system components that needed a further increase or decrease after an initial perturbation post-KO to reach an optimal level. The unrestored profile sought to categorize system components that were immediately adjusted to an optimal level post-KO and required no further adjustments during adaptation. Metabolomics (b), transcriptomics (c), and fluxomics (d) data for each replicate of each KO lineage were binned into each of the six profiles (Pearson’s R, R > 0.88). See section “Component profiles reveal systematic variations” for a discussion of trends that were found based on these six profiles

Fig. 4

Suboptimal pathway usage limits allocation of carbon to biomass precursors. Toy network schematic of flux distribution in Ref (a) and in uKO (b). A reaction knockout is highlighted in red. The flux distribution in eKO could be categorized as c changed flux distribution (i.e., the pathway usage was changed) or d changed flux capacity (i.e., the same pathways was used but at a higher level, see Methods). Four examples of changed flux distribution and changed flux capacity for f gnd, g sdhCB, h ptsHIcrr, and i pgi lineages. f flux was initially re-routed through the ED pathway after removing the gnd gene The ED pathway has a net yield of one ATP, NADH, and NADPH per molecule of glucose, whereas glycolysis has a net yield of two ATP and NADH. Instead, the evolved gnd endpoints limited the use of the PPP and increased the flux capacity through the higher energy and redox equivalent producing pathway of glycolysis. g flux was initially re-routed through the TCA cycle in uSdhCB by diverting flux through the anaplerotic reactions phosphoenolpyruvate carboxylase (PPC) and inverting the direction of flux through malate dehydrogenase (MDH). The eSdhCB re-inverted the direction of malate dehydrogenase towards production of nadh or quinone reduction, and downregulated flux through the rest of the TCA cycle. h A significant portion of flux was bifurcated between the methylglyoxal pathway and lower glycolysis in uPtsHIcrr in response to elevated levels of dihydroxyacetone phosphate (DHAP) and depletion of lower glycolytic intermediates that inhibit the activity of methylglyoxal synthase^, . The flux through the methylglyoxal pathway was essentially eliminated in endpoints 2 and 4, and significantly decreased in replicates 1 and 3, in order to utilize the less toxic and more energy and redox producing lower glycolytic pathway. i The abnormally high levels of flux directed through the oxidative Pentose Phosphate Pathway (oxPPP) in uPgi was initially re-routed through the ED pathway. Several evolved pgi endpoints retained the flux through the ED pathway to varying degrees, but most re-distributed flux through GND, and all increased the flux capacity through the non-oxidative Pentose Phosphate Pathway (non-oxPPP). Green and orange colored reaction lines in f–i correspond to the grouping of changed flux distribution or changed flux capacity shown in the bar plot in h. Color bars for all flux values are shown next to their corresponding reaction(s)

Fig. 5

Mapping between network components and annotated regulation. a An algorithm for determining agreement and disagreement between system components categories and annotated biochemical pathways and regulation. 1. The algorithm inputs include the component profiles, the network components, and the network interactions. 2. An on/off Boolean interaction network that describes the biochemical and/or regulatory relationship between two components is constructed. 3. The component categories and on/off interaction between each component can then be determined. 4. For components that were not directly measured, a consensus category and confidence score can be determined. b Example of metabolite-mediated transcription factor activation between tyr-L, TyrR, and aroF. c Example of an unresolved discrepancy involving Fur regulation. d Example of transcription factor hierarchy between cAMP-CRP and Cra

Fig. 6

Overview of mutation statistics. See Supplementary Data 6 for detailed statistics of each category and categories not shown. a The type of mutation. Mutations include amplification (AMP), deletion (DEL), insertions (INS), mobile element aided insertions or deletions (MOB), single nucleotide polymorphism (SNP). b The location of the mutation. Locations include coding regions, regions associated with cryptic prophages, intergenic regions, regions two coding genes not classified as an intergenic region (intergenic/intergenic), and repetitive elements (REP or RIP). c The class of mutation. Classes include frameshifts, frameshifts resulted in a truncated CDS, missense, non-frameshifts, peptide truncations, and other unclassified mutations. d The functional or structural category of the mutated gene. Categories are based on the “parent class” as found in the EcoCyc database

Fig. 7

A model of biological systems adaptation following the KO of key metabolic enzymes. (i) Suboptimal pathway usage limited allocation of carbon to biomass precursors. (ii) Perturbed metabolite levels triggered transcription regulatory network (TRN) responses in the uKOs. (iii) Activation of the TRN revealed a hierarchy of regulation involving competing and overlapping regulatory interactions between various system components including DNA, RNA, and proteins. (iv) Mutations selected during adaptive evolution changed many regulatory networks, and also introduced innovations that targeted specific pathway or metabolite imbalances