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).

Species	GLM	RF	GAM	Mean AUC
B. hortorum	0.676	0.570	0.394	0.546
B. humilis	0.891	0.801	0.764	0.819
B. hypnorum	0.930	0.676	0.621	0.743
B. lapidarius	0.754	0.691	0.677	0.707
B. mendax	0.776	0.704	0.664	0.715
B. mesomelas	0.892	0.672	0.703	0.756
B. monticola	0.759	0.686	0.625	0.690
B. mucidus	0.834	0.754	0.507	0.699
B. pascuorum	0.903	0.882	0.861	0.882
B. pratorum	0.693	0.570	0.641	0.635
B. pyrenaeus	0.854	0.800	0.808	0.821
B. ruderarius	0.747	0.668	0.692	0.702
B. sichelii	0.837	0.722	0.676	0.745
B. soroeensis	0.569	0.429	0.422	0.473
B. wurfleinii	0.743	0.571	0.629	0.647

2.3 Species loss under climate change

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).

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].

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.