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.