Bee Specimen Data.

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).

Changes in Aggregate Species Richness over Time.

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.

Changes in Relative Abundance.

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.

Ecological Traits Associated with Species Responses.

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.

Potential Bias Due to Unknown Sampling Effort over Time.

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.