Introduction

Climate change is undeniably one of the main driving forces behind the ongoing reorganisation of biological communities (Parmesan and Yohe 2003; Walther et al. 2002) and may put many species at risk of extinction in the near future (Urban 2015). Numerous studies have provided evidence that climate change affects biodiversity in several different ways, including phenological shifts (Cohen et al. 2018; Parmesan and Yohe 2003), changes in abundance (Parmesan 2006; Sturm et al. 2001) and evolutionary responses (Parmesan 2006; Rank and Dahlhoff 2002). Another, widespread, consequence of climate change is a shift in the distribution of species towards previously colder areas, i.e. poleward and uphill range shifts (Barton et al. 2016; Chen et al. 2011; Hickling et al. 2005; La Sorte and Thompson 2007; Parmesan 2006; Perry et al. 2005). At the local scale, it results in a gradual replacement of cold-adapted species by warm-adapted ones. This species turnover is often revealed by an increase in the community temperature index (CTI), a community weighted mean of species’ temperature preferences (Devictor et al. 2008). While climate change plays a considerable role in reshaping communities, so far, the main cause of biodiversity decline is the direct effect of human activity, especially habitat loss (Sala et al. 2000). As climate and land-use changes are expected to interact to shape future species’ distribution (Jetz et al. 2007; Marshall et al. 2018), it is imperative to examine how they together contributed to long-term biodiversity changes.

Mountain regions are particularly at risk with regard to the effect of climate change (Dirnböck et al. 2011; La Sorte and Jetz 2010). They often host specialised and endemic species, adapted to a short growing season and low average annual temperature, which can be at risk of extinction in a warmer climate (Laurance et al. 2011). Along altitudinal gradients, we typically observe an altitudinal shift of species’ distributions (Grytnes et al. 2014; Kelly and Goulden 2008) with variable lag behind the actual climate change (Alexander et al. 2018). This can have far-reaching consequences on species and community persistence since asynchronous elevation shifts can locally reshape competitive (Alexander et al. 2015) and trophic interactions (Pyke et al. 2016). Moreover, distribution shifts may be impossible for species already at their upper elevation limit, making mountainous species highly vulnerable to climate change (Elsen and Tingley 2015). Interestingly, some studies demonstrated that land-use change, alone (Ameztegui et al. 2016) or in interaction with climate change (Guo et al. 2018), also contributes to observed altitudinal shifts.

Investigating the joint effect of climate and land-use change on biodiversity is not straightforward because it generally requires monitoring long-term changes in species richness and community composition over large environmental gradients. Historical ecological data can be a valuable source of baseline data to study how biological communities have changed over time (Vellend et al. 2013), especially in combination with information on local climatic and land-use trends. Fortunately, there exist high-quality distributional data of bumblebee species in Scandinavia from the 1940s–1960s (Løken 1973), including in mountainous regions. Bumblebees are important providers of pollination services and are known in other regions to respond to climate (Kerr et al. 2015) and land-use change (Goulson et al. 2010). However, most studies reporting decline in bumblebees were carried out in lowland, intensively managed agricultural systems (e.g. Bommarco et al. 2012; Dupont et al. 2011; but see e.g. Ploquin et al. 2013). Thus, restricting analyses to subalpine areas allows describing changes along an altitudinal gradient, with a relatively natural land cover where subtle effects of land-use and climate change are not masked by the intensification of agricultural practices.

For the present study, we take advantage of a historical high-quality data set of bumblebee distributions (Løken 1973), and re-visited pairs of low and high altitude sites in Norway. We report changes in species richness as well as in CTI to quantify the change in the relative proportion of warm- and cold-adapted species. While an increase in CTI is expected as a result of a warmer climate, it is unclear how species richness could have changed in this 50-year period. Climate change may have contributed to an increase in the overall species pool because of the colonisation of new species moving northward (e.g. Martinet et al. 2015), but may have also caused the extinction of climate-sensitive species, especially cold-adapted species at their lower altitudinal margin. At high altitude, the balance between the colonisation of species previously restricted to lowlands and the local extinction of cold-adapted species that shifted even higher—if possible—is hardly predictable. Local changes in species richness also depend on the relative speed of these processes. For example, we may observe a temporary increase in species richness if there is a time lag in extinctions while colonisations occur at a faster rate. On top of that, land-cover change may have exacerbated or counteracted the effect of climate change. To disentangle more precisely the different processes, we used historical and present aerial photographs, as well as local meteorological data, to test the effect of site-specific characteristics on the variation in species richness and CTI change. We hypothesised that the change in species richness was greatest in sites where the magnitude of climate and land-cover change was maximum. Similarly, we expected CTI to have increased more in the sites that have experienced the largest warming. Altogether, these analyses will improve the comprehension of the processes driving biodiversity changes in subalpine areas in a context of global change.

Methods

Study design and bumblebee survey

Løken (1973) surveyed bumblebees at a large number of sites in Norway during the mid-1900s, and compiled distribution records based on her own and other collections. The records from Løken’s collections have later been digitized, and are avaliable though Artskart (https://artskart.artsdatabanken.no/). It is not possible to identify the exact sampling effort from Løken’s publication. Instead, we extracted the dates that bumblebee specimens had been collected, from the digitized records in Artskart. This gives us the minimum number of days that a certain site was visited (Supplementary material, Appendix S1). In a few sites (Øyer high, Øyer low, Lom low) that had been visited extensively as a part of a study focusing on flower visitation patterns rather than systematically surveying for bumblebee species in general (Løken 1949), we excluded days when no new bumblebee species had been collected. The majority (10) of sites had been visited in a single year, but eight sites had been visited more than one year (Appendix S1). The median number of days a certain site had been visited was 3.

We selected 18 sites to be re-surveyed in 2012 (Fig. 1). A criterion for selection was the presence of the plant Aconitum septentrionale, which was the subject of a parallel study, and occurs mainly in nutrient-rich, wet meadows and sub-alpine forests with high soil pH. However, since we surveyed bumblebees within a large area with a 1 km radius (see below), this specific criterion is unlikely to have any larger impact on the bumblebee species composition except for the occurrence of Bombus consobrinus which is specialized on A. sepentrionale. We selected sites in order to cover a large latitudinal and longitudinal gradient (Lat.: 60.42–62.62°, ca. 240 km; long.: 7.24–10.61°, ca. 180 km), and to form 9 pairs of sites where each pair consisted of one site at a high altitude—approximately 1000 m.a.s.l. (min = 730 m, max = 1000 m) and one site at a low altitude—approximately 500 m.a.s.l. (min = 440 m, max = 765 m) in each pair (Fig. 1). We also made sure that each site contained a variety of flower-rich habitats such as different types of grassland, road verges, semi-open forest etc. The median distance between two sites in a pair was 12.3 km (min = 2.6 km, max = 44.4 km), and the median altitudinal difference was 352 m (min = 235 m, max = 470 m). The high-altitude sites were located within the upper part of the subalpine forest zone, close to the forest line which is around 1000 m.a.s.l. in this region (Moen 1999). We identified the locations and checked the suitability of sites to be re-surveyed in the field in July 2011.

Fig. 1
figure 1

Map of surveyed sites in Norway

Full size image

In some cases, it was difficult to identify the exact location of the site for the historical records, which typically consisted of the name of a farm or village, rather than any exact coordinates. The identification of sites was however facilitated by the fact that the altitude was often stated. To account for this uncertainty, we took a conservative approach, and surveyed bumblebees in multiple flower-rich habitats within a circular area with a radius of 1 km around the most likely location for the historical records. This means that we are less likely to under-estimate current bumblebee species richness than the historical surveys. Prior to the field visits, potentially flower-rich habitats were identified from aerial photos and by traveling along roads within the 1 km circle by car. Each site was then visited twice between July 6 and July 27 2012. In the field, two experienced field entomologists visited as many different flower-rich habitats within the 1 km circle as possible during 2 h (excluding time for transport between habitats) and recorded the abundance of all bumblebee species encountered. The bumblebee surveys were carried out by experienced entomologists with expertise not only in bumblebee identification but also with large field experience of their ecology. This ensured that the search efforts were as efficient as possible, as these persons were trained to identify suitable habitats of different bumblebee species. In addition, flower-visiting bumblebees on Aconitum sepentrionale was observed during 45 min in each site, during the same day as the surveys. The data from the flower visitation and the 2-h surveys were pooled for the analyses. Most species could be identified in the field, but in a few cases bumblebee specimens were collected and identified later in the laboratory. To allow comparison with historical data, we discarded abundances in further analyses and used species occurrences only. Data are available from the Figshare repository: https://doi.org/10.6084/m9.figshare.6833531.v1.

Many of the rarest species of the historical data were not observed in 2012 (see “Results” section and Table 1). This might potentially be caused by a failure to detect those species due to an inadequate sampling protocol relative to their low detection probability. To account for this potential bias, we repeated all analyses (see below) after removal of the seven species that were present in less than three sites in the historical data, therefore focusing only on the common species. Moreover, because simple counts of the number of observed species typically underestimate the true number of species, we aimed to test whether our results were consistent when using a non-parametric estimator of species richness instead of our raw count of species richness. Therefore, we repeated all our analyses of species richness replacing, in each site, the number of species encountered in the present survey by the Chao 1 abundance-based estimator of species richness (Chao 1987) computed with the “vegan” R package (Oksanen et al. 2015). Note that we could not derive richness estimators from the historical data because they require either abundance data or multiple surveys. In addition, we also computed accumulation curves of total species richness for historical and present survey, using the “iNEXT” R package (Hsieh et al. 2016), in order to ensure that our survey protocol captured the entire species richness in the region.

Table 1 Species detected in the historical and present surveys, and number of sites (low and high altitude) occupied by each species
Full size table

Climate and land-cover change

We obtained total daily precipitation and mean daily temperature from the Norwegian Meteorological Institute as interpolated raster grids with a 1 km × 1 km spatial resolution, based on the national network of weather stations (www.met.no, data available at: ftp://ftp.met.no/projects/klimagrid). Daily values from 1960 to 2012 were extracted at each survey site. As measures of climate change, we used the slopes of the linear trends of the mean annual temperature and the mean annual precipitation during this period. Temperature increased in all sites [mean slope = 0.026 °C/year (+ 1.46 °C over the entire period), min. = 0.003 °C/year (+ 0.17 °C), max. = 0.040 °C/year (+ 2.24 °C)], while precipitation increased in all sites except one [mean slope = 0.225 mm/year (+ 12.6 mm), min. = − 0.0007 mm/year (− 0.04 mm), max. = 0.560 mm/year (+ 31.4 mm)].

To assess the direction and intensity of land-cover change that occurred during the ca. 50 years separating the two surveys, we used aerial photographs from the Norwegian Mapping Authority taken in 1958—1968 and 2008–2013. Using ArcGIS 10.3 (® ESRI), we defined a buffer of 1 km radius around each survey site and extracted for each period the proportion of land in the buffer covered by: arable land, forest, grassland, wetland, built-up area, clear-cut forest or open (heath or scrub) area above the tree-line. There was high variability in land-cover among sites (Supplementary material, Appendix S2-A). Some were predominantly covered by forest (e.g.: Oppdal, high altitude, in the 1960s: 75% forest), others by grassland (e.g.: Geilo, high altitude, in the 1960s: 95% grassland) or had a more even distribution of land-cover types (e.g.: Lom, low altitude, nowadays: 40% arable land, 21% forest, 31% grassland). We calculated an index of dissimilarity, based on the Euclidian distance, between the historical and present land-cover. This index could thus range from 0 for sites whose land-cover did not change to 1 for sites whose land-cover has been entirely replaced (mean land-cover change across all surveyed sites: 0.23, min. = 0.023, max. = 0.55). This index therefore represented the magnitude of land-cover change between the two periods irrespective of the direction of this change. Additionally, we computed a principal component analysis (PCA) based on the proportion of each land-cover category (Appendix S2-B). For each site, we quantified the direction of land-cover change by calculating the difference along the two first PCA axes between the historical and present periods. These two principal components account for 48.8% of the total variation in land-cover.

Data analyses

From the survey data, we calculated for each site and for the historical and present surveys the species richness (the number of bumblebee species identified) and the community temperature index (CTI). CTI is a measure of the relative proportion of cold- and warm-adapted species in a community and has been commonly used to describe the species turnover that occurs as a result of climate change (Devictor et al. 2012). CTI was calculated as the mean species temperature index (STI) of all co-occurring species (Devictor et al. 2008). STI represents the average temperature experienced by a species within its geographical range. STI for all bumblebee species was extracted from Rasmont et al. (2015), where temperature data were based on 1970-2000 climate. We compared the estimated annual rate of CTI change—assuming that our historical data represent community composition in 1960—to the observed temperature trend to assess whether communities track or lag behind the actual climate change. In addition, for each species, we averaged the altitude, latitude and longitude of the sites where the species has been encountered in the historical and present surveys. This way, we aimed to describe shifts in species’ distributions that may have occurred over the 50-year period. We tested whether these mean values differed significantly between surveys using Welch’s t tests. Note, however, that since we recorded the same limited number of sites in both surveys, this analysis is unable to reveal broad patterns of range shifts and is only relevant to describe the distribution of species among our 18 surveyed sites.

In a first step, we described how species richness and CTI varied geographically and how they changed between the historical and present surveys. We used linear mixed models with either species richness or CTI as response variables, and the period of survey (historical or present), the altitude (low or high), the longitude and latitude of sites and the interaction between altitude and the period of survey as explanatory variables. Since surveyed sites were grouped by pairs of high and low altitude sites, we included the identity of site pairs as a random effect. Mixed models were computed in R using the “lme4″ package (Bates et al. 2015), with the significance of fixed effects assessed using the “lmerTest” package (Kuznetsova et al. 2017).

In a second step, we aimed to test whether the physical characteristics of the sites (altitude, magnitude of land-cover and climate change) could explain the magnitude and direction of the change in species richness and CTI between historical and present surveys. For this purpose, we computed linear models with, for each surveyed site, the net change between surveys in species richness (log-transformed) or in CTI as response variable. Here, we log-transformed species richness to standardise its temporal change whatever the level of species richness observed in the site. Explanatory variables were the magnitude of land-cover change (Euclidian distance), land-cover change along the first and second principal components, temperature trend, precipitation trend, and the altitude (low or high). Due to the relatively low number of sites, we could not test the interaction between land-cover change, climate change and altitude. In addition, in order to control for spatial autocorrelation between sites, we included a spatial correlation structure as a function of longitude and latitude. We tested several types of correlation structures and selected a Gaussian correlation structure because it yielded the lowest AICc. Models with correlation structures were fitted using generalized least squares with the “nlme” R package (Pinheiro et al. 2018). We additionally performed a post hoc analysis with the “lsmeans” R package (Lenth 2018) to test if the changes in species richness at low and high altitude were significantly different from 0.

Results

Twenty-three bumblebee species were identified in the historical survey. Only 16 were found in the present survey, and no new species were detected (Table 1). Accumulation curves show that we reached the plateau of species richness in the present survey, while more species probably remained undetected in the historical survey (Appendix S3). Among the 16 species that persisted, only one (B. cingulatus) showed a significant distribution change among our sites (P = 0.043), namely a northward latitudinal shift by 0.69°, i.e. ca. 77 km (Appendix S4). CTI had increased by 0.014 °C/year on average, corresponding to a lag of 0.012 °C/year compared to the actual temperature trend (0.026 °C/year on average). Looking at each sites individually, 11 out of 18 clearly lagged behind the observed climate change (Appendix S5).

Species richness and CTI were both on average higher at low altitude (Table 2 and Fig. 2, P = 0.017 and P < 0.001 respectively). In addition, there was an increase of CTI from east to west (Table 2, P = 0.047). Due to large variation between sites (Appendix S6), there was no consistent change in species richness between the historical and present surveys (Table 2 and Fig. 2, P = 0.428). In contrast, we observed that CTI increased over time both at low and high altitudes (Table 2 and Fig. 2, P = 0.007).

Table 2 Results of the linear mixed models that tested the effect of survey period, altitude and geography on species richness and community temperature index
Full size table
Fig. 2
figure 2

Variation of species richness (left) and community temperature index (CTI, right) between historical and present survey and between altitudes. Values are mean ± s.e.m. for low (blue and dotted lines) and high altitudes (red and solid lines). Species richness and CTI were both higher at low altitude. While species richness did not change between survey periods, CTI increased over time (see Table 2)

Full size image

The variation in the relative change of species richness among sites could to a large extent be explained (R2 = 0.78) by the combination of the altitude (P = 0.003) and the magnitude in the changes in land cover (P = 0.005), temperature (P = 0.021) and precipitation (P = 0.008) (Table 3). The direction of land-cover change along both PCA axes had, however, no effect on the change in species richness (Table 3, P = 0.766 and P = 0.963). We observed on average a positive change in species richness in high-altitude sites (+ 25.9%, P = 0.029), but a stable species richness (or slight decrease) in low altitude sites (− 20.1%, P = 0.074). There was a positive relationship between the change in species richness and the magnitude of land-cover change. Finally, the change in species richness was negatively correlated to the trends in temperature and precipitation (Fig. 3). The variation in CTI change could not, however, be explained by any of the site-specific variables (Table 3, all P > 0.08).

Table 3 Results of the linear models explaining the change in log-transformed species richness and in community temperature index (CTI) between the historical and present surveys
Full size table
Fig. 3
figure 3

Partial regression plots showing that the relative change in species richness was positively related to the magnitude of land-cover change, negatively related to the increase in temperature and precipitation, and was larger at high altitude. Full model results are shown in Table 3. The y-axis can be interpreted as the proportional change in species richness between the present and the historical survey. Each plot shows the predicted values when all other variables are kept at their median. Note that post hoc tests revealed that the change in species richness was significantly positive at high altitude (P = 0.03), while only on the verge of significance at low altitude (P = 0.07)

Full size image

Removal of the seven rarest species did not change the results qualitatively (Appendix S7). The only difference was that the negative relationship between the change in CTI and land-cover change on the second PCA axis became significant (P = 0.039). This shows that our results are robust to a possible failure to detect species with a low detection probability. Similarly, using Chao estimates of species richness instead of the raw species counts in the present survey did not affect the results concerning the drivers of the change in species richness (Appendix S8). The only difference here is that the difference in species richness between low and high altitude was not significant anymore.

Discussion

We observed changes in the composition and diversity of communities consistent with a combined effect of climate and land-cover change on bumblebee assemblages. The general increase in community temperature index (CTI) observed between the present survey and the 1960s reflects the gradual replacement of cold-adapted species by warm-adapted species, though lagging behind climate change in most locations. A temporal increase in CTI has been observed in different organisms including birds (Devictor et al. 2012; Lindström et al. 2013), butterflies (Devictor et al. 2012) and plants (Savage and Vellend 2015). This trend was expected under the hypothesis that species track their climatic niche by shifting their distribution, leading to a turnover towards more warm-adapted species at the local scale. In this regard, a latitudinal shift of bumblebee distributions was predicted (Rasmont et al. 2015) and has been consistently described worldwide (Kerr et al. 2015). Moreover, some species previously restricted to more southern areas have recently invaded northern Scandinavia (Martinet et al. 2015), which exemplifies the shift towards warmer communities in historically cold areas. However, we did not detect any new species compared to the historical inventory and there was no evidence for a latitudinal shift in the species we observed, suggesting that the increase in CTI was not primarily driven by the colonisation of species expanding their northern range margin.

Generally, we showed that climate change had an overall negative impact on bumblebee communities. Indeed, we observed a decrease in species richness in the sites that had the largest increase in temperature or precipitation. A possible explanation is that higher temperature had a direct effect on bumblebee mortality, for example because of more frequent heat waves or droughts (Iserbyt and Rasmont 2013; Rasmont and Iserbyt 2013). Since the precipitation trend was calculated for the entire year, an increase in precipitation may correspond to a longer period of snow cover in spring or fewer active days without rainfall in summer, two climatic trends that may be detrimental for bumblebees (Sanderson et al. 2015). Alternatively, we may observe an indirect effect of climate change through a changed composition of the vegetation. There is plenty of evidence for shifts in the composition of alpine vegetation communities in response to climate change (Cannone et al. 2007; Dullinger et al. 2012; Gottfried et al. 2012). Such a change in vegetation may have caused local extinctions of bumblebee species that were not compensated by colonisations of new species. Climate change can also affect vegetation and bumblebees asynchronously in such a way that it generates a phenological mismatch between bumblebees and their floral resources (Bartomeus et al. 2011; Ogilvie et al. 2017; Pyke et al. 2016). We noted, however, that temperature and precipitation trends showed a strong geographical pattern characterised by a lower increase in precipitation and a higher temperature increase towards the east (Appendix S9). It remains thus possible that their observed effects on species richness were in reality driven by another, spatially correlated, variable.

Despite the fact that we demonstrated an adverse effect of climate change on bumblebee species richness, overall we did not detect a decline of local richness. This is partly because it was counterbalanced by an increase of species richness at high altitude. This is in accordance with an altitudinal shift of low-altitude species, a typical consequence of climate change (Chen et al. 2011; Walther et al. 2002), while species of higher altitude were able to persist. A possible process is that climate change made high-altitude areas more suitable as they became warmer (Franzén and Öckinger 2011), while human disturbance, usually higher at low altitude, would have driven extinctions in lower areas. An altitudinal shift of bumblebee distributions has been documented in several areas of the world (Biella et al. 2017; Kerr et al. 2015; Ploquin et al. 2013; Pyke et al. 2016). More specifically, Ploquin et al. (2013) described in bumblebee communities of northern Spain a pattern similar to the one we observed, with a range retraction at low altitude and a net gain of species at higher altitude. Similarly, Biella et al. (2017) showed an uphill shift in the lower range of Bombus alpinus in the Alps, but that was not compensated by a similar shift at high altitude because of the persistence of glaciers and harsh weather conditions. Finally, Pyke et al. (2016) recorded altitudinal shifts in the distribution of North-American bumblebees, although they also noted that it resulted in a reduction of synchrony with plants’ phenology. None of the species we observed exhibited a significant upward shift here (Appendix S4), although a majority (10 out of 16 persisting species) showed a tendency to be detected on average at higher altitude than in the historical data. Altogether, this suggests that elevation shifts are common in bumblebees inhabiting montane habitats.

Another process that buffered the loss of local species richness caused by climate change was the effect of local land-cover change. Indeed, we observed that the diversity of bumblebee species seemed to benefit from large land-cover change, whatever the direction of this change, for reasons that are unclear so far. This positive relationship between the magnitude of land-cover change and the change in species richness is intriguing because there is no consistent trend towards more natural or heterogeneous habitats that could explain this result (Andersson et al. 2013). In a management perspective, understanding the link between land-cover change and the modification of bumblebee communities may help maintaining or improving bumblebee-friendly land-use to compensate the detrimental effect of climate change.

Long-term monitoring data distributed over large environmental gradients are rare, but are invaluable for describing biodiversity trends in relation to slow-acting processes. They are especially important when biodiversity changes may not always occur concomitantly with the underlying drivers because of extinction debts or temporal lags in climate tracking (Jackson and Sax 2010; Talluto et al. 2017). Here, by replicating inventories from 50-70 years ago, we were able to show evidence for an effect of both climate and land-cover change on subalpine bumblebee assemblages. These factors are likely to interact with each other to shape future species’ distributions and community composition in a large range of organisms (Guo et al. 2018; Jetz et al. 2007; Marshall et al. 2018). Although we were not able to exactly repeat the sampling procedure from the original surveys due to lack of information, our methodology with directed searching in a variety of habitat types, by experienced field entomologists was especially designed to maximize the likelihood of finding even the rarer species. The observation that species accumulation curves reached the plateau of species richness in the present survey, indicates that this strategy was likely successful. The fact that the results remained qualitatively the same both when we used estimated rather than observed species richness, and when we omitted the seven rarest species from the analyses indicates that our conclusions are robust to these differences in sampling effort.

In mountains, land-cover is partially determined by climatic factors, since temperature constrains the upper limit of forest cover, making it challenging to disentangle direct and indirect effects of climate change (Grace et al. 2002; Guo et al. 2018). Another important aspect, that we could not test here, is the synchrony in the response to climate change of interacting species (HilleRisLambers et al. 2013). In the present case, flowering plants that constitute an essential food resource for bumblebees must be able to shift simultaneously their distribution uphill so that plant-bumblebee interactions can be maintained in the long term (Pyke et al. 2016). Generally, elevation shifts caused by climate change have the potential to perturb biotic interactions, by facilitating the invasion of novel competitors previously restricted to lower altitudes (Alexander et al. 2015) or by disrupting trophic interactions (HilleRisLambers et al. 2013). Moreover, in low mountain areas, where the potential for elevation shifts is limited and where the available area rapidly decreases as species move up, failure to track climate change towards higher altitude may cause the regional extinction of the most sensitive species (Elsen and Tingley 2015; Wilson et al. 2005). Here, despite no overall change at the local scale, 7 out of 23 species found in the 1960s were not observed anymore. These species are known to be normally flying in July, when our survey was conducted (Løken 1973; Ødegaard et al. 2015), suggesting that they have gone extinct from the area. These species were not among the most cold-adapted, contrary to we may have expected in a context of climate change. However, they were already rare in the historical survey and generally habitat specialists (Ødegaard et al. 2015)—three are typical alpine species (B. flavidus, B. hyperboreus, B. pyrrhopygus) and three other species are strictly specialised in open lowland habitats (B. distingueus, B. humilis, B. ruderarius) – which could have made them sensitive to any type of disturbance. This may indicate that the environmental changes occurring in subalpine areas have overall a negative impact on biodiversity that needs to be investigated further. A deeper understanding of the interactions between climate, topography, land-use change and species’ ecology may help to predicting the response of subalpine communities to global changes.