Edited by May R. Berenbaum, University of Illinois at Urbana–Champaign, Urbana, IL, and approved February 1, 2013 (received for review October 24, 2012)
Pollinators such as bees are essential to the functioning of terrestrial ecosystems. However, despite concerns about a global pollinator crisis, long-term data on the status of bee species are limited. We present a long-term study of relative rates of change for an entire regional bee fauna in the northeastern United States, based on >30,000 museum records representing 438 species. Over a 140-y period, aggregate native species richness weakly decreased, but richness declines were significant only for the genus Bombus. Of 187 native species analyzed individually, only three declined steeply, all of these in the genus Bombus. However, there were large shifts in community composition, as indicated by 56% of species showing significant changes in relative abundance over time. Traits associated with a declining relative abundance include small dietary and phenological breadth and large body size. In addition, species with lower latitudinal range boundaries are increasing in relative abundance, a finding that may represent a response to climate change. We show that despite marked increases in human population density and large changes in anthropogenic land use, aggregate native species richness declines were modest outside of the genus Bombus. At the same time, we find that certain ecological traits are associated with declines in relative abundance. These results should help target conservation efforts focused on maintaining native bee abundance and diversity and therefore the important ecosystems services that they provide.
Pollination is an essential ecosystem function because 87% of the world’s angiosperm species are pollinated by animals (1), including most of the leading global food crops (2). Bees (Hymenoptera: Apoidea: Anthophila) are regarded as the most important pollinators, both for their efficiency and their ubiquity (3). However, despite concerns about pollinator declines and a global pollinator crisis (4), long-term data on the status of bee populations are scarce (5). Thus, a recent US National Academy of Sciences report concluded that “for most pollinator species, the paucity of long-term data and the incomplete knowledge of even basic taxonomy and ecology make definitive assessment of status exceedingly difficult” (6).
Heretofore, most studies reporting bee population declines have been focused on the bumble bee genus, Bombus. Some Bombus species are declining sharply in North America (7⇓–9) and elsewhere (10, 11), although others remain numerous (7) or are expanding their ranges (10). Furthermore, Bombus may not be representative of the world’s 442 other bee genera because they may have been impacted by recent pathogen introductions from managed Bombus colonies (7). Much less is known about the status of other genera, which account for >95% of the ∼20,000 described species of bees worldwide (12). The only long-term analysis to date that included species other than Bombus showed significant declines in estimated species richness in the United Kingdom and the Netherlands when comparing pre- versus post-1980 records (13). However, both study areas are characterized by particularly intensive land use and may not be representative of changes in the status of bees in other parts of the world. Thus, the existence of a widespread crisis in pollinator declines, as often portrayed in the media and elsewhere (4), rests on data of limited taxonomic or geographic scope.
Environmental change affects species differentially, creating “losers” that decline with increased human activity, but also “winners” that thrive in human-altered environments (14). Although there are likely winners and losers among bees, the identity of these species is largely unknown (6). In general, certain life-history traits are predicted to make species more vulnerable (15), such as having a small niche breadth (e.g., a specialized diet; 16). However, predictions for some other traits such as body size have resulted in contrasting predictions (17, 18). In bees, additional traits such as nest site location and brood parasitism or sociality also determine a species’ response to environmental change (19⇓–21).
Here we present a long-term study of relative rates of change for all 47 northeastern North American bee genera, comprising 438 species. To achieve a long-term (140-y) dataset, we databased, identified, and filtered >30,000 bee specimens from major collections of leading northeastern North American museums.
We first binned the 30,138 specimens into 10 time periods, each containing a similar number of independent records. For the non-Bombus species, we found that the number of rarefied bee species per time period has declined by 15%, but the trend is not significant (permutation test P = 0.07; Fig. 1A and Fig. S1A). Modest richness declines are in accordance with the few field data available (22, 23). For Bombus, species richness declined by 30% over the 140-y period (permutation test P = 0.01; Fig. 1B and Fig. S1B). Our result confirms previous studies documenting North American declines in Bombus species richness (7⇓–9). Last, we identified 20 exotic species (i.e., not native to the United States) in our study area (Table S1) and found that the number of exotic species collected increased by a factor of 9 over time (permutation test P = 0.01; Fig. 1C).
Trends in species richness over time. Number of bee species (± SE) in a standard number of independent specimen records per time period. Dashed line indicates a nonsignificant trend and solid line a significant trend. (A) All native bee species excluding Bombus (rarefied to 1,000 specimens). (B) Genus Bombus (rarefied to 400 specimens). (C) Exotic bee species (rarefied to 1,000 specimens).
A limitation of our study, and of previous published studies of long-term trends in pollinators (13), is that sampling effort is unknown, particularly in the earlier time periods. Although the rarefaction techniques we use partially correct for this, rarefied richness may reflect changes in species composition and dominance as well as changes in species richness; thus, our results should be interpreted as combining these two metrics. For example, one possible interpretation is that some species are becoming more common over time whereas other species are becoming less so.
In a second, species-level analysis, only three of 187 species, all in the genus Bombus, exhibited a rapid and recent population collapse [generalized linear model (GLM) using presence-absence as response variable; P < 0.05]. Two of these species, Bombus affinis and Bombus pensylvanicus, have recently been identified as being of conservation concern based on data independent of ours (7), whereas Bombus ashtoni, a social parasite of B. affinis and the related B. terricola, was also identified as declining in our study. One additional species, Macropis patellata, stood out as being of conservation concern, having shown a gradual, historical (pre-1950) decline, and has not been collected in the region since 1991 (Table S1). In contrast to the paucity of steep declines, there have been large shifts in relative abundance over time, as indicated by the fact that 29% of the species decreased significantly, whereas 27% had significant increases (GLM using year as the predictor and the proportion of records per year as the response variable, with the response weighted by the total number of specimens collected in that year; P < 0.05; Table S1). For example, although we did not detect a rapid and drastic decline for Bombus terricola, which has been previously identified as of conservation concern (7), we do find a significant decline in its relative abundance. Although our data show steep declines for this species in the southern part of its range, recent records demonstrate persistence, especially in northern areas within the core of the species’ historic range. Many other species with significant declining trends are common and are still collected regularly, but nevertheless should be monitored because such slow declines could be taken as an early warning signal of imperilment. Last, nine of the 87 rare species (defined as having 10 < n < 30 specimens; all of these species were excluded in our species-level analysis because of inadequate sample size) have not been recorded from the past 10 years (these are listed in Table S1) (24). Furthermore, our failure to detect particular species in the most recent period is conservative because the greatest collection effort took place during this period (Fig. 1).
There are several important caveats to our species- and genus-level results, which stem from limitations of our dataset. First, small sample sizes prevented us from statistically analyzing trends for the rarest species, which are potentially of the greatest conservation concern. Second, our analyses measure relative changes in abundance with respect to other species in the collections, not absolute changes. Third, some groups (e.g., Nomada) are particularly difficult to identify; the most taxonomically problematic species are noted in Table S1.
Investigation of ecological traits that differ between increasing and decreasing species can provide a more mechanistic understanding of how biodiversity is affected by environmental changes. In our study system, during the past 140 y, human population density more than doubled, from 140 people/km2 in 1900 to 325 people/km2 in 2000, with associated increases in human land use (25). We first assessed the phylogenetic signal (λ) of species relative rates of change itself, but found that although they do differ across bee genera, the strength of this signal at the genus level is not strong (26; λ = 0.24, P < 0.001; Fig. 2). The best-fit phylogenetic generalized least squares model (PGLS) based on the Akaike Information Criterion shows that species with a small dietary breadth, narrow phenological breadth, and large body size are more likely to be in decline (dietary breadth estimate = 0.008 ± 0.002, P = 0.001; phenological breadth estimate = 0.00009 ± 0.00004, P = 0.01; body size estimate = −0.003 ± 0.001, P = 0.01; Fig. 3 A–C). The body size result holds when Bombus, which are unusually large and also may be declining for particular reasons (7), are excluded from the analysis (body size estimate = −0.003 ± 0.001, P = 0.06). In addition, species with lower northern latitude range boundaries are increasing in relative abundance within our study area (northern latitude estimate = −0.002 ± 0.0004, P < 0.001; Fig. 3D), a finding that may represent a positive response to climate change on the part of these species (27). Mean April temperature has increased by more than 1° during the past 40 y in our study region and climate-associated phenological advances for bees have been already demonstrated (28); climate change may also be affecting bee population abundances and distribution. Brood parasitic species are predicted to be more vulnerable to decline because they are dependent on particular host bees (19, 29; but see ref. 30). However, we found no differences between parasitic and nonparasitic species in terms of their rates of decline. We likewise did not detect any influence of sociality (31; Fig. S2A) or nest site location (Fig. S2B) on rate of decline. Overall, ecological traits explained a small proportion of the variation (F = 13.61 on 180 DF, P < 0.0001; λ = 0, r2 = 0.22), indicating that a diversity of responses exists among species with those traits. Last, species that have increased significantly in relative abundance over the 140-y period include some known to tolerate human disturbance (32) as well as exotic species (Table S1).
Estimates of relative rate of change over the 140-y interval for different bee genera. Estimates are derived from logistic regression models. Negative values indicate declining trends. Boxplot width is proportional to the number of species sampled for each genus.
Relationships among species traits, phylogeny, and rates of change in relative abundance over time. The three most speciose genera are colored to demonstrate how traits are often shared among species within a genus: Bombus (blue), Andrena (red), and Lasioglossum (green). (A) Dietary breadth (floral specialization); boxplot width is proportional to number of species and individual species values are indicated along each side. (B) Body size measured as the intertegular distance (mm). (C) Phenological breadth (days of adult activity per year). (D) Northern distributional limit (degrees latitude).
We obtained our dataset by entering information from pinned bee specimens from the American Museum of Natural History and from collections maintained by the University of Connecticut, Cornell University, Rutgers University, Connecticut Agricultural Station, University of New Hampshire, University of Massachusetts, Vermont State Bee Database, New York State Museum, and the Bohart Museum of Entomology (University of California, Davis). To focus the geographic extent of our study area, we used records ranging from 38° to 45° N latitude and −85° to −70° W longitude (Fig. S3). This region has the most extensive historical collections of bees in the New World because of its early settlement and large number of universities. For each specimen, species identifications were made or verified by an expert taxonomist associated with our project (J.S.A. or J.G.), and full label data were captured into the Planetary Biodiversity Inventory database (33). Once databased, records were filtered, cleaned, and standardized in the following ways. We retained only specimens for which we had data on the locality and the year. Because the method used to collect bees in the field can affect the taxa collected (34), we excluded specimens known to be collected by methods other than hand netting. To ensure independence of samples, we used only one specimen of a given species from a given collection event, defined by unique combinations of collector, date, and location. Although this required removing ∼40% of the records from our analysis, it was an important step to minimize bias in our analyses. We removed taxa for which the species-level taxonomy is still unresolved (e.g., in the genus Nomada). We also excluded the honey bee Apis mellifera, because it is a managed species. For museum collections where complete databasing of all species was not possible, we avoided temporal bias in our dataset by deciding a priori which genera were most important to include (e.g., which had the smallest sample sizes in our existing dataset), and then fully databasing these genera across all time periods. We also excluded Bombus specimens held by Cornell University to ensure that our data were independent from those used in a previous study of Bombus declines (7). The data used in our final analysis consisted of 30,138 specimens, collected by at least 1,550 different known collector teams (with 3,708 specimens being collected by unknown collectors), in 11,295 different collection events.
Our final dataset included 438 species from 47 genera (Table S1). The frequency distribution of these is predictably skewed (35), such that only 80 species have more than 100 independent records, and 187 species more than 30 independent records. Fifty-four species are represented by a single specimen; however, because the number of singletons present show no trend by time period, they should not have a strong influence on measured changes in species richness over time (Fig. S4A). Species accumulation curves indicate that our databasing efforts captured nearly all of the species diversity present in the study area (36; Fig. S4B; see SI Text, Investigation of Potential Bias for investigation of potential bias in our data).
We compiled information on six ecological traits for the complete set of 187 species for which we ran individual species-level analyses (SI Text, Ecological Traits).
To investigate changes in bee communities over time, we first binned the data such that they had a similar number of independent specimen records per bin. We created the bins using quantiles from the overall cumulative distribution function such that each bin represented an equal number of quantiles. Hence, each bin comprises a similar number of records, but different numbers of years (36). This approach allowed us to maximize the sample size of each bin; it also minimized the influence of periods with low collection effort. To assure that the number of bins chosen did not affect our conclusions, we performed a sensitivity analysis on the number of bins used. Results were similar across the range of 3–10 bins (Fig. S1). To estimate species richness changes over time, we first rarefied all bins to a common number of specimens (1,000 independently collected individuals for all non-Bombus species and 400 for Bombus, representing 82% and 79% of the specimens in the smallest bin, respectively), and then calculated the mean ± SE number of species per time period (36). The resampling method allowed us to standardize sampling effort across time periods and gave us SE estimates, thus correcting for the fact that true sampling efforts in museum collections are unknown. This method provided a relative rather than an absolute measure of species’ abundances, because an apparent decline in some species could in fact be due to an increase in abundance of another species. As such, although it is standard in the ecological literature to refer to this measure as “species richness,” it actually measures both evenness and richness. To estimate statistical significance, we used a permutation test to randomly reorder the time periods and then calculated the correlation between time period and species richness across 1,000 such permutations. Bin size was not equal and hence we were unable to evaluate whether the rates of change reported are linear or take some other form. The reported P value is equal to the fraction of these permutations that had higher or lower correlations compared with the correlation we observed using the chronological time period sequence. Rarefaction analysis was done with the package vegan in R (37).
Given the different collection, identification, and databasing efforts that have been focused on Bombus (as indicated by the nonparametric multidimensional scaling plot; Fig. S5), we analyzed this genus separately. All other genera were analyzed altogether. Exotic species were excluded from analyses, except when explicitly noted.
To analyze relative trends in individual species over time, we selected all native bee species with more than 30 records (mean number of specimens per species = 143.2 ±11.7 SE) and investigated how the probability of finding a given species in the collection changes over time. For this analysis, we used a general linear model with a binomial distribution and a logit link. This analysis did not use resampling methods. Instead, each year was weighted by the total number of specimens collected that year to account for differential sampling effort among years. The logit link transforms the response to probabilities constrained between 0 and 1, and the directionality and magnitude of change is given by the model estimate. For species showing overdispersion, a quasibinomial distribution was used.
The analysis described is designed to find gradual, long-term declines. However, it has low power for detecting recent, steep declines. Thus, we used a second model for the latter purpose. The second model uses species presence/absence information for each year. This model is uninformative for species collected in most years, but can detect abrupt disappearances/extirpations. Again, Bombus was analyzed separately, using as the outcome variable the proportion of each Bombus species with respect to all Bombus collected in a given year. Likewise, exotic bee species were analyzed separately, using as the outcome the proportion of each exotic species with respect to all non-Bombus species collected in a given year. All other taxa were analyzed together.
To correct for phylogenetic nonindependence among traits, we inferred a phylogenetic tree (SI Text, Phylogenetic Tree; Fig. S6). We used this tree to determine whether ecological traits were correlated with relative change of abundance over time using a PGLS approach (38) as implemented in the R package Caper (39). Our outcome variable was the rate of change in relative abundance over time as estimated by the slope of the individual species models. Both the outcome variable and all traits included in the models show a significant phylogenetic signal (λ) when tested alone. Phylogenetic signal was calculated using packages Phytools (40) and Geiger (41) in R (Table S2). The variance-covariance between rate of change in abundance and each trait was then estimated simultaneously with λ and fitted to the generalized least squares models (42). The best model based on AIC is presented. Note that when the full model λ is estimated to be 0 (as the one we report), the estimate from the PGLS model is not identical to that of a standard linear model.
All museum data—including ours—are subject to potential biases if collectors change their methodology over time. We discussed this issue with data providers and the main contemporary collectors to ensure a proper interpretation of the data. For example, our dataset includes records from a diverse array of land use types across all time periods, but early collectors (19th and early 20th century) may have been more prone to creating synoptic collections in which each species is represented by only one specimen, whereas contemporary collectors are more likely to be conducting ecological studies in which their goal is to collect each species in proportion to its abundance. Although there is no definitive way to resolve this issue, we were able to explore and minimize potential collecting bias. First, we filtered the data to use only one specimen of each species per collection event, thus largely removing the “redundant” specimens of common species. Second, the fact that we had 1,550 collectors and that 88% of these each contributed fewer than 20 specimens makes systematic bias in collector behavior over time less likely. Third, if the bias hypothesized existed, it would bias the results toward finding steeper declines in species richness than in fact occur. This is a bias against our main finding of relatively weak and nonsignificant richness declines, at least for taxa other than Bombus.
We thank all those who collected the bees used in our analyses and the museums for access to their collections, J. Pickering for maintaining www.discoverlife.org, R. T. Schuh for access to the Planetary Biodiversity Inventory database, and D. Sol and D. Cariveau for comments on the manuscript. Data capture was supported by National Science Foundation Division of Biological Infrastructure Grant 0956388 (to J.S.A.), with additional support from Robert G. Goelet and a state wildlife grant [09DEP10012AA (to D.L.W.)]. This work was also supported by a Postdoctoral Fellowship from the Spanish Education Ministry Grant EX2009-1017 (to I.B.) and a Rutgers University Pre-Tenure Career Enhancement Award (to R.W., I.B., and J.S.A.).
Author contributions: I.B., J.S.A., and R.W. designed research; I.B., J.S.A., J.G., and D.L.W. performed research; B.N.D. and S.M.H. contributed new reagents/analytic tools; I.B. analyzed data; and I.B. and R.W. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
Data deposition: The data reported in this paper have been deposited in the DRYAD repository, http://dx.doi.org/10.5061/dryad.0nj49.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1218503110/-/DCSupplemental.
