Large variability in response to projected climate and land‐use changes among European bumblebee species

Bumblebees (Bombus ssp.) are among the most important wild pollinators, but many species have suffered from range declines. Land‐use change, agricultural intensification, and the associated loss of habitat have been identified as drivers of the observed dynamics, amplifying pressures from a changing climate. However, these drivers are still underrepresented in continental‐scale species distribution modeling. Here, we project the potential distribution of 47 European bumblebee species in 2050 and 2080 from existing European‐scale distribution maps, based on a set of climate and land‐use futures simulated through a regional integrated assessment model and consistent with the RCP–SSP scenario framework. We compare projections including (1) dynamic climate and constant land use (CLIM); (2) constant climate and dynamic land use (LU); and (3) dynamic climate and dynamic land use (COMB) to disentangle the effects of land use and climate change on future habitat suitability, providing the first rigorous continental‐scale assessment of linked climate–land‐use futures for bumblebees. We find that direct climate impacts, although variable across species, dominate responses for most species, especially under high‐end climate change scenarios (up to 99% range loss). Land‐use impacts are highly variable across species and scenarios, ranging from severe losses (up to 75% loss) to considerable gains (up to 68% gain) of suitable habitat extent. Rare species thereby tend to be disproportionally affected by both climate and land‐use change. COMB projections reveal that land use may amplify, attenuate, or offset changes to suitable habitat extent expected from climate impact depending on species and scenario. Especially in low‐end climate change scenarios, land use has the potential to become a game changer in determining the direction and magnitude of range changes, indicating substantial potential for targeted conservation management.


| INTRODUC TI ON
The global food system greatly benefits from insect pollination that provides an estimated economic value between €9.7 bn (Vallecillo et al., 2019) and ~€14 bn in the European Union  and ~€153 bn globally (Gallai et al., 2009). While the area of insectpollinated crops is on the rise (Aizen et al., 2008), a decline of wild pollinator abundance and diversity has been observed across scales (Goulson et al., 2008;Potts et al., 2010;Powney et al., 2019), driven by agricultural intensification Grab et al., 2019;Newbold et al., 2015;Potts et al., 2016;Woodcock et al., 2017) and a changing climate (Kerr et al., 2015;Soroye et al., 2020). Yield losses due to pollinator decline have been observed (Deguines et al., 2014;Garibaldi et al., 2011), thus threatening future food supply (Potts et al., 2016) and the quality of human diet, as the crops concerned provide high shares of proteins, vitamins, and minerals (Eilers et al., 2011;Klatt et al., 2014).
Bumblebees (Bombus ssp.) are among the most important wild pollinators in temperate regions due to their relatively high abundance and coevolved pollinator associations with many native plant species and their generalist life-history strategy (Goulson et al., 2008). Moreover, bumblebees are able to pollinate at relatively low temperatures and in bad weather conditions, thus being more efficient than managed honey bees (Goulson, 2010). Their ability to 'buzz-pollinate' makes them the most important pollinators for certain flowering crops (e.g., tomatoes) that rely on this form of insect pollination (De Luca & Vallejo-Marín, 2013). However, bumblebees are sensitive toward several environmental changes, such as the loss of floral resources and increased pesticide usage associated with agricultural intensification (Carvell et al., 2006;Rundlöf et al., 2015), resulting in declining populations across spatial scales. Potts et al. (2015) report declining population trends for roughly half of the 68 European bumblebee species (with 20 species remaining stable and only nine species increasing); a finding repeated for more than half of the 23 species present in the United Kingdom (with seven species increasing) (Comont & Dickinson, 2020). Similar trends can be observed in North and South America (Bartomeus et al., 2013;Cameron et al., 2011;Krechemer & Marchioro, 2020). Concerns have been raised that such trends will continue or even amplify in the future with ongoing climate and land-use change (Marshall et al., 2018;Rasmont et al., 2015;Sirois-Delisle & Kerr, 2018). An improved understanding of how individual environmental pressures, and particularly the interactions among them, affect bumblebee populations and bumblebee habitat is essential to develop conservation strategies that counteract declines in bumblebee distribution.
A rich body of literature exists demonstrating the impacts of a changing climate on bumblebee populations and ecological traits across geographic regions (e.g., Casey et al., 2015;Gérard et al., 2020;Kerr et al., 2015;Marshall, Perdijk, et al., 2020). On the European scale, Rasmont et al. (2015) showed a decrease in climatically suitable areas for 49-55 European bumblebee species in 2100, depending on the climate scenario considered. Similarly, Sirois-Delisle and Kerr (2018), who assessed potential range changes of 31 North American bumblebee species under four representative concentration pathway (RCP) climate scenarios, identified range losses for 15-30 species depending on climate scenario and assumptions about dispersal ability. Recently, Soroye et al. (2020) identified a higher risk of local extinction rates for North American and European bumblebees, if climate conditions exceed historically observed ranges. However, continental-scale studies often do not account for the effects of land-use change (e.g., agricultural intensification, urbanization, or expansion of cropland and grazing land) on bumblebee habitat. When included, it is only through ad hoc scenarios and with little thematic detail (Marshall et al., 2018;Soroye et al., 2020). As a result, important drivers of the future distribution of bumblebee populations and associated pollination services are underrepresented in current assessments.
The knowledge regarding the impacts of land-use change and agricultural management on bumblebee ecology and habitat (e.g., fewer nests sites in intensive agricultural landscapes or reduced colony growth and queen production due to pesticide application) is well established from the plot-to landscape scale (e.g., Aguirre-Gutiérrez et al., 2015;Goulson et al., 2010;Senapathi et al., 2015;Whitehorn et al., 2012). However, the effects of these drivers on large-scale bumblebee distributions and their interaction with climate change are less understood. Land-use change, where included, was found to have non-significant (Kerr et al., 2015) or minor negative effects (Soroye et al., 2020) as a driver of shifts in historical distribution patterns at continental scale and is to date hardly accounted for in projections of future distribution at the continental scale Sirois-Delisle & Kerr, 2018). Marshall et al. (2018) were the first to include dynamic land-use changes in European-scale, predictive bumblebee distribution modeling, finding significant changes in range dynamics upon accounting for land-use variables in their models. Yet, the land-use changes they considered were limited, and had relatively low importance compared to climate change impacts; a finding they attributed partly to the difficulty of capturing land-use and land-management changes at a level of detail relevant for bumblebee ecology (Marshall, Beckers, et al., 2020).
Partly due to such methodological limitations, there have been no attempts yet to quantify potential future changes of bumblebee habitat area in a framework where different levels of climate change and various plausible land-use futures are combined to a coherent set of linked climate-land-use scenarios. Moreover, there has been little focus on the variability among individual species' responses to these interacting environmental pressures. However, an improved understanding of how individual species react to different levels of environmental change is crucial to identify locations that require more detailed analyses on the interactions between land use and bumblebees and to design targeted land-management policies that support as many species as possible under a changing climate. Here, we thus combine land-use change projections for Europe from a regional integrated assessment platform running a full set of socioeconomic and climate scenarios (the IMPRESSIONS Integrated Assessment Platform (IAP2); Harrison et al., 2019) with established species distribution models (SDMs) for 47 European bumblebee species developed by Polce et al. (2018). We do so in order to (1) quantify changes in bumblebee habitat suitability due to isolated climate and land-use change effects at the individual species level and (2) assess the impact of different land-use change futures on bumblebee habitat suitability under different levels of projected climate change. We specifically focus on the heterogeneous impact across bumblebee species to evaluate at which levels of climate change and for which species the choice of land use has the potential to become a game changer in a changing climate.

| Overview
We combined established knowledge on the distribution of 47 European bumblebee species (Table S1) with projections of climate and land-use change to quantify the species-specific response to a set of climate and land-use scenarios at the European scale ( Figure 1). Baseline bumblebee distributions were based on SDMs at 10 × 10 km spatial resolution (Polce et al., 2018) and adjusted to be consistent with land-use data from the projections. Future climate and land-use data were taken from CMIP5 climate model simulations (Taylor et al., 2012) 'Neill et al., 2017;van Vuuren et al., 2011), with socioeconomic developments adjusted to the European context (Kok et al., 2019). We applied the SDMs to the time periods representative for the years 2050 and 2080 to obtain estimates of the future distributions of the 47 bumblebee species. We established an expert-based mapping between the Coordination of Information on the Environment (CORINE) land cover and the IAP2 land-use classes to represent major changes in European land management in the SDMs. All extrapolations were repeated for a series of modeling experiments where climate predictors (CLIM), land-use predictors (LU), or both (COMB) were allowed to dynamically change. The details of the individual steps in the analysis are described in the following sections.

| Baseline bumblebee distribution
The potential distribution of 47 bumblebee species (Table S1) and SDM model coefficients were made available by Polce et al. (2018). The distribution maps indicate for each species (1) the probability of occurrence and (2) presence/absence at a spatial resolution of 10 x 10 km, based on a maximum entropy (MaxEnt) modeling approach (Phillips et al., 2006) that establishes a functional relationship between species occurrence and 22 environmental predictor variables (Tables 2 and 3). Species occurrence data for the years 1991-2012 underlying the MaxEnt models originate from the Atlas Hymenoptera  and consisted of validated presence-only bumblebee records, gathered from different data donors in Europe (see Rasmont et al., 2015 andSupplementary Material 7 in Polce et al., 2018 for further details), collated during the 'Status and trends of European pollinators' project . The environmental predictors were based on CORINE land cover for the year 2006 for land-related variables (Bossard et al., 2000), 'E-OBS' gridded meteorological data for bioclimatic variables (Cornes et al., 2018), and F I G U R E 1 Overview of the analysis. MaxEnt models from Polce et al. (2018) are adjusted through the exchange of land-related predictor variables based on IAP2 land use and projected to 2050 and 2080 under various combinations of climate (WorldClim v1.4 RCP scenarios) and land-use change (IMPRESSIONS RCP-SSP core scenarios; Table 1). Net habitat changes per species are calculated as the difference between suitable habitat area in the future maps and suitable habitat area in the baseline map (based on mapped IAP2 land use for around 2010 and climate 1991-2012). All net habitat changes are expressed as percentage of the baseline habitat extent to provide a comparable measure across species [Colour figure can be viewed at wileyonlinelibrary.com] the 'EU-DEM' for topography (EEA, 2017). Details of the modeling process are described in Polce et al. (2018). To establish a consistent baseline for the extrapolations that were based on IAP2 land-use projections, we adjusted the baseline distribution maps using the IAP2-simulated baseline land-use areas (approximating the year 2010; Figure 1; Table 2). This simulated baseline is used in the IAP2 as a coherent outcome of varied input data, and used here as a consistent starting point from which scenario outcomes develop; a consistency that would not be possible with other (observational) baseline land-cover data. The same methodology as for retrieving the extrapolations was applied (details in the following sections), including the preprocessing of IAP2 outputs, the mapping between CORINE and IAP2 legends, and the reconstruction of the MaxEnt models from the model coefficients. The exchange of land-use predictor variables resulted in some deviations compared to the original distribution maps ( Figure S1).

| Climate and land-use data
The extrapolations of the SDMs were based on two major data sources (Tables 2 and 3). Land-cover variables (i.e., fractions of CORINE classes; Table 2) were obtained for a range of scenarios (Table 1)

SSP1
Europe of sustainable development, characterized by political stability, low level of inequality, resource-saving lifestyles, and a successful transition toward renewable energy. Moderate economic growth and CO 2 neutrality by 2050.

SSP3
Highly unequal Europe characterized by political instability, resource-intensive lifestyles, severe ecosystem failures, and substantial increases in energy and food prices. Low economic growth and high carbon intensity.
Europe with increasing concentration of power in a small business and political elite. Characterized by growing inequalities and shift toward green technologies, driven by economic crisis and climate change impacts. High economic growth.

SSP5
Europe of rapid technological progress and economic growth that is almost fully based on further exploitation of fossil fuels due to a lack of environmental concern. Very energy-intensive lifestyle and severe environmental degradation.
TA B L E 2 Overview of land-related environmental predictors to model the distribution of bumblebees and mapping between CORINE and IAP2 legends. 'Original data source' refers to data used in Polce et al. (2018). 'Baseline data source' refers to data used to derive baseline bumblebee distributions consistent with IAP2 data 2080 (2071-2100) by mapping the fractions of IAP2 land-use output to CORINE classes as indicated in Table 2 and described in the following section. The IAP2 is an interactive web-based platform that integrates a series of interlinked meta-models representing urban development, water resources, flooding, forests and agriculture, and biodiversity to assess climate change impacts, vulnerability, and adaptation (Harrison et al., 2016. The IAP has been applied and evaluated in a large number of studies including sensitivity and uncertainty analyses (e.g., Brown et al., 2015;Fronzek et al., 2019;Harrison et al., 2015Harrison et al., , 2016Harrison et al., , 2019Holman et al., 2017;Kebede et al., 2015). A full model description and the online model itself are available at http://www.impre ssion s-proje ct.eu/show/IAP2_14855. Land use is modeled on a 10 arcmin grid across Europe and the allocation is based on both biophysical conditions (e.g., soil type, climate suitability) and socioeconomic aspects (e.g., prices, support rules, costs), resulting in multiple land uses proportionally for each cell (Holman & Harrison, 2012). We used the core scenarios of the IMPRESSIONS project (Kok et al., 2015; Table 2 the European extent is 10 arcmin (WGS84; EPSG:4326). To obtain land uses for the 10 × 10 km grid, we overlaid polygons projected to WGS84 with the IAP2 maps and aggregated the IAP2 land uses weighted by the 10 arcmin cell contribution to these polygons. In this way, we aimed to minimize distortions in the spatial configuration of land-use fractions from the IAP2 in the 10 × 10 km grid.

| Mapping between IAP2 land use and CORINE land cover
CORINE land-cover classes as listed in Table 2 have been used as predictor variables in the SDMs by Polce et al. (2018). As land-cover projections based on the CORINE legend do not exist, we use projections from the IAP2 and map these onto the CORINE legend ( Table 2). The IAP2 integrates a set of meta-models to represent the land-use sector (Holman & Harrison, 2012)  areas'). Such classes represent a mixture of land cover/use, which is impossible to map from the IAP2 classes without additional information on the spatial configuration within grid cells. In consequence, the quality of the mapping differs across categories ( Figure S2). However, all our analysis is based on the relative changes from the mapped CORINE land cover at around 2010 to the mapped CORINE land cover in the projections and thus we present a consistent set of modeling experiments.

| Reconstruction of SDMs and modeling experiments
The

| Assessment of land-use and climate impacts
To assess the impacts of land-use and climate change on the distribution of bumblebees, we calculated 'net habitat changes' from the difference of gains and losses of suitable habitat area in 2050 and 2080 compared to the baseline. We assume, following Polce et al.  (Table S1), net habitat changes are expressed as percentage of the baseline habitat extent to provide a comparable measure across species.  it is the rare species that lose larger percentages of their baseline habitat extent.

| Land-use impact on bumblebee habitat suitability
The  (Table S2). Rare species tend to benefit from land-use changes following the SSP3, SSP4, and SSP5 storylines, while SSP1, which is dominated by agricultural intensification and loss of forest, induces severe habitat losses (up to −75% in 2080). In contrast, common species show habitat increases mostly in scenarios associated with SSP1 land-use changes and habitat losses in SSP3-, SSP4-, and SSP5related scenarios, indicating a higher adaptive capacity to agriculturally dominated landscapes.

| Comparison of isolated land-use and climate effects
In Figure  scenarios, but habitat gains in all land-use scenarios but LU-RCP2.6-SSP1. In consequence, the effect of including dynamic land use in the extrapolations may amplify, attenuate, or offset changes to suitable habitat area expected from climate impacts.

| Land-use modulation within climate scenarios
Given a certain level of climate change (RCP2.6, RCP4.5, or RCP8.5), including dynamic land-use change (COMB modeling experiment) has varying impacts on the species' response across land-use scenarios (Table 5). In a world consistent with RCP2.6 climate in 2080, land use following the SSP1 storyline results in an average additional habitat change of −7.50%, which is largely driven by large losses of individual species (Table S4). In contrast, SSP4 land-use changes

F I G U R E 4
Relation between the maximum change in habitat extent in the CLIM scenarios (n = 3) and LU scenarios (n = 7) in 2080. The maximum change in habitat extent (can be both net gain and loss) is considered as a measure of the potential magnitude of isolated climate or land-use impact across the scenarios, respectively. Each dot represents a bumblebee species. The color gradient is based on the habitat extent in the baseline maps (~2010), indicating the rareness/commonness of a species (Table S1) In the COMB-RCP2.6-SSP1 scenario 22 species lose additional habitat compared to the respective CLIM-RCP2.6 scenario, while 25 species gain habitat. In COMB-RCP2.6-SSP4 19 species lose habitat, while 28 gain habitat; but, more importantly, it is different species that are affected in the one or other direction. SSP1 land use leads to strong additional habitat losses for rare species (−31.21% for the lower quartile) while common species tend to benefit (+2.35% and +1.15% for the third and upper quartiles). SSP4 land use, in contrast, attenuates habitat losses for rare species (+11.37% for the lower quartile) and common species experience additional losses (−1.66% for third and −0.83% for the upper quartile). For 35 out of 47 species, the choice of land-use scenario has an opposite effect on the net change to suitable habitat area (i.e., habitat loss in one land-use scenario contrasts with habitat gain in the other) and for 24 species it has the potential to change the sign of the combined signal compared to the CLIM modeling experiment ( Figure 5; Table S2).
Differences between land-use scenarios do also exist in RCP4.5 and RCP8.5 worlds in 2080, although the importance of land-use decreases with higher levels of climate change. The average modulating effect of including dynamic land use is −3.43%, −2.54%, and −0.69% for SSP1, SSP3, and SSP4 in a RCP4.5 climate, and −4.72% and −1.69% for SSP3 and SSP5 in a RCP8.5 climate. For all landuse scenarios, the majority of species experiences additional habitat losses, with COMB-RCP4.5-SSP4 and COMB-RCP8.5-SSP5 representing the scenarios with the largest number of species still likely to gain from land use (20 and 22, respectively). Rare species benefit most from SSP4 in a RCP4.5 climate, and SSP5 land use tends to be more beneficial than SSP3 land use in a RCP8.5 world (Table 5).
However, land-use impacts and interactions between land use and climate are-in contrast to RCP2.6-more important than isolated climate impacts for only a few species (Tables S3 and S4).

| DISCUSS ION
In this study, we extrapolate the potential distribution of 47 European bumblebee species to 2050 and 2080, considering various levels of climate and land-use change, consistent with the RCP-SSP scenario framework. By contrasting projections with dynamic climate or land use (while keeping the other constant at baseline conditions), we isolate changes to the potential future habitat suitability due to climate and land use, at the species level. An additional modeling experiment with both dynamic climate and dynamic land use allows us to assess the modulating effect of including land use in the SDMs within a given climate. To our knowledge, this is the first implementation of the RCP-SSP scenario framework in bumblebee distribution modeling at the European scale and, hence, it is the first assessment of species dynamics across an established, representative range of potential futures. While we use bumblebees as a target species here, our methodology is also applicable to other taxa.   (Escalante et al., 2013;Marshall, Beckers, et al., 2020). our study due to differences among species and environmental conditions, but North American and European species have been previously shown to respond similarly to climate change (Kerr et al., 2015;Soroye et al., 2020).
Isolated land-use impacts on bumblebee habitat suitability (=LU modeling experiment) can reach similar magnitude than isolated climate impacts (=CLIM modeling experiment; Figure 4), although overall changes in habitat extent due to climate change are larger.
Relatively small land-use effects in large-scale assessments and at coarse spatial resolutions have been reported before (e.g., Kerr et al., 2015;Luoto et al., 2007;Marshall et al., 2018). However, the small effects might be to some extent an artifact of the coarse thematic resolution of land-use categories in current assessments that are not able to sufficiently characterize bumblebee habitat (Marshall, Beckers, et al., 2020), as well as of limited land-use change scenarios that do not explore a representative range of future conditions.
Nonetheless, we find a large variability of responses to land-use changes across individual species (from up to >65% gain to <−75% loss), reflecting the variation in habitat requirements and ecological traits across bumblebee species (Cariveau & Winfree, 2015;Winfree et al., 2011). Different land-use storylines result in heterogeneous habitat change responses across species, with SSP1-and SSP4related scenarios revealing the most contrasting patterns between rare and common species (see Figure 3). While the SSP1 storyline results in widespread expansion of arable land and pasture, both are contracting in the SSP4 scenario. The rare species, which tend to be habitat specialists, may be disproportionally threatened by this strong expansion of agricultural land in SSP1 scenarios (Cariveau & Winfree, 2015;Persson et al., 2015;Rasmont et al., 2015), while common species may better adapt or even benefit from additional resources in agricultural landscapes (Westphal et al., 2003). On the other hand, SSP4 scenarios seem to be beneficial for rare species due to higher availability of semi-natural land while common species might be not able to further expand their existing habitat due to abandoned agricultural land turning to forests in these scenarios, which is not considered a preferred bumblebee habitat (Winfree et al., 2011). While the SSPs represent narrative storylines rather than predictive pathways (O'Neill et al., 2017), previous uncertainty analysis suggests that, in the context of the IAP model used here, SSPs 3 and 4 may produce land-use outcomes reasonably central within overall scenario uncertainty . This does not necessarily make these outcomes more likely, but suggests that they could be representative of a range of scenario conditions.

| Game changer land use in a RCP2.6 world?
The COMB simulation (=dynamic climate and land use) suggests that the choice of the land-use scenario becomes particularly important under low-level climate impacts, that is, the RCP2.6 scenario. A majority of the species show substantial discrepancies in response to the two land-use scenarios in a RCP2.6 world, with some species losing >50% of their climatically suitable habitat under one scenario (see Figure 5). In consequence, some of the land-use change impacts on bumblebees seem currently obscured by a focus on climate change and high-end climate scenarios in existing assessments (Titeux et al., 2016  , spatial mismatches are unavoidable ( Figure S2) and evaluation of the IAP2 land-use projections based on more recent CORINE maps is not possible. We approach these issues by proportionally mapping broad IAP2 land-use classes to the spatial pattern of the CORINE baseline map and assessing changes in comparison to this mapped baseline only. While this approach provides an internally consistent assessment of land-use impact on bumblebee habitat suitability, we advise against using our results as absolute effects, for example, to derive recommendations for conservation management. Instead, the results emphasize the heterogeneous response of individual species to a set of combined climate and land-use changes that requires further attention in future research.
The CORINE predictor variables feeding into our SDMs include little information about land-use intensity (e.g., crop diversification or fertilizer and pesticide application) and landscape configuration (e.g., the availability of floral resources or semi-natural habitat patches within agricultural landscapes). Similarly, the representation of bumblebee habitat with broad land-use/land-cover classes may be overly simplistic at a spatial resolution of 10 × 10 km. Both simplifications have been shown to be important for bumblebee population dynamics and distribution (e.g., Luoto et al., 2007;Marshall, Beckers, et al., 2020). In consequence, a more detailed description of bumblebee habitat would affect some of the results of our analysis. For example, scenarios following the SSP1 storyline lead to a strong increase in agricultural land which translates to severe losses of suitable habitat for several species, while scenarios following the SSP4 storyline are characterized by the abandonment of agricultural land leading to smaller losses or habitat gains Kok et al., 2019). However, these results may be biased to some extent as in a more sustainable world (=SSP1) agricultural systems may substantially differ from the ones in SSP4 in terms of land-use intensity (Kok et al., 2019;Mitter et al., 2020), resulting in different habitat qualities provided by agricultural landscapes (Persson et al., 2015).
Including such effects, however, would require a more mechanistic representation of key ecological processes, which, in turn, rely on far higher-resolution data and modeling than currently available to overcome the limitations of SDM approaches (Singer et al., 2016). A key element for a better understanding of the distribution of bumblebee species as a function of land use and climate is set by the EU pollinators' initiative (European Commission, 2018). Action 1 of this initiative, a European pollinator monitoring scheme is proposed (Potts et al., 2021) based on more than 2000 sites where insects, including bumblebee species, will be systematically monitored. It will provide a valuable source of data to test the assumptions and predictions made in this study. Nevertheless, our analysis provides insights into the modulating effects of land use for bumblebee habitat suitability, especially in a world with limited climate change, where a wise management of land can become a game changer in the conservation of key pollinator species.

ACK N OWLED G EM ENTS
The research in this paper was financially supported by the Baden-Württemberg Foundation ELITE program and the Helmholtz Association.
PW is supported by the Alexander von Humboldt Foundation. Open Access funding enabled and organized by Projekt DEAL.

D I SCL A I M ER
The views expressed in the article are personal and do not necessarily reflect an official position of the European Commission.

DATA AVA I L A B I L I T Y S TAT E M E N T
Land-use and bumblebee projections are made available at https:// doi.org/10.5281/zenodo.4557528.