The importance of species interactions in eco-evolutionary community dynamics under climate change
Modeling frameworkWe consider S species distributed in L distinct habitat patches. The patches form a linear latitudinal chain going around the globe, with dispersal between adjacent patches (Fig. 1). The state variables are species’ local densities and local temperature optima (the temperature at which species achieve maximum intrinsic population growth). This temperature optimum is a trait whose evolution is governed by quantitative genetics18,19,20,21,22: each species, in every patch, has a normally distributed temperature optimum with a given mean and variance. The variance is the sum of a genetic and an environmental contribution. The genetic component is given via the infinitesimal model23,24, whereby a very large number of loci each contribute a small additive effect to the trait. This has two consequences. First, a single round of random mating restores the normal shape of the trait distribution, even if it is distorted by selection or migration. Second, the phenotypic variance is unchanged by these processes, with only the mean being affected25 (we apply a reduction in genetic variance at very low population densities to prevent such species from evolving rapidly; see the Supplementary Information [SI], Section 3.4). Consequently, despite selection and the mixing of phenotypes from neighboring patches, each species retains a normally-shaped phenotypic distribution with the same phenotypic variance across all patches—but the mean temperature optimum may evolve locally and can therefore differ across patches (Fig. 1).Fig. 1: Illustration of our modeling framework.There are several patches hosting local communities, arranged linearly along a latitudinal gradient. Patch color represents the local average temperature, with warmer colors corresponding to higher temperatures. The graph depicts the community of a single patch, with four species present. They are represented by the colored areas showing the distributions of their temperature optima, with the area under each curve equal to the population density of the corresponding species. The green species is highlighted for purposes of illustration. Each species has migrants to adjacent patches (independent of local adaptedness), as well as immigrants from them (arrows from and to the green species; the distributions with dashed lines show the trait distributions of the green species’ immigrant individuals). The purple line is the intrinsic growth rate of a phenotype in the patch, as a function of its local temperature optimum (this optimum differs across patches, which is why the immigrants are slightly maladapted to the temperature of the focal patch.) Both local population densities and local adaptedness are changed by the constant interplay of temperature-dependent intrinsic growth, competition with other species in the same patch, immigration to or emigration from neighboring patches, and (in certain realizations of the model) pressure from consumer species.Full size imageSpecies in our setup may either be resources or consumers. Their local dynamics are governed by the following processes. First, within each patch, we allow for migration to and from adjacent patches (changing both local population densities and also local adaptedness, due to the mixing of immigrant individuals with local ones). Second, each species’ intrinsic rate of increase is temperature-dependent, influenced by how well their temperature optima match local temperatures (Fig. 2a). For consumers, metabolic loss and mortality always result in negative intrinsic growth, which must be compensated by sufficient consumption to maintain their populations. Third, there is a local competition between resource species, which can be thought of as exploitative competition for a set of shared substitutable lower-level resources26. Consumers, when present, compete only indirectly via their shared resource species. Fourth, each consumer has feeding links to five of the resource species (pending their presence in patches where the consumer is also present), which are randomly determined but always include the one resource which matches the consumer’s initial mean temperature optimum. Feeding rates follow a Holling type II functional response. Consumers experience growth from consumption, and resource species experience loss due to being consumed.Fig. 2: Temperature optima and climate curves.a Different growth rates at various temperatures. Colors show species with different mean temperature optima, with warmer colors corresponding to more warm-adapted species. The curves show the maximum growth rate achieved when a phenotype matches the local temperature, and how the growth rate decreases with an increased mismatch between a phenotype and local temperature, for each species. The dashed line shows zero growth: below this point, the given phenotype of a species mismatches the local temperature to the extent that it is too maladapted to be able to grow. Note the tradeoff between the width and height of the growth curves, with more warm-tolerant species having larger maximum growth at the cost of being viable for only a narrower range of temperatures62,63. b Temperature changes over time. After an initial establishment phase of 4000 years during which the pre-climate change community dynamics stabilize, temperatures start increasing at t = 0 for 300 years (vertical dotted line, indicating the end of climate change). Colors show temperature change at different locations along the spatial gradient, with warmer colors indicating lower latitudes. The magnitude and latitudinal dependence of the temperature change is based on region-specific predictions by 2100 CE, in combination with estimates giving an approximate increase by 2300 CE, for the IPCC intermediate emission scenario27.Full size imageFollowing the previous methodology, we derive our equations in the weak selection limit22 (see also the Discussion). We have multiple selection forces acting on the different components of our model. Species respond to local climate (frequency-independent directional selection, unless a species is at the local environmental optimum), to consumers and resources (frequency-dependent selection), and competitors (also frequency-dependent selection, possibly complicated by the temperature-dependence of the competition coefficients mediating frequency dependence). These different modes of selection do not depend on the parameterization of evolution and dispersal, which instead are used to adjust the relative importance of these processes.Communities are initiated with 50 species per trophic level, subdividing the latitudinal gradient into 50 distinct patches going from pole to equator (results are qualitatively unchanged by increasing either the number of species or the number of patches; SI, Section 5.9–5.10). We assume that climate is symmetric around the equator; thus, only the pole-to-equator region needs to be modeled explicitly (SI, Section 3.5). The temperature increase is based on predictions from the IPCC intermediate emission scenario27 and corresponds to predictions for the north pole to the equator. The modeled temperature increase is represented by annual averages and the increase is thus smooth. Species are initially equally spaced, and adapted to the centers of their ranges. We then integrate the model for 6500 years, with three main phases: (1) an establishment period from t = −4000 to t = 0 years, during which local temperatures are constant; (2) climate change, between t = 0 and t = 300 years, during which local temperatures increase in a latitude-specific way (Fig. 2b); and (3) the post-climate change period from t = 300 to t = 2500 years, where temperatures remain constant again at their elevated values.To explore the influence and importance of dispersal, evolution, and interspecific interactions, we considered the fully factorial combination of high and low average dispersal rates, high and low average available genetic variance (determining the speed and extent of species’ evolutionary responses), and four different ecological models. These were: (1) the baseline model with a single trophic level and constant, patch- and temperature-independent competition between species; (2) two trophic levels and constant competition; (3) single trophic level with temperature-dependent competition (where resource species compete more if they have similar temperature optima); and (4) two trophic levels as well as temperature-dependent competition. Trophic interactions can strongly influence diversity in a community, either by apparent competition28 or by acting as extra regulating agents boosting prey coexistence29. Temperature-dependent competition means that the strength of interaction between two phenotypes decreases with an increasing difference in their temperature optima. Importantly, while differences in temperature adaptation may influence competition, they do not influence trophic interactions.The combination of high and low genetic variance and dispersal rates, and four model setups, gives a total of 2 × 2 × 4 = 16 scenarios. For each of them, some parameters (competition coefficients, tradeoff parameters, genetic variances, dispersal rates, consumer attack rates, and handling times; SI, Section 6) were randomly drawn from pre-specified distributions. We, therefore, obtained 100 replicates for each of these 16 scenarios. While replicates differed in the precise identity of the species which survived or went extinct, they varied little in the overall patterns they produced.We use the results from these numerical experiments to explore patterns of (1) local species diversity (alpha diversity), (2) regional trends, including species range breadths and turnover (beta diversity), (3) global (gamma) diversity, and global changes in community composition induced by climate change. In addition, we also calculated the interspecific community-wide trait lag (the difference between the community’s density-weighted mean temperature optima and the current temperature) as a function of the community-wide weighted trait dispersion (centralized variance in species’ density-weighted mean temperature optima; see Methods). The response capacity is the ability of the biotic community to close this trait lag over time30 (SI, Section 4). Integrating trait lag through time31 gives an overall measure of different communities’ ability to cope with changing climate over this time period; furthermore, this measure is comparable across communities. The integrated trait lag summarizes, in a single functional metric, the performance and adaptability of a community over space and time. The reason it is related to performance is that species that on average live more often under temperatures closer to their optima (creating lower trait lags) will perform better than species whose temperature optima are far off from local conditions in space and/or time. Thus, a lower trait lag (higher response capacity) may also be related to other ecosystem functions, such as better carbon uptake which in turn has the potential to feedback to global temperatures32.Overview of resultsWe use our framework to explore the effect of species interactions on local, regional, and global biodiversity patterns, under various degrees of dispersal and available genetic variance. For simplicity, we focus on the dynamics of the resource species, which are present in all scenarios. Results for consumers, when present, are in the SI (Section 5.8). First, we display a snapshot of species’ movement across the landscape with time; before, during, and after climate change. Then we proceed with analyzing local patterns, followed by regional trends, and finally, global trends.Snapshots from the time series of species’ range distributions reveal useful information about species’ movement and coexistence (Fig. 3). Regardless of model setup and parameterization, there is a northward shift in species’ ranges: tropical species expand into temperate regions and temperate species into polar regions. This is accompanied by a visible decline in the number of species globally, with the northernmost species affected most. The models do differ in the predicted degree of range overlap: trophic interactions and temperature-dependent competition both lead to broadly overlapping ranges, enhancing local coexistence (the overlap in spatial distribution is particularly pronounced with high available genetic variance). Without these interactions, species ranges overlap to a substantially lower degree, diminishing local diversity. Below we investigate whether these patterns, observed for a single realization of the dynamics for each scenario, play out more generally as well.Fig. 3: Species’ range shift through time, along a latitudinal gradient ranging from polar to tropical climates (ordinate).Species distributions are shown by colored curves, with the height of each curve representing local density in a single replicate (abscissa; note the different scales in the panels), with the color indicating the species’ initial (i.e., at t = 0) temperature adaptation. The model was run with only 10 species, for better visibility. The color of each species indicates its temperature adaptation at the start of the climate change period, with warmer colors belonging to species with a higher temperature optimum associated with higher latitudes. Rows correspond to a specific combination of genetic variance and dispersal ability of species, columns show species densities at different times (t = 0 start of climate change, t = 300 end of climate change, t = 2500 end of simulations). Each panel corresponds to a different model setup; a the baseline model, b an added trophic level of consumers, c temperature-dependent competition coefficients, and d the combined influence of consumers and temperature-dependent competition.Full size imageLocal trendsTrophic interactions and temperature-dependent competition indeed result in elevated local species richness levels (Fig. 4). The fostering of local coexistence by trophic interactions and temperature-dependent competition is in line with general ecological expectations. Predation pressure can enhance diversity by providing additional mechanisms of density regulation and thus prey coexistence through predator partitioning28,29. In turn, temperature-dependent competition means species can reduce interspecific competition by evolving locally suboptimal mean temperature optima22, compared with the baseline model’s fixed competition coefficients. Hence with temperature-dependent competition, the advantages of being sufficiently different from other locally present species can outweigh the disadvantages of being somewhat maladapted to the local temperatures. If competition is not temperature-dependent, interspecific competition is at a fixed level independent of the temperature optima of each species. An important question is how local diversity is affected when the two processes act simultaneously. In fact, any synergy between their effects is very weak, and is even slightly negative when both the available genetic variance and dispersal abilities are high (Fig. 4, top row).Fig. 4: Local species richness of communities over time, from the start of climate change to the end of the simulation, averaged over replicates.Values are given in 100-year steps. At each point in time, the figure shows the mean number of species per patch over the landscape (points) and their standard deviation (shaded region, extending one standard deviation both up- and downwards from the mean). Panel rows show different parameterizations (all four combinations of high and low genetic variance and dispersal ability); columns represent various model setups (the baseline model; an added trophic level of consumers; temperature-dependent competition coefficients; and the combined influence of consumers and temperature-dependent competition). Dotted vertical lines indicate the time at which climate change ends.Full size imageRegional trendsWe see a strong tendency for poleward movement of species when looking at the altered distributions of species over the spatial landscape (Fig. 3). Indeed, looking at the effects of climate change on the fraction of patches occupied by species over the landscape reveals that initially cold-adapted species lose suitable habitat during climate change, and even afterwards (Fig. 5). For the northernmost species, this always eventuate to the point where all habitat is lost, resulting in their extinction. This pattern holds universally in every model setup and parameterization. Only initially warm-adapted species can expand their ranges, and even they only do so under highly restrictive conditions, requiring both good dispersal ability and available genetic variance as well as consumer pressure (Fig. 5, top row, second and third panel).Fig. 5: Range breadth of each species expressed as the percentage of the whole landscape they occupy (ordinate) at three different time stamps (colors).The mean (points) and plus/minus one standard deviation range (colored bands) are shown over replicates. Numbers along the abscissa represent species, with initially more warm-adapted species corresponding to higher values. The range breadth of each species is shown at three time stamps: at the start of climate change (t = 0, blue), the end of climate change (t = 300, green), and at the end of our simulations (t = 2500, yellow). Panel layout as in Fig. 4.Full size imageOne can also look at larger regional changes in species richness, dividing the landscape into three equal parts: the top third (polar region), the middle third (temperate region), and the bottom third (tropical region). Region-wise exploration of changes in species richness (Fig. 6) shows that the species richness of the polar region is highly volatile. It often experiences the greatest losses; however, with high dispersal ability and temperature-dependent competition, the regional richness can remain substantial and even increase compared to its starting level (Fig. 6, first and third rows, last two columns). Of course, change in regional species richness is a result of species dispersing to new patches and regions as well as of local extinctions. Since the initially most cold-adapted species lose their habitat and go extinct, altered regional species richness is connected to having altered community compositions along the spatial gradient. All regions experience turnover in species composition (SI, Section 5.1), but in general, the polar region experiences the largest turnover, where the final communities are at least 50% and sometimes more than 80% dissimilar to the community state right before the onset of climate change—a result in agreement with previous studies as well7,33.Fig. 6: Relative change in global species richness from the community state at the onset of climate change (ordinate) over time (abscissa), averaged over replicates and given in 100-year steps (points).Black points correspond to species richness over the whole landscape; the blue points to richness in the top third of all patches (the polar region), green points to the middle third (temperate region), and yellow points to the last third (tropical region). Panel layout as in Fig. 4; dotted horizontal lines highlight the point of no net change in global species richness.Full size imageGlobal trendsHence, the identity of the species undergoing global extinction is not random, but strongly biased towards initially cold-adapted species. On a global scale, these extinctions cause decreased richness, and the model predicts large global biodiversity losses for all scenarios (Fig. 6). These continue during the post-climate change period with stable temperatures, indicating a substantial extinction debt which has been previously demonstrated34. Temperature-dependent competition reduces the number of global losses compared to the baseline and trophic models.A further elucidating global pattern is revealed by analyzing the relationship between the time-integrated temperature trait lag and community-wide trait dispersion (Fig. 7). There is an overall negative correlation between the two, but more importantly, within each scenario (unique combination of model and parameterization) a negative relationship is evident. Furthermore, the slopes are very similar: the main difference between scenarios is in their mean trait lag and trait dispersion values (note that the panels do not share axis value ranges). The negative trend reveals the positive effect of more varied temperature tolerance strategies among the species on the community’s ability to respond to climate change. This is analogous to Fisher’s fundamental theorem35, stating that the speed of the evolution of fitness r is proportional to its variance: dr/dt ~ var(r). More concretely, this relationship is also predicted by trait-driver theory, a mathematical framework that focuses explicitly on linking spatiotemporal variation in environmental drivers to the resulting trait distributions30. Communities generated by different models reveal differences in the magnitude of this relationship: trait dispersion is much higher in models with temperature-dependent competition (essentially, niche differentiation with respect to temperature), resulting in lower trait lag. The temperature-dependent competition also separates communities based on their spatial dispersal ability, with faster dispersal corresponding to greater trait dispersion and thus lower trait lag. Interestingly, trophic interactions tend to erode the relationship between trait lag and trait dispersion slightly (R2 values are lower in communities with trophic interactions, both with and without temperature-dependent competition). We have additionally explored the relationship between species richness and trait dispersion, finding a positive relationship between the two (SI, Section 4.1).Fig. 7: The ability of communities in four different models (panels) to track local climatic conditions (ordinate), against observed variation in traits within those communities (abscissa).Larger values along the ordinate indicate that species’ temperature optima are lagging behind local temperatures, meaning a low ability of communities to track local climate conditions. Both quantities are averaged over the landscape and time from the beginning to the end of the climate change period, yielding a single number for every community (points). The greater the average local diversity of mean temperature optima in a community, the closer it is able to match the prevailing temperature conditions. Species’ dispersal ability and available genetic variance (colors) are clustered along this relationship.Full size image More
