Predicting bee community responses to land-use changes: Effects of geographic and taxonomic biases

Land-use change and intensification threaten bee populations worldwide, imperilling pollination services. Global models are needed to better characterise, project, and mitigate bees' responses to these human impacts. The available data are, however, geographically and taxonomically unrepresentative; most data are from North America and Western Europe, overrepresenting bumblebees and raising concerns that model results may not be generalizable to other regions and taxa. To assess whether the geographic and taxonomic biases of data could undermine effectiveness of models for conservation policy, we have collated from the published literature a global dataset of bee diversity at sites facing land-use change and intensification, and assess whether bee responses to these pressures vary across 11 regions (Western, Northern, Eastern and Southern Europe; North, Central and South America; Australia and New Zealand; South East Asia; Middle and Southern Africa) and between bumblebees and other bees. Our analyses highlight strong regionally-based responses of total abundance, species richness and Simpson's diversity to land use, caused by variation in the sensitivity of species and potentially in the nature of threats. These results suggest that global extrapolation of models based on geographically and taxonomically restricted data may underestimate the true uncertainty, increasing the risk of ecological surprises.

impacts 9 . There has therefore been much research into responses of bee communities to human impacts such as land-use change and intensification.
A number of syntheses have attempted to identify general trends in the response of bees to human impacts 5,10 . However, their datasets have often been geographically limited, with the majority of data arising from North America and Western Europe 11 . The geographic patterns of bee decline and diversity are not understood sufficiently well to ensure that such generalisations are valid 11,12 . If species' responses to disturbance vary among regions, geographically-restricted models will be inadequate to support broad conclusions. The consequences of basing management strategies on extrapolations from such models could be severe, as many under-studied regions have a high economic dependency upon animal-pollinated crops 11,13 and may generally have limited governmental capacity to adapt to environmental changes 14 .
Geographic variation in bee community responses could arise because differences in land-use history and practices mean that the threats facing assemblages differ across regions. Species subject to very recent disturbance may be more vulnerable, whereas extinction filters [15][16][17] may have already removed many susceptible species from landscapes where the intensification of farming started already decades ago, such as in temperate European agricultural landscapes. Extinction debt may make matters worse still, if the full impact of land-use changes is not yet evident 18,19 . In addition, differences in landscape context across regions can influence species' responses. For instance, Winfree et al. 5 found that habitat loss and fragmentation significantly affected bee communities, but only in areas where little natural habitat still remained.
Bee community responses may also vary regionally because community composition varies geographically. Taxa can differ in their intrinsic susceptibility to land-use change and intensification, through having different functional response traits [20][21][22] , the distribution of which within a community can affect resilience to pressures 23 . A geographic bias towards North America and Western Europe has also resulted in a taxonomic bias; for instance, bumblebees (Apidae: Bombus) are particularly diverse in these areas, whereas large areas of the world have no native bumblebee species (e.g., most of Africa and Australasia). In addition, bumblebees are large, often abundant species with long flight seasons and relatively slow flight, making them fairly easy to sample and, in many cases, to identify. Bumblebees may be more or less sensitive than other bees due to their ecological traits and habitat requirements 24 , which have been shown to influence responses to human impacts and vulnerability to decline 25,26 . In addition, bumblebees have shown clearer declines than other bees in North America 25 and some European countries 27 , so they may be atypical of broader bee diversity.
We compiled a global dataset of bee diversity from published sources of bee assemblages in sites differing in pressures such as land use, and used this to explore whether models of responses to human impacts are robust against geographic and taxonomic biases. Specifically, we hypothesized that bee responses to land-use pressures should vary significantly with region and with taxonomic group (i.e., bumblebees or other bees) and so models and projections will not be transferable across regions and taxa. Improved understanding in this area will help to clarify whether knowledge based on a few regions and taxa is sufficient to underpin policy decisions as well as highlight systems for future study.

Methods
Data Collation. Data were sought from the literature where bee species abundance and/or occurrence were reported for multiple sites. Suitable papers were identified by searching Web of Science at various times from 2011 to 2015, as well as searching journal alerts and assessing references cited in reviews. Papers were further considered if more than one site was sampled for bee diversity using the same sampling method in the same season and geographic coordinates of each site were available. Papers were prioritised if their data were collected from February 2000 onwards, so that biodiversity data could be matched with remote-sensed data from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS). Data were supplemented with sources found through the PREDICTS project (www.predicts.org.uk), which aims to develop global statistical models of how local biodiversity responds to human impacts 28 . The database presented here is not a comprehensive compilation of published sources on occurrence and abundance of bee species across sites differing in land use or intensity, because of regional differences in the ability to retrieve information about potential sources and because most researchers we contacted did not make their data available. The dataset will, however, still be useful for researchers wishing to study land-use impacts on this important taxonomic group.
Where possible we extracted site-level records of bee species (Hymenoptera: Apoidea) occurrence and abundance from suitable papers, along with data for other taxonomic groups if available. Raw data were usually not included within the papers or supplementary files, so the papers' corresponding authors were asked for these data. Relevant data were available from 69 papers, hereafter referred to as 'sources' (Table 1). Each source contains one or more studies, where a study is defined as the set of samples within the same country that were taken using the same methodology. By defining studies in this way, we reduce the impact of broad-scale biogeographic differences in diversity and avoid the confounding effects of methodological differences: within, but not between, studies, diversity data can be compared among sites in a straightforward fashion. Differences in sampling effort within a study were corrected for when necessary by dividing abundance by the sampling effort unit. This assumes a linear relationship between abundance and sampling effort; generalised additive models suggested that this assumption was appropriate (gamm4 package 29 , see Supplementary Data S1 for details). Within each study, we recorded any blocked or split-plot design. The major land-use class and use intensity at each site were assessed based on information in the associated paper, using the scheme described in Hudson et al. 28 (reproduced in  Supplementary Table S1). Briefly, land use was classified as primary vegetation (native vegetation not known to have ever been completely destroyed), secondary vegetation (where the primary vegetation has been completely destroyed; this can include naturally recovering, actively restored, or semi-natural sites), cropland (planted with herbaceous crops), plantation forest (planted with crop trees or shrubs), pasture (regularly or permanently grazed by livestock) or urban (areas with human habitation, where vegetation is predominantly managed for civic or Continued personal amenity). Use intensity was classified according to a three point scale: low, medium and high intensity. For instance, high-intensity cropland would be monocultures with many signs of intensification such as large fields with high levels of external inputs, irrigation and mechanisation; medium intensity cropland would only show some, but not all, features of higher intensity cropland; low-intensity would refer to small fields with mixed crops and little to no external inputs, irrigation or mechanisation. In one data source, information on the use intensity was unavailable at the site-level, so information at the landscape level was used.
The dataset contained 111 studies from 69 sources and 3211 within-study sites (Table 1). This amounted to 195,357 species diversity measurements (i.e., bee taxa and other taxa, Table 1), including 107,176 measurements of bee diversity (a single measurement being, for example, the abundance of a given species at a given site; see Supplementary Data S2 for species list).

Analysis.
For this analysis, we did not include studies that recorded only particular target species (for instance, studies that were only interested in the abundance of a single species across sites), so that site-level diversity measures would be meaningful. The final dataset for the analysis included 101,524 diversity records from 837 bee species at 2421 sites from across the globe (North America: 239 sites; Central America: 103; South America: 176; Western Europe: 1211; Northern Europe: 325; Eastern Europe: 64; Southern Europe: 50; Middle and Southern Africa: 39; South Eastern Asia: 31; Australia and New Zealand: 183). In this reduced dataset, many combinations of land use and use intensity had too few sites to permit robust modelling. The data were therefore aggregated to give a variable of combined Land Use and Intensity (LUI) with the following levels: primary vegetation, secondary vegetation, low-intensity cropland, medium-intensity cropland, high-intensity cropland, pasture, plantation forest and urban. All LUI levels had at least 170 sites, except for plantation forest and urban areas, which were scarce in the dataset with only 105 and 94 sites respectively. Sites were also classified by region and subregion (according to United Nations classifications), with Middle and Southern Africa combined into a single category to increase the sample size.
For each site, we calculated three measures of bee community diversity as our response variables: total abundance, within-sample species richness and Simpson's diversity. Simpson's diversity was calculated as:  Table 1. Data sources and sample sizes. mMLE = largest Maximum Linear Extent (in meters) of any site in the source. MLE is the maximum distance between sampling points within a site, e.g. the length of a transect or the distance between pan traps. nr = not reported. Numbers of taxa are the numbers of unique taxa for which diversity measurements are given (so, if diversity measurements are available only for all bees combined, this would count as one taxon). The percentage of bee species with a known binomial name is also given (% binomial). Note that the figures here represent available data as curated by the PREDICTS team; these will not necessarily match figures in the original papers. + Data were used in the presented analysis. † Data will be incorporated into the PREDICTS database (which will be made openly available). ‡ Data are available from the referenced paper. For all other datasets, please contact the corresponding author of that paper directly.
where P i is the proportion of individuals belonging to species i. We use Simpson's diversity as it stabilises faster than species richness and other diversity measures as specimens accumulate 30 .
As total abundance measurements are not necessarily integers (e.g. densities and effort-corrected measures), use of the Poisson error structure was not possible, so total abundance was ln + 1 transformed before modelling to normalise residuals and equalise variance. Total abundance and Simpson's diversity were modelled using Gaussian errors (model-checking showed that these treatments were appropriate). Species richness was modelled with Poisson error distribution and log-link function; there was evidence of significant overdispersion in these models so an observation-level random effect was included to account for this (i.e., a Poisson-lognormal model) 31 .
All analyses were carried out using R 3.1.0 32 . We constructed models for each response variable, using mixed-effects models (lme4 package 33 ) to account for non-independence of data due to differences in collectors ('source'), sampling methodologies and biogeographic source pools ('study') and the spatial structure of sites ('block'); the initial random-effects structure was therefore block nested within study within source. The initial fixed-effects structure of models included LUI, subregion and their interaction. Subregion is treated as a fixed rather than random effect as we are interested in testing the effect, rather than simply estimating the variance associated with geographic subregion. We test differences in responses to LUI among subregions rather than assessing how responses vary with the latitude and longitude of sites, as subregions represent political differences in land-use patterns and data availability, as well as to some extent reflecting biogeographical differences in community composition.
The best random-effects structure was assessed using likelihood ratio tests 34 , with models fit using Restricted Maximum Likelihood for total abundance and Simpson's diversity, and Maximum Likelihood for species richness. We then attempted to simplify the fixed-effects structure using backwards stepwise model simplification and likelihood ratio tests, with models fit using Maximum Likelihood [34][35][36] . Significance of terms in the minimum adequate models were assessed using Type II Wald Chi Square Tests 37 . However, to better appreciate the uncertainty in the models 38 , if the interaction between LUI and subregion remained in the minimum adequate model, we also constructed the following models: additive model (with LUI and subregion included as additive effects); LUI only (univariate model); and subregion only (univariate). We then compared the explanatory power and predictive error of the interactive model with these simpler alternatives.
Explanatory power was calculated using the MuMIn package in R 39 , as the marginal and conditional R 2 glmm values: i.e., the variance explained by fixed effects alone and by fixed and random effects combined, respectively 40 . Predictive error was calculated as the Mean Squared Error (MSE) from ten-fold cross validation, where the model was iteratively fit to nine-tenths of the data (training set), and validated on the final tenth (validation set); we did this by randomly assigning sites into ten approximately equal-sized groups 41 . As the data are structured, the training data may not be fully independent of the validation data 42 , but any bias in prediction error that this causes will apply equally to all models being compared as the random effect structures are identical. In addition, some combinations of explanatory variables only occur in few studies or sources; splitting the dataset by these higher-level strata would mean that some combinations would rarely appear in the training data, leading to overestimates of predictive error. MSE was decomposed into measures of bias and variance, which give an indication of the accuracy and precision of predictions respectively 43 (See Supplementary Methods for details).
The dataset was then subset to include only studies where both bumblebees (Apidae: Bombus) and other bees were sampled (bumblebees contributed over 19% of the bee abundance records); this resulted in 1636 sites from 47 studies. We calculated the site-level diversity measures separately for each group and fitted the initial model with a three-way interaction between LUI, subregion and taxonomic group (Bombus or otherwise). The initial random structure was as above, but included a site-level random effect to account for multiple samples (bumblebees and other bees) being taken from the same site. As above, we first attempted to simplify the initial model (both in terms of random effects and then fixed effects) and, if the initial three-way interaction remained in the model, compared the explanatory power and predictive error with simpler models, where responses to LUI were permitted to vary with subregion (LUI, subregion and their interaction) or with taxonomic group (LUI, taxonomic group and their interaction).
To further understand heterogeneity in community response to LUI, planned comparisons were performed (multcomp package 44 ). Within each subregion (and each taxonomic group, if assessed), we tested for differences between natural vegetation (primary vegetation) and all other land uses; between semi-natural vegetation (secondary vegetation) and all other land uses (except primary); whether low-intensity cropland differed from medium-intensity cropland; and whether medium-intensity cropland differed from high-intensity cropland. To avoid rank-deficiency, LUI and subregion were collapsed into a single factor in these models. Not all comparisons were possible in all subregions. Multiple comparisons were corrected for using the False Discovery Rate method to adjust significance values 45,46 .
An alpha value of 0.05 was used in all tests for significance. Spatial autocorrelation was assessed in residuals of minimum adequate models using Moran's I, for each study in turn (spdep package 47,48 ). As multiple tests are carried out, we expect 5% of these to be significant by chance so we additionally test whether the proportion of studies showing autocorrelation exceeds this expected proportion (using a one-sided Chi squared test).

Results
For total abundance, Simpson's diversity and species richness, the minimum adequate models were those in which responses to LUI were free to vary among geographic subregions. These models also always had the greatest explanatory power and were always among the models having the lowest predictive error (Fig. 1). Overall, however, explanatory power of fixed effects alone was fairly low, with most variation instead being attributed to random effects, which is not surprising given the variation in sampling methodology and effort among studies and sources.
For total abundance, the interaction between LUI and subregion explained a highly significant amount of variation (χ 2 = 133.15, df = 27, p < 0.0001) (Fig. 2), which resulted in this model having the lowest predictive error compared to simpler models. The interaction between LUI and subregion was also significant for Simpson's diversity (χ 2 = 66.48, df = 27, p < 0.0001) and species richness (χ 2 = 96.41, df = 27, p < 0.05), although the predictive error was not much lower than for models based on subregion alone. In all cases, the interactive model had slightly higher bias than some simpler models, but the lowest variance (See Supplementary Table S2).
In South America, bees were more sensitive to agricultural land uses: medium-intensity cropland maintained significantly lower Simpson's diversity than secondary vegetation (z = − 5.15, p = 0.00029). Urbanization had differing effects between subregions, with increased species richness (z = 5.29, p = 0.00022) in Middle and Southern Africa, but no strong effect detected elsewhere.
When the dataset was split by taxon (Bombus vs. others), the best models for each response variable according to likelihood ratio tests included significant three way interactions between LUI, subregion and taxon (total abundance: χ 2 = 217.9, df = 13, p < 0.0001; Simpson's diversity: χ 2 = 27.62, df = 13, p = 0.0102; species richness: χ 2 = 76.08, df = 13, p < 0.0001). These models also had the greatest explanatory power and lowest predictive error, compared with simpler models (Fig. 3); for total abundance, the lower predictive error was driven by both lower Figure 2. Predicted means of total (logged) abundance of bees for different land-use classes in each subregion, with 95% confidence intervals. Also shown are significant results of multiple comparisons, testing differences between natural (Primary vegetation) and semi-natural land uses (Secondary vegetation) to humandominated land uses, and differences between low, medium and high intensity cropland (*p < 0.05 , * *p < 0.01 , * * * p < 0.001). bias and variance, while for Simpson's diversity and species richness, the slightly higher bias was offset by lower variation (See See Supplementary Table S2).
Planned comparisons provided more detail into these context-and taxon-specific differences in response to human impacts. In Western Europe, bumblebees and other species respond similarly in terms of total abundance, which was higher in cropland than in secondary vegetation (e.g. in low-intensity cropland, Bombus: z = 6.85, p < 0.0001; other species: z = 9.33, p < 0.0001). However, bumblebee species richness tended to be higher in low-intensity cropland than secondary vegetation (z = 4.68, p = 0.004), which was not true for other bees (z = 1.96, n.s.). Increasing agricultural intensity also resulted in a decline in species richness, but the response was consistent across taxa: low-intensity cropland maintained higher richness than medium-intensity cropland both in bumblebees (z = − 4.30, p = 0.016) and other bees (z = − 3.75, p = 0.042).
Bumblebees also responded differently from other bee species in North America. Bombus species richness was lower in secondary vegetation relative to primary vegetation (z = − 3.93, p = 0.027), but this was not true for other bees (z = 0.86, n.s.); similarly, total abundance was slightly lower for bumblebees (z = 3.64, p = 0.064) but not other species (z = 0.45, n.s.). However, other genera in North America appeared to be relatively more sensitive to medium-intensity cropland, with reduced species richness relative to both primary vegetation (z = − 3.90, p = 0.027) and secondary vegetation (z = − 3.33, p = 0.017). Other genera were also more sensitive to urban areas, with reduced species richness in urban sites relative to primary vegetation (z = − 3.99, p = 0.027), while bumblebees showed no significant response (z = − 2.5, n.s.).
The number of studies showing significant autocorrelation was not significantly higher than the 5% expected by chance (See See Supplementary Table S3 for details).

Discussion
Bees are facing declines across the globe as a result of changing and intensifying land use 4,17,49 . Detailed statistical models that relate bee diversity to drivers of change have the potential to inform mitigation and conservation efforts and to help safeguard food security. However, the transferability of models based on restricted data to other regions and taxa is not guaranteed 11,50 . If responses to threats are context dependent, extrapolation from well known study systems could carry significant risks for biodiversity and food security. The areas where food production is most highly dependent upon animal pollination are also those for which the fewest data are available 11,13,51 , due to a lack of infrastructure and funding in many areas of the world 11 . These same areas are often poorly buffered against the disruption of ecosystem service provision from whatever cause, meaning that effects of any ecological surprises on human well-being could be more severe here than elsewhere. We have shown that bee community responses to land-use change and agricultural intensification can indeed be highly context-dependent, but whether this impacts the transferability of models depends on the facet of diversity that is of interest.
The response of total abundance, Simpson's diversity and species richness of bee communities to land use and intensity (LUI) varied significantly with geographic region, in line with our hypotheses and with previous work in tropical regions 52 . For all response variables, the greatest predictive ability could only be achieved by allowing regional variation in responses to LUI; at the very least, it was necessary to allow regional variation in baseline diversity. This suggests that conclusions based on geographically restricted data cannot reliably be generalized to other regions. Indeed, bee community responses to agricultural intensification varied between regions; only in Western Europe was there an evident decline in diversity with increasing use-intensity of cropland, in line with previous suggestions that agricultural land-use intensity is more important in temperate than in tropical or subtropical systems 10 . The most negative impact of agriculture, however, was seen in South America, where Simpson's diversity was significantly lower than in secondary vegetation; this is congruent with a previous meta-analysis by Gibson et al. 52 that focussed on tropical areas. The effects of urbanisation likewise depended on the subregionwith increased abundance in Africa but few effects seen elsewhere-but these inferences were based on relatively few data. More data from more regions are needed to better understand the impact of urbanisation on bee communities and associated ecosystem services 53 .
It is likely that the geographic variation in responses is in part due to differences in community composition 11 , as we found that taxonomic biases towards bumblebees, which frequently dominate datasets geographically limited to North America and Western Europe, can mask the responses of other species. In Western Europe, for instance, bumblebees had higher species diversity in low-intensity cropland than in secondary vegetation, while other bee species did not show the same effect. Bumblebees have longer flight distances than many smaller bees, so may be better able to persist in more human-dominated land-uses, where foraging resources tend to be further from nesting sites 22 , and can benefit from mass-flowering crops such as oilseed rape 54 .
The effect of taxonomic group on responses to LUI also differed between subregions, suggesting that other factors may also affect generalities. For instance, geographic variation in the nature of threats may be important. Although our land-use intensity classification is applied in an equivalent fashion across regions, it remains extremely coarse. For example, high-intensity cropland may be more intensive in Western Europe than in South America, with regards to some pressures (e.g., pesticide load 55 ) but potentially not others (e.g. spatial extent of monocultures). Such variation in agricultural intensification among regions (even within the same land-use intensity class) could in part be driving observed regional differences in biodiversity responses. More detailed data on different aspects of land-use intensity, such as pesticide load and fertilizer application rates, as well as data on the landscape structure, would enable a more robust and precise analysis of how responses vary across regions. This limitation still highlights, however, that models mostly underpinned by data from regions with a long and intensive history of cultivation are unlikely to provide meaningful inferences for many other regions of the world.
Variation in response among regions could also be driven by differences in community composition and therefore in the distribution of traits that may confer resistance or resilience to human impacts. Previous work has shown that trait-based models of species distributions are only transferable-even within a subregion-when land cover is similar 56 . Transferability across subregions is likely to be even more difficult: variability in the sensitivity of bee communities will be influenced by a complex interaction between the trait distribution (and phylogeny) across communities and variation in the threats they face. For instance, in a global analysis, species that reproduce socially were more vulnerable to isolation and pesticide use than solitary species, but were less sensitive to tillage and agricultural intensity than solitary species 20 ; however, bees that reproduce solitarily are more common in temperate areas than the tropics 57 , while the distribution of these pressures also vary regionally 55 . While we only assessed how a single aspect of community composition may influence results (bumblebees vs other bees), further work into phylogenetic patterns of sensitivity may help to disentangle these two mechanisms that may be driving regional variation in responses. Another important extension to our work would be to explore the interaction of multiple threatening processes, as the pressures faced by bee communities can vary regionally 58 . For instance, competition with introduced species and fragmentation are likely to be more important drivers of native bee diversity in the Neotropics than in temperate regions 59 . While it was not the focus of this work, a spatially-explicit analysis of latitudinal gradients in vulnerability to land-use pressures may be an interesting avenue for further research, potentially highlighting other factors of the environment or community structure that could contribute to geographic variation in sensitivity. For example, species richness of bees peaks at approximately 35° latitude, in dry, Mediterranean climates 60 , rather than in the tropics (as is the case for many other groups 61 ) and this variation in baseline diversity may alter both actual and detected responses to human impacts.
Although our dataset includes over 2000 sites from five continents, it is not a comprehensive compilation of published sources and is still both geographically and taxonomically biased. Africa and Asia in particular are still poorly represented and as a consequence we may still be underestimating the uncertainty in bee responses to land use in these regions. Even biomes that have high bee diversity are underrepresented; for example, only six studies were in the Mediterranean biome although bee species richness tends to peak at this latitude 60 . In addition, the explanatory power of fixed effects was fairly low, as most variation in diversity is explained by methodological differences between studies and sources in most models. Nonetheless, our analysis has important implications for pollinator research and conservation action. We show that results based on geographically and taxonomically restricted datasets may not be transferable to other regions. Responses vary across regions due to a combination of differences in the inherent vulnerability of species and variation in the nature of threats. The provision of pollination services can be influenced by the abundance 62,63 , species diversity 64 and species richness 64-66 of bee communities, although the relative importance of each facet of diversity appears to vary with study system 62,64 . Therefore, if we are to safeguard pollinators and the services they provide, research effort to enhance the representativeness (if not the amount) of available data will be needed to make context-dependent recommendations and to better understand the state of pollination services worldwide.