Introduction

Urbanization and conversion of land for agriculture is altering landscapes globally, with levels of urban land cover by 2030 projected to be triple what they were at the turn of this century (Seto et al. 2012). The corresponding loss, fragmentation and degradation of natural habitats is expected to have devastating effects on biodiversity (McDonald et al. 2013). The proliferation of impervious surfaces including roads, the establishment and spread of non-native plant species, climate change and pollution are likely to worsen these effects. Urbanization is unlikely to affect all insect taxa equally, however, and developed areas have been shown to be of great value to some taxa (Hall et al. 2016; Turo and Gardiner 2019; Wenzel et al. 2020). Bees appear to be relatively tolerant of urban environments where they can be equally or even more diverse than in non-urban land use types (Baldock et al. 2015; Kaluza et al. 2016; Lerman and Milam 2016; Banaszak-Cibicka et al. 2018). Such findings highlight the potential role urban areas can play in conserving this group of insects.

Various local and landscape factors are likely to affect urban bee diversity (Matteson and Langellotto 2010; Williams and Winfree 2013; Lerman and Milam 2016; Ayers and Rehan 2021; Turo and Gardiner 2019; Wenzel et al. 2020). Several studies have reported a positive relationship between bee diversity and the amount of natural habitat in the surrounding landscape (Kremen et al. 2002; Steffan-Dewenter et al. 2002; Carré et al. 2009; Le Féon et al. 2010). In regions that were extensively forested prior to human colonization, a substantial proportion of the bee fauna consists of forest-dependent species. For example, Smith et al. (2021) estimated that about one third of bee species in the northeastern United States are forest-associated while a similar proportion is associated with anthropogenic habitat. The remaining species were categorized as habitat generalists. The researchers found forest area to have a positive effect on the abundance and richness of forest-associated bees but no effect on habitat generalists. Harrison et al. (2018) concluded that forest bees cannot persist in anthropogenic habitats, and are being replaced in agricultural and urban environments by species adapted to open habitats. The composition of bee communities in urban environments is thus likely to vary greatly from urban to rural environments.

A key challenge in research on the effects of landscape context on biodiversity is to determine the scale(s) at which species are influenced by various landscape factors, i.e., the ‘scale of effect’ (Fahrig 2013; Holland et al. 2005). Because the scale of effect varies widely among taxa and response variables, and because significant landscape effects can be missed if the wrong scale is considered (Holland et al. 2005), it is important to test a wide range of scales to fully explore the effects of landscape context (Jackson and Fahrig 2015). Many bee taxa, especially small species, limit their foraging to within a few hundred meters of the nest (Gathmann and Tscharntke 2002; Zurbuchen et al. 2010), suggesting these species may respond to smaller scales. Indeed, several studies have shown that bees respond more strongly to smaller scales (Saturni et al. 2016) and that the scale of effect is proportional to body size, a correlate of foraging distance (Benjamin et al. 2014; Tscheulin et al. 2011).

Here we present the results from a study aimed at understanding how landscape context affects bee diversity in residential habitats in an area of the southeastern United States experiencing rapid urbanization. We were particularly interested in understanding how these insects are affected by the amount of development and forest cover at different spatial scales. The urban-forest interface was of particular interest. We predicted that bees would be negatively affected by development and positively influenced by the amount of forest on the landscape. We further predicted that these effects would be strongest at smaller scales.

Methods

Study sites

This study took place in northeast Georgia, USA, a region that was largely deforested for cotton production beginning in the mid 1800s. After the collapse of the cotton industry in the early twentieth century, much of the region became reforested through natural succession and timber production. This trend reversed in recent decades, as more forests were lost to development, particularly the creation of residential communities (Miller 2012). Most sites were in and around Athens, Clarke County. Athens is the county seat of Clarke County. With a 2020 population of 130,081, it is the 6th largest city in Georgia and the 217th largest city in the United States. Athens is currently growing at a rate of 0.82% annually and its population has increased by 12.67% since the most recent census, which recorded a population of 115,452 in 2010. Spanning over 118 square miles, Athens has a population density of 1118 people per square mile. (World population review https://worldpopulationreview.com/us-cities/athens-ga-population).

We obtained permission from homeowners to sample bees on 34 and 50 residential properties in 2019 and 2020, respectively (Fig. 1; Online Appendix C for Lat., Long.). Sites included several land cover types (Fig. 1) including development, agriculture and forest. Forests accounted for an average of 38% (range 0–90%) and 44% (range 6–75%) of land area surrounding our sampling sites at the spatial extents 0.2025 km2 and 14.98 km2, respectively. These spatial extents correspond to scales 2 and 12 used in the analysis (see “Results”, Online Appendix A). All properties were at least 1600 m apart to avoid spatial autocorrelation, and by design, represented a continuum from older properties within the city of Athens to peripheral suburban neighborhoods established on former agricultural land or forests.

Fig. 1
figure 1

Map of the study area showing sampling locations for 2019 and 2020

Full size image

Bee sampling

Bees were trapped using a set of three colored plastic pan traps (white, yellow, and blue) filled with soapy water (Dawn dishwashing soap). Although pan traps are known to capture smaller bees more effectively than larger bees (Cane et al. 2000; Roulston et al. 2007), they are a highly standardized and efficient method allowing simultaneous and consistent sampling of a large number of sites. The traps were placed in areas with direct sunlight, at least 10 m from the nearest mature tree, and were arranged in a straight line with 1 m separation. Wire stands were used to hold the traps in place about 30 cm above the ground. Sampling took place 2 days per week from 19 June to 6 August in 2019 and from 30 March to 26 June in 2020. The contents from the three bowls at each sampling site were combined into a single jar and returned to the laboratory. Insects were strained from the pooled water samples and stored in ethanol until bees could be sorted, pinned and identified. They were identified by MDU using an established reference collection and a variety of printed and online resources (Gibbs 2011; Gibbs et al. 2013; Mitchell 1960, 1962; https://www.discoverlife.org). Voucher specimens are retained at the University of Georgia Natural History Museum.

Bee data and analyses

All statistical analyses were conducted in R version 4.0.3 (R Core Team 2020) unless otherwise stated. For each year, we calculated bee diversity at each site using Hill numbers, a mathematically unified family of diversity indices (Chao et al. 2014), using the hill_taxa function of the hillR package (Li 2018). Hill numbers allow for comparisons between species richness estimates allotting more weight to rare (q = 0), common (q = 1), and dominant (q = 2) species (Chao et al. 2014). Hereafter, Hill number estimates will be referred to as species richness (q = 0), Shannon diversity (q = 1), and Simpson diversity (q = 2) as recommended in a recent review (Roswell et al. 2021).

True species richness was estimated at each site with the Chao1 estimator (Chao 1984; Colwell and Coddington 1994) using the chao1 function of the rareNMtests package (Cayuela and Gotelli 2014). Our bee data (total abundance, richness and diversity indices) were tested for spatial autocorrelation by calculating Moran’s I (Dormann et al. 2007) using the ape package’s Moran.I function (Paradis and Schliep 2019).

Landscape analysis

We quantified landscape composition at 11 spatial extents (scales), ranging from ~ 0.20 to 2.2 km in radius (i.e., 0.20 to 14.98 km2) (Online Appendix A). These scales were chosen as they encompass the extent of documented foraging ranges of bee species (Taki et al. 2007; Winfree et al. 2007; Watson et al. 2011). At each site the percent of the landscape occupied by each cover type for each spatial scale was calculated using the most recent USGS National Land Cover Database (NLCD) 2016 (Dewitz 2019). Land cover categories considered included total forest cover (i.e., “AllF”: NLCD classes 41, 42, 43, 90), agriculture (i.e., “AllAg”: NLCD classes 81, 82), and development (i.e., “AllDev”: NLCD classes 21, 22, 23, 24).

Further descriptions of landscape data are at: National Land Cover Database. (accessed 9/1/2021 date).

We also calculated land cover type diversity, which we calculated in FRAGSTATS v4 (McGarigal et al. 2002) for each site and scale. Landscape diversity (LD) represents Shannon’s diversity based on the number and evenness of land cover classes from the NLCD data (forest, development and agriculture). Thus, sites with higher LD values will be surrounded by a greater variety of land uses. Importantly, LD does not consider which land use types are present, or their relative abundances. So, for example, a site could have high LD but all of that diversity could come from highly modified, anthropogenic land uses.

Scale of effect

We determined the ‘scale of effect’ (Jackson and Fahrig 2012) for each landscape variable and bee response variable combination by identifying the spatial scales (Online Appendix A) most highly correlated with each response variable (Holland and Yang 2016). We used Pearson correlation coefficients for response variables with normal distributions and Spearman correlation coefficients for those with non-normal distributions. For each land cover type, the scale(s) with the highest and lowest correlation coefficients (to include negative correlations) were selected for initial analyses.

Linear modeling

Generalized linear models on estimates of total bee diversity (Hill numbers 0, 1 and 2), were constructed separately for 2019 and 2020 using the MASS package’s (Venables and Ripley 2002) stepAIC function. For each Hill number, the scale of effect analysis described above determined which scale(s) to include for each landscape parameter, with the scale which had the maximum and minimum correlation coefficient values included in each Hill number’s initial model. All landscape parameters were included as fixed effects. A Gaussian distribution was used for q = 0, q = 1 and q = 2. Model parameter variance inflation factors (VIF) were confirmed to be under 4, removing the parameter with the highest VIF until all met this criterion. Model assumptions were checked using the DHARMa package’s (Hartig 2019) simulateResiduals function.

Generalized linear models were also constructed to model the responses of individual species abundances to land covers for species that occurred in at least 10% of the sites and with a total count of more than 25 individuals (Brin et al. 2016). Sampling year was included as a fixed effect for each model to assess if sampling year had an effect on each species. A full initial model containing year and the full set of landscape scale parameters (Table 2 column 2) for each species was created with MASS’s glm.nb function, to fit negative binomial generalized linear models with a log link appropriate for abundance count data. A stepwise procedure to select the best model for each species was performed via the stepAIC function.

Community analysis

Sites were classified by primary cover types: agriculture, development, and forest. As we were interested in the effect of forest cover, the primary cover was defined as the cover type with the plurality (largest percentage of land use cover) at the spatial scale with the strongest correlation with forest. Pluralities often were less than 50% of cover, with multiple cover types having similar quantities.

We assessed how bee communities differed among land cover types by performing non-metric multidimensional scaling (NMDS) ordination of communities using Bray–Curtis dissimilarity. The species matrices used for community analyses were Hellinger transformed to reduce over-emphasis of rare and abundant species. Permutational anovas (PERMANOVAs) compared Bray–Curtis dissimilarities of sites to determine if bee communities differed statistically by primary land cover type (agriculture, forest, or development. Community ordinations by sample location were made with vegan’s (Oksanen et al. 2020) metamds function using Bray–Curtis dissimilarities with axes set to k = 3. Significant land covers (p ≤ 0.05) were vectorized onto the resulting ordinations via non-parametric linear regression using the vegan package’s envfit function. To gain further insight into the contributions of each bee species to community dissimilarities among cover types we performed indicator analysis (multipatt function from the indicspecies package (Cáceres and Legendre 2009) and SIMPER analysis using vegan’s SIMPER function on the Hellinger transformed species matrix.

Results

Bee community diversity

A total of 4938 bees representing 111 species (Online Appendix B), from 28 genera in five families were captured in pan traps during the 2 years. The species represented 20% of those reported by Schlueter (2020) (https://native-bees-of-georgia.ggc.edu) for the state of Georgia. Spatial autocorrelation between sample locations was assessed using Moran’s I for bee abundance and diversity and bees were found to be not significantly autocorrelated (in all cases, P > 0.05). Hill numbers for sample sites within each year (Online Appendix C) revealed generally higher species richness and diversity during 2020 when the sampling interval included spring months and the number of sample sites was doubled. The number of species, Chao1 estimator and % coverage were 49, 59 and 83% (for 2019); 99, 125 and 79% (for 2020); 111, 138 and 80% (for combined years).

Bee community diversity was significantly affected by land cover type (Table 1) and scale often influenced direction of effect. Development correlated positively with bee diversity (q = 1) during both years at a small scale (2 and 4). Agriculture correlated positively with q = 1 and q = 2 diversities during 2019 at a large scale (12), but negatively with species richness at a small scale (5). Bee diversity was negatively affected by forest cover at small scales in both years, but we detected a positive effect of forest at large scales (12) for all diversity indices in 2020. Landscape diversity correlated negatively with all diversity indices in 2019, but was not a significant factor in the 2020 models.

Table 1 Linear model coefficient, T value and significance level for each predictor land cover variable (scale, L = Large, S = Small) identified in the final model for community diversity indices (Hill numbers)
Full size table

Individual species models

Twenty-two species met the selection criteria for individual species abundance modeling (more than 25 individuals collected in 10 or more sampling events) and all model assumptions (Table 2). The negative binomial distribution was used for all species models. At the small landscape scale, forest coverage had a consistently significant, but negative, effect on Lasioglossum species (four out of nine, L. illinoense, L. tegulare/puteulanum, L. trigeminum, L. zephyrum,). Overall, forest cover was identified as a model parameter for 13 species. Forest coverage’s effect switched from mostly negative to mostly positive as the landscape scale increased. Development had a primarily positive effect, especially at small scales. Landscape diversity had a mostly negative effect on individual bee species abundances in this study. In contrast to other studies, landscape diversity was negatively associated with bee abundance possibly due to the difference in setting (semi-natural vs. modified, urban-suburban).

Table 2 Linear model coefficient, z value and significance level for each predictor variable [land cover (scale, L = Large, S = Small)] identified in the final model for species that met the count and incident criteria (bee species occurring in at least 10% of the sites and with a total count of more than 25 individuals)
Full size table

Community ordination analysis

NMDS ordinations revealed clustering of communities according to land cover types (Fig. 2). PERMANOVA tests confirmed significantly different communities between primary land cover types for 2019 (F = 1.92; P = 0.03, Fig. 2) and 2020 (F = 1.84; P = 0.003, Fig. 2). SIMPER and indicator species analyses identified species that contributed to the dissimilarities between groups and cover types and that typified each group (Table 4, Online Appendix E). For the 2019 data, indicator species analysis identified L. callidum, L. trigeminum, and L. truncatum as significantly associated with agriculture P < 0.05 and SIMPER analysis revealed that the three species contributed heavily to differences observed among the three main cover types. Similarly, indicator analysis identified B. bimaculatus and L. coreopsis as significantly associated with both Forest and Agriculture (P < 0.05).

Fig. 2
figure 2

Species ordination and envfit vectoring for 2019 and 2020 bee species, comparing communities by primary land cover type (agriculture, development and forest. Points represent the bee community at each sample location. Ellipses represent the standard deviation around the centroid of each cover community cluster. Vectors (FS_12_AllAg = Agriculture (Large Scale), FS_2 (or 5)_AllAg = Agriculture (Small Scale)) represent land covers significantly (P ≤ 0.05) influencing differentiation in communities with direction indicating the direction of positive association and length proportional to the strength of the effect. Communities differed significantly by cover types

Full size image

Discussion

The results from our linear modelling reveal the importance of both landscape context and spatial scale on bee diversity at the urban–rural interface and urban-forest interface. For example, we found bee diversity and the abundances of individual species to correlate negatively with forest cover at small spatial scales, i.e., 0.2025 km2 but positively at larger spatial scales, i.e., 14.98 km2. This likely points to the enriching effect of a more blended landscape (Fig. 2) relative to any one predominant landscape, including predominantly forest covered sites. Thus, as the landscape scale increases and forest cover can be included along with other beneficial cover types, the correlation with pollinator richness switches to a positive one. This suggests that more open areas may support higher bee diversity locally but that forest cover still promotes regional bee diversity, perhaps by supporting forest species that may not persist in open habitats (Harrison et al. 2018; Smith et al. 2021). It is noteworthy that the positive relationship between forest cover and bee richness at large scales was detected only in 2020. This may reflect the fact that sampling began much earlier in 2020 compared to 2019. Because forest bee diversity peaks during the spring months in our region (Ulyshen et al. 2022), it is not surprising that the benefits of forests to bee diversity were detected in 2020 but not 2019. Two species previously showing strong positive associations with forests, Ceratina strenua Smith and Osmia georgica Cresson, are cavity nesters, suggesting forests may provide important nesting resources for species with similar nesting behavior. The role of nesting resources may be an important driver of community structure (Fortuin and Gandhi 2021).

With regard to bee community composition there are significant differences in community composition for primarily forested areas relative to primarily agricultural or developed sites (PERMANOVA P < 0.05). This implies that including forest cover in the urban landscape matters for enriching regional pollinator communities but that forest-associated bees may not be particularly sensitive to the amount of forest cover. When looking at cover types (Fig. 2, SIMPER Analysis Online Appendix E) we can see distinct communities in blended Forest-Agriculture land cover and Forest-Development land cover, likely indicating a diversity bolstering effect of blended urban landscapes.

Previous efforts to sample bees from the forest interior support the conclusion that forests provide important habitats to these insects in the eastern United States. For example, Ulyshen et al. (2010) examined the vertical distribution of bees in a deciduous temperate forest in Ogelthorpe Co., GA and reported 71 species from five families. Thirty-one of those species were collected in our urban/suburban study, representing almost half of the forest bee species collected in that previous study. More recently, Traylor et al. (2022) caught 111 species of bees in forests within the boundaries of Clarke Co., Georgia, the same area investigated in the current study. Similarly, Urban-Mead et al. (2021) emphasize the importance of forest canopy for bee conservation reporting that 20% of New York’s 417 bee species were active in the springtime canopy, especially female wood-nesting and female social bees. We thus conclude that despite the negative correlations between bee diversity and forest cover at small scales, forests are important to regional bee diversity at the landscape level. The value of the urban forest canopy to bee conservation would be an interesting future avenue of research.

Compared to forests, developed areas sit at the opposite end of the disturbance spectrum. Although development is often associated with negative effects on bee diversity it had a positive effect on six of eight individual bee species abundances where development appeared in generalized linear model selection. Such findings underscore the variability in habitat associations among native bee species and are consistent with the idea that some taxa prefer open habitats (Harrison et al. 2018; Smith et al. 2021). With respect to the level of disturbance, agricultural areas fall between forests and developed areas. Similar to the results for forests, we found agriculture to correlate positively with bee diversity at larger spatial scales, and negatively at smaller spatial scales. Five of 12 bee species abundances where agriculture was significant in generalized linear models were positively affected, Bombus impatiens, Halictus ligatus/poeyi, Lasioglossum callidum, Lasioglossum imitatum, Lasioglossum tegulare/puteulanum. Landscape diversity (diversity of land cover types) generally correlated negatively with bee diversity and individual bee species abundances at larger scales while positive correlations were sometimes detected at smaller scales. Such findings may reflect the negative effects of development on bee communities at larger scales.

Opportunities for native bee conservation in human-dominated habitats can be better identified through increased understanding of the effects of anthropogenic activities at multiple spatial scales (Winfree et al. 2007; Bennett and Lovell 2019; Hall et al. 2016; Turo and Gardiner 2019; Wenzel et al. 2020). Wenzel et al. 2020 reviewed 140 studies to identify the drivers of urban pollinator populations and pollination. Positive responses were often associated with moderate levels of urbanization of rural, mostly agricultural land. Pollinator responses were commonly highly trait- specific. Cavity nesters and generalist species usually benefited from urbanization more than ground nesters and specialists.

A study of bee fauna and floral abundance in suburban yards in MA, USA identified 111 species active between May and September and emphasized the importance of dandelion (Taraxacum officinale), clover (Trifolium spp.) and other spontaneous lawn flowers as supplemental resources supporting pollinators (Lerman and Milam 2016). Research has demonstrated local level effects of increasing floral area and density suggesting that habitat gardening for bees can predictably increase abundance and diversity in urban settings (e.g., Frankie et al. 2009). Educational efforts in support of this concept have greatly increased in recent years (e.g., Griffin and Braman 2018). Pollinator plant preference for ornamental herbaceous and woody annual and perennial plants was the focus of some recent research efforts to assist consumers in plant selection, (e.g., Harris et al. 2016; Braman and Quick 2018; Erickson et al. 2019). Turo and Gardiner (2019) discuss the challenges associated with pollinator conservation in public green spaces. They advocate optimizing pollinators’ needs in a manner that both is economically feasible and respects societal norms and values when designing urban bee spaces, an inclusive process that should involve a diverse team of stakeholders.

Hamblin et al. (2018), determined that bee species richness and large bee abundance benefited from floral density, but cautioned that simply adding floral resources to otherwise hot, impervious sites was unlikely to evenly restore wild bee communities. Measuring conservation progress by species counts alone also does not address phylogenetic diversity and evolutionary history that may have important implications for ecosystem function and pollination services (Grab et al. 2019). While local site variables influence pollinators and their services, understanding influences of landscape variables at multiple spatial scales can inform regional urban planning to protect pollinators (Grab et al. 2019). Previous studies have found pollinator diversity to be both negatively and positively influenced by forest, agriculture, etc., but few have looked at species-specific traits which may differ from general measures of diversity.

When delving more deeply into the particular species within our surveyed bee pollinator communities, a handful of species appear to be driving the dissimilarities (Table 3, Online Appendix E). Researchers interested in these or other particular species can use our modeling results (Table 2) to determine what cover types and what scales were important drivers of species abundance and diversity in temperate suburban Oak/Pine landscape settings. Importantly, landscape scale, especially for forest influenced the direction of effect. We predicted that bees would be negatively affected by development and positively influenced by the amount of forest on the landscape. In fact some bees responded positively to development. We further predicted that these effects would be strongest at smaller scales. Forest at the small scale was always negative even for known forest associates that respond positively to forest at higher landscape scales, such as Osmia georgica. Bees that responded positively to forest habitat also included Andrena morrisonella, Augochlorella aurata, Bombus griseocollis, Lasioglossum callidum, Lasioglossum hitchensi, Lasioglossum illinoense, Lasioglossum tegulare/puteulanum and Melissodes communis. Those responding positively to more open environments included Bombus impatiens, Calliopsis andreniformis, Eucera (Peponapsis) pruinosa, Halictus ligatus/poeyi, Lasioglossum callidum, Lasioglossum hitchensi, Lasioglossum imitatum and Lasioglossum tegulare/puteulanum. Those responding positively to both included Lasioglossum callidum, Lasioglossum hitchensi and Lasioglossum tegulare/puteulanum suggesting an urban-forest interface association that merits further research.

Table 3 Indicator species by Community Cover type groups
Full size table

Conclusions

It is clear from this study that residential landscapes can support high bee diversity and that this diversity is sensitive to landscape context at different scales. Although development appears to have a negative effect on bee diversity overall, some bee species are favored by the open conditions characteristic of developed areas. Our study area was historically forested and forest remnants appear to be valuable habitats for many species and are thus important to regional bee diversity. NMDS ordination and PERMANOVA confirmed significantly different bee communities between sites with differing predominant land cover types (agriculture, forest and development). SIMPER and indicator species analysis revealed which species contributed heavily to the observed patterns and helped to determine group distinctions. Our results challenge the wisdom of overly generalized interpretation of landscape trends for urban pollinator communities, and show that bees respond to a complex assortment of landscape characteristics and this is driven by species-specific relationships with the land cover variables. These complexities highlight a need for further studying the interacting ecological drivers of pollinator communities at the urban-forest interface.