Methods

Study site

Our study sites were located in-and-around the Pennsylvania State Fruit Research and Extension Center in Adams County, Pennsylvania, USA (39.935226, -77.254530). The area has an average yearly rainfall of 112 cm, average summer temperature ranging from 16 °C to 28 °C, and average winter temperatures of -5 °C to 5 °C (Biddinger et al., 2018). This landscape is hilly with well-drained soils and the broader area is approximately 56% broadleaf forest fragments, 25% pastureland, 9% developed areas, and 8% commercial orchards (Biddinger et al., 2018). All orchards were managed under growers’ choice conventional pest management programs that use pesticide classes including insect growth regulators, anthranilic diamide, tetramic acid, microbials, and neonicotinoid insecticides (Biddinger et al., 2018). We sampled bees at 8 locations adjacent to 4 different active apple orchards. Sampling locations were within 150 m of orchards and 250 m of a forest fragment (Figure 1) which have diverse plant and pollinator communities (Kammerer et al., 2016). Often orchards rely, in part, on managed honey bee colonies for pollination, which have the potential to negatively impact native bee populations (Mallinger et al., 2017). However, our sampling sites did not have managed honey bee hives within 2 km and growers managing the adjacent orchards had not rented honey bees for at least 15 years.
Our bee monitoring traps were located within perennial wildflower strips approximately 50 m x 10 m in size that were sown between 2-3 years before the beginning of our study. Wildflower sites used in this study were established and managed using the specific planting guidelines developed by the Pennsylvania USDA-NRCS and the Xerces Society for Invertebrate Conservation (NRCS, 2011). They were sown with 21 species of native forbs and grasses sourced from a local native seed supplier (Ernst Conservation Seed, Meadville, PA 16335). All wildflower sites were mowed once a year and received spot sprays of common selective herbicides as needed.

Bee collections

We trapped bees continuously from April to October from 2014 to 2019 using Blue Vane traps (BanfieldBio Inc., Woodinville WA). A previous study in this region showed that Blue Vane traps collect a higher abundance and total richness of bees than colored bowl traps, also called pan traps (Joshi et al., 2015). Although the overall community composition of bees captured in Blue Vane traps was different from bowl traps, nearly all species were more likely to be captured in Blue Vane traps over bowls, except some Andrena and Lasioglossumspecies (Joshi et al., 2015). In our study, Blue Vane traps were filled with about 7 cm of 60% ethylene glycol (Supertech® Wal-Mart Stores, Inc., Bentonville, AR), hung from posts about 1.5 m off the ground. At each of our 8 locations, we placed 2 traps 25 m apart (Figure 1). Traps were left outside continuously from April to October every year and traps were replaced each year in case wear over time decreases their attractiveness. Each week, all specimens were removed and the ethylene glycol was replaced. Bee specimens were separated from other insects collected in the traps and stored in 70% alcohol until they were washed, pinned, and labeled. All bees were identified to the species level except 14 individuals were removed because of uncertain species-level identification. For bee identification, we used published dichotomous keys (Mitchell, 1960, 1962; Michener et al., 1994, Michener, 2000) and various interactive bee identification guides available at Discover Life (http://www.discoverlife.org). Species identifications were conducted by David Biddinger (Pennsylvania State University), Robert Jean (Senior Entomologist, Environmental Solutions and Innovations, Inc.), Jason Gibbs (University of Manitoba), and Sam Droege (United States Geological Survey). All specimens from this study are stored at the Pennsylvania State Fruit Research and Extension Center, Biglerville, PA, or the Frost Insect Museum at Pennsylvania State University, University Park, PA.

Data and statistical analyses

All analyses were conducted with R version 4.0.3. To explore how thorough our sampling of bee biodiversity was we created a species accumulation curve based on the number of individuals sampled. We did this using the ‘specaccum’ function in the vegan package (Oksanen et al., 2013).
We calculated biodiversity metrics for each of our 8 sampling locations using the combined data from both traps at each site. For month analyses (within years) we calculated biodiversity metrics for each month and then averaged across all years. We then averaged together sites closer than 900 m resulting in 4 replicates. Therefore, the interpretation of each replicate is the average biodiversity value per site for a single month. For year analyses (across years), we first summed data across all months and then calculated metrics for each year and then averaged to get 4 replicates per year. The interpretation of these data is then the total biodiversity value per site for a single year.
Using community abundance data, we measured total abundance, richness, and diversity (inverse Simpson’s) within and across years using the ‘specnumber’ and ‘diversity’ functions in the vegan package (Oksanen et al., 2013). We also used community abundance data to measure differences in community composition using a Bray-Curtis dissimilarity matrix. To measure phylogenetic structure, we use a genus-level molecular phylogeny from Hedtke et al. (2013). We made the phylogeny ultrametric with the ‘force.ultrametric’ function in the phytools package (Revel, 2012) using the non-Negative Least-squares method. We then amended our species below the genus-level using the ‘genus.to.species.tree’ function which creates bifurcating subtrees among species in each genus and then binds them at a random place along the terminal edge. With this approach, the relationships among species below the genus level are created at random. We only had to drop 2 species from these analyses because their genera (Cemolobus, Triepeolus ) were not in the tree, which only had 7 and 2 individuals in the dataset, respectively.
Using the community data and this phylogeny, we calculated mean pairwise distance (MPD) among all taxa. This metric is a measure of the average evolutionary distance between all pairs of species in a community and was calculated with the ‘mpd.query’ function in the PhyloMeasures package (Tsirogiannis and Sandel, 2016). Because raw values of MPD can be impacted by species richness of samples, we used a standardized effect size of MPD that is based on the difference between the observed measure of MPD and a random null model of a community with the same number of species (Tsirogiannis and Sandel, 2016). This value is also standardized by variance making the metric an effect size in standard deviation units. Negative values of this metric indicate a community is more clustered (less average evolutionary distance among pairs of species) than a random community with the same richness, and positive values indicate a more even community, also known as overdispersed. The resolution of phylogenies below the species level has little-to-no impact on MPD calculations. Qian and Jin (2021) found that MPD values are nearly identical when calculated with a genus-level phylogeny with species amended (like ours) compared to a fully resolved phylogenetic tree. And in our case, repeating the process of randomly adding species below the genus level resulted in nearly perfectly correlated measures of MPD (r = 0.99). However, trees without species-level resolution do not provide reliable estimates of phylogenetic signals of species’ traits (Davies and Kraft, 2012), so we did not include those analyses in the current study.
We modeled the changes in bee biodiversity within and across years by fitting general additive models using the ‘gam’ function in the mgcv package (Wood, 2017). These allowed for non-linear fits to the data. For all tests, we used 5 knots which allowed for sufficient curviness to represent observed patterns and to produce linear relationships between observed and fitted values. The amount of variation explained by these models was typically similar to ANOVA models fit to the same data. We report percent change effect sizes among extreme values of months and years as the difference in means between a pair of values (e.g. the mean from April minus the mean for August) then divided that difference by the overall mean for that variable. This gives a standardized effect size that can be compared among variables. We used perMANOVA to test differences in composition among groups using the ‘escalc’ function and visualized results with non-metric multidimensional scaling fit with the ‘metaMDS’ function, both from the vegan package (Oksanen et al., 2013).
To assess how phenological patterns differ among bee families and species, we calculated a metric for seasonality (time of year when the species showed the highest abundance) and phenological breadth (a measure of the amount of time bees are active as adults). Phenology measures were calculated for species with 30 or more individuals (a total of 40 species) as this is enough to reliably estimate phenological breadth (Bartomeus et al., 2013). We measured seasonality as the median day of year (Julian date) of capture across individuals collected for a species. Our measure of breadth was the difference between the 10th and 90th percentile of capture dates. To control for the different numbers of individuals across species, we randomly drew 30 data points for each species repeated 500 times. For each of these subsamples, we calculated the seasonality and breadth statistics and then averaged those 500 values to get our final statistics (Bartomeus et al., 2013).
To evaluate how bee families and species differ in their changes in abundance over time, we calculated a metric of change over time for each species as the slope of the linear relationship between abundance and year. We standardized abundance data for each species to have a mean of 0 and a standard deviation of 1 to allow comparisons across species. We multiplied the model coefficients by five (the change in years in our sampling) to get a predicted level of change in standard deviation units over the course of the study.