Modeling genome-wide enzyme evolution predicts strong epistasis underlying catalytic turnover rates

Abstract

Systems biology describes cellular phenotypes as properties that emerge from the complex interactions of individual system components. Little is known about how these interactions have affected the evolution of metabolic enzymes. Here, we combine genome-scale metabolic modeling with population genetics models to simulate the evolution of enzyme turnover numbers (k_cats) from a theoretical ancestor with inefficient enzymes. This systems view of biochemical evolution reveals strong epistatic interactions between metabolic genes that shape evolutionary trajectories and influence the magnitude of evolved k_cats. Diminishing returns epistasis prevents enzymes from developing higher k_cats in all reactions and keeps the organism far from the potential fitness optimum. Multifunctional enzymes cause synergistic epistasis that slows down adaptation. The resulting fitness landscape allows k_cat evolution to be convergent. Predicted k_cat parameters show a significant correlation with experimental data, validating our modeling approach. Our analysis reveals how evolutionary forces shape modern k_cats and the whole of metabolism.

Fig. 1

The MCMC algorithm used for simulating genome-scale k_cat evolution. A single iteration of the algorithm proceeds as follows: (I) A mutation in the k_cat of a random reaction of a single cell in the population is introduced. The original growth rate µ(k_cat) and the novel growth rate µ(k_cat,mutated) are predicted by solving the respective MOMENT problems (see Methods). (II) The probability of fixation for the novel mutation is calculated with a population genetics model based on µ(k_cat), µ(k_cat,mutated), and the population size N. Fixation of the novel change in k_cat is then decided based on this probability. If fixation fails, the mutation is discarded. A typical simulation run includes around 10⁸ of the described iterations

Fig. 2

Evolutionary trajectories exhibit convergence and diminishing returns epistasis. a The growth rate of the population against the number of simulated mutations. The black line shows the average growth rate across replicates and its standard error. The replicates showed 4880 fixation events on average. See Supplementary Figure 3 for a fit of the analytical model presented in Supplementary Note 1. b The selection coefficient s (defined as the change in growth rate relative to the novel growth rate) plotted against the cumulative number of simulated mutational events for all fixed mutations. The inset shows s against the background growth rate in which a mutation occurred for the three reactions that had the most changes in k_cat fixed. These reactions are: ATPS4rpp: ATP synthase (orange), GLUDy: Glutamate dehydrogenase (NADP) (green), and PPC: Phosphoenolpyruvate carboxylase (blue)

Fig. 3

Multifunctional enzymes cause synergistic epistasis in k_cat evolution. a A multifunctional enzyme with two distinct active sites catalyzes two reactions in the same linear fitness-relevant pathway. b Mutations that increase either k_cat individually cannot be used to reduce protein cost of the pathway and thus exhibit synergistic epistasis

Fig. 4

Comparison between k_cat predictions for evolution in diverse environments and experimental data. a Distribution of turnover rates in in vitro (n = 188), in vivo (n = 210), and simulated data (n = 276). Simulated data are only shown for non-constrained reactions that contribute to growth. The arrow indicates reactions that were essentially left unchanged by the simulation, indicating that they were not used in most environments. b Comparison between experimental in vitro data and simulated data for reactions contributing to growth (n = 83). c Comparison between experimental in vivo data (k_app,max) and simulated data for reactions contributing to growth (n = 210). The outliers in the upper left suggest that these reactions are rarely used in the environmental conditions that we model. Horizontal error bars in b and c show the standard deviation across three simulated replicates. The p-values in b and c are based on Pearson’s R to test for significant correlation with a two-sided t-test (see Methods). See Supplementary Figure 13 for sensitivity analysis against assumptions about the ancestral state and Supplementary Figure 15 for sensitivity analysis against reaction stoichiometries