1Background

Climate change has affected distribution and phenology of a wide range of taxa [1] and is forecasted to cause high rates of species extinctions [2]. In order to track climate change, species are expected to shift their range toward higher latitudes and elevations. Considering these movements, species-specific responses to climate change have the potential to create temporal or spatial mismatch between interacting species [3]-[5]. Due to biotic interdependences within trophic networks, loss or addition of species in local plant and animal assemblages may cause cascading changes [6],[7]. Studies forecasting the effect of climate change on species turnover within communities generally overlook the potential effect of biotic interactions that link species fate to each other [8],[9]. Hence, the consideration of the effect of species interactions represents a significant challenge for forecasting the future consequences of climate change [9],[10]. A first step to understand how interactions may change in a warmer future is to investigate shift in the functional structure of communities [11].

Pollination is a key ecosystem service and an essential component for the maintenance of biodiversity in ecosystems [12]. From the coevolution of pollinator and flowering plant, intricate mutualistic networks emerged within communities [13]. Functional complementarity suggests that a higher diversity of pollinators contributes to increased pollination success of the plants or, in turn, that a higher diversity of flowers may better sustain a diverse guild of pollinators [14]. However, climate change is expected to affect the co-occurrences of plant and pollinator species in space due to climate tracking or local extinction and thus disrupt such networks. Species-specific range shift of plants [15] and pollinators [16]-[20] suggests that distribution overlaps between plants and pollinators will not remain the same in the future [21]. In contrast to highly specialized plant species (e.g., [22]), the pollination success of a vast majority of plant species is determined by the availability of a functionally appropriate flower visitor, with adequate features such as body size or proboscis length [23],[24]. The floral morphologies in plant communities are thus tightly linked to the functional structure of the available pollinators through complementarity [14],[25]-[27]. In this context, Lavorel et al. [28] suggested that the effect of climate change on pollinator communities might affect the frequency of functional groups in communities more than composition per se and that this will drive change in plant community structure.

Among the world’s ecosystems, mountains and their unique biota are predicted to be disproportionally exposed to climate change with potential shifts in communities functional structure [8],[29],[30]. Bumblebees are the most efficient and active pollinators in alpine ecosystems [31], due to their capacity to generate heat [32] and fly over long distances [33]. Studying shift in bumblebee communities and functional structure informs on future plant-pollinator interaction in mountains. Bumblebees have evolved functional characters, like abundant hairs or muscles generating heat prior to flying enabling them to forage in cold environments [34],[35]. These adaptations to cold make them potentially more sensitive to climate change. They also display proboscis of varying length which mediates the type of flowers they interact with [36]. Bumblebees with longer proboscides visit deeper flowers than those with short proboscides [36], but long-tongued species are currently restricted to lower elevation [37]. Climate change may greatly disrupt the established structure of plant-pollinator interactions through change in the proboscides of available pollinator in communities.

Stacked species distribution models (S-SDMs) allow assessing the potential impact of climate change on biological communities [38],[39]. SDMs relate field occurrences to a series of environmental variables, in order to model a species realized environmental niche and its associated potential distribution in the landscape. S-SDMs then assemble communities issued from individual SDM predictions [40],[41]. Maps issued from SDMs are often binarized using threshold values based on different optimization techniques maximizing for instance sensitivity and specificity (reciever operating characteristic (ROC), see Elith et al. [42] for more information), enabling the conversion of probabilities into binary maps. An alternative is to use resampling from probabilities using a binomial distribution [41]. This probabilistic approach allows modeling realized endpoint of assembly (i.e., one possible assemblage out of other possible ones) accounting for stochasticity [39],[41] may also account for uncertainty in the modelled response curves [38],[41].

Here, we used high-resolution environmental maps of proximal environmental predictors derived from weather stations, three IPCC-based climate change scenarios, and a novel probabilistic ensemble S-SDM approach to assess comprehensively the possible effect of climate change on the structure of communities of an entire regional bumblebee pollinator fauna by the end of the 21st century (2070–2099 time period average). More specifically, 1) to assess the potential consequence of climate change on species and communities, we estimated the percentage of suitable habitat loss for each species compared to the current distribution of suitable habitats. 2) We compared the current functional structure of bumblebee communities with regard to average proboscis length, with the forecasted future communities under climate change.

2Results

2.1 Overall diversity

During the summers of 2009 and 2010, we sampled and identified 2,020 individuals across 149 sites, resulting in a total of 30 bumblebee species in the study area. The species from the Bombus sensu stricto subgenus occurring in the area (B. terrestris, B. lucorum, and B. cryptarum) have very similar workers, and identification of these workers to the species level can prove difficult for older individuals [43]. To avoid modeling bias resulting from poor taxonomic differentiation, we did not consider this group in further analyses. The Psithyrus subgenus (eight species) was also removed from the analysis as they parasitize the other Bombus species and are not frequent pollinators as they do not produce workers. Finally, five species—B. alpinus, B. gerstaeckeri, B. jonellus, B. subterraneus, and B. veteranus—had a too low number of occurrences (minimal number of occurrences needed = 15) to model their distributions.

2.2 SDM predictive power and community evaluation

The average area under the curve (AUC) value across all models ranged from 0.53 to 0.89 (average = 0.7, Table 1). For each species and from the 20 repetitions of model runs, the ensemble forecasting method conserved only the models with AUC values above 0.7, keeping only the models with sufficient predictive power to be used for projections following the recommendation of the biomod2 package [44]. The prediction success of bumblebee communities (success of predicting the right species in right places) ranged from 0.37% to 0.62% and was the lowest at middle elevation corresponding to where communities have the highest species richness [37], thus where the species pool is also larger (Figure 1).

Table 1 Predictive power of the models of the 15 bumblebee species measured with the area under the curve (AUC) of the ROC plot
Full size table
Figure 1
figure 1

Current and future (scenario A2, period 2085) species richness obtained after the 250 draws. The blue shows low species richness, whereas the red shows high species richness. The grey line shows the 1,000 m isoline.

Full size image

2.3 Species loss under climate change

No species was projected to lose 100% of its suitable habitat by 2085 for any scenario, but some species were forecasted to suffer a drastic reduction of their range. In particular, the species that occupy the coldest environment at high elevation are projected to lose more than 80% of their occupied range—B. sichelii (−92%) and B. pyrenaeus (−96%), B. mucidus (−83%), B. mesomelas (−91%), and B. mendax (−90%)—under the most extreme A2 scenario (Table 2). The loss is roughly similar between scenarios except for RCP3PD with limited range loss for those species (Table 2). In contrast, several species are forecasted to increase the size of potentially suitable habitat. This is especially the case for three species from the warm-adapted subgenus Thoracobombus, B. pascuorum (+97%), B. hortorum (+30%), and B. hypnorum (+405%). Overall, some groups were predicted to suffer greater loss of potential surfaces compared to other. The short-tongued clade Pyrobombus + Rhodobombus + Kallobombus is forecasted to lose a higher proportion of habitat surface (B. lapidarius and B. hypnorum excepted) compared to the long-tongued clade Thoracobombus + Megabombus, which tends to maintain its current range size with the exception of B. mucidus.

Table 2 Percentage of future to current occupied surfaces for all the species and the 15 bumblebee species individually represented following the three scenarios
Full size table

2.4 Functional structure under climate change

Under all climate change scenarios except the milder RCP3PD, the S-SDMs suggest an increase of the mean proboscis length in communities at high elevation, attenuating the current decrease in mean proboscis length along elevation (Figures 2 and 3). In a warmer future, bumblebee communities at the highest elevations will be more species rich and are expected to be composed of species with longer proboscides, currently restricted to lower elevation (Figure 1). Fewer differences in pollinator functional structure should be observed between lowland and high elevations at least with regard to bumblebees (Figure 4, for the variance in the 3 scenarios see Additional file 1: Figure S1).

Figure 2
figure 2

Current and future (scenario A2, period 2085) proboscis mean obtained after the 250 draws. The blue shows low mean proboscis length score whereas the red shows high mean proboscis length score. The grey line shows the 1,000 m isoline.

Full size image
Figure 3
figure 3

Change in suitable surface for the 15 Bombus species for the current period and the three different scenarios corresponding to the end of the century (2085). The numbers on the y axis correspond to the sum of pixels of 50 m × 50 m. The variability within each group corresponds to the 250 threshold resampling.

Full size image
Figure 4
figure 4

Relationship between mean proboscis length score and elevation for the three climate change scenarios: (a) A2, (b) A1B1, and (c) RCP3PD. The grey surface depicts the maximum and minimum current mean proboscis length projected over the study area. The blue solid and dashed lines indicate the future mean proboscis length together with uncertainties.

Full size image

3Discussion

Bumblebees have evolved functional characters that allow them not only to tolerate cold temperatures but to show optimal efficiency at high elevations [45]. This specialization to cold environmental conditions should in turn render them more sensitive to climate change. Our results indicate that in the western Swiss Alps at least five bumblebee species could lose more than 80% of their range by 2085 and will risk local extinction if warming continues. This is particularly the case of some high-elevation specialist species like B. monticola, B. sichelii, and B. pyrenaeus which are potentially not tolerant to high temperature stress or are less efficient compared to more competitive bumblebee species occupying warmer conditions [46]. Our forecasted shift in elevation range of bumblebee species parallels results of a recent colonization of northern Europe by bumblebee species linked to climate change [47]. In addition, we showed that under climate change species with longer proboscides will be able to colonize higher elevation. In turn, the average and variance in proboscis length in higher elevation communities will increase potentially affecting plant-pollinator interactions for instance by reducing pollination limitation at high elevation [48].

3.1 Species gain and loss under climate change

The tolerance to climatic conditions is phylogenetically conserved in bumblebees [37]. One clade exhibits tolerances toward much lower temperature compared to its sister group. The present study indicates that species from the short proboscis group occurring overall at higher elevation will suffer a higher loss of occupied surfaces compared to species with a longer proboscis. Some groups such as the Thoracobombus (B. humilis, B. pascuorum, B. ruderarius) clade may even increase their occupied range under climate change. Even if some species show a decrease of more than 90% of their range, no species included in the study was predicted to become extinct by 2085. However, as warming is expected to continue increasing beyond the time frame considered in this study, this ongoing range contraction may lead to the local disappearance of cold-adapted bumblebee species. In addition, we only considered more frequent species (>15 occurrences). It is likely that the rarest species with highly specialized niche will suffer the greatest loss. In particular, B. alpinus was only found in two sites above 2,500 m. The habitat of those already rare species is likely to shrink under climate change. Altogether, species forecasts are in accordance, with observed trend of expansion for some bumblebee species [49] and decrease in others [34], and suggest a species-specific response to climate change.

3.2 Homogenization of proboscis diversity under climate change along elevation

Plant-pollinator networks are strongly modulated by functional traits [24]. The variation in proboscis length of the different pollinator species mediates pollination networks and climate change may remodulate how plants and pollinators interact with each other [35],[48]. In addition to abiotic filtering, pollination filtering as proposed by Sargent and Ackerly [50] was suggested to currently exclude plant species from high elevation and maintain a pollination-generalist flora at high elevation [51]. However, climate change may remove such limitation. The average proboscis length of bumblebee communities decreases with elevation, because species from the long-tongued clade of bumblebees are overall less tolerant to low temperatures and thus restricted to lower elevation [37]. In contrast, bumblebee communities at high elevation are overall composed of species with shorter proboscides that are less able to forage flowers with long spurs [37]. The same pattern was also recorded for plant species where it has been shown that simple actinomorphic flowers were more frequent at higher elevation [52]. Our results suggest that the bumblebee species now occurring at lower elevation and having longer proboscides will track climate change and reach the current alpine landscape and increase the diversity of pollination service. In turn, this functional shift of pollinator assemblages may facilitate the establishment of plant species with more specialized pollination types, which are currently restricted to lower elevation [51].

3.3 Caveats of the current modeling approach

Pollinators influence plant species distribution and community assembly, but the reverse is also expected. Here, we did not consider the floral resource influence on bumblebee community structure. However, climate change may affect the local survival of bumblebee species in interaction with availability of floral resources. Aldridge et al. [53] showed that climate change may induce longer periods of low flowering abundance in the middle of the summer season that could negatively affect pollinators that are active throughout the season such as bumblebees. In parallel, shifts in flowering peaks within habitats might create mismatches between floral resources and demand by pollinators [53]. In addition, bumblebee species with long proboscides tend to have a more restricted trophic niche compared to species with short proboscides [54]. As a consequence, long-tongued bumblebee species may be more sensitive to global change perturbation compared to more generalist species. Using historical data, Dupont et al. [55] showed a consistent and dramatic decline in species richness and abundances of specialized bumblebees (with long proboscis) throughout the flowering season of red clover, while short-tongued species were largely unaffected. Future models should account for feedback processes between plants and pollinators under climate change. Finally, the potentially impact of climate changes and especially weather means and extremes could lead to individualistic responses which are very hard to predict as demonstrated with alpine butterflies [56].

4Conclusions

Climate change may affect species distribution and community composition either directly through abiotic tolerances to environmental conditions or because of change in biotic interactions. Here, we show that in addition to causing the decrease of several bumblebee species, climate change likely affects the functional structure of communities. In turn, this may impact plant pollination network along elevation under warmer climate. Here, we investigated the potential effect of climate change on pollination as one example of biotic interactions, but further studies are required to investigate other types of biotic interactions and achieve an ecosystem-wide assessment of climate change impacts.

5Methods

5.1 Study area and species sampling

The study area is located in the western Swiss Alps, at elevations ranging from 1,000 to 3,210 m above sea level, with temperate vegetation belts, i.e., montane, subalpine, alpine, and nival. Lowland sites (below 1,000 m) have not been sampled as they are mostly forested and intensively cultivated areas with a strong human influence on the plants and insect communities. The sampling plot selection was made following a balanced random-stratified sampling design, based on elevation, slope, and aspect, and considering only regions outside of forested areas [57]. During the months of June to September in 2009 and 2010, 149 bumblebee communities were sampled. Sites were visited every 3 weeks during hours of high bumblebee activity (10:00–17:00) and only under good weather conditions (i.e., little wind, sunny, and temperatures above 15°C). Visiting the site several times allows a robust sampling of the communities independently of the phenology. Each site was visited for 45 min, during which bumblebees were sampled in a 50 m × 50 m area. During the summers 2009 and 2010, we sampled and identified 2,020 individuals. Among these, we found 28 species. We considered only the non-parasitic species that produces colonies and whose worker greatly contributed to plant pollination. We also excluded species with less than 15 occurrences. We ranked the species according to proboscis length as in Pellissier et al. [37]: 1, very short proboscis (6.5–7.25 mm); 2, short proboscis (7.25–8 mm); 3, medium proboscis (8–8.75 mm); 4, long proboscis (8.75–9.5 mm); and 5, very long proboscis (9.5 mm and above).

5.2 Environmental variables

Recent studies have highlighted the need to produce more accurate thermal predictors to model the distribution of mountain flora and small fauna [58],[59]. Using the same approach as Zimmermann and Kienast [60], we modeled temperature at a resolution of 25 m from 188 temperature stations in Switzerland. We calculated the mean temperature for each month of the growing season (June to August) from 1981 to 2009 according to monthly lapse rates. Using the digital elevation model and the monthly lapse rates, we then projected the mean monthly temperatures for the whole study area. We assessed the local variation, adding the residuals interpolated over the study area and thus decoupling partly the temperature from elevation. The 1981–2009 period corresponds to the new reference period used by the Center for Climate System Modeling (http://www.c2sm.ethz.ch/) and is related to new regional climatic warming scenarios. We calculated the amount of solar radiation received in each month of the growing season (June, July, and August) in each pixel with the spatial analyst tool in ArcGIS 10 [61]. To fit the climatic layers with the sampling sites, we aggregated the temperature layers at a resolution of 50 m using an averaging function before extracting the data corresponding to each sampling point.

5.3 Climate change scenario

We used three different climate projection scenarios (A2, A1B, RCP3PD) averaged for three different models coupling regional and global climatic models (ARPEGE-ALADIN, ECHAM5-REMO, HadCM3Q0-CLM) using three averaged time periods (2021–2050, 2045–2074, 2070–2099) developed in the Swiss Climate Change Scenario CH2011 project from the Center for Climate Systems Modeling (http://www.c2sm.ethz.ch/). The Swiss Climate Change Scenarios CH2011 provide state-of-the-art projections of how climate may change over the 21st century. They are based on new generations of climate models with higher resolution combining global climate models (e.g., ECHAM5) and regional climate models (REMO). They also provide anomalies for every weather station in Switzerland [62]. These anomalies have been interpolated following the same techniques used for temperature and added to the actual maps. We used a total of 188 temperature stations for the interpolations of the anomalies.

5.4 Community modeling

We modeled individually the species distributions using presence and absence data with two variables: temperature average and solar radiations. We used three different modeling techniques implemented in the biomod2 package [44] and following Engler et al. [8]: a generalized linear model (GLM, [63]), a generalized additive model (GAM, [64]), and a random forest model (RF). For each technique, we used the default settings in biomod2, as these were optimized for SDMs [44]. We randomly split the dataset into 20 partitions, successively calibrating the models using 80% of the data and sequentially predicting the species distribution and communities based on the remaining 20% of the data. We measured the SDM predictive power for each species as the AUC of a ROC [65] using a 20-fold split sample procedure (modeling on 80% of the data, predicting on 20%). We applied an ensemble forecasting framework by averaging all single model projections, by AUC [66]. To build prediction, we used an the ensemble forecasting method from the biomod2 package, keeping only the models with AUC values above 0.7 and averaging the remaining models [44].

We used the probabilities predicted by SDMs as actual probabilities following Pellissier et al. [41]. Species may display different tolerances to environmental conditions while still governed by stochastic processes. To obtain predictions of species distributions as binary presence or absence from the species probabilities of occurrence, we used the rbinom function in R. This function draws presences and absences from the probabilities out of a binomial distribution, as previously used in Dubuis et al. [67] and Pellissier et al. [41]. We repeated this process 250 times and produced the same number of predicted potential communities.

We projected the community models over the study area in the current period as well as considering the three climate change scenarios and the three time periods. We subsequently applied a mask representing anthropized areas (cities, roads) and forests as a habitat that was not covered by the current sampling. Species were assumed to have unlimited dispersal capabilities, which is a fair assumption considering the bumblebee flight ability. We calculated the number of pixels occupied by each species in the current climate as well as under the nine climate change scenarios for all runs. We computed the ratio between current and future potential surfaces. Finally, we computed the average of all runs to obtain the main tendency. We evaluated whether some lineages were more likely to benefit or lose from climate warming. Based on the projected communities over each pixel of the study area in the current period as well as considering the three climate change scenarios, we measured the proboscis length average in each pixel and for all runs. Proboscis length for each species was obtained from Pellissier et al. [41]. We related the current and future average and variance in proboscis length to the elevation gradient to evaluate the effect of climate change on the bumblebee functional structure from lowland to the alpine belt.

Additional file