Phenological data
We assembled digitally available specimen records from 220 herbaria across the United States (Note S1). Species names were harmonized using the Taxonomic Name Resolution Service 54. We excluded: conspecific specimens collected within the same GPS locations and the same dates (i.e., duplicates); specimens not explicitly recorded as bearing flowers; lacking GPS coordinates, dates of collection, or species-level identification; and specimens collected outside of the United States. For subsequent analysis, we selected species represented by at least 300 specimens to ensure that our model was computationally tractable and that we had sufficient sample sizes for estimating temperature responses in space and time. The net result was a sample of 1,038,027 specimens that included 1,605 species (Fig. S3). We used day of year (‘DOY’) of collection of each specimen as a proxy for flowering date. Because flowering spanned year-ends for many species, we accounted for the DOY discontinuity between December 31st and January 1st using an Azimuthal correction, whereby DOYs from the year prior become negative values 29.
Climatic data
Temperature conditions preceding and leading up to anthesis can mediate flowering time through their effects on developmental rates of preceding phenophases or by cuing floral development and anthesis. Accordingly, we used mean surface temperatures averaged over a standard period of three months 18,21,50,55 leading up to (and including) the mean flowering month for each species (hereafter ‘TMEAN’) as a predictor. For each collection site, we obtained monthly TMEAN time series (January 1901 – December 2021) at a 4-km2 spatial resolution from the Parameter-elevation Regressions on Independent Slopes Model (PRISM Climate Group, Oregon State University, http://prism.oregonstate.edu). We characterized each collection site by its long-term mean temperature (hereafter ‘TMEANNormal’), averaging observed TMEAN across all years between 1901 and 2021. Annual deviations from long-term TMEAN conditions (hereafter ‘TMEANAnomaly’) at each site and in each year were calculated by subtracting the TMEANNormal from the observed TMEAN conditions in the year of collection. Positive and negative TMEANAnomaly values respectively reflect warmer-than-average and colder-than-average years. TMEANNormal and TMEANAnomaly were uncorrelated irrespective of the latitudinal and elevational range spanned by a species (median r = −0.04), thus representing independent axes of climatic variation (Fig. S4). TMEANNormal spanned a wider temperature range than TMEANAnomaly for most species, with respective median 90% ranges of 9.9 °C and 3.7 °C (Fig. S5). Species occurring in cold climates tended to show later mean flowering dates than species occupying warmer regions (Fig. S6a); consequently, average TMEANNormal values were well above 0°C leading up to the mean flowering dates of all species in our data (Fig. S6b)
To assess how sensitivities varied across climatic gradients (see Analyses, below), we first characterized long-term precipitation and temperature at each site of collection using a Principal Component Analysis (PCA), with mean annual temperature normal (MATNormal), mean annual precipitation normal (PPTNormal), temperature seasonality, and precipitation seasonality as input features. We obtained precipitation (hereafter ‘PPT’) data from PRISM and calculated PPT and temperature seasonality for each collection site as the difference between the months with the highest and lowest PPT and mean temperature normal, respectively. We made PPT seasonality proportional to local levels of precipitation by dividing differences in maximum versus minimum monthly precipitation normal by PPTNormal at each site. The PCA identified 2 principal components accounting for more variance than its input features, jointly explaining 78% of observed variation. PC1 was associated with increasing PPT seasonality (36%), decreasing temperature seasonality (31%) and increasing MATNormal (28%) (Fig. S7). PC2 represented a gradient of decreasing PPTNormal (74%) and increasing temperature seasonality (22%).
Analyses
Estimating flowering time sensitivity to temperature over space and time—We estimated flowering time sensitivity to TMEANNormal and TMEANAnomaly using a Bayesian mixed-effect model. The model fitted species-specific intercepts and slopes and treated them as random effects sampled from ‘community-level’ distributions (defined by among-species mean and standard deviation of intercepts and slopes). This hierarchical structure improved estimation of parameters by using information and estimates from all species in the data. In turn, the Bayesian inference framework allowed for estimation of the correlations between TMEAN sensitivities over space and time and their differences for each species while propagating parameter uncertainty.
We used DOY for each observation i as a response, assuming a normal distribution with mean µi and species-specific standard deviation σsp:
We used weakly informative priors, with wide, 0-centered normal distributions for intercepts, slopes, and rate parameters for exponential distributions (used to obtain species-specific variances). For the variance-covariance matrix Σ, we used a Lewandowski-Kurowicka-Joe (LKJ) Cholesky covariance prior, with ŋ = 1 to allow for high correlations among parameters. Posterior distributions were obtained using Hamiltonian Monte Carlo (HMC) in Stan (code provided in Note S2) as implemented in R v.4.2.1 using the ‘rstan’ package v.2.21.2 56. We implemented a non-centered parameterization to improve sampling of the parameter space. Sampling was done using three MCMC chains with a training period of 1000 iterations and sampling of 4000 iterations. All Sspace, Stime, and Sspace − Stime estimates had Gelman-Rubin statistics (‘R-hat’) of less than 1.002, and visual examination of trace plots confirmed convergence.
Fitting the model on simulated data (Note S3), which emulated the average range of TMEAN conditions and the signal-to-noise ratio of DOY vs. TMEAN observed within species in our data, confirmed that our model could accurately recover the parameters of interest (Stime, Sspace, and Sspace − Stime ) for a range of sample and effect sizes (Figs. S8—10). Moreover, we found that apparent plasticity (Stime) and apparent adaptation (Sspace − Stime) could be estimated with similar degrees of precision (Fig. S11).
Because our model did not include an explicit temporal predictor, it may appear to ignore widespread trends in phenology and temperature reported in recent decades (19). However, additional simulation analyses (Note S4) showed that our model does account for temporal trends in phenology among species that experience trends in TMEANAnomaly over time and are responsive to TMEANAnomaly (i.e., non-zero Stime) (Fig. S12a). To evaluate the model’s implicit assumption that trends in TMEANAnomaly cause observed trends in phenology, we used the herbarium dataset to determine empirically whether observed temporal trends in TMEANAnomaly and a species’ Stime indeed explain observed trends in DOY. We recovered the same patterns observed in the simulation (Fig. S12b), suggesting that phenology and TMEANAnomaly trends are causally related. Moreover, detrending DOY and TMEANAnomaly prior to fitting the model did not affect our results, suggesting that omitting time as a covariate was unlikely to bias our results (Fig. S13).
Finally, we evaluated the impact on our estimates of choosing alternative reference periods to calculate TMEANNormal (i.e., 1901–2020 vs. 1901–1930, 1931–1960, 1961–1990, 1991–2020) (Note S5, Figs. S14–16). These analyses confirmed that period selection was unlikely to have affected our results.
Exploring biological assumptions—Herbarium specimens rarely are collected repeatedly at the same location across years. Accordingly, we found few repeated collections over time and in close enough proximity to represent single populations. Because of this, we estimated Sspace and Stime using statistical methods different from 9 and 23 (Note S6). Nevertheless, the interpretation of our model relied on the same simplifying assumptions: spatial slopes reflect variation among populations along a temperature gradient, temporal slopes reflect plasticity, plasticity does not vary within and among populations, and the temporal and spatial relationships between phenology and climate are not biased by confounding factors.
We evaluated the plausibility of many of these assumptions. Sspace likely represented phenological variation among populations because conspecific specimens were collected over vast regions spanning median 90% latitudinal and longitudinal ranges of 938 km and 1250 km, respectively. In turn, Stime likely reflected the effects of plasticity and not adaptation: analyses including only long-lived perennials (unlikely to show microevolutionary changes over short time periods) yielded very similar results to those presented below (Fig. S2); moreover, detrending DOY and TMEANAnomaly prior to fitting the model—which may account for temporal confounds or microevolution 57—yielded nearly identical estimates (Fig. S13). Furthermore, we generated a single estimate of Stime per species, thus assuming uniform plastic responses within and among populations. This assumption was supported by the observation that, for a large majority of species, Stime did not vary along geographic gradients of long-term TMEAN, long-term PPT, TMEAN seasonality, PPT seasonality, or the joint gradients described by PC1 and PC2 (Fig. S17). Cumulative precipitation and photoperiod are unlikely to confound Sspace and Stime: accounting for cumulative PPT yielded nearly identical estimates in single-species models (Fig. S18), and an analysis of 120 species collected withing geographic ranges restricted to narrower latitudinal bands (≤1°)—and limited geographic variation in photoperiod—yielded results very similar to those based on the entire dataset (Fig. S19). Finally, we detected no biases in Sspace or Stime due to differences in sample size among species (Fig. S20), spatial autocorrelation (Fig. S21), phylogeny (Fig. S22), non-linear phenology-temperature relationships (Fig. S23), or difference range size among species (Fig. S24).
Categorizing sensitivity patterns—To assess the prevalence of apparent plasticity and adaptation among species, we categorized each species’ Sspace versus Stime patterns as consistent with the effects of plasticity alone (Figs. 1a,b), adaptation alone (Figs. 1c,d), the joint effects of plasticity and adaptation (co- or counter-gradient adaptation; Figs. 1e–h), or neither. Classifications were based on the proportion of the posterior probability distribution of Stime and Sspace − Stime lying in the direction of their maximum a posteriori (MAP) estimate (i.e., their “probability of direction”, henceforth ‘PD’). PD is bound by 0.5 (maximum uncertainty about the effect of the predictor) and 1 (certainty of an effect in the direction of the MAP estimate). We subjectively considered apparent plasticity (Stime) and adaptation (Sspace − Stime) as significant when their PD was ≥ 0.95 (Table 1). Apparent plasticity and adaptation showed similar levels of estimation uncertainty both empirically (SD = 0.87 ± 0.34 d/°C for Stime; SD = 0.93 ± 0.32 d/°C for Sspace − Stime) and in preliminary simulation analyses (Note S3), suggesting sensitivity patterns were not substantially more likely to be classified as consistent with plasticity than with adaptation (and vice versa).
Table 1—Criteria for classifying the sensitivity pattern of each species as consistent with the role of plasticity only, adaptation only, the joint effects of plasticity and adaptation in a co- or counter-gradient adaptation pattern, or neither adaptation nor plasticity. The probability that Stime or Sspace − Stime differed from 0 in the direction of its maximum a posteriori (MAP) estimate (i.e., their probability of direction) was obtained from the posterior distribution of these parameters for each species.
Biological Process
|
Empirical Sensitivity Pattern
|
Plasticity only
|
- Probability of direction for Stime ≥ 0.95
- Probability of direction for Sspace − Stime < 0.95
|
Adaptation only
|
- Probability of direction for Sspace − Stime ≥ 0.95
- Probability of direction for Stime < 0.95
|
Plasticity and Adaptation
|
Co-gradient
|
- Probability of direction for Stime ≥ 0.95
- Probability of direction for Sspace − Stime ≥ 0.95
- Sspace and Stime have the same direction
- |Sspace| > |Stime|
|
Counter-gradient
|
- Probability of direction for Stime ≥ 0.95
- Probability of direction for Sspace − Stime ≥ 0.95
Case 1:
- Sspace and Stime have opposite direction
Case 2:
- Sspace and Stime have the same direction
- |Sspace| < |Stime|
|
Neither
|
- Probability of direction for Stime < 0.95
- Probability of direction for Sspace − Stime < 0.95
|
Plasticity and adaptation across phenological niches, local climates, and ecoregions—To assess how apparent plasticity and adaptation varied with native climate and phenological niche among species, we first calculated the mean flowering DOY and the mean coordinates along the climate gradients described by PC1 and PC2 among specimens of each species. We then fit two generalized additive models (GAMs) using Stime or Sspace − Stime as responses—assumed to be normally distributed—and a three-variable tensor-product smooth of mean flowering DOY, mean PC1, and mean PC2 as a predictor. This design allowed us to assess how native climate and phenological niche jointly determined the apparent roles of plasticity and adaptation while accounting for possible interactions and non-linearities. Because Stime and Sspace − Stime are estimated quantities, we accounted for parameter uncertainty by weighting each observation by the inverse of its posterior variance (i.e., its precision).
Additionally, we assessed the relative contributions of apparent plasticity and adaptation throughout the season within ecoregions of the contiguous United States. To do so, we identified the Level II Ecoregion—as classified by the USA Environmental Protection Agency (EPA) 58,59—within which each specimen was collected. We used Level II Ecoregions because they provide sufficient ecological detail to distinguish regional floras while encompassing areas broad enough for each to capture multiple species in our data. To avoid inflating species overlap among regions or the influence of species that were rarely sampled within an ecoregion, we arbitrarily considered a species as present within an ecoregion if at least 10% of its collections occurred within it. We then retained only ecoregions represented by a minimum of 8 species. Under this scheme, the median species was classified as occurring within 2 ecoregions (range = 1–7), the median ecoregion was represented by 156 species (range = 17–956 for Atlantic Highlands and Western Cordilleras, respectively), and pairs of ecoregions shared, on average, 4% of their species (range = 0–39%; Fig. S25). Of the 120 ecoregion pairs examined, 57 shared less than 1% of species, 100 shared less than 10% of species, and 114 shared less than 20% of species.
Once species × ecoregion combinations were identified (n = 3,570), we fitted two GAMs including apparent plasticity (Stime) or apparent adaptation (Sspace − Stime) as a response, ecoregion as a categorical predictor, mean flowering DOY as continuous predictor, and a mean flowering DOY × ecoregion spline assessing the ecoregion-specific effects of mean DOY on apparent plasticity or adaptation. Again, we accounted for parameter uncertainty by weighting each observation according to the precision of its corresponding apparent plasticity or adaptation estimate. Collection locations in different ecoregions differed substantially in their long-term climatic conditions (Fig. S26). However, we assumed no intraspecific variation in Stime across ecoregions an assumption partially supported by the observation that Stime did not tend to vary along climatic gradients within species (Fig. S16). All GAMs were implemented using the ‘mgcv’ package v.1.8-40 in R 60,61 .