Eco-evolutionary dynamics lead to functionally robust and redundant communities

Microbial communities are functionally stable and taxonomically variable: species abundances fluctuate over time and space, while the functional composition is robust and reproducible. These observations imply functional redundancy: the same function is performed by many species, so that one may assemble communities with different species but the same functional composition. The clarity of this observation does not parallel with a theoretical understanding of its origin. Here we study the eco-evolutionary dynamics of communities interacting through competition and cross-feeding. We show that the eco-evolutionary trajectories rapidly converge to a “functional attractor”, characterized by a functional composition uniquely determined by environmental conditions. The taxonomic composition instead follows non-reproducible dynamics, constrained by the conservation of the functional composition. Our framework provides a deep theoretical foundation to the empirical observations of functional robustness and redundancy.

Introduction 1 two survives when competing. We will use the word 'strain' to identify a group of individuals with equal resource 66 preferences and intrinsic fitness. 67 Resource dynamics are described explicitly. Resources are introduced in the system with a resource-specific rate h i and consumed by the individuals present in the community. Their concentrations decrease because of consumption but also vary due to cross-feeding. A fraction 1 − � of resources consumed by each individual is used for growth, while a fraction � is transformed into different resources and released again in the environment [8]. The cross-feeding matrix, with elements D ij , specifies the relative rates of resource transformation (see Materials and Methods for its parameterization). The per-capita growth rate of a strain µ is a function g µ (c) of the resource concentration c, which, in turn, depends dynamically on the population abundance because of consumption and cross-feeding. We consider the following choice The In our framework, both resource preferences and intrinsic fitness values are subject to mutations and evolution. We 77 consider different implementations of the mutational steps (e.g., including different scenarios for the relative rate of 78 Horizontal Gene Tranfer, see Materials and Methods) which, however, do not affect the results. The timescales between 79 two successive successful mutations is comparable with the ecological timescale, set by the ecological dynamics. We 80 assume that a mutation of the resource preference always corresponds to a mutation of the intrinsic fitness, which 81 is drawn at random from a fixed distribution with width �. The parameter � sets the typical difference of intrinsic 82 fitness values between two individuals. We focus on the case of small fitness differences, and we extensively explore 83 the effect on the eco-evolutionary trajectories of increasing the value of �.  The final community structure is remarkably simple if, instead of analyzing strain abundances, we focus on its 88 functional composition. We define functional occurrence F i as the fraction of individuals able to grow on i (i.e., with 89 a i = 1). After a short transient, the functional occurrences and the total biomass N relax to the respective stationary 90 values F * i and N * , which are very reproducible across different realizations.

91
Two phases characterize therefore the eco-evolutionary dynamics. The first one is an initial-condition-dependent 92 transient, where the community structure is mainly shaped by rapid invasions. In the second phase, conversely, the 93 community has converged to a stable functional composition, which we will refer to as "functional attractor" in the 94 following, and slowly evolves reaching the final strain-level equilibrium.

95
The total biomass converges to a constant value during the second phase of the eco-evolutionary dynamics, which 96 affects therefore only the relative abundance of strains. The sequence of invasions and extinctions of strains is 97 determined by the interplay of fitness differences � µ and niche differences, related to the dissimilarity of the resource 98 preferences. Importantly, the trajectories of strain abundances are effectively restrained to occur on the low(er)-99 dimensional space determined by the constraint enforced through the functional occurrences F * i . In this second phase 100 the community has thus reached a "functional maturity" and the subsequent evolution only affects strain composition 101 while leaving unaltered the functional one.

102
The stability and reproducibility of the functional attractor suggest that it is possible to predict analytically its 103 properties. We considered a toy model of the eco-evolutionary dynamics which aims at mimicking the effective 104 exploration of the phenotypic space performed by mutations. In particular, we consider only the ecological dynamics, 105 initialized with an infinitely large species pool, which encompasses all the possible strains (i.e., the 2 R possible resource 106 preferences). A similar approach has been considered to study a simpler version of the model [22] (corresponding to 107 the limit χ � 1 and no cross-feeding). The toy model further postulates a timescale separation between resource and 108 population dynamics [22, 23], which is not assumed in the full eco-evolutionary dynamics.

109
The consumer-resource-crossfeeding model with infinite pool of diversity and no intrinsic fitness differences can be analytically solved. In the Materials and Methods we show that the stationary functional occurrences F * i and the total biomass N * are given by and . ( The parameter h ef f i is the effective resource inflow in the system, which is given by the combination of resources that 110 are externally supplied and the ones produced via cross-feeding. This quantity is in simple linear relation to the inflow 111 rate of externally provided resource h i through the cross-feeding matrix D (see Materials and Methods).

112
The analytical calculations are based on many simplifying assumptions (infinitely large pool of diversity, not explicit 113 resource dynamics, absence of fitness differences) which do not strictly hold for the more complex setting of the eco-114 evolutionary model. Nevertheless, Figure 2 shows that the predictions of eq. 2 and eq. 3 accurately describe the 115 outcomes of the eco-evolutionary dynamics.

116
Resources can be partitioned in two groups based on their effective influx rate h ef f i . If the influx rate is larger than 117 a critical value h ef f c , then the ability to metabolize that resource is a "core" function, shared by all the individuals in the community (i.e., F * i = 1). The value of h ef f c depends on both the spread of the effective influx rate (the variability 119 among the h i ) and the metabolic cost χ. The higher the metabolic cost and the variability, the higher the critical 120 influx rate threshold h ef f c and, consequently, the fewer the core resources.
attractor of the eco-evolutionary dynamics. While intrinsic fitness differences are small, they are not negligible as 156 they determine the taxonomic composition within the functional attractor. Variation of the intrinsic fitness leads to 157 functional redundancy, high taxonomic variability with conserved functional profile.

158
The assumption of small intrinsic fitness differences is critical for the observations of functional robustness and 159 redundancy. Increasing the magnitude of fitness differences affects also the functional profile. Typical intrinsic fitness 160 differences of 1% do not alter substantially the functional composition (see Figure S1). Larger differences (of the order 161 of 10%) disrupt the structure of the functional attractor: the functional composition is determined by the resource 162 preferences of the individuals with the largest intrinsic fitness and the functional profile becomes largely decoupled 163 from the resource input. Differences of 0.1% and smaller are indistinguishable from the analytical prediction and the 164 functional profile closely matches the one predicted by the resource influx rates.

165
The existence of these different regimes, where the functional composition is or is not affected by intrinsic fitness 166 differences, strictly relates to the identification of the limiting factors shaping the communities. As mentioned before, 167 in fact, the intrinsic differences could be due to abiotic factors, but also to limiting factors (such as phages) other responsible for intrinsic fitness differences (e.g. phages).

173
The metabolic trade-off is an essential ingredient of our framework. We implemented it as a fitness cost that is . This property is likely to hold more generally and not to be 189 restricted to consumer-resource systems or microbial communities. We expect that a similar approach could be 190 developed to study mutualistic communities or pathogen dynamics.
We consider a consumer-resource model in presence of cross-feeding [8], which describes the dynamics of population abundances n σ (for σ ∈ S) and resources concentration c i (for i ∈ R). Changes in population abundance are defined by where δ σ is a death term and η σ is the efficiency of the conversion of energy into biomass. E g iσ is the energy flux used for strain σ to grow from metabolite i. We can similarly define the energy flux E in iσ into a cell from resource i and the energy released in the environment by the cell E out iσ in the form of other metabolites obtained from from metabolites of type i. The associated dynamics of resource concentration c i is defined by where w i defines the conversion between energy and concentration of resource i. The function h i (c i ) specifies the 194 dynamics of resource concentration in absence of consumers. 195 We assume that energy fluxes used for growth are a constant fraction 1 − � of the total ones: energy fluxed from secreted metabolites is then given by E out iσ = � � j∈R D ij E in jσ . The cross-feeding matrix element 197 D ij defines energy conversion between resource j and resource i. Energy conservation implies � i D ij = 1.

198
The energy flux E in iσ takes the form where r i (c i ) is a non-decreasing function of the concentration of resource i and ν i is the maximal intake rate of resource

202
We consider a metabolic trade-off by assuming for death rates and yield the expression Without loss of generality, we can set the timescale τ = 1 equal for all strains, as differences in τ can be reabsorbed In the eco-evolutionary simulations, we always consider resources and populations changing over a similar timescale.

207
To make analytical progress we approximate the full dynamics with the effective one obtained by assuming timescale 208 separation -i.e. resource concentrations equilibrate faster than the changes in population abundances. We underline 209 that we assume the separation of timescales only as an approximation, for the purpose of predicting analytically the 210 outcomes of the numerical simulations, which are always obtained with explicit resource dynamics.

211
In this case, one can effectively describe the dynamics of populations as where is a Lyapunov function. With our choice for the metabolic trade-off (7), such functional can be conveniently rewritten as We then introduce the total population size N = � σ∈S n σ and define the functional abundances F i as which correspond to the fraction of individuals that are able to metabolize resource i. Interestingly, and surprisingly, when � σ = 0, the Lyapunov function can then be written as function of N and {F } alone: The fact that the Lyapunov function depends only on the total biomass and the functional profile already suggest, 212 even if it does not imply, that functional abundances are the relevant variable for the study of community composition.

213
By minimizing the function over F i in [0, 1] one obtains where the total biomass is the solution of These equation can be solve iteratively, starting from F i = 1 ∀i and N = � j∈R h ef f j /(1 + χR).

214
In the case with no intrinsic fitness differences (� σ = 0), the equilibrium solutions are identified by equations 13 215 and 14. For a given system, a fraction of resources will be core resources, i.e. shared by everyone F * i = 1. These core Eco-Evolutionary dynamics 218 The mutation probability of a preference of resource i in strain µ depends on whether µ consumes or not i. The 219 rate U −,i at which a mutantμ stops consuming resource i (the parent has a i = 1 and the mutant a i = 0) is constant, 220 independent of i, and equal to U − . The rate at which a mutant starts consuming a resource i (the parent has a i = 0 221 and the mutant a i = 1) equals to U +,i = U + (P h F i + P dn ). The quantity P h is the probability that an addition 222 happens because of horizontal gene transfer, while P dn = 1 − P h the probability of "de novo" mutations. The rate 223 of horizontal transfer is proportional to the frequency F i of that allele in the population, while the rate of a de-novo 224 mutation is independent of i.

225
The rate at which the resource preference i mutates in strain µ is then equal to where b µ is the per-capita birth rate on strain µ, which is equal to In theory one could expect a new mutant to have abundance 1. The initial phase of its dynamics is then dominated by demographic stochasticity, with many mutants going to extinction despite having a positive (average) growth rate.
In our framework, we do not consider this effect of demographic stochasticity explicitly, but we include it effectively.
Since the initial abundance of the mutantμ is a small fraction of the total population, its stochastic dynamics can be approximated by a stochastic exponential growth. In this regime, the per-capita birth rate of the mutant is given where c * i is the concentration of resource i prior to the mutant arrival. The per-capita death rate of the mutant reads (1 + χ � i aμ i )(1 − �μ). Under the assumption of a stochastic exponential growth the survival probability is given by The strain intrinsic fitness values � σ are independently drawn from a Gaussian distribution with mean 0 and standard 226 deviation �.

227
By calculating all these quantities for all possible mutations of all existing strains, one obtain the rate of invasion W inv iµ of a mutantμ which is obtained by changing the resource preference of strain µ for resource i. The rate of invasion W inv iµ reads where the mutantμ differ from µ in the resource preference i. 228 We simulate the eco-evolutionary dynamics as a sequence of discrete small time steps Δt. After a step of integration 229 of equations 4 and 5 we update the values of W inv iµ , as they depend on strain abundances, and checked whether a 230 mutant appeared. Each mutant, identified by the parent strain µ and a resource i, has probability W inv iµ Δt to invade.

231
If such an event occurs, the new mutant is introduced with an initial relative density equal to 10 −5 . The results presented in this paper were achieved using generic parameters, whose details can affect the distribution of taxa or relaxation time but not the macroscopic observables that characterize the functional attractor. In order to quantify the convergence to the functional attractor, we measure the discrepancy between the functional composition of the community during its eco-evolutionary trajectory and the functional composition predicted by equations 13 and 14. As a measure of the discrepancy, we consider the Kullback-Leibler divergence between the normalized functional profiles The  shows that the value of K i has no effect on the functional composition of evolved communities.

247
The intrinsic fitness on any new mutant was drawn from a Gaussian distribution with mean zero and variance � 2 , 248 independently of the fitness of the parent. In the main text we considered � = 0.001. Fig S1 explores  fitness. For intrinsic fitness differences with a width of the order of 10 −2 , the functional composition converges to the 252 analytical prediction, which becomes more and more accurate as fitness differences decrease.

253
The strength of cross-feeding � has no effect on convergence to the functional attractor (see Fig. S2). The cross- 259 Figure S3 shows that the outcomes of the evolutionary trajectories are independent of frequencies of the different 260 mutation steps. We varied the (average) total mutation rate U tot = (U + + U − )/2, the ratio U − /U + between mutation 261 leading to deletions of resource preferences (with rate U − ) and the ones leading to additions (U + ), and the ratio ratio U − /U + , affected the evolutionary trajectories and speed of adaptation, none of these parameters affected the 264 convergence of the functional composition to the predicted attractor.

265
Both the eco-evolutionary simulations and the analytical approximation are based on the assumption that the metabolic cost is linear in the consumed resources, as expressed in equation 7. In general, one could assume a non-linear tradeoff [24] that takes the form where g(z) is an arbitrary non-linear, monotonically increasing, function. We considered the outcomes of the evolution- show the evolution of the composition of a system with a dynamically varying noise .