Mutualistic interactions reshuffle the effects of climate change on plants across the tree of life.

Climatically induced local species extinctions may trigger coextinction cascades, thus driving many more species to extinction than originally predicted by species distribution models. Using seven pollination networks across Europe that include the phylogeny and life history traits of plants, we show a substantial variability across networks in climatically predicted plant extinction-and particularly the subsequent coextinction-rates, with much higher values in Mediterranean than Eurosiberian networks. While geographic location best predicts the probability of a plant species to be driven to extinction by climate change, subsequent coextinctions are best predicted by the local network of interactions. These coextinctions not only increase the total number of plant species being driven to extinction but also add a bias in the way the major taxonomic and functional groups are pruned.

INTRODUCTION

Climate change affects the physiology, abundance, and distribution of species (–), and recent evidence suggests that it may also affect species interactions (, ). A remaining challenge, however, is to be able to predict these effects when progressing from pairwise interactions all the way to species-rich networks (). Extinctions induced by climate change may trigger coextinction cascades—groups of species disappearing as a consequence of the extinction of species they depend on—thus driving many more species to extinction than originally predicted (). The magnitude of these coextinctions may depend on the web of mutual dependencies among species and how this web affects community persistence (, ). Species distribution models (), however, have traditionally treated species independently from each other. The gap between research on climate change and ecological networks is only now beginning to be reduced (, –). A first approach has added species interactions (oftentimes inferred through species co-occurrences) as an additional filter to improve species distribution forecasts (–). A second approach, in turn, has started to explore the effects of climate change on network structure (, ) and robustness (). For example, mutualistic networks, such as those between plants and their pollinators or seed dispersers, have been found to be more sensitive to climatically projected plant (rather than animal) extinctions (). It remains to be seen, however, to what degree the factors predicting direct extinctions are decoupled from those predicting subsequent coextinctions and whether the presence of coextinction results in different pruning of the taxonomic and functional groups. This information is important if we want to understand the consequences of these coextinctions for the functioning of these networks and their potential to adapt, which very much hinge on their functional and phylogenetic diversity. Here, we use climatically informed local plant extinction sequences in different communities and examine how subsequent coextinction cascades erode plant taxonomic, phylogenetic, and functional diversity. Our network approach is similar to that of previous work using species distribution models in that our model includes neither potential adaptation nor immigration ().

Our study is based on a dataset of seven European pollinator networks compiled from the web-of-life dataset (see Materials and Methods), for which phylogenetic relationships (figs. S1 and S2) and information about plants’ functional traits (database S1) were gathered from various sources. Specifically, the seven networks belong to Mediterranean () and Eurosiberian () biogeographic regions (Fig. 1), ranging from 52 to 797 species (plants and pollinators), and contain information about the following plant traits: clonal reproduction (yes/no), seed dispersal (wind/animal/others), flower shape limits pollen exchange by animals (yes/no), seven classes of flower color, and five types of Raunkiaer life-form (see Materials and Methods and database S1). These are functional traits central to plant life history. For example, the seed bank may allow a species to persist across unfavorable years, seed dispersal may allow colonizing distant habitats, and floral shape may condition the degree of dependence on pollinators.

Fig. 1

Climate warming across Europe and geographic locations of the pollination networks.

Color codes the predicted increase in maximum temperature of the warmest month between 2020 and 2080 (one of the five variables considered in the species distribution models used here). The warmer the color, the higher the increase in temperature between these time horizons. The points indicate the location of the seven networks, which are plotted in the insets with green and orange nodes representing plant and insect species, respectively. The different colors of the dots serve to identify each network in the following figures. The two darker colors in the South indicate the two Mediterranean networks, while the rest belong to the Eurosiberian region.

We collected geographic distribution data for the 244 distinct plant species of the seven pollinator networks (see Materials and Methods and database S1). Subsequently, we considered the projected climatic conditions for time horizons 2050 and 2080 [through the main text and figures, we are considering the A1b socioeconomic scenario (SRES) using the RCA30 regional climate model (RCM) driven by the ECHAM5 global circulation model (GCM) and ensembles of six different species distribution models (see Materials and Methods). However, we have also explored other scenarios and combinations (see Materials and Methods) with qualitatively similar results (figs. S3 and S4 and tables S1 and S2).] (see Materials and Methods, Fig. 1, and fig. S5). For each time horizon, we ran species distribution models to assign climatic suitability to each plant species at each one of the seven locations where it is present. From this information, we estimate an extinction probability based on the degree to which the predicted future local climate at the network location is similar to the plant’s climatic range, measured across its current Euro-Mediterranean distribution (see Materials and Methods and database S2). At this stage, several plant species may have gone locally extinct due to the direct effects of climate change. As a consequence, their interactions in the pollination network are no longer present. Since a species experiencing a reduction in the number of partners it depends on may have reduced fitness, we further imposed a probability of coextinction as follows. Each plant and animal species that has faced a reduction on its number of partners is assigned a coextinction probability proportional to the fraction of interactions lost [for three of the networks containing weighted information, we estimated a weighted version of the coextinction probability, which is proportional to the sum of strengths of the interactions lost. Overall, both measures of coextinction probability were significantly correlated (see Materials and Methods and fig. S6)]. We assume that those species left without any resource will go locally extinct. This procedure is repeated recursively so that coextinctions cascade through the network until no more species are affected, independently for each time horizon. This enables us to disentangle the direct impact of climate change on extinction rates versus the indirect effects mediated by the depletion of biotic interactions. We keep track of the identity of the plant species that go extinct, so we can later on assess to what degree the extant species represent a biased sample of the original communities (i.e., before the climatically induced extinctions).

RESULTS

Figure 2 shows the fraction of plant species lost through the two time horizons for each of the seven pollinator networks. Each panel compares the direct loss of species induced by climate versus the same amount plus the subsequent coextinction cascades for 2050 and 2080. The difference between both amounts (the slope of the lines in Fig. 2), therefore, shows how many more species are predicted to be lost, as one takes into account their mutual dependencies.

Fig. 2

Plant species extinctions and subsequent coextinction cascades for each pollination network across the two time horizons.

Each panel compares the direct loss of species induced by climate change versus the same amount plus the subsequent coextinction cascades for 2050 and 2080. As in Fig. 1, the different colors identify the specific networks, with the darker and lighter colors representing the Mediterranean and Eurosiberian networks, respectively. For visualization purposes, different points are slightly displaced across the x axis when they overlap. Figure shows the average and SD fraction of species lost across 1000 replicates.

As Fig. 2 demonstrates, there is quite a broad range of biogeographical variability in the rate of species extinctions. The two Mediterranean networks (dark blue and brown circles in Fig. 2) are the ones experiencing a higher rate of extinction. The difference is much higher, however, when considering the subsequent coextinctions. This suggests that the Mediterranean networks may be much more vulnerable to climate change than the Eurosiberian ones, albeit our sample size is too small to make a general claim. This would be expected primarily not only because climate change is assumed to be greater in Mediterranean zones () but also because Mediterranean networks contain plant species with narrower distribution ranges (t = 19.1, df = 277, P < 0.001, two sample t test) and therefore a higher probability of extinction as predicted by species distribution models. This higher extinction probability pushes the network closer to the threshold of collapse as identified by previous papers on the fragility of ecological networks (–).

Besides the apparent major differences between Mediterranean and Eurosiberian networks, there is also some variability in coextinction rates across the Eurosiberian networks for the time horizon of 2080. This observation leads us to the more general question of assessing what are the best correlates predicting the fate of species. To answer this question, we compared the explanatory power of several models that predict climatically induced species extinctions on one hand and extinctions resulting from the combined effect of climate and subsequent coextinctions on the other hand based on phylogenetic information [phylogenetic signal was quantified by optimizing Pagel’s λ, which varies between 0 and 1 and quantifies the extent to which resemblance between species reflects their evolutionary history as approximated by Brownian motion ().], geographic location (defined by latitude, longitude, and their interaction), and network identity (i.e., any other property of each network such as, for example, its topology). We then compared candidate models using Akaike’s information criterion (AIC; see Materials and Methods and tables S3 to S7).

Our statistical analyses reveal that while geographic location is the best predictor of the probability of climatically driven extinctions (table S3), network identity is a substantially better predictor of the ultimate fate of a species through the combined effect of climatically induced extinctions and subsequent coextinctions (ΔAIC > 4.27 in the two time horizons; table S4). This provides strong evidence that information about the network of mutual dependencies among species is key when assessing the response of communities—rather than that of isolated species—to climate change.

Once we have described the rate of species lost and the differential predictors of both climatically induced extinctions and those extinctions plus subsequent coextinctions, we next partition the vulnerability to each of these two processes across plant species. Quite a large fraction of plant species that have a low probability of becoming extinct through climate change show a high probability of being driven to coextinction (Fig. 3). This means that the effect of species interactions and their induced coextinctions is not only increasing the number of species disappearing but also inducing a bias in the identity of those species. This may result in a different pruning of the phylogenetic and functional trees than the one resulting from merely climatically induced extinctions (Fig. 4). Specifically, plant orders with a higher probability of disappearing by 2080 through the combined effects of extinction and coextinction include Malvales, Dipsacales, Santalales, Oxalidales, and Asparagales, represented in the networks by 16, 5, 3, 1, and 25 species, respectively (Fig. 5A). In line with this, although there is no phylogenetic signal in the overall pattern of plant coextinctions when looking simultaneously at all the networks (table S5), there is a phylogenetic signal in four of the networks when analyzed separately in 2050; this signal, however, remains significant in only one of the networks by 2080 (fig. S7), which would suggest that our ability to predict the vulnerability of species to climate change based on phylogenetic information decreases with the length of the temporal horizon.

Fig. 3

Coextinctions may target different plant species than climatically induced primary extinctions.

The figure shows the probability of coextinction for each plant species at each network versus its probability of being climatically driven to extinction at time horizon 2080. Species toward the right of the red x = y isocline show a higher vulnerability to climate, while those toward the left are more susceptible to the loss of biotic interactions. Details of the simulations are the same as in Fig. 2.

Fig. 4

The pruning of plants’ phylogenetic and functional trees through extinctions and subsequent coextinctions.

Phylogenetic tree (A) and functional similarity dendrogram (B) of all plant species of the seven pollination networks. The inner color bar circle represents the probability of direct, climatically induced extinction. The outer circle, in turn, represents the overall probability of being driven to extinction by either climate or the subsequent coextinction cascades. Time horizon is 2080. Details of the simulations are the same as in Figs. 1 to 3. The total number of plant species is 244, and the extinction and coextinction probabilities are the average of 1000 replicates.

Fig. 5

Plant vulnerability across orders and functional groups.

The figure shows the ranking of plant orders (A) and functional groups (B) to the overall risk of extinction and subsequent coextinctions for time horizon 2080. Figure represents the average and SD across 1000 replicates of climatically induced extinctions (gray) and those plus subsequent coextinctions (black). The different functional groups are as follows: group 1, open flowers, not clonal, wind dispersed; group 2, closed flowers, not clonal, wind dispersed; group 3, open flowers, not clonal, animal dispersed; group 4, open flowers, clonal, animal dispersed; group 5, closed flowers, not clonal, yellow, neither animal nor wind dispersed; group 6, closed flowers, not clonal, other colors, neither animal nor wind dispersed; group 7, open flowers, clonal, geophytes, neither animal nor wind dispersed; group 8, open flowers, not clonal, therophytes, neither animal nor wind dispersed; group 9, open flowers, not clonal, chamaerophytes, neither animal nor wind dispersed; group 10, open flowers, not clonal, geophytes, neither animal nor wind dispersed; group 11, closed flowers, clonal, hemicryptophytes, neither animal nor wind dispersed; group 12, closed flowers, clonal, hemicryptophytes, wind dispersed; group 13, open flowers, clonal, hemicryptophytes; group 14, open flowers, not clonal, hemicryptophytes; group 15, closed flowers, clonal, geophytes; group 16, open flowers, clonal, chamaephytes.

As for the case of taxonomic groups, the combined effects of extinctions and subsequent coextinctions can target different functional groups and therefore have contrasting effects on the functioning of local communities and their ability to cope with further environmental change (Fig. 5B). One example is functional group 10 constituted by plant species that have open flowers, are not clonal, have an underground storage organ (geophytes), and are neither animal nor wind dispersed. Another example is group 15 constituted by plant species showing open flowers, clonal reproduction, and perennating buds at the soil surface (hemicryptophytes). This group has an almost null probability to be climatically driven to extinction but a probability close to 0.3 to disappearing through subsequent coextinctions (Fig. 5B).

Because coextinctions target different phylogenetic and functional groups and given the fact that we observe a highly significant phylogenetic signal across all functional traits (P < 0.001 in all cases), one could conclude that the functional dendrogram (Fig. 4B) merely depicts phylogenetic relationships (Fig. 4A). In line with this, the effect of locality and network identity on coextinction probabilities based on functional similarity is qualitatively identical to analyses including phylogeny (tables S5 and S6). However, a cophenetic correlation (c = 0.045, P = 0.040) and a Mantel test (Spearman r = 0.052, P = 0.029) indicate a very weak albeit significant association between these descriptors. Therefore, our incorporation of plants’ functional traits provides additional information to that provided by the phylogeny ().

Although information about the pollinators is much too scarce to run species distribution models, we can perform a similar analysis on the insect coextinction cascades resulting from plant extinctions and coextinctions. After building a phylogeny for all the insects present in the seven pollinator networks (fig. S8), we repeated our model selection analyses for the best predictor of insect coextinction. In strong agreement with the patterns found in plants, our results show that the local network of interactions provides a better predictor of the joint risk of extinction and coextinction than biogeographic information at both time horizons (table S7). In addition, there is a significant phylogenetic signal in the combined extinctions and coextinctions for the time horizon 2050 (λ = 0.24; AIC is lower than the equivalent when λ is forced to 0; table S7). This signal, however, is absent in 2080, which again seems to suggest a decrease in our ability to predict species vulnerability based on phylogenetic information with the length of the temporal horizon.

DISCUSSION

Our approach certainly makes some simplifications that are worth considering. To begin with, one should be careful when estimating climatically induced extinction probabilities from the decrease in climatic suitability predicted by species distribution models. For example, species may have different climatic niches than the ones observed, which would distort our estimation of extinction risk. In addition, our climatic scenario predicts the extinction of species at a particular site but does not consider the possible arrival of new species at this site. Long-lived herbs may stay alive for even decades without recruitment, and populations may be more resilient than predicted because of evolutionary responses and trophic flexibility (, –). Similarly, our approach does not consider pollinator-mediated benefits such as those arising when one plant sharing a pollinator with another plant benefits from the extinction of the latter. Adding these levels of realism could reduce the indirect effects of climate change (, , ). Therefore, one has to look at our scheme as the worst-case scenario. On the other hand, our framework is missing climatically induced extinctions of insects, which would add to the simulated coextinction cascades, thus pushing even more plants toward extinction. Last, it is not clear to what extent our results would apply to other types of interactions such as in food webs. Given the similar collapse of food webs once species extinctions reach a critical point (, ), however, one would anticipate a similar role of coextinctions as the one reported here for pollination networks (when considering the response of consumers following the extinction of their resources).

We have shown that considering the network of interdependencies among species is relevant when predicting extinction risk in the face of climate change. Incorporating species interactions not only increases the pool of species most likely being driven to extinction but also changes the way extant species are selected from the evolutionary and functional trees with potential implications for the functioning and robustness of the resulting communities.

MATERIALS AND METHODS

Species distributions

We considered the following plant-pollinator networks from the repository www.web-of-life.es: M_PL_006, M_PL_007, M_PL_009, M_PL_017, M_PL_018 (corresponding to Eurosiberian regions), and M_PL_015 and M_PL_016 (corresponding to Mediterranean regions). Species distribution maps (presence/absence maps) were collected from the following sources: Flora dels Països Catalans (by digitizing the different volumes), Atlas Florae Europaeae, and Hulten’s Flora (–). For the species not represented in these atlases, occurrence records were collected from the online databases GBIF (www.gbif.org) and Anthos (www.anthos.es). All distribution maps were reported on a Common European Chorological Grid Reference System (CGRS) grid (50 km by 50 km; datumD) using the Europe Lambert Conformal Conic coordinate system, edited with image editing software (ImageMagick and GIMP), and processed with geographic information systems (ArcMap). Given that GBIF and Anthos only provide presence data, we generated pseudo-absence records in the unoccupied grid cells.

Climatic information

Species distribution models were calibrated using current climatic information from the WorldClim database () at a 10–arc min resolution. Current climate was represented by five different bioclimatic variables: temperature seasonality (within-year SD × 100; bio 4), maximum temperature of the warmest month (bio 5), minimum temperature of the coldest month (bio 6), precipitation of the wettest month (bio 13), and precipitation of the driest month (bio 14). These variables were previously shown to present little multicollinearity in the study area (). Species distributions were projected according to different future climatic scenarios represented by several RCM runs generated by the ENSEMBLES EU project, which downscaled the very coarse resolution GCMs’ climatic data obtained for the fourth assessment report of the Intergovernmental Panel on Climate Change (). Specifically, we used three different RCMs (HadRM3, RCA3, and RACMO2) (–) downscaled from three different GCMs (HadCM3, ECHAM5, and CCSM3). We also considered three different SRESs available for these models (A1b, A2, and B1) (). All combinations of RCM/GCM/SRES scenarios were interpolated to the same 10–arc min resolution.

Species distribution modeling

We calibrated the distribution models with presence and absence records provided by the atlases. The sampling of the species ranges provided by GBIF and Anthos may rather be incomplete, and our initial presence–pseudo-absence datasets may therefore include many false absences, which may, in turn, induce an underestimation of these species climatic niches (, ). To compensate for this potential underestimation, we discarded pseudo-absences that fell inside these species rectilinear surface envelope (, ) using the sre function of the biomod2 package in R and calibrated distribution models with the remaining ones. Given that the species distribution modeling accuracy may vary depending on the statistical models used and the quality of the input data (), we used a consensus method (, ) by considering six different probabilistic models implemented in the biomod2 package in R: generalized additive model, generalized boosting model, artificial neural network, multivariate adaptive regression splines, random forest, and Maxent. For each model run, we used 70% of the original data for the model calibration and kept 30% for model evaluation. This procedure was repeated 20 times to ensure that the model predictive accuracy was not affected by the random-splitting strategy. Moreover, we evaluated this accuracy according to two different metrics: the true skill statistics (TSS) () and the area under the receiver operating characteristic curve (ROC) (). For each species, we ran a total of 240 models using six different statistical models, two evaluation metrics, and 20 repetitions.

Last, we projected each species’ presence probabilities under the current and future climatic conditions of its respective pollination network location. To summarize the different ensembles of projections (240 projections per time step, per RCM, and per species), we discarded models exhibiting TSS and ROC scores below 0.6 and 0.8, respectively, and computed a mean of the occurrence probabilities of the remaining ones, weighted by their respective model predictive accuracy (i.e., TSS scores). Overall, the different models exhibited high evaluation scores (mean TSS score, 0.82 ± 0.09; mean ROC score, 0.95 ± 0.05), and atlas and GBIF species exhibited similar ones (mean TSS score, 0.81 ± 0.08 for atlas species and 0.87 ± 0.11 for GBIF species; mean ROC score, 0.95 ± 0.04 for atlas species and 0.95 ± 0.06 for GBIF species). We therefore obtained one probability of occurrence per species, per network, per time step, and per RCM.

Extinction simulations

We used the occurrence probability for each species provided by the species distribution models as described in the previous section. For each network, scenario, and time horizon (current, 2050, and 2080, respectively), the decrease in occurrence probability for species i was given by , where is the occurrence probability at present time and is the occurrence probability at time horizon h. For those species with an undetermined occurrence probability, the decrease in occurrence was assigned to the value of the maximum decrement of occurrence within its community for that particular time horizon and scenario.

We then followed by simulating species extinctions. Simulations were performed following a stochastic process, where the probability of a plant species being driven to extinction is proportional to its decrease in occurrence. This stochastic process is independent for each network, scenario, and time horizon. That is, calculations are both from the present time to 2050 and from the present time to 2080. We ran 1000 replicates. This first stage leads to a number of extinguished plant species and a number of directly coextinguished insect species. We proceeded by considering plant coextinctions. The probability of coextinction is proportional to the fraction of interactions lost through the simulation. For each network, scenario, and time horizon, the probability of coextinction for species i at iteration t is given by , where is the number of interactions of species i at time step t (after removing extinct and coextinct species) and is the number of interactions at previous time step. This reflects the fact that a given species can experience a reduction in its number of interactions (and so, have a nonzero probability of becoming coextinct) at several times as coextinction cascades travel through the network. This scheme ensures that the probability of coextinction for each species is the same than if all the interactions were lost at once. The process was executed recursively at each time step until there were no further losses of interactions. For the networks containing weighted information describing the frequency of visits by each pollinator (three of seven), we calculated a weighted version of our probability of coextinction (estimated as the relative fraction of interaction strength contained in the extinct interactions) and compared it to the previous one to assess whether similar trends would be expected for more realistic coextinction scenarios.

Phylogenetic tree

We built the plant phylogeny encompassing 244 species using phylomatic, which combines taxonomic information with a backbone tree to obtain the phylogeny of a given set of plant species (figs. S1 and S2) (). We used their conservative tree with relations up to the family level and pseudobranch lengths in million years and resolved relationships within Asteraceae encompassing 56 species following (). Species in multiple networks were included as soft polytomies in analyses because they exhibit multiple estimates of extinction and coextinction probabilities (276 tip data). We then estimated the phylogenetic signal of categorical functional traits with a function designed ad hoc by one of us (E.L. Rezende, phylo.signal.disc). This function randomizes the tip data and compares the minimum number of character state changes inferred with parsimony with a null model (). An equivalent phylogeny comprising 1122 tip data was built for the pollinators combining phylogenetic and taxonomic information. Briefly, a backbone phylogeny of insect orders was built on the basis of the relationships described in two large-scale analyses (, ). Then, relationships between families within each order were resolved on the basis of different sources: Coleoptera (), Lepidotpera (, ), Hemiptera (), Hymenoptera (), and Diptera ().

Functional diversity tree

The functional tree of plants was constructed by applying the Sorensen distance () (S7 coefficient of Gower and Legendre, function dist.ktab in library ade4) to the plant functional traits. This creates the distance matrix that is used to generate the tree. For that purpose, we used hierarchical clustering to build the most reliable dendrogram of all species in functional trait space using the complete linkage method (function hclust) (). Species with unknown information for any of the traits were discarded, which reduced the tree to 179 species of the 244 original ones. As with the phylogeny, species in multiple networks were included more than once and emerged as polytomies in the functional tree (207 tip data).

Statistical analysis

The phylogenetic and functional structure of extinction and coextinction patterns were analyzed with generalized linear models. We included the phylogenetic or functional correlation matrix with off-diagonals multiplied by a parameter λ, which fitted to the residual distribution of the model and quantified the strength of association between related or functionally similar species (). We compared three models. The first model included only the intercept, the second model also added spatial information embedded in the latitude, longitude, and the interaction, and the third model also included network identity as a categorical factor. We then used Akaike weights to estimate the relative evidence favoring each model.