Wild bees as winners and losers: Relative impacts of landscape composition, quality, and climate.

Wild bees, like many other taxa, are threatened by land-use and climate change, which, in turn, jeopardizes pollination of crops and wild plants. Understanding how land-use and climate factors interact is critical to predicting and managing pollinator populations and ensuring adequate pollination services, but most studies have evaluated either land-use or climate effects, not both. Furthermore, bee species are incredibly variable, spanning an array of behavioral, physiological, and life-history traits that can increase or decrease resilience to land-use or climate change. Thus, there are likely bee species that benefit, while others suffer, from changing climate and land use, but few studies have documented taxon-specific trends. To address these critical knowledge gaps, we analyzed a long-term dataset of wild bee occurrences from Maryland, Delaware, and Washington DC, USA, examining how different bee genera and functional groups respond to landscape composition, quality, and climate factors. Despite a large body of literature documenting land-use effects on wild bees, in this study, climate factors emerged as the main drivers of wild-bee abundance and richness. For wild-bee communities in spring and summer/fall, temperature and precipitation were more important predictors than landscape composition, landscape quality, or topography. However, relationships varied substantially between wild-bee genera and functional groups. In the Northeast USA, past trends and future predictions show a changing climate with warmer winters, more intense precipitation in winter and spring, and longer growing seasons with higher maximum temperatures. In almost all of our analyses, these conditions were associated with lower abundance of wild bees. Wild-bee richness results were more mixed, including neutral and positive relationships with predicted temperature and precipitation patterns. Thus, in this region and undoubtedly more broadly, changing climate poses a significant threat to wild-bee communities.

INTRODUCTION

Wild bees provide critical pollination services for agriculture and natural ecosystems but are threatened by land‐use and climate change (Cardoso & Gonçalves, 2018; Kerr et al., 2015; Potts et al., 2010, 2016; Settele et al., 2016). However, the impacts of land‐use and climate change are manifold, and bee species span an array of behavioral, physiological, and life‐history traits that can increase or decrease resilience to land‐use or climate change (Benjamin et al., 2014; De Palma et al., 2015; Harrison et al., 2018). Indeed, within a given region, populations of some species are declining while others are stable or increasing (Bartomeus et al., 2013a; Powney et al., 2019; Soroye et al., 2020). Understanding how land use, climate, and species characteristics interact is critical to predicting and managing pollinator populations, but these studies require datasets spanning diverse sampling locations, years, and pollinator communities.

Anthropogenic land‐use change can have complex and cascading effects on bee communities (Potts et al., 2010). Developed and agricultural landscapes typically have lower wild‐bee richness and abundance than semi‐natural landscapes (Kennedy et al., 2013). This pattern likely results from reduced floral and nesting resources for wild bees in anthropogenic habitats, although other factors like pesticide use and predation risk also contribute to spatial variation in bee communities (Goulson et al., 2018; Roulston & Goodell, 2011). In addition to reducing bee abundance, habitat loss and fragmentation are associated with reduced bee body size, phylogenetic diversity, and pollination services provided by wild bees (Grab et al., 2019; Renauld et al., 2016; Warzecha et al., 2016). Compared with natural habitats, plant–pollinator networks in urban areas generally have fewer plant–insect interactions and increased generalism of floral visitors (Geslin et al., 2013). This shift to highly generalized networks with less redundancy may mean pollination services in urban areas are more vulnerable to future disturbance events (Memmott et al., 2004), although some habitats within cities support a high diversity of flowering plants and pollinating insects (Baldock et al., 2019).

Concurrent with land‐use change, climate change is driving increasingly unpredictable weather patterns (Lynch et al., 2016; Thibeault & Seth, 2014), which could influence communities of wild bees in several ways (Forrest, 2016; Rafferty, 2017). Changes in seasonal weather patterns coupled with increased frequency of extreme weather events could alter availability of floral resources (Phillips et al., 2018), cause phenological mismatches between plants and their pollinators (Bartomeus et al., 2011, 2013a; Kudo & Ida, 2013), shift weather‐dependent activity patterns and foraging behavior (Straka et al., 2014), and potentially increase pathogen spread and dominance of non‐native species (Settele et al., 2016). Heat and drought stress can reduce flower number and nectar volume and concentration, with significant implications for plant–pollinator interactions (Descamps et al., 2018; Gallagher & Campbell, 2017; Mu et al., 2015). In montane meadows in Colorado USA, climate change altered timing and abundance of floral resources, which was positively correlated with interannual abundance of bumble bees (Aldridge et al., 2011; Ogilvie et al., 2017). In low and mid‐elevation environments, however, considerable uncertainty remains about how wild bees are responding to temperature and precipitation drivers.

Both within and among taxa, wild bees vary considerably in how they respond to land‐use and climate change (Cariveau & Winfree, 2015). Of studies included in a recent review, 42% found a negative relationship between pollinators and land‐use change while 45% were neutral and 13% were positive (Winfree et al., 2011). Bumble bee (Bombus spp.) ranges in North America and Europe have contracted because they failed to move northward as southern sections of their range became too warm (Kerr et al., 2015; Sirois‐Delisle & Kerr, 2018). However, other taxa, such as Ceratina australensis, an arid‐adapted, Australian species, should expand their range under future climate conditions, in this case, due to increasing area of arid habitat (Dew et al., 2019). Physiological differences may be particularly important in determining the response of species to altered climate. For example, compared with wild‐bee species that overwinter as larvae or pupae, species that overwinter as adults have higher metabolic rates during diapause, and, when exposed to warmer winter temperatures, lose more weight, and emerge earlier (Fründ et al., 2013; Schenk et al., 2018; Slominski & Burkle, 2019).

To implement effective conservation strategies, we need to move from considering single stressors to quantifying multiple, potentially interacting pressures on wild‐bee communities (González‐Varo et al., 2013; Goulson et al., 2015; Vanbergen & The Insect Pollinators Initiative, 2013). Most studies have evaluated only land‐use or climate effects, not both (but see Oliver et al., 2015). To address these critical knowledge gaps, we analyzed a long‐term dataset of wild‐bee occurrences from Maryland, Delaware, and Washington DC, USA, examining how different bee species and communities respond to land‐use and climate factors. We asked the following research questions: (1) What is the relative importance of land use and climate in shaping wild‐bee abundance and species richness? (2) Do functional traits explain differential responses of wild‐bee species to land use and climate? and (3) Does life‐cycle phenology explain differential responses of wild‐bee species to land use and climate? Since both land use and climate influence wild‐bee communities, we expected some variables from both predictor sets would correlate with bee abundance and richness. Specifically, we hypothesized that bees with a short season of activity, solitary life history (or eusocial species in a solitary phase), and univoltine lifecycle would be most affected by loss of floral resources associated with more intensive land use. In our study region, these traits typify the bee community active in spring; thus, if our hypothesis is supported, spring bees will have a stronger negative response to anthropogenic land uses than bee species foraging in the summer or fall (autumn). Similarly, we predicted lower abundance of spring bees overwintering as adults following warm winters. In laboratory studies, compared with species overwintering as larvae, species overwintering as adults have higher mortality with warmer overwintering temperature (Fründ et al., 2013; Schenk et al., 2018; Slominski & Burkle, 2019).

METHODS

United States Geological Survey, Native Bee Inventory, and Monitoring Lab (BIML) dataset

To investigate wild bee responses to land use and climate, we cleaned and analyzed a dataset of wild‐bee occurrences from Maryland, Delaware, and Washington DC USA. From 2002 to 2016, the United States Geological Survey's Native Bee Inventory and Monitoring Lab (BIML) collected, identified, and archived over 100,000 wild bee observations across this region (Droege & Sellers, 2017; Figure S1). To prepare BIML occurrence data for analyses, we (1) extracted only records generated through pan trap sampling methods, (2) cleaned and verified species binomials, time, and locality information, and (3) generated identifier variables, sampling method, and effort (Kammerer et al., 2020a, 2020b). To focus on unmanaged, wild bee species, we removed all occurrences of honey bees, Apis mellifera (L.).

Starting with occurrence‐level BIML data (Kammerer et al., 2020a, 2020b), we filtered data to years and sites suitable for our analyses. For species‐richness analyses, we retained all wild‐bee occurrences (Kammerer et al., 2020a, 2020b, dataset #1), but for abundance models, we only used occurrences with known sampling effort (e.g., number of traps and length of sampling; Kammerer et al., 2020a, 2020b, dataset #2). Before 2003, while developing their monitoring program, BIML employed a variety of sampling methods. In 2014, they made a substantial change to their standard method using propylene glycol as a preservative, longer sampling times, and larger traps. To reduce variation in sampling method, we only included data from 2003 to 2013, and incorporated remaining variation in sampling method (trap color and volume) in our analyses (see Kammerer et al., 2020b for a detailed description of BIML sampling methodology over time).

Study region

Our study region encompasses natural, agricultural (arable), and developed land, with approximately 48, 28, 19%, respectively, in each category (USDA NASS, 2016, see Table S1 for our classification of USDA NASS land‐use classes as agricultural, natural, or developed). Deciduous forest is the primary type of natural land (USDA NASS, 2016), although the Eastern Shore (eastern Maryland and Delaware) also has large areas of wetland and coastal habitats. Historically, forest in our study region was dominated by oak and chestnut, but more recently, Maryland, Washington DC, and Delaware were classified within the compositionally diverse Mesophytic forest region, with western Maryland demarcated as part of the Appalachian‐oak section (Dyer, 2006). Agricultural land cover in this region is primarily field crops (corn, soybeans, and winter wheat), hay, or pastureland (USDA NASS, 2016).

Like the broader study region, BIML sampling locations represent varying landscape composition and local habitat types. The BIML dataset includes highly natural, agricultural, and urban landscapes, with landscape composition ranging from 0 to 97% natural, 0 to 83% arable, and 0 to 99% developed land (USDA NASS, 2016) within 1750 m of each sampling location (Figure S2, see Section 2.4.1 below). At locations of BIML sampling, local land cover was similarly variable. Approximately 39%, 16%, and 45% of BIML sampling locations (summarized by unique site and year combinations) were natural, agricultural, and developed habitat, respectively (USDA NASS, 2016). Developed/open space, deciduous forest, and developed/low intensity (USDA NASS, 2016) were the most common types of local habitat, representing 30, 17, and 11%, respectively, of all site‐years in the BIML dataset (Figure S3). Developed/open space, deciduous forest, and grass/pasture were the most intensively sampled (Table S2). Due to substantial variation across sites, we included landscape composition and the land‐cover class at the location of bee sampling in all our analyses (see Section 2.7). Excepting some counties with increasing urban land, our study region experienced relatively modest land‐use change during the study period (Yang et al., 2018), but the varied local and landscape conditions in the BIML dataset allowed us to quantify the effect of more extreme future change in land use on wild‐bee richness and abundance.

Bee phenology and functional traits

To examine seasonal patterns in bee responses to land use, we divided BIML data into “spring” and “summer/fall” seasons. While many previous studies examined land‐use and climate effects on season‐total bee richness or abundance, it is clear that flowering plant and bee communities, and their interactions, vary dramatically with season (CaraDonna et al., 2017; Harrison et al., 2018; Leong et al., 2016; Mandelik et al., 2012; Russo et al., 2013). In the eastern US, some wild‐bee species emerge very early in the season and finish provisioning their offspring by mid to late May while others are present from mid‐spring until fall, some with multiple generations per season (Adamson et al., 2012; Figure 1). Similarly, the quality and quantity of floral resources in habitats are not static within a growing season (Ogilvie & Forrest, 2017). In our study region, forests contain important early resources (e.g., flowering trees and understory plants) for wild bees, but have few flowering plants after canopy closure (Kammerer et al., 2016).

FIGURE 1

Phenology of the most common wild‐bee genera in the mid‐Atlantic USA (genera representing least 1.5% of total abundance). For more interpretable visualization, we show the proportional abundance of each genus in 34 seven‐day intervals

To define seasons relevant for wild bees, we characterized bee phenology based on relative abundance over time. We simplified the BIML‐recorded sampling date by translating sampling date into 1‐week intervals. Then, we defined seasonal groups by (1) visually examining relative abundance of genera (n = 11 genera) that represented at least 1.5% of the total abundance of the BIML dataset (Figure 1) and (2) visualizing community composition across sampling dates with an NMDS ordination analysis. Using the ampvis2 package in R (Andersen et al., 2018; R Core Team, 2019) with a Bray–Curtis distance measure, we conducted a NMDS ordination of relative abundance of all bee species. We found composition of bee communities sampled from June to September was extremely similar, with some overlap between June–September and October/November (Figure S4). But, bee communities sampled in March and April were entirely distinct from those active in June–November. Sampling dates in early May were more similar to March and April, while communities sampled in late May clustered with June–November (see https://land‐4‐bees.github.io/DroegeLandscapeClimate/ for an interactive version of Figure S4 to differentiate early vs. late May). Based on these temporal clusters and the need to maintain adequate sample sizes in each season, we split the sampling dates into two seasons, “spring” and “summer/fall.” We defined “spring” as any transect where the sampling period midpoint was before 15 May. We considered the rest of the year “summer/fall.” For statistical analyses, to combine multiple transects sampled at the same site, we calculated mean abundance day−1 trap−1 per site in each season.

We separated the bee community into functional groups by characterizing the nesting location, sociality, body size, and native vs. exotic status of each species. For each non‐parasitic species that represented at least 0.01% of the total abundance in the full 2002–2016 BIML dataset (n = 173 species), we summarized trait information from existing literature (Bartomeus et al., 2013b; Danforth, 2015; Michener, 2007). We were not able to identify the nesting location of some species (n = 39), so removed them from trait‐based analyses, leaving only species with known nesting locations (n = 134). We extracted intertegular distance of each species from published datasets (Bartomeus et al., 2013b; Kammerer et al., 2016) and calculated body mass from bee intertegular distance with the pollimetry package 1.0.1 in R (Kendall et al., 2019; R Core Team, 2019). These data are available in Table S3 and the associated Dryad archive (see “Data Availability Statement”).

Landscape characteristics

To quantify landscape characteristics surrounding each bee‐sampling site, we calculated landscape composition and landscape indices of floral resources, nesting resources, and insecticide toxic load. For all landscape metrics, we used land‐cover data from the Cropland Data Layer (CDL) generated by the United States Department of Agriculture, National Agricultural Statistics Service (USDA NASS, 2016).

Landscape composition

We calculated landscape composition surrounding each site by grouping CDL classes into “developed,” “natural,” and “arable” categories (Table S1). We calculated proportional area of each land‐cover category at 250, 500, and 1750 m radii from the center of our sampling site, corresponding, respectively, to the mean foraging range for small‐bodied species, all species, and large‐bodied species in our region (Kammerer, Biddinger, Joshi, et al., 2016). In preliminary analyses, we found that landscape composition at 250‐ and 500‐m radii were highly correlated with 1750‐m radius (Figure S5), so used a 1750‐m radius for all analyses presented here. When possible, we calculated landscape indices for each site using land‐cover data matching the year of bee sampling. However, the CDL was not available for our study area from 2003 to 2007, so for sites sampled in this period, we used land‐cover data from 2008.

Nesting and floral‐resource indices

We generated floral‐ and nesting‐quality indices using the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) crop pollination model, version 3.5.0 (Lonsdorf et al., 2009; Natural Capital Project, 2018). To calculate landscape‐scale availability of floral and nesting resources, the InVEST model requires estimates of nesting and floral quality for each land‐cover class. We adopted estimates of nesting and floral quality from Koh et al. (2016), who conducted an expert‐opinion survey and validated their survey results with field data from over 500 sites across the United States. We assumed that “Background,” “Clouds/NoData,” and “Undefined/NonAg” CDL classes did not have floral and nesting resources, nor any insecticide toxic load (see description below).

To calculate foraging resources and nesting quality for a specific site, the InVEST model requires a bee‐community table specifying phenology, nesting location, and foraging range of taxa included in the model (Tables S3 and S4). We assembled these parameters from BIML and published data. For each wild‐bee species (n = 45) that represented at least 0.5% of total abundance, we calculated the proportion of specimens captured in each season (spring, summer, and fall). For rarer species, we did not have enough observations to characterize species‐level phenology. Based on our genus‐level phenology (Figure 1), we defined spring as before mid‐May, summer as mid‐May to mid‐August, and fall as after mid‐August. We had to define “fall” as a separate season for InVEST calculations because Koh et al.’s (2016) resource coefficients include fall‐specific values for habitat floral resource. We utilized the information on nesting location and intertegular span described above and estimated foraging ranges from intertegular span using a published equation (Greenleaf et al., 2007). Then, to reduce the number of taxa for computational efficiency, we calculated an abundance‐weighted average of seasonal activity and foraging range of all species sharing a nesting location (ground, stem, wood, cavity/ground, and cavity/stem, Table S5).

Using land‐cover maps (1750‐m radius CDL), resource values for each land‐cover type, and the bee‐community table, we ran the InVEST model for each site‐year. To formulate indices of total floral and nesting resources, we averaged InVEST resource maps over all taxa, weighting according to the relative abundance of each taxon. For nesting resources, prior to averaging we also applied a distance‐weighting function based on the foraging range of each taxa, as in Lonsdorf et al. (2009). We used an exponential decline function to weight nesting contribution of each cell as a function of its distance from the focal cell (sampling site). Then, we summed all pixels within two times the foraging range to calculate total accessible nesting resources at each BIML site. Floral resource maps were already distance‐weighted within InVEST calculations, so we did not apply an additional distance‐weighting procedure.

Insecticide toxic‐load index

We quantified insecticide toxic load at each BIML sampling site with an insecticide index presented by Douglas et al. (2020), Douglas et al. (2020). This index was derived from state‐specific data on identity and quantity of insecticides per ha applied to agricultural crops, scaled by toxicity to the honey bee (Apis mellifera) to estimate bee lethal doses per ha. For this analysis, we used the Douglas, Soba, et al. (2020), Douglas, Sponsler, et al. (2020) insecticide toxic load from each CDL land‐cover class, with separate indices for each state in the BIML dataset. Necessary state‐level insecticide data do not exist for Washington DC, so for BIML sites in Washington, DC, we used insecticide data from comparable land‐cover classes in Maryland. For all BIML sites (sampled from 2003–2013), we used insecticide data matching the year of bee sampling.

Our implementation of the Douglas, Soba, et al. (2020), Douglas, Sponsler, et al. (2020) index integrated toxic load of insecticides at sampling locations and in the surrounding landscape. For each sampling location in the BIML dataset, we converted CDL land cover to insecticide toxic load using Douglas, Soba, et al's (Douglas, Soba, et al., 2020) reclassification tables. Then, as for nesting resources, we applied a distance‐weighting function and calculated total insecticide toxic load at the sampling location. To distance‐weight insecticide index maps, we used the mean foraging range, weighted by relative abundance, of all 45 species in the bee community table (see Section 2.4.2).

Climate and topography

We represented annual climate by generating a standard set of bioclimatic (BIOCLIM) variables used in species distribution models from 4 km gridded PRISM daily temperature and precipitation data (Busby, 1991; PRISM Climate Group, 2019). The BIOCLIM indices include monthly maximum, minimum, and mean, annual total, and seasonality of temperature and precipitation variables, among others (Table 1). We calculated these variables for the year of wild‐bee sampling (denoted “Year0” in the figures and tables), one year before sampling (denoted “Year‐1”), and a 15‐year climate normal (2001–2015, denoted “15Year”).

TABLE 1

Predictor variables included in random‐forest analyses. All climate and weather variables were calculated from the 4 km PRISM daily data(PRISM Climate Group, 2019) for the year of bee sampling, 1 year prior to sampling, and 15‐year means (2001–2015; U.S. Geological Survey, 2014; USDA Douglas, Soba, et al., 2020, Douglas, Sponsler, et al., 2020; Koh et al., 2016; NASS, 2016)

Predictor variables	Unit	Source
Climate/Weather
Annual mean temperature	°C	BIOCLIM a
Mean diurnal range	°C	BIOCLIM a
Isothermality (temperature evenness)	percent	BIOCLIM a
Temperature seasonality (standard deviation × 100)	°C	BIOCLIM a
Maximum temperature of warmest month	°C	BIOCLIM a
Minimum temperature of coldest month	°C	BIOCLIM a
Temperature annual range	°C	BIOCLIM a
Mean temperature (wettest, driest, warmest, and coldest quarter)	°C	BIOCLIM a
Growing degree days (base 4.5 and 10°C)	°C	USDA NRCS a
Annual maximum and minimum temperature	°C	USDA NRCS a
Length of frost‐free growing season (−2.2 and 0°C threshold)	days	USDA NRCS a
Date of last spring frost (−2.2 and 0°C threshold)	date	USDA NRCS a
Date of first fall frost (−2.2 and 0°C threshold)	date	USDA NRCS a
Mean temperature in July	°C	USDA NRCS a
Monthly temperature range	°C	USDA NRCS a
Annual precipitation	mm	BIOCLIM a
Precipitation of wettest and driest month	mm	BIOCLIM a
Precipitation seasonality (coefficient of variation)	mm	BIOCLIM a
Precipitation (wettest, driest, warmest, and coldest quarter)	mm	BIOCLIM a
Growing season total precipitation	mm	USDA NRCS a
Mean number of days between 0.254 cm or greater rain events	days	USDA NRCS a
Mean number of days between 0.76 cm or greater rain events	days	USDA NRCS a
Monthly precipitation range	cm	USDA NRCS a
Topography
Elevation	m	USGS b
Slope	percent	USGS b
Winter and summer solar radiation	watt hours/m 2	USGS b
Topographic convergence index	NA	USGS b
Aspect (degrees NS, EW)	degrees	USGS b
Landscape quality
Percent developed land within 1750 m radius	percent	USDA NASS c
Percent agriculture within 1750 m radius	percent	USDA NASS c
Percent natural land within 1750 m radius	percent	USDA NASS c
Floral resource quality (spring, summer, and fall)	NA	Koh et al d
Nesting resource quality (for spring, summer, and fall bees)	NA	Koh et al d
Combined resource quality (spring, summer, and fall)	NA	Koh et al d
Local habitat type	NA	USDA NASS c
Insecticide toxic load		Douglas, Soba, et al. (2020), Douglas, Sponsler, et al. (2020)
Abundance of Apis mellifera in pan trap samples	number of individuals	BIML data
Sampling method
Volume and color of pan traps	NA	BIML data

PRISM Climate Group. http://prism.oregonstate.edu (2019).

U.S. Geological Survey. National Elevation Dataset. https://gdg.sc.egov.usda.gov/ (2014).

USDA NASS. USDA National Agricultural Statistics Service Cropland Data Layer. http://nassgeodata.gmu.edu/CropScape/ (2016).

Koh et al. Modeling the status, trends, and impacts of wild bee abundance in the United States. Proc. Natl. Acad. Sci. 113, 140–145 (2016).

In our study region, mean foraging distance of furthest flying, large‐bodied bees is approximately 1750 m (Kammerer, Biddinger, Joshi, et al., 2016), so we assumed that virtually all species would respond to variation in temperature and precipitation at a smaller spatial extent than the 4 km PRISM data (Wang et al., 2017). To account for small‐scale variability, we calculated several topographic variables including northerly and westerly components of aspect, slope, and curvature (Table 1). We derived all topographic variables from 30‐m gridded elevation data (U.S. Geological Survey, 2014) using GRASS GIS software (GRASS Development Team, 2018).

Species‐richness rarefaction

To compare species richness at sites sampled with uneven effort, we adjusted observed species richness with coverage‐based rarefaction and extrapolation using the iNEXT package in R (Hsieh et al., 2016, 2018; R Core Team, 2019). Coverage‐based rarefaction and extrapolation is preferred over traditional rarefaction methods because it adjusts species richness based on completeness of the sample (coverage), rather than sample size, which reduces bias when comparing low‐ and high‐diversity communities (Chao & Jost, 2012). We estimated species richness at the mean coverage level of all site‐years in each random‐forest dataset (see below). Using sample‐coverage‐based rarefaction, species richness can be robustly estimated at sample sizes up to twice the observed number of occurrences (Chao & Jost, 2012; Hsieh et al., 2016). We removed site‐years that exceeded this threshold, as resulting species‐richness estimates are subject to large prediction bias.

Statistical analyses

For each of (1) spring and summer/fall total abundance and richness, (2) spring genus‐specific abundance, and (3) summer/fall functional group abundance, we examined wild‐bee responses to land use (local habitat type, landscape composition, and resource indices), climate, and topography (Table 1) with random‐forest machine‐learning models. We selected functional groups comparisons for summer/fall bees based on the results of a fourth‐corner analysis (see Supplemental Materials and Figure S6). Specifically, for summer/fall bees, we fit random‐forest models for ground versus cavity/stem nesting, eusocial versus solitary species, and native versus non‐native species. For spring bees, we examined the three most abundant genera in lieu of trait‐based groups as the fourth‐corner analysis did not detect strong trait–environment interactions. Our trait groups varied in mean abundance, so to quantitatively compare between trait groups, we used z‐score normalized bee abundance day−1 trap−1or species‐richness response variables (set to a mean of zero and unit‐variance) before tuning and executing random‐forest analyses. To evaluate the performance of specific random‐forest models, we also fit models with untransformed abundance or richness data and include these results in supplemental materials (Table S6).

We conducted all random‐forest analyses in R 3.6.0 (R Core Team, 2019). We used the caret package 6.0–84 to construct random regression forests with the “ranger” function (Kuhn, 2008, 2019; Wright & Ziegler, 2017). For each random‐forest model, we selected the optimal number of trees (1000–5000 trees incremented by 500) and the optimal number of variables included at each tree split (25–55 variables incremented by 5) with a grid search. We defined the optimal model as the model that minimizes root mean‐squared error. We employed 10‐fold cross validation performed three times to assess model performance. To interpret random‐forest results, we calculated permutation‐based variable importance scores with ranger package 0.11.2. Below, we present variable importance results for all variables which scored 80 or higher (out of a maximum of 100) in any random‐forest model (Figure 2). Variable importance scores for all variables are available in supplemental materials (Table S7). We also generated accumulated local‐effects plots with iml package 0.9.0 (Molnar et al., 2018; Wright & Ziegler, 2017). We present accumulated local‐effects plots instead of partial dependence plots as they are more robust for datasets with correlated predictors (Molnar, 2019).

FIGURE 2

Relative importance of climate, topography, landscape quality and land‐use variables in predicting season‐total abundance or richness (left panel), genus (middle panel), or trait‐specific (right panel) abundance of wild bees. Dark blue indicates the most important variables for each random‐forest analysis. Climate and weather variables were calculated from the year of bee sampling (‘Year0’), 1 year before sampling (‘Year‐1’), and 15‐year climate normal for 2001–2016 (‘15Year’)

RESULTS

Wild‐bee community characterized by long‐term monitoring data

After all cleaning and quality checks, the BIML dataset we used for richness analyses included 65,949 occurrences from 917 sites and 976 site‐years. The abundance dataset (occurrences with recorded sampling effort) contained 54,694 occurrences from 895 sites and 952 site‐years. The richness and abundance datasets were similarly speciose, with 299 and 288 species, respectively, which represent 40 of 169 bee genera found in North and Central America (Michener et al., 1994).

Of the 134 most abundant, non‐parasitic species of wild bees in the BIML dataset (Table S3), 70% excavate nests in the ground, 19.5% nest in pithy stems of plants or other pre‐existing cavities, 6% nest in pre‐existing cavities in or near the ground, like rodent burrows, under organic material, or rock crevices, and 4.5% excavate nests in dead or rotting wood. Just over half of the wild‐bee species in this dataset are solitary (60%), 30% are eusocial, and 10% facultatively or primitively social. Only six of the 134 most abundant, non‐parasitic species are not native to North America, Andrena wilkella (Kirby 1802), Anthidium manicatum (Linnaeus, 1758), Anthidium oblongatum (Illiger, 1806), Megachile rotundata (Fabricius, 1787), Osmia cornifrons (Radoszkowski, 1887), and O. taurus (Smith, 1873), five of which are cavity or stem nesting species from the Megachilidae family. We also observed substantial differences in phenology among wild‐bee genera represented in the BIML dataset (Figure 1). Three genera (Andrena, Osmia, and Nomada) were highly abundant early in the season, but much less common after mid to late May.

Seasonal richness and abundance of wild bees

Weather conditions in the year prior to sampling and 15‐year average climate were better predictors of seasonal abundance and richness of wild bees than landscape composition, landscape indices, or topography (Figure 2). Specifically, for spring bee communities, we were able to explain 25% and 17% of the variation in wild‐bee richness and abundance, respectively (Table 2). The maximum temperature in the warmest month (July) and annual maximum temperature in the year before sampling were the best predictors of species richness, with a positive relationship between temperature and richness (Figures 2 and 3a,b). In contrast, for bee abundance, trap color was the most important variable, and temperature was not predictive (Figures 2 and 3a,b).

TABLE 2

Random‐forest model performance predicting abundance or richness of wild bees. Response variables were z‐score normalized prior to analysis. The R‐squared, root mean squared error (RMSE), and mean absolute error (MAE) were calculated using 10‐fold cross validation performed three times and are presented as mean ± standard deviation. Parameters mtry and ntree are the optimal values (for number of variables and number of trees, respectively) for each random‐forest model as determined by grid‐search parameter tuning

Wild‐bee taxa	Season	Response variable	RMSE	R 2	MAE	mtry	ntrees
All species	Spring	Abundance day−1 trap−1	0.89 ± 0.33	0.17 ± 0.10	0.59 ± 0.13	15	5000
	Summer/fall		0.90 ± 0.32	0.13 ± 0.09	0.59 ± 0.09	25	1000
	Spring	Richness	0.88 ± 0.16	0.25 ± 0.17	0.67 ± 0.09	35	5000
	Summer/fall		0.86 ± 0.08	0.27 ± 0.08	0.66 ± 0.05	35	2000
Andrena	Spring	Abundance day−1 trap−1	0.76 ± 0.56	0.28 ± 0.20	0.37 ± 0.14	35	2000
Nomada			0.93 ± 0.28	0.10 ± 0.10	0.59 ± 0.12	15	1000
Osmia			0.89 ± 0.28	0.21 ± 0.15	0.54 ± 0.11	15	5000
Cavity/stem nesters	Summer/fall		0.85 ± 0.36	0.23 ± 0.15	0.45 ± 0.09	55	2000
Ground nesters			0.87 ± 0.40	0.17 ± 0.11	0.53 ± 0.10	50	3000
Eusocial species			0.89 ± 0.15	0.22 ± 0.10	0.57 ± 0.06	40	3000
Solitary species			0.89 ± 0.27	0.17 ± 0.08	0.51 ± 0.09	35	5000
Native species			0.90 ± 0.33	0.13 ± 0.08	0.57 ± 0.08	30	5000
Non‐native species			0.83 ± 0.40	0.23 ± 0.17	0.42 ± 0.09	55	3000

FIGURE 3

Relationships between the most important landscape and climate variables and seasonal abundance and richness of wild‐bee communities. The abundance values shown are z‐score normalized (mean of zero and unit‐variance) to enable comparisons between different random‐forest models. We defined “spring” as prior to May 15th and “summer/fall” as after May 15th, based on the phenology of wild‐bees in our region (Figure 1)

For summer/fall bee communities, our models explained 27% and 13% of the variation in richness and abundance, respectively (Table 2). Mean temperature in the warmest quarter (i.e., June–August) in the year before sampling was the strongest predictor of wild‐bee abundance (Figure 2). In contrast to the effects of temperature on spring bees, abundance of bees in summer and fall decreased with warmer temperatures the previous year; specifically, when mean temperature the previous summer exceeded about 24°C (Figure 3c,d). For wild‐bee richness in the summer and fall, mean number of days between precipitation events and season‐total precipitation were the best predictors (Figure 2). If the average time between precipitation events was greater than about 6.5 days, the number of summer bee species decreased dramatically (Figure 3c,d).

Genus‐specific responses of spring‐active bees

Analyzing subsets of the spring bee community revealed notable differences among genera. Andrena abundance was strongly predicted by mean temperature of the winter immediately preceding bee sampling, with fewer bees when average winter temperature was greater than about 1.5°C (Figures 2 and 4a,b). Like Andrena, Osmia abundance was lower with warmer winter temperatures but also dropped substantially when annual precipitation exceeded about 875 mm (Figures 2 and 4a,b). In summer and fall, we estimated higher abundance of solitary bees above 875 mm of annual precipitation, but for spring Osmia we found the opposite trend. Our random‐forest models for Andrena and Osmia explained approximately 4%–11% more of the variation in spring‐bee abundance than the season‐total model (28% and 21%, respectively, Table 2), but the model for Nomada performed worse than the season‐total model, so we chose not to interpret the results for Nomada.

FIGURE 4

Relationships between the most important landscape and climate variables and abundance of contrasting wild‐bee genera and functional groups. The abundance values shown are z‐score normalized (mean of zero and unit‐variance) to enable comparisons between different random‐forest models. To compare Andrena and Osmia (panels a and b), we used bee occurrences sampled in the spring (prior to May 15th). For native vs. non‐native (panels c and d) and eusocial vs solitary (panels e and f) comparisons, we utilized bees collected in the summer and fall (after May 15th)

Functional group responses of summer/fall active bees

In the summer/fall bee community, we found several notable contrasts between functional groups differing in native status and sociality (Figure 2). We documented a contrast between native and exotic species, with a strong positive relationship between exotic species abundance and percent‐developed land; this relationship was not apparent for native species. Estimated abundance of exotic species (summer/fall) increased four times across a gradient from 25% to 75% developed land cover (Figure 4, middle). Exotic‐bee abundance was constant when land use was less than 25% or greater than 75% developed land. In contrast, native‐bee abundance did not change with amount of developed land in the landscape but decreased with mean summer temperatures greater than approximately 23°C (Figure 4c,d).

Solitary species in summer and fall strongly differed from eusocial species in their response to land use and weather (Figure 4e,f). We found a threshold in the 15‐year season‐total precipitation, at approximately 875 mm, above which we estimated higher abundance of solitary bees (Figure 4e,f). Eusocial species had no relationship with annual precipitation. Abundance of eusocial species declined with more natural land surrounding sampling sites while solitary species had a very weak, positive response to the proportion of natural land around sites (Figure 4e,f).

Explaining approximately 17%–23% of the variation in wild‐bee abundance, the random‐forest models based on functional groups of wild bees performed better than the model for all taxa, except the model for native species (R 2 = 0.13), which was comparable to the all‐taxa model (Table 2). Almost all species in the BIML dataset were native, so it is not surprising that the native‐only model was very similar to the model for all species. Generally, compared to functional groups with fewer species (cavity/stem nesting, eusocial, and exotic species), we observed higher unexplained variation when modelling abundance of functional groups with many species, like ground‐nesting, solitary or native species (Table 2). This indicates that these large functional groups likely obscured relevant species or genera‐specific responses to landscape or climate that might be revealed by modelling narrower taxonomic or functional groups.

DISCUSSION

Despite a large body of literature documenting land‐use effects on wild bees (Kennedy et al., 2013; Ricketts et al., 2008), we found that climate factors were the main drivers of wild‐bee abundance and richness, although total variance explained was low across all models. Almost all our models in both seasons (spring or summer/fall) indicated that temperature and precipitation were more important predictors of wild‐bee communities than landscape composition, landscape quality, or topography; however, relationships varied substantially between wild‐bee genera and functional groups. Generally, our results are consistent with previous research showing bee abundance fluctuates in response to varying temperature and precipitation (Papanikolaou et al., 2017b). We did not find that the negative effect of temperature was offset by semi‐natural habitat in the landscape (Papanikolaou et al., 2017a), but the landscapes surrounding our study sites contained substantially more semi‐natural habitat than the previous study. The lack of correlation of landscape composition and resource quality with wild bee abundance and diversity is surprising, but it is likely that our measures of these factors (using general land‐use metrics or the floral/nesting index) did not accurately reflect local conditions, and/or flowering‐plant communities are strongly influenced by climate as well. Remote sensing approaches may improve the predictive abilities of floral and nesting metrics by more accurately measuring vegetation characteristics at local and landscape scales (Galbraith et al., 2015).

For spring bees (such as Osmia), precipitation was the most important driver of abundance, with higher precipitation associated with lower abundance. Increased spring precipitation may prevent females from finding suitable nest sites or, for individuals with established nests, reduce the number of suitable foraging days (Drummond et al., 2017; Tuell & Isaacs, 2010). Some studies report Osmia can forage in rainy weather (Vicens & Bosch, 2000), but others show decreased or no foraging during overcast or rainy conditions (Drummond & Stubbs, 1997; McKinney & Park, 2012). Individuals that continue to forage in rainy conditions are more likely to return with smaller pollen loads, take longer to complete foraging bouts, or fail to return to their nests (Tuell & Isaacs, 2010), likely due to unfavorable flight energetics or reduced sensory abilities (Lawson & Rands, 2019). Additionally, precipitation influences spring plant‐community composition and floral‐resource availability. Finally, increased precipitation could reduce nesting success, as found in recent study on Anthidium vigintipunctatum, a cavity‐nesting species (Megachilidae; Vitale et al., 2020). In all of these scenarios, reproductive output would be reduced, resulting in lower abundance the following year. Comparative behavioral studies would be extremely valuable to understand why increased precipitation is more detrimental to Osmia compared with Andrena species. From available literature and our data, we found no clear explanation for this trend.

For summer/fall bee community, temperature was the main factor dictating bee abundance, with lower abundance, particularly of native species, following years with hot summers. We observed lower abundance of bees in summer/fall when mean temperature the previous summer exceeded about 24°C, which may represent summer temperatures too far above bees’ thermal optima of 24.9 ± 4.7°C (Kühsel & Blüthgen, 2015). With high temperatures, bees risk overheating during flight, so commonly cease flight for short periods, seek shade and cooler microclimates, or elevate convective cooling by flying faster (Corbet & Huang, 2016; Willmer & Stone, 2004). These behaviors are effective in mitigating heat risk to the individual, but reduce amounts of pollen and nectar collected for nest provisions, likely leading to lower reproductive output.

Notably, eusocial and solitary species had different responses to increased precipitation. In addition to fewer species of all summer/fall bees with longer temporal gaps between precipitation events, we documented fewer solitary bees at locations with less rain. We hypothesize that summer/fall solitary species may be more sensitive to drought conditions than eusocial species. In studies of one solitary species, Andrena vaga, individuals began foraging earlier and had a shorter lifespan in warm, dry years than during cold, wet years, which may be due to heat stress of foraging females (Straka et al., 2014). The pattern we observed of more summer/fall solitary bees at locations with higher precipitation contrasts with our finding of fewer spring bees at rainier locations. This contrast may be due to the typical frequency of precipitation in each season. In our study region, on average, spring months (March–May) have 11 days with measurable precipitation (greater than 0.0254 cm) each month while summer and fall (June–October) months have 9 days (NOAA National Centers for Environmental Information, 2018). Thus, spring bees, with their shorter phenological window, may already be foraging‐limited by precipitation, which is then further limited by increased precipitation.

For summer/fall eusocial species, above approximately 50% natural land, we observed a decline in the number of bees, which might be related to the composition and mobility of eusocial versus solitary‐bee communities. In summer/fall, 34 out of 40 eusocial species and 96% of eusocial individuals were in the family Halictidae, mostly small‐bodied Lasioglossum species, with reduced foraging ranges (132 m, weighted by relative abundance) relative to the solitary species (344 m). In our study region, forest is the primary type of natural land (USDA NASS, 2016), and in summer provides few floral resources. In landscapes with high forest cover, small‐bodied eusocial species might not be able to access other habitats to collect adequate pollen and nectar resources. Of the five most‐abundant, eusocial species in our dataset, three (Lasioglossum pilosum, Lasioglossum tegulare, and Lasioglossum bruneri) have been reported elsewhere to avoid forest habitat, two were neutral toward forest habitat, and none preferred forests (Collado et al., 2019).

We documented a strong pattern of higher abundance of exotic bees in developed areas and predict that the prevalence of exotic bees in our study region will continue to increase. Studies in other cities in eastern US and Canada found similar results (Matteson et al., 2008; Normandin et al., 2017), while crops and natural habitat in Eastern Canada supported only five exotic species (Grixti & Packer, 2006; Sheffield et al., 2003). Most of the urban exotic species in our study region are cavity nesters (see Section 3.1 and Table S3), and increased abundance of exotic bees in cities has not been linked to proximity to ports, urban warming, or abundance of exotic floral resources, leaving increased availability of cavity nest sites as a likely driver (Fitch et al., 2019). Between 2001 and 2016 in our study region, total impervious cover associated with developed land increased by 13.4%, which represents an additional 18,426 ha of impermeable surfaces (Yang et al., 2018). Impervious cover in some counties increased by as much as 30.5% (2180 ha; Yang et al., 2018). In urbanizing landscapes like our study region, exotic‐bee abundance will almost certainly continue to increase, with unknown consequences for native‐bee communities (Fitch et al., 2019; Russo, 2016).

Warm winters were challenging for several key taxa of wild bees, specifically spring bees that overwinter as adults. Both Andrena and Osmia overwinter as adults, and were less numerous following years with warmer winters. This result coincides with studies showing that, with earlier spring onset, taxa overwintering as adults had higher pre‐emergence weight loss and mortality and shorter life span post‐emergence (Fründ et al., 2013; Schenk et al., 2018; Slominski & Burkle, 2019). Conversely, species that overwinter as pre‐pupae were more responsive to spring temperature (Slominski & Burkle, 2019) and microclimate, which mediates nest‐site selection (Wilson et al., 2020). This pattern of reduced Andrena and Osmia abundance with warmer winters is concerning for several reasons. Both genera are important pollinators of spring‐blooming crops (Adamson et al., 2012; Park et al., 2016; Russo et al., 2017). In Northeastern USA, winter temperatures have warmed three times faster than summer (Thibeault & Seth, 2014), advancing shrub leaf‐out and blooming by 1.2–1.6 days per decade (Ault et al., 2015). If Andrena and Osmia populations are threatened by increasingly warm winters, it could force management changes in some tree‐fruit cropping systems that rely heavily on these taxa (Blitzer et al., 2016; Russo et al., 2017). Additionally, in the future, warmer winters may coincide with wetter springs (Easterling et al., 2017; Lynch et al., 2016; Thibeault & Seth, 2014), which would pose a particular challenge to Osmia and may jeopardize efforts to develop species in this genus as managed pollinators for tree‐fruit production in the Eastern USA (Sedivy & Dorn, 2014).

There are some important limitations to this study. We did not consider the effect of landscape configuration, because, for wild bees generally composition of landscapes is more important than configuration (Holzschuh et al., 2010; Kennedy et al., 2013; Steckel et al., 2014). But, some wild‐bee species are affected by configuration of the surrounding landscape (Hass et al., 2018; Hopfenmüller et al., 2014; Martin et al., 2019), so future work should consider both landscape composition and configuration interactions with climate. All of the specimens in the BIML dataset we used were passively collected with pan traps. Pan traps are known to under‐sample large‐bodied bee taxa, like Bombus and Xylocopa (Joshi et al., 2015; Roulston et al., 2007), so these genera are not well represented in our analyses. Additionally, sampling locations were not consistent throughout the study period, with many sites sampled fewer than 3 years. This design increased the number of sites and habitat types included but also likely increased intra‐annual variation in bee abundance and richness. The sampling method employed by the BIML also changed over time. We addressed this by selecting a subset of the BIML data with more consistent sampling methods and including all available sampling method information (trap color and volume) in our analyses. Differences in sampling locations and method likely contributed to the substantial variation in bee abundance and richness unexplained by our models (Table 2). Overall, our analysis of this large, observational dataset cannot determine specific mechanisms of land‐use and climate effects on wild bees. Rather, we report the relative importance of climate and land‐use drivers, and present wild‐bee responses to specific weather and landscape variables to generate hypotheses to be tested with more targeted field studies and experimental research.

CONCLUSION

Here, we showed that climate patterns, specifically spring precipitation and average summer and winter temperatures, were key drivers of interannual and spatial variation in abundance of wild bees in temperate ecosystems. In the Northeast USA, past trends and future predictions show a changing climate with warmer winters, more intense precipitation in winter and spring, and longer growing seasons with higher maximum temperatures (Easterling et al., 2017; Lynch et al., 2016; Thibeault & Seth, 2014). In almost all our analyses, these conditions were associated with lower abundance of wild bees. Wild‐bee richness results were more mixed, including neutral and positive relationships with temperature and precipitation patterns predicted to increase in the future. In the Northeast USA, combined with continued urbanization, changing climate imposes a significant threat to wild‐bee communities.

The relationship between climate conditions and wild‐bee abundance and richness deserves more research attention. We especially recommend research to elucidate the mechanisms underlying these variable relationships and implications for fluctuating wild‐bee abundance for pollination service provisioning. A more mechanistic understanding of direct and indirect effects of temperature and precipitation on wild bees, and how these interact with land use, is crucial to inform climate‐resilient conservation of bee populations. By including climate variables, landscape pollination models and decision‐support tools would likely more accurately predict interannual variation in wild‐bee abundance and effects on pollination services for crops and wild‐plant communities.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Supplementary Material

Click here for additional data file.