Edited by F. Stuart Chapin, University of Alaska, Fairbanks, AK, and approved November 7, 2011 (received for review September 21, 2011)
The phenology of many ecological processes is modulated by temperature, making them potentially sensitive to climate change. Mutualistic interactions may be especially vulnerable because of the potential for phenological mismatching if the species involved do not respond similarly to changes in temperature. Here we present an analysis of climate-associated shifts in the phenology of wild bees, the most important pollinators worldwide, and compare these shifts to published studies of bee-pollinated plants over the same time period. We report that over the past 130 y, the phenology of 10 bee species from northeastern North America has advanced by a mean of 10.4 ± 1.3 d. Most of this advance has taken place since 1970, paralleling global temperature increases. When the best available data are used to estimate analogous rates of advance for plants, these rates are not distinguishable from those of bees, suggesting that bee emergence is keeping pace with shifts in host-plant flowering, at least among the generalist species that we investigated.
Climate warming over the past 50 y is associated with phenological advances in a wide variety of organisms including plants, birds, and insects (1–3). Responses to climate warming are particularly important to understand for species that provide critical ecological functions such as pollinators. Furthermore, many ecological functions result from interactions among species, and because not all species respond to climate warming in the same manner, this could potentially lead to phenological mismatches that result in the loss of function (4–6). Alternatively, the interacting species may be buffered against climate variation if they have evolved similar responses to environmental variation (7). Here we present an analysis of climate-associated shifts in the phenology of wild bee pollinators, and compare the rates of advance for bees to those of bee-pollinated plants from the same region.
Bees (Hymenoptera: Apoidea: Anthophila) are the primary animal pollinators in most ecosystems (8). However, there is only one study of climate-associated phenological shifts in bees, and this focused on a single managed species, the honey bee (9). The honey bee (Apis mellifera L.) represents a special case relative to the ∼19,700 described species of bees existing worldwide (10), first, because it is a domesticated species, and second, because it is one of the minority of perennial bee species, meaning that adults remain active over the winter and regulate hive temperatures in temperate latitudes. In contrast, most wild bee species outside the tropics have annual cycles that include an obligatory larval or adult diapause before spring emergence. The development of bees and the environmental triggers regulating seasonal activity are largely unknown, and the few species studied show complex responses to both winter and spring temperatures (11, 12). Thus, although we would predict phenological shifts in bee activity due to climate change, the directionality and magnitude of these shifts are difficult to predict. In contrast to the minimal information available for bees, there is a large literature reporting phenological advances in plants (e.g., 13, 14). We used the published literature on native, bee-pollinated plants from northeastern North America to compare with rates of phenological advance between plants and bee pollinators.
To evaluate long-term phenological trends in wild bees from northeastern North America (Fig. S1), we used museum data dating back to the 1880s. Ten bee species that emerge in early spring were selected for study, because spring-active taxa are known to be good indicators of response to climate change (2). Additionally, species were selected to encompass a range of natural history traits, including both cavity and soil nesters, and both solitary and eusocial species. All 10 species are generalists that visit a wide range of spring-blooming flowers. A primary challenge in investigating long-term phenological shifts is obtaining reliable historical data. In contrast to other animal taxa such as birds and butterflies, for bees there are no long-term standardized monitoring schemes that could provide historical data. We used data from the contemporary period (2000–2010) in conjunction with data we recorded from museum specimens dating back to the 1880s to examine a 130-y period for which no other form of data on bee phenology is available (the final analysis comprises 3,447 records). Museum specimen records indicate that a species was in flight on the collection date, thus representing the span of activity for a given species in a given year. Such complete distributions may be more robust to sampling bias compared with records of the earliest activity in each year, and have been used in studies of climate-induced phenological change for other taxa (e.g., 15, 16).
To evaluate the rate of change in bee phenology over time, we used a general linear mixed model with the day of the year that the specimen was collected as the outcome, year, latitude, and bee sex as predictors, and bee species as a random factor. The overall model including all bee species showed a significant negative slope for collection date as a function of year (Fig. 1A; year estimate ± SE = −0.08 ± 0.01 d⋅y−1, P < 0.001; latitude estimate = 5.42 ± 0.26 d⋅degree latitude−1, P < 0.001; sex estimate = −18.13 ± 0.76 d, P < 0.001, pseudo-R2 = 0.45), indicating that spring-active bees are advancing their phenology at a rate of 0.8 d per decade. Over the entire time period studied, from 1880 to 2010, the regression predicts a mean advance of 10.4 ± 1.3 d. Most of this change has occurred during the last 40 y: The slope of bee advance between 1970 and 2010 is more than twice as steep (Fig. 1A; year estimate = −0.18 ± 0.05 d⋅y−1, P < 0.001; latitude estimate = 5.90 ± 0.37 d⋅degree latitude−1, P < 0.001; sex estimate = −18.38 ± 0.93 d, P < 0.001, pseudo-R2 = 0.43) as the overall slope, implying that 69% of the total advance, or roughly 7.2 ± 1.9 d, has occurred since 1970. When we analyzed the 1870–1970 period alone, the slope is indistinguishable from 0 (year estimate = −0.01 ± 0.02 d⋅y−1, P = 0.60; latitude estimate = 3.17 ± 0.48 d⋅degree latitude−1, P < 0.001; sex estimate = −16.68 ± 1.31 d, P < 0.001, pseudo-R2 = 0.40).
Rates of change in bee phenology and temperature over time. Each point represents a bee specimen used in the analysis. Raw data are shown without correcting for covariates (see text for details). The collection day for each bee specimen, where January first is day 1 (A), and the mean temperature during the month of April (in °C) at the location where the specimen was collected (B). The overall trend (black lines) and the trend from 1970 to 2010 (red lines) are shown. (C) The correlation between collection day and temperature. The overall trend is shown (black line). (D–G) Representative species for the genera in our analysis: A. miserabilis, C. inaequalis, O. lignaria, and B. impatiens queen. (Photos by J.S.A.)
The directional changes we found in bee phenology are likely due to recent climate warming, as has been demonstrated for other organisms (17, 18). Across all of the locations where bee specimens were collected, there is a positive association between mean April temperature and year, when corrected for latitude and longitude (Fig. 1B; year estimate = 0.01 ± 0.001 d⋅y−1, P < 0.001; latitude estimate = −1.30 ± 0.02 d⋅degree latitude−1, P < 0.001; overall model R2 = 0.70). As is the case for bee phenology, the slope of this relationship is greater for the last 40 y (year estimate = 0.035 ± 0.003 d⋅y−1, P < 0.001; latitude estimate = −1.30 ± 0.02 d⋅degree latitude−1, P < 0.001; longitude estimate = 0.06 ± 0.02 d⋅degree longitude−1, P < 0.001; overall model R2 = 0.74; see also ref. 19). Importantly, when bee specimen collection day was regressed directly against April temperature, there was a significant positive association (Fig. 1C; temperature estimate = −3.61 ± 0.15 d⋅°C−1, P < 0.001; sex estimate = −18.57 ± 0.81 d, P < 0.001, pseudo-R2 = 0.46). This supports the idea that climate change is the main factor explaining the observed phenological advances.
When bee specimen collection day was regressed against year for each bee species individually, all 10 species showed negative slopes. In four species this trend was significant (P < 0.05), whereas three other species showed a strong but nonsignificant trend (P < 0.1) (Table 1). Latitude and sex were significant for all species, with bees flying later further north and males of solitary bees flying before females (Table 1). No interaction was detected between sex and year, indicating that flight seasons are shifting earlier at similar rates for male and female solitary bees. Protandry is well-known in solitary bees (20). Thus, detecting both protandry and known latitudinal patterns indicates the appropriateness of our data and methodology for detecting unknown, climate-associated phenological shifts. Last, when bee specimen collection day was regressed directly against April temperature, there was a significant positive association in all species (Table S1). Although in principle bee phenology could respond to climate-associated changes in precipitation as well as temperature, only A. crataegi, the latest-emerging species, showed a model improvement when precipitation was added to the model (Table S1).
Bee phenology models: The best models for each bee species analyzed
Bee natural history traits including nest type and sociality had only weak associations with bee phenological patterns (SI Text, Influence of Bee Natural History Traits in Determining Rates of Advance for Bees). The only significant pattern to emerge was that the bee species that emerge earlier in the season (Fig. S2) are advancing their phenology more quickly (P < 0.001, R2 = 0.59; Fig. S3). This same pattern has been observed in plants (14, 21).
Eighty-seven percent of the world's angiosperm species are pollinated to some degree by animals (22), and recent work has emphasized the potential for negative consequences if phenological mismatches occur between plants and pollinators (23, 24). However, data on changes in pollinator phenology and on corresponding changes in the plants they pollinate are scarce. Two short-term observational studies provide evidence that a focal plant species and its pollinators may experience some mismatch when spring advances (25, 26), but this does not have to be the case (27). Here we used long-term data to compare phenological shifts for 10 bee species to shifts in 106 native plant species that are visited by these same bee species.
Four published studies report long-term data on phenological shifts in entomophilous plants within the broad geographical area of our bee dataset, and we used all four in our analyses. First, Primack and co-workers (15, 28) investigated rates of advance for multiple plant species in Massachusetts, using herbarium records collected between 1885 and 2003 and comparing these to bloom dates for the same species as observed in 2003. Second, Bradley et al. (29) recorded the date of first flowering of 24 plant species common in our study area from 1936 to 1999. They used data collected in Wisconsin, which is geographically outside but at a similar latitude as our study region. Third, Cook et al. (30) recorded the first flowering of 11 plant species in New York State from 1931 to 2008. Finally, Abu-Asab et al. (31) recorded the phenology of 44 plant species in Washington, DC from 1970 to 1999, which corresponds to the period of the greatest temperature changes, both in our bee dataset (Fig. 1B) and in larger climate studies from our region (19).
From each published study, we extracted data only for those plant species that are native to our study region, flower during the flight period of the bee species we studied, and are known to be visited by at least some of these same bee species (Table S2). We compared rates of advance for these plant species to the rates of advance for our 10 bee species as analyzed over the same time period. In all four cases, the mean advance of plants was not different from the mean advance of bees computed over the same time period (all t tests, P > 0.20; Fig. 2). In two studies (28, 31), plants were advancing faster than bees by 23% and 26%, whereas in the other two studies (29, 30), bees advanced faster than plants by 3% and 33%. Our results suggest that there is no clear pattern of divergence between bees and plants. However, it is difficult to know how much of this variability is due to real biological differences in phenology, as opposed to sampling error or methodological differences among studies. For example, in the three plant species that were reported by more than one study (29, 30), the measured rates of phenological advance differ by a mean of 0.14 d⋅y−1 (Fig. S4), which is greater than the mean rate of advance reported by most studies.
Comparison of bee and plant studies. Mean and 95% confidence interval (CI) for the rate of phenological advance for all bee species in combination (filled circles), compared with the mean and 95% CI for plants (hollow circles). The values for bees differ among comparisons because we used only the bee data from the same period over which plants were studied in each comparison. A dotted line is drawn between the bee and plant mean for purposes of visual comparison. From left to right, the plant studies are refs. 28–31.
The lack of a statistically significant difference does not demonstrate that the two groups have similar responses, especially given the large variation in the responses (Fig. 2). Thus, we used an equivalence test (32) to determine the size of the difference between the rates of bee and plant advance within which the observed rates can be considered statistically equivalent. These values ranged from 0.06 to 0.09 d⋅y−1 for the three plant studies covering long-term time periods to 0.17 d⋅y−1 for the plant study covering only the most recent period (Fig. S5). Those equivalence intervals suggest a maximum mismatch of 4–6 d over the ∼70- to 100-y span of the three long-term studies. However, for the most recent period, during which the greatest temperature changes occurred, it suggests a maximum mismatch of 5 d over only 30 y. The observed mismatches are likely to be of little biological significance, given that poor weather can easily prevent plant–pollinator interactions for several days in a row in early spring in our study region. However, if climate warming continues at its recent rapid rate, more significant phenological mismatches could occur in the future.
Overall, results indicate that phenological changes in bees have paralleled changes in the plants that they visit. First, both bees and plants responded to the temperature increases of 1971–1999 (Fig. 1B) by more than doubling their rates of phenological advance (Fig. 2), suggesting a parallel response to climate change. Second, the mean difference between plant and bee responses is small in most cases [Primack et al. (28), 0.024 d⋅y−1; Cook et al. (30), −0.027 d⋅y−1; Bradley et al. (29), −0.002 d⋅y−1; Abu-Asab et al. (31), −0.080 d⋅y−1; Fig. S5]. Third, the directionality of the difference varies among studies (Fig. 2 and Fig. S5). Interestingly, the two plant studies conducted inside cities showed greater advances, suggesting that the potential urban heat-island effect merits further study (33). Fourth, both bees and plants showed greater rates of advance for species that are active earlier in the season, suggesting some correspondence between taxa (14). Fifth, most species of plants and bees are generalists in terms of their interactions (34) (Table S2); thus, most species are not tied to the phenology of a small number of partners. Finally, for most species of plants and bees, an individual is in flower or in flight for at least several weeks. This period is longer than the maximum phenological mismatch predicted to occur by the equivalence test when taking into account the large variability observed. Overall, we conclude that phenological mismatch probably has not occurred already, but that it could occur in the future, as rates of temperature warming increase.
There are several important caveats to our findings. First, our entire approach assumes generalization on the part of both plants and pollinators, as we compared entire distributions rather than specific pairwise interactions. Analogous studies of specialist species are needed, particularly in ecosystem types that are known to be dominated by oligolectic bee species (35, 36) or where precipitation rather than temperature is the main driver for the timing of phenological activities (37, 38). Second, because the most significant climate warming has occurred during the last 40 y, caution is needed when extrapolating from existing data to predict future trends. There is a need to develop more mechanistic, physiological models for bees as well as plants if we want to make predictions beyond the range of our associational data.
Before beginning analysis, we selected the bee species for study according to several criteria. First, we chose species known to emerge in early spring for comparison with spring-flowering plants, because early spring-emerging species are the most responsive to climate change (14). Second, we selected species for which we could find sufficient historical specimen records in the principal northeastern North American entomology museums (see SI Text, Data Capture Protocol for a list of institutions accessed and data capture protocol). Finally, we included species with contrasting natural histories (Fig. 1 D–G), including solitary ground nesters (Andrena miserabilis Cresson, A. carlini Cockerell, A. crataegi Robertson, and Colletes inaequalis Say), primitively eusocial ground nesters (Bombus impatiens Cresson and B. bimaculatus Cresson), and solitary above-ground cavity nesters (Osmia bucephala Cresson, O. pumila Cresson, O. atriventris Cresson, and O. lignaria lignaria Say). For bumble bees, we restricted our analyses to queens. Bumble bee queens of the selected species emerge from diapause in early spring, slightly after the earliest-emerging solitary bees, and are therefore the appropriate caste to use for phenological analyses (see SI Text, Bumble Bee Analysis for details on other bumble bee castes). All 10 study species are widely distributed over eastern North America and are known to be important pollinators for both wild plants and crops (39–41).
To limit the geographic extent of our study area, we used records ranging from 36° to 50° N latitude and 85° to 70° W longitude (Fig. S1). These geographic limits were selected so as to minimize the extent to which underlying geographic variation might complicate an understanding of phenological effects, while simultaneously using in the analysis 80% of our original data records. In particular, Schwartz and Reiter (42) analyzed plant phenology across North America showing that despite differences between main geographical areas, the Northeast region presents a consistent response. To assure 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. The data used in the final analysis consisted of 3,447 specimens, collected by 763 collectors in 3,277 different collection events. Only 8 collectors contributed more than 50 specimens, and the top collector contributed only 232 specimens, making it unlikely that the behavior of a particular collector is biasing the sample.
Monthly climate data (monthly minimum and maximum temperature and overall precipitation) corresponding to the location and year where each specimen was collected were gathered from two sources. Both datasets consist of monthly climate surfaces for North America, generated by using weather station data along with spatial interpolation. We used data from the Canada Forest Service (CFS; http://cfs.nrcan.gc.ca/projects/3?lang=en_CA) (43) where possible, because it is spatially continuous, whereas the other data sources are provided at a coarser resolution. CFS uses thin plate spline smoothing algorithms (ANUSPLIN) to interpolate the data. The predictive accuracy of the historical surfaces was validated against independent climate station data, and mean errors for both temperature and precipitation were generally small (±0.5 °C for temperature) (43). In cases where CFS data were not available, we obtained the three climatic variables from PRISM (parameter-elevation regressions on independent slopes model; http://www.prism.oregonstate.edu) (44). PRISM is available at a 2.5-arc⋅min (∼4 × 4-km) resolution, which although coarser than the CFS data, is adequate for the large-scale questions we are asking. All data were processed in GRASS software (45). The CFS and PRISM datasets are highly correlated (R2 across our collection points = 0.92); thus, the use of the two datasets only minimally affects our results. CFS and PRISM both report minimum and maximum monthly temperature for a given location. These values were highly correlated for our locations (R2 = 0.83); thus, in our analyses, we used the mean of the maximum and minimum monthly temperature, and refer to this value hereafter as the mean monthly temperature. There were no data available for some data points, including all pre-1900, and hence temperature analysis had a slightly smaller sample size of 3,339 data points, compared with 3,447 data points for the analysis of bee specimen collection date versus year.
We conducted a joint analysis of all bee species combined by using R package nlme (46). This analysis used collection day as an outcome variable and the collection year, sex, longitude, latitude, and the interactions between collection year and sex as predictors. Genus and species nested within genus were included as random factors. Predictor variables were not strongly correlated (mean of absolute values of Pearson correlations = 0.22; maximum = −0.43). Diagnostic plots were examined to check for heteroscedasticity, as well as to ensure the normality of errors (47). We present here the best model based on Akaike's information criteria (AIC) (48). An R2 statistic cannot be computed on models including random factors; thus, we present pseudo-R2 values, which correspond to the squared correlation between the fitted and observed values. We performed a separate analysis using only the most recent 40 y of data, during which our study region has experienced accelerated climate warming (19). Because the choice of breakpoint can influence results, we did a sensitivity analysis on the choice of breakpoint year. Results were qualitatively similar when different breakpoints ranging from 1960 to 1980 were used. In all cases, the slope of the more recent period was steeper than the overall period, with the most recent slope ranging from −0.12 to −0.24 (Fig. S6). We report the 1970 breakpoint in the results to be consistent with the consensus of climatologists from our region (19). A variance component analysis of this complete model performed with R package ape indicates at which taxonomic level (species or genus) most of the variability in the response exists. This helps us understand how taxonomic relationships, which reflect phylogenetic relationships, might affect phenological responses (49). In a second analysis, we analyzed each bee species separately to investigate trends in phenology over time using linear models. We used the specimen collection day as the outcome variable and the fixed predictors mentioned above. To investigate differences in response by natural history group, slopes were compared between nesting substrate and sociality categories using t tests.
Whereas the temporal trends documented above represent a comprehensive analysis of multiple climatic factors, temperature is likely the main environmental driver of phenological shifts in plants and animals. Temperature data from our collection locations indicate a mean temperature increase of ∼1.4 °C since 1970 (Fig. 1B), similar to increases reported for other parts of northeastern North America (28, 31). Thus, we built an additional linear model by regressing the collection day directly against the climate data. Sex was included as a covariate, but not latitude, because it was highly correlated with temperature (Pearson correlation = 0.83). In the studied bees, emergence is known to depend primarily on early spring temperatures, corresponding primarily to April for the latitudes in our study region (Fig. S2; see ref. 12 for O. lignaria). However, plants may respond to other climatic cues, such as the accumulation of degree-days. The high correlation among monthly temperatures across sites prevents us from including other months in the model (i.e., Pearson correlation of March–April mean temperature = 0.80); thus, we present the regressions with April, the month that shows most explanatory power in all species based on AIC. We also investigated whether the model's explanatory power was improved by including April precipitation as an additional predictor. Throughout our analyses we used 95% confidence intervals, and thus assessed significance at the 0.05 level.
Data on the flowering phenology of plants occurring in the study area were obtained from four papers (28–31). Plant species were used only if they are native to our study region and start flowering between March and May, which corresponds to the activity period for the bees in our dataset. This information was gathered from the US Department of Agriculture PLANTS webpage (http://plants.usda.gov) and local plant guides (e.g., 50). Second, we only included plant species that are known to be visited by at least one of the 10 bee species that we analyzed, or in some cases for which the plant genus is known to be visited by at least one of our bee species (Table S2). To determine bee–plant associations, we used our own databases of more than 1,800 pairwise records, as well as the Flower-Visiting Insect Database (http://flowervisitors.info), which records pairwise records for 725 plant species and 463 bee species. Additional information was gathered from the published literature (51, 52).
Analyses were conducted in the following way. Data from Primack et al. (28) were reanalyzed following the same technique used by the authors after selecting the subset of plants meeting our criteria (114 observations of 27 species). For the other three studies (29–31), we were able to extract the rate of advance for each plant species that we selected, in units of d⋅y−1 that flowering period advanced. Overall, we analyzed 11 individual responses for Cook et al. (30), 24 for Bradley et al. (29), and 44 for Abu-Asab et al. (31). The bee analysis was performed in each case for the same time period that each plant study reported. In the case of Cook et al. (30) and Abu-Asab et al. (31), the median starting point of all plant responses (1,936 and 1,971, respectively) was used as the early threshold for the bee data. Our methodology for the bee data are most similar to that of Primack et al. (28), in that we both analyze the phenophase (i.e., specimens collected throughout the entire period of flight or throughout the peak bloom of the plant species), whereas the other three plant studies are based on the first reported flower in a given year.
We used t tests to determine whether plant and bee slopes, measured in units of days of advance per year, were significantly different. To further determine the difference in advancing rates that can be considered statistically similar, we used an equivalence test (32, 53). In general, when performing an equivalence test, the researcher defines an a priori minimum difference between datasets, the equivalence interval, which is then assumed as the null hypothesis (H0 = the distributions differ by more than the chosen equivalence interval). However, in our case, there is no biological basis for defining a meaningful difference; thus any threshold chosen would be arbitrary. Instead, we took an alternative approach and calculated the maximum equivalence interval that would be required to demonstrate that bee and plant slopes were statistically equivalent.
We thank all those who collected the bees used in our analyses, and the museums for access to their collections. We thank S. Droege and collaborators, Y. Mandelik, and L. Richardson for contributions to the contemporary dataset, and F. Gadallah for cleaning and error-checking the Bombus dataset. J. Gibbs conducted and supervised the databasing at Cornell University. Data capture at the American Museum of Natural History (AMNH), the University of Connecticut, Rutgers University, and Cornell University was supported by National Science Foundation Biological Infrastructure (DBI) Grant 0956388 (to J.S.A.), with additional support at AMNH from Robert G. Goelet and at the University of Connecticut by state wildlife Grant 09DEP10012AA (to D.W.). I.B. was funded by postdoctoral Fellowship EX2009-1017 from the Spanish Education Ministry 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., D.W., B.N.D., S.C., and S.K. performed research; I.B. and R.W. 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.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1115559108/-/DCSupplemental.
