Genome-driven evolutionary game theory helps understand the rise of metabolic interdependencies in microbial communities

Abstract

Metabolite exchanges in microbial communities give rise to ecological interactions that govern ecosystem diversity and stability. It is unclear, however, how the rise of these interactions varies across metabolites and organisms. Here we address this question by integrating genome-scale models of metabolism with evolutionary game theory. Specifically, we use microbial fitness values estimated by metabolic models to infer evolutionarily stable interactions in multi-species microbial "games". We first validate our approach using a well-characterized yeast cheater-cooperator system. We next perform over 80,000 in silico experiments to infer how metabolic interdependencies mediated by amino acid leakage in Escherichia coli vary across 189 amino acid pairs. While most pairs display shared patterns of inter-species interactions, multiple deviations are caused by pleiotropy and epistasis in metabolism. Furthermore, simulated invasion experiments reveal possible paths to obligate cross-feeding. Our study provides genomically driven insight into the rise of ecological interactions, with implications for microbiome research and synthetic ecology.

Fig. 1

General scheme for the proposed genome-driven evolutionary game theory approach. Annotated genomes of community members are used to construct genome-scale metabolic models. For each possible pair of genotypes in the community, constraint-based analysis tools for metabolic models, such as flux balance analysis, are used to estimate the fitness (or “payoff”) of each genotype as they engage in a specific metabolite-mediated interaction. These payoffs form the payoff matrix of the game. Based on this payoff matrix, we identify all pure strategy Nash equilibria of the game, using an automated pipeline (NashEq Finder, see “Methods”). The payoff matrix also allows one to model evolutionary dynamics (i.e., how genotype frequencies change over time) and to determine which of the identified Nash equilibria are evolutionarily stable (see “Methods”). Supplementary Figs. 11 and 12 provide a more specific representation of this scheme for the presented case studies in this paper

Fig. 2

Metabolic dependencies in populations of S. cerevisiae growing on sucrose. a Metabolic interactions between producer (wild-type, WT) and non-producer (mutant, MT) genotypes of S. cerevisiae growing on sucrose. Here, e represents the percentage of glucose/fructose that diffuses away and serves as a public good. b Nash equilibria and c the equilibrium frequency of WT for the community shown in a as a function of the capture efficiency of the glucose/fructose and the invertase production cost (that latter was implicitly modeled by changing the stoichiometric coefficient of ATP in the sucrose hydrolysis reaction, indicated by x). An alternative in silico formulation of the energetic cost of invertase production that reproduces exactly the setup used in the experiment by Gore et al. (based on histidine auxotrophy) proved to be qualitatively equivalent to the analysis presented here (see Supplementary Note 2 for details). The equilibrium frequency of WT in c was obtained from in silico invasion experiments (see “Methods”) for two cases of a small fraction of MT invading a resident population of WT and vice versa. This analysis demonstrated that the equilibrium frequency of WT is the same in both cases (results are shown here for only one case). d Metabolic interactions between WT and MT when additional glucose is provided in the growth medium (see Supplementary Methods for details of implementation). e Nash equilibria and f the equilibrium frequency of WT in the presence of glucose in the growth medium. The entire Snowdrift game region and part of the Mutually Beneficial region in b are replaced by the Prisoner’s Dilemma game, in e, which is consistent with previous reports and serves as an additional verification of our modeling approach. This is because in the presence of an external supply of glucose, MT is less dependent on WT, leading to an increase in the average fitness of MT

Fig. 3

Metabolic dependencies in populations of E. coli secreting one amino acid. a Possible genotypes in populations of E. coli leaking an amino acid include a prototrophic wild-type strain (WT) self-synthesizing a leaky amino acid and a mutant strain (MT) lacking the gene(s) for the biosynthesis of this amino acid. b The identified Nash equilibria for various leakiness levels (as a percentage of an in silico determined maximum: see Supplementary Methods) across all 20 amino acids. Amino acids are shown here by using their standard three-letter code in the order of increasing in silico growth cost (see also Supplementary Fig. 1). c Experimentally reported leakiness levels of amino acids averaged over three different data sets^,. Values in each data set were normalized to their maximum (see Supplementary Data 2 for values of data). Error bars show standard deviation over the three data sets. d The equilibrium frequency of WT as a function of the leakiness level and amino acid type. In silico invasion experiments for two cases of MT invading WT and WT invading MT revealed that the equilibrium frequencies are insensitive to the initial frequencies (results are shown here for only one case). e Predicted selection coefficients across all amino acids and leakiness levels vs. the experimentally reported ones for E. coli . Empty circles and error bars show, respectively, the average and the range (i.e., minimum and maximum) of the computed selection coefficients across all amino acids and leakiness levels. The blue line and the shaded region around it show the line fitted to experimentally measured selection coefficients and their range (from ref. , amino acid-deficient regime), respectively

Fig. 4

Equilibrium metabolic dependencies in populations of E. coli with two leaky amino acids. a Genotypes involved include a prototrophic strain self-synthesizing two leaky amino acids (i.e., 11), two single-mutant strains each is auxotrophic for one amino acid but synthesizing and leaking the other (i.e., 01 and 10), and a mutant strain auxotrophic for both leaky amino acids (i.e., 00). Here, “1” and “0” denote the presence or absence of biosynthesis pathways (genes) for an amino acid, respectively. b The identified Nash equilibria of two-player games (i.e., pairwise interactions) for all amino acid pairs across different leakiness levels, zoomed in for two selected pairs including c lysine and isoleucine, d glutamate and leucine. A sample payoff matrix of the game is shown in c for a leakiness level of 15% for both lysine and isoleucine. Non-viable equilibria signify associations between genotypes that are non-viable leading to community collapse. Nash equilibria of three- and four-player games for a selected number of amino acid pairs are also given in Supplementary Figs. 2–7

Fig. 5

Metabolic pathway interactions can shape inter-species metabolic dependencies. The identified Nash equilibria for a arginine and glutamate and b glycine and threonine. Sample payoff matrices are given for the region of sustainable leakiness levels by 11 genotypes. Unidirectional dependency is the only possible equilibrium in this region for (arginine, glutamate) while cross-feeding is the only possible equilibrium for glycine and threonine. Metabolic maps show the metabolic pathways involved in the biosynthesis of the respective amino acids

Fig. 6

Impact of the initial genotype frequencies on the evolutionary emergence of metabolic dependencies in populations of E. coli with two leaky amino acids. Here, we have shown the results of targeted in silico invasion experiments for a representative amino acid pair (lysine and isoleucine) (see also Fig. 4c). a 00, 01, and 10 simultaneously originate from 11 (i.e., WT) through genome streamlining and invade an existing population of 11 genotypes. b A small population of 10 and 00 invade a resident population of the 11 and 01. This simulates the second step of a two-step process for the loss of the leaky functions hypothesized in ref. (see the main text for details and also Supplementary Fig. 8). c A small population of 01 and 10 invades a resident population of 11. This models an alternative scenario for the two-step loss of leaky functions leading to stable cross-feeders: Two partial producer mutant genotypes (01 and 10) simultaneously originate from 11, followed by the rise of the 00 genotype from 01 and/or 10 in a later stage. As shown here, cross-feeders can evolutionarily stabilize and coexist with 11 genotypes in the first step. Further analysis showed that cross-feeders are also resistant to invasion by 00 genotypes arising in the second step (see Supplementary Fig. 9). Dynamic plots in a–c show the sample evolutionary dynamics of the system for selected equal leakiness levels for both lysine and isoleucine. Pie charts show the equilibrium frequencies of each genotype starting from the initial genotype frequencies shown in each panel. These equilibrium frequencies are given only for the sustainable leakiness region (green region in Fig. 4c)