Plasticity and not adaptation is the primary source of temperature-mediated variation in flowering phenology in North America

Phenology varies widely over space and time because of its sensitivity to climate. However, whether phenological variation is primarily generated by rapid organismal responses (plasticity) or local adaptation remains unresolved. Here we used 1,038,027 herbarium specimens representing 1,605 species from the continental United States to measure flowering-time sensitivity to temperature over time ( S time ) and space ( S space ). By comparing these estimates, we inferred how adaptation and plasticity historically influenced phenology along temperature gradients and how their contributions vary among species with different phenology and native climates and among ecoregions differing in species composition. Parameters S space and S time were positively correlated ( r = 0.87), of similar magnitude and more frequently consistent with plasticity than adaptation. Apparent plasticity and adaptation generated earlier flowering in spring, limited responsiveness in late summer and delayed flowering in autumn in response to temperature increases. Nonetheless, ecoregions differed in the relative contributions of adaptation and plasticity, from consistently greater importance of plasticity (for example, southeastern United States plains) to their nearly equal importance throughout the season (for example, Western Sierra Madre Piedmont). Our results support the hypothesis that plasticity is the primary driver of flowering-time variation along temperature gradients, with local adaptation having a widespread but comparatively limited role. The timing of life-cycle events (phenology) determines the environmental conditions that organisms encounter throughout development and often mediates their fitness 1 . Phenology usually is cued by seasonally and interannually variable climatic factors—such as temperature—that enable individuals to adjust growth and reproduction plastically in response to fluctuating environmental conditions 1,2 . Phenology also varies within species as a result of evolutionary adaptation to local environments, which may

Phenology varies widely over space and time because of its sensitivity to climate.However, whether phenological variation is primarily generated by rapid organismal responses (plasticity) or local adaptation remains unresolved.Here we used 1,038,027 herbarium specimens representing 1,605 species from the continental United States to measure flowering-time sensitivity to temperature over time (S time ) and space (S space ).By comparing these estimates, we inferred how adaptation and plasticity historically influenced phenology along temperature gradients and how their contributions vary among species with different phenology and native climates and among ecoregions differing in species composition.Parameters S space and S time were positively correlated (r = 0.87), of similar magnitude and more frequently consistent with plasticity than adaptation.Apparent plasticity and adaptation generated earlier flowering in spring, limited responsiveness in late summer and delayed flowering in autumn in response to temperature increases.Nonetheless, ecoregions differed in the relative contributions of adaptation and plasticity, from consistently greater importance of plasticity (for example, southeastern United States plains) to their nearly equal importance throughout the season (for example, Western Sierra Madre Piedmont).Our results support the hypothesis that plasticity is the primary driver of flowering-time variation along temperature gradients, with local adaptation having a widespread but comparatively limited role.
The timing of life-cycle events (phenology) determines the environmental conditions that organisms encounter throughout development and often mediates their fitness 1 .Phenology usually is cued by seasonally and interannually variable climatic factors-such as temperature-that enable individuals to adjust growth and reproduction plastically in response to fluctuating environmental conditions 1,2 .Phenology also varies within species as a result of evolutionary adaptation to local environments, which may select for different mean phenological timings among or within populations in space and time [3][4][5][6] .Although both plasticity and adaptation alter phenology, their relative contributions rarely have been measured within the same system largely because doing so requires experiments or spatiotemporally extensive genetic sampling [7][8][9] (but see ref. 6).Accordingly, most studies have highlighted either plasticity or adaptation as mechanisms of phenological variation attributable to environmental change 7 but their relative importance across species and ecological contexts remains unresolved.Elucidating https://doi.org/10.1038/s41559-023-02304-5across a climatic gradient, we should observe negligible sensitivity to temporal climatic variation (that is, no plasticity; S time = 0) and a substantial difference between the slopes of the temporal and spatial relationships (S space − S time ≠ 0 attributable to adaptation along the gradient; Fig. 1a,b).Alternatively, a phenologically plastic species whose populations are not locally adapted along the gradient should show sensitivity to interannual climatic variation (that is, S time ≠ 0) and no differences between temporal and spatial slopes (S space − S time = 0; Fig. 1c,d), implying that variation along the gradient can be attributed to plastic responses (that is, S space = S time ).When both adaptation and plasticity drive phenological variation along the climate gradient (that is, S time ≠ 0 and S space − S time ≠ 0), the resulting empirical pattern should depend on the relative direction of plastic and adaptive responses.Specifically, when adaptation operates in the same direction as plasticity (that is, co-gradient adaptation), we should observe a greater spatial than temporal sensitivity (for example, S time < 0 and S space − S time < 0 implies that S space < S time , so S space is more negative; Fig. 1e,f).In turn, when adaptation operates in the opposite direction to plasticity (that is, counter-gradient adaptation 15,16 ), we should observe a lesser spatial sensitivity or one of opposite direction to the temporal relationship (for example, S time < 0 and S space − S time > 0 implies that S space > S time , so S space is either less steep or positive; Fig. 1g,h).Finally, if a species shows no plasticity or local adaptation along a climate gradient, we would expect negligible temporal and spatial sensitivities (Fig. 1i,j).
Phenological sensitivity to temperature often varies among species occurring in different regions or that initiate phenological events at the degree to which species have phenologically responded to historical climatic variation through plasticity or adaptation could provide important context for predicting whether organismal responses may be sufficient-or evolutionary change necessary-to maintain development synchronized with suitable climatic conditions in a warming world 8 .
Phillimore and others 9 proposed that the relative and joint contributions of plasticity and local adaptation to spatial variation in phenology within a species can be estimated from the difference between the slopes of spatial and temporal phenology-climate relationships.This proposition rests on several observations.The effects of interannual climatic variation on phenology generally reflect plastic responses, especially among long-lived species less liable to experience micro-evolutionary changes from year to year 10 .Phenological variation over space also can be caused by phenotypic plasticity for which, for example, growing-degree day thresholds that trigger life-cycle events occur on different dates across sites 11 .However, among populations, local adaptation also can generate phenological variation along climatic gradients 12,13 .Therefore, assuming no confounding factors and in the absence of substantial variation in phenological plasticity within and among populations, phenological variation along spatial climate gradients should reflect the joint effects of plasticity and adaptation 14 .
Given these observations and assumptions, plasticity and adaptation can generate five empirical patterns of sensitivity to temporal climatic variation (hereafter, S time ) and to spatial climatic variation (hereafter, S space ) (Fig. 1).First, if a species does not show phenological plasticity but population-level phenological means are locally adapted  e,f, Plasticity and adaptation operating in the same direction (for example, both negative) (e) should result in a clear temporal relationship and a spatial relationship of substantially greater magnitude (f).g,h, In contrast, plasticity and adaptation operating in opposite directions (for example, plasticity negative, adaptation positive) (g) should result in a clear temporal relationship and a spatial relationship of substantially lesser magnitude (or having a different sign altogether) (h).i,j, Species exhibiting no plasticity or adaptation along the gradient (i) would generate negligible temporal and spatial slopes (j).Orange lines in a, c, e and g illustrate phenological responses of spatially separated populations to temporal temperature variation, which spans a narrower temperature range than spatial temperature variation across the entire species range (segmented red lines).The biological processes in a, c, e and g generate the empirical patterns in b, d, f and h.In turn, the empirical patterns imply the processes that generated them.See 'Exploring assumptions' for an overview of the assumptions of this approach and the degree to which they were met by our data.For examples of species exhibiting each of these patterns, see Supplementary Fig. 1. https://doi.org/10.1038/s41559-023-02304-5 different times throughout the growing season [17][18][19][20][21][22][23][24] .However, comparisons of phenological sensitivity to climate over space and time-which are necessary to evaluate the apparent contributions of plasticity and adaptation (Fig. 1)-across species differing in phenology and occupying different climates require spatiotemporally extensive datasets and therefore remain rare.Herbaria provide abundant and increasingly available data to conduct these analyses at unprecedented taxonomic, temporal and spatial scales 21,[25][26][27][28][29][30] .However, few studies have separately estimated sensitivity to spatial versus temporal climate variation using specimens (but see refs.28,31-36) and none has leveraged their unique scope to determine the ecological contexts in which plasticity or adaptation might contribute more strongly to spatial variation in phenology.
Here we analysed a dataset of over a million flowering specimens from 1,605 species across the continental United States to compare phenological sensitivities to spatial and temporal variation in temperature (S space and S time , respectively).For each species, we assessed whether its empirical sensitivity patterns were consistent with the effects of plasticity, adaptation or both along temperature gradients (Fig. 1).Additionally, we evaluated how apparent temperature-related plasticity and adaptation of flowering time varied among species with different native climates, phenological niches and occurring within different regional floras.Together, our analyses identified ecological contexts in which plasticity or adaptation appears to have most strongly influenced spatial phenological variation, providing the most taxonomically and geographically extensive assessment of temperature-mediated variation in flowering time among North American angiosperms conducted to date.

Plasticity versus adaptation as determinants of phenology
S space and S time of 93% and 79% of species, respectively, differed from 0 with at least 95% probability.S space and S time agreed in direction for 94% of species and estimates of both S time and S space were negative for 89% and 91% of species, indicating earlier flowering across increasingly warmer locations and in warmer-than-average years (Fig. 2a).
Both apparent plasticity and adaptation were associated with clinal variation in flowering time along temperature gradients, with plasticity playing a predominant role among species.S space and S time were highly positively correlated and their magnitude tended to correspond one-to-one for many species (Fig. 2b).Therefore, flowering shifts in warmer-than-average years typically had similar direction and magnitude (d °C −1 ) to those observed across increasingly warmer locations, consistent with a scenario of plasticity as the cause of phenological variation along spatial temperature gradients (Fig. 1c,d and Table 1).
More species showed sensitivity patterns consistent with plasticity (79%) than with adaptation (45%) (Fig. 1 and a detailed classification scheme in Table 1).Apparent plasticity explained ~52% of the variance in flowering-time clines along temperature gradients among species (straight black line in Fig. 2b).Of the species, 41% showed sensitivity  Patterns were classified 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 S time or S space − S time differed from 0 in the direction of its maximum a posteriori (MAP) estimate (that is, their probability of direction) was obtained from the posterior distribution of these parameters for each species. https://doi.org/10.1038/s41559-023-02304-5 patterns consistent with plasticity as the sole driver of phenological variation across gradients.In contrast, only 7% of species showed sensitivity patterns consistent solely with adaptation (Fig. 1a,b).Of the species, 38% showed both apparent local adaptation and evidence of plasticity.Among these, a greater proportion showed flowering advances (and co-gradient patterns; 27%) than flowering delays (and counter-gradient patterns; 10%) resulting from apparent adaptation along temperature gradients (Fig. 2b).Of the species, 14% showed patterns that were consistent neither with temperature-related plasticity nor with adaptation.These patterns remained consistent when analysing only long-lived species (whose responses to yearly temperature anomalies are certain to be plastic) (Extended Data Fig. 1).

Plasticity and adaptation across ecological contexts
Apparent plasticity (S time ) varied substantially among species with different phenological niches and across local climates (R 2 = 0.55; Fig. 3a,c).Species flowering during late winter and spring tended to show flowering advances in warmer-than-average years.Such advances decreased in magnitude throughout the season, typically reversing to flowering delays during late summer and autumn (Fig. 3a,c).The timing of the transition from positive values was consistent throughout PC1 (Fig. 3a) but occurred much earlier in arid regions with high temperature seasonality along PC2 (Fig. 3c).Apparent adaptation (S space − S time ) also varied with phenological niche and native climate (R 2 = 0.47; Fig. 3b,d).Apparent adaptation varied from negative to positive values throughout the growing season, indicating a transition from flowering advances to delays attributable to local adaptation.Such transitions occurred much earlier in cool, thermally seasonal regions (that is, the low range of PC1) (Fig. 3b).Apparent adaptation also varied throughout the growing season along PC2, with transition from advances to delays under warmer conditions occurring earlier in regions with high precipitation (Fig. 3d).
These patterns were mirrored at the regional level: throughout the season, average apparent plasticity and adaptation among species transitioned from generating flowering advances to generating delays in response to higher temperatures in all sampled ecoregions (R 2 for S time = 0.44; R 2 for S space − S time = 0.35; Fig. 4).This transition invariably occurred during the summer months.The magnitude of apparent adaptation tended to be lower than that of apparent plasticity during most of spring and early summer for all ecoregions.Their difference tended to be less among species flowering during early spring and the magnitude of adaptation was often greater among species flowering during late summer and early autumn (Fig. 4a-n).Nonetheless, we detected regional differences in the relative contributions of apparent adaptation and plasticity among species throughout the season.For example, apparent adaptation and plasticity had similar magnitudes within the Western Sierra Madre Piedmont (Fig. 4g).In contrast, mean apparent plasticity was consistently greater than adaptation among species in the southeastern United States plains (Fig. 4j).The difference in magnitude between apparent plasticity and adaptation was greatest among early-to mid-summer flowering species in the Western Cordilleras and cold deserts (Fig. 4b,c).

Discussion
This study provides evidence that, for 1,605 North American plant species, phenotypic plasticity historically has been the primary mechanism generating flowering-time variation along temperature gradients.Nonetheless, apparent adaptation and plasticity jointly generated phenological variation in many species.Both apparent plasticity and adaptation consistently generated flowering advances in spring, lesser advances during summer and flowering delays during early autumn and this pattern was consistent across climates and ecoregions.Whether phenological reaction norms to historical climatic conditions will remain adaptive under future climatic regimes is unclear 10 .Nonetheless, these results suggest that plasticity historically has enabled flowering phenology to respond quickly to a wide range of temperature conditions among North American angiosperms, with adaptation frequently playing an important but context-dependent role.

Plasticity causes clinal variation in flowering time
Extensive research has documented phenological plasticity to spatial climatic variation in plants [37][38][39][40] that can result in clinal phenological variation 41 even among short-lived taxa 11 .Our study extends these results by showing that the predominance of plasticity over adaptation associated with temperature-related variation in phenology over space appears to be the norm among North American species.
The greater importance of plasticity found in this study does not contradict the well-established role of phenological adaptation in space and time 40 , which can mediate rapid temporal shifts in phenology 5 or facilitate ecological invasions 6,42 .Indeed, 45% of species in our data showed evidence of adaptation-driven phenological variation along temperature gradients (Fig. 2b).It is also possible that we did not detect nonlinear or patchy adaptation patterns or that the contributions of apparent adaptation and plasticity may be different in regions under-represented in our data (for example, the Great Plains and prairies; Extended Data Fig. 2).Crucially, we only assessed the apparent contributions of plasticity and adaptation to observed variation in flowering time over temperature gradients, so our results do not rule out the possibility that adaptation is the primary driver of phenological variation along gradients of different climatic variables.Finally, determining the exact environmental conditions within microsites where herbarium specimens were collected is impossible because continental-scale climate products have relatively coarse spatial resolution and because specimen coordinates typically are inexact.Climatic variation at the microsite level could confound our estimates of S space and our assessment of the prevalence of local adaptation if, for example, different populations along the gradient occupied https://doi.org/10.1038/s41559-023-02304-5 distinct microsites that maintained temperatures more constant than apparent when looking at coarser pixel-level averages.However, to our knowledge, such microsite sorting across species ranges has only been reported at their trailing edges where climate is most limiting 43 .Nonetheless, these potential complexities underscore the ultimate need for molecular or quantitative genetic studies to corroborate the broad correlational patterns outlined in this study.
Still, the strong correlation between S space and S time has important implications for phenoclimatic research.For example, it suggests that temperature-related variation in flowering time among conspecific populations is a good proxy of responsiveness to interannual temperature variation.Therefore, space-for-time substitutions might be viable approaches for quantifying plastic flowering responsiveness to temperature in North American angiosperms, for most of which we lack long-term phenological records 26,44 .Specifically, the match between S space and S time shows that substituting space for time might reveal the direction and approximate magnitude on flowering sensitivity to temperature over time within species or relative differences in sensitivity among species.However, co-gradient adaptation frequently generated spatial sensitivities of greater magnitudes than those over time, demonstrating that S space might overestimate S time in many species.
Our results also indicate that plasticity may have generated phenological variation across a temperature range (a median range of 13.7 °C) exceeding the degree of warming forecasted for most regions in coming decades.However, such historical plastic flowering shifts over space will not necessarily be mirrored by temporal shifts within populations as warming trends continue.For example, historical temperature cues may become uncorrelated from the factors mediating the fitness consequences of phenology, rendering plastic reaction norms maladaptive 10 .Plastic phenological shifts associated with warming also may be constrained by physiology 45 or by other competing cueing mechanisms such as photoperiod or winter chilling that may be disrupted by phenological shifts associated with higher temperatures [46][47][48] .These complexities highlight the need for research on the fitness consequences of recent and ongoing phenological shifts 49,50 and on the inter-related mechanisms underpinning associations between multiple abiotic cues (for example, chilling, forcing, photoperiod and resources) and seasonal development beyond model systems 48,51 .

Plasticity and adaptation vary across ecological contexts
Sensitivities transitioned from flowering advances under warming in spring to reduced or no responsiveness during summer and even flowering delays in early autumn (Figs. 3 and 4).This pattern implies that temperature trends will probably drive changes to the structure of the flowering season during spring and autumn under global change but that other environmental factors might play predominant roles during summer.
These results support studies showing decreases in phenological sensitivity to temperature among species throughout the season in temperate biomes 18,21,52,53 and others showing flowering delays among autumn-flowering species or lengthening of the growing and flowering seasons under warming 23,[54][55][56] .While we cannot unambiguously identify the causes of this pattern, studies have shown that warming typically advances phenology during spring due to accelerated developmental rates, while phenophases occurring during autumn are cued directly by seasonal cooling [57][58][59] .This difference would explain why autumn-flowering species showed phenological delays under warming (that is, autumn cooling occurs later in warmer-than-average years) or why the transition from advances to delays was more pronounced within cool regions with high temperature seasonality (that is, those showing more pronounced cooling during autumn; Fig. 3).Regardless of its causes, our study corroborates that transitions from spring flowering advances to autumn delays because of climatic warming are consistent across thousands of species and diverse climate zones and biomes in the continental United States.
Likewise, apparent adaptation throughout the season typically transitioned from generating mean flowering advances to generating delays along temperature gradients.Our results are consistent with https://doi.org/10.1038/s41559-023-02304-5those reported by Delgado and others 23 , who found changes in the direction of apparent plasticity and adaptation throughout the growing season for multiple trophic levels (that is, saprotrophs, primary producers and primary and secondary consumers) in Eastern Europe.That changes in apparent plastic and adaptive responses to warming throughout the year might be robust across different phenophases, taxa, trophic levels or climatic regimes across the temperate zone may reflect shared cueing mechanisms or selective pressures for different phenological events occurring during the same seasons 56 , with factors other than temperature (for example, resources or photoperiod) probably driving phenological variation for developmental events in summer.Additionally, the greater prevalence of co-gradient adaptation as opposed to counter-gradient adaptation suggests that adaptation typically operates to generate greater variation in phenology along temperature gradients than is generated by plasticity alone.

Conclusions
Our findings indicate that phenotypic plasticity is the predominant historical mechanism of spatial phenological variation across a wide range of temperature conditions in the continental United States; adaptation plays more context-specific roles.Whether and how species-level attributes such as functional traits and life history may mediate these relative contributions or whether historical responses will tend to be adaptive under non-analogue climatic conditions remain open questions and important directions for future research.Our results outline broad correlational patterns whose verification will require direct measurements of plasticity and adaptation across species and climate regions.Nonetheless, our data-across many biomes and thousands of species-confirmed patterns of plastic and adaptive phenological advances in spring and delays in autumn in response to warming observed in detailed empirical studies, highlighting the increasing utility of biological collections for studying plant responses to global change at vast taxonomic and spatiotemporal scales.

Specimen data
We assembled specimen records from 220 herbaria made available digitally through 16 consortia from Mexico, the United States and Canada (accessed during July and August of 2022; Supplementary Note 1).We retained only specimens explicitly recorded as bearing flowers, which we determined by summarizing all unique entries in the DarwinCore 'reproductiveCondition' column and identifying those that unambiguously indicated presence of flowers.After harmonizing species names using the taxonomic name resolution service 60 , we removed specimens lacking species-level identification, GPS coordinates or dates of collection.To match the spatial and temporal coverage of the climate data (see 'Climatic data'), we retained only specimens collected from 1896 to 2020 within the United States.We considered as duplicates any conspecific specimens collected within 111 m (0.001 of a decimal degree) of one another on the same date.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.This filtering yielded a sample of 1,038,047 specimens from 1,605 species (Extended Data Fig. 2) (see ref. 61 for additional methodological detail).
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 31 December and 1 January using an azimuthal correction, whereby DOYs from the year before 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 cueing floral development and anthesis.Accordingly, we used mean surface temperatures averaged over a standard period of 3 months 18,21,53,62 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 1896 to December 2020) at a 16 km 2 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, TMEAN Normal ), averaging observed TMEAN across all years between 1896 and 2020.Annual deviations from long-term TMEAN conditions (hereafter, TMEAN Anomaly ) at each site and in each year were calculated by subtracting the TMEAN Normal from the observed TMEAN conditions in the year of collection.Positive and negative TMEAN Anomaly values, respectively, reflect warmer-than-average and colder-than-average years.TMEAN Normal and TMEAN Anomaly were uncorrelated irrespective of the latitudinal and elevational range spanned by a species (median r = −0.04),thus representing independent axes of climatic variation (Supplementary Fig. 2).TMEAN Normal spanned a wider temperature range than TMEAN Anomaly for most species, with respective median ranges of 13.7 and 5.4 °C (Supplementary Fig. 3).Species occurring in cold climates tended to show later mean flowering dates than species occupying warmer regions (Supplementary Fig. 4a); consequently, average TMEAN Normal values were well above 0 °C leading up to the mean flowering dates of all species in our data (Supplementary Fig. 4b).
To assess how sensitivities varied across climatic gradients (see 'Analyses'), we first characterized long-term precipitation and temperature at each site of collection using a principal component analysis (PCA), with mean annual temperature normal (MAT Normal ), mean annual precipitation normal (PPT Normal ), 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 PPT Normal at each site.The PCA identified two principal components (PCs) 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 MAT Normal (28%) (Extended Data Fig. 2).PC2 represented a gradient of decreasing PPT Normal (74%) and increasing temperature seasonality (22%).

Analyses
Estimating apparent plasticity and adaptation.We estimated flowering-time sensitivity to TMEAN Normal and TMEAN Anomaly using a Bayesian mixed-effects 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 modelled µ i as a linear function of TMEAN Normal (TMEANnorm i ) and TMEAN Anomaly (TMEANanom i ) for each observation i: https://doi.org/10.1038/s41559-023-02304-5 For each species (sp), the model yielded intercepts representing mean flowering dates (α sp ), sensitivities (that is, regression slopes) for TMEAN normal (S space sp ) and sensitivities for TMEAN anomaly (S time sp ).
To assess the correlation between S space and S time , we modelled community-level distributions for intercepts and slopes as generated by a multivariate normal distribution with a vector of hypermeans µ and a variance-covariance matrix Σ: We also calculated the difference between sensitivity types (S space sp − S time sp ) as a derived quantity within the model, which we interpreted as the degree of apparent local adaptation in DOY observed across the TMEAN normal gradient (Fig. 1), with negative and positive values, respectively, indicating advances and delays in flowering DOY across warmer locations.
We used weakly informative priors, with wide, 0-centred 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 Cholesky covariance prior, with ŋ = 1 to allow for high correlations among parameters.Posterior distributions were obtained using Hamiltonian Monte Carlo in Stan (code provided in Supplementary Note 2) as implemented in R v.4.2.1 using the rstan package v.2.21.2 (ref.63).We implemented a non-centred parameterization to improve sampling of the parameter space.Sampling was done using three MCMC chains with a training period of 1,000 iterations and sampling of 4,000 iterations.All S space , S time and S space − S time estimates had Gelman-Rubin statistics (R-hat) of <1.002 and visual examination of trace plots confirmed convergence.
Fitting the model on simulated data (Supplementary Note 3), which emulated the average range of TMEAN conditions and the signal-to-noise ratio of DOY versus TMEAN observed within species in our data, confirmed that our model could accurately recover the parameters of interest (S time , S space and S space − S time ) for a range of sample and effect sizes (Supplementary Note 3 and Supplementary Figs.5-7).Moreover, we found that apparent plasticity (S time ) and apparent adaptation (S space − S time ) could be estimated with similar degrees of precision (Supplementary Fig. 8).
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.However, additional simulation analyses (Supplementary Note 4) showed that our model does account for temporal trends in phenology among species that experience trends in TMEAN Anomaly over time and that are responsive to TMEAN Anomaly (that is, non-zero S time ) (Supplementary Fig. 9a).To evaluate the model's implicit assumption that trends in TMEAN Anomaly cause observed trends in phenology, we used the herbarium dataset to determine empirically whether observed temporal trends in TMEAN Anomaly and a species' S time indeed explain observed trends in DOY.We recovered the same patterns observed in the simulation (Supplementary Fig. 9b), suggesting that phenology and TMEAN Anomaly trends are causally related.Moreover, detrending DOY and TMEAN Anomaly before fitting the model did not affect our results, suggesting that omitting time as a covariate was unlikely to bias our results (Extended Data Fig. 3).
Finally, we evaluated the impact on our estimates of choosing alternative reference periods to calculate TMEAN Normal (that is, 1901-2020  versus 1901-1930, 1931-1960, 1961-1990 and 1991-2020) (Supplementary Note 5 and Supplementary Figs.10-12).These analyses confirmed that period selection was unlikely to have affected our results.
Exploring 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 S space and S time using statistical methods different from ref. 9 and ref. 23 (Supplementary Note 6).Nevertheless, the interpretation of our model relied on the same simplifying assumptions: spatial slopes reflect variation in DOY 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.Parameter S space probably represented phenological variation among populations because conspecific specimens were collected over vast regions spanning median latitudinal and longitudinal ranges of 1,356 and 1,819 km (removing outliers), respectively.In turn, S time probably reflected the effects of plasticity and not adaptation: analyses including only long-lived perennials (unlikely to show micro-evolutionary changes over short periods) yielded very similar results to those presented below (Extended Data Fig. 1); moreover, detrending DOY and TMEAN Anomaly before fitting the model-which may account for temporal confounds or micro-evolution 64 -yielded nearly identical estimates (Extended Data Fig. 3).Furthermore, we generated a single estimate of S time per species, thus assuming uniform plastic responses within and among populations.This assumption was supported by the observation that, for most species, S time 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 (Extended Data Fig. 4).Cumulative precipitation and photoperiod are unlikely to confound S space and S time : accounting for cumulative PPT yielded nearly identical estimates in single-species models (Extended Data Fig. 5) and an analysis of 120 species collected withing geographic ranges restricted to narrower latitudinal bands (≤1°)-and therefore to limited geographically driven variation in photoperiod-yielded results very similar to those based on the entire dataset (Extended Data Fig. 6).Finally, we detected no biases in S space or S time due to differences in sample size among species (Extended Data Fig. 7a,b), phylogeny (Extended Data Fig. 7c,d), spatial autocorrelation (Extended Data Fig. 7e,f), nonlinear phenology-temperature relationships (Extended Data Fig. 8) or difference in range size among species (Extended Data Fig. 9).
Although herbarium data have many spatial and temporal collection biases and limitations-including preferential collection near roads and urban areas and decreases in collection intensity in recent decades 65 -such biases are probably not severe in our data (Supplementary Notes 7 and 8 and Supplementary Figs.13-20).Our estimates of S space , S time and S space − S time were robust to inclusion in our models of factors such as urbanization (Supplementary Fig. 14) and proximity to major roads (Supplementary Figs. 17 and 18) and showed no evidence of various forms of temporal non-independence (Supplementary Fig. 20).Collector preferences can result in over-representation of certain taxa or traits among specimens 65 .While we cannot rule out these biases in our data, our study encompassed species from 106 families and 740 genera, capturing vast functional, evolutionary and life history diversity.Therefore, we consider it unlikely that our results were driven by over-representation of taxa or traits.Finally, some herbaria obscure location data for endangered or heavily poached species.However, since we only included georeferenced specimens from well-represented species-of which only 12 (0.7% of the total) are listed as endangered by the United States Department of Agriculture 66 -it is unlikely that our species list includes many such taxa.
Categorizing sensitivity patterns.To assess the prevalence of apparent plasticity and adaptation among species, we categorized each species' S space versus S time patterns as consistent with the effects of plasticity alone (Fig. 1a,b), adaptation alone (Fig. 1c,d), the joint effects of plasticity and adaptation (co-or counter-gradient adaptation; Fig. 1e-h) or neither.Classifications were based on the proportion of the posterior probability distribution of S time and S space − S time lying in the https://doi.org/10.1038/s41559-023-02304-5direction of their MAP estimate (that is, 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 (S time ) and adaptation (S space − S time ) as significant when their PD was ≥0.95 (Table 1).Apparent plasticity and adaptation showed similar levels of estimation uncertainty both empirically (s.d.= 0.87 ± 0.34 d °C −1 for S time ; s.d.= 0.93 ± 0.32 d °C −1 for S space − S time ) and in simulation analyses (Supplementary Note 3), suggesting sensitivity patterns were not substantially more likely to be classified as consistent with plasticity than with adaptation (and vice versa) due to estimation uncertainty.

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 S time or S space − S time 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 nonlinearities.Because S time and S space − S time are estimates, we accounted for parameter uncertainty by weighting each observation by the inverse of its posterior variance (that is, 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 67,68within 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 eight species.Under this scheme, the median species was classified as occurring within two 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%; Supplementary Fig. 21).Of the 120 ecoregion pairs examined, 57 shared <1% of species, 100 shared <10% of species and 114 shared <20% of species.
Once species × ecoregion combinations were identified (n = 3,570), we fitted two GAMs including apparent plasticity (S time ) or apparent adaptation (S space − S time ) 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 by the precision of its corresponding apparent plasticity or adaptation estimate.Collection locations in different ecoregions differed substantially in their long-term climatic conditions (Extended Data Fig. 10).However, we assumed no intraspecific variation in S time across ecoregions, an assumption partially supported by the observation that S time did not tend to vary along climatic gradients within species (Extended Data Fig. 4).All GAMs were implemented using the mgcv package v. 1.8-40 in R (refs.69,70).

Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.https://doi.org/10.1038/s41559-023-02304-5Data Fig. 1 | Distributions of and relationship between S space and S time among 201 long-lived species in the continental United States.Light blue and red shaded regions in (a) respectively correspond to the kernel-density distributions of S space and S time among the 201 species included in this analysis.The solid black line in (b) indicates a 1:1 relationship corresponding to perfect agreement between sensitivity types.The solid curved line indicates the line of best fit obtained from a Generalized Additive Model (GAM) of S time vs. S space , with the shaded area around it denoting the standard error of the predicted mean value.Each point in (b) represents a species whose x, y coordinates are given by the maximum a posteriori (MAP) estimates for S space and S time , respectively.Point shapes and colours in (b) indicate whether sensitivity patterns were consistent with plasticity or adaptation as the sole drivers of flowering time variation along the temperature gradient, with both plasticity and adaptation having statistically significant effects in a co-or counter-gradient adaptation pattern, or not showing statistically significant adaptation nor plasticity.The straight, solid black line in (b) indicates a 1:1 relationship (that is, S space = S time ), whereas the curved solid line shows the observed relationship estimated from a generalized additive model (GAM).The shaded region along the curved solid line in (b) corresponds to the standard error of the predicted value of S time .The percent of species showing each pattern is shown in the legend in parenthesis.The 95% credible interval for the correlation between S space and S time is provided as a text inset in (b).The subset of 201 species was selected based on growth form data from the United States Department of Agriculture Plant Database (USDA Plant Database, https://plants.usda.gov).We downloaded all records of growth habit information available through the search tool and subset the resulting dataset to contain only species represented among the 1,605 species included in the analyses presented in the main text.We then retained only flowering specimens from species classified as 'Tree' (n = 5), 'Shrub' (n = 164), 'Subshrub' (n = 27) or 'Vine' (n = 5), which yielded a dataset of 201 species.Using this subset dataset, we ran the model presented in the main text, obtaining estimates of sensitivity to TMEAN Normal and TMEAN Anomaly and of their difference for each species, as well as an estimate of their correlation accounting for parameter uncertainty.The resulting patterns closely mirrored those of the larger dataset, with a high correlation and agreement in magnitude between S space and S time and similar relative frequencies among species for each sensitivity pattern.Extended Data Fig. 3 | Estimates with and without detrending DOY and TMEAN Anomaly .Iler et al. 64 showed that shared temporal trends between DOY and temperature can generate spurious relationships between these variables that often disappear when the phenological and temperature time series are detrended prior to estimating their relationship.Alternatively, a non-spurious but trended relationship between DOY and temperature might reflect the effects of adaptation to directional changes in temperature, at least in short-lived species.Therefore, relationships between phenology and temperature that persist following detrending are more likely to reflect phenological plasticity.Accordingly, we assessed whether estimates of sensitivity to TMEAN Anomaly (S time ) presented in the main text could be confounded by temporal trends in DOY and TMEAN Anomaly .To do so, we first ran single-species linear regressions using DOY or TMEAN Anomaly as responses and year as a single predictor, storing the resulting residuals as detrended versions of both responses.Then, for each species, we ran two linear models of DOY against TMEAN Normal and TMEAN Anomaly : one with observed DOY and TMEAN Anomaly and another with detrended DOY and TMEAN Anomaly .Trended and detrended estimates of sensitivity to TMEAN Anomaly were very highly correlated among species, suggesting that TMEAN sensitivity estimates presented in the main text do not reflect the confounding effect of shared temporal trends.Similarly, detrending DOY and TMEAN Anomaly did not substantially alter estimates of S space and S space − S time .a-c, In each panel, points represent the combinations of trended or detrended estimates of S space , S time , or S space − S time for each species in the data, whereas diagonal black lines correspond to 1:1 to relationships denoting perfect agreement between trended and detrended estimates.Solid blue lines in each panel indicate the observed relationship between trended and detrended estimates, with the shaded region around the trend line (nearly imperceptible due to the large sample size) indicating the standard error of the predicted value.

Extended Data
Extended Data Fig. 4 | Variation in S time across geographic climatic gradients for 1,605 species across the coterminous United States.Distribution of interaction terms between TMEAN Anomaly and long-term climatic conditions within sites of specimen collection, including (a) TMEAN Normal , (b) PPT Normal, (c) TMEAN Seasonality, (d) PPT Seasonality, (e) the gradient of increasing temperature and precipitation seasonality described by PC1 and (f) the gradient of decreasing precipitation and increasing temperature seasonality described by PC2.The interaction coefficients for all variables were obtained from single-species models including flowering DOY as a response, the focal longterm climatic variable and TMEAN Anomaly as a predictor and an interaction term between them.Long-term climatic variables were standardized (mean = 0, SD = 1) before fitting the models.Accordingly, the interaction terms quantify the change in the slope of TMEAN Anomaly vs. DOY (S time , in days/°C) for an increase of 1 SD in the long-term climatic variable.The values for the 25 th , 50 th and 75 th percentiles of each among-species distribution are indicated in each panel text insets as well as the proportion of species for which the interaction coefficient had a p-value greater than 0.01 (based on a two-sided t-test).Within a species, phenological sensitivity to temperature can vary among portions of the range with different long-term climatic conditions.Therefore, differences between S space and S time presented in the main text may result from variation in the DOY-TMEAN Anomaly slope across the geographic climatic gradients and not from the effects of adaptation across the gradient.Despite this, we found no evidence of pervasive variation in S time along various geographic climatic gradients, suggesting intraspecific variation in phenological sensitivity is unlikely to generate the patterns reported in the main text.
Extended Data Fig. 5 | Consequence of including precipitation in models estimating S space and S time .Relationship between estimates of flowering sensitivity to TMEAN Normal (S space ) (a, c) and TMEAN Anomaly (S time ) (b, d) with or without accounting for the effects of cumulative precipitation normal and anomaly (x-axis and y-axis, respectively) during the same 3-month periods used to calculate TMEAN Normal and TMEAN Anomaly for each species (see Methods in main text).Panels (a) and (b) show the relationships between estimates from temperature-only models with those obtained from models including PPT normal and net PPT anomaly for the focal 3-month period in the year of collection.In turn, panels (c) and (d) show the same relationship but with estimates from a model including PPT normal and PPT anomaly proportional to the long-term average for that period (that is, divided by the PPT normal).Proportional anomalies were included to account for differences in the biological significance that the same amount of precipitation might have in chronically dry compared to chronically wet locations.The method developed by Phillimore et al. 9 assumes that the variables causing phenological variation along spatial temperature gradients are correctly identified and included in the model.Although temperature has been found to be a predominant environmental cue inducing flowering in temperate biomes, other variables such as precipitation, or those that emerge from the interaction between temperature and precipitation, such as snow cover or water stress, routinely have been implicated in phenological variation in many North American species.Therefore, it is possible that differences in spatial vs. temporal patterns of temperaturerelated phenological variation might stem from the confounding effects of phenologically important variables not included in our models.Estimates of phenologically important variables such as the timing of snowmelt or the onset of drought conditions in xeric environments are not available at the temporal and spatial scales spanned by our data.However, most of these variables are highly correlated with precipitation and temperature over space and time and including both in phenoclimatic models might account for the effects of predictors other than temperature and precipitation.Accordingly, we assessed whether estimated S space and S time changed when accounting for the effects of long-term cumulative precipitation (PPT normal) and PPT anomalies in the year of collection, separately assessing the effects of both net PPT anomalies and of anomalies scaled proportionally to long-term means (PPT normal) for the focal 3-month period.Estimates of phenology-temperature relationships in space and time did not change substantially when including precipitation variables, resulting in a very high correlation between estimates from temperature-only models and those from models including precipitation (r = 0.95 or 0.96).Therefore, the estimates presented in the main text are unlikely to be biased by the omission of precipitation during the months leading up to flowering.
Extended Data Fig. 6 | Distribution of and relationship between S space and S time among species sampled within narrow latitudinal bands.These analyses included 157 species with 200 or more specimens collected within a latitudinal band of 1° (~111 km) in the continental United States (analogous to Fig. 2 of the main text).Light blue and red shaded regions in (a) respectively correspond to the kernel-density distributions of S space and S time among the 157 species included in the analysis.The solid black line in (b) indicates a 1:1 relationship corresponding to perfect agreement between the two types of sensitivity.The solid curved line indicates the line of best fit obtained from a Generalized Additive Model (GAM) of S time vs. S space , with the shaded area around it denoting the standard error of the predicted mean value.Each point in (b) represents a species whose x, y coordinates are given by the maximum a posteriori (MAP) estimates for S space and S time , respectively.Point shapes and colours in (b) indicate whether sensitivity patterns were consistent with plasticity or adaptation as the sole drivers of flowering time variation along the temperature gradient, with both plasticity and adaptation having significant effects in a co-or countergradient adaptation pattern, or not showing statistically significant adaptation nor plasticity.The straight, solid black line in (b) indicates a 1:1 relationship (that is, S space = S time ), whereas the curved solid line shows the observed relationship estimated from a generalized additive model (GAM).The shaded region along the curved solid line in (b) corresponds to the standard error of the predicted value of S time .The percent of species showing each pattern is shown in the legend in parenthesis.The 95% credible interval for the correlation between S space and S time is provided as a text inset in (b).Both temperature and photoperiod are known to be the predominant environmental cues controlling both vegetative and reproductive phenology among plants in temperature regions.Therefore, across latitudinal ranges such as those spanned by most species in our data (median latitudinal range = ca.12.2°), it is possible that differences in S time and S space (for example, geographic temperature gradients) might reflect the confounding influence of latitudinal shifts in photoperiod on our estimates of sensitivity to TMEAN Normal .To account for this possibility, we identified 157 species in our data that were well sampled (200 or more specimens) within narrow latitudinal bands (≤1°).Using this subset of species and including only specimens from such 1° bands, we ran the model presented in the main text, obtaining estimates of S time and S space and their difference for each species and an estimate of their correlation accounting for parameter uncertainty.The results did not qualitatively differ from those presented in the main text, with a high correlation between S space and S time and similar relative frequencies of each sensitivity pattern among species.
Extended Data Fig. 7 | Effects of sample size differences, spatial autocorrelation and phylogeny on estimates of S space and S time .Comparison of S space and S time estimates obtained by (a) homogenizing sample sizes among species, (b) accounting for spatial autocorrelation among observations and (c) accounting for phylogenetic relationships among species against estimates generated ignoring these factors (as those presented in the main text).In (a, b), we fit the model presented in the main text using a thinned dataset were each species was represented by 300 specimens, comparing its output to that of the model in the main text.In (c, d), we compared the results of models omitting or accounting for phylogenetic relationships.We selected a random subset of 300 species from which to generte a phylogeny, thinning these data to include only 300 specimens for each species (to make the model computationally tractable).S space and S time estimates that did not account for phylogeny were obtained using the model described in the main text.In turn, the model accounting for phylogeny included a prior for the covariance structure of species-specific parameters consisting of the evolutionary distance between each pair of species as estimated from a phylogenetic hypothesis and a model of trait divergence among species.The phylogenetic tree (or hypothesis) was generated using the R package 'v.PhyloMaker' version 0.1.0 71and generated a phylogeny resovled to the genus level.Using this tree, we then calculated the variance-covariance phylogenetic matrix predicted by a Brownian model of trait evolution using the R package 'ape' version 5.6-2 72 .Finally, both models were implemented using the 'brms' package version 2.18.0 73 .Finally, in (e, f) we compared estimates obtained from models ignoring or accounting for spatial autocorrelation of the residuals.All S space and S time estimates were obtained using single-species models, but those accounting for spatial autocorrelation included a covariance structure for the residuals determined by the geographic distance between each pair of points.All models were fitted using the'nlme' package version 3.1 74 in R. Estimates of S space and S time obtained accounting for or ignoring spatial autocorrelation were nearly indentical across species.Across panels, the x-axes show the estimates obtained when omitting the focal factor (sample size, phylogeny, or spatial autocorrelation), whereas the y-axes show estimates obtained when accounting for it.Solid black lines represent a 1:1 line, representing perfect agreement in magnitude and direction between estimates.S space and S time estimates obtained ignoring sample size differences, phylogeny and spatial autocorrelation where highly correlated to estimates obtained from models accounting for these factors.Accordingly, we consider it unlikely that omitting these factors could have biased our results.
Extended Data Fig. 8 | Assessing evidence for nonlinear phenologytemperature relationships.Comparison of R 2 values obtained using 10-fold cross-validation of models of flowering DOY versus TMEAN Normal and TMEAN Anomaly obtained from (a) linear regressions assuming linear relationships between phenology and temperature or (b) generalized additive models (GAMs) accounting for potential nonlinear relationships.The shaded region in each panel represents the among-species kernel distribution of cross-validated R 2 values obtained using each model type (linear regression or GAM).The mean and SD of R 2 values each are presented as text insets in each panel.The model that generated the sensitivity estimates presented in the main text assumed linear relationships between flowering dates and TMEAN Normal and TMEAN Anomaly .To verify whether such an assumption was warranted for our data, we compared the predictive ability of single-species models assuming linear relationships between phenology and temperature (fitted using linear regression) and models accounting for possible nonlinear relationships (fitted using Generalized Additive Models).We reasoned that if omitted nonlinear relationships between flowering time and temperature were pervasive in our data and potentially biased our results, then models accounting for nonlinear relationships would tend to perform better than linear regressions among species in our data.We used 10-fold cross-validation to compare the out-of-sample performance (quantified through R 2 values) of linear regressions and GAMs.For each model type (linear regression or GAM), this procedure randomly split the observations for each species into 10 groups, each of which was omitted from a model estimated from the remaining 9 groups.The performance of each of these models was then assessed against the observations omitted in fitting the model, generating 10 out-of-sample R 2 values for each model type (linear or GAM) per species.We then compared the distribution of mean cross-validated R 2 values obtained from linear models and GAMs to assess whether nonlinear models explained additional variance.
Extended Data Fig. 9 | Effects of geographic range on apparent plasticity and adaptation.Relationships between the latitudinal and longitudinal range of specimens of a species and estimates of apparent plasticity (S time ) and apparent adaptation (S space -S time ).Latitudinal ranges in (a, b) and longitudinal ranges in (c, d) were obtained by first removing the extreme 1% of observations among observations for each species.In a-d, blue lines in each panel correspond to bestfit lines obtained using generalized additive models (GAMs), with blue ribbons showing the standard error of the predicted value of the response for each value of the predictors.R 2 are provided as text insets in each panel.Although apparent plasticity and adaptation showed marginally greater magnitude among species with narrower latitudinal and longitudinal range, these relationships explained a very small proportion of the variance.Therefore, we conclude that it is unlikely that differences in latitudinal or longitudinal range size could confound the results presented in the main text.GAMs using apparent plasticity or apparent adaptation as a response and including both latitudinal and longitudinal range as predictors also explained a marginal proportion of the variance (R 2 = 0.10 and R 2 = 0.05, respectively).
Extended Data Fig. 10 | Climatic space captured among specimen collection locations across ecoregions.a-n, Variation in long-term climatic conditions among sites of specimen collection occurring within different Level II ecoregions throughout the contiguous United States.Variation in long-term conditions was calculated using principal components (PCs).PC1 represents a gradient of increasing precipitation seasonality, decreasing temperature seasonality and increasing long-term mean annual temperature.In turn, PC2 represents a gradient of decreasing long-term mean annual precipitation and increasing temperature seasonality (see 'Climatic data').

Fig. 1 |
Fig. 1 | Spatial and temporal relationships between flowering time and temperature resulting from plasticity and adaptation.a,b, Local adaptation acting as the sole driver of flowering time along the gradient (that is, no phenological plasticity) (a) should result in a negligible temporal relationship and a substantial difference between temporal and spatial slopes (b).c,d, In contrast, plasticity acting as the sole driver of flowering-time variation along the gradient (that is, no adaptation) (c) should result in a substantial temporal relationship and negligible differences between spatial and temporal slopes (d).Local adaptation and plasticity jointly influencing flowering time should result in different empirical patterns depending on the direction of their effects.e,f,Plasticity and adaptation operating in the same direction (for example, both negative) (e) should result in a clear temporal relationship and a spatial relationship of substantially greater magnitude (f).g,h, In contrast, plasticity and adaptation operating in opposite directions (for example, plasticity negative,

Fig. 2 |
Fig. 2 | Distributions of and relationship between S space and S time among 1,605North American angiosperms.a, Shaded regions correspond to the kernel density distributions of S time (red) and S space (blue) among species.b, Each point represents a species whose x, y coordinates are given by the MAP estimates for S space and S time , respectively.Colours in b indicate whether sensitivity patterns were consistent with plasticity (green) or adaptation (magenta) as the sole drivers of flowering-time variation along the temperature gradient, with both plasticity and adaptation in a co-or counter-gradient adaptation pattern (blue,

SFig. 3 |
Fig. 3 | Variation in apparent plasticity (S time ) and apparent adaptation (S space − S time ) attributable to differences in phenological niche and native climate among species.a,b, PC1 represents a climate gradient of increasing precipitation seasonality, decreasing temperature seasonality and increasing mean annual temperature.c,d, In turn, PC2 corresponds to a gradient of decreasing mean annual precipitation and increasing temperature seasonality.The colour gradients in each panel represent the predicted magnitude of S time or S space − S time (d °C −1 ) for a combination of mean flowering DOY and PC1 or PC2 values.The predicted surfaces represented by the colour gradients were obtained using three-variable tensor smooths in a GAM framework.In each panel, the value of the third variable (the one not plotted) was fixed at its mean.

Fig. 4 |
Fig. 4 | Variation in apparent plasticity and apparent adaptation among species with varying phenological niches across ecoregions of the United States.a-n, Shaded regions in each panel represent the 95% confidence

Fig. 2 |
Sampling intensity and long-term climatic conditions across collection sites in the continental United States.Pixels correspond to 20 × 20-km grid cells, with their colour representing (a) the total number of specimens collected and (b) their mean PC1 and PC2 values.PC1 represents a gradient of increasing precipitation seasonality, decreasing temperature seasonality and increasing long-term mean annual temperature.In turn, PC2 represents a gradient of decreasing long-term mean annual precipitation and increasing temperature seasonality (see 'Climatic data').