Evolution of diversity in metabolic strategies
Abstract
Understanding the origin and maintenance of biodiversity is a fundamental problem. Many theoretical approaches have been investigating ecological interactions, such as competition, as potential drivers of diversification. Classical consumerresource models predict that the number of coexisting species should not exceed the number of distinct resources, a phenomenon known as the competitive exclusion principle. It has recently been argued that including physiological tradeoffs in consumerresource models can lead to violations of this principle and to ecological coexistence of very high numbers of species. Here, we show that these results crucially depend on the functional form of the tradeoff. We investigate the evolutionary dynamics of resource use constrained by tradeoffs and show that if the tradeoffs are nonlinear, the system either does not diversify or diversifies into a number of coexisting species that do not exceed the number of resources. In particular, very high diversity can only be observed for linear tradeoffs.
Introduction
Life on Earth is spectacularly diverse (May, 1988). For example, one study in the early 2000s found that the number of species of fungi is, by a conservative estimate, ca. 1.5 million (Hawksworth, 2001), which was subsequently revised to be between 2.2 and 3.8 million species (Hawksworth and Lücking, 2017). Microbes are by far the most diverse form of life. They constitute approximately 70–90% of all species (Larsen et al., 2017). Perhaps even more astonishing than the number of species is the fact that all of them came from a single common ancestor (Darwin, 1859; Steel and Penny, 2010; Theobald, 2010). To understand the fundamental mechanisms behind such diversification is one of the most relevant problems addressed by the scientific community (Mayr and Mayr, 1963; Coyne, 1992; Rice and Hostert, 1993; Higashi et al., 1999; Dieckmann and Doebeli, 1999; Gavrilets and Waxman, 2002; de Aguiar et al., 2009; Doebeli, 2011).
Recently, ecological interactions, such as competition, have received a lot of attention as potentially very strong drivers of diversification and speciation. A widely used class of models in which this phenomenon can be observed is based on classical LotkaVolterra competition models, which are augmented by assuming that the carrying capacity is a (typically unimodal) function of a continuous phenotype, and that the strength of competition between two phenotypes is measured by a competition kernel, which is typically assumed to be a (symmetric) function of the distance between the competing phenotypes, with a maximum at distance 0 (so that the strength of competition decreases with increasing phenotypic distance).
These assumptions are biologically plausible, and such models have been widely used to provide insights into evolutionary diversification due to competition (Dieckmann and Doebeli, 1999; Doebeli and Ispolatov, 2010; Doebeli and Ispolatov, 2017). However, these models are not derived mechanistically from underlying resource dynamics, and in fact it is known that the commonly used Gaussian functions for the carrying capacity and the competition kernel are not compatible with resourceconsumer models (Abrams, 1986; Ackermann and Doebeli, 2004). A more mechanistic approach is desirable.
Recently, a MacArthur consumerresource model (Macarthur and Levins, 1967) was studied in an ecological context with a view toward explaining the existence of very high levels of diversity (Posfai et al., 2017; Erez et al., 2020). The authors consider different species competing for $p$ interchangeable resources, each supplied at a constant rate (Posfai et al., 2017) or periodically repleted after being used (Erez et al., 2020). A consumer species is characterized by an uptake strategy, $\alpha =({\alpha}_{1},\mathrm{\dots},{\alpha}_{p})$, where the $j$ th component ${\alpha}_{j}\ge 0$ represents the amount of cellular metabolism allocated to the uptake of the $j$ th resource. The rate of consumption of the $j$ th resource and thus its contribution to the growth rate is assumed to be proportional to ${\alpha}_{j}$. The total amount of cellular metabolism available for resource uptake is limited, and hence it is natural to assume a tradeoff between the uptake rates of different resources. In general mathematical terms, a tradeoff is typically given by a function $T(\alpha )=T({\alpha}_{1},\mathrm{\dots},{\alpha}_{p})$ that is increasing in each of the arguments ${\alpha}_{j}$, and such that the only permissible allocation strategies α are those satisfying $T(\alpha )\le E$, where $E$ is a constant. The analysis is then typically restricted to the subspace of strategies defined by $T(\alpha )=E$ (because $T$ is increasing in each ${\alpha}_{j}$). It was shown in Posfai et al., 2017; Erez et al., 2020 that, under the assumption of a linear tradeoff, ${\sum}_{j}^{p}{\alpha}_{j}=E$, very high levels of diversity, that is, many different species with different αstrategies, can coexist. This is a very interesting finding because it violates the competitive exclusion principle (Hardin, 1960), according to which at most $p$ different species should be able to stably coexist on $p$ different resources. Such high levels of diversity emerging from simple consumerresource models could help solve the paradox of the plankton (Hutchinson, 1961) from an ecological perspective.
However, metabolic tradeoffs are not necessarily linear, and in fact there is reason to believe that they almost never are. Nature owes its complexity and diversity to the nonlinearity of the underlying physical and chemical processes. In particular, the nonlinearity of tradeoffs is an essentially inevitable consequence of the general nonlinearity of chemical kinetics. The rate and mass action equilibrium of even a simple bimolecular reaction are in general nonlinear functions of the concentrations of reactants. Linear approximations are commonly used when the concentrations of certain reactants are vastly exceeding the concentrations of others, or when the binding is so strong that the dissociation constant of a complex is much less than typical concentrations of its constituents. However, while the concentrations of enzymes in bacteria (which are probably the most realistic prototype for models of Posfai et al., 2017; Erez et al., 2020) are generally below those of their substrates, the difference is often only few or 10fold, which is insufficient to approximate the enzymatic kinetics by functions that are linear in enzymatic concentrations. For example, a detailed study (Bennett et al., 2009) of the model microbe Escherichia coli revealed that out of 103 metabolites, 35 have concentrations above 1 mM, but the concentrations of 46 metabolites are in tens or single micromole digits, including two metabolites with concentrations below 1 μM. Supporting this, BIONUMBERS (Milo et al., 2010) estimate the typical metabolite concentration in an E. coli bacterium as 32 μM. At the same time, BIONUMBERS provide the evidence for concentrations of important E. coli glycolysis enzymes in tens and even hundreds of μM, and hence the difference between metabolite and enzyme concentrations generally does not seem to be large enough to justify linear approximations.
Another argument for the prevalence of nonlinearity in tradeoffs is based on the oligomerization of more than half of all metabolic enzymes (Marianayagam et al., 2004). The dissociation constants of dimer or oligomer enzymes is often comparable to the concentrations of its monomer units to make the dimerization sensitive to environmental conditions and use it as a regulator of enzymatic activity (Ali and Imperiali, 2005 Traut, 1994). Thus, doubling the concentration of an oligomer requires more (in case of heterooligomer) or less (in case of homooligomer) than doubling the concentrations of its monomers, and hence the metabolic costs of the former in terms of the metabolic costs of the latter are nonlinear.
Since metabolic tradeoffs can often be expected to be nonlinear, here we generalize the models of Posfai et al., 2017; Erez et al., 2020 by incorporating nonlinear tradeoffs in resource use. Specifically, we consider energy budgets of the form
where γ and $E$ are positive constants.
In addition, we incorporate evolutionary dynamics into the ecological models of Posfai et al., 2017; Erez et al., 2020, which allows us to investigate not only the conditions under which diversity can be maintained, but also the evolution of diversity from a single ancestral species. We show that in the resulting evolutionary model, coexistence of more than $p$ species only emerges for the (structurally unstable) linear case $\gamma =1$. Using adaptive dynamics and numerical simulations, we show that regardless of the value of γ, an initially monomorphic population always evolves to an attractive fixed point (also called ‘singular point’), after which two generic scenarios are possible: (i) if $\gamma <1$, the population branches and diversifies, with the maximal number of coexisting species equal to the number of resources $p$, a state in which each species is a complete specialist on exactly one of the resources; (ii) if $\gamma >1$, an initially monomorphic population also evolves to a singular point, but subsequently does not diversify and instead remains a monomorphic generalist.
To make the argument for the relevance of nonlinear tradeoffs even more solid, we prove that an omnipresent nonlinearity in the dependence of nutrient uptake rates on α can be transformed into the nonlinearity of tradeoff (Equation 1), and vice versa. Thus, a nonlinearity in either the tradeoff or the metabolic rates is sufficient to bring the diversity down to the competitive exclusion limit. We also show that the two scenarios (of either a generalist or $p$ specialists) emerge as a result of purely ecological dynamics in a system initially populated with multiple species with different uptake strategies α that satisfy (Equation 1).
Overall, our results show that very high levels of diversity do not evolve in the consumerresource model considered here in a realistic scenario where tradeoffs in resource preference or the resource uptake rates are nonlinear.
Model and results
We consider a population competing for $p$ substitutable resources in wellmixed environments. A phenotypic species α is characterized by its metabolic allocation strategy $\alpha =({\alpha}_{1},\mathrm{\dots},{\alpha}_{p})$, where ${\alpha}_{j}$ is the per capita rate at which individuals of species α take up the $j$ th nutrient. Various coexisting species are distinguished by their specific α’s. From a physiological perspective, ${\alpha}_{j}$ is proportional to the amount of metabolic effort allocated by the individuals of species α to capture nutrient $j$. Intrinsic limitations on metabolic activities impose a restriction on the total amount of nutrient uptake. For simplicity, we assume that this intrinsic limitation leads to a tradeoff in the components ${\alpha}_{j}$ of the form (Equation 1). (Note that we also assume ${\alpha}_{j}\ge 0$ for all $j$.) Throughout, we will set the scaling parameter $E=1$. (See Appendix 1 for a more general treatment, in which the exponent γ can differ for different directions ${\alpha}_{j}$ in phenotype space.)
Following Posfai et al., 2017, we denote by ${c}_{j}(t)$ the concentration of resource $j$ at time $t$, and we assume that the amount of resource $j$ available for uptake per individual (e.g., the amount of resource bound to the outer membrane of a microbial cell) is given by a monotonously increasing function ${r}_{j}({c}_{j})$. Specifically, we assume this function to be of Monod type, ${r}_{j}({c}_{j})={c}_{j}/({K}_{j}+{c}_{j})$. Thus, the rate of uptake of resource $j$ by an individual consumer with uptake strategy α is ${\alpha}_{j}{r}_{j}({c}_{j})$.
Chemostat conditions
We assume that resources are supplied to the system at a constant rate defined by the supply vector $s=({s}_{1},\mathrm{\dots},{s}_{p})$, so that resource $j$ is supplied at a constant total rate s_{j} and decays at a rate ${\mu}_{j}$ (Posfai et al., 2017). This generates the following system of equations for the ecological dynamics of the concentrations c_{j}, $j=1,\mathrm{\dots},p$:
Here, ${n}_{\alpha}(t)$ is the population density of species α at time $t$, so that ${\sum}_{\alpha}{n}_{\alpha}(t){\alpha}_{j}$ is the total amount of metabolic activity invested into uptake of resource $j$ (the sum runs over all species α present in the community). We further assume that the cellular per capita birth rate of species α is equal to the amount of nutrient absorbed by each individual. The dynamics of the population density ${n}_{\alpha}$ then becomes
where δ is the per capita death rate, which is assumed to be the same for all consumers.
The evolutionary dynamics of the the traits ${\alpha}_{j}$ can be solved analytically only for a simplified system in which the resource decay (dilution) rates ${\mu}_{j}$ are set to 0. This assumption, also made in Posfai et al., 2017, corresponds to rapid consumption of almost all resource. In Appendix 1, we derive the adaptive dynamics for the allocation strategies, that is, for the traits ${\alpha}_{j}$ (Metz et al., 1992; Dieckmann and Law, 1996; Dieckmann and Doebeli, 1999; ; Hui et al., 2018; Doebeli, 2011; Geritz et al., 1997). We show that with vanishing decay rates, there is a unique singular point
Calculations of the Jacobian of the adaptive dynamics (an indicator of convergence stability of a fixed point) and of the Hessian of the invasion fitness function (which distinguishes whether the fixed point is an evolutionary endpoint or a branching point) yield the following conclusions: Regardless of the value of γ, the singular point ${\alpha}^{*}$ is always convergent stable, so that the system approaches ${\alpha}^{*}$ from any initial condition. If $\gamma >1$, the singular point ${\alpha}^{*}$ is also evolutionarily stable and hence represents the evolutionary endpoint. In particular, no diversification takes place. On the other hand, if $\gamma <1$, the singular point is evolutionarily unstable and hence is an evolutionary branching point. In particular, if $\gamma <1$, the system will diversify into a number of coexisting consumer species. If $\gamma =1$ (linear tradeoff), the fitness Hessian is 0, representing evolutionary neutrality.
To check our analytical approximations and to investigate the details of diversification after convergence to the evolutionary branching point, we performed numerical simulations of evolving populations consisting of multiple phenotypic strains. The simulations were performed without the simplifying assumption of zero resource degradation (dilution) rates; further details of the numerical simulations are presented in Appendix 1.
In the figures below we show evolving populations as circles with radii proportional to the square root of population size ${n}_{\alpha}$ in threedimensional strategy space $({\alpha}_{1},{\alpha}_{2},{\alpha}_{3})$, viewed orthogonally to the simplex plane ${\sum}_{i=1}^{3}{\alpha}_{i}=1$. With the constraint ${\sum}_{i=1}^{3}{\alpha}_{i}^{\gamma}=1$, the coordinates of each population are $({\alpha}_{1}^{\gamma},{\alpha}_{2}^{\gamma},{\alpha}_{3}^{\gamma})$. In the following numerical examples, we considered a symmetric supply of resources ${s}_{i}=1$ and a slow resource degradation, ${\mu}_{i}{K}_{i}=0.1$.
We first consider scenarios with linear tradeoffs, $\gamma =1$. Figure 1 shows the evolution of a population (shown in blue circles) whose individuals die at constant rate $\delta =1$ (corresponding videos of the simulations can be accessed through the links provided in the figure legends). The black circle represents the singular point that is calculated in the limit of low degradation of nutrients, given by Equation 4. Figure 1(a) shows the initial monomorphic population far from the singular point. An intermediate time of the evolutionary process is shown in Figure 1(b), in which the population remains monomorphic and is approaching the singular point ${\alpha}^{*}$. For $\gamma =1$, the singular point is neutral evolutionarily (all eigenvalues of the Hessian of the invasion fitness function are 0 due to the linearity of the tradeoff), and once the population converges to the singular point, it starts to diversify ‘diffusively’, as anticipated in Posfai et al., 2017: neutrality of selection results in communities consisting of a large number of species. Thus, the high diversity observed in this case is an evolutionary consequence of the selective neutrality caused by a linear enzymatic tradeoff.
The situation changes for nonlinear tradeoffs, $\gamma \ne 1$, which generates two very different evolutionary regimes depending on whether $\gamma >1$ or $\gamma <1$ (even when the deviation of γ from one is small). Figure 2(a–c) shows an example of the evolutionary dynamics for $\gamma =1.1$.
The dynamics starts with an initial monomorphic population far from the singular point, as shown in Figure 2(a). As in the linear case, and as predicted by the analytical theory, the monomorphic population converges toward the singular point Figure 2(b). However, because $\gamma >1$ the singular point is evolutionarily stable, and no diversification occurs (apart from mutationselection balance around the singular point). Instead, when the population reaches the singular point, evolution comes to a halt, and all individuals are generalists, that is, use all resources to some extent (as determined by the location of the singular point), as depicted in Figure 2(c).
On the other hand, Figure 3(a–c) shows the evolutionary process for a community with $\gamma =0.9$. The initial configuration is shown in Figure 3(a). As in the previous examples, the initial phase of evolution ends with the population converging to the singular point ${\alpha}^{*}$. However, in this case, the singular point is an evolutionary branching point giving rise to the emergence of distinct and diverging phenotypic clusters (Figure 3(b)). The final state of the evolutionary process is shown in Figure 3(c): there are three coexisting phenotypic clusters, each being a specialist in exactly one of the resources. Our numerical simulations indicate that the results shown in Figures 1–3 are general and robust: nonneutral diversification occurs only for $\gamma <1$ and typically leads to coexistence of $p$ specialists. In fact, the results easily generalize to situations in which the exponent γ in the tradeoff function may be different for different directions in phenotype space, that is, for different ${\alpha}_{j}$. As we show in Appendix 1, evolutionary branching along a direction ${\alpha}_{j}$ in phenotype space can occur if the corresponding exponent ${\gamma}_{j}<1$. Appendix 1—figure 2 and Appendix 1—figure 3 in Appendix 1 illustrate scenarios in which only a subset of the phenotypic directions ${\alpha}_{j}$ are branching directions along which evolutionary diversification occurs. In such a case, the number of distinct species resulting from the evolutionary process is less than $p$.
Finally, we note that our results for the effects of nonlinear tradeoffs on evolutionary dynamics have corresponding results in purely ecological scenarios, such as those studied in Posfai et al., 2017. We simulated ecological time scales by seeding the system with a set of for example randomly chosen phenotypes throughout phenotype space and running the population dynamics with the mutational process turned off. Again, as shown in Appendix 1—figure 4, nonlinear tradeoffs have a profound effect on the number of surviving species in such ecological simulations, with many species coexisting when $\gamma =1$, as reported in Posfai et al., 2017, but with typically only $p$ species surviving when $\gamma <1$ and only very few species surviving in the close vicinity of the singular point when $\gamma >1$.
Serial dilution conditions
Serial dilution conditions are defined as a sequence of explicitly nonstationary inoculation and growth events (Erez et al., 2020), which mimics seasonality or batch culture experiments (e.g., Lenski and Travisano, 1994). Each growth phase starts with the introduction of a diluted collection of species from a previous batch
into a fresh batch of resources with a given composition ${c}_{j}(0)$. In each batch, the species densities ${n}_{\alpha}(0)$ increase with time as
while resources are depleted:
Unlike in the chemostat model, the death of individuals and the decay of resources are ignored ($\delta =0$ and $\mu =0$). Each event ends at time ${t}_{fin}$ when all resources are almost completely depleted,
and the process is repeated.
Due to the explicit nonstationarity of such serial dilution processes, one of the main assumptions of our adaptive dynamics analysis, the stationarity of resident populations, is not satisfied. Nevertheless, our numerical simulations show that the conclusions drawn for the chemostat case also hold for the serial dilution conditions, to the point that the simulation snapshots are visually indistinguishable from those shown in Figure 2 and Figure 3. However, in the videos, which can be found here , it is possible to see the oscillating population density, caused by the serial dilution protocol.
Specifically, we simulated the serial dilution for three limits considered in Erez et al., 2020, ${c}_{j}(0)=10K$, ${c}_{j}(0)=K$, and ${c}_{j}(0)=0.1K$ for ${\rho}_{0}={10}^{3}$ and ${c}_{fin}={10}^{8}$. All other parameters were the same as used in Figure 1, Figure 2, Figure 3 and corresponding videos.
In all three cases ${c}_{j}(0)\gg K$, ${c}_{j}(0)\sim K$, and ${c}_{j}(0)\ll K$, we observed that for $\gamma >1$, the monomorphic population converges toward the singular point ${\alpha}^{*}$ (Figure 2(b)) and video files here . The singular point is evolutionarily stable, hence, as shown in Figure 2(c), no subsequent diversification occurs (apart from narrow mutationselection spreading around the singular point).
On the contrary, Figure 3(a–c) and videos accessible here show the evolutionary process for a community with $\gamma <1$. The initial configuration is shown in Figure 3(a). As in the previous examples, in the initial phase the monomorphic population evolves close to the singular point ${\alpha}^{*}$. However, in this case, the singular point is again an evolutionary branching point giving rise to the emergence of distinct and diverging phenotypic clusters (Figure 3(b)). The final state of the evolutionary process is shown in Figure 3(c): there are three coexisting phenotypic clusters, each being a specialist on one of the resources.
In addition, purely ecological (i.e., mutationless) simulations performed similarly to what is described above and in Erez et al., 2020 resulted in similar outcomes as in the chemostat model. In a system initially filled with many (200) species, only a few species survive after a fairly short transitory time. When $\gamma >1$, one or a few species remain very close to the singular point ${\alpha}^{*}$, while for $\gamma <1$, typically $p$ specialist species remain in the system. The videos of pure ecological simulations can be seen here.
Once evolution has come to its steady state, resulting in a single generalist species when $\gamma >1$ or $p$ specialist species when $\gamma <1$, each species is represented by a ‘cloud’ of phenotypes α (Panel C in Figure 2 and Figure 3). Such a cloud is formed by a competition between the deterministic selection gradient that acts toward the center of the cloud, making the survival of peripheral species less likely, and the stochastic mutational process that broadens the distribution of strains in all directions, ‘reseeding’ new strains everywhere in the cloud, including its periphery. This is analogous to the classical mutationselection balance occurring with stabilizing selection. When stabilizing selection is relatively weak, which occurs when the tradeoff is only weakly nonlinear (with γ close to 1), the dispersion of phenotypes around the centers of clouds is larger. Technically, this can be concluded from the factor $1\gamma $ in Equation A20. Thus, we make a potentially testable predictions that a weaker nonlinearity in tradeoffs or uptake rates should result in broader distributions of corresponding phenotypes within specialist or generalist species. At the same time, the perspective of mutationselection balance makes it easier to see the difference between the neutral evolutionary scenario of linear tradeoffs and the weakly nonlinear case: While in the former case the distribution of strains will be uniform across the simplex (constrained only by the ‘convex envelope’ condition; Posfai et al., 2017), the nonlinear tradeoffs lead to distinct species with welllocalized distributions of phenotypes for any $\gamma \ne 1$. In an analogy with critical phenomena in physics, correlations typically decay exponentially, except at critical points, where they are exceptionally longranged. Such an ‘anomalous’ behavior requires careful tuning of parameters to get exactly to the critical point, unless the system is ‘selforganized critical’. A similar situation appears to be the case with linear tradeoffs. It could be possible in principle that a system possesses carefully adjusted metabolic parameters so that for a range of uptake rates, the tradeoffs in enzyme concentrations are linear. Yet there appears to be no evolutionary reason for ‘selforganization’ to such a state, and the accidental cancelation of all nonlinearities is very unlikely.
Discussion
To understand the origin and maintenance of diversity is a fundamental question in science. In particular, the mechanisms of diversification due to ecological interactions still generate lively debates.
Recently, tradeoffs in the rates of uptake of different resources were suggested as a mechanism to generate large amounts of diversity (Posfai et al., 2017; Erez et al., 2020), possibly solving the ‘paradox of the plankton’ (Hutchinson, 1961), and violating the competitive exclusion principle (Hardin, 1960), which states that the number of coexisting species should not exceed the number of resources. It has been shown that enzymatic allocation strategies that are plastic instead of fixed, so that individuals can change their allocation (while maintaining a linear tradeoff under a fixed allocation budget) in response to resource availability during their lifetime, tend to reduce the amount of diversity maintained in the ecological communities (PaccianiMori et al., 2020). Perhaps this is not surprising, since more plastic strategies tend to be able to be more generalist as well. As in Posfai et al., 2017; Erez et al., 2020, here, we consider the case of nonplastic strategies, in which each individual is defined by its allocation vector α, but assuming a more general, nonlinear form of tradeoffs. Moreover, we investigate evolutionary rather than just ecological dynamics to determine the conditions under which evolutionary diversification can occur. There are no true jacksofall trades in biology and tradeoffs are a ubiquitous assumption in evolutionary thinking and modeling. However, the cellular and physiological mechanisms that underly such tradeoffs are typically very complicated and the result of biochemical interactions between many different metabolic pathways. Attempts have been made to understand tradeoffs more mechanistically, particularly in microbes (Litchman et al., 2015), but higherlevel modeling efforts most often still require a mostly phenomenological approach to incorporating tradeoffs. In this paper we assumed that each of $p$ resources is available to each microbial organism at a certain rate that depends on the resource concentration in the system. The microbe in turn is described phenotypically by the metabolic allocation strategy that defines its uptake of the available resources.
Without tradeoffs, and everything else being equal, the best strategy would be to allocate an infinite amount (or at least the maximal amount possible) of metabolic activity to every resource, a scenario that is generally unrealistic biologically. Rather, tradeoffs inherent to cell metabolism prevent such strategies. Formally, tradeoffs are given by one or more equations (or more generally inequalities) that the phenotypes of individuals have to satisfy.
In our simplistic models, tradeoffs are determined by the parameter γ, which essentially describes the curvature of the tradeoff function, with the linear tradeoff $\gamma =1$ being the threshold between concave ($\gamma <1$) and convex ($\gamma >1$) tradeoffs. Formally, linear tradeoffs are the simplest case, but there is no a priori general reason why tradeoffs should be linear. Our results show that generically, diversity only evolves with concave tradeoffs, and the number of coexisting species never exceeds the number of resources. Only in the structurally unstable linear case ($\gamma =1$), it is possible for very high levels of diversity to evolve due to the cessation of selection at the evolutionary equilibrium. Any value of $\gamma \ne 1$ precludes high amounts of diversity. Extensive numerical explorations revealed that these results are robust and qualitatively independent of particular parameter choices, such as the number of resources or the dynamics of resource input.
Furthermore, in Appendix 1 we show that the originally nonlinear tradeoffs can be made linear by redefining uptake rates ${\alpha}_{i}$ (Equation A8), thus ‘transferring’ the nonlinearity to the nutrient uptake and the birth rate functions (Equation A8). But a metabolic and nutrient uptake rate is itself a linear function in the enzyme concentration only when the concentration of the substrate vastly exceeds the enzyme concentration. A good example is the wellknown MichaelisMenten approximation, which is identical to the formula used in Posfai et al., 2017; Erez et al., 2020 for the dependence of nutrient uptake on enzyme allocation α. While such linear approximations have been successfully applied in chemical kinetics for over a century, often without questioning their formal validity, the effect of linearization on ecological and evolutionary properties turns out to be very significant. The MichaelisMenten kinetics is valid when the formation of enzymesubstrate complexes does not reduce the concentration of free substrate. Yet the intracellular concentration of enzymes in bacteria are often comparable to or are just few or 10fold smaller than those of their substrates (Bennett et al., 2009; Milo et al., 2010). In Appendix 1 we sketch a derivation of kinetics of an enzymatic reaction in the general case assuming the steadiness of the concentration of the enzymesubstrate complex, but without the assumption that the enzyme concentration is negligible compared to that of the substrate. It follows that enzymatic reaction rates are generally sublinear in the concentrations of enzymes, which is intuitively clear from considering the rate saturation in the limit of infinite enzyme concentrations. However, sublinear rates are not the only possible deviation from linearity: the formation of enzyme oligomers (Marianayagam et al., 2004, Traut, 1994; Ali and Imperiali, 2005; Traut, 1994) and spatially organized complexes (Schmitt and An, 2017) are controlled by intrinsically nonlinear (superlinear in case of homooligomers) mass action equilibria, thus making the enzymatic rates generally sigmoid functions (Ricard and Noat, 1986) of the amount of enzyme. Again, it follows that the physiological costs of the production of individual enzymes are typically nonlinear.
There are also more direct ways to demonstrate the ubiquity of nonlinear dependences of metabolic rates or fluxes $f$ on enzyme concentrations α, for which quantities known as reaction elasticities or flux control coefficients are normally defined as doublelogarithmic derivative, $d\mathrm{ln}[(f(\alpha ))]/d\mathrm{ln}(\alpha )$. For example, for a general power law $f(\alpha )\equiv C{\alpha}^{\gamma}$ that we used to define metabolic tradeoffs (or uptake rates, see Appendix 1), the loglog derivative is equal to γ, the nonlinearity parameter. For the tradeoffs used in Posfai et al., 2017; Erez et al., 2020, this derivate is always 1. However, it is not surprising that realistic assessments of such coefficients (e.g. Loder et al., 2016; Giersch, 1995; Sun and Qian, 2002; Saavedra et al., 2005; Rohwer et al., 2000; Rutkis et al., 2013; Schmidt et al., 2016; van der Vlag et al., 1995; and many other references) produce values that rarely come close to 1, and hence that the measured dependencies of metabolic fluxes on enzyme concentrations are significantly nonlinear. For an easier parametrization of these nonlinearities, it was suggested to express rates of complex enzymatic reactions as products of powerlaw functions of concentrations of enzymes and substrates (Savageau, 1969). This idea, originally suggested more than 50 years ago, has since developed a substantial following, which once again indicates the necessity to account for nonlinearity in the kinetics of enzymatic pathways. All this indicates that reaction elasticities and flux control coefficients are typically distinct from one, which is essentially the main raison d’être for those quantities and for the science of metabolic engineering itself.
Whether sub or superlinear, any deviation of the growth rates from the linear form (Equation 3) and Equation A3 results in a revalidation of the competitive exclusion limit, similarly to nonlinearity in tradeoffs. This serves as another indication that linear tradeoffs in metabolic rates is a biologically unrealistic and exceptional case, while generic nonlinearities do not generate high levels of diversity, and instead the outcomes are in line with classical results about the evolution of resource generalists vs. resource specialists (Ma and Levin, 2006).
It is well known that the shape of tradeoff curves is, in general, an important component in adaptive dynamics models (Kisdi, 2006; Kisdi, 2015). In particular, studies of evolution of cooperation (e.g. Damore and Gore, 2012; Archetti and Scheuring, 2012) have stressed that the outcome of evolution is conditional on the curvature of the public good and cost functions and provided numerous biochemical reasons for nonlinearity of metabolic rates in enzyme concentrations. Here, we have shown the importance of the tradeoff curvature for the evolution and maintenance of diversity in a general consumerresource model. Of course, many potentially important ingredients that could yet lead to high or low diversity in these models were not considered in the present work. For example, dynamic and optimal metabolic strategies (PaccianiMori et al., 2020) and crossfeeding have recently been suggested as factors that could potentially enable such diversity (Goyal and Maslov, 2018), while ‘soft constraints’ that allow random deviations of metabolic strategies from the exact tradeoff constraint were reported in Cui et al., 2020 to reduce the diversity even below the competitive exclusion limit. It will be interesting to consider these model extensions with nonlinear tradeoffs.
Furthermore, it is possible that nonequilibrium ecological dynamics can allow for the maintenance of excess diversity. While this is not the case for externally imposed batch culture dynamics, as reported in the present paper, we have recently shown, using a different ecological model (Doebeli et al., 2021), that endogenous nonstationary ‘boombust’ population dynamics can lead to a significant increase in diversity above the saturation limit expected with equilibrium population dynamics. Together with many experimental results reporting nonstationarity and apparent chaoticity of the population dynamics of actual plankton species, this leads to the conjecture that rather than the neutral evolutionary regime predicted in Posfai et al., 2017, nonstationary population dynamics induced by competition and predation (and perhaps external factors) may be more important in explaining high levels of diversity in natural systems.
Appendix 1
Ecological and evolutionary dynamics
We assume as in Posfai et al., 2017 that metabolic reactions occur on a much faster time scale than cellular division, so that resource concentrations are always at their ecological equilibrium values ${c}_{j}^{*}$ determined as solutions of equations:
with $d{c}_{j}/dt=0$. (Note that these equilibrium resource concentrations are determined by the current populations sizes ${n}_{\alpha}(t)$.) In practice, a faster time scale can be achieved by multiplying the righthand side of Equation A1 by a large dimensionless constant. The cellular per capita birth rate ${g}_{\alpha}$ of species α is proportional to the amount of nutrient absorbed by each individual,
The dynamics of the population size ${n}_{\alpha}$ then becomes
To derive the evolutionary dynamics for the allocation strategies, that is, for the trait α, we follow the adaptive dynamics approach, a powerful tool to study gradual evolutionary diversification due to frequencydependent ecological interactions (Metz et al., 1992; Dieckmann and Law, 1996; Doebeli, 2011; Dieckmann and Doebeli, 1999; Doebeli and Dieckmann, 2003; Geritz et al., 1997; Hui et al., 2018). In particular, adaptive dynamics can generate the paradigmatic phenomenon of evolutionary branching (Metz et al., 1992; Geritz et al., 1997; Doebeli and Dieckmann, 2000; Hui et al., 2018; Doebeli, 2011; Dieckmann and Doebeli, 1999), during which a population that evolves in a continuous phenotype space first converges to a fitness minimum (evolutionary branching point) and then splits into two (or more) diverging phenotypic branches. We start with considering a monomorphic resident population at its ecological equilibrium ${n}_{\alpha}^{*}$, which is defined as the population size for which the equilibrium resource levels ${c}_{j}^{*}$ are such that ${g}_{\alpha}({c}_{1}^{*},\mathrm{\dots},{c}_{p}^{*})=\delta $ (note again that the $({c}_{1}^{*},\mathrm{\dots},{c}_{p}^{*})$ implicitly depend on α). The invasion fitness of a rare mutant ${\alpha}^{\prime}$ is then the per capita growth rate of the mutant ${\alpha}^{\prime}$ at the resource levels defined by the resident:
To derive the adaptive dynamics, we consider the selection gradient $q(\alpha )=({q}_{1}(\alpha ),\mathrm{\dots},{q}_{p}(\alpha ))$, with components
$q(\alpha )$ defines a $p$dimensional dynamical system in unrestricted αspace,
The speed of evolution of α is proportional to the current ecological equilibrium population size ${n}_{\alpha}^{*}$ because the number of mutations occurring at any given point in time is proportional to ${n}_{\alpha}^{*}$. The parameter σ describes both the per capita rate and effective size of mutations. Without loss of generality, we set $\sigma =1$.
To take the enzymatic tradeoff into account, the unconstrained adaptive dynamics (Equation A6) needs to be restricted to the surface in αspace that is defined by the tradeoff ${\sum}_{j=1}^{p}{\alpha}_{j}^{\gamma}=E$, where $E$ is a positive number (Ito and Sasaki, 2016). An illustrative example is the one in which the nutrients come from three different resources. The tradeoff ${\alpha}_{1}^{\gamma}+{\alpha}_{2}^{\gamma}+{\alpha}_{3}^{\gamma}=E$ defines a surface in αspace containing all strategies. The curvature of each surface is determined by γ. Appendix 1—figure 1a shows an example of the surface defined by the tradeoff for the case that $\gamma >1$ while Appendix 1—figure 1b and c show the curvature for the case where $\gamma =1$ and $\gamma <1$, respectively. The blue star and the orange diamond illustrate possible position of the strategies in αspace. The individuals with strategy indicated by the blue star uptake nutrients only from resource s_{1}, while the individuals with strategy indicated by the orange diamond uptake nutrients from all three resources.
Equilibrium points of the adaptive dynamics, the socalled singular points, are resting points ${\alpha}^{*}$ of the resulting dynamical system in phenotype space. Given a singular point ${\alpha}^{*}$, two stability concepts are important. First, there is stability in the usual sense of converging to ${\alpha}^{*}$ from nearby initial conditions, which is measured by the Jacobian matrix of the functions defining the adaptive dynamics, evaluated at ${\alpha}^{*}$. Second, evolutionary stability is measured by the Hessian of the invasion fitness function $f({\alpha}^{*},{\alpha}^{\prime})$ with respect to the mutant trait, evaluated at the singular point ${\alpha}^{*}$, and taken along the constraint surface (Ito and Sasaki, 2016). A negative definite Hessian (all eigenvalues negative) means that the singular point is a maximum of invasion fitness and no branching occurs. Alternatively, a singular point is called an evolutionary branching point if it is both convergent stable with regard to the Jacobian and evolutionarily unstable with regard to the Hessian. Thus, a singular point is a branching point if all eigenvalues of the Jacobian have negative real parts, and if the Hessian matrix is not negative definite.
Singular points and their convergence and evolutionary stability
The case where the decay rates, ${\mu}_{j}$, are zero for any $j$ admits an analytical solution. We consider allocation strategies $\alpha =({\alpha}_{1},\mathrm{\dots},{\alpha}_{p})$ as in the main text, but here we assume a more general tradeoff function:
It turns out to be convenient to reparametrize the strategy space as follows:
for $j=1,\mathrm{\dots},p$. This simplifies the tradeoff expression to
Because the ${\beta}_{j}$ increase monotonically with ${\alpha}_{j}$, the adaptive dynamic properties in terms of convergence and evolutionary stability of singular points are the same for $\alpha =({\alpha}_{1},\mathrm{\dots},{\alpha}_{p})$ and $\beta =({\beta}_{1},\mathrm{\dots},{\beta}_{p})$ phenotypes. However, the tradeoff in β, Equation A9, is linear, which simplifies the analysis.
In terms of β, the per capita rate of use of resource $j$ of an individual with phenotype β is
We assume that nutrients are supplied to the system at a constant rate given by the vector $s=({s}_{1},{s}_{2},\mathrm{\dots},{s}_{p})$, where s_{j} is the supply rate of the $j$ th resource. We consider the low degradation rate regime, that is, ${\mu}_{j}\to 0$ for all $j$ in Equation A1. Setting the righthand sides of Equations A1, A3 equal to zero and taking into account that the sum in Equation A1 consists of a single term, we obtain for the equilibrium density of a population monomorphic in β
The invasion fitness of a rare mutant with uptake strategy ${\beta}^{\prime}$ in a resident β at ecological equilibrium ${n}_{\beta}^{*}$ becomes
To derive the adaptive dynamics of β, we calculate the selection gradient $q(\beta )=({q}_{1}(\beta ),\mathrm{\dots},{q}_{p}(\beta ))$ and project it onto the linear constraint space:
Here, the term
is the component of the selection gradient (Equation A13) that is orthogonal to the tradeoff hyperplane (note that $(1/\sqrt{p},\mathrm{\dots},1/\sqrt{p})$ is a unit vector orthogonal to the tradeoff hyperplane).
If we set the mutational parameter $\sigma =1$, the adaptive dynamics of β becomes
Note that we only need $p1$ equations due to the (linear) tradeoff. It is easy to see that Equation A16 has a unique fixed point, that is, there is a unique singular point for the adaptive dynamics given by
for $j=1,\mathrm{\dots},p$. In terms of the original trait α, Equation A17 is (Equation 4) in the main text.
To check for convergence stability of ${\beta}^{*}=({\beta}_{1}^{*},\mathrm{\dots},{\beta}_{p}^{*})$, we have to calculate the Jacobian matrix $J$ of the righthand side of Equation A16, evaluated at the singular point ${\beta}^{*}$. It is easy to see that the $jk$ th element of $J$ is
Thus, $J$ is of the form
where ${J}_{d}$ is a diagonal matrix with element $J}_{jj}=\frac{{s}_{j}}{{\gamma}_{j}{\beta}_{j}^{\ast 2}$ and $A$ is a matrix whose elements in the $k$ th column are all identical and equal to $\frac{1}{p}\frac{{s}_{k}}{{\gamma}_{k}{\beta}_{k}^{*2}}$. This implies that the matrix $A$ maps any vector in phenotype space to a vector that is orthogonal to the tradeoff hyperplane (i.e., to a multiple of the vector $(1,\mathrm{\dots},1)$). If $\mathrm{\Delta}\beta =(\mathrm{\Delta}{\beta}_{1},\mathrm{\dots},\mathrm{\Delta}{\beta}_{p})$ is any vector of deviations from the singular point, it follows that the projection of $J\mathrm{\Delta}\beta $ onto the tradeoff hyperplane is the same as the projection of the vector ${J}_{d}\mathrm{\Delta}\beta $. Since all eigenvalues of ${J}_{d}$ are real and negative, it follows that the singular point ${\beta}^{*}$ is a local attractor, that is, convergent stable, regardless of the exponents ${\gamma}_{j}$, $j=1,\mathrm{\dots},p$.
For evolutionary stability, we have to calculate the Hessian matrix $H$ of second derivatives of the invasion fitness function, Equation A12 with respect to the mutant trait ${\beta}^{\prime}$ and evaluated at the singular trait value ${\beta}^{*}$. The $jk$ th element of $H$ is
where ν is a constant:
Thus, $H$ is diagonal (due to the transformation from α to β), and $H$ is negative definite, that is, all eigenvalues are negative, if and only if ${\gamma}_{j}>1$ for all $j=1,\mathrm{\dots},p$. Because the tradeoff hyperplane is linear in β, it follows that any index $j$ with ${\gamma}_{j}<1$ provides a branching direction ${\beta}_{j}$, that is, a direction in phenotype space along which evolutionary diversification is possible. More precisely, any direction in βspace (other than orthogonal to the tradeoff surface) along which the unconstrained Hessian (Equation A20) has a minimum corresponds to a direction on the tradeoff surface along which diversification is possible.
The results presented in the main text now follow from the above by setting ${\gamma}_{j}=\gamma $ for $j=1,\mathrm{\dots},p$. But the above analysis also suggests that with suitably chosen ${\gamma}_{j}$, it is possible to generate evolutionary branching in some directions, but not in others. This is illustrated in Appendix 1—figure 2. Our analysis and numerical procedure are applicable to an evolving system of populations with any number of resources. To facilitate visualization, in the following, we consider just three resources, so that because of the constraint, each population is characterized by two independent parameters ${\alpha}_{i}$. Appendix 1—figure 2 illustrates diversification in the direction of ${\alpha}_{3}$ (${\gamma}_{3}<1$) without diversification in the directions ${\alpha}_{1}$ and ${\alpha}_{2}$ (${\gamma}_{1},{\gamma}_{2}>1$).
Appendix 1—figure 3 illustrates diversification in the directions ${\alpha}_{1}$ and ${\alpha}_{2}$ (${\gamma}_{1},{\gamma}_{2}<1$), with no diversification in ${\alpha}_{3}$ (${\gamma}_{3}>1$).
Numerical procedures
In the chemostat simulations, we numerically integrate the system of population dynamics Equation A3 for $M$ populations using a simple Euler update ($M=1$ at the beginning of the simulations). After each integration step, the populations that fall below a small ‘extinction’ threshold density (normally ${n}_{min}={10}^{6}$) are removed from the system. The resource concentrations c_{i} and uptake rates r_{i} are considered relaxed to their steady states for a given set of populations $\{{n}_{\alpha}\}$,
where ${\varphi}_{i}={\sum}_{\alpha}{n}_{\alpha}{\alpha}_{i}$.
To simulate serial dilutions, we numerically integrate Equations 6, 7 for $M$ populations and $p$ resources using also the Euler update ($M=1$ at the beginning of the simulations). After each integration step, the populations that fall below a small ‘extinction’ threshold density (normally ${n}_{min}={10}^{6}$) are removed from the system. Once the resources are depleted so that the condition (Equation 8) is satisfied, the populations of all existing species are rescaled according to Equation 5 and the resource concentrations are reset to ${c}_{j}(0)$.
To mimic mutations in both simulation setups, periodically (typically once every $\mathrm{\Delta}{t}_{mutation}=1$ time units) a mutant is split from an ancestor, which is randomly chosen with probability proportional to its total birth rate. The mutant’s phenotype is randomly offset from the ancestral phenotype along the constraint surface. The offset distance is drawn from a uniform distribution in the interval $[m,m]$. Unless otherwise noted, $m=0.005$. The mutant population is set to be 10% of the ancestral one, and the ancestor population is reduced by 10%. In addition to mutations, periodically (typically once every $\mathrm{\Delta}{t}_{merge}=100\mathrm{\Delta}{t}_{mutation}$ time units) populations that are within a distance $m$ of each other are merged (preserving their phenotypic center of mass) and their population sizes added. Periodic repetition of mutation and merging procedures preserves the phenotypic variance necessary for evolution while limiting computational complexity.
This produces clouds or clusters of αvalues in phenotype space (see figures), with each phenotype α representing a monomorphic population of individuals with that phenotype α. Somewhat imprecisely, we refer to a distinct cluster of phenotypes as a species. The clusters move in phenotype space due to extinction and merging of phenotypes, and due to creation of new phenotypes by mutation. This movement represents evolution and occurs along the constraint surface. A diversification event occurs when a cluster corresponding to the diversifying species spontaneously splits into two or more clusters that diverge from each other and move apart.
Maintenance of diversity in ecological time scales
Here, we briefly show how nonlinear tradeoffs affect diversity on ecological time scale. To do this, we initiate the simulations with a set of for example randomly chosen phenotypes throughout phenotype space and then run the systems with the mutational process turned off. In Appendix 1—figure 4a, we show the initial configuration used for three different scenarios with different exponents γ of the tradeoff (here, we again assume that ${\gamma}_{j}=\gamma $ for all $j$). The functional form of the tradeoffs has a profound effect on the number of surviving species, with many species coexisting when $\gamma =1$, as reported in Posfai et al., 2017 (Appendix 1—figure 4b), but with typically only $p$ species surviving when $\gamma <1$ (Appendix 1—figure 4c) and only very few species surviving in the close vicinity of the singular point when $\gamma >1$ (Appendix 1—figure 4d).
Videos of ecological simulation of serial dilution scenario can be found here.
Nonlinear metabolic rates
Consider the consumption or transformation of a resource substance $c$ into a downstream metabolic product $m$ using a specific enzyme α:
The assumption that the concentration of the complex $c\alpha $ instantaneously relaxes to its steady state defined by the current concentrations of $c$ and α constitutes the MichaelisMenten approximation,
Assuming mass action kinetics and denoting by $c$ and α the total (bound in the complex plus free) concentrations of the resource and enzyme, and dropping the traditional symbols for concentrations, one gets a quadratic equation for the concentration ψ of the $c\alpha $ complex:
with the solution
Here, $\chi \equiv ({k}_{1}+{k}_{2})/{k}_{1}$ is the dissociation constant for the complex. The more common form of the MichaelisMenten approximation,
is obtained in the limit when the substrate concentration is much larger than that of the enzyme. (We note that this form is linear in α and is identical to the product of a Monod function and the corresponding α used in Equation 3). Since
the resource uptake and the growth rates are always sublinear in the concentration of the enzyme α.
Data availability
All data generated or analysed during this study are obtained through the codes which have been deposited in https://github.com/jaros007/Codes_for_Evolution_of_diversity_in_metabolic_strategies (copy archived at https://archive.softwareheritage.org/swh:1:rev:d0a9ad7ca4459a1cc221b7bf1d1d311733400f0a).
References

Character displacement and niche shift analyzed using consumerresource models of competitionTheoretical Population Biology 29:107–160.https://doi.org/10.1016/00405809(86)900079

Protein oligomerization: how and whyBioorganic & Medicinal Chemistry 13:5013–5020.https://doi.org/10.1016/j.bmc.2005.05.037

Review: game theory of public goods in oneshot social dilemmas without assortmentJournal of Theoretical Biology 299:9–20.https://doi.org/10.1016/j.jtbi.2011.06.018

Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coliNature Chemical Biology 5:593–599.https://doi.org/10.1038/nchembio.186

Effect of resource dynamics on species packing in diverse ecosystemsPhysical Review Letters 125:048101.https://doi.org/10.1103/PhysRevLett.125.048101

Understanding microbial cooperationJournal of Theoretical Biology 299:31–41.https://doi.org/10.1016/j.jtbi.2011.03.008

BookOn the Origin of Species by Means of Natural SelectionMurray, London: The royal society.

The dynamical theory of coevolution: a derivation from stochastic ecological processesJournal of Mathematical Biology 34:579–612.https://doi.org/10.1007/BF02409751

BookAdaptive Diversification (MPB48)Princeton, United States: Princeton University Press.

Evolutionary branching and sympatric speciation caused by different types of ecological interactionsThe American Naturalist 156:S77–S101.https://doi.org/10.1086/303417

Continuously stable strategies as evolutionary branching pointsJournal of Theoretical Biology 266:529–535.https://doi.org/10.1016/j.jtbi.2010.06.036

Diversity and coevolutionary dynamics in HighDimensional phenotype spacesThe American Naturalist 189:105–120.https://doi.org/10.1086/689891

Dynamics of adaptation and evolutionary branchingPhysical Review Letters 78:2024–2027.https://doi.org/10.1103/PhysRevLett.78.2024

Diversity, stability, and reproducibility in stochastically assembled microbial ecosystemsPhysical Review Letters 120:158102.https://doi.org/10.1103/PhysRevLett.120.158102

The competitive exclusion principleScience 131:1292–1297.https://doi.org/10.1126/science.131.3409.1292

The magnitude of fungal diversity: the 1.5 million species estimate revisitedMycological Research 105:1422–1432.https://doi.org/10.1017/S0953756201004725

Fungal diversity revisited: 2.2 to 3.8 million speciesMicrobiology Spectrum 5:2016.https://doi.org/10.1128/microbiolspec.FUNK00522016

Modelling coevolution in ecological networks with adaptive dynamicsMathematical Methods in the Applied Sciences 41:8407–8422.https://doi.org/10.1002/mma.4612

Evolutionary branching under multidimensional evolutionary constraintsJournal of Theoretical Biology 407:409–428.https://doi.org/10.1016/j.jtbi.2016.07.011

Tradeoff geometries and the adaptive dynamics of two coevolving speciesEvolutionary Ecology Research 8:959–973.

Construction of multiple tradeoffs to obtain arbitrary singularities of adaptive dynamicsJournal of Mathematical Biology 70:1093–1117.https://doi.org/10.1007/s0028501407885

Inordinate fondness multiplied and redistributed: the number of species on earth and the new pie of lifeThe Quarterly Review of Biology 92:229–265.https://doi.org/10.1086/693564

The evolution of resource adaptation: how generalist and specialist consumers evolveBulletin of Mathematical Biology 68:1111–1123.https://doi.org/10.1007/s1153800690966

The limiting similarity, convergence, and divergence of coexisting speciesThe American Naturalist 101:377–385.https://doi.org/10.1086/282505

The power of two: protein dimerization in biologyTrends in Biochemical Sciences 29:618–625.https://doi.org/10.1016/j.tibs.2004.09.006

How many species are there on earth?Science 241:1441–1449.https://doi.org/10.1126/science.241.4872.1441

How should we define 'fitness' for general ecological scenarios?Trends in Ecology & Evolution 7:198–202.https://doi.org/10.1016/01695347(92)90073K

BioNumbersthe database of key numbers in molecular and cell biologyNucleic Acids Research 38:D750–D753.https://doi.org/10.1093/nar/gkp889

Dynamic metabolic adaptation can promote species coexistence in competitive microbial communitiesPLOS Computational Biology 16:e1007896.https://doi.org/10.1371/journal.pcbi.1007896

Metabolic TradeOffs promote diversity in a model ecosystemPhysical Review Letters 118:028103.https://doi.org/10.1103/PhysRevLett.118.028103

Catalytic efficiency, kinetic cooperativity of oligomeric enzymes and evolutionJournal of Theoretical Biology 123:431–451.https://doi.org/10.1016/s00225193(86)802120

Laboratory experiments of speciation: what have we learned in 40 years?Evolution 47:1637–1653.

Understanding glucose transport by the bacterial phosphoenolpyruvate:glycose phosphotransferase system on the basis of kinetic measurements in vitroJournal of Biological Chemistry 275:34909–34921.https://doi.org/10.1074/jbc.M002461200

Biochemical systems analysis. I. some mathematical properties of the rate law for the component enzymatic reactionsJournal of Theoretical Biology 25:365–369.https://doi.org/10.1016/s00225193(69)800263

The quantitative and conditiondependent Escherichia coli proteomeNature Biotechnology 34:104–110.https://doi.org/10.1038/nbt.3418

Flux control analysis for biphenyl metabolism by Rhodococcus pyridinovorans r04Biotechnology Letters 24:1525–1529.

Dissociation of enzyme oligomers: a mechanism for allosteric regulationCritical Reviews in Biochemistry and Molecular Biology 29:125–163.https://doi.org/10.3109/10409239409086799

Control of glucose metabolism by the enzymes of the glucose phosphotransferase system in Salmonella typhimuriumEuropean Journal of Biochemistry 230:170–182.https://doi.org/10.1111/j.14321033.1995.0170i.x
Decision letter

Wenying ShouReviewing Editor; University College London, United Kingdom

Aleksandra M WalczakSenior Editor; École Normale Supérieure, France

Wenying ShouReviewer; University College London, United Kingdom
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
This manuscript generalizes a recentlypublished framework by the Wingreen group which considers adaptation/evolutionary dynamics in a consumer resource model when a tradeoff on metabolic rates is imposed. The original model involved a linear tradeoff: e.g. the faster a microbe can eat nutrient 1, the slower it can eat nutrient 2, so that a sumrule is imposed on the maximal rates. This paper considers the mathematics of the effect of nonlinear tradeoffs, in the form of a sum of rates raised to a power γ, with γ not just equaling 1 (γ=1 represents the regime of linear tradeoff). The authors find that, keeping all other model features the same, diversity is lost when γ is not 1. That is, nonlinear tradeoff destroys diversity.
Decision letter after peer review:
[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]
Thank you for submitting your work entitled "Evolution of diversity in metabolic strategies" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Wenying Shou as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by a Senior Editor .
We are sorry to say that, after consultation with the reviewers, we have decided that your work will not be considered further for publication by eLife.
Although one reviewer was positive, the overwhelming sentiment was that the work does not sufficiently advance the field.Reviewer #1:
Previous theoretical work argued that among species that compete for resources, physiological tradeoff (e.g. consuming more of food 1 leads to consuming less of food 2) can give rise to species diversity that greatly exceeds the number of resources, even in a wellmixed environment not conducive for species diversity. If previous work were to be general, it would be exciting because this offers a clue to the puzzle scientists have been trying to solve for a long time: what supports high species diversity? Caetano et al show that the finding from previous work only holds under a very restrictive condition (tradeoff function being linear). When that condition is violated (which can frequently occur in biology), we end up with either a single generalist species, or specialists each specializing on a single resource. Thus, in general, the total number of species cannot be larger than the total number of limiting resources in a wellmixed environment, as posited by the competitive exclusion principle. In short, we are back at where we were.
The authors started by introducing the adaptive dynamics framework which has been used to study evolutionary diversification due to frequencydependent selection. The introduction is not adequate in giving nonspecialists an intuitive feeling about how evolutionary branching point from fitness minimum works. Of course, making theoretical conclusions accessible is not trivial, and may not be achievable. However, authors can try harder by including supplementary figures.
The authors then introduced work from the Wingreen group: R number of resources in a wellmixed environment, and each consumer species has an uptake strategy for each of the R resources. The total uptake for each species is fixed, creating a tradeoff: more uptake of one resource reduces the uptake of another. Previous work showed that when the tradeoff function is linear (e.g. the sum of [uptake of each resource)^{γ}] = a constant where the exponent γ=1, then many species with different uptake strategies can coexist.
The authors showed that the conclusion from previous work is rather restricted in its scope: large diversity only exists when the exponent γ in the tradeoff function is 1 (i.e. linear). When that exponent is greater than 1, no diversification occurs (saving for mutationselection balance) and all individuals are phenotypically similar generalists. When the exponent is less than 1 (i.e. concave), the initial convergence to an unstable steady state later evolutionarily diverges into specialists, each specializing on a resource.
Authors also tested their conclusions within the ecological framework of the Wingreen group, and in different scenarios (chemostat versus serial dilutions), and reached the same conclusions.
The paper is relatively easy to read (in the realm of theoretical papers). This work reminds me of the work that has been done in the field of the evolution of cooperation. For example, public goods games often assume that the effect of the public good is a linear function of the number of contributions, an assumption that is often violated in biology. Depending on whether this function is linear or nonlinear, one can get very different outcomes in cooperator/noncooperator coexistence (e.g. Archetti and Scheuring, JTB 299:920; Damore and Gore, JTB, https://doi.org/10.1016/j.jtbi.2011.03.008). Authors may want to add a discussion on that.
I wonder whether this paper should be added as a Research Advance to the original paper from the Wingreen paper published in eLife.
Reviewer #2:
This paper deals with the diversification of metabolic strategies in an evolving population. The authors consider a consumerresource model under different metabolic tradeoffs (sublinear, linear, and superlinear). They show that the linear case is a marginal scenario that corresponds to high diversity as a consequence of neutral evolution. Both the sub and superlinear cases lead to the coexistence of fewer species than resources, as expected by the competitive exclusion principle.
The manuscript is well written and easy to follow. The derivation using adaptive dynamics is interesting and the results are robust. I am mainly concerned by the premises of the work.
l 70 "This is a very interesting finding because it violates the competitive exclusion principle". This is not strictly correct (and I think this is a central point). To violate the competitive exclusion principle, more species than resources should *stably* coexist. The stability requirement is essential. Otherwise, the principle can be easily falsified by a neutral model: in presence of even 1 single resource, an arbitrary number of ecologically equivalent species coexist (neutrally). Neutral coexistence is, as well known, structurally unstable: arbitrary small differences (which break the ecological equivalence) drive many species to extinction (restoring the bound on diversity given by the number of resources).
In the model by Posfai et al. (2017) coexistence is in fact only neutral. There is a manifold of fixed points and stability is marginal (several eigenvalues of the community matrix are equal to zero, e.g. see https://arxiv.org/pdf/2002.04358 ). The fixed point of their dynamics (abundance of different species) depends on the initial conditions. The high levels of diversity are, as a consequence, structurally unstable. This can be shown in multiple ways: introducing an (arbitrarily small) species variability in the tradeoff (E depends on species identity), introducing variability in the dilution rate d (appearing in eq 3), or, as done in the paper, by altering the functional form of the tradeoff.
This is a central point. It explains why many species are observed in the model by Posfai et al. And it explains why the result is extremely (infinitely) sensitive to the parameterization. These ecological considerations are mirrored in the fact that for gamma = 1 the evolutionary dynamics are neutral.
My main concern about this work is whether it has a sufficient degree of novelty and interest. As mentioned in the public review, the results are robust. But, after a close analysis of the model by Posfai et al., they are not unexpected. The manuscript, as it stands, mostly demonstrates the weaknesses of the paper by Posfai et al.: fixed points with more species than resources exist, but they are only marginally stable (https://arxiv.org/pdf/2002.04358 ). Structural instability is a direct consequence of this fact. Nonlinear tradeoffs are just one (among many) ways to show that the results are infinitely sensitive to the parameter choices.
Reviewer #3:
There has been considerable recent interest in understanding the high degree of diversity observed in microbial communities. From a theoretical perspective, this has led to a resurgence of interest in resourcecompetition models. Several recent papers have studied the effects of tradeoffs on total enzyme budgets within these models. An interesting observation is that with exact tradeoffs, communities can selforganize to a state with an arbitrarily large number of species coexisting. One assumption of these models is that the total "cost" of enzymes is a linear function. The current work relaxes this assumption, and shows that this state of arbitrarily high coexistence relies on the linearity of the cost function.
Strengths: This study is rigorous, clearly presented, and the conclusions are mathematically sound. The authors analyze both chemostat and serial dilution systems.
Weaknesses: The results are qualitatively as expected from previous studies of the role of inexact tradeoffs, and are more limited. The nonlinear tradeoffs explored here are essentially equivalent to the unequal enzyme budgets explored in prior work. Indeed, these nonlinearities can be directly mapped to unequal budgets: for example, a nonlinearity that favors expression of a single enzyme is directly equivalent to a larger enzyme budget for species that produce only a single enzyme. Previous studies showed that unequal enzyme budgets lead to a loss of diversity, as is found in this work. Moreover, these prior studies found that even if tradeoffs are not exact, the slow loss of diversity due to inexact tradeoffs can be offset by invasion of new strategies and can therefore still lead to a large number of coexisting species.
The likely impact of this work on the field is modest, given that those who are already experts in the field will recognize that nonlinear tradeoffs are equivalent to unequal enzyme budgets. Moreover, the current study does not actually provide any specific support for nonlinear tradeoffs other than a few remarks in the Introduction.
The work would be of greater general interest if the biological evidence for nonlinearities in enzyme costs were carefully examined; mechanistic insights on how enzyme budget nonlinearities may arise in nature would be of significant utility to the field. However, this would require a substantial additional undertaking and would shift the focus of the work from the specific implications of nonlinearities in resourcecompetition models. An alternative would be to publish the current study in a more specialized journal, with a more theoretical focus.
It would also be helpful to provide a quantitative assessment of the sensitivity of diversity to the degree of nonlinearity. It is clear that any nonlinearity (or inexactness of tradeoffs) leads to loss of diversity at long times. However, a small rate of invasion by new strategies can still lead to a diverse stationary state of the population. Given a certain degree of nonlinearity, how much invasion is required to maintain diversity? The adaptive dynamics calculations performed by the authors do not address this point because new strategies are only introduced if these are more fit than the residents. The question of diversity requires introducing invaders that may be slightly less fit, but still manage to survive due to demographic noise.
[Editors’ note: further revisions were suggested prior to acceptance, as described below.]
Thank you for submitting your article "Evolution of diversity in metabolic strategies" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Wenying Shou as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Aleksandra Walczak as the Senior Editor.
The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.
Essential Revisions:
Please consult the Reviewers' comments and address these in your revision.
Reviewer #1:
I remain supportive, especially if authors can discuss empirical measurements of biological tradeoffs and whether in the natural environment the linearity assumption might break down.
Reviewer #3:
The authors have made several cosmetic changes which have improved the clarity of their manuscript. However, the revised version does not substantially address my main concerns. The real question is whether the current manuscript makes a substantial contribution to the topic of microbial diversity. The focus of the paper is a critique of a model of resource competition with tradeoffs. It is certainly legitimate to be critical of existing models. However, I believe the readers of eLife already appreciate the adage "all models are wrong, but some are useful". The authors have focused their attention on the first part of the adage, arguing that because growth functions will not be exactly linear the model is "wrong". But it's not news that the model is "wrong" (see above), the question is whether the model might still be a useful starting point for understanding diversity? What seems to me to be missing in the discussion, both in the original studies of Ref. 3 and 5 and the current manuscript, is quantification of how "wrong" the initial model is, and whether this undermines its utility. This is why I suggested that the authors carefully examine the "biological evidence for nonlinearities in enzyme costs". Their revised manuscript adds some sentences on this point, but in a nonquantitative way: the authors continue to make the mathematical point that the model is "wrong", but have not taken up the challenge of addressing whether it is or is not "useful". Yes, it is mathematically correct as the authors state that bimolecular reactions are strictly nonlinear in the reactants. But for typical enzyme concentrations in the μM range and typical metabolites in the 0.11 millimolar range, these nonlinearities are in the 0.11% range. From a biological perspective, a linear function might therefore still be a useful starting point. I've also read the references the authors cite on other sources of nonlinearity – they are equally nonquantitative. For example, the review by Marianayagam et al. states (without citations) "In its simplest form, oligomerization functions as a general mechanism for sensing protein concentration. An increase in protein concentration above the oligomerization threshold can be the stimulus for enzyme activation; similarly, enzyme deactivation will apply when cellular levels of the enzyme fall." I absolutely agree that for enzymes that need to oligomerize to function, this implies a mathematically nonlinear processing rate as a function of enzyme concentration. However, again for enzyme levels in the μM range and oligomerization dissociation constants in the commonly observed 110 picomolar range, the nonlinearities are again ~0.11%. Despite the revisions on this and other points in the reviews, in the end I am left still wondering whether the original model is "useful". My conclusion is that the current manuscript will be primarily of interest to researchers whose focus is on the mathematics of resourcecompetition models, and would therefore be appropriate for a more mathematically focused journal.
Reviewer #4:
This is a wellreasoned and wellpresented paper, which considers adaptation/evolutionary dynamics in a consumer resource model when a tradeoff on metabolic rates is imposed. This finding will be of interest to the ecology community, and is well presented. Moreover, these consumerresource models have relevance beyond the microbial ecology scenarios outlined here with applications in biology and beyond. Though the mathematics of their specific suggested model is wellanalyzed, it stands to reason that when considering realistic systems, other timescales can be involved. As a result, the interplay between those timescales and other model parameters could quite possibly lead again to stable coexistence, as in the case of gamma=1.
The main caveat of this paper is that one does not need to resort to nonlinear tradeoffs to destroy diversity in the original Wingreen models, since simply giving different organisms different enzyme budgets would suffice  a plausible enough scenario. In the case of unequal enzyme budgets, by adding other (plausible) diversity stabilizing mechanisms it is possible to retain the coexistence state. Similarly, for values gamma sufficiently close to 1, I would expect that other realistic effects could again lead to a state resembling gamma=1. I would have liked to see this paper attempt to understand how diversity is maintained away from gamma=1, rather than just show us a way to break it.
The reality of the field (if I may) is that none of us have the answer to the socalled Paradox of the Plankton. There are several competing theories out there, and each of them has its merits and imperfections. Simply outlining an expected and obvious lack of generalizability of a particular aspect of an otherwise solid and informative theory is less constructive than this paper could be. Especially considering it is well known that consumerresource models (e.g. Wingreen) support coexistence only on equal enzyme budgets unless other diversity stabilizing mechanisms are included.
What I think that the authors should focus on, is that neither the Wingreen papers nor recent papers from other groups claim that their particular framework captures all features of diversity maintenance. Instead, we are all trying, as a community, to piece together a coherent theory, stumbling along as we inevitably do. For example, the Wingreen papers consider violation of the linear constraint (i.e. unequal enzyme budgets) in several stochastic and deterministic cases, outlining their relevance. There are several timescales at play here and one can think how their interplay with particular model parameters influences the resulting steady state. This is likely the case with regards to this manuscript. Already, if one considers small deviations from γ = 1, one introduces a new timescale, to compete with other model timescales. Moreover, the "mutation" dynamics (Eq. 14) introduce one or more timescales, which are conveniently swept under the rug with a "without loss of generality" (l463), a statement without actual backing as far as I understood. Furthermore, the authors do not sufficiently acknowledge that their suggested "mutation" dynamics are a very specific simplification of true dynamics. Let us agree that true dynamics again include many timescales and confounding factors and that therefore this manuscript is no less finetuned than the theory it challenges. What the authors essentially do is, impose specific dynamics and then demonstrate that these specific dynamics fail to maintain diversity at . It seems likely that their "mutation" dynamics can be plausibly modified so that diversity is regained. After all, natural ecosystems do involve resource competition and the plankton are diverse.
Reading the other referee reports, it appears that some referees do not consider this paper sufficiently surprising to be published in eLife. However, one referee suggested that it be published as a Research Advance, to the original Wingreen paper published in eLife. I think that if the authors sufficiently improve this manuscript, the Research Advance track might make sense. Specifically, what I suggest, is for the authors to readjust their focus. Their results are well argued but simply show a known weakness in the theory – a weakness that can overcome. Why not argue them in a way that seeks to expand the field? The authors raise an interesting question, what needs to be added to the plainvanilla version of consumerresource models so that diversity is regained despite (slightly) nonlinear tradeoffs and a specific form of "mutation" dynamics? More elaborate model variants already incorporate other diversity stabilizing mechanisms which might well maintain diversity at nonunity values of γ, e.g. a recent preprint by Huang group in Stanford https://www.biorxiv.org/content/10.1101/2021.05.13.444061v1. Should the authors answer this question, and salvage diversity from nonlinear tradeoffs, I think this manuscript will be much improved. Exploring an example whereby diversity is reinstated in the consumerresource framework (despite nonlinear tradeoffs) would demonstrate how in fact, suitably adjusted, consumerresource models can be used to capture such competing ecological forces. I believe that by steering their manuscript to open new avenues of thought rather than closing avenues of thought, it would promote future inquiries in the field, hopefully ultimately leading to a deeper understanding. With this suitable addition and adjustment of the manuscript, I hope that the Editor and other referees would then agree that it would pass the threshold of contribution to be included as a Research Advance.
2. If indeed this paper is to be considered to be published in eLife as a Research Advance following the Wingreen eLife paper, it is a good idea that the authors change their notation to match precisely that parent paper's notations. I understand this is an unusual request but I do believe that it would serve future readers best – to smoothly carry over notations from one paper to its immediate followup in the same journal. As it is, the notation differences between the two papers are small and this is not a big request. Future readers will thank you.
https://doi.org/10.7554/eLife.67764.sa1Author response
[Editors’ note: The authors appealed the original decision. What follows is the authors’ response to the first round of review.]
Reviewer #1:
The authors started by introducing the adaptive dynamics framework which has been used to study evolutionary diversification due to frequencydependent selection. The introduction is not adequate in giving nonspecialists an intuitive feeling about how evolutionary branching point from fitness minimum works. Of course, making theoretical conclusions accessible is not trivial, and may not be achievable. However, authors can try harder by including supplementary figures.
We rearranged the text, moving the general description of the adaptive dynamics to
the Appendix and focusing more on the predictions derived from both the numerical simulation of the complete model and the adaptive dynamics treatment of its approximate version. We also tried to better explain the convergence and evolutionary stability of the fixed point. This made the corresponding part of the Appendix a fairly complete verbal description of the standard adaptive dynamics procedure, illustrated by the equations pertinent to our model. However, given that the adaptive dynamics framework is almost 30 years old, and that many of the original and review papers on adaptive dynamics are already cited in the manuscript (there are thousands of publications based on adaptive dynamics), we stayed short of retelling a more detailed description of the general adaptive dynamics procedure.
The authors then introduced work from the Wingreen group: R number of resources in a wellmixed environment, and each consumer species has an uptake strategy for each of the R resources. The total uptake for each species is fixed, creating a tradeoff: more uptake of one resource reduces the uptake of another. Previous work showed that when the tradeoff function is linear (e.g. the sum of [uptake of each resource)^{γ}] = a constant where the exponent γ=1, then many species with different uptake strategies can coexist.
The authors showed that the conclusion from previous work is rather restricted in its scope: large diversity only exists when the exponent γ in the tradeoff function is 1 (i.e. linear). When that exponent is greater than 1, no diversification occurs (saving for mutationselection balance) and all individuals are phenotypically similar generalists. When the exponent is less than 1 (i.e. concave), the initial convergence to an unstable steady state later evolutionarily diverges into specialists, each specializing on a resource.
Authors also tested their conclusions within the ecological framework of the Wingreen group, and in different scenarios (chemostat versus serial dilutions), and reached the same conclusions.
The paper is relatively easy to read (in the realm of theoretical papers). This work reminds me of the work that has been done in the field of the evolution of cooperation. For example, public goods games often assume that the effect of the public good is a linear function of the number of contributions, an assumption that is often violated in biology. Depending on whether this function is linear or nonlinear, one can get very different outcomes in cooperator/noncooperator coexistence (e.g. Archetti and Scheuring, JTB 299:920; Damore and Gore, JTB, https://doi.org/10.1016/j.jtbi.2011.03.008). Authors may want to add a discussion on that.
We thank the Reviewer for bringing up the connection between these two fields. Not
only the curvature of the benefit and cost functions play essential roles in the evolution of cooperation, the biochemical justifications of nonlinearity in these functions are equally applicable to the case of nonlinear uptake rates and tradeoffs used in our work. We added the suggested references and a short explanation to the Discussion.
I wonder whether this paper should be added as a Research Advance to the original paper from the Wingreen paper published in eLife.
We agree that this would be the best place to publish our work.
Reviewer #2:
My main concern about this work is whether it has a sufficient degree of novelty and interest. As mentioned in the public review, the results are robust. But, after a close analysis of the model by Posfai et al., they are not unexpected. The manuscript, as it stands, mostly demonstrates the weaknesses of the paper by Posfai et al.: fixed points with more species than resources exist, but they are only marginally stable
(https://arxiv.org/pdf/2002.04358). Structural instability is a direct consequence of this fact. Nonlinear tradeoffs are just one (among many) ways to show that the results are infinitely sensitive to the parameter choices.
As we stated at the beginning of this letter, we did our best to find any publication that
considers the models suggested in3,5, reports their structural instability, and presents the evolutionary and ecological results in the two generic scenarios of concave and convex tradeoffs. We have not found any, and therefore we believe that such a publication is in order. We agree that there could be many other mechanisms that break the neutral evolutionary scenario and reduce the number of coexisting species to the competitive exclusion limit. Our goal in this work is to consider the interesting metabolic models developed in (3, 5) preserving their mechanistic definitions, yet supplying them with more realistic functional forms of tradeoffs and nutrient uptake rates.
Reviewer #3:
The work would be of greater general interest if the biological evidence for nonlinearities in enzyme costs were carefully examined; mechanistic insights on how enzyme budget nonlinearities may arise in nature would be of significant utility to the field. However, this would require a substantial additional undertaking and would shift the focus of the work from the specific implications of nonlinearities in resourcecompetition models. An alternative would be to publish the current study in a more specialized journal, with a more theoretical focus.
In short, the fundamental reason for nonlinearities in tradeoffs (enzyme costs) and uptake rates, which we showed to be interchangeable, is the nonlinear dependence of chemical kinetics on the concentrations of reactants. A more detailed explanation is presented at the beginning of this letter. We added paragraphs to the Introduction and Discussion that present and discuss the evidence for nonlinearity.
It would also be helpful to provide a quantitative assessment of the sensitivity of diversity to the degree of nonlinearity. It is clear that any nonlinearity (or inexactness of tradeoffs) leads to loss of diversity at long times. However, a small rate of invasion by new strategies can still lead to a diverse stationary state of the population. Given a certain degree of nonlinearity, how much invasion is required to maintain diversity? The adaptive dynamics calculations performed by the authors do not address this point because new strategies are only introduced if these are more fit than the residents. The question of diversity requires introducing invaders that may be slightly less fit, but still manage to survive due to demographic noise.
$\frac{\beta \delta}{\beta \phantom{\rule{0.222em}{0ex}}}$
To give credit to adaptive dynamics, it does provide some estimates to how advantageous or disadvantageous ecological conditions are for an invader or a mutant. A selection gradient (e.g. (2122)) is an indicator of the steepness of the fitness landscape, i.e. how quickly growth rate advantages or disadvantages change with the phenotypic separation from the resident population. Likewise, the Hessian of the invasion fitness (28) serves as an estimate of the secondorder effects of phenotypic separation on the relative fitness. For example, the presence of 1− ƴ factor in (28) indicates that the strength of deviation from neutral selection increases in the leading order linearly with the deviation of the tradeoff exponent ƴ from one.
Some information about the fate of an invader can be obtained from the videos corresponding to Figure 3 in the Appendix. Those simulations are initiated with a dense species packing and the majority of species quickly go extinct. To show that the effect is pronounced even for slight deviation from nonlinearity, we chose the values for the tradeoff exponent close to one, ƴ = 0.9 and ƴ = 1.1.
Studying the extinction rates of invaders under the conditions of the models (3, 5) could be an interesting separate project. In the present manuscript we prefer to focus on the steady state of evolving diversity.
References
[1] J. A. Borghans, R. J. De Boer, and L. A. Segel. Extending the quasisteady state approximation by changing variables. Bulletin of mathematical biology, 58(1):43–63, 1996.
[2] A. Ciliberto, F. Capuani, and J. J. Tyson. Modeling networks of coupled enzymatic
reactions using the total quasisteady state approximation. PLoSComput Biol, 3(3):e45,
2007.
[3] A. Erez, J. G. Lopez, B. G. Weiner, Y. Meir, and N. S. Wingreen. Nutrient levels and
tradeoffs control diversity in a serial dilution ecosystem. eLife, 9:e57790, sep 2020.
ISSN 2050084X.
[4] C. M. Hill, R. D.Waightm, andW. G. Bardsley. Does any enzyme follow the michaelismenten equation? Molecular and cellular biochemistry, 15(3):173–178, 1977.
[5] A. Posfai, T. Taillefumier, and N. S. Wingreen. Metabolic tradeoffs promote diversity
in a model ecosystem. Phys. Rev. Lett., 118:028103, Jan 2017.
[Editors’ note: what follows is the authors’ response to the second round of review.]
Essential Revisions:
Please consult the Reviewers' comments and address these in your revision.
Reviewer #1:
I remain supportive, especially if authors can discuss empirical measurements of biological tradeoffs and whether in the natural environment the linearity assumption might break down.
Following the recommendation of Reviewer 1, in the new version of the manuscript we expanded the Introduction, Results and Discussion sections (see text in blue in the accompanying document). Assessing the empirical evidence for nonlinearity in tradeoffs and uptake rates, we focused on similarity in concentrations of metabolic enzymes and substrates, and on measured reaction elasticities and flux control coefficients as the most direct quantitative indicators of the “degree” of nonlinearity. We also emphasized that our work leads to empirically testable predictions, at least in principle: microbial ecosystems with tradeoffstructures that are close to linear should tend to be more diverse than those with highly nonlinear tradeoff structures (see new paragraph in the Results section). Finally, we propose possible ways to find a solution to the paradox of plankton based on nonstationary endogenous population dynamics (see new paragraph in the Discussion section).
Reviewer #3:
The authors have made several cosmetic changes which have improved the clarity of their manuscript. However, the revised version does not substantially address my main concerns. The real question is whether the current manuscript makes a substantial contribution to the topic of microbial diversity. The focus of the paper is a critique of a model of resource competition with tradeoffs. It is certainly legitimate to be critical of existing models. However, I believe the readers of eLife already appreciate the adage "all models are wrong, but some are useful". The authors have focused their attention on the first part of the adage, arguing that because growth functions will not be exactly linear the model is "wrong". But it's not news that the model is "wrong" (see above), the question is whether the model might still be a useful starting point for understanding diversity? What seems to me to be missing in the discussion, both in the original studies of Ref. 3 and 5 and the current manuscript, is quantification of how "wrong" the initial model is, and whether this undermines its utility. This is why I suggested that the authors carefully examine the "biological evidence for nonlinearities in enzyme costs". Their revised manuscript adds some sentences on this point, but in a nonquantitative way: the authors continue to make the mathematical point that the model is "wrong", but have not taken up the challenge of addressing whether it is or is not "useful". Yes, it is mathematically correct as the authors state that bimolecular reactions are strictly nonlinear in the reactants. But for typical enzyme concentrations in the μM range and typical metabolites in the 0.11 millimolar range, these nonlinearities are in the 0.11% range. From a biological perspective, a linear function might therefore still be a useful starting point. I've also read the references the authors cite on other sources of nonlinearity – they are equally nonquantitative. For example, the review by Marianayagam et al. states (without citations) "In its simplest form, oligomerization functions as a general mechanism for sensing protein concentration. An increase in protein concentration above the oligomerization threshold can be the stimulus for enzyme activation; similarly, enzyme deactivation will apply when cellular levels of the enzyme fall." I absolutely agree that for enzymes that need to oligomerize to function, this implies a mathematically nonlinear processing rate as a function of enzyme concentration. However, again for enzyme levels in the μM range and oligomerization dissociation constants in the commonly observed 110 picomolar range, the nonlinearities are again ~0.11%. Despite the revisions on this and other points in the reviews, in the end I am left still wondering whether the original model is "useful". My conclusion is that the current manuscript will be primarily of interest to researchers whose focus is on the mathematics of resourcecompetition models, and would therefore be appropriate for a more mathematically focused journal.
We also appreciate George Box’s ”All models are wrong, but some are useful” quote and frequently use it in teaching to motivate students to aspire to design useful models.
However, there are different kinds of ”wrong” in modeling. There is Box’s wrong, which
applies to all models, as they necessarily always represent a simplification of reality. But
there is also a more technical, and arguably more important model failure, which occurs
when a model is structurally unstable, which essentially means that if a particular model parameter is changed by a very small amount, the main quantitative and qualitative behaviour of the model changes completely. This is the case for Wingreen’s model with regard to the model parameter, which was set to 1 in Wingreen’s paper. As we show in our contribution, the model outcomes completely change when ƴ≠ 1 (no matter how large or small the deviation from 1 is, see below!). Thus, Wingreen’s model is wrong in the sense that it is structurally unstable. To suggest that our model is equally ”wrong” misses the point of structural instability.
Regarding usefulness, we maintain that generally assuming linearity of tradeoff appears
to be wrong biologically. That is, we argue that Wingreen’s model is too restrictive to be useful. Assuming that George Box’s statement was primarily related to his main field, statistics, we further support this claim by statistical arguments, as follows.
While we agree that there are metabolites, particular in fastgrowing bacteria, with concentrations well in the millimolar range, many common metabolite have concentrations in the range of tens and even single digits of micromoles. For example, a detailed study (2) of the canonical model microbe E. coli during exponential growth in glucose, glycerol, or acetate as the carbon source revealed that out of 103 metabolites, 35 have concentrations above 1mM, but the concentrations of 46 metabolites are in tens or single micromole digits, including 2 metabolites with concentrations below 1 μM (as summarized in Table 1 in that paper). The database BIONUMBERS7 estimates the typical metabolite concentration in an E. coli bacterium as 32 μM, and provides evidence for concentrations of important E. coli glycolysis enzymes in tens and even hundreds of μM. Hence, while it looks natural that concentrations of metabolic substrates are larger than those of the corresponding enzymes, the difference is commonly few to tenfold rather than 23 orders of magnitude, as suggested by Reviewer 3.
A similar mischaracterization appears in the Reviewer’s estimate of the degree of nonlinearity caused by the disparity between the “oligomerization dissociation constants in the commonly observed 110 picomolar range” (”picomolar” is apparently rarely heard of, perhaps a typo in the Reviewer comments?), and concentrations of their monomer constituents. In reality, the difference between these quantities does not seem to be nearly as pronounced as the reviewer claims, as can e.g. be seen from the following summarizing quote in the review (1), page 5014: “The association between subunits can vary in strength and duration. Some proteins are found only, or primarily, in the oligomeric state. These proteins generally have dissociation constants in the nanomolar range. Others have a weak tendency to associate, with oligomerization dependent on environmental conditions, such as concentration, temperature, and pH. Such proteins often have higher Kd values in the μM or even millimolar range. Still other proteins oligomerize dynamically in response to a stimulus, such as a change in nucleotide binding, nucleotide hydrolysis or phosphorylation state. Such a change can have a dramatic effect on the affinity of the subunits for one another, often by orders of magnitude.”
Another review (15) discusses enzymes that are known to regulate their in vivo catalytic activity by dimeric and oligomeric associationdissociation, which rules out saturated binding between monomers and a corresponding linear dependence between the concentrations of monomers and complexes.
In addition to the comparison between enzyme and substrate concentrations and enzyme dissociation constants (tricky to measure in vivo even today, e.g.5), there is a more direct way to show omnipresence of nonlinearity in metabolic reactions and fluxes. The Reviewer’s assessment that ”these nonlinearities are in the 0.11% range” probably refers to the definition of the degree of nonlinearity of a function f(x) as a doublelogarithmic derivative, d ln[(f(x))]/d ln(x). For a general powerlaw function, such as f(x) ≡ Cx that we used to define tradeoffs, the loglog derivative is equal to, the nonlinearity parameter. For example, the loglog derivatives with respect to enzyme concentration of the tradeoff used in8 are equal to one. This derivative is equal to the local value of an exponent from an approximation (also local) of the function f(x) by a power law, and has long been used in metabolic analysis, being called an “elasticity coefficient” when f(x) is a particular reaction rate and x is the concentration of enzyme, or a flux control coefficient when f(x) is a flux through a metabolic pathway and x is again the concentration of enzyme. Given that we are more interested in the effect of concentration of enzymes on the growth rate, the latter definition seems more relevant.
As with metabolic control analysis and metabolic engineering, the mere existence of
these terms indicates that their values are often distinct from one. Indeed, such is the case in publications that present experimentally measured values of these coefficients, or quantities that are derived from those coefficients. For example, flux control coefficients of metabolic enzymes and plots of fluxes vs. enzyme concentrations shown in Figure 6 in (6), Table 6 in (4), Table 2 in (14), Figure 2 in (11), Figure 1 in (9), Table 5 in (10), Figure 2B in (13), Table 3 in (16) reveal that these coefficients rarely come close to 1, and that the dependencies of metabolic fluxes on enzyme concentrations are significantly nonlinear. (We found these and many other examples using Google Scholar to search for Flux Control Metabolic Coefficients Bacteria.) Indeed, back in the 60s Michael Savageau suggested to parametrize rates of complex enzymatic reactions as products of powerlaw functions of concentrations of enzymes and substrates12. This idea developed a substantial following, which once again indicates the necessity to account for nonlinearity in kinetics of enzymatic pathways. Overall, the 0.11% range for nonlinearities seems doubtful.
In conclusion, it seems highly questionable to simply assume that all biologically relevant tradeoffs are close to linear. Rather, it seems very plausible that many of them are highly nonlinear. To argue that readers of eLife would only need to know what happens in generic consumerresource models when tradeoffs are linear therefore does not make sense to us. On the contrary, we think that readers of Wingreen’s papers, including their eLife paper, should be made aware of the structural instability of their results. We emphasize that because of the structural instability, tradeoffs would have to be exactly linear to arrive at Wingreen’s results, as any slight deviation from linearity destroys the neutrality necessary for excess diversity in Wingreen’s model, and hence results in completely different outcomes. This can be seen in (Author response image 1A) which shows that even for ƴ = 0.99, i.e., for ƴ within 1% of the linear case, there is no excess diversity. The case for ƴ = 1.01 is shown in (Author response image 1B) and the full evolutionary dynamics for this case can be found at https://figshare.com/s/f65ed0bf9b4305e9018f.
In the updated text, we substantially edited the Introduction and Discussion to point to
evidence in favour of the importance of nonlinear kinetics and tradeoffs, and to reflect the above considerations. We also added a paragraph to the Results section discussing systems with different degrees of nonlinearity and the expectations regarding the amount of diversity that can be derived based on the degree of nonlinearity in tradeoff structures.
Reviewer #4:
The reality of the field (if I may) is that none of us have the answer to the socalled Paradox of the Plankton. There are several competing theories out there, and each of them has its merits and imperfections. Simply outlining an expected and obvious lack of generalizability of a particular aspect of an otherwise solid and informative theory is less constructive than this paper could be. Especially considering it is well known that consumerresource models (e.g. Wingreen) support coexistence only on equal enzyme budgets unless other diversity stabilizing mechanisms are included.
Again, this seems strange to us. It is exactly papers like ours that go some way to clarify the situation that a general model with nonlinear tradeoff does not “support coexistence only on equal enzyme budgets ” (assuming that it refers to excessive diversity), but the Reviewer seems to think that all this is already well known. As we show, Wingreen’s models are structurally unstable, and therefore not ”solid and informative”, as the Reviewer claims. To us, it seems important that this is pointed out in the literature.
What I think that the authors should focus on, is that neither the Wingreen papers nor recent papers from other groups claim that their particular framework captures all features of diversity maintenance. Instead, we are all trying, as a community, to piece together a coherent theory, stumbling along as we inevitably do. For example, the Wingreen papers consider violation of the linear constraint (i.e. unequal enzyme budgets) in several stochastic and deterministic cases, outlining their relevance.
As far as we can tell, both papers coming from Wingreen’s group regarding this topic assume linear tradeoffs. We are unsure why the reviewer claims otherwise. We agree that we are all trying, and that all models and approaches have their own benefits and shortcomings. However, if a model has a particularly severe shortcoming, such that the
results simply don’t withstand scrutiny, this needs to be pointed out. In our opinion, this is the case for the structural instability due to assuming linear tradeoffs in Wingreen’s models.
There are several timescales at play here and one can think how their interplay with particular model parameters influences the resulting steady state. This is likely the case with regards to this manuscript. Already, if one considers small deviations from γ=1, one introduces a new timescale, to compete with other model timescales. Moreover, the "mutation" dynamics (Eq. 14) introduce one or more timescales, which are conveniently swept under the rug with a "without loss of generality" (l463), a statement without actual backing as far as I understood. Furthermore, the authors do not sufficiently acknowledge that their suggested "mutation" dynamics are a very specific simplification of true dynamics. Let us agree that true dynamics again include many timescales and confounding factors and that therefore this manuscript is no less finetuned than the theory it challenges. What the authors essentially do is, impose specific dynamics and then demonstrate that these specific dynamics fail to maintain diversity at . It seems likely that their "mutation" dynamics can be plausibly modified so that diversity is regained. After all, natural ecosystems do involve resource competition and the plankton are diverse.
Introducing ƴ≠1 does not introduce new timescales. It is true that our evolutionary framework does introduce a new (evolutionary) timescale that is dictated by mutations. However, this merely serves to generalize the model, and it does not affect the basic finding of the effect of nonlinear tradeoffs on diversity. This can be seen by considering purely ecological dynamics, as happens in our models either at the evolutionary end state (when there is no evolutionary change anymore), or by simply setting mutations to 0, as we have done in the supplementary material. The resulting ecological dynamics show exactly the same effect of nonlinear tradeoffs on diversity as the more general evolutionary model.
Thus, it is not the evolutionary time scale that causes these effects. As an aside: what we
meant by the phrase “Without loss of generality, we set σ = 1” quoted by the Reviewer, is that the constant σ, which essentially describes the rate and size of mutations, can be absorbed by defining a new time t’ = σt, which simply affects the frame rate in our videos, rather than what is in those frames.
Reading the other referee reports, it appears that some referees do not consider this paper sufficiently surprising to be published in eLife. However, one referee suggested that it be published as a Research Advance, to the original Wingreen paper published in eLife. I think that if the authors sufficiently improve this manuscript, the Research Advance track might make sense. Specifically, what I suggest, is for the authors to readjust their focus. Their results are well argued but simply show a known weakness in the theory – a weakness that can overcome. Why not argue them in a way that seeks to expand the field? The authors raise an interesting question, what needs to be added to the plainvanilla version of consumerresource models so that diversity is regained despite (slightly) nonlinear tradeoffs and a specific form of "mutation" dynamics? More elaborate model variants already incorporate other diversity stabilizing mechanisms which might well maintain diversity at nonunity values of γ, e.g. a recent preprint by Huang group in Stanford https://www.biorxiv.org/content/10.1101/2021.05.13.444061v1. Should the authors answer this question, and salvage diversity from nonlinear tradeoffs, I think this manuscript will be much improved. Exploring an example whereby diversity is reinstated in the consumerresource framework (despite nonlinear tradeoffs) would demonstrate how in fact, suitably adjusted, consumerresource models can be used to capture such competing ecological forces. I believe that by steering their manuscript to open new avenues of thought rather than closing avenues of thought, it would promote future inquiries in the field, hopefully ultimately leading to a deeper understanding. With this suitable addition and adjustment of the manuscript, I hope that the Editor and other referees would then agree that it would pass the threshold of contribution to be included as a Research Advance.
Once again, we are a bit perplexed by these comments. We tried to be succinct in our presentation of the fact that Wingreen’s model is structurally unstable, and that with
nonlinear tradeoffs, no excess diversity can be maintained. Despite succinctness, this resulted in a fullfledged paper. Now the Reviewer wants us to write an additional paper about possible mechanisms that can lead to excess diversity with nonlinear tradeoffs. This seems strange. Shouldn’t it be one paper at a time? Moreover, despite what the Reviewer may think, it is not clear at all at this point what possible mechanisms could indeed lead to excess diversity with nonlinear tradeoffs. One possibility would be to consider nonequilibrium ecological dynamics. This has been done in the eLife paper from the Wingreen group to which our paper refers. In that paper, the authors considered seasonal ”batch culture” dynamics rather than equilibrium consumerresource dynamics. However, they still always assumed linear tradeoffs! Note that this is a paper all of its own, and not merely an additional section added on to their original PRL paper…how could it be otherwise? We note that such batch culture dynamics do not lead to excess diversity with nonlinear tradeoffs, as we have shown in our manuscript.
We are equally keen to find a solution to the Paradox of Plankton, and it is possible that
nonequilibrium ecological dynamics can allow for the maintenance of excess diversity, as we have recently shown using a different ecological model3. More precisely, we show that endogenous nonstationary “boombust” population dynamics can lead to a fewfold increase in diversity above the saturation limit expected with equilibrium population dynamics. This together with many experimental results reporting nonstationarity and apparent chaoticity of the population dynamics of actual plankton species makes us believe that the key to explain the astounding diversity of species is not the neutral evolutionary regime predicted in 8, but rather nonstationary population dynamics induced by competition and predation and perhaps external factors. We think that the neutral evolutionary state and linear tradeoff to 9 which an evolving system may accidentally “selforganize” is very fragile and at most temporal and cannot play a major role in establishing and maintaining (even less likely!) excess diversity.
However, it seems obvious that to extend investigations of nonequilibrium ecological
dynamics with nonlinear tradeoffs to Wingreen’s consumerresource model would need a new research project that would require its own space, and cannot be added to a paper about the structural instability of Wingreen’s model. The preprint cited by the Reviewer also proves the point: such investigations must be done in their own right, not tacked on to something else. Clearly, we must leave the exploration of excess diversity with nonlinear tradeoffs to the future.
We have added a paragraph to the Discussion where we bring up the possibility of exploring the effect of nonlinear tradeoffs in the presence of nonequilibrium ecological dynamics.
2. If indeed this paper is to be considered to be published in eLife as a Research Advance following the Wingreen eLife paper, it is a good idea that the authors change their notation to match precisely that parent paper's notations. I understand this is an unusual request but I do believe that it would serve future readers best – to smoothly carry over notations from one paper to its immediate followup in the same journal. As it is, the notation differences between the two papers are small and this is not a big request. Future readers will thank you.
There seem to be four discrepancies between our and Wingreen’s notation:
the number of resources p in 8 vs R (us), resource input rates in 8 vs S (us), death rate _ in8 vs d (us), and a dedicated summation index _ in sums over all species in 8, which is absent in our paper, because it is not needed. We have fixed the first three discrepancies, making our notation identical to that in 8. As for the summation index, we think that our notation is more succinct, less cumbersome, and actually not contradicting that in 8, so we kept it.
References
[1] M. H. Ali and B. Imperiali. Protein oligomerization: how and why. Bioorganic and
medicinal chemistry, 13(17):5013–5020, 2005.
[2] B. D. Bennett, E. H. Kimball, M. Gao, R. Osterhout, S. J. Van Dien, and J. D. Rabinowitz.
Absolute metabolite concentrations and implied enzyme active site occupancy
in Escherichia coli. Nature chemical biology, 5(8):593–599, 2009.
[3] M. Doebeli, E. C. Jaque, and Y. Ispolatov. Boombust population dynamics increase
diversity in evolving competitive communities. Communications biology, 4(1):1–8,
2021.
[4] C. Giersch. Determining elasticities from multiple measurements of flux rates and
metabolite concentrations: application of the multiple modulation method to a reconstituted pathway. European journal of biochemistry, 227(12):194–201, 1995.
[5] I. Jarmoskaite, I. AlSadhan, P. P. Vaidyanathan, and D. Herschlag. How to measure and evaluate binding affinities. eLife, 9:e57264, 2020.
[6] A. J. Loder, Y. Han, A. B. Hawkins, H. Lian, G. L. Lipscomb, G. J. Schut, M. W. Keller, M. W. Adams, and R. M. Kelly. Reaction kinetic analysis of the 3hydroxypropionate/4ydroxybutyrate co2 fixation cycle in extremely thermoacidophilic archaea. Metabolic engineering, 38:446–463, 2016.
[7] R. Milo, P. Jorgensen, U. Moran, G. Weber, and M. Springer. Bionumbersthe database
of key numbers in molecular and cell biology. Nucleic acids research, 38(suppl 1):
D750–D753, 2010.
[8] A. Posfai, T. Taillefumier, and N. S. Wingreen. Metabolic tradeoffs promote diversity
in a model ecosystem. Phys. Rev. Lett., 118:028103, Jan 2017.
[9] J. M. Rohwer, N. D. Meadow, S. Roseman, H. V. Westerhoff, and P. W. Postma. Understanding glucose transport by the bacterial phosphoenolpyruvate: glycose phosphotransferase system on the basis of kinetic measurements in vitro. Journal of Biological Chemistry, 275(45):34909–34921, 2000.
[10] R. Rutkis, U. Kalnenieks, E. Stalidzans, and D. A. Fell. Kinetic modelling of the
zymomonasmobilisentner–doudoroff pathway: insights into control and functionality.
Microbiology, 159(Pt 12):2674–2689, 2013.
[11] E. Saavedra, R. Encalada, E. Pineda, R. JassoCh´avez, and R. MorenoS´anchez. Glycolysis in entamoeba histolytica: biochemical characterization of recombinant glycolytic enzymes and flux control analysis. The FEBS journal, 272(7):1767–1783, 2005.
[12] M. A. Savageau. Biochemical systems analysis: I. some mathematical properties of the rate law for the component enzymatic reactions. Journal of theoretical biology, 25(3): 365–369, 1969.
[13] A. Schmidt, K. Kochanowski, S. Vedelaar, E. Ahrn´e, B. Volkmer, L. Callipo,
K. Knoops, M. Bauer, R. Aebersold, and M. Heinemann. The quantitative and
conditiondependent Escherichia coli proteome. Nature biotechnology, 34(1):104–110,
2016.
[14] Y. Sun and S. Qian. Flux control analysis for biphenyl metabolism by rhodococcus
pyridinovorans r04. Biotechnology letters, 24(18):1525–1529, 2002.
[15] T. W. Traut. Dissociation of enzyme oligomers: a mechanism for allosteric regulation.
Critical reviews in biochemistry and molecular biology, 29(2):125–163, 1994.
[16] J. Van Der Vlag, R. Van’t Hof, K. Van Dam, and P. W. Postma. Control of glucose
metabolism by the enzymes of the glucose phosphotransferase system in Salmonella
typhimurium. European journal of biochemistry, 230(1):170–182, 1995.
https://doi.org/10.7554/eLife.67764.sa2Article and author information
Author details
Funding
FONDECYT (1200708)
 Yaroslav Ispolatov
NSERC (Discovery Grant 219930)
 Michael Doebeli
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Senior Editor
 Aleksandra M Walczak, École Normale Supérieure, France
Reviewing Editor
 Wenying Shou, University College London, United Kingdom
Reviewer
 Wenying Shou, University College London, United Kingdom
Publication history
 Preprint posted: October 21, 2020 (view preprint)
 Received: February 22, 2021
 Accepted: August 5, 2021
 Accepted Manuscript published: August 5, 2021 (version 1)
 Version of Record published: September 9, 2021 (version 2)
Copyright
© 2021, Caetano et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,199
 Page views

 152
 Downloads

 3
 Citations
Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Evolutionary Biology
 Genetics and Genomics
The genetic variants introduced into the ancestors of modern humans from interbreeding with Neanderthals have been suggested to contribute an unexpected extent to complex human traits. However, testing this hypothesis has been challenging due to the idiosyncratic population genetic properties of introgressed variants. We developed rigorous methods to assess the contribution of introgressed Neanderthal variants to heritable trait variation relative to that of modern human variants. We applied these methods to analyze 235,592 introgressed Neanderthal variants and 96 distinct phenotypes measured in about 300,000 unrelated white British individuals in the UK Biobank. Introgressed Neanderthal variants have a significant contribution to trait variation consistent with the polygenic architecture of complex phenotypes (contributing 0.12% of heritable variation averaged across phenotypes). However, the contribution of introgressed variants tends to be significantly depleted relative to modern human variants matched for allele frequency and linkage disequilibrium (about 59% depletion on average), consistent with purifying selection on introgressed variants. Different from previous studies (McArthur 2021), we find no evidence for elevated heritability across the phenotypes examined. We identified 348 independent significant associations of introgressed Neanderthal variants with 64 phenotypes . Previous work (Skov 2020) has suggested that a majority of such associations are likely driven by statistical association with nearby modern human variants that are the true causal variants. We therefore developed a customized statistical finemapping methodology for introgressed variants that led us to identify 112 regions (at a false discovery proportion of 16%) across 47 phenotypes containing 4,303 unique genetic variants where introgressed variants are highly likely to have a phenotypic effect. Examination of these variants reveal their substantial impact on genes that are important for the immune system, development, and metabolism. Our results provide the first rigorous basis for understanding how Neanderthal introgression modulates complex trait variation in presentday humans.

 Developmental Biology
 Evolutionary Biology
The gill skeleton of cartilaginous fishes (sharks, skates, rays, and holocephalans) exhibits a striking anterior–posterior polarity, with a series of fine appendages called branchial rays projecting from the posterior margin of the gill arch cartilages. We previously demonstrated in the skate (Leucoraja erinacea) that branchial rays derive from a posterior domain of pharyngeal arch mesenchyme that is responsive to Sonic hedgehog (Shh) signaling from a distal gill arch epithelial ridge (GAER) signaling centre. However, how branchial ray progenitors are specified exclusively within posterior gill arch mesenchyme is not known. Here, we show that genes encoding several Wnt ligands are expressed in the ectoderm immediately adjacent to the skate GAER, and that these Wnt signals are transduced largely in the anterior arch environment. Using pharmacological manipulation, we show that inhibition of Wnt signalling results in an anterior expansion of Shh signal transduction in developing skate gill arches, and in the formation of ectopic anterior branchial ray cartilages. Our findings demonstrate that ectodermal Wnt signalling contributes to gill arch skeletal polarity in skate by restricting Shh signal transduction and chondrogenesis to the posterior arch environment and highlights the importance of signalling interactions at embryonic tissue boundaries for cell fate determination in vertebrate pharyngeal arches.