The importance of species interactions in eco-evolutionary community dynamics under climate change

Eco-evolutionary dynamics are essential in shaping the biological response of communities to ongoing climate change. Here we develop a spatially explicit eco-evolutionary framework which features more detailed species interactions, integrating evolution and dispersal. We include species interactions within and between trophic levels, and additionally, we incorporate the feature that species’ interspecific competition might change due to increasing temperatures and affect the impact of climate change on ecological communities. Our modeling framework captures previously reported ecological responses to climate change, and also reveals two key results. First, interactions between trophic levels as well as temperature-dependent competition within a trophic level mitigate the negative impact of climate change on biodiversity, emphasizing the importance of understanding biotic interactions in shaping climate change impact. Second, our trait-based perspective reveals a strong positive relationship between the within-community variation in preferred temperatures and the capacity to respond to climate change. Temperature-dependent competition consistently results both in higher trait variation and more responsive communities to altered climatic conditions. Our study demonstrates the importance of species interactions in an eco-evolutionary setting, further expanding our knowledge of the interplay between ecological and evolutionary processes.


Quantitative genetic recursion in space
Let S species be distributed across L habitat patches, with possible migration in between. Each species i is characterized by its density N k i and trait distribution p k i (z) in patch k. Here z measures a unidimensional quantitative trait of interest; N k i p k i (z)dz is then the density of species i's individuals in patch k whose phenotype values fall between z and z + dz. By definition, at any moment of time. The mean trait of species i in patch k is Our starting point is the framework of quantitative genetic recursion 1-5 in the weak selection limit 6 . We extend this approach by adding migration. The basic equation giving the change in the density distribution of species i in patch k reads where the tilde denotes values after selection and migration, but before reproduction. The first term on the right hand side describes selection via the fitness function W k i (z). The second term is immigration from all other patches (therefore M kl i is the dispersal rate from patch l to patch k, with M kk i = 0 for all i and k). The final term is the loss of species i's trait z in patch k due to emigration from the focal patch.
For the underlying genetics, we assume a large number of loci each contributing a very small additive effect to the trait 7,8 . In this infinitesimal model, all trait distributions p k i (z) are normal with a variance σ 2 i that does not change in response to selection: The change in mean trait from one generation to the next is given by the breeder's equation 9 : where the hat denotes values in the next generation, andμ k i denotes the trait mean after selection and migration, but before reproduction. After reproduction,N k i =Ñ k i , andp k i is given by Eq. 4 withμ k i calculated from the breeder's equation 5 .
In the weak selection limit, the fitness W k i (z) of species i's phenotype z in patch k can be written W k i (z) = 1 + sr k i (z), (6) where s is a small parameter, and r k i (z) is the per capita growth rate of species i's phenotype z in patch k, defined by ecological interactions. The migration rates are similarly written as To obtain the dynamics of the densities, we integrate Eq. 3 across z: Using Eq. 1 andÑ k i =N k i : In the weak selection limit (Eqs. 6 and 7), which, using Eq. 1, is written aŝ Subtracting N k i and dividing by s leads tô With appropriate scaling 5 , this can be written in differential equation form when s is sufficiently small: To obtain the equation for the change in trait means, we writeμ k i by rearranging Eq. 3: We use Eq. 9 to expand the denominator: We multiply both sides by z and integrate. Using Eq. 2, we get In the weak selection limit (Eqs. 6 and 7), this reads Using Eqs. 1 and 2, we get Simplifying by N k i and factoring out s from both the numerator and denominator: Taylor expanding around s = 0 to first order, we get where O(s 2 ) means terms of order s 2 or higher. Simplification yields Substituting this expression into the breeder's equation (Eq. 5): After dividing both sides by s and performing the s → 0 limit with appropriate scaling 5 , we get the differential equation approximation

Per capita growth rates and species-level dynamics
We model a trophic community in space, with L distinct habitat patches. Each patch experiences an average temperature, which is increasing through time due to climate change. The evolvable trait z is identified with temperature tolerance: an individual with trait z achieves maximal intrinsic growth in a habitat with temperature z. Local population dynamics, excluding migration, are governed by four processes: temperature-dependent intrinsic growth, intra-and interspecific competition, growth due to consumption, and loss due to being consumed. The per capita growth rate of of species i's phenotype z in patch k is written where, in patch k, r k 0,i (z) is the intrinsic growth rate of species i's phenotype z, a k ij (z, z ) are competition coefficients between species i's phenotype z and species j's phenotype z , i is species i's resource conversion efficiency, and F k ij is the feeding rate of species i on j. It has the form where q i is species i's attack rate, W ij is the adjacency matrix of the feeding network (W ij = 1 if i eats j and 0 otherwise), ω ij is the relative consumption rate of i on j (we assume an equal split across resources, so that ω ij = [number of resources] −1 for all consumers), and H i is species i's handling time.
Substituting the growth rate in Eq. 24 into Eqs. 13 and 23 yields for the population densities, and for the trait means. Rearranging, we get By Eq. 1, the integral in the third term of Eq. 28 is simply 1. In turn, the integral in the third term of Eq. 29 is zero, because the integrand is a product of an odd and an even function in z − µ k i . This leads to Introducing the definitions Eqs. 30-31 read with F k ij given by Eq. 25.

Model parameterization
Here we present the detailed derivation and parameterization of each model component.

Temperature-dependent intrinsic growth rates
Following Amarasekare & Johnson 10 , we model the temperature-dependence of r k 0,i (z) using the Gaussian form where κ i is a mortality rate for species i, A k i − κ i is the maximum intrinsic growth that can be achieved for species i in patch k, T k is the temperature in patch k, and w k i is the species-and patch-specific width of the growth curve.
The w k i are not constant, but co-evolve with species' mean temperature optima µ k i . Within reasonable limits (such that the tolerance widths remain positive), this can be modeled using the following linear relationship: where a w and b w are positive constant parameters. This creates a trend of increasingly broader curves for lower temperature optima, corresponding to the observed correlation between tolerance width and latitude 11,12 . For determining A k i , we use the observation that there is a tradeoff between the width and maximum height of temperature optima. This phenomenon is brought about by higher annual temperature variability at higher latitudes: since reaction norms are nonlinear functions of temperature, the optimal response to the mean temperature is not the same as the mean optimal response, due to Jensen's inequality 10 . This also means that the distance between experienced and optimal temperatures is larger at higher latitudes, where climate tends to be more variable 13,14 . However, our model does not take annual temperature fluctuations into account, which means that the observed optimum shift cannot evolve as in the work of Amarasekare & Johnson 10 . Instead, we imposed this tradeoff by assigning a smaller A k i to wider temperature adaptations at low µ k i values (Eq. 39), mimicking the fact that populations at higher latitudes are adapted to more variable environments and thus also living at temperatures further away from those that would yield maximum growth in a stable environment. The tradeoff is implemented as where i is a parameter modulating the shape of the tradeoff. Substituting the parameter specifications of Eqs. 39-40 into Eq. 38 gives This form of the tradeoff, and its particular parameterization (Supplementary Table 1), were chosen so the resulting growth curves would qualitatively mimick those found empirically 13,14 . Supplementary Figure 1 shows what the growth rates of Eq. 41 look like, in six species equally spaced along a temperature gradient, both for consumers and resources.
To obtain the species-level parameters b k i and g k i from this growth rate, we substitute Eqs. 4 and 41 into Eqs. 32 and 34: The integrals can be evaluated directly, leading to the final forms

Competition coefficients
We only model competition between resource species. Consumers also compete but only indirectly, via shared resources. This means that the competition kernel a k ij (z, z ) is zero unless both i and . Colors indicate species, with cooler shades standing for more cold-adapted ones. The dashed line indicates zero growth; therefore the r k 0,i of consumer species is always a mortality rate, which must be compensated by consumption. For resources, r k 0,i is positive unless the local temperature is too far from the optimum for that species. In this figure, κ i = 0.1 for all species, i = 0.1 for all consumers, and i = 0.1 for all resources. In our simulations, to make sure that the model's predicted community patterns are not an artifact of the above precise tradeoff shapes, the i were randomized in a ±10% range of these values (Supplementary Table 1).
j refer to resource species. When they do, we use two different kernel forms. One only depends on species identity, and is independent of trait or patch value: Applying Eqs. 33 and 35 to obtain α ij and β ij , we get In the first equation, the integrals evaluate to 1 because of Eq. 4, while in the second equation, integration with respect to z involves the product of an odd and an even function in z − µ k i and is therefore zero. The constant a ij values are sampled randomly (Supplementary Table 1), in a way that ensures intraspecific competition being always stronger than interspecific competition.
The second form of the competition kernel we use depends on the difference between phenotypic values: where η is the competition width, determining how distant two phenotypes must be for competition to be significantly reduced between them. Here α k ij and β k ij can be calculated by direct integration using Eqs. 33 and 35: This form of the competition kernel facilitates trait divergence even within a single patch, since it may be profitable for a species to specialize on a suboptimal temperature to thus weaken interspecific competition and coexist locally. Biologically, this mechanism is explained by having microhabitats within each patch, with higher and lower local temperatures than the patch average. A species with a higher or lower temperature optimum than the average can exploit these microhabitats better than the locally dominant species, and therefore persist in the patch as a whole-though it will not be able to achieve as much growth due to the relatively limited availability of such microhabitats.

Feeding network
We model two strict trophic levels, with S resource and S consumer species. The bipartite network of feeding connections W ij (with W ij = 1 if i eats j and 0 otherwise) is generated as follows. First, both resources and consumers are labeled consecutively, based on their initial temperature adaptations: resource 1 / consumer 1 are the most cold-adapted, and resource S / consumer S the most warm-adapted. Next, we always put a feeding link between consumer i and resource i. Finally, each consumer is randomly linked to four other resource species. This yields a feeding network where every consumer is connected to five resources altogether.

Genetic and environmental variances
Genetic variances V G,i are randomly drawn for each species, while the environmental variances V E,i are fixed (Supplementary Table 2). We assume that there is no epistatic or dominance variance, which means that the total phenotypic variance, σ 2 i , is the sum of the additive genetic and environmental variances: The heritability of species i's trait value is, by definition, the ratio of genetic to total phenotypic variance 9 , so we have Our model has no demographic stochasticity. This means we are assuming that even "low" population densities are sufficiently high to neglect genetic and ecological drift. This, however, can introduce problems with artifactual evolutionary rescue, where an unrealistically tiny population evolves to the point where it can finally grow, thus surviving and establishing itself even though it was bound to go extinct. To control for this artifact, we introduced a reduction in species' genetic variances whenever their local population densities dropped below the threshold N c . Below this threshold, genetic variances were multiplied by a factor This is a smoothed step function: it is 0 for x < 0, increases monotonically to one at x = 1, and stays equal to 1 for x > 1 (Supplementary Figure 2). It is constructed to be twice continuously differentiable for all x. This means that heritabilities were also modified below N c , being instead given by the patch-dependent expression

Spatial structure and dispersal rates
We discretize latitudinal position into L equidistant patches: patch 1 is at the pole, patch L at the equator, with the rest linearly spaced in between. Migration happens between adjacent patches with a species-specific rate d i : To the approximation that the two hemispheres of Earth are symmetric with respect to their temperature profiles, the pole-to-equator spatial profile could be repeated from equator to the other pole again (in reverse), and then wrapped around the whole planet. Due to this symmetry and periodicity, the range going from pole to equator already contains all information, and is the only part that needs to be modeled, provided one implements appropriate boundary conditions to account for the periodicity. Specifically, for patches 1 and L, the following replacements are made compared to Eqs. 36-37:

Climate model
Let x denote latitudinal position, measured as the latitudinal distance from the north pole. The temperature at position x and time t is given by T (x, t). The local temperature T k in patch k is the one at the latitude corresponding to patch k. The components of the climate change model are as follows: 1. Temperature is initially constant in each patch, for a burn-in period from t 0 = −4000 to 0 years. The initial average temperature T 0 (x) is T min at the pole, T max at the equator (Supplementary Table 1), and linearly extrapolated in between: where X is the distance from pole to equator, measured in the same units as x.
2. Climate change starts at t = 0, after the burn-in period from t = t 0 to t = 0. In each patch, the qualitative shape of the relative temperature rise is given by the smoothed step function Q(t/t E ) (Eq. 53), where t E is the end of the climate change period.
3. Climate change stops after t E = 300 years.
4. The extent of the temperature increase depends on latitude. This latitudinal trend is predicted with high confidence for the northern hemisphere for the next century, and is observed in all of IPCC's emission scenarios 15 . We assume that the latitude-dependence of the final temperature increase (at t E ) is linear: where C max is the maximum temperature increase (at the pole) and C min is the minimum increase (at the equator). The values of C max and C min (Supplementary Table 1) are based on region-specific predictions of increase in temperature by 2100 CE, in combination with estimates giving approximate increase by 2300 CE, for the IPCC intermediate emission scenario 15 .
Considering all the above, the time-dependence of the temperature profile T (x, t) is written compactly as

Initial conditions
Both for resource species 1, 2, . . . , S and for consumer species 1, 2, . . . , S, the initial temperature adaptations µ k i (t 0 ) are determined by the initial temperature profile such that species 1 is adapted to the polar temperatures, species S to equatorial temperatures, and all other species are linearly spaced in between: at all spatial locations where the species is initially present. Initial population densities are centered at the location to which the species is best adapted, and are normally distributed around this point in space:

Community response capacity
Let there be S species distributed across L patches, evolving according to our model. We denote the relative density of species i in patch k at time t by n k i (t): and the community-weighted mean trait in patch k at time t byμ k (t): We then define the trait lag in patch k, A k (t), as where T k (t) is the local temperature in patch k at time t. This quantity measures the extent to which traits in a local community match their environment, on average. A k (t) may be further averaged across patches and time, yielding The time averaging moves from the onset to the end of climate change, i.e., between years 0 and t E = 300. In turn, we define the density-weighted dispersion of mean trait values, V k (t), as The more different species are in their trait means (with relatively more common species getting more weight for being different), the larger this value becomes. Again, we can average this quantity across time and patches to obtain A and V can be calculated for each community and plotted against one another, as was done in Figure 7 in the main text.

Trait dispersion and species richness
Regardless of the model analyzed, higher species richness results in higher trait dispersion overall (Supplementary Figure 4). However, this global trend is sometimes locally reversed for particular parameterizations-namely, when genetic variance is low and competition is temperature-dependent. This creates a version of Simpson's paradox 16 , where the aggregated data show different trends than disaggregated pieces of it. The extent of trait dispersion strongly depends on the mode of competition between resource species: if competition is temperature-dependent (Supplementary Figure 4, bottom two panels), trait dispersion is about an order of magnitude higher than otherwise (top two panels; note the difference in scale along the abscissas). This is in line with the fact that, with temperaturedependent competition, species can coexist by being sufficiently different from one another, by adapting to locally suboptimal temperatures which nevertheless reduces competition to tolerable levels. Without temperature-dependent competition, whether species adapt to have similar temperature optima has no consequence to how strongly they compete. Thus, if a set of species can coexist, they can do so at the local trait optimum. The presence of many species at the optimum naturally results in reduced trait dispersion (Eq. 69).
With temperature-dependent competition, high spatial dispersal ability results in higher trait dispersion (Supplementary Figure 4, bottom panels, blue and yellow). High genetic variance and low dispersal reduce both richness and dispersion in all models: the lack of dispersal means that each local community is more independent from its neighbors, at which point the locally readily adaptable species organize themselves into a stable community relatively quickly. This is known to lead to lower richness than equivalent communities without the ability to adapt-though these lower-richness communities are also more robust against environmental perturbations 5 . However, high genetic variance and dispersal ability versus low variance and dispersal shift position depending on whether competition is temperature-dependent.

Community turnover
Species turnover in time is measured by the Jaccard distance. Let S(t) denote the set of species present at time t, and |S(t)| the number of elements in the set. The Jaccard distance J(t 0 , t) between the community state at some reference time t 0 and a later time t reads where ∩ is the intersection and ∪ the union of two sets. This formula is agnostic about the spatial extent of a community: it can be applied to each patch separately, or the community as a whole. An in-between approach is to obtain regional turnover by applying Eq. 71 to the top third of all L patches (the polar region), the middle third (temperate region), and bottom third (tropical region). Supplementary Figure 5 shows how regional species compositions change in  time, with the community state at the start of climate change playing the role of the reference community at t 0 .

Beta diversity and rank-abundance curves
Beta diversity, the turnover in species diversity when moving from one local site to another within a larger regional community, is calculated following Whittaker's approach 17 , where β = γ/α.
Here α accounts for local species diversity, and γ is the total species richness over all patches. The former, α, is calculated as in Figure 4 of the main text, with species mean per patch over the landscape. In turn, we obtain β as in the global patterns of Figure 6 of the main text. Supplementary Figure 6 shows how beta diversity changes over time for the resource species.
To display species' relative abundances, we obtained rank-abundance curves for resource species 7. The calculation is based on species' mean abundances over all patches and replicates, at the end of our simulations (t = 2500). Species are then ordered from highest to lowest abundance.

Generalist and specialist consumers
In our model, a consumer always feeds on the resource species with an initial corresponding temperature adaptation, and four additional randomly chosen species (Section 3.3). Here we relax this assumption and show result for generalist and specialist consumers. In the specialist case, each consumer feed on one single resource species (the one with corresponding initial temperature adaptation). In the generalist case, each consumer feeds on their matching resource, and is also randomly linked to S/2 − 1 other resources (rounded if necessary). For S even, this will lead to the feeding network having a connectance of 1/2.  Figure 5: Jaccard distance between the regional community composition of resource species at the start of climate change and subsequent times, in steps of 100 years (points; dashed lines connect them for better visibility). A value of 0% means that the same set of species is present as at the start, while 100% indicates a complete turnover in regional species composition. The Jaccard distance is averaged over patches (merged into regions, indicated by color) and replicates. Columns show different ecological models, rows various parameter combinations of average dispersal ability and available genetic variance.  Figure 6: Beta diversity as measured by the ratio of regional to local species diversity, at the start of climate change and subsequent times in steps of 100 years (points; dashed lines connect them for better visibility). A high beta diversity indicates that the difference between local and global species diversity is large. A low beta diversity indicates that species are evenly distributed across the landscape, with no great difference in local and global species diversity. The beta diversity is averaged over replicates. Columns show different ecological models, rows various parameter combinations of average dispersal ability and available genetic variance.
The patterns of our original setup are broadly retained with generalist and specialist consumers ( Supplementary Figures 8-13). The difference is that greater resource specificity leads to increased diversity, both on the global and local scales. This is in line with the intuition that specialist predators help maintain resource diversity, since now each resource species is additionally regulated by an independent factor. However, even when predators are generalists, they have a positive effect on local and global resource diversity patterns-albeit this positive effect is diminished in comparison with the specialist case. different from the effect of just one or the other mechanism in isolation. We can further elucidate this pattern with a 2-way ANOVA, and obtaining confidence intervals for the effect sizes of trophic interactions alone, temperature-dependent competition alone, and their interaction, using Tukey's test. We do so separately for each combination of parameter settings (high and low genetic variance and dispersal ability) in Supplementary Figures 14-17 below.

Zero dispersal
We have also repeated our analyses setting all dispersal rates to zero. In this case, we find narrower range breadths and lower local species diversity than in the presence of dispersal ( Supplementary Figures 18-21). In terms of global species richness, global losses are mitigated by temperature-dependent competition also without dispersal. Trophic interactions do so as well, but only when available genetic variance is limited (slow evolution). By contrast, local diversity is promoted by trophic interactions alone but mostly for high levels of available genetic variance. Both for local and global diversity, temperature-dependent competition has a stronger positive effect on diversity than trophic interactions alone. Isolating patches by setting dispersal to zero allows us to contrast our results with some earlier ones, obtained for trophic dynamics under climate change but without dispersal. Similar to the findings of Osmond et al. 18 , consumer presence in our model reduces resource densities (Figure 3). With consumer presence but without dispersal, we observe slightly boosted local resource richness and increased ranges, in comparison to what we see without consumers (Supplementary Figure 18-19). Contrasting with Osmond et al. 18 , we observe a similar but marginal effect of consumer presence, with slightly fewer global losses when evolution is slow (Supplementary Figure 20).
Cortez & Yamamichi 19 focus on how evolution in prey, predators, or both, affect responses of predator populations to an increase in mortality. Both synergistic and antagonistic effects between prey and predator evolution appear, but compared to their scenario with no evolution at all, evolution of one or both species increases the predator mortality extinction threshold in most cases. What we observe is that higher levels of available genetic variance (leading to faster evolution) actually slightly decreases consumer persistence (Supplementary Figure 21). While this appears to contrast with the observed increase in mortality thresholds by Cortez & Yamamichi 19 , they in fact have shown that consumer density can temporarily increase as a function of increasing mortality (hydra effect). For this reason, the outcome of increasing the speed of evolution is not immediately obvious in our multispecies setting. The fact that we find a slight reduction in consumer persistence is therefore not necessarily inconsistent with their findings.

The effect of dispersal kernels
For the main results of this study, we used a simple, symmetric dispersal kernel in which species migrate with a given rate to the adjacent patches in both a southward and northward direction (Eq. 55). We also repeated our simulations using a Gaussian, symmetric dispersal kernel given by exp(−x 2 /(2σ 2 )), where x is the distance between two patches and σ = 0.02 2/π is chosen to keep the total rate of migration the same as with the original, nearest-neighbor dispersal. Our qualitative conclusions are robust to this change of dispersal mode ( Supplementary  Figures 22-24). We additionally explored asymmetric dispersal, by modifying Eq. 55 to introduce a northward or southward bias in movement: 5d i in the latter case, where d i is the original dispersal rate of the symmetric dispersal. Overall, asymmetry does not influence the community patterns much, independent of whether we have a northward bias (Supplementary Figures 25-27) or a southward bias (Supplementary Figures 28-30).

Competitive equivalence
We also performed a robustness analysis, contrasting our results with a model where species are competitively equivalent. This means that if species also had the same intrinsic growth rates, they would have equal fitness under all circumstances, and so neutral ecological drift would be the only process changing relative abundances. To create such a model, we have re-run our simulations with all intra-and interspecific interactions equal (species still retain the ability to adapt differently to temperature). We used a value of 0.2 for species' intra-and interspecific competition, since this value is the mean of the range of intraspecific competition coefficients used otherwise (Supplementary Table 1). As inclusion of both a second trophic level and temperature-dependent competition would break competitive equivalence, we only simulated our baseline model with these equal competition coefficients. While the results are qualitatively similar to our original baseline model ( Supplementary  Figures 31-33), an important difference is that competitive equivalence leads to a 25-50% reduction in local species diversity. In comparison with our other models, with a second trophic level and/or temperature-dependent competition, local diversity is approximately 50-80% lower. However, global species richness is slightly higher in the temperate and tropical regions, in particular when evolution is slow. Overall, temperature-dependent competition mitigates global extinctions significantly (global losses of 20-25%), when compared with the baseline model with competitive equivalence (global losses of 40-50%).    in the main text assume 50 species per trophic level. This number is sufficient for all relevant community patterns to stabilize, such that adding more species has no more substantial effect. To show this, we display species ranges, local richness, and global richness for S = 30 species per trophic level ( Supplementary Figures 38-46). The qualitative patterns are identical to those in Figures 4-7 of the main text, and Supplementary Figures 5 and 37.

Increasing the number of habitat patches from 50 to 100
Beyond a certain density of patches between pole and equator, results no longer change by adding even more patches. Our choice of L = 50 patches is sufficient for achieving this kind of convergence. To show this, Supplementary Figures 47, 48, and 49 below reproduce the main text's Figures 4, 5, and 6, respectively-but with 100 instead of 50 equally spaced patches. We generated only 20 instead of 100 replicates; nevertheless, the results are visibly unchanged by this doubling of the number of patches.