Bee community response to local and landscape factors along an urban-rural gradient

Over a third of the world’s crops require insect pollination, and reliance on pollination services for food continues to rise as human populations increase. Furthermore, as interest in urban agriculture has grown, so has a need for studies of urban pollinator ecology and pollination. Analyzing pollinator assemblages along a rural-urban gradient provides powerful mechanistic insight into how urbanization impacts pollinators. Yet, studies examining pollinators along urban-rural gradients are limited and results vary. Since pollinators vary tremendously in life history characteristics and respond to urbanization differently, studies from different regions would improve our understanding of pollinator response to urbanization. This study documents different bee assemblages along a high-plains semi-arid urban-rural gradient in Denver, Colorado, USA. Percent impervious surface was used to define the extent of urbanization at 12 sites and local and landscape characteristics were estimated using field assessments and geospatial analysis. Wild bees were collected and the relationships between urbanization and bee communities were explored using linear modeling. Overall, bee abundance and diversity decreased with increasing urbanization, suggesting that urban areas negatively impact bee communities. However, all bee guilds responded positively to local floral richness and negatively to the degree of landscape urbanization, suggesting that different types of bees responded similarly to urbanization. These findings suggest that providing a greater diversity of floral resources is key to mitigating the negative impacts of urbanization on pollinator communities.


Introduction
Pollination services provided by bees and other insects is fundamental to the production and regeneration of wild plants and agricultural crops. Over a third of the world's cropsincluding fruits, vegetables, nuts, spices, and oilseed-require insect pollination (Klein et al. 2007), and our reliance on pollination services to promote these crops continues to rise due to increasing demands from growing human populations.
Between 1961 and 2006, dependence on animal pollination increased by 50% in developing nations and by 62% in developed countries (Aizen and Harder 2009). As of 2010, insect pollination was estimated at US$212 billion globally, representing 9.5% of the total value of agricultural production (van Engelsdorp and Meixner 2010), and in the US alone, bee pollination is responsible for $14 billion of agricultural production (Morse and Calderone 2000;Cane 2005). Furthermore, as urban populations expand, community gardening activities increase, demanding greater pollination services (Frail 2010;Palmer 2018). Demand for pollinatordependent crops has tripled during the last half century and pollinator losses could have serious implications for plant reproduction, animal populations that depend on those plants, and also world food security (Majewski 2016;Novais et al. 2016;Aizen and Harder 2009). Thus, understanding the biogeography of insect pollinators in a variety of environments (e.g., urban, rural, and agricultural) will become increasingly important as pollination needs increase.
Recently, bees have received substantial attention from the media and scientific community due to global population Electronic supplementary material The online version of this article (https://doi.org/10.1007/s11252-020-00956-w) contains supplementary material, which is available to authorized users. declines. Particularly, honeybee declines are well documented in North America and Europe, with 49.5% managed honeybee colony losses in North America and 25.5% colony losses in Europe between 1961(FAOSTAT 2009). Moreover, according to the Xerces Society (2017), 57 different wild bee species have been identified as endangered, threatened, or atrisk in North America, with over 30 bee species classified as critically imperiled or possibly extinct. These declines are due to a range of factors, including land-use change and intensification, climate change, pesticide applications, lack of food and nesting resources, and the spread of exotic species and diseases (Vanbergen and the Insect Pollinators Initiative 2013;Potts et al. 2010;Potts et al. 2016;Grunewald 2010). Interactions between these pressures likely drive much of the decline in pollinators (Vanbergen and the Insect Pollinators Initiative 2013). For example, converting natural land to agricultural fields reduces food and nesting resources and increases exposure to pesticides. These impacts result in declining nutrition, which may increase a bee's susceptibility to pathogens, thus decreasing foraging activity, reducing homing capabilities, crippling reproductive efforts, and curtailing its lifespan (Vanbergen and the Insect Pollinators Initiative 2013;Gill et al. 2012;Feltham et al. 2014;Gill and Raine 2014;Henry et al. 2012;Fischer et al. 2014;Whitehorn et al. 2012;Rundlof et al. 2015).
Bees and other pollinators are significantly impacted by anthropogenic disturbances such as land use modification (Ahrne et al. 2009;Leong et al. 2014;Theodorou et al. 2017) as evidenced in the extensive literature describing bees and other pollinator assemblages in both rural and urban contexts. However, less is known about variation in the degree to which bees are affected by these changes and research examining pollinator assemblages along spatial continua, such as urban-rural gradients, is limited and results vary. Some scholars have found that certain pollinators, such as bees, are positively associated with urbanization (Martins et al. 2017;Theodorou et al. 2017). Others have found that bee communities decline with increasing urbanization (Geslin et al. 2016;Lagucki et al. 2017;Ahrne et al. 2009;Bates et al. 2011;Choate et al. 2018;Verboven et al. 2014). Variation in these findings may be due to differences in pollinator assemblages or regional variation in bee communities. For example, specific bee guilds, or a species subgroup that achieves similar functions in the community and has similar resource requirements (de M Santos et al. 2013), may respond differently to increasing urbanization (Banaszak-Cibicka and Zmihorski 2012;Fortel et al. 2014). Banaszak-Cibicka and Zmihorski (2012) examined the differences in bee species richness and community composition along an urban-peri-urban gradient in Poznan, Poland and found urban landscapes can "filter" certain bee guilds according to ecological traits, favoring small-bodied eusocial bees that begin activity later in the season. Another study from Grand Lyon, France found that cavity-nesting bees and long-tongued bees occurred in areas with intermediate to high proportions of impervious surfaces (Fortel et al. 2014). They found no strong effects of urbanization on the occurrence of species according to body size or social behavior. Similarly, others have found that cavity nesting may be an important ecological trait governing bee prevalence in urban environments (Pardee and Philpott 2014;Tonietto et al. 2011;Fitch et al. 2019). Given geographic differences in landscapes, habitats, and urban practices, guild-specific responses of bee communities to urbanization likely vary across different regions or biomes. Thus, exploring bee community response to urbanization in more regions, and especially in under surveyed areas, is critical to building a better understanding of how urbanization impacts bees and the critical ecosystem services they provide.
The objectives of this study were to understand how pollinators are influenced by urbanization in an understudied highplains semi-arid metropolitan region. The Denver metropolitan area of Colorado is one of the fastest growing urban regions in the U.S. and by examining bee abundance and diversity along the city's urban-rural gradient, we may determine if bee species with different life history characteristics are more resilient to urbanization, and how to best mitigate the impacts of urbanization for more susceptible bees. Specifically, we asked: 1) how do bee abundance and diversity respond to urbanization, 2) how does urbanization impact bees from different foraging, nesting, and social guilds, and 3) what mechanisms drive bee diversity and community composition across the urban-rural gradient? We hypothesized that bee diversity and abundance in the Denver region would be greatest in rural settings, intermediate in suburban areas, and lowest in urban landscapes. This hypothesis was based on research suggesting that while urban areas have the capacity to support certain bee species, bee diversity and abundance diminish with increasing urban intensity (Bates et al. 2011). Second, we hypothesized that bee communities in areas with greater urbanization would be dominated by smaller bees with more ecologically generalized characteristics (e.g., polylectic and eusocial) than rural areas which likely harbor bee communities with more specialized traits (e.g., oligolectic). Finally, we hypothesized that bee abundance and diversity would respond negatively to increasing impervious surfaces and positively to increasing floral resources.

Study area
We conducted this study at 12 sites across the Denver metropolitan area (39.7392°N, 104.9903°W). The area is relatively high elevation (1610 m above sea level) and in a semi-arid high-plains environment 19 km east of the Rocky Mountains.
The estimated population of over 2.7 million residents as of 2017 makes it the 19th most populous metropolitan area in the United States (U.S. Census Bureau Population Division 2017). Denver's 1.4% growth rate is ranked seventh in the nation, adding nearly 80,000 people between 2017 and 2018 (Murray 2018). The study area offers a typical mosaic of paved surfaces, residential, commercial, and industrial areas, parks, and gardens containing diverse ornamental and native plants. The region's native semi-arid ecosystems transition from short grass steppe east of Denver, where tree growth is generally limited, to riparian zones and mesic spaces to a gamble oak-juniper woodland up the foothills of the Rocky Mountains to the west of Denver (Sims et al. 1978;Decker 2007;Hinners and Hjelmroos-Koski 2009). Average annual rainfall measures between 35 and 40 cm, with most precipitation occurring between May and July. Mean daily temperatures reach lows of 0°C in January to highs of 23°C in July (National Oceanic and Atmospheric Administration 2018). Approximately half of the study sites consisted of exotic landscaping, including turf grass, shade trees, and ornamental flower plantings, often requiring supplementary water input during the summer months (Hinners and Hjelmroos-Koski 2009).

Defining the urban-rural gradient
For a gradient analysis it is important to define metrics that identify "urban," "suburban," and "rural," spaces. Comparable studies have opted to use land use/land cover classification systems to identify urban, agricultural, natural, or other land use variables within the study system (Leong et al. 2014;Verboven et al. 2014). Others have chosen to use impervious surface data as a proxy for increasing urban intensity (Fortel et al. 2014;Geslin et al. 2013Geslin et al. , 2016Lagucki et al. 2017;Choate et al. 2018). We visualized urban landcovers surrounding our sites using ArcGIS (version 10.5; Environmental Systems Research Institute, Redlands, CA, USA). Land use data were downloaded from the National Land Cover Database (Homer et al. 2015). Adhering to the methods by Jha and Kremen (2013) and Quistberg et al. (2016), urban land cover types were divided into seven categories according to levels of resistance to bee nesting: 1) low intensity urban (<50% impervious surface; NLCD values 21 and 22), 2) medium intensity urban (50-79% impervious surface area; NLCD 23), 3) high intensity urban (80% impervious surface area; NLCD 24), 4) forest (evergreen and deciduous; NLCD 41 and 42); 5) undeveloped open space (NLCD 31, 52 and 71), 6) agriculture (NLCD 81 and 82), and 7) wetland and open water (NLCD 11,90,and 95).. Given that we were also interested in rural areas which may have relatively low degrees of urbanization, we also chose to classify the relative degree of urbanization within our 12 study sites by using percent impervious data directly using surface raster data with 30 m resolution downloaded from the National Land Cover Database (Xian et al. 2011).
We extracted urban landcovers and impervious surface from 2 km radii surrounding sites as many bee species forage locally and seldom travel beyond a 2 km radius from their nest in urban settings (Garbuzov et al. 2015); sites were also located approximately 4 km or more from the other sites. Sites containing an average of 50% or more of imperviousness were classified as "urban," sites with an average of 25% to <50% impermeability as "suburban," and sites with an average of <25% impermeability as "rural" (Fig. 1). Four sites within each category were sampled. Urban sites include U01 -U04. U01 is a university campus located near downtown Denver. U02 and U03 are urban community gardens, and U04 incorporates a small plot behind a commercial building. Suburban sites include S05 -S08. S05 is located near the intersection of two major roads and is near the center of a commercial district. S06 and S08 are both residential, and the S07 site is in the lawn of a large suburban church. The rural sites are R09 -R12. R09 and R10 are both located on the outskirts of suburban Denver; the first is a community garden and the second is in a residential subdivision. R11 is a state park north of the city, and R12 is a privately-owned horse ranch outside of town.

Local habitat characterization
For each of the 12 sites, local and landscape habitat characteristics were determined to identify the human-environmental factors that drive changes in pollinator assemblages along the gradient. Local habitat characteristics measured or identified were based on Bates et al. (2011). Richness of flowering forbs, woody plants, and trees, both in and out of bloom, were sampled once in May and once in July, using 25-m radius circles centered around each study site. Plants were identified in the field whenever possible using flowers or, when the plant was out of bloom, using leaf and stem morphology. When plants could not be identified in the field, floral samples or photographs were taken for identification in the laboratory. Floral abundance was estimated using a logarithmic scale (McFrederick and LeBuhn 2006) where floral "units" are composed of a single flower, umbel, head, spike, or capitulum (Carvell et al. 2007). Though this is a coarse measure of floral resources, we maintain that two sampling periods allowed us to identify most floral species present at each site, regardless of phenological stage. Percent canopy cover was also estimated within 25 m of study site centers (Matteson et al. 2013), and we assessed wind exposure by estimating the percent of land within a 50-m radius not protected by tall shrubs, fences, hedgerows, trees, or buildings (Bates et al. 2011). Floral abundance (given its coarse estimation), percent canopy cover, and wind exposure are all included in our published data (Online Resource 1) but were not used in statistical analyses in order to reduce the number of potential predictors in our models. Instead we focus on floral richness as predictors of bee species richness in subsequent analyses.

Bee sampling
Bees were captured using pan traps and sweep netting based on The Bee Inventory Plot methods (2018) in the summer of 2017. Westphal et al. (2008) found that pan trapping is the most efficient method of collecting bees across a broad geographical range. However, pan trapping can lead to biased estimates of bee community composition as it disproportionately attracts small generalist bees such as those in the genus Halictus (Cane et al. 2000) and does a poor job of collecting larger social bees such as honeybees and bumble bees, as well as specialist species such as Colletes (Roulston et al. 2007;Wilson et al. 2008). Despite its shortcomings, pan trapping remains one of the most common sampling methods for bee studies. Sweep netting and/or hand searching are often implemented in conjunction with pan trapping to overcome biases.
Bee surveys were conducted approximately once every 4 weeks from May to August on sunny days with light wind and no precipitation. For pan traps, we suspended two feet from the ground three different colored (white, yellow, and blue) plastic soufflé cups (Solo brand) at each of the sites for 24 h once a month. The cups were half-filled with a solution of 90% water and 10% All Free Clear liquid laundry detergent (to break the water surface tension). Bees not attracted to traps (e.g., pollen specialists) were captured within 25 m of each site using a long-handled sweep net for 30 min once a month. Sweep netting occurred 1 to 5 days before or after pan traps were placed to avoid interference between sampling methods. Because bee activity may vary throughout the day (Rader et al. 2013;Stanghellini et al. 2002), sites were visited in varying order each month such that sites were sampled both in the morning and in the afternoon. All honeybees, except the first captured specimen, were recorded and released after the 30 min of sweep netting. All recorded specimens were included in analyses.
Specimens were strained, rinsed in water, and temporarily stored in denatured isopropyl alcohol and pinned for identification. Because Agapostemon angelicus and A. texanus females are morphologically identical (Roberts 1972), these Fig. 1 Study sites in the greater Denver, CO metropolitan area. Urban sites include U01, U02, U03, and U04. Suburban sites include S05, S06, S07, and S08. Rural sites include R09, R10, R11, and R12 species were analyzed together as a single species aggregate (McIntyre and Hostetler 2001). All specimens were labeled and sorted to species or morphospecies and are housed in the University of Colorado Museum of Natural History Entomology Section. Ecological characteristics for each species, such as mean body length, sociality (eusocial, solitary, or parasitic), preferred nesting substrate (soil, cavity, or hive), pollen specificity (oligolectic or polylectic), and phenology (i.e., time of first activity) were determined from primary literature, catalogs, and other references (Wilson and Carril 2016;Michener 2000;Quistberg et al. 2016;Ascher and Pickering 2012;Gibbs 2011;Packer 1993;McGinley 1986;McGinley 2003;LaBerge and Ribble 1975;Bouseman and LaBerge 1978;Scott et al. 2011;Brooks 1988;Wright et al. 2017;Cockerell and Sumner 1931;Cockerell 1903). Species origin (whether a species is native or exotic to North America) was also identified per Scott et al. (2011). Oligolectic species were identified as those that forage on a single plant family or a subset of flowering genera; no attempt was made to evaluate the extent of oligolecty. Body size was quantified according to the full body length of each specimen, and collected bees were classified as small (<7 mm), medium (7-12 mm), and large (>12 mm) (Quistberg et al. 2016). Characteristics were left as undetermined when species-specific ecological data were incomplete. There were 14 species with undetermined sociality, a n d 11 s p e c i e s w i t h u n d e t e r m i n e d p h e n o l o g y (Online Resource 2); these were excluded from further analyses.

Bee community composition
All analyses were conducted in R version 3.5.1 (R Core Team 2013). Bee abundance was calculated as the total number of recorded bees, and species richness as the total number of bee species documented across the season at each study site. Species diversity was calculated using Shannon's H′ because species richness may increase with greater diversity of land uses within each study site (Matteson and Langellotto 2010;Verboven et al. 2014). Bee abundance data were logtransformed to generate a normal distribution. We used oneway Analysis of Variance (ANOVA) to determine whether bee abundance, richness, and diversity differed between community types.
To determine our sample effectiveness of the bee community, we combined all bee captures across sites and explored regional species richness by fitting an accumulation curve and estimating asymptotic species richness. Using the vegan package (Oksanen et al. 2018) for (R Core Team 2020) we generated a species accumulation curve using the random method with 100 permutations to estimate the total number of unseen species in the study area. Total expected species richness was assessed using the Chao estimator because bee abundance data contains many singletons and it is a relatively unbiased estimator for species-rich assemblages (Fortel et al. 2014;Walther and Morand 1998). We then ran an asymptotic model on the species accumulation curve using the "fitspecaccum" function in vegan to predict the number of sites needed to be within 80% (as a reasonable goal for estimated community coverage) of the estimated Chao. The "fitspecaccum" function is adequate for non-linear species accumulation models such as in asymptotic models.
Differences in local floral richness among sites was tested with a one-way ANOVA. Landscape variables were found to be highly correlated with one another. Therefore, we used a principal component analysis (PCA) to incorporate the variables into a single variable that explains most of the variation in the original variables. The PCA yielded both positive and negative values and explained 56% of the variance. We removed negative values, which then explained 61% of the variance; the remaining values were combined to create a single component called "urbanization". A table describing the landscape characteristics and the component urbanization value for each site can be found in Online Resource 3.
To visualize dissimilarities in bee community composition between urban, suburban, and rural sites, non-metric multidimensional scaling (NMDS) was used (vegan package) based on Bray-Curtis dissimilarities, which accounts for the incidence of unique species between two sites relative to the total of all species sampled (Carper et al. 2014). This measure has been shown to have a strong relationship with natural ecological distance when tested with simulated data (Faith et al. 1987). To evaluate the relative fit of the NMDS, final stress measurements of the ordination were recorded, with values <0.05 indicating robust support for the ordination (Carper et al. 2014). In addition, community ellipses representing the 95% confidence interval around each community centroid were plotted, and permutational analysis of variance (Adonis function, with 999 permutations) was used to test the significance of whether the centroids differed from one another in the NMDS To visualize how bee communities respond to different environmental factors, we overlaid the "urbanization" and flower diversity variables as vectors on the NMDS, and tested the significance of these factors using "envfit." This function performs F-tests on vector directions to determine if a variable is significantly driving differences between communities in the ordination space.

Bee community response to urbanization
Linear regression was used to explore the relationships between total bee abundance, species richness, and diversity and urbanization. Percent imperviousness was used as the predictor variable. Correlation coefficients were used to test the strength of relationships, and one-way ANOVA was used to test for model significance. To explore guild-specific responses within the bee community to urbanization, separate stepwise linear regressions (step function in R) were used with the abundance of bees with different life history characteristics (i.e. mean body length, sociality, preferred nesting substrate, pollen specificity, and phenology) as responses. Stepwise linear regression is a method of fitting regression models by adding or subtracting predictive variables to find the most parsimonious model that reports the most predictive factors driving bee community dynamics (Efroymson 1960;Hocking 1976;Draper and Smith 1981;SAS Institute 1989). Thus, in each step a variable is considered for addition into or subtraction from the set of explanatory variables based on Akaike Information Criterion (AIC). We included the amount of "urbanization" in the landscape and local floral richness as fixed effects and conducted forward and backward steps in model selection. Categories containing bees with indeterminate characteristics were excluded from the analysis. If a guild type was recorded at half the sites or less, it was also removed from the analysis. All abundances were log-transformed to meet the assumptions of normality.

Results
In all, 1611 bees belonging to five families (Andrenidae, Apidae, Colletidae, Halictidae, and Megachilidae), 23 genera, and 81 species, were recorded throughout the field season. Most bees collected were native. Exotic species accounted for 9% of collected bees and included Apis mellifera, Megachile rotundata, and Hyleaus leptocephalus. The most abundant species, accounting for 52% of all specimens collected, was Halictus tripartitus (Online Resource 4): a small, eusocial flower generalist that emerges in early spring. Other abundant species included A. mellifera, Lasioglossum semicaeurelum, Agapostemon texanus/angelicus, and L. sisymbrii. All species are polylectic, collecting pollen from a variety of flowers. The most represented families included Halictidae (85% of individuals, 41% of species) and Apidae (7% of individuals, 20% of species). Thirty-seven species were represented by a single individual.
While rural sites had over four times more bees with nearly double the species richness, and more than twice the diversity as suburban and urban sites, we found no significant differences between site categories for either measure (F 2,9 = 3.47, p = 0.077; F 2,9 = 2.02, p = 0.190; and F 2,9 = 3.61, p = 0.071, respectively, Fig. 2, Online Resource 5). This was likely driven by relatively small sample sizes. According to the Chao estimate (Fig. 3) the study area harbored approximately 129 (± 20) species; thus, we recorded 63% of the total species predicted for the region. To obtain 80% (103 species) of the total estimated species richness, we would have needed to sample approximately 40 sites (Fig. 3). Still, given that Scott et al. (2011) documented 156 species for Denver County, these results suggest that our sampling estimates are consistent with documented local species richness.
NMDS indicated little overlap between urban and suburban communities, and no overlap with rural bee communities (Fig. 4). The final stress of the two-dimensional NMDS was 0.176, indicating acceptable confidence in the ordination distances of each community with significantly different centroids (PERMANOVA: F 2,11 = 1.37, p = 0.04). In addition, results from smooth surface fitting suggest that separation in communities was likely driven more by urbanization (envfit: R 2 = 0.70, p = 0.003) and less so by floral richness (envfit: R 2 = 0.27, p = 0.236) across the gradient from urban to rural sites (Fig. 4). Rural sites had the greatest floral richness with 113 recorded species. Suburban sites had the lowest floral richness with 44 species, and urban sites had intermediate floral richness with 64 species, though again differences were not significant (F 1,2 = 1.46, p = 0.283). Conversely, rural sites had the lowest percent imperviousness with a mean of 7.7%, suburban sites had an intermediate amount of imperviousness with a mean of 35.7%, and urban sites had the highest amount of imperviousness with a mean of 62.0% (F 1,2 = 38.45, p ≤ 0.001; Online Resource 6).
Most bees across the study area were small (75% of individuals, 52% of species), social (66% of individuals, 14% of species), ground nesting (94% of individuals, 69% of species), and pollen generalists (94%, of individuals 73% of species). Parasitic bees composed a very small fraction of the bee community for the study area (1% of individuals, 5% of species), and pollen specialists were also uncommon (5% of individuals, 22% of species). Stepwise linear regression suggested guild-specific responses within the bee community to urbanization and floral resources, though not all guilds responded to either metrics (see Table 1). In all cases where models indicated significant responses to the explanatory variables, responses were always positive for floral richness and always negative for urbanization and no bee guild was positively associated with urbanization. In general, small bees declined with increasing urbanization (Fig. 6a) while large bees increased with increasing local floral diversity. Ground-nesting bees were not impacted by urbanization, while hive-nesting bees responded significantly to floral richness but not urbanization. Polylectic bees declined significantly with increasing urbanization (Fig. 6b), though oligolectic bees did not respond to either explanatory variable. Neither social nor solitary bees responded to urbanization. Only early-season bees responded to urbanization (Fig. 6c). There was also no response from either native or exotic species to urbanization, although interestingly exotic species increased with floral richness.

Discussion
These results demonstrate that while the Denver metro area has a diverse bee community, both bee abundance and diversity decrease with increasing urbanization along the rural to urban gradient. Moreover, multiple guilds of bees, including small, polylectic, and early season bees, responded negatively to increasing urbanization and we found no positive relationship for any guild with increasing impervious surface cover. Conversely, floral richness had positive effects on bee communities, and both metrics likely drive variation in communities across the Denver rural-urban gradient. While these results suggest that urbanization in Denver does have negative effects on bee communities, it also suggests that increasing floral richness could help mitigate the impacts of urbanization on bees in general, and that the benefits of increased floral richness may not necessarily be guild specific.
Contrary to our hypothesized guild specific-responses within the bee community to urbanization, our results showed that all bee guilds either responded negatively to urbanization or not at all. For instance, the most dominant bee guild collected in this study was small bees which were highly abundant in suburban and rural locations, but much less abundant in urban sites. Small bees have short foraging distances (Gathmann and Tscharntke 2002;Greenleaf et al. 2007) and may be more sensitive to highly dispersed floral resources and physical barriers, such as busy roads and tall buildings that can prevent movement between resource patches (Ahrne et al. 2009;Banaszak-Cibicka and Zmihorski 2012;Geslin et al. 2016). Instead, medium and large bees are likely less impacted by habitat fragmentation (Fortel et al. 2014;Steffan-Dewenter and Tscharntke 1999;Westphal et al. 2006) and may have higher performance in urban areas with high impervious surfaces because of their increased flight mobility. Similarly, while generalist and cavity-nesting bees have been shown to be greater in urban areas compared to specialists and groundnesting bees (Banaszak-Cibicka and Zmihorski 2012; Hernandez et al. 2009;Deguines et al. 2016;Frankie et al. 2005;Geslin et al. 2013;Jedrzejewska-Szmek and Zych 2013;Pardee and Philpott 2014;Tonietto et al. 2011;Fitch et al. 2019), our results suggest that they too are susceptible to urbanization. We found very few cavity-nesting bees, potentially as a result of the study region; differences between this study and other findings could indicate that bee community responses to urbanization are region specific. The current study occurred in a high-plains semi-arid climate, which may be less supportive of cavity-nesting bee populations due to the short length of the growing season, low availability of Fig. 3 Species accumulation curve combining captures from all study sites. Suggested regional species richness was approximately 37% higher than sample species richness. The dotted line represents the asymptotic model on the species accumulation curve, and the vertical dashed line identifies the approximate number of sites needed to sample within 80% of the Chao estimate Fig. 2 Log-transformed bee abundance a, richness b, and diversity c across urban, suburban, and rural communities. Bars represent standard error of the mean wood nesting resources, or elevational constraints. Thus, disparities between results may also be due to regional climate differences between the study areas. Given these regionspecific responses of bee guilds to urbanization, more studies, such as this one, are needed to explore potential mechanisms driving community responses to urbanization.
Although overall bee abundance and diversity decreased with increasing urbanization, we also found that floral resources had a positive effect on bee communities, particularly in urban spaces. For example, the U03 site, an urban area with high floral richness, supported substantially more bee species compared to other urban and suburban sites, suggesting that local floral richness, more than the surrounding urban landscape, is important in sustaining a diverse bee community. These results are consistent with other research. Ahrne et al. (2009), demonstrated that bumble bee abundance and species richness was positively associated with local floral abundance and negatively associated with increasing urbanization in Stockholm, Sweden. They concluded that landscape changes associated with urbanization negatively affect bumble bee diversity, while local factors, such as floral abundance, positively affect bumble bee abundance. Similar patterns were found for wild bees and other pollinators in Boulder, CO (Kearns and Oliveras 2009), Paris, France (Geslin et al. 2013;Geslin et al. 2016), New York City, NY (Matteson et al. 2008), Grand Lyon, France (Fortel et al. 2014), and Toledo, OH (Lagucki et al. 2017). While we sampled floral diversity only twice (in May and July) with a relatively coarse metric, that we detected  Given the lack of cavity-nesting and oligolectic bees, we were unable to include them in the analyses significant relationships suggests that floral diversity is a strong predictor of bee diversity. Among our sites, the U03 site stood out as a community garden and others agree that community gardens have great potential to support urban bees and other pollinators (e.g., Ahrne et al. 2009;Gotlieb et al. 2011;Samnegård et al. 2011;Hinners et al. 2012). For example, Kaluza et al. (2016) found that nectar forager abundances were higher in community gardens compared to plantations and adjacent natural forests in Queensland, Australia, likely due to a constant nectar availability throughout the season. On the other hand, while the U03 site contained a rich bee community, the U02 site, another small urban community garden, maintained substantially fewer individuals and species, possibly because the garden is embedded in an impervious landscape. Given the greater percentage of medium bees captured at the U02 site, it was likely too physically remote for many smaller bee species to access. Similarly, Matteson et al. (2008) found substantially reduced native bee richness in New York City community gardens compared to the surrounding natural landscape. Therefore, urban gardens should be accessible and produce floral resources throughout the season to be most beneficial to bees.
Notably, the U04 site maintained higher abundances and greater diversity of bees relative to the U02 and U01 urban sites. The small lawn space at U04 contained many ruderal species compared to the surrounding turf lawns. Thus, U04 was likely highly attractive to the local bee community in an area with sparse resources. This has significant implications for the Denver metropolitan area where ruderal sites are disappearing with the incursion of housing construction and infrastructural development. Similarly, Jedrzejewska-Szmek and Zych (2013) found that even small patches of natural vegetation in urban spaces may support appreciable pollinator diversity. Additionally, Larson et al. (2014) found that weeds grown in turf lawns were an important food source, and may even act as resource corridors or steppingstones for certain bee species in Lexington, KY. The importance of maintaining diverse local resource patches for bees and other pollinators is well documented (e.g., Ahrne et al. 2009;Kearns and Oliveras 2009;McFrederick and LeBuhn 2006;Potts et al. 2003;Threlfall et al. 2015). Therefore, patches of remnant vegetation, if they are accessible, may be critical to sustain urban bee populations.
There was a greater abundance of exotic bees, namely honeybees, in suburban spaces relative to urban and rural areas, possibly due to an increase in the amount of exotic and ornamental vegetation planted in these spaces (Thompson et al. 2003;Frankie et al. 2005;Matteson et al. 2008) or an influx of managed honeybee hives. There is increasing interest in maintaining honeybee (Apis mellifera) hives in Denver and other cities to help supplement pollination in community gardens (Cane 2005;Peters 2012). However, wild bees may be just as or more relevant for suburban pollination. Lowenstein et al. (2015) found that most bees visiting experimental crop plants in Chicago were wild. To our knowledge, there is no vegetation in Denver that relies exclusively on honeybee pollination. It may prove beneficial to increase efforts to improve wild bee populations in the city to enhance urban pollination and concurrently reduce reliance on a single pollinator species.
Ultimately, the implications of this study suggest some simple recommendations for mitigating the negative impacts of urbanization. Because bees responded positively to floral richness, planting a wide variety of bee plants should be an effective strategy to better manage bee communities along the urban-rural gradient. This was especially evident in the U03 site, a florally rich urban site. Additionally, it has been shown that wild bees prefer native vegetation over exotic plant species (Frankie et al. 2005) and that they require various substrates for nesting (Cane 2005). Therefore, other simple adjustments could include planting more native pollinatorfriendly plants in Denver's parks and open spaces as well as leaving patches of fallen wood, leaves, and bare ground for nesting substrate. Prioritizing bare ground may provide the greatest benefit given the preponderance of ground-nesting bees in the Denver metro area. These simple changes can furthermore benefit other wildlife populations, promote biodiversity, and enhance ecosystem services in the city (Wratten et al. 2012;Hipólito et al. 2016). Furthermore, because bees responded negatively to urbanization, it should prove beneficial to incorporate more natural park space in citywide management. For example, improvements may include planting native shortgrass-steppe grassland species in parking lot islands and along road verges instead of mulching and exotic shrub plantings. Finally, actively involving Denver residents in pollinator habitat restoration efforts would advance local knowledge on the importance of maintaining urban pollinators and encourage residents to manage lawns more sustainably. However, such initiatives will demand cooperation among policy makers, city planners, conservation scientists, and the community in general to effectively move forward with urban pollinator conservation (Hall et al. 2016). For example, revising development plans to compel endorsement from municipal governments and residents when considering changes to the built environment can help advance pollinator conservation in the city. Given that 67% of the world's population will live in cities by 2050 (United Nations 2014), addressing pollinator declines and their impacts on food resilience and human well-being should be a common goal to all.
identifications and Peter Anthamatten at The University of Colorado Denver for his insight, support, and compositional expertise. We also thank the Auraria Higher Education Center, Duane McClanahan, Denver Urban Gardens, Christy Shires, Ellen Haroutunian, the Jay family, St. Andrews United Methodist Church, Barr Lake State Park, and Gwen Ka'awaloa for allowing us to survey bee populations on their properties. Finally, we would like to thank Micah Birdshire, Julie Borja, Maida Pearce, Becca Jay and family, Audrey Yanos and family, Ellen Haroutunian, Paige Alexander, and Sarah Lambert for helping with field collections.
Funding information This study was funded by The Garden Club of America Centennial Pollinator Fellowship, the American Association of Geographers Biogeography Specialty Group, the University of Colorado Denver Graduate School and the College of Liberal Arts and Sciences.

Compliance with ethical standards
Grant support Kristen Birdshire has received research grants from The Garden Club of America Centennial Pollinator Fellowship, the American Association of Geographers Biogeography Specialty Group, the University of Colorado Denver Graduate School and the College of Liberal Arts and Sciences.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.