Spatial mismatch between wild bee diversity hotspots and protected areas

Wild bees are critical for multiple ecosystem functions but are currently threatened. Understanding the determinants of the spatial distribution of wild bee diversity is a major research gap for their conservation. We modeled wild bee α and β taxonomic and functional diversity in Switzerland to uncover countrywide diversity patterns and determine the extent to which they provide complementary information, assess the importance of the different drivers structuring wild bee diversity, identify hotspots of wild bee diversity, and determine the overlap between diversity hotspots and the network of protected areas. We used site-level occurrence and trait data from 547 wild bee species across 3343 plots and calculated community attributes, including taxonomic diversity metrics, community mean trait values, and functional diversity metrics. We modeled their distribution with predictors describing gradients of climate, resource availability (vegetation), and anthropogenic inﬂuence (i.e., land-use types and beekeeping intensity). Wild bee diversity changed along gradients of climate and resource availability; high-elevation areas had lower functional and taxonomic α diversity, and xeric areas harbored more diverse bee communities. Functional and taxonomic β diversities diverged from this pattern, with high elevations hosting unique species and trait combinations. The proportion of diversity hotspots included in protected areas depended on the biodiversity facet, but most diversity hotspots occurred in unprotected land. Climate and resource availability gradients drove spatial patterns of wild bee diversity, resulting in lower overall diversity at higher elevations, but simultaneously greater taxonomic and functional uniqueness. This spatial mismatch among distinct biodiversity facets and the degree of overlap with protected areas is a challenge to wild bee conservation, especially in the face of global change, and calls for better integrating unprotected land. The application of spatial predictive models represents a valuable tool to aid the future development of protected areas and achieve wild bee conservation goals.


INTRODUCTION
Wild bees are of great ecological, economic, and social importance (Potts et al., 2016), but they are threatened by ongoing global change (Goulson et al., 2015). Drivers of wild bee declines are known (Goulson et al., 2015), and evidence shows that some species have suffered extirpations and contractions of their ranges (Cameron et al., 2011). This has led to a loss of species and a potential erosion of functional diversity (Pradervand et al., 2014), challenging wild bee resilience to future changes (Lavorel et al., 2013). Therefore, efficient protection of wild bee taxonomic and functional diversity is urgently needed. However, knowledge of wild bee diversity is severely constrained by existing taxonomic bottlenecks and a lack of comprehensive data sets on wild bee occurrence and traits (Nieto et al., 2014) (but see Woodcock et al. [2014] and Polce et al. [2018]). These constraints apply in Europe, one of the continents on which bees have been extensively studied. Europe hosts around 10% of the world bee diversity (approximately 2000 species), but over 56% of the species are considered data deficient (Nieto et al., 2014). Overall, these data shortages hamper the development of conservation plans (Di Marco et al., 2017;Guisan et al., 2013) because an understanding of the drivers shaping diversity patterns and the spatial distribution of diversity is required to designate protected areas for wild bee conservation (Chowdhury et al., 2023). Information on biodiversity in all its complexity should inform conservation planning (Devictor et al., 2010;Villalta et al., 2022). Taxonomic metrics, such as species richness and species diversity, are widely used to evaluate the importance of drivers of wild bee diversity, including land-use change (Ekroos et al., 2020), disturbance (Winfree et al., 2009), pollution (Moroń et al., 2012), and climate (Bystriakova et al., 2018). Approaches based on functional traits, that is, the phenotypic attributes of an individual that determine its fitness (Violle et al., 2007), are increasingly being used in combination with taxonomic metrics to improve understanding of the variation in species assemblages along ecological gradients (Coutinho et al., 2021) and to predict the consequences of the variation in diversity for ecosystem functions and services (e.g., Fründ et al., 2013). Comprehensive national bee monitoring and functional trait data sets are becoming available (e.g., Fournier et al., 2020;Woodcock et al., 2014) for specific regions, enabling the study of functional trait gradients at higher resolutions. Phylogenetic diversity, a metric of the shared evolutionary history among species (Faith, 1992), represents another biodiversity facet. Unfortunately, a comprehensive phylogeny for bees is still lacking (Hedkte et al., 2013). Examining α and β diversity provides complementary information. Local α diversity typically provides a measure of diversity (i.e., taxonomic, functional, phylogenetic) at the spatial scale at which studies are conducted (Socolar et al., 2016). Conversely, β diversity represents a measurement of the compositional differences among local species assemblages. It is often an important metric in the understanding of diversity change resulting from, for example, environmental gradients and in the study of biodiversity loss and homogenization (Socolar et al., 2016). Among other applications, β diversity can help reveal areas harboring unique species assemblages. Thus, the combination of α and β taxonomic and functional diversity represents an important metric with which to examine diversity change that can be used to inform conservation planning (Socolar et al., 2016).
The determinants of wild bee diversity have been explored. Nonetheless, existing work has limitations regarding spatial coverage, diversity facets included, environmental predictors tested, and completeness of bee diversity data. Recent analyses at the global scale confirm the bimodal latitudinal bee richness gradient proposed by Michener (1979). Specifically, climate and resource availability (i.e., vegetation communities) have arisen as the main drivers of bee diversity worldwide (Bystriakova et al., 2018;Orr et al., 2021) and regionally (Dzekashu et al., 2022;Sponsler et al., 2022). Moreover, anthropogenic stress in the form of land-use change negatively affects wild bee α and β diversity for Palearctic bees at the continental (Bystriakova et al., 2018) andsmaller scales (e.g., De Palma et al., 2017;Moroń et al., 2012). For instance, highly intensified agricultural and urban areas tend to reduce wild bee diversity by selecting species with specific traits that allow them to persist and thrive in those environments, which decreases α and β diversity (Collado et al., 2019;Fournier et al., 2020;Villatla et al., 2022). Urban areas may harbor a greater wild bee diversity than highly intensified agricultural areas (Baldock, 2020) and have higher β diversity (Fournier et al., 2020;Villalta et al., 2022). Finally, anthropogenic stress through the influence of managed pollinators (e.g., honeybees), which may compete with wild bees, has rarely been tested at large spatial scales, but evidence shows a gradual replacement of wild bees by honeybees across the Mediterranean basin (Herrera, 2020). How the determinants of bee distribution act at smaller spatial scales (e.g., countrywide) remains largely unknown, limiting knowledge of wild bee ecology and the capacity to manage and conserve their biodiversity.
Switzerland hosts a relatively high bee diversity compared with other European countries and exhibits large gradients in climate and land-use intensity; therefore, it is a good location to study bee distribution patterns. Bee richness in Switzerland is estimated at 633 species (Ascher & Pickering, 2020), 45% of which are threatened (Duelli, 1994). National surveys of bees have a good spatial coverage, which enables spatial modeling approaches at relatively high resolutions (e.g., Vitasse et al., 2021). The quality and quantity of these data are unparalleled in Europe; few other countries have long-term surveys of bee occurrence and none cover such sharp environmental gradients. Switzerland has developed a network of protected areas aimed at the conservation of the country's biodiversity. As in other countries, protected areas have mostly been designed for plants and vertebrates, and their effectiveness for insect conservation is largely unknown (Chowdhury et al., 2023). For instance, how the protected areas cover the geographic range of wild bees and their diversity hotspots has not yet been assessed We studied bee diversity spatial patterns by mapping the distribution of taxonomic and functional attributes of bee communities based on a unique data set containing records, taken at an unprecedented spatial resolution, of 547 wild bee species in Switzerland. We modeled wild bee taxonomic and functional diversity with sets of predictors that represent gradients of climate, resource availability (i.e., vegetation communities), and anthropogenic stress (i.e., land-use composition and beekeeping intensity). We aimed to uncover countrywide taxonomic and functional α and β diversity patterns and determine the extent to which they provide complementary information; assess the importance of the different drivers structuring wild bee α and β taxonomic and functional diversity; identify hotspots of α and β taxonomic and functional diversity; and quantify the degree of overlap between diversity hotspots and the network of Swiss protected areas.

Study location
Switzerland has pronounced elevational gradients, a broad range of climatic conditions, and a variety of land-use types. Thus, sharp environmental and land-use gradients are characteristic of Swiss landscapes, making them ideal for studying the drivers of spatial patterns in community-level taxonomic and functional diversity.

Bee occurrence, trait values, and community metrics
Occurrence data were provided by the Swiss Biological Records Centre (http://www.cscf.ch, data accessed 12 April 2020). These data originate from community-level surveys performed in 100 × 100 m plots in 2015-2020 in the context of a project focused on updating the red list of Swiss bees (Müller & Praz, unpublished data) (Appendix S1). In total, 6200 plots were surveyed. All quadrats were sampled several times per year by seasoned specialists in bee taxonomy. The number of sampling campaigns per year varied as a function of elevation; higher elevation plots were sampled fewer times because of the shorter vegetation period. Most of the data were taken in a standardized way, but a minority of the samples come from diverse projects that, in some cases, used different methodologies. We included these samples to get maximal cover of the territory. To minimize the risk of varying sampling intensity and ensure the comparability of all the samples, we removed all plots with fewer than 5 species to avoid including undersampled localities (2730 plots). We also removed plots in the highest 2% of quantiles of species richness to avoid potential oversampling (128 plots). In the end, we used data from 3343 plots, containing 52,092 records, that is, wild bee species occurrences (mean [SD] = 15.58 [9.81] wild bee species per plot), and a total of 547 wild bee species (Appendix S2). Our database included more than 98% of the species predicted to be present in Switzerland based on species accumulation curves (Appendix S3). All the data we used were projected on the same 100 × 100 m grid covering all of Switzerland (see Appendix S4 for methodological framework). Our entire data set can be accessed on GBIF (https://doi.org/10.15468/ksfmzj) (Praz et al., 2022).
We selected 8 functionally relevant traits from the European trait database (compiler: Stuart Roberts; pollinator loss module of the EU-FP6 ALARM project, version 01.2017, Table 1) for which sufficient data were available (Fournier et al., 2020). Specifically, we used the following traits and definitions (Table  1): intertegular distance (ITD), feeding specialization, tongue length, beginning of adult flying period (hereafter phenology start), duration of flying period (hereafter phenology duration), nesting mode (hereafter belowground), parasitic behavior (hereafter cleptoparasite), and sociality (hereafter solitary). Details on the traits are in Appendices S2 and S5.
For taxonomic α diversity, we calculated species richness and the Shannon diversity index. For functional α diversity metrics, we first standardized all trait values by subtracting the mean and dividing by the standard deviation to ensure they had the same unit. We then calculated the community-weighted mean value (CWM) of each trait, which can be used to identify major patterns in trait variation in wild bee communities (Woodcock et al., 2014). Moreover, we calculated 3 complementary indices of functional diversity: functional richness based on the trait onion peeling index, evenness based on the trait evenness distribution index, and dispersion, following Fontana et al. (2016) and Laliberté and Legendre (2010). We further calculated β diversity metrics in terms of the local community contribution to taxonomic and functional β diversity (LCBD) with Legendre and De Cáceres's (2013) approach. This approach provides the ecological uniqueness (e.g., taxonomic and functional) of the sites in terms of community composition. Large LCBD values indicate sites that have unusual species combinations, strongly differentiated species compositions, or both relative to sites with low LCBD values (Legendre, 2014). Further information on the definitions and calculation of the metrics is in Appendix S6.
The raw species records are protected by a code of conduct common to all Swiss national data centers. Those data can be ordered via infospecies (https://www.infospecies.ch/ fr/donnees/) according to this deontology. The data set is also available on GBIF.org (https://doi.org/10.15468/ksfmzj) (5-km grid squares) according to a nationally agreed ethical framework. The bee trait database is managed and maintained by Stuart Roberts and enquiries and requests for data can be made to stuart.roberts@cantab.net. Trait data used in this article will be made available in ENVIDAT once the full database from Stuart Roberts is published in Oracle for Research.

Environmental predictors
To infer climate gradients, we considered the main trends in climatic conditions across Switzerland with the 19 Bioclim variables of the CHELSA database set at 1 × 1 km resolution (Karger et al., 2017), which provide information about biologically relevant aspects of climate for 1979−2013. Using these data, we first ran a principal component analysis (PCA) with 100,000 randomly sampled cells (Appendix S7). We then projected the remaining cells onto the PCA. The first 4 PCA axes represented the main trends in climate (i.e., temperature, precipitation, climate seasonality, and temperature range) (92% of the total variation [Appendix S7]) and were selected for further analyses. The resulting maps were downscaled (from 1 × 1 km cells to 100 × 100 m cells) to match the resolution of the other data sets. Details on the PCA axes are in Appendix S8, and Appendix S9 contains maps corresponding to the PCA results.
To infer resource gradients, we mapped the major trends in vegetation across Switzerland with plant occurrence data from the biodiversity monitoring program of Switzerland (InfoFlora). Their data pertain to 500 plots distributed evenly across Switzerland that were surveyed from 2015 to 2019 (detailed plant survey information available from https:// biodiversitymonitoring.ch/index.php/en/). In total, the data set includes 1727 species, representing about half of the Swiss flora.
We ran a PCA of the plant occurrence data to reduce dimensionality of the data set and capture the main trends in vegetation structure and composition. We selected the first 4 PCA axes, representing 43% of the total variation in plant occurrence, for further analyses (PCA1 = 24%, PCA2 = 12%, PCA3 = 4%, PCA4 = 3% [Appendices S10 & S11]). All other axes explained <2% of the variation. We then modeled the distribution of the 4 PCA axes. This step was meant to avoid gaps in the data sets and to transfer the vegetation information onto the same grid as those used for bees. Specifically, each plant PCA axis was modeled individually with the raw 19 CHELSA variables as descriptors (Appendices S7 & S8). The random forest algorithm was trained on 80% of the data and evaluated on the remaining 20%, which were stratified according to the response variable (function createDataPartition in R package caret 6.0-86 [Kuhn, 2008]). Model training and parameter tuning were done using 3 times 3-fold cross-validation (function train in R package caret). The best model was chosen based on the root mean squared error (RMSE), mean absolute error (MAE), and R 2 measured with the trained data set.
The 4 plant PCA axes represented the main trends in plant community composition across Switzerland. The first PCA axis was highly correlated with climate (>0.7) and thus excluded from the analyses. The other 3 axes (PCA2-PCA4) represented gradients of plant communities: from those dominated by coniferous trees to other communities (named midelevation coniferous), from dry to wet vegetation (named dry-wet), and from woody to open plant communities (named forest). (Details on interpretation are in Appendix S10, and results are in Appendices S11 & S12. ) We used land-use and land-cover data from Swiss Area Statistics set at 100 × 100 m resolution (Altwegg & Weibel, 2015) and collected from 2015 to 2018. For modeling, we focused on 3 of the 4 principal domains: urban, agricultural, and forest (details in Appendices S13 & S14). We calculated the percentage of urban, agricultural, and forest cells in 200-, 500-, 1000-, and 2500-m radii around each raster cell center (Appendix S15). We initially considered land-use intensity an important driver of bee diversity based on a country-scale land-use intensity map (Meier et al., 2020) with habitat type and environmental data. Due to high collinearity with several predictors (Appendix S16), however, we did not use land-use intensity in the final analyses.
To assess beekeeping intensity, we used annual data on the spatial distribution of beekeeping locations and the number of hives in Switzerland, which were obtained from the cantonal veterinary offices, from 2012 to 2018. Data were only available for 2012−2014 for the canton of Basel and for 2013−2018 for the canton of Vaud. The data from each veterinary office were checked separately, and only records of beekeeping locations with reliable coordinates were included. We then calculated the number of honeybee hives in 200-, 500-, 1000-, and 2500-m radii around each 100 × 100 m raster cell center (Appendix S17).
Data on the swissTLMRegio protected areas are available from https://opendata.swiss/en/dataset/swisstlmregioschutzgebiete. The DHM25 is available from https:// www.swisstopo.admin.ch/en/geodata/height/dhm25. html#technische_details. Land-use data for Switzerland can be obtained by request from the Federal Office of Statistics (FOS) (www.bfs.admin.ch). Plant data are available from infospecies (https://www.infospecies.ch/fr/donnees/). Data on beekeeping in Switzerland can be obtained from the Cantonal Offices with a confidentiality agreement.

Elevation and protected areas
We used the digital elevation model of Switzerland (DHM25; see Federal office of Topography swisstopo, 2021) to obtain elevation values and to study the elevation patterns of the predictors (Appendix S18). Swiss protected areas include different types of spatial objects that vary in the degree of anthropogenic influence; some of them are subject to strict levels of protection and cover around 27% of the surface (Appendix S19). Hence, we classified the protected areas in 2 groups. First, we considered only those areas that are strictly protected. These protected areas include only the biotopes of national importance (dry grasslands, fens, bogs, amphibian reproduction sites, floodplains), the Swiss National Park, Ramsar areas, forest reserves, and private forest and nature reserves (i.e., owned by Pro Natura), which cover around 5.6% of the Swiss surface (Table S4; Delarze et al., 2016). Second, we considered protected areas with less protection and variable land-use intensity and anthropogenic interference, referred as protected multipleuse areas (Appendix S19). These areas cover approximately 21% of the Swiss surface. We also combined the 2 classes of protected areas into a single category of all protected areas.

Variable selection
We selected variables that had Pearson intercorrelations <0.7 (Appendix S16). For the variables calculated using 200-, 500-, 1000-, and 2500-m buffers (i.e., beekeeping, land use, and land-use intensity), we also ran preliminary analyses including all variables to assess which neighborhood windows were the best descriptors of the bee diversity metrics. The variables calculated at 2500 m around the raster cell center were the best predictors.
After the variable selection, we retained climate (4 PCA axes), vegetation (3 PCA axes), beehive density (number of beehives in a 2500-m radius), and land use (percentage of urban, agricultural, and forest in a 2500-m radius) for further analyses. These variables represented gradients in climate, resource availability, and habitat amount and disturbance, respectively. All variables were scaled and centered prior to analyses.

Statistical analyses
All analyses were conducted in R 4.0.3 (R Core Team, 2020). Flat violin plots were created using the raincloud package in R (Allen et al., 2019). All α and β taxonomic and functional diversity metrics (species richness, Shannon diversity, functional richness, trait evenness, dispersion, LCBD taxonomic, and LCBD functional) and the CWM of the 8 studied traits were modeled using the selected variables (see "Variable selection"). Model calibration, parameter tuning, and model performance and selection were done following the same steps as for the vegetation structure and composition (see "Environmental predictors"). Three algorithms were tested: generalized linear models, neural networks (Venables & Ripley, 2002), and random forests (Liaw & Wiener, 2002). Model performance on the test data was estimated with 3 metrics: RMSE, MAE, and R 2 . Random forest analyses produced the best results for each diversity metric: lowest RMSE and MAE, and the highest R 2 (random forest R 2 was always >0.9 [Appendix S20] and was used to produce countryscale predictions of diversity patterns [Appendices S20-S22]). The importance of each descriptor in the random forest models was estimated as the averaged difference across all trees (before and after permuting a variable) in RMSE computed on the out-of-bag data and normalized by the standard error (function varImp, package caret). We used partial dependence plots (PDPs) to assess the changes in diversity metrics while limiting the influence of other descriptors (Friedman, 2001;Greenwell, 2017). A PDP is especially useful for visualizing the relationships revealed through complex machine learning algorithms, such as random forests (Greenwell, 2017). In particular, PDPs help visualize the relationship between 1 specific predictor and the response variable while accounting for the average effect of all other predictors in the model.
We identified hotspots of wild bee α and β taxonomic diversity and functional diversity with the maps of predicted diversity metrics. Specifically, we retained only those cells containing diversity values from the upper 10th percentile of their distribution.

Wild bee spatial diversity patterns in Switzerland
We found countrywide patterns in wild bee taxonomic and functional diversity in Switzerland (Figure 1; Appendix S23). Specifically, wild bee diversity metrics were strongly structured along elevational gradients (Figure 1a-e; Appendix S23); taxonomic and functional α diversity decreased at higher elevations (>2000 m asl) (Appendix S24). Community attributes showed high Pearson's correlations not only between species richness and Shannon diversity (r = 0.86) (Figure 1h), but also between species richness and functional richness (r = 0.83). Although functional evenness generally declined with elevation, there was also a local maximum at approximately 1300 m asl (Appendix S24). The predicted CWM of the 8 selected traits also showed important shifts as elevation increased (Appendices S23 & S25). Specifically, high-elevation wild bee communities had a lower proportion of belowground-nesting, cleptoparasitic, and solitary bee species (Appendix S25). The proportion of short-tongued species (Appendix S25) and the proportion of generalist species decreased (Appendix S25), whereas body size increased (Appendix S25). Taxonomic LCBD and functional LCBD, which provided an indication of the ecological uniqueness of the cells, were correlated (r = 0.63) and negatively associated with α diversity metrics, such as richness (r = −0.71), Shannon diversity (r = −0.83), and dispersion (r = −0.80). Taxonomic and functional uniqueness increased as elevation increased, resulting in a negative correlation with taxonomic and functional richness. This pattern was stronger for taxonomic than for functional LCBD (Figure 1f,g).

Drivers of wild bee diversity
We assessed the importance of each predictor individually from the 4 categories considered, that is, climate, vegetation, land use, and beekeeping. Climatic variables, represented by the first 4 axes of the climate PCA, were always among those FIGURE 2 Importance of environmental descriptors (y-axis) as predictors of α and β taxonomic and functional community attributes of wild bees in Switzerland: (a) species richness, (b) Shannon diversity index, (c) functional richness, (d) functional evenness, (e) functional dispersion, (f) local community contributions to taxonomic β diversity (LCBD taxonomic), and (g) local community contributions to functional β diversity (LCBD functional) (the longer the bar, the better the predictor). Importance values divided by maximum value to obtain a comparable range from 0 to 1. Climate and vegetation variables represent the principal component analysis (PCA) axes (PCA 1−4 for climate and PCA 2−4 for vegetation), representing 17% of the variation (details in Methods and Appendices S7-S12). explaining the highest proportion of variance for all responses, followed generally by vegetation, land-use, and beekeeping metrics (Figure 2; Appendix S26). For climate, temperature range (PCA4) (Appendix S8) explained a large proportion of the variation in taxonomic and functional α diversity metrics (Figure 2a-e), whereas temperature and precipitation gradients explained a large part of the variation in the taxonomic and functional LCBD β diversity (Figure 2f,g; Appendix S8) and in most CWMs, including the proportion of belowground-nesting, cleptoparasitic, and solitary bee species, ITD, phenology dura-tion, and tongue length (Appendix S26). Higher temperature values (lower PCA scores) enhanced wild bee taxonomic diversity and functional richness (Figure 3 climate PC1), but the patterns were less clear for the remaining α and β functional diversity metrics (Figure 3 climate PC1) and for the CWM of the studied traits (Appendix S27). Xeric conditions boosted species and functional richness, Shannon diversity, and functional LCBD β diversity (Figure 3 climate PC2), particularly in the lowlands of southwestern (canton of Wallis) and southeastern (canton of Graubünden) Switzerland. Finally, wild bee diversity metrics generally decreased in areas with higher temperature ranges (Figure 3 climate PC4).
Wild bee diversity metrics and traits were structured along the 3 PCA axes depicting different gradients of vegetation change. First, changes in the plant communities based on their drought tolerance (PCA2 dry-wet) (see also Appendix S11) influenced several wild bee diversity metrics. Plant communities that were more drought tolerant boosted wild bee species richness, functional richness, Shannon diversity, and taxonomic β diversity, whereas functional evenness and dispersion peaked at intermediate values. Furthermore, drylands had more solitary, small, and belowground-nesting bees (Appendix S27). Second, changes from open vegetation to forest had a small effect on wild bee diversity metrics, with the exception of functional evenness, which peaked at intermediate values. Finally, plant PCA2, representing plant communities in coniferous forest at midelevation (Appendix S11), mostly affected functional evenness and functional β diversity (Figure 3).
Wild bee diversity metrics were influenced differently by the 3 types of land use (i.e., agricultural, forest, urban) (Figure 2; Appendix S26). Wild bee species and functional richness were lowest in landscapes with about 50% agricultural area (Figure 3). Landscapes containing higher proportions of agricultural land and assemblages with shorter tongues and lower proportions of belowground-nesting and solitary wild bees were associated with reduced functional evenness ( Figure 3) and taxonomic LCBD β diversity (Appendix S27). Wild bee species and functional richness increased with larger proportions of forest in the landscape (Figure 3). Conversely, Shannon diversity, functional evenness, and functional dispersion peaked in landscapes with intermediate levels of forest ( Figure 3). Furthermore, in landscapes increasingly dominated by forest, wild bee communities had higher proportions of belowground-nesting and solitary species and species with longer tongues (Appendix S27). Finally, increasingly urbanized landscapes enhanced wild bee species and functional richness, but reduced Shannon diversity, functional evenness, and functional dispersion ( Figure 3). Strikingly, the local community contributions to functional β diversity (LCBD functional) increased when landscapes became more urbanized (Figure 3), suggesting original combinations of traits, whereas taxonomic LCBD remained largely unaffected ( Figure 3). Furthermore, urban areas triggered changes in trait values. For instance, the proportion of belowgroundnesting bees and feeding specialization decreased as urban area increased, whereas the proportion of cleptoparasitic and solitary bees and the phenology duration increased (Appendix S27).
Beekeeping was the predictor explaining the lowest percentage of variance for all the responses modeled ( Figure 2); some metrics were completely unaffected, such as dispersion ( Figure 3) and phenology duration (Appendix S27). Nonetheless, species and functional richness and Shannon diversity increased as beekeeping increased and then stabilized at high values ( Figure 3). Interestingly, trait evenness, taxonomic β diversity, and feeding specialization decreased as the number of beehives increased (Figure 3; Appendix S27).

Identification of hotspots
Predicted species richness, Shannon diversity, functional richness, functional evenness, and functional dispersion were greater in urban areas than in forest and agricultural areas (Appendix S28), whereas taxonomic and functional LCBD were lower in urban areas (Appendix S29). At mid and high elevations, urban land use represented a very small proportion of the total land use (Appendices S30 & S31) and there were no major settlements. But these elevations had other types of urban land use (Appendix S13) that hampered comparisons. Further, at low elevations, agricultural and urban areas had a similar species richness (Figure 4). Agricultural areas harbored larger Shannon diversity and functional dispersion (Figure 4), whereas urban areas had a higher wild bee species richness, functional richness, and functional evenness. In urban areas, β diversity was also higher (Figure 4).

Protected areas
Strictly protected areas and protected multiple-use areas had different elevational distribution patterns. Protected multiple-use areas exhibited a clear positive elevational gradient (Figure 5b), whereas strictly protected areas had a unimodal distribution along elevation (peaking at intermediate elevations of 1000−2000 m asl) (Figure 5d).
The majority of α and β diversity hotspots (cells with upper 10% diversity values) were predicted in unprotected areas (Figures 5-7; Appendix S36). This was particularly prominent for α diversity metrics (species richness, Shannon diversity, functional richness, and dispersion), typically found in the lowlands. Over 75% of the hotspots were predicted in unprotected land ( Figure 5). In the case of β diversity hotspots, the overlap with protected areas was substantially higher (Figures 5a &  7a). For instance, nearly 50% of the taxonomic LCBD hotspots, typically found in highlands, were predicted in protected land (Figure 7a).
The overlap between hotspots of the different diversity facets of wild bee communities varied substantially depending on the degree of protection ( Figure 5). We found α diversity hotspots were better represented in the strictly protected areas (Figure 5c), with the proportion of protected cells containing the highest diversity metrics (e.g., 1-10% of the diversity gradient) (Figure 5c), larger than the national average (approximately 5.6%). Conversely, β diversity hotspots were better represented in protected multiple-use areas (Figure 5a). Hotspots of α functional metrics were equally distributed in strictly protected areas (∼49%; Figure 7b) and protected multiple-use areas (Figure 7b). However, for the remaining metrics, diversity hotspots were always in higher proportions in protected multiple-use areas, especially for β diversity hotspots (Figure 7b). Finally, regarding the occurrences of the bee species in the community plots, 7 species were only sampled in strictly protected areas and 9 species in protected multiple-use areas (Appendix S36). FIGURE 4 Differences in the α and β taxonomic and functional metrics for wild bees in Switzerland among the 3 main land-use types (agricultural, urban, and forest) at low elevation (197−1000 m asl): (a) species richness, (b) Shannon diversity index, (c) functional richness, (d) functional evenness, (e) functional dispersion, (f) local community contributions to taxonomic β diversity (LCBD taxonomic), and (g) local community contributions to functional β diversity (LCBD functional) (vertical lines, 95% confidence interval of the median; horizontal lines, median; data shown are from 10th to 90th percentiles to facilitate comparison of box graphs).
(h) Land-use composition at low elevations. Additional plots of the α and β taxonomic and functional diversity metrics and of the community-weighted means of the 8 wild bee traits at low, medium, and high elevations are in Appendices S28-S34.

DISCUSSION
Wild bee diversity metrics followed clear elevational gradients. Species richness, functional richness, Shannon diversity, and functional divergence declined monotonically as elevation increased; taxonomic and functional β diversity increased as elevation increased; and functional evenness exhibited a humpshaped relationship. Elevational gradients represent one of the most prominent ecological gradients, in which plant assemblages (e.g., Pellissier et al., 2010) and other taxonomic groups (e.g., grasshoppers  and bumblebees [Sponsler et al., 2022]) sift. Elevational gradients select for specific phenotypes, which generates multiple patterns, including midelevation peaks (hump shapes), monotonic declines, and unimodal distributions, depending on the group and context (Rahbek, 1995). The number of wild bee species is typically lower in the highlands of mountain slopes because changes in climatic variables, a shortened phenological season, and reduced productivity and resource availability constrain species persistence (Rahbek, 1995).
Nonetheless, high-elevation communities tend to be ecologically unique because they are composed of species that FIGURE 5 Proportion of taxonomic and functional diversities included in protected areas that (a) include areas with variable management intensity and anthropogenic influence and relatively lower protective measures and (c) include only those with strict protection measures (Appendix S19). For each diversity metric, each cell is ranked from most to least diverse based on quantile values. The lower the x-axis value, the more diverse the cells (horizontal black line, proportion of protected cells in all of Switzerland [protected multiple-use areas approximately 21%; strictly protected areas approximately 5.6% of the surface]; taxonomic α diversity metrics: species richness and Shannon diversity; functional α diversity metrics: functional richness, functional evenness, and functional dispersion; taxonomic β diversity metrics: local community contributions to taxonomic diversity [LCBD taxonomic]; functional β diversity metrics: local community contributions to functional diversity [LCBD functional]). The graph format is based on Devictor et al. (2010). Proportion of protected cells along an elevation gradient in (b) strictly protected areas and (d) protected multiple-use areas (gray on maps, distribution of protected areas in Switzerland). In Switzerland, wild bee (e) α taxonomic and functional diversity hotspots and their overlap and (b) taxonomic and functional β diversity hotspots and their overlap. Hotspot cells contain upper 10th percentile of the distribution of all β and α diversity metrics. display specific adaptations to survive at high elevations and consequently contribute disproportionately to large-scale taxonomic and functional β diversity (e.g., invertebrates [Fontana et al., 2020], lichens [Nanda et al., 2021]). In the case of wild bee communities in temperate latitudes, high-elevation areas contain multiple wild bee species, including several highelevation specialist bumblebee species, which have been studied in Switzerland (Pellissier et al., 2013). High-elevation bumblebees are highly sensitive to heat stress (Pellissier et al., 2013), which might extend to other high-elevation bees. Hence, global warming threatens to trigger drastic declines of these taxa due to rising temperatures and, potentially, increased competition with migrating bees from lower elevation (Pradervand et al., 2014). Such a change has been observed for plant communities (Vitasse et al., 2021) and herbivorous insects (e.g., grasshoppers; Descombes et al., 2021).
Consistent with the xeric hypothesis (Michener, 1979), wild bee diversity thrived in warm, xeric areas, such as dry grasslands. At a global scale, wild bee diversity is highest in these climatic regions (e.g., in Mediterranean-type regions) (Orr et al., 2021). We found that this also applies to much smaller spatial scales in the Alps. Xeric areas might recreate conditions in which the clade originated in Gondwana (Michener, 1979), and bees may therefore possess adaptations to such environments (Michener, 1979). Xeric areas tend to have large plant species pools and to promote specialized interactions (Minckley et al., 2000). In the Alps, dry grasslands are well-known biodiversity hotspots but are increasingly threatened by climate and land-use change (Boch et al., 2019).  We found urban and agricultural land uses promoted distinct diversity metrics. Specifically, although Shannon diversity and functional dispersion increased as urbanization increased, species richness, functional richness, and functional evenness increased as agriculture increased. Furthermore, our results showed lower β diversity in cities than in forest and agricultural areas. This is in line with recent research showing that urban features (e.g., technological innovations, gardening techniques) tend to be similar among different cities (Alberti, 2015), thereby leading to a convergence and homogenization in the selected phenotypes and species (Groffman et al., 2014). While cities can harbor diverse bee communities (Theodorou et al., 2020), they simultaneously exert a strong filter for certain traits (Fournier et al., 2020) and might even represent a subset of species from agricultural areas (Banaszak-Cibicka &Żmihorski, 2020).
We did not find clear evidence that beekeeping negatively influences bee diversity, even though Switzerland has one of the highest hive densities in Europe (around 4 hives/km 2 in 2019 [Charriere et al., 2018]). This might be due to the elevational structure of beekeeping, which is far more intense at low elevations, where bee richness is also highest. Beekeeping might be in areas with high resource availability. Our data could not be used to consider the temporal and spatial variability in urban beekeeping, including the movement of hives uphill during the summer season. In any case, wild bee declines might lag and be gradual, thus requiring temporal data, as observed in the Mediterranean basin (Herrera, 2020). This is concerning because beekeeping is changing fast in Switzerland. For instance, urban beekeeping doubled in Swiss cities from 2012 to 2018 (Casanelles-Abella & Moretti, 2022).
Our study comes with limitations. First, we did not model individual species occurrences, but rather focused on community attributes. Although species distribution models (SDMs) provide species-specific responses to environmental gradients, they require a sufficient number of observations to be reliably applied , which in our case would have excluded many rare species. Conversely, modeling community attributes makes full use of community data and helps highlight general diversity trends. Nonetheless, wild bee SDMs are a necessary next step to further understand bee ecology and forecast changes in distribution ranges due to global change drivers (e.g., climate, land use), as has been done for a large number of taxa (Isaac et al., 2020). Second, we assessed resource availability gradients through vegetation patterns because we could not use direct metrics of plant diversity or nectar or pollen productivity, both important drivers of bee diversity at different scales (Orr et al., 2021). Developing high-resolution countrywide maps of plant diversity and nectar and pollen productivity, including spontaneous and horticultural non-native plants, will be key to improving predictions of wild bee diversity and occurrence. Third, we did not explore wild bee phylogenetic metrics, which might have different patterns than functional and taxonomic diversities (Devictor et al., 2010) due to a lack of a comprehensive bee phylogeny (Hedtke et al., 2013). Fourth, we did not consider the extent of overlap among wild bee diversity, ecosystem functions (especially pollination), and ecosystem services. Finally, we did not explore the spatial configuration and connectivity of protected areas, which might play a key role in bee biodiversity (Schüep et al., 2010).
Our results call for the implementation of multiple-faceted approaches to inform biodiversity conservation planning and strategies. The spatial mismatch in the distribution of the different diversity facets and scales shows the necessity of considering multiple dimensions to study and protect biodiversity in all its complexity (Pollock et al., 2020). In that regard, strictly protected areas and protected multiple-use areas provide complementary protection on different biodiversity dimensions. On the one hand, strictly protected areas, which have highly regulated anthropogenic intervention and that contain habitats with high plant richness (Delazare et al., 2016), serve as highly valuable habitat for bee α diversity. On the other hand, protected multiple-use areas, which are mostly located in higher elevations, also represent an important infrastructure to protect bee ecological uniqueness, although they experience more anthropogenic influence.
Our study points out the role of unprotected land in wild bee conservation. Although α diversity hotspots were well represented in strictly protected areas, these areas cover a small surface and some of them are under stress due to global change (Boch et al., 2019). For instance, in Switzerland the habitat quality of dry grasslands is decreasing at an alarming rate (Boch et al., 2019). Furthermore, protected multipleuse areas are often found at high elevation, areas that are expected to be substantially, negatively affected by changing temperatures. High-elevation areas may also have their plant assemblages modified and be colonized by lowland species, which will likely negatively affect high-elevation specialist bees (Pradervand et al., 2014). Therefore, the existing network of unprotected ecological infrastructure (different natural, seminatural, and unique habitats) in anthropogenic landscapes (e.g., agricultural areas, cities) or in seminatural areas represents a key element for further safeguarding wild bee diversity, particularly if properly integrated and planned within the strictly protected areas. Unprotected areas of value for wild bee conservation include a myriad of land-cover types, such as flower strips, urban gardens, and ruderal sites, that are often quite sensitive to globalchange drivers, such as land-use change, global warming, and invasion by non-native species. It will be important to quantify the specific contributions of this network of unprotected areas to wild bee diversity.
The Post-2020 Global Biodiversity Framework draft stresses the importance of protected and conservation areas to mitigate ongoing biodiversity loss and sets an ambitious goal of protecting 30% of all land by 2030 (Convention for Biological Diversity, 2020). Nonetheless, without an evaluation of the performance of protected areas, these actions may ultimately have a limited impact on reducing biodiversity declines. Protecting insect diversity is challenging because the spatial distribution of insect multifaceted diversity remains elusive for many groups (Chowdhury et al., 2023). In that regard, our study provides a practical example of the already stated importance of spatial predictive models to overcome the distributional gaps and contribute to the planning and evaluation of biodiversity conservation measures (Pollock et al., 2020). To achieve the global conservation goals, countrywide prioritization of existing and future protected areas emerges as a major issue due to the potential spatial mismatch among biodiversity facets and the uncertainty of the conservation value of some protected areas under future conditions. Solving these problems represents a pressing issue in the face of future global change.

ACKNOWLEDGMENTS
We acknowledge financial support from the DFG (FO 1420/1-1, project FunShift) and the WISNA program from the German Federal Ministry of Education and Research to B.F. J.C.A. was supported by the Swiss National Science Foundation (project 31BD30_172467) within the program ERA-Net BiodivERsA project BioVEINS: Connectivity of Green and Blue Infrastructures: Living Veins for Biodiverse and Healthy Cities (H2020 BiodivERsA32015104) and by the Swiss Federal Office for the Environment (FOEN) in the framework of the project City4Bees (contract 16.0101.PJ/S284-0366). We also acknowledge Göhner Stiftung for supporting this project (2019-2917/1.1). Finally, we thank Pro Natura for providing the maps on their forest and nature reserves. The bee data set was assembled with funding from the Federal Office for the Environment to info fauna (Neuchatel, Switzerland) within a project coordinated by C. Praz and A. Müller. We thank K. Bollmann for his input on Swiss protected areas.
Open access funding provided by ETH-Bereich Forschungsanstalten.