The interplay of climate and land use change affects the distribution of EU bumblebees

Bumblebees in Europe have been in steady decline since the 1900s. This decline is expected to continue with climate change as the main driver. However, at the local scale, land use and land cover (LULC) change strongly affects the occurrence of bumblebees. At present, LULC change is rarely included in models of future distributions of species. This study's objective is to compare the roles of dynamic LULC change and climate change on the projected distribution patterns of 48 European bumblebee species for three change scenarios until 2100 at the scales of Europe, and Belgium, Netherlands and Luxembourg (BENELUX). We compared three types of models: (1) only climate covariates, (2) climate and static LULC covariates and (3) climate and dynamic LULC covariates. The climate and LULC change scenarios used in the models include, extreme growth applied strategy (GRAS), business as might be usual and sustainable European development goals. We analysed model performance, range gain/loss and the shift in range limits for all bumblebees. Overall, model performance improved with the introduction of LULC covariates. Dynamic models projected less range loss and gain than climate‐only projections, and greater range loss and gain than static models. Overall, there is considerable variation in species responses and effects were most pronounced at the BENELUX scale. The majority of species were predicted to lose considerable range, particularly under the extreme growth scenario (GRAS; overall mean: 64% ± 34). Model simulations project a number of local extinctions and considerable range loss at the BENELUX scale (overall mean: 56% ± 39). Therefore, we recommend species‐specific modelling to understand how LULC and climate interact in future modelling. The efficacy of dynamic LULC change should improve with higher thematic and spatial resolution. Nevertheless, current broad scale representations of change in major land use classes impact modelled future distribution patterns.

performance, range gain/loss and the shift in range limits for all bumblebees. Overall, model performance improved with the introduction of LULC covariates. Dynamic models projected less range loss and gain than climate-only projections, and greater range loss and gain than static models. Overall, there is considerable variation in species responses and effects were most pronounced at the BENELUX scale. The majority of species were predicted to lose considerable range, particularly under the extreme growth scenario (GRAS; overall mean: 64% AE 34). Model simulations project a number of local extinctions and considerable range loss at the BENELUX scale (overall mean: 56% AE 39). Therefore, we recommend species-specific modelling to understand how LULC and climate interact in future modelling. The efficacy of dynamic LULC change should improve with higher thematic and spatial resolution.
Nevertheless, current broad scale representations of change in major land use classes impact modelled future distribution patterns.

K E Y W O R D S
biodiversity loss, dynamic, future, land use change scenarios, pollinators, projections, species distribution models (SDMs), wild bees

| INTRODUCTION
Recent scientific consensus suggests that we are facing a sixth mass extinction event, correlated strongly to anthropogenic factors (Ceballos et al., 2015). To avoid the dramatic loss of biodiversity and associated ecosystem services, immediate and thorough conservation efforts are required (Barnosky et al., 2011). An important role of biodiversity conservation research is to understand and estimate potential changes in biodiversity alongside changing abiotic and biotic conditions (Elith, Kearney, & Phillips, 2010;Porfirio et al., 2014).
In an effort to understand these effects experts have produced scenarios of climate, and land use and land cover (LULC) change.
Land use and land cover change scenarios use potential climate change, policy decisions and strategies to represent socioeconomic developments which will inevitably shift land use and management (Rounsevell, Ewert, Reginster, Leemans & Carter, 2005;Van Vuuren et al., 2011;Verburg, Rounsevell & Veldkamp, 2006). Scientists have developed scenarios with the goal to evaluate the impact of environmental changes on biodiversity (Spangenberg et al., 2012). Their role in biodiversity analyses is to allow the production of dynamic land use variables which better reflect future habitat suitability for a species and may be useful to explain additional drivers of distributional changes alongside climate change. There is strong consensus that both climate and LULC change are important in driving the observed patterns of biodiversity declines (Luoto, Virkkala & Heikkinen, 2007;Ostberg, Schaphoff, Lucht & Gerten, 2015). Historically, LULC change has been the dominant cause of observed biodiversity changes and researchers expect that it will remain an ongoing threat to worldwide biodiversity (Millennium Ecosystem Assessment, 2005, Ostberg et al., 2015. Climate and land use change underlie a multitude of environmental pressures that may have a greater joint impact on biodiversity than when operating in isolation (Clavero, Villero & Brotons, 2011;Mantyka-Pringle, Martin & Rhodes, 2012). Therefore, models which exclude LULC change from modelling biodiversity in the future neglect a significant factor in potential drivers of species distribution change, even if these projections are coarse and at broad spatial scales.
Species distribution models (SDMs) represent a powerful tool for understanding patterns in biodiversity. They combine species occurrence data with environmental conditions to estimate the distribution of species in space and time . Often used to project species distributions into unsampled areas, or areas of possible invasion, they also project species distributions into the future (Franklin, 2010). The majority of future distribution models include only climate change variables and do not include LULC variables or use only LULC variables based on current conditions (static; Bellard, Bertelsmeier, Leadley, Thuiller & Courchamp, 2012;Titeux et al., 2016). At broad spatial scales, climate is expected to be the main constraint to species distributions, but at finer resolutions, the effect of LULC covariates increase; landscape-specific features that provide nesting and feeding resources occur at this finer scale (Luoto et al., 2007;Rahbek et al., 2007;Thuiller, Ara ujo & Lavorel, 2004). Therefore, improved estimations of biodiversity change require detailed land use change scenarios (Titeux et al., 2016).
Even though studies recommend the inclusion of LULC variables to avoid producing unrealistic projections, few studies have used dynamic LULC covariates to model biodiversity patterns in the future. Reasons for this is that projections of LULC change are rarely available or only at coarse resolution and with few land use classes (Titeux et al., 2016). However, climate predictions offer similar limitations with resolution and parameters often not directly relevant to the habitat suitability of species. Interestingly, the studies that explicitly include dynamic LULC variables in the SDM process show considerable variation in the effect this has on species distribution patterns, specifically range change (Barbet-Massin, Thuiller & Jiguet, 2012;Chytr y et al., 2012;Ficetola et al., 2010;Martin, Van Dyck, Dendoncker & Titeux, 2013;Riordan & Rundel, 2014;Sohl, 2014;Wisz et al., 2008). The variation is most likely due to differences in species, spatial scale and explanatory variables included in these studies. Likewise, the performance of SDMs usually depends strongly on the modelling framework used, the species modelled, the distribution, quality and quantity of collection data, and the resolution of the species occurrence data and covariates (Aguirre-Gutierrez et al., 2013;Bellard et al., 2012;Harris et al., 2013;Warren & Seifert, 2011). Testing the effect of dynamic LULC covariates with multiple species, different resolutions and covariates is essential to understand their role in SDMs (Martin et al., 2013).  Spangenberg et al., 2012). We expect to observe differences in the projected distributions produced by climate-only models vs. models which include LULC. We expect that the differences between static and dynamic LULC models will be less pronounced and species-specific, and will likely depend on the spatial scale and resolution at which the LULC covariates are projected (Luoto et al., 2007;Martin et al., 2013). Overall, we aim to illustrate the bias associated with using climate change-only scenarios when modelling bumblebees that land use change will undoubtedly affect. We also aim to show how presently available dynamic LULC projections affect the modelled distributions for multiple species. Following this important step, we discuss the extent to which our results provide improvements to land use change scenarios in development and the conservation implications of using such SDMs.

| Target species
Our study group is the genus Bombus, for which we have detailed, long-term, biogeographical records for most of Europe, and which has shown significant decline in the last one hundred years (Biesmeijer et al., 2006;Carvalheiro et al., 2013;Kerr et al., 2015;Rasmont et al., 2005). Forty-eight European bumblebee species were included in the analysis (see Table S1). The species modelled share similar life histories, but exhibit vastly different ranges and distributions in Europe (Rasmont et al., 2015). According to the IUCN Red List of threatened species, Bombus in Europe includes species of all threat levels (Nieto et al., 2014). Climate change impacts have been modelled for the genus Bombus at the European scale, projecting severe declines and northerly shifts for the majority of the species (Rasmont et al., 2015). However, loss of habitat for feeding and nesting resources has been cited as a major driver of past Bombus decline (Biesmeijer et al., 2006;Carvalheiro et al., 2013;Goulson et al., 2010;Williams & Osborne, 2009). Therefore, climate might not necessarily be the only significant driver of change for this group over the next one hundred years. Furthermore, the distribution patterns of wild bee species are reported to be affected by change in major land use classes, particularly the presence of arable land (Aguirre-Guti errez et al., 2015;Senapathi et al., 2015).

| Species presence data
This study includes bumblebee collection records from 22 European countries and multiple sources including professional and amateur scientists (see Fig. S1). The data were collated as part of the EU FP7 project STEP (Potts et al., 2011), and is aggregated and available to view on the Atlas Hymenoptera webpage (Rasmont & Iserbyt, 2013).
We used records from 1970 until 2000, as these represent the 'current' period of climate data, which we used to train the species distribution models. We had 462,636 records available to use.

| Spatial extent and resolution
The spatial extent was limited to the extent of the ALARM projections of European land use, which in turn limited the species collection records available to use (see Fig. S1). Europe in the context of this study is defined as the European Union without Ireland, Romania, Bulgaria, Canary Islands and Cyprus, and including Norway and Switzerland. We created 5 9 5 km, 10 9 10 km and 20 9 20 km European grids for training the SDMs to project onto the BENELUX (Belgium, Netherlands and Luxembourg) region. We also created a 50 9 50 km European grid for training the SDMs to project onto the original spatial extent of Europe. All map projections use the European terrestrial references system 1989 (ETRS89).

| Climate and Land Use Data
Variables of current climatic conditions were produced from monthly interpolated rainfall and temperature data from 1971 to 2000, at a 10 0 resolution (Fronzek, Carter & Jylh€ a, 2012;Mitchell, Carter, Jones, Hulme & New, 2004). We considered 14 climate variables for the modelling process (see Table S2). However, because climate variables are often strongly correlated. Including all climate variables in the models would have added redundant information. Therefore, to avoid collinearities, we conducted a selection according to Pearson correlation coefficients (<0.7; Dormann et al., 2013). When two variables were highly correlated, we selected the variable that we estimated to have the greatest ecological relevance to Bombus species.
We selected total annual growing degree-days (>5°c), which was correlated with other temperature variables, because it is linked to the presence of wildflowers and flowering crops, both important food sources for bumblebees. Furthermore, we chose water balance, which was correlated with the majority of other precipitation variables because it is representative not only of total precipitation, but has a direct link with temperature, making it an important influence for terrestrial vegetation (Gerten, Schaphoff, Haberlandt, Lucht & Sitch, 2004). Five climate variables were used as explanatory covariates in the model: average precipitation of the wettest month; total annual number of growing degree-days above 5°C; mean diurnal range (mean of monthly difference between daily maximum and minimum temperatures); annual temperature range (maximum temperature of warmest month-minimum temperature of coldest month); and annual water balance (mean monthly precipitation minus the monthly potential evapotranspiration; Gerten et al., 2004).
Each of the five climate variables was aggregated to the 50 9 50 km and 20 9 20 km grids, and downscaled to the 10 9 10 km and 5 9 5 km grids using bilinear interpolation (Randin et al., 2009). All spatial analyses were conducted using Rstatistics • 'Growth Applied Strategy' (GRAS)-IPCC A1FI; mean projected temperature rise in Europe at 2100 is 5.6°C; a maximum change scenario driven by policies of deregulation and economic growth.
scenario; mean projected temperature rise in Europe at 2100 is 3.0°C; a moderate change scenario driven by economic, social and environmental policies, related to stabilizing atmospheric greenhouse gases emissions and stopping the loss of biodiversity.
Current land use was obtained from the Coordination of Information on the Environment (CORINE) Land Cover at 250 9 250 m resolution (Bossard, Feranec & Otahel, 2000). The CORINE classes were reclassified as six classes to match the future projections. We removed the class 'others' from our analysis because it represents diverse land use types and was inexplicable in an ecologically relevant context for bumblebee species. Future land use was obtained from the ALARM EU project downscaled to 250 9 250 m for each of the three scenarios for 2050 and 2100 (Dendoncker, Bogaert & Rounsevell, 2006;Spangenberg et al., 2012). At each grid resolution, we determined the percentage cover for each land use covariate.
The final five land use layers were: percentage cover arable land; percentage cover forest; percentage cover grassland; percentage cover permanent crops; and percentage cover urban.
The role of the covariates will be tested in three ways using three variable sets in the models: (1) Dynamic climate-only models, suggesting that only climate variables matter in the future distribution of bumblebee species.
(2) Static land use and dynamic climate, suggesting that land use variables are important in delimiting species habitat suitability, but that their future change will be driven only by climate change and changes in land use are redundant. (3) Dynamic climate and dynamic land use, suggesting that future distribution patterns will be dependent on the interaction between changing climate and changing land use.

| Species distribution modelling
We used a SDM approach to compare the role of dynamic land use data in the future distribution patterns of bumblebees. We modelled the distribution of 48 species using R (R Core Team, 2012) with the biomod2 package (version 3.3-3; Thuiller, Georges, & Engler, 2015).
We chose an ensemble modelling approach, which creates a consensus of the predictions of multiple algorithms and is an established method to account for projection variability (Thuiller, 2014). Even small differences between algorithms can lead to different projections of future distribution change. Ensemble modelling aims to limit the many uncertainties of forecast modelling and has become increasingly used in studies of biodiversity change (Thuiller, 2014).
We chose three algorithms to include in the ensemble model based on their previous performances with analogous collection data for a similar insect species group (Aguirre- Gutierrez et al., 2013). The three algorithms chosen were a generalized linear model, GLM (Nelder & Wedderburn, 1972), with linear and quadratic effects, and stepwise selection based on the Akaike Information Criteria (AIC); a generalized boosted model, GBM (Friedman, 2001), with 3,000 trees and five cross-validation folds; and Maximum Entropy Modelling Models for each species were trained at multiple resolutions at the European scale; 5 9 5 km, 10 9 10 km, 20 9 20 km and 50 9 50 km. We had 462,636 records available to use; these were aggregated as unique species occurrences for each grid cell resolution. The number of occurrences per resolution is as follows: 67030 at 5 9 5 km, 49146 at 10 9 10 km, 30104 at 20 9 20 km and 21,162 at 50 9 50 km. We modelled 48 species (see Table S1) with at least 50 unique records, and for which there are no ongoing taxonomic debates surrounding their species definition (see Rasmont et al., 2015). A number of occurrences in the database were not point level GPS coordinates, but were recorded as UTM grids of varying sizes. To be confident in the spatial accuracy of collection records we removed occurrences that were recorded as UTM grids larger than 1 9 1 km. As the sampling methods were diverse and nonsystematic, there are likely spatial biases amongst the records.
To deal with this potential spatial autocorrelation between closely sampled locations we selected a subset of points per species. A random starting observation was selected and all points in adjacent grid cells removed; this was then repeated for all remaining points. This produced a more even spread of observations and minimized the effects of heavy sampling at particular locations.
As true absences were not available (it is not possible to accurately say that a bee species is not present during sampling) we generated randomly distributed pseudo-absences for GBM and GLM and selected a background sample for MAXENT (Elith et al., 2011;Phillips et al., 2009). We used target-group sampling to select our background points (Mateo, Croat, Felic ısimo & Muñoz, 2010;Phillips et al., 2009). We specified that the background samples and pseudoabsences could only be selected from areas where other bumblebees have been recorded since 1970. This approach is more objective than taking the background and pseudo-absence samples from sites that have not been sampled, accounting for potential sampling bias (Elith et al., 2011;Phillips et al., 2009) and providing more accurate results (Mateo et al., 2010). To account for within algorithm variation we trained the models 10 times for each of the 48 species, the three algorithms, the three model hypotheses, and the four grid resolutions. This resulted in 360 models per species. We used a bootstrap approach where random subsets of 80% of the data were used for model training and the remaining 20% to produce Area Under the Curve (AUC) values to test model performance (Bahn & Mcgill, 2013;Jim enez-Valverde & Lobo, 2007). For each covariate included in the model, we calculated variable contribution as the change in correlation between the covariates and the response with and without the selected variable (Thuiller et al., 2015). We then produced an ensemble model for each of the three model hypotheses, creating a median representation of the predictions of the 10 runs and three algorithms together. We chose the median value as it is less sensitive to extreme values than the mean.
We projected the models trained at 5 9 5 km, 10 9 10 km and 20 9 20 km, onto BENELUX. BENELUX comprises no novel conditions under the scenarios (i.e., there are no conditions in BENELUX in 2100 that do not already occur within Europe). Therefore, no forecasting into unknown ecological space occurred (Fig. S2). We also projected the data trained at 50 9 50 km onto the entire European study area. For each species we produced habitat suitability maps of the median ensemble predicted distribution. One map was produced for each of the three model types at 2050, and 2100 under the three change scenarios at the 4 grid resolutions. Habitat suitability maps were converted to binary presence absence maps using the values under which specificity and sensitivity is optimized (Thuiller et al., 2015). [SLM]), we created separate mixed effects models for each of the three metrics for both Europe and BENELUX projections. We included species as a random effect, as we were interested in how changes in distribution of the species vary across the different model types, periods and scenarios, and not in the inherent variation between species. Furthermore, to determine if our results are related to the structure of the data we also included the current range of the species as a covariate.

| Statistical analysis
Due to large numbers of zeros both range loss and range gain at the BENELUX scale were analysed with two separate mixed models: Bernouli distributed models of the probability of gain or loss and a linear mixed effects model of values given range loss/gain were projected.
Finally, in addition to presenting results for bumblebees as a group, we chose two species, Bombus argillaceus (Scopoli, 1763; increasing in range) and B. veteranus (Fabricius, 1793; decreasing in range), to look more closely at the difference between model projections with and without LULC covariates. We chose these two species as they are at opposite end of the spectrum of climate risk, both had high model performance values, both have a large number of collection records within Europe and we believe them to be representative of two futures, i.e. considerable range gain and considerable range loss, respectively (Rasmont et al., 2015). The current distribution of B. argillaceus is in Southern and South Eastern Europe as well as Western Asia (Rasmont & Iserbyt, 2013). In previous climate-only models of future conditions B.
argillaceus was projected to increase its range considerably in Western Europe (Rasmont et al., 2015). Bombus veteranus exhibits an already patchy distribution in the plains of Northern Europe and has already declined in Belgium, shifting from an abundant species to one which is barely present (Rasmont & Iserbyt, 2013). Under future climate-only projections B. veteranus is expected to decrease in range considerably (Rasmont et al., 2015).

| Model training fit and variable contribution
For models trained on the current period, we assessed model fit using AUC scores. An AUC value below 0.5 indicates a model fit that is not better than random, values above indicate enhanced model fit.
We used AUC values to compare the change in model fit per species with LULC vs. a COM (Figure 1) We also compared the variable contributions of the different explanatory covariates included in the models (Figure 2). Climatic variables are the most important in explaining the current distribution of all species. The total annual number of growing degree-days was included amongst the four most important variables for 44 of the species modelled. The most important LULC covariate is the percentage cover of arable land but the percentage cover of forest is also important for a number of species (Figure 2). Overall LULC variables contribute 15% of total variable importance.

| The future of bumblebees projected at the BENELUX scale
Of the distribution change metrics analysed, the largest discrepancies were found in the projected range loss (Figure 3a,b). There is considerable variability between species and between scenarios but model type has a significant effect on the projections of whether species will lose range and how much range will be lost (Table 1). Overall species are more likely to lose range under DLMs than both COMs and SLMs (p < .001 and .002; Table 2). However, given range loss occurs (i.e. excluding species that showed no range loss) then greater loss is projected by COMs than both SLMs and DLMs (1.3%; p < .001; Table 2). However, this relationship is highly variable and species specific. Under COMs 11 species show greater mean range loss averaged across scenario and resolution, however, five species show greater range loss under DLMs (Figure 3a). The relationship between projected range loss of SLMs and DLMs, while not significant at the BENELUX scale, (Table 2) also appears to be species specific, with some species below the equal projection line, indicating greater range loss under DLMs (Figure 3b). There are no significant interactions between model type and other explanatory variables, suggesting a consistent effect of model type across scenarios, periods and resolutions (Table 1).
Model type, period, scenario and resolution at which the modelling occurred significantly influence the probability of range gain (Table 1). Only 50% of species were projected to gain any range at all within BENELUX by 2100 (Figure 4a,b). The odds of range gain are significantly higher for DLM projections than for COM and SLM (p < .0001; Table 2). When range gain occurs there is no significant difference between COMs and DLMs, however, both projected significantly higher loss than SLMs (1.4 and 1.2%, p < .0001 & .03; Table 2). This can be visualized in Figure 4a, where variation F I G U R E 1 Area under the curve (AUC) statistics for median-ensemble-model performance visualized per species. Black squares represent models with only climate covariates and grey triangles models with land use land cover (LULC) covariates and climate covariates. Groupings represent Climatic risk as calculated by the Climate Risk Atlas for Bumblebees (Rasmont et al., 2015). Potential risk ( Period and scenario at which the modelling occurred significantly influence the directional shift of the distribution centroid (p < .001; Table 1 and Figure 1). Model type did not significantly affect the projected shift.

European scale
At the European scale with lower spatial resolution (50 9 50 km), SLMs project significantly less range loss than both COMs and DLMs (2.9% and 1.7%; p = <.001 and .02, Table 2). Overall, all 48 species are projected to lose at least some range and the relationships between the different model types shows a strong linear correlation, but with considerable deviation from the assumption of the projections being equal (Figure 5a,b). Eighteen species are projected to lose greater range under COMs whilst fourteen species are projected to lose greater range under DLMs (Figure 5a). The relationship between DLMs and SLMs is clearer with a higher number of species below the equal protection line than above, which supports the significant effect found in the mixed models (1.21%, p < .01; Figure 5b and Table 1).
At the European scale greater range gain is projected by COMs than SLMs and DLMs (2% and 1.6%; p < .001; Table 2). DLMs project greater range gain than SLMs (1.2%, p value = .01; Table 2). This relationship is visible in Figure 6a with the majority of species considerably above the equal projection line. The same pattern is observed for SLMs and DLMs, with 12 species below the equal projection line. The majority of species only illustrate modest range gain, and the differences between model types are emphasized when range gain is high (Figure 6a Finally, the current size of the distribution was also included in some best models, at the European scale more widespread species lose less and gain more range (for full details of all models see Table S3 and Figs. S3-S10).

Range size present 9
Resolution p-values: .01 ≤ p ≤ .05 = *, .001 ≤ p ≤ 0.01 = ** and <.001 = *** The most parsimonious models as chosen by Bayesian information criteria (BIC) for the percentage range loss, percentage range gain, and shift in the distributional centroid for 48 bumblebee species at European and BENELUX scales. The significance of each term included in the model is shown. The symbol "-" represents a variable not included in the best model. The random term for all models was '1 | species.' For a detailed version of the table see Supplementary Table S3. of Europe (Figure 9). A large amount of the projected range loss is in areas with novel climatic conditions, making the predictions unreliable.
Bombus veteranus is one of the numerous European bumblebee species projected to lose a large part of its suitable habitat under climate change; it is therefore representative of the majority of p-values: .01 ≤ p ≤ .05 = *, .001 ≤ p ≤ .01 = ** and <.001 = *** Showing the fixed effect and the significance of the best models as chosen by Bayesian information criteria BIC. Null hypothesis tested: that the difference between contrasts is equal to 0. Values are averaged over other explanatory variables included in the model (see Table S1.)

| DISCUSSION
This study shows that incorporating dynamic LULC change scenarios, even those with low precision and few classes, can have significant effects on the projected distributions of bumblebee species. Species can only occur in a location at any time when a series of critical conditions are met, including both suitable climate and land use and land cover types that allow them to feed, grow, survive and reproduce.
Therefore, it is surprising that the use of climate change projections is commonplace, whereas LULC change projections are rarely used in species forecasting (Titeux et al., 2016). We tested the effect of dynamic LULC variables on projecting future distribution changes for showing that using only climate covariates may over-represent the species range in the present (Luoto et al., 2007;Sohl, 2014;Stanton, Pearson, Horning, Ersts & Res ßit Akc ßakaya, 2012). This is likely to misrepresent species range change as well as the shift of species range limits. The importance of LULC change is dependent on whether habitat requirements, namely nesting and feeding resources (Busch, 2006), can be adequately captured by the relationship between these six land use covariates and the climate change covariates. Therefore, we saw variation for bumblebees as they differ significantly in their landscape requirements (Goulson et al., 2010;Persson, Rundl€ of, Clough & Smith, 2015 (Kerr et al., 2015).

| Models including dynamic LULC compared to static LULC models
Including static LULC change in SDMs is based on the incorrect assumption that LULC will not change in the future or that this change is negligible in comparison to climate change (Stanton et al., 2012). In this study, loss and gain of suitable habitat was more likely with dynamic LULC covariates. The distribution patterns of DLMs represent more variable suitable habitat conditions in time than SLMs under equivalent climate change, resulting in greater projected range loss and gain. However, this pattern varied between species and was more discernable for some over others. This variability is supported by other studies; including dynamic LULC covariates previously led to more accurate model predictions for invasive bullfrogs (Ficetola et al., 2010) and central European plants (Chytr y et al., 2012), but not so for a European butterfly species (Martin et al., 2013 4.4 | LULC-inclusive models for forecasting and guiding conservation efforts The importance of including LULC projections depends on the goals and desired outcomes of the modelling process. As a tool, SDMs explore unknown areas and periods where conditions may be suitable for species occurrence, observe the role of environmental covariates and influence conservation management (Franklin, 2010).
However, due to limitations in data availability and modelling methods their value to conservation and ability to predict occurrence should not be overestimated (Lobo, 2016) New scenarios should emphasize a relevance to biodiversity and land use management, for example, separating between natural-grassland and agricultural-grassland, and intensive and less intensive farming systems. While the incidence of and change in forest and arable land cover correlates with habitat suitability, this is an indirect effect. The coarseness of these classifications does not provide an adequate foundation to extract causal information or infer on the importance of land use management (Thuiller et al., 2004). Moreover, national and international policies, such as the CAP in Europe, tend not to change land cover per se (grassland remains grassland), but the management level and thus biodiversity value. For example, arable land cover is the most important LULC covariate for the majority of bumblebees as defined by the correlative variable importance values (see Table S3). However, the ecological significance of this importance could relate to agricultural intensification, pesticide use, availability of floral resources, or most likely, a combination of these factors.

| Differences between the data sources
Among the 48 bumblebees modelled there are examples of polytypic species representing significant intraspecific variation (Rasmont, 1983). For example, SDMs at subspecies level for B. terrestris performed differently from aggregated models with all subspecies as a single unit (Lecocq, Rasmont, Harpke & Schweiger, 2016). We did not utilize this variation; we modelled the habitat requirements of each species using all available records. Occurrence points were selected to represent the highest resolution possible to model at 5 9 5 km resolution, and many points were removed. However, due to the nature of the data and the multitude of sources it is still likely that some point records were not accurately recorded, though we expect this number to be minimal (Duputi e, Zimmermann & Chuine, 2014).
There were distinctions between the resolution of the climate and land use sources in the past and in the future. Due to the coarse nature of Atmosphere-Ocean General Circulation Models (AOGCMs) used to calculate the climate-change covariates they do not represent accurately fine scale effects (Fronzek et al., 2012). This means at the 10 9 10 and 5 9 5 km resolutions that fine-scale topographic effects of climate may be lost. This may result in a more homogenous representation of climate at these resolutions, which may overrepresent range size and connectivity. However, this is representative of climate data often used in studies of this type to model in the future, and in general climate is more homogenous than land use, particularly at the BENELUX scale. To understand in detail the climate effects on biodiversity, fine scale climate change projections are required. The land-use change maps were downscaled to match the availability of current LULC data at European scale. However, the downscaling algorithm is likely to produce some clustering for the future LULC covariates . Therefore, we focused on percentage cover covariates and it was not possible to include covariates of connectivity and edge effects, as they would not be comparable to present conditions. Furthermore, the land-use change models were created in congruence with climate variables; this means that present and future comparisons are valid at the different modelled resolutions .
Finally, there are many methods for SDM and changes to the modelling algorithms, covariates and resolutions can affect the results (Aguirre- Gutierrez et al., 2013;Warren & Seifert, 2011). We chose to use simplified algorithms in an ensemble modelling approach to account for this variation (Thuiller, 2014).

| CONCLUDIN G REMARKS
This work represents a detailed analysis of the effect of dynamic LULC scenarios at different scales on the projected distributions of multiple species. We show species dependent responses to the effect of dynamic LULC, which demonstrates that these types of scenarios can play a significant role in projecting species distributions under climate change. Climate variables alone, whilst driving habitat suitability, are unlikely to project accurately the detailed distribution patterns of all species. Therefore, we advocate repeated use and testing of these available scenarios with multiple species. However, new scenarios and projections of LULC change at finer spatial and thematic resolutions that indicate management practices will be necessary to better assess biodiversity change in an uncertain future.

ACKNOWLEDG EMENTS
We acknowledge funding from the BELSPO (BELSPO; BR/132/A1/ BELBEES). We specifically thank Leopoldo Castro, Centre Suisse de Cartographie de la Faune (CSCF/SZKF), Bj€ orn Cederberg, Gilles Mah e, Aulo Manino, for providing data and we thank all others who provided their data in the context of the STEP Project (www.step-project.net).
We also thank Steve Donovan for his help in proofreading the manuscript. Finally, we thank the two anonymous reviewers who provided