Evolutionary coexistence in a fluctuating environment by specialization on resource level

Microbial communities in fluctuating environments, such as oceans or the human gut, contain a wealth of diversity. This diversity contributes to the stability of communities and the functions they have in their hosts and ecosystems. To improve stability and increase production of beneficial compounds, we need to understand the underlying mechanisms causing this diversity. When nutrient levels fluctuate over time, one possibly relevant mechanism is coexistence between specialists on low and specialists on high nutrient levels. The relevance of this process is supported by the observations of coexistence in the laboratory, and by simple models, which show that negative frequency dependence of two such specialists can stabilize coexistence. However, as microbial populations are often large and fast growing, they evolve rapidly. Our aim is to determine what happens when species can evolve; whether evolutionary branching can create diversity or whether evolution will destabilize coexistence. We derive an analytical expression of the invasion fitness in fluctuating environments and use adaptive dynamics techniques to find that evolutionarily stable coexistence requires a special type of trade‐off between growth at low and high nutrients. We do not find support for the necessary evolutionary trade‐off in data available for the bacterium Escherichia coli and the yeast Saccharomyces cerevisiae on glucose. However, this type of data is scarce and might exist for other species or in different conditions. Moreover, we do find evidence for evolutionarily stable coexistence of the two species together. Since we find this coexistence in the scarce data that are available, we predict that specialization on resource level is a relevant mechanism for species diversity in microbial communities in fluctuating environments in natural settings.

| 623 WORTEL stabilization of diversity from frequency-dependent interactions (Chesson, 2000). An example of an environment that received a lot of attention lately and where microbial species and strains coexist over long timescales is the human gut microbiota (Garud et al., 2019;Zhao et al., 2019). A main property of this environment is that there are fluctuating nutrient conditions. Even in seemingly simple environments, with as only complexity fluctuating resources, long-term coexistence of different microbial strains is very widespread (Good et al., 2017). Several explanations have been given for this coexistence, including crossfeeding (Doebeli, 2002;D'Souza et al., 2018;Goyal et al., 2021;Smith et al., 2019), specialization on different resources (Estrela et al., 2021;Wickman et al., 2019) and specialization on resource level (Stewart & Levin, 1973). The first observation of coexistence of two microbial strains in an environment with a fluctuating resource level is by Levin (1972). With a mathematical toy model, Stewart & Levin (1973) showed that, when one species takes up nutrients faster at low resource levels and the other at high, coexistence is possible. An extension of the model with fixed cycle lengths and lag phases also shows possibilities of coexistence (Smith, 2011). In a different experiment, Turner et al. (1996) observed coexistence and found evidence of a trade-off on growth at different resource levels, but a mathematical model based on Monod growth, described in the same article, showed that this was not enough to sustain coexistence. Closer analysis showed that indeed, in addition to the resource-level specialization trade-off, there was cross-feeding between the strains.
The above studies are all restricted to analysing strains with fixed properties. To observe coexistence of cohabiting cells, the coexistence should have some evolutionary stability (Egas et al., 2004), that is, adaptation should not lead to the loss of the coexistence.
We therefore aim to study the relevance of specialization on resource level for evolutionary stable coexistence in fluctuating environments. Thereby we link ecological studies of coexistence of specialists on low-and high-resource concentrations to studies of evolutionary dynamics in microbial evolution. The adaptive dynamics framework is well suited to study these evolutionary dynamics, because this framework links ecological dynamics, such as the fluctuations and changes in population sizes, with evolving traits (Dieckmann & Law, 1996;Geritz et al., 1998;Metz et al., 1992).
Adaptive dynamics has been used to study microbial population dynamics of long-term evolution studies (Friesen et al., 2004).
In this paper, we present an analytical expression for the invasion fitness in serial dilution experiments, which is typically defined as the relative fitness over a full fluctuation cycle (called the selection coefficient Chevin, 2011;Lenski et al., 1991;Manhart & Shakhnovich, 2018). Since this concerns the common way to do laboratory evolution experiments with micro-organisms, this result can directly be applied to experimental design and interpretation. We also supply a method to use experimental data together with computational models to obtain an explicit trade-off function. We show that certain trade-off functions can and some cannot lead to evolutionary stable coexistence. We combine the obtained trade-off data with the trade-off necessary for evolutionary coexistence and postulate that evolutionary coexistence is most likely in communities of different species rather than different phenotypes of the same species.

| MODEL DE SCRIP TION
The model describes an environment where microbes grow and resources are consumed. At the moment that the resource is completely consumed, cells are diluted to a fixed population size and new resource is added (for simulations, a cut-off at very low-substrate levels is used, after which the species densities hardly change, and the case of substrate remaining at the transfer is discussed in Appendix S1: Section A.3.2). This setup relates to the conditions that are commonly used for experimental evolution. We focus only on the exponential and nutrient-limited growth phases and ignore lag phase and stationary phase. We assume Monod-growth, because this captures the main features of growth at high-and low-substrate concentrations, and is also the most complex description for which we can derive an analytical equation of the invasion fitness. Moreover, growth properties from more detailed models can be fit well to a Monod equations, making sure that our results will be confirmable with those models. This leads to the following dynamical system: where x is the biomass concentration, s the substrate concentration, K S the Monod constant (the substrate concentration at with half maximal growth is reached) which we also refer to as the saturation constant and max the maximal growth rate. Y X∕S is the numerical yield of the species x on the substrate s. To enhance readability for the derivations, we will use a notation where we write s for s(t), x for x(t), for max , k for K S and y for Y X∕S . When we simulate multiple strains, we can extend these equations (see Appendix S1: Section A.7).

| Selection in fluctuating environments
We calculate the invasion fitness as the ratio of the Malthusian parameters (Lenski et al., 1991;Wiser et al., 2013) of the mutant and the resident minus one, such that positive values correspond to a mutant that can invade, and negative values to one that cannot (see Appendix S1: Section A.2). We obtain the following invasion fitness: where mutant trait values are denoted with an * and initial and final concentrations of a single transfer are denoted by subscripts 0 and ∞ .
This invasion fitness can directly be applied to experimental design and replace numerical simulations to estimate the relative selection on different growth traits (Vasi et al., 1994). Some direct observations from the invasion fitness are: (1) There is no selection on the numerical yield-although the resident yield does influence whether other types can invade-and (2) higher initial substrate leads to stronger selection on the maximal growth rate, while higher initial biomass leads to stronger selection on the Monod constant. The derivation of this equation, calculation of the selection gradient and more insights on the parameter dependencies can be found in Appendix S1: Sections A.3 and A.5. (which corresponds to a low saturation constant K S ) benefits growth at low-substrate concentrations. The frequency-dependent selection can lead to stable coexistence of two fixed phenotypes (without considering evolution), which we will further refer to as ecological coexistence (Egas et al., 2004).

| Evolutionary coexistence
We study evolutionary changes, which could potentially destabilize ecological coexistence, with evolutionary stable strategies (ESSs), that is a strategy ( , k, y) that cannot be invaded by any other strategy (∄ ( * , k * , y * ): > 0). In adaptive dynamics, ecological and evolutionary timescales are separated. In this system, the ecological timescale coincides with the dynamics over the course of many transfers, and the evolutionary timescales with the introduction of new phenotypes. Our system includes a third timescale, the fluctuation timescale, which encompasses the changes within a transfer (see Figure 1a). This extra timescale makes the ecological steady state a stable fluctuation, instead of a, in adaptive dynamics more common, steady state.
Using the invasion fitness, we can investigate the evolutionary dynamics using techniques from adaptive dynamics (Geritz et al., 1997(Geritz et al., , 1998Metz et al., 1996;Stefan et al., 1999): a pairwise invasibility plot (PIP) shows which mutants can invade in which resident populations. By mirroring such a PIP in the diagonal and overlaying it with the original, a mutual invasibility plot (MIP) is created, showing the area of ecological coexistence. When we include the directional selection on the traits of both residents in the area of coexistence, we create a trait evolution plot (TEP). In this way, it can be seen whether there is an attracting state in which two residents coexist (i.e. when the arrows are pointing towards the coexistence state) (Geritz et al., 1998). In a conceptual way, we show that when the trade-off between maximal growth rate and affinity is completely convex (Figure 2a), there is only a small area of mutual invasibility between two strains with different properties (Figure 2c). In this area, there is ecological coexistence, but, as the arrows show, evolution will destabilize this coexistence and the evolutionary outcome will be a single stable strategy (see Appendix S1: Figure A.11 for the derivation of the arrows). However, a trade-off with a concave part ( Figure 2d) can lead to a larger area of mutual invasibility with an evolutionarily steady coalition (i.e. a state with more than one phenotype that cannot be invaded (Metz, 2012), also called an ESS stable evolutionarily singular coalition (Geritz et al., 1998); Figure 2f and Appendix S1: Figure A.11). When the coexistence is stable even under evolutionary changes, we define it as evolutionary coexistence (termed 'evolutionarily stable coexistence' by Li et al., 2020).
There is no evolutionary branching, which means that this coexistence cannot arise by small mutational steps from a single ancestor, and has to arise from large mutational steps or by migration.

| Case study: Trade-off in growth parameters for Escherichia coli and Saccharomyces cerevisiae
We have included observations from E. coli and S. cerevisiae to study whether the trade-offs between growth at different resource levels of organisms fall in the category of a 'convex' trade-off or a 'partly concave' trade-off, and to study the resulting evolutionary dynamics.
We have collected experimental data from the literature. A previous study on this trade-off in E. coli already supplies an equation for this trade-off curve (Wirtz, 2002) (see Appendix S1: Section A.4.1). A study where chimeric glucose transporters are built into S. cerevisiae gives data for cells that mostly similar, but differ in their affinity and maximal growth rate (Elbing et al., 2004). We have used that data to generate an equation for the trade-off curve (see Appendix S1: Since the direct experimental data are limited, we have found a method to derive the trade-off from computational models ( Figure 3). The method includes starting with a metabolic or growth model, which is based on measurements of enzymes and cellular properties, and optimizing investment in different enzymes/modules at different substrate concentrations, which simulates evolutionary adaptations to different conditions. Next, the investments are fixed and adapted models are simulated to obtain the growth at a range of substrate concentrations. A Monod curve is fitted to obtain a maximal growth rate ( max ) and saturation constant (K S ) for every adapted model. Finally, an evolutionary trade-off curve is fitted through the max -K S pairs.
We have used the above methodology to create an evolutionary trade-off for E. coli, using a mechanistic computational model of central carbon metabolism that is based on known stoichiometries, | 625 WORTEL measurements of enzyme properties and measured growth rates (Wortel et al., 2018) (see Appendix S1: Section A.4.2). This model can be adapted to anaerobic conditions; therefore, we obtain two different curves. For S. cerevisiae, we use a coarse-grained metabolic model (Wortel et al., 2016) that uses measured enzyme weights, cell composition and is fitted to measured growth rates and the switching point to overflow metabolism. With the methodology outlined in Figure 3, we obtain a trade-off for this model (see Appendix S1: Section A.4.4). As S. cerevisiae already uses anaerobic fermentation at high-substrate concentrations when oxygen is present, the anaerobic and aerobic trade-off curves are very similar (data not shown), and only one generic curve is used.
In some cases, the trade-off is related to transport kinetics (Elbing et al., 2004), but mostly it is a property of the whole cell and surfaces even with a single type of transporter. The trade-offs for the same species resulting from different sources lead to different trade-offs (see Figure 4a). However, the shape of the curves does show some consistency, and resembles the 'convex' trade-off In the ESS, no other type can invade (inset in Figure 4b). We can overlay this PIP with the PIP reflected in the main diagonal (Geritz et al., 1998), to see when two types can mutually invade (green area in 4c), showing the area of ecological coexistence.
The obtained analytical expression (Equation 2) for the invasion fitness allows us to investigate many different scenarios. However, for all trade-offs, from experimental and computational data, and for different initial substrate and biomass concentrations, the general picture is similar: There is a small area of ecological coexistence, but analysis of the shape of the mutual invasibility plot shows that this coexistence is not evolutionary stable and the only attracting point is an ESS (see Appendix S1: Section A.6). This means that two phenotypes adapted to different conditions, such as a high-nutrient The dynamics of serial transfers. On fast 'fluctuation timescales', the phenotype densities and metabolite concentrations are continuously changing. A stable pattern can arise over ecological timescales, where the concentrations at the beginning of each transfer remain constant. Over evolutionary timescales, phenotypes can evolve properties such as maximal growth rate and affinity.
Invading mutants either displace one of the phenotypes and form a new 'ecologically stable' coexistence (as in the figure) or take over the whole community. (b) Negative frequency-dependent selection between low-and high-substrate (S) specialists (arising from e.g. different protein allocations; investment in ribosomes for higher max and in transport protein for higher affinity (lower K S )). (c) When low-substrate specialists (LSS) are abundant, the substrate concentration decreases slower but is completely depleted faster than when high-substrate specialists (

F I G U R E 3
Methodology to obtain evolutionary trade-offs from computational models.
Step 1: Starting point is a detailed metabolic model or a coarse-grained growth model that is based on experimental data (e.g. enzyme properties, fluxes and growth measurements).
Step 2: The next step is to optimize the model for growth rate at different substrate concentrations by adjusting the investment in enzymes or modules. Three examples are shown: red/arrow up (optimized at high substrate; investment mainly in ribosome), purple/equal symbol (optimized at intermediate substrate) and blue/arrow down (optimized at low substrate; investment mainly in transporters).
Step 3: Then, the investments remain fixed, but the growth rate of the different optimized models is calculated at different substrate concentrations and a Monod curve with a max and K S is fitted for each model.
Step 4: The max is plotted versus 1 ∕ K S for each model (the procedure of steps 2 and 3 is repeated for more substrate concentrations than the three that are shown as examples) and then an evolutionary trade-off curve, which described a relation between max and K S , is fitted through all of those points (solid multicoloured line).

| 627
WORTEL and low-nutrient condition, could coexist when placed together in a fluctuating environment. However, when they adapt in the same fluctuating regime, they will evolve in such a way that one of the phenotypes will be lost, and a single phenotype will remain.

| Evolutionary coexistence in multiple species consortia
While we have not predicted evolutionary coexistence for strains within a single species, specialization on resource level could explain coexistence between species. To study evolutionary species coexistence, we need one species to be in an ESS when the other species has a certain phenotype, and then, that phenotype has also to be an ESS in the ESS phenotype of the first species. Here, we present a method to find such a set of phenotypes. This method allows us to study evolutionary coexistence analytically, thereby being able to check many different scenarios, while still guaranteeing that the found solutions will also be a solution of the complete two-species dynamical system. We can simplify the two trade-offs into one (that possibly is 'partly concave') to find this set of ESS phenotypes, because a species with a higher affinity or higher maximal growth rate for the same value of the other property will always outcompete the species with a lower value, and therefore, we can ignore the inferior phenotype.
As a proof of principle, we combine the two trade-offs that we found to intersect, namely the E. coli trade-off derived from the computational model under anaerobic conditions and the S. cerevisiae trade-off derived under experimental conditions. Although the experimental conditions were aerobic, we confirmed with the computational model that for S. cerevisiae, the effect of oxygen is not very strong. Therefore, we deem this artificial selection of trade-off curves reasonable for a proof of principle and illustration of the method and its applicability for biologically realistic estimates. First, we construct a new trade-off that combines the two (see Figure 5a and Appendix S1: Section A.4.7). The intersection if the trade-off curves is in a fitness minimum, therefore crossings of this intersection do not happen. Such crossings are not biologically relevant in the current setting as they would imply a mutation from one species to the other, but could be relevant when studying closer relatives and/or longer evolutionary timescales.
Allowing for large mutations does not affect our conclusions except that the evolutionarily steady coalition is that reached from any initial (set of) phenotype(s).
The analytical approach allows us to screen different conditions (initial substrate and biomass concentrations) and find scenarios of coexistence (trait evolution plots resembling Figure 2d).
When we construct the pairwise invasibility plot for the twospecies trade-off, there are three equilibria, one local minimum at the species intersection and two local maxima, one for each species (Figure 5b). The single species maxima are both ESSs and will be reached when we start with a single species (Figure 5c and Appendix S1: Figure A.14). There is an area of mutual invasibility between the species, shown as the green area in Figure 5c. If both species would be adapted without the other species under the same fluctuation regime, they would coexist when one migrates into the other community. Moreover, using the analysis of trait evolution plots (Appendix S1: Figure A.11), we can show that there is an evolutionary attracting point in the species coexistence area (Stefan et al., 1999). With simulations (Appendix S1: Figure   A.14), we can find this point and we show numerically that in this point, the invasion fitness for all other types is negative (inset in Figure 5c), meaning that both phenotypes are in an ESS. This means that after the species coexist together, they will continue to F I G U R E 4 Evolutionary trade-offs and ecological coexistence for different phenotypes within a species. (a) Trade-offs between ( max ) and the substrate affinity (1 ∕ K S ) derived from experimental data or experimentally parameterized models (see Appendix S1: Section A.4).
(b) Pairwise invasibility plot for the experimentally derived trade-off for E. coli at the conditions of the LTEE (initial glucose = 0.1375 mM, initial cells are 0.5 × 10 6 cells/mL) (Vasi et al., 1994). Black area depicts where the mutant can invade and the axes run from low-substrate specialists (low saturation constant and low maximal growth rate) to the high-substrate specialists (high saturation constant and high maximal growth rate). The inset shows that the equilibrium is invasion stable. The equilibrium is also attracting and therefore an evolutionarily stable strategy (ESS). (c) Mutual invasibility plot. By reflecting the PIP of (b) over the main diagonal and overlaying it with the PIP, we can determine the area of mutual invasibility and therefore ecological coexistence (green area) (Geritz et al., 1998).

(a) (b) (c)
evolve, but reach an evolutionary steady coalition. There is a range of initial biomass and resource concentrations that supports this coexistence (Appendix S1: Section A.6.1).

| DISCUSS ION
Our main result is that specialization on resource level is likely to support species coexistence in fluctuating resource conditions. Nutrient fluctuations are abundant in nature and are often mimicked in laboratory evolution studies with a serial transfer regime.
In this paper, we derive an analytical expression for the invasion fitness in such regimes, develop a method to use computational models to obtain trade-off curves between growth parameters and apply adaptive dynamics to study the evolutionary dynamics. From the limited available data, we do not find any evidence for evolutionary stable coexistence by specialization on resource level within a species, making this less likely to be an explanation for the observed diversity in laboratory evolution studies (Good et al., 2017). However, it could be a factor contributing to coexistence, such as in Turner et al. (1996), where this specialization is combined with cross-feeding. We show with a proof of principle that the methods can in this case be applied to species coexistence, and although we do not have data on naturally co-occurring species, from our results, we postulate that evolutionarily stable coexistence by specialization on resource level is much more likely to occur.
Although we have not been able to show it rigorously, the evolutionary stable coexistence depends on the shape of the trade-off curve (see Figure 2). Combining the two species leads to a 'partly concave' trade-off curve (see Figure 5). In theory, a single species could have such a shape trade-off curve. We have not found any evidence for this, but as our data here are limited, phenotype and fitness because of the frequency-dependent selection (Rueffler et al., 2004). The distinction between species and strains is made here on the bases of the trade-off curves.
Strains are variants that are evolutionarily constrained by the same trade-off curve, while different species follow different trade-off curves. This could be interpreted as species being further apart in phenotype space than strains. We do acknowledge that other properties will also affect the trade-off curve, and trade-off curves themselves can evolve over longer timescales (Roff & Fairbairn, 2007).
Our results give indications which conditions are most promising to lead to coexistence. Evolutionary coexistence between species is most likely when dilution is not too extreme, and relatively high cell densities remain before nutrients are replenished (see Appendix S1: Section A.6.1). This is because high dilution leads to a long growth phase, which benefits a fast growth strategy, making it difficult for low-substrate specialists to survive (this was also found in Abreu et al., 2019). Natural habitats such as the mammalian gut satisfy such conditions. When experimental data are available for species and conditions of interest, our methods can be used for a thorough parameter scan of substrate and biomass levels and fluctuation F I G U R E 5 Evolutionary coexistence between species. (a) Combined trade-off of two species between the maximal growth rate ( max ) and the affinity (1 ∕ K S ). (b) In the PIP (with initial biomass x 0 = 10 6 cells/mL and initial glucose s 0 = 500 mM), we can see there are three equilibria, of which the middle one is a local minimum corresponding to the switching point of one species to the other (grey dashes lines; those cannot be crossed by mutations). There are two equilibria that are attracting and locally invasion stable, which correspond to the ESS for each species. (c) Three different simulations are plotted onto the area of coexistence, created with the mutual invasibility plot (see Appendix S1: Figure A.14 for the simulations). When initially one species is present, a single species ESS is reached. However, when we start with two species within the area of coexistence, an evolutionarily steady coalition (ESC) is reached. As shown in the inset, in this ESC, all other types have a negative invasion fitness, and therefore, this ESC is evolutionarily stable. | 629 WORTEL lengths that support coexistence (see Appendix S1: Sections A.3.2 and A.6.1).
We have used a trade-off between maximal growth rate and substrate affinity to determine the possible trait space of the organisms. Data and detailed models that can infer these properties are not abundant; therefore, we focussed on E. coli and S. cerevisiae, where data are still limited and not completely consistent (as can be observed from the difference between the curves obtained from different data sets in Figure 4a, indicating that more study is needed on these trade-offs.). We accept that these might not be the most biologically relevant species for coexistence studies, but if we can already find possibilities for coexistence in this very limited data set, it is likely to occur for some combinations that cohabit in nature.
We have shown some support for the trade-off in these species; however, Velicer & Lenski (1999) set out to measure this trade-off experimentally in E. coli and found more contra-than supporting evidence. However, as they also note, perhaps their experiment was not long enough to actually reach the trade-off, and general growth performance mutations on the supplied substrate were selected for.
Moreover, E. coli might not be the organism where this phenomenon is most important, because it tends to use active transport that usually already coincides with a high affinity. This trade-off might be more relevant for other micro-organisms, for CO 2 uptake by cyanobacteria and green algae (Ji et al., 2017) or plants (see Appendix S1 in Beardmore et al., 2011).
Our results could be verified with more mechanistic models.
Our basic model uses Monod growth, which is a rather simple description of microbial growth. However, it is also the most complicated one for which we can derive an analytical expression for the invasion fitness, which is a new result obtained here. Moreover, because of the analytical results which make it easy to check many scenarios, our results can be used as a basis for the conditions under which more detailed models can be simulated. For the more detailed models used to derive the trade-offs in this paper (Wortel et al., 2016(Wortel et al., , 2018, we do expect to find the same results, because the important properties are captured in the fit of the trade-off curve. The only differences would arise from the error of the fit and from the way we would include mutations, as we have only studied adaptation on the Pareto front, that is, with strategies optimal in the max -K S domain. This does lead to two interesting directions for further investigation: Simulation with more detailed growth models, from a Hill equation to a metabolic network or a genome scale model with uptake kinetics (Mahadevan et al., 2002), and including mutants that are not directly perfectly adapted, which could be done with (agent based) simulations (e.g. Virtual Microbes; van Dijk et al., 2019).
In Figure 5c, we see bistability in the evolutionary sense, that is, the reached ESS depends on the initial traits of the species. Only when the initial traits of the two species are in the coexistence area, stable coexistence will evolve, and otherwise, the outcome will be a single species. However, the proposed model can, besides coexistence, also show ecological bistability, that is, the final outcome depends on the initial ratio of the strains. Bistability was also observed in another study with a trade-off between maximal growth rate and lag phase (Manhart & Shakhnovich, 2018). The conditions for bistability of single-strain states include mutual non-invasibility, and can therefore be studied with Equation 2. Bistability is only possible if the strain specialized on high substrate has a low numerical yield, such that when that strain is abundant, nutrients remain high, favouring itself and supplying the necessary positive feedback. However, when we look at the measured parameters, bistability is not observed (there is no overlay of non-invasibility areas in figure A.9), and therefore, initial ratios are not important for the outcomes presented here. Note that bistability is not possible when instead of the initial biomass, the dilution is kept constant, because then the yield does not affect the invasion fitness (see Appendix S1: Section A.3.1).
We have only discussed the coexistence of two strains. Previous studies have shown that theoretically, more than two species could coexist, albeit mostly though numerical simulations or only showing neutral coexistence (Manhart & Shakhnovich, 2018;Smith, 2011;Stewart & Levin, 1973 the results. If nutrients are not completely depleted, this will benefit the high-substrate specialist, and therefore, coexistence will be less likely (see Appendix S1: Section A.3.2). Second, in the coexistence, the high-substrate specialist is often more abundant than the low-substrate specialist. If there would be a relation between the equilibrium ratio for the different strategies and the location of the ESC, this could simplify finding the ESC, which is difficult because an analytical solution for the fitness with more than one resident cannot be found. However, we have not been able to find this relation.
Third, it would be interesting to apply our methods to other growth and community dynamics, such as lag phases, death in stationary phase and cross-feeding. It might be possible to assign the relative contributions to coexistence when combining different properties in one model. Last, it would be ideal to have experimental verification of these results. For that, first, models and theory can be used to identify an experimental system that would be suitable, that is, where we expect evolutionarily stable coexistence under some (initial) conditions. Then, experiments can be designed to determine if only resource-level specialization achieves coexistence, or if it is one of the contributors. As long as a suitable experimental system is lacking, virtual organisms could be a first test.
In conclusion, we have shown that we can use adaptive dynamics and trade-offs based on experimental data and computational models to show the biologically relevance of an observed and suggested mechanism behind diversity of microbial communities in fluctuating environments. As more quantitative data become available and more mechanistic models can be parameterized, we can investigate how common it is that species could coexist by specialization on resource level. Moreover, we provide an analytical expression for the invasion fitness in fluctuating environments and techniques to investigate evolutionary stable diversity within and between species, which can be applied more broadly.

ACK N O WLE D G E M ENTS
I would like to thank Han Peters for discussion on the mathematical derivations, Hans Metz for discussions on applying adaptive dynamics techniques, Guilhem Doulcier for the online availability of his lecture notes and code on adaptive dynamics analysis with Python,

Martijn Egas for thorough reading, advice and discussions and Juan
Bonachela, Bob Planqué and Huub Hoefsloot for helpful discussions.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The author declares no conflict of interests.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.14158.