Plant and animal occurrence data

We compiled occurrence data for plant and animal species recorded in eight quantitative pollination and five quantitative seed-dispersal networks from central Europe (Supplementary Data 1). Range maps of plant distributions were compiled from published distribution maps40,41,42, occurrence data from the Global Biodiversity Information Facility (GBIF), national and regional floristic databases and further maps from the floristic literature (see bibliographic details in Index holmiensis43,44,45,46). Contiguous large areas of plant occurrence were generalized as range polygons; spatially isolated occurrences were digitized as single-point locations. Distributional data on most wild bee taxa were extracted from GBIF and from a database hosted at the University of Mons (Belgium) (Atlas Hymenoptera47,48,49,50,51,52,53). Most of the Andrena bee data were originally derived from maps associated with the Warncke collection (Biocentre of the Upper Austrian Museum, Linz, Austria54). Data on distribution of Colletes bees were provided by Michael Kuhlmann. Distributional data for European butterflies were provided by the database LepiDiv55, an updated version of the database used for the Distribution Atlas of Butterflies in Europe56. Data for hoverflies were provided by a database hosted at the University of Novi Sad57. Range maps of bird distributions were compiled from a database of global avian distribution maps58. All distribution data were gridded to match a European CGRS grid (3024 equal-area cells with a resolution of 50 × 50 km). We omitted areas from eastern Europe with a low sampling intensity for insect pollinators (see Supplementary Fig. 3 for the spatial extent of the used grid). Data aggregation at the 50 × 50 km grid further mitigates effects of low sampling intensity, that is, an overestimation of species’ absences for some insect taxa.

Climatic niche estimation and species distribution models

We computed the current range size of a species as the number of grid cells with observed presences within the European CGRS grid; in the case of avian migrants, only presences within the breeding range were considered. We omitted species from the analyses for which no or deficient occurrence data were digitally available; the number of occupied grid cells per species ranged from 19 to 2,917 (median=815 occurrences across the 50 × 50 km grid, Supplementary Data 2). Species’ realized climatic niches were quantified as a function of their occurrences and four variables of current climate59 (annual mean temperature, temperature annual range, annual precipitation and precipitation seasonality, sampled for the used grid). We used two alternative methods for quantifying the current climatic niche breadth of a species (see Supplementary Fig. 4 for a data overview). First, we used the hypervolume method60, which calculates the realized climatic niche breadth by using a multidimensional kernel density estimation procedure to estimate an n-dimensional hypervolume from a set of species’ occurrences and the respective climatic variables. To calculate the climatic hypervolume, we z-standardized climatic variables and used a Silverman bandwidth estimator and a 0% quantile threshold60. Second, we used the outlying mean index (OMI) that quantifies the distance between realized and background climate conditions61. In contrast to other multivariate methods for niche quantification (such as canonical correspondence or redundancy analysis), it makes no assumption about the shape of the response curve to the environment (that is, climate) and is not influenced by species richness62. Along each ordination axis, the OMI calculates niche breadth as variances based on the climatic conditions at the localities of species’ occurrences. Based on the first and second OMI ordination axes (cumulative inertia of both axes: 89%), we defined species’ realized niche breadth as the geometric mean of variances along these two axes. Both methods resulted in similar estimates of realized climatic niche breadth (n=709 species, r=0.92) and were positively associated with the current range size of a species (hypervolume, r=0.54; OMI niche breadth, r=0.67). We also tested whether the spatial extent of the analysis may have affected estimates of climatic niche breadth. We found that the OMI climatic niche breadth, derived from occurrences across the entire Palearctic for plants and birds, was closely correlated to that at the European scale (n=346 species, r=0.83). Occurrence data beyond Europe were not available for insect pollinators and all taxa were therefore analysed at the European scale.

We quantified probabilities of occurrence from species’ recorded presences and absences with species distribution models based on boosted regression trees63, using a cross-validation approach to estimate the optimal number of trees (number of initial trees=10, tree complexity=2, learning rate=0.01). To evaluate model performance, we calculated the area under the receiver-operator curve (AUC) and the True Skill Statistic (TSS)64 where the sum of the model’s sensitivity and specificity was maximal (TSS max ). AUC and TSS max were calculated as the arithmetic mean of 10 random splits of data into 75% used for model calibration and 25% for model testing. Arithmetic means (±1 s.d.) of AUC/TSS max values were: 0.85±0.046/0.58±0.093 (bees), 0.84±0.057/0.57±0.120 (butterflies), 0.79±0.057/0.49±0.107 (hoverflies), 0.94±0.031/0.77±0.074 (birds), and 0.95±0.026/0.80±0.073 (plants), indicating good to very good (pollinators) or excellent (birds, plants) model performance under current conditions (Supplementary Data 2). We used the full set of current occurrences of each species for calculating model projections under current and future conditions. Future climate projections were obtained from two general circulation models from the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC) (CCSM4, MIROC5 (refs 59, 65)) using two scenarios of RCPs assuming an average increase of 2.85±0.62 °C (RCP 6.0) or 4.02±0.80 °C (RCP 8.5) in mean annual temperature for the geographic area covered (see Supplementary Fig. 3 for the projected changes in all climatic variables). We quantified the potential vulnerability of a species to projected climate change by changes in climatic suitability for each grid cell covered by the species’ current European range, defined by the difference between the probabilities of occurrence under current (years 1950–2000) versus projected future conditions for 2070 (averaged for 2061–2080). For each species, changes in climatic suitability were summarized using the median change across all grid cells of the species’ current European range; projections were averaged between the two circulation models and calculated separately for each RCP scenario (see Supplementary Fig. 4 for a data overview). Restricting the vulnerability quantification to a species’ current distribution may overestimate extinction risk because areas outside the current range, which may become suitable in the future, are not considered. However, this approach accurately quantifies species’ exposure to projected changes in climatic conditions and avoids several simplifying assumptions, for example, about species’ dispersal ability, non-analogue climates or novel biotic interactions that are particularly problematic in range projections beyond current distributions17.

Changes in climatic suitability were very closely correlated between the two RCP scenarios (n=709 species, r=0.97), but were only weakly related to current climatic niche breadth (n=709 species, climatic hypervolume, r<0.5 for both RCP scenarios; OMI climatic niche breadth, r<0.4 for both RCP scenarios). Changes in climatic suitability at the regional scale (that is, within the nine grid cells adjacent to each network’s study region) were closely correlated to range-wide suitability changes (n=709 species, r=0.74 for both RCP scenarios). We used range-wide suitability changes because they are more representative for a species’ vulnerability across its range and are less sensitive to local climatic projections.

Mutualistic pollination and seed-dispersal networks

We matched data on climatic niche breadth and vulnerability with empirical data of biotic specialization derived from eight quantitative pollination and five quantitative seed-dispersal networks, each recorded within a different region in central Europe (Supplementary Data 1). We included these networks in the analyses because they report comprehensive data on interaction frequencies between species pairs, consistently recorded by direct plant observations over several months (Supplementary Data 1). Interaction frequencies equal the number of visits of an animal to a plant species. To cover all potential interaction partners of a species within the regional context, we summed interaction frequencies over time and space within each region as most networks were recorded on repeated visits at several localities within each region. Interaction data from other animal taxa than bees, butterflies, hoverflies and birds were excluded from the original networks. As it is generally the case for mutualistic networks, we lack information on the actual functional dependences of pollinators and seed dispersers on their foraging plants and vice versa and, thus, assume that interaction frequency is closely associated with the reciprocal functional importance of the interacting species18. Overall, the 13 pollination and seed-dispersal networks were characterized by a low to intermediate degree of specialization; complementary specialization (H 2 ') ranged from 0.29 to 0.47 (mean=0.38, see Supplementary Data 1 for a compilation of other standard network metrics for the 13 mutualistic networks).

In total, we were able to derive independent data of biotic specialization, climatic niche breadth and vulnerability to climate change for 363 insect (196 bees, 70 butterflies, 97 hoverflies; pollinators hereafter), 51 bird (seed dispersers hereafter) and 295 plant species (Supplementary Data 2). These species comprised most of the species of the respective study taxa (bees, butterflies, hoverflies, birds and plants) in the original networks (mean across networks: 88% (range: 74–100%), see Supplementary Data 1). Including species with deficient occurrence data and/or other animal taxa for which no distribution data were available (for example, other Diptera or Coleoptera that had been sampled for a few of the original pollination networks) resulted in qualitatively identical estimates of biotic specialization for the analysed plant and animal species.

Network analysis and simulations of species coextinctions

We measured plant and animal specialization in two ways based on the number and uniqueness of interaction partners within each network. First, we calculated the effective number of partners. It equals the diversity eH of interaction frequencies per species (based on the Shannon index H) and is equivalent to the number of partners if each link was equally common66. This metric is positively related to the total interaction frequency of plants (log–log scale, n=591, r=0.68) and animals (n=1009, r=0.76). On average, plants had a higher diversity of partners than animals (mean eH±1 s.e., 4.20±0.17 versus 2.61±0.10). Second, we calculated complementary specialization d’ as the deviation between the observed interactions and partner selection according to species’ total interaction frequencies67. The metric ranges from 0 for a generalist species (sharing partners with many others) to 1 for a fully specialized species (with a unique set of partners). d’ was unrelated to total interaction frequency (plants, r=0.03; animals, r=0.06) and was similar for plants and animals (0.34±0.01 versus 0.32±0.01).

We tested the statistical associations between species’ biotic specialization within each network and climatic hypervolume (square-root transformed), OMI climatic niche breadth (geometric mean of variances along the first two OMI ordination axes) and projected changes in climatic suitability (median change across a species’ current European range) with linear mixed-effects models; error distributions of all models did not deviate from normality. We fitted statistical models including main and interaction effects of the respective climatic predictor and trophic level (animal versus plant) on biotic specialization. In all models, we accounted for random variation due to network identity, plant and animal taxonomy on the model intercepts (taxonomic levels: family, genus, species). To account for variation in the performance of species distribution models and, thus, for the uncertainty of projected changes in climatic suitability, we weighted the respective linear mixed-effects models with the accuracy of the respective species distribution model as given by the TSS max value64 for each species (Supplementary Data 2). In the interest of comparability, we z-transformed realized climatic niche breadth (climatic hypervolume, OMI climatic niche breadth) and changes in climatic suitability before the statistical analyses. In order to examine network-specific relationships between biotic specialization and climatic variables, we additionally tested effects of climatic predictors on biotic specialization in models accounting for both random intercepts and slopes of network identity, separately for plants and animals.

Based on the projected impacts of climate change on plant and animal species, we modelled effects of climate change on each network. We simulated secondary species extinction as a consequence of the sequential loss of plant and animal species, respectively11,12. The order of species loss followed the projected changes in climatic suitability; thus, we first removed the species experiencing the largest decline in climatic suitability across the current European range, followed by the removal of the species with the second largest decline until the least vulnerable species had been removed. In the simulation model, species became secondarily extinct once they had lost at least 25, 50 or 75% of their interaction events in respect to the original network, assuming that interaction frequencies are proportional to the functional dependences of animals on plants and vice versa18. These thresholds for secondary extinction are probably more realistic than the assumption that all interaction events have to be lost before a species goes extinct15. We further assumed that species could reallocate lost interactions to other persisting species in the network accounting for the flexibility of partner choice in interaction networks13,16. We simulated two different rewiring scenarios under the assumption that no new species will enter the network. First, we reallocated a varying proportion of removed interactions to all persisting partners (constrained rewiring), proportional to the relative interaction frequencies of each species. Second, we reallocated lost interactions to all persisting species (unconstrained rewiring), relative to species’ total interaction frequencies. Thus, the first scenario assumes that species will be constrained in their interactions to their current partners, whereas the second scenario assumes that species are able to freely establish new links under future conditions. We varied the flexibility of species to reallocate lost interactions to persisting species between 0% (no reallocation) and 100% (reallocation of all interactions). In the scenario where all partners are interchangeable (100% flexibility), all interaction events must be lost to trigger secondary extinction as all lost interaction events are reallocated to persisting species. In a scenario of unconstrained rewiring and 100% flexibility, thus, all species need to go extinct to trigger secondary extinction. As this is a very unlikely scenario, it was omitted from our simulations.

For each network and simulation scenario, we quantified network sensitivity to plant and animal species extinction by the area above the secondary extinction curve; the metric ranges from 0 (no species go secondarily extinct) to 1 (all species go secondarily extinct after removing a single species) and is conversely equivalent to network robustness (the area under the extinction curve12). We computed the difference in network sensitivity to plant and animal extinction for each RCP projection and each coextinction, rewiring and species’ flexibility scenario. We further compared network sensitivity for the different simulation scenarios between projected climate change and 200 random sequences of species loss; here network sensitivity was averaged across iterations of random extinction for each simulation scenario. We tested whether the risk of secondary animal versus secondary plant extinction differed between climate change and random extinction using two-sided, pair-wise t-tests.

Data availability