We did two sets of 1000 simulations each where we first populated an ecologically plausible ‘virtual Earth’ with open, trophically structured communities that we then subjected to rapid and catastrophic global environmental-change trajectories, consisting of either a planetary heating (for example, following volcanic eruptions26), or cooling (as in a ‘nuclear winter’15,27). For each trajectory, we replicated the environmental-change phase using two alternative scenarios. In the first, we assumed that species go extinct only due to environmental tolerances being exceeded, i.e., when environmental conditions become unsuitable for their survival. In the second scenario, we considered the additional contribution of co-extinctions to diversity loss, which we modelled as extinction cascades through simulated food webs, triggered by the primary extinctions driven by environmental tolerances alone (as in the first scenario). Henceforth, we refer to those two scenarios as the ‘environmental tolerance’ and ‘co-extinction’ scenarios, respectively. By keeping track of global diversity at each step of the temperature trajectories in both scenarios, we were able to obtain two measures of robustness of planetary life (with respect to each scenario’s assumptions) as the area under the curve of diversity versus time (temperature). Comparing those measures allowed us to disentangle the direct effect of environmental tolerances from that of co-extinctions, and hence obtain, for the first time, a quantitative assessment of the relative contribution of co-extinctions to the global biodiversity crisis. We provide details on the different model components in the following sections. All the code implementing our model and permitting full replication of the study (with line-by-line detailed comments) is is freely available at https://github.com/giovannistrona/co-extinctions.

Overview

We ran two sets of 1000 simulations for planetary heating and cooling trajectories, respectively. For each simulation, we generated a set of 10000 species having ecologically plausible tolerance limits, trophic level, and specificity in the use of resources (see Data calibration). We also generated a subset of 100 species as extremophiles, with large tolerance limits (see Generating virtual species) — we refer to those species as ‘tardigrades’, and we included them in our simulations to test explicitly how much assuming species’ tolerances to extreme conditions as a measure of their extinction risk in a changing environment14,25 departs from ecological realism. We attributed to each species a phenotype consisting of a set of functional traits determining a consumer’s accessibility to a given resource (see Functional traits). Those contributed much to the ecological realism of our simulations, playing a fundamental role in the construction of food webs (see Building local networks), and in the dispersal and colonization processes (see Species dispersal/colonization).

In each simulation, we generated a random, spatially explicit set of localities having realistic climatic features consistent with their relative geographic position on the virtual Earth. We then extracted random samples of species from the 10000 set and ‘dropped’ them into each locality, retaining in each only those species having compatibility with the local (locality) climate conditions (see Generating virtual localities and populating them with communities). The size of the samples of candidate species was the same for each locality, but varied among simulations, being randomly sampled with uniform probability between 100 and 1000. This provided variation in local/global diversity, which enabled us to explore the effect of varying species richness on community robustness to environmental change across simulations. We then arranged each initial pool of species in a structured food web (Section 7), and finally filtered them by excluding all species without trophic links (Section 6). To increase the ecological and biogeographical realism of the virtual Earth, we simulated a large number of dispersal events/colonization attempts, randomly sampling for each simulation with uniform probability between 103 and 105 events.

After the preliminary dispersal/colonization phase, we subjected the virtual Earth to ‘rapid’ global environmental change. We did not focus on a precise time scale for change, we instead assumed this to be much smaller than an evolutionary time scale by excluding speciation processes from our simulations. We focused instead on the magnitude of the temperature change as our temporal frame of reference. We modelled this in two opposing ways as a progressive, linear monotonic (1) increase, and (2) decrease in temperature. In particular, we shifted the upper and lower temperature limits of each locality by a random value of increase or decrease of temperature (in °C) sampled from a normal distribution with mean 0.01 and standard deviation 0.0025 for 5000 steps (with a final average global increase/decrease of ~ 50 °C that would be enough to annihilate all life directly, with the possible exception of some extremophile species). We also took into account a potential latitudinal effect on the speed of environmental change, as currently observed with current global warming where a nearly linear increase in temperature change has been observed in the last decades from 60° latitude upward, leading to a twofold increase in temperatures at the North Pole relative to the Equator31. For this, we increased the magnitude of temperature change at latitudes (y) > 60° north/south as (y-60)/30.

At each step of the environmental-change trajectory, we removed from each locality all species with temperature-tolerance limits no longer compatible with the changed conditions (see Measuring environmental compatibility). This single mechanism defined species loss in the environmental-tolerance scenario. In the co-extinction scenario, in addition to the primary extinctions caused by climate change at each step, we also accounted for the loss of consumers driven to extinction by the depletion of their resources. In so doing, we explored various assumptions regarding the minimum amount of resources ensuring the survival of a consumer, and the ability of the food web to rearrange interactions following species loss (see Modelling co-extinctions).

Data calibration

We derived a proxy for the distribution of trophic levels in natural systems by merging into a single list all the trophic levels of individual members of a large set of food webs32. We computed trophic level for each species as the shortest path length (i.e., steps through network links) from the target species to a basal resource33. We also associated a ‘specificity’ to each individual’s trophic level, computed as the fraction of resources consumed by the target species over the total number of items in the food web to which it belonged. We obtained information on species tolerances to temperature limits from various datasets (we focused on terrestrial species — see Results and Discussion for justification). In particular, we used data for 458 species of endotherms34, and 239 species of terrestrial ectotherms35. For plants, we could only find data for cold tolerance, but not for upper temperature limits; thus, we generated an original dataset of plant temperature tolerances based on species distribution. For this, we extracted plant occurrence data from the Global Biodiversity Information Facility36, limiting our search to records from “observation” or “literature occurrence”. This provided > 7.2 million occurrence records of > 25000 plant species. We then restricted our search to species having > 100 occurrences (on land), which reduced the set to 4445 species. We then combined these occurrences with Worldclim data37 to identify tolerance limits. To limit the number of outlier errors, we took the lower 95% bootstrapped confidence bounds of the minimum temperatures of the coldest month across all occurrence points, and the upper 95% confidence bounds of the maximum temperatures of the hottest month across all occurrence points, as the lower and upper tolerance limits, respectively. The purpose of the total list is not to provide species-specific tolerances; rather it provides a proxy for the ‘natural’ distribution of temperature tolerances in plants at the global scale.

Generating virtual species

In each simulation, we generated a species pool (100000 different species) to populate localities. For this, we:

(1). extracted an element from the list of trophic levels/specificities; (2). if the extracted trophic level = 0 (i.e., indicating a basal resource not consuming any other species in the food web), combined that value with the tolerance limits of a randomly extracted item from the Global Biodiversity Information Facility/Worldclim-derived plant dataset; (3). if the extracted trophic level > 0, combined it with the tolerance limits of a randomly extracted item from the endotherm list with a probability = 0.001, and with an item extracted at random from the ectotherm list with a probability = 0.999 (taking into account the predicted difference of around three orders of magnitude between endo- and ectotherm species diversity at the global scale38).

For 0.1% of the species (i.e., 100), we used different, ‘tardigrade-like’ tolerance limits: a random integer between −50 and −80, and a random integer between 50 and 100 for the planetary cooling and heating) trajectories, respectively22,23. We attributed to these virtual tardigrades a random trophic level between 1 and 2, consistent with tardigrades’ role as grazers and micro-predators in natural systems24, and a corresponding specificity value (see previous section). To each non-basal species (i.e., those with trophic level > 0), we assigned a real number sampled with uniform probability from 0 to 1 to indicate the relative position of the species within its trophic level; for example, this permits an apex predator to consume another apex predator. Finally, we provided each species with a set of functional traits, defining its ‘phenotype’, that we used to assess the accessibility of a given consumer to a candidate resource, as well as functional similarity between species (see next section).

Functional traits

The attribution of a functional phenotype to species in the simulations play an important role both to arrange species into food webs, as well as in species’ colonization processes. For each simulation, we implemented the functional traits and the related ecological mechanisms according to the following procedure:

(1). We first defined an arbitrary set of 26 different functional traits, with each trait being identified by a letter from A to Z. (2). We then populated an adjacency matrix with rows and columns corresponding to functional traits with real numbers sampled uniformly from −1 to 1. The value in each cell x ij indicates the degree to which the ith trait enables a consumer to use a resource having the jth trait. Positive values indicate that having the trait i increases the consumer’s accessibility to the resource with trait j, while negative values indicate that having the trait j protects (to some extent) the candidate resource from being consumed by a species with trait i. (3). We attributed a random set of letters of a size sampled with uniform probability from 1 to 9 to each species in the global diversity pool, as a proxy for phenotype. (4). We then quantified the ability of a consumer to use a potential resource (ca ij ) by summing all the x ij entries in the adjacency matrix for each ith trait of the consumer and each jth trait of the candidate resource. We rescaled this value from 0 to 1 by applying the formula:

$$\bar{c{a}_{ij}}=1-(c{a}_{max}-c{a}_{ij})/(c{a}_{max}-c{a}_{min})$$

with ca min and ca max being estimates of the minimum and maximum possible accessibility of a consumer to a resource. Since those values depend on the trait adjacency matrix, they are estimated as the minimum and maximum observed ca ij value over a large number (1 million) of pairwise comparisons between randomly generated phenotypes.

Measuring environmental compatibility

To assign species to a given locality (while populating the virtual Earth at the beginning of each simulation), and to determine the likelihood of success of colonization attempts as well as that of local species’ extinction events following environmental change, we focused on the compatibility between a species’ lower and upper thermal tolerance limits (sp_t, sp_T), and minimum and maximum temperatures of the target locality (loc_t, loc_T). In particular, we computed the minimum distance between the species’ tolerance limits and the local condition as:

$$Tol\_d=min[(loc\_t-sp\_t),(sp\_T-loc\_T)]$$

We then set the extinction probability to 0 in cases where the Tol_d value was ≥ 5, and to 1 in cases where Tol_d was < 0. Otherwise, we set the extinction probability to 1/(1 + Tol_d). Then, we evaluated the probability of a species going extinct at a given moment in the simulation by integrating over the curve of extinction probability versus time. We obtained this by computing Tol_d at 100 equally spaced time steps between time 0 and the target time, according to the corresponding temperature observed in the locality at each step. Time 0 corresponds to the moment a species established itself in the target locality, i.e., either the beginning of the environmental change simulation for native species, or the time of successful colonization for alien species. To derive the environmental conditions at a given moment (needed to compute the extinction probability at that moment, and hence compute the area under the curve), we invoked a linear relationship between time and temperature change as a good approximation, with a 1 °C of change corresponding to 1 time unit.

Generating virtual localities and populating them with communities

We collated minimum and maximum annual temperature data from Worldclim37, and interpolated them at a resolution of 1° × 1° latitude/longitude. Then we used the list of virtual species to ‘populate’ a set of localities extracted at random from all 1° × 1° Worldclim cells. For each target locality, we sampled at random n species from the species pool, with n for each virtual Earth generated being a random integer sampled with uniform probability in [100, 1000]. This gave different outcomes of average local species richness, and hence explores the sensitivity of our results to varying species diversity (see Figs S1, S2). Among the n candidates, we assigned species to the target locality on the basis of their niche compatibility with local climatic conditions (see previous section).

If fewer than 5 species among the candidates were attributed to a target locality, or if basal species contributed to < 20% of local diversity, we discarded the locality and replaced it with another random locality. Otherwise, we attempted to arrange species in the pool in a structured food web (see next section for details). If this was not possible, we discarded the locality, and replaced it with another candidate one. Otherwise, after we generated the food web, we excluded all the species not having trophic links from the local species pool. These theoretical communities allowed us to test the relevance of co-extinction processes in biodiversity loss.

We reiterated the procedure until we populated a random target number of localities sampled (for each simulation) with a uniform probability between 100 and 500. As in the case of global diversity, this allowed us to test the sensitivity of our results not only with respect to sample size, but also to various spatially explicit processes we included in our model, which are affected in turn by pairwise distance between localities, and hence by the density of localities in our simulated planet.

Building local networks

We used trophic levels (TL) and specificities (see Data calibration) to build a food web for each locality. Here, we assumed all i species in the target locality were potential consumers of all j species in the locality (with i ≠ j, and the order of the jth species randomized for every ith species). We then assigned species j as a resource to consumer i if TL j < TL i < (TL j + 1) and with a probability given by trait compatibility between the target consumer and the candidate resource \((\overline{c{a}_{ij}})\). For each target consumer, we stopped the comparisons with potential resources when the number of assigned links was equal to the product of the consumer’s specificity and the total number of species in the locality (in the following, we refer to this as the ‘expected number of resources’).

In cases where the expected number of resources was not associated to a target consumer after this first step, we attempted to connect the consumer to each one of the species within the same trophic level having a smaller intra-level trophic score (see Generating virtual species), again with a probability dictated by resource-consumer trait compatibility \((\overline{c{a}_{ij}})\). We presented candidate resources to a target consumer in random order, with the process terminating in cases where the expected number of resources was reached. If even this second step did not result in the expected number of resources to a target consumer, and the consumer belonged to a trophic level > 1, we attempted to link a target consumer to each one of the resources of two trophic levels lower. As in the previous steps, we established links with a probability given by resource-consumer trait compatibility, presented potential resources to target consumer in random order, and terminated the process when the maximum number of expected resources was reached. This three-step procedure gives more realism to the food webs by increasing their connectivity and clustering, and by reducing their average path length.

Once all species in the locality were evaluated as potential consumers, we identified the set of basal species as those with a trophic level = 0. We then filtered the network by excluding from each target locality all species for which a path to basal species did not exist. This prevented the generation of unrealistic food webs where, for example, large carnivores have access to herbivore prey, but the latter have no access to food. Finally, we assigned weights to the food-web links on the basis of functional-trait compatibility between the target resource-consumer pair \((\overline{c{a}_{ij}})\).

Species dispersal/colonization

To make the virtual communities connected within a biogeographical context, we implemented a mechanism of species dispersal/colonization from one locality to another (with slight differences between the environmental-tolerance and co-extinction scenarios — see below). We assumed that the likelihood of successful colonization is a function of the distance between two localities, of the structure of the target community (with more mature/structured communities being less open to colonization than degraded communities), and of the ability of the colonizer to establish itself in the new locality. The latter depends on several factors, particularly the degree of tolerance of the colonizer to local climatic conditions, the colonizer’s trophic level, and the colonizer’s functional traits.

We simulated a single dispersal/colonization attempt in the model as follows:

(1). A random pair of localities i and j (with i ≠ j) is sampled with uniform probability; (2). A random species is then sampled with uniform probability from the source locality i; this disperses successfully to locality j with probability d ij −1, with d ij representing the shortest Euclidean distance (d) between two localities; (3). In cases of successful dispersal, colonization is attempted. Colonization fails regardless of the colonizer’s identity with a probability given by 1 - p i (the target locality’s invasion susceptibility), which is equal to the connectance of the local food web in the co-extinction scenario (i.e., the ratio between the number of links in a food web to its squared number of nodes, with the latter being the theoretical maximum number of possible pairwise interactions)20, and to the ratio of initial diversity (i.e., the number of species present in the target locality at the beginning of the simulation) to the diversity at the time of the colonization attempt in the tolerance scenario. (4). In the environmental-tolerance scenario, the colonizer outcompetes species of the same trophic level that are less tolerant to local environmental conditions (see Section 5), with a probability given by their trait similarity. We calculated this probability as 2.0 × M/T, where T is the total number of traits in both species’ phenotypes, and M is the number of matches (this formula yields 1 if the sequences are identical, and 0 if they have nothing in common). Thus, trait similarity identifies the potential invader’s competitors, and then colonization succeeds only if the invader is better adapted to the novel, rapidly changing conditions than local species. We acknowledge that ‘adaptation’ here is limited to thermal tolerance, and that there could be several other ecological aspects making a native species ‘better’ adapted to its native environment than an alien one, but in our simplified world, thermal tolerance was in fact the main element determining species distribution, making our choice ecologically reasonable. In cases where the potential colonizer has no competitors in the target locality, it succeeds in colonization with a probability given by its compatibility with local environmental conditions. That is, invaders with no local competitors succeed in colonization with the only constraint of tolerance to local temperature range. This completes the explanation for the colonization mechanism in the environmental-tolerance scenario; all the following steps in this section refer to the co-extinction scenario only. (5). In the co-extinction scenario, if colonization has not failed due to local invasion susceptibility alone, the potential colonizer attempts to ‘enter’ the food web. If the colonizer is not a basal resource (i.e., trophic level > 1), it will need to find some resources to consume in order to establish itself in the target locality. For this, it will attempt to replace each local consumer having the same trophic level in each one of its resource-consumer interactions. For this, for each link in the network involving a consumer of the same trophic level as the potential colonizer, the local consumer will be replaced by the colonizer if the colonizer has an identical or better compatibility with local environmental conditions than the local species, and with a probability given by the rescaled trait compatibility \((\overline{c{a}_{ij}})\) between the colonizer and the target resource. In cases of successful colonization, the old trophic link is removed from the food web, while the new one between the colonizer and the local resource is added. The interaction weight is inherited from the original link between the local resource and the out-competed local consumer. This choice is motivated by the fact that the weight of a trophic link in our modelling framework represents resource availability. The availability of local resources is (at least at the beginning of the invasion) not affected by the switch in consumers between the previous species and the invader. Thus, an alien species replacing ≥ 1 local species in some or all of their trophic interactions will have access to the same amount of resources previously used by the native, now-outcompeted consumer/s, and hence be assigned trophic links having weight equal to that of the original interactions. (6). In cases where a non-basal colonizer has successfully established, the model accounts for the possibility that the colonizer becomes an additional resource for other species. Thus, for each resource consumer link in the food web, if the colonizer has the same trophic level of target resource, is accessible to consumers (i.e., with probability \(\overline{c{a}_{ij}}\)), and is well-adapted to local climate (i.e., with probability equal to climatic suitability), then the colonizer might be used as an additional resource. A new link is therefore added to the food web, connecting the colonizer to the target local consumer, with an interaction weight equal to a random fraction of the original link between local resource and consumer. Such a link is then maintained in the network. (7). If the potential colonizer is a basal resource, then it will have to outcompete, to some degree, another local basal resource, under the assumption that the system is at equilibrium carrying capacity. In the real world, such a process could lead to loss of local diversity, especially in cases where the local, outcompeted resource supports several species unable to survive by establishing novel interactions with the colonizer. Despite being ecologically interesting, analyzing this process further is beyond the scope of our study. To avoid situations where our results are biased by diversity loss triggered by biological invasions, we conservatively allowed basal colonizers to replace local basal species as resources in resource-consumer interactions, but we did not consider normal outcompeting events. As in the colonization processes we described above, the colonizer must have trait compatibility with the target local consumer (which also has ecological meaning, since it might indicate that the colonizer’s phenotype is in fact functionally similar to that of the local resource it is outcompeting).

According to these mechanisms, a successful colonization event can either: (a) lead to an increase in local diversity if the colonizer replaces ≥ 1 local species in only some of their resource-consumer interactions without fully outcompeting any of them (or, in the environmental-tolerance scenario, when a colonizer establishes itself by exploiting an empty niche), (b) leave local diversity unchanged even if the colonizer completely outcompetes a single local species, or (c) reduce diversity if the colonizer entirely outcompetes > 1 local species. In the environmental-tolerance scenario, diversity can also increase when a colonizer establishes itself by exploiting an empty niche. In the co-extinction scenario, when the colonizer is a non-basal resource that becomes an additional resource to local species, the potential increase in diversity is also paired to an increase in overall network interaction strength. However, this does not necessarily produce an increase in network connectance, since this decreases more rapidly with the addition of nodes than it increases with the addition of links.

Community rescue via exogenous recruitment

In addition to colonization events, we also implemented an explicit mechanism of recruitment from closer localities as a means to replenish depleted resource stocks. At each step of the simulation, every target locality i receives exogenous recruits from any other source locality (j) with a probability = 0.01. When recruitment happens, the interaction strength of each link of the food web in the ith locality is increased (up to the original interaction strength observed at the beginning of the simulation) by d ij −1 times its original value (with d ij representing the shortest linear distance between two localities), if the resource species in the target link is also present in the jth (i.e., the source) locality. As we discussed above, the interaction weight in our model is a proxy for resource availability, and hence, population size of the focal resource. The magnitude of recruitment is proportional to distance in that we expect more recruitment from closer localities than from more distant ones. The inverse relationship between distance and increase in interaction weight (ideally representing an increase in the target resource’s population size) offers a simple way to account for this aspect in our simulations.

Species adaptability to change

Although in our simulations we assumed environmental change happens on a temporal scale much smaller than evolutionary processes, we still provided species with the ability to develop some thermal adaptation — an ability observed in nature in both invertebrates39 and vertebrates40. We assumed that in each simulation step, each species shifted its thermal limits (by a random value sampled from a normal distribution with mean 0.75 and standard deviation 0.25) with a particular probability. We sampled this probability at random (with uniform probability) between 0.0 and 0.0001, and set it at the beginning of the simulation to remain stable throughout each simulation. Using a variable probability across individual simulations allowed us to test the sensitivity of our results to overall species adaptability.

Modelling co-extinctions

We modelled co-extinctions following Säterberg et al.12. We initiated simulations with the removal of species obliterated by mismatches of environmental tolerance as described above. This possibly triggered co-extinctions in consumers that experienced an initial loss of resources exceeding a certain fraction. We quantified the ‘amount’ of resources as the cumulative weights of the links from a target consumer to all of its resources, and we explored various thresholds for co-extinctions (in terms of fraction of lost interactions) by sampling in each simulation a real number at random (with uniform probability) between 0 and 1. We permitted the ‘rewiring’ of lost interactions; that is, we allowed resources made available by the extinction of some of their original consumers to be redistributed among extant consumers14. We constrained the redistribution among consumers already using the target resource, proportionally to the current use of that resource. We weighted the amount of redistributed resources by a ‘reallocation ratio’. As in the case of the co-extinction threshold, we explored different reallocation ratios by sampling in each simulation a real number at random (with uniform probability) between 0 and 1. We reiterated the three steps of the co-extinction process — (1) removing lost species, (2) evaluating lost resources, and (3) including the interaction rewiring — until there were no new (i.e., cascaded) co-extinctions, or when all species in the network had gone extinct.

Exploring food-web robustness to environmental-change trajectories

To explore the robustness of local food webs towards either climate heating or cooling, we ‘disassembled’ 1000 food webs sampled at random between all food webs in all virtual Earths (after the pre-dispersal/colonization phase to use the same food webs we subjected to environmental change). Disassembly consisted of removing species progressively from the least to the most tolerant to warm or cold temperatures alternatively. Following each species removal, we modelled co-extinctions according to the same procedure we used in our main simulations (see previous section). Finally, to assess robustness of the target food web, we built curves tracking the overall decline of local species diversity (accounting for both primary extinctions triggered directly by climate change, and co-extinctions) following progressive temperature change. We also identified approximate upper and lower boundaries for robustness by replicating the same procedure under different criteria of species removal16. To approximate the upper boundary of food-web robustness, we progressively removed species from those supporting the lowest number of consumers to those supporting the most consumers. To approximate the lower robustness boundary, we removed species in reverse order, starting from those supporting the most consumers, and then moving progressively towards species supporting the fewest or no consumers. We also generated reference curves for each food web by removing species in random order.