Differences in spatial versus temporal reaction norms for spring and autumn phenological events.

For species to stay temporally tuned to their environment, they use cues such as the accumulation of degree-days. The relationships between the timing of a phenological event in a population and its environmental cue can be described by a population-level reaction norm. Variation in reaction norms along environmental gradients may either intensify the environmental effects on timing (cogradient variation) or attenuate the effects (countergradient variation). To resolve spatial and seasonal variation in species' response, we use a unique dataset of 91 taxa and 178 phenological events observed across a network of 472 monitoring sites, spread across the nations of the former Soviet Union. We show that compared to local rates of advancement of phenological events with the advancement of temperature-related cues (i.e., variation within site over years), spatial variation in reaction norms tend to accentuate responses in spring (cogradient variation) and attenuate them in autumn (countergradient variation). As a result, among-population variation in the timing of events is greater in spring and less in autumn than if all populations followed the same reaction norm regardless of location. Despite such signs of local adaptation, overall phenotypic plasticity was not sufficient for phenological events to keep exact pace with their cues-the earlier the year, the more did the timing of the phenological event lag behind the timing of the cue. Overall, these patterns suggest that differences in the spatial versus temporal reaction norms will affect species' response to climate change in opposite ways in spring and autumn.

To stay tuned to their environment, species need to respond to both short- and long-term variation in climatic conditions. In temperate regions, favorable abiotic conditions, key resources, and major enemies may all occur early in a warm year, whereas they may occur late in a cold year. Coinciding with such factors may thus come with pronounced effects on individual fitness and population-level performance (1–4). As phenological traits also show substantial variability within and among populations, they can be subject to selection in nature (5–7), potentially resulting in patterns of local adaptation (8–10).

At present, the rapid rate of global change is causing shifts in species phenology across the globe (11–13). Of acute interest is the extent to which different events are shifting in unison or not, sometimes creating seasonal mismatches and functionally disruptive asynchrony (3, 14–16). If much of the temporal and spatial variation in seasonal timing is a product of phenotypic plasticity, then changes can be instant, and sustained synchrony among interaction partners will depend on the extent to which different species react similarly to short-term variation in climatic conditions. If geographic variation in phenology reflects local adaptive evolutionary differentiation, then, in the short term, as climate changes, phenological interactions may be disrupted due to the lag as adaptation tries to catch up (17–19). By assuming that space can substitute time, it is possible to make inference about the role that adaptation to climate may play. How well species stay in synchrony will then depend on the extent to which local selective forces act similarly or differently on different species and events.

Local adaptation in phenology may take two forms. 1) The magnitude of phenological change might vary along environmental gradients in ways that intensify the environmental effects on phenological traits, a process known as cogradient variation (Fig. 1). In such a case, the covariance between the genetic influences on phenological traits and the environmental influences is positive. Under this scenario, the effect of environmental variation over space and time will be larger than if all populations were to follow the same reaction norm regardless of location. 2) Genotypes might counteract environmental effects, thereby diminishing the change in mean trait expression across the environmental gradient. In such a case, the effect of environmental variation over space and time will be smaller than if all populations were to follow the same reaction norm regardless of location. This latter scenario, termed countergradient variation, occurs when genetic and environmental influences on phenotypic traits oppose one another (Fig. 1) (20, 21).

Fig. 1.

Schematic illustration showing slopes of phenology on temperature. Adapted with permission from ref. 30. A corresponds to phenological plasticity with respect to temperature and no local adaptation. B reveals phenological plasticity with respect to temperature plus cogradient local adaptation. C reveals phenological plasticity with respect to temperature plus countergradient local adaptation. For each scenario, we have included two examples of events showing this type of pattern in our data. For the exact climatic cues related to these biotic events, see . In each plot, the red lines correspond to the within-population reaction norms through time (i.e., temporal slopes within locations), and the blue line corresponds to the between-population reaction norm (i.e., spatial slopes). If all populations respond alike, then the same reaction norm will apply across all locations, and individuals will respond in the same way to the cue no matter where they were, and no matter whether we examine responses within or between locations. If this was the case, then the reaction norm would be the same within (red lines) and between locations, and the blue and the red slopes would be parallel (i.e., their slopes identical). This scenario is depicted in A. What we use as our estimate of local adaptation is the difference between the two, i.e., whether the slope of reaction norms within populations differs from that across populations. If the temporal slopes are estimated at a relatively short time scale (as compared to the generation length of the focal organisms), then we can assume that within-location variation in the timing of the event reflects phenotypic responses alone, not evolutionary change over time. This component is then, per definition, due to phenotypic plasticity as such, i.e., to how individuals of a constant genetic makeup respond to annual variation in their environment. By comparison, the spatial slope (i.e., the blue line) is a sum of two parts: first, it reflects the mean of how individuals of a constant genetic makeup respond to annual variation in their environment, i.e., the temporal reaction norm defined above. These means are shown by the red dots in A–C. However, second, if populations differentiate across sites, then we will see variation in their response to long-term conditions, with an added element in the spatial slope reflecting mean plasticity plus local adaptation. Therefore, if the spatial slope differs from the temporal slope, this reveals local adaptation (see for further details). Such local adaptation in phenological response may take two forms. 1) The magnitude of phenological change might vary along environmental gradients in ways that intensify the environmental effects on phenological traits, a process known as cogradient variation (Fig. 1). In such a case, the covariance between the genetic influences on phenological traits and the environmental influences is positive. Under this scenario, variation in the environmental cue over space and time will cause larger variation in phenological timing than if all populations were to follow the same reaction norm regardless of location. 2) Genotypes might counteract environmental effects, thereby diminishing the change in mean trait expression across the environmental gradient. In such a case, the effect of variation in the environmental cue over space and time will be smaller than if all populations were to follow the same reaction norm regardless of location. This latter scenario, termed countergradient variation, occurs when genetic and environmental influences on phenotypic traits oppose one another (C).

For phenology, the overall prevalence of co- versus countergradient patterns is crucial, as it will dictate the extent to which local adaptation will either accentuate or attenuate phenological responses to temporal shifts in climate (10). Across environmental gradients in space, the relative prevalence of counter- versus cogradient variation in spring versus autumn will critically modify how climatic variation affects the length of the activity period of the entire ecological community. Overall, geographic variation in the activity period will be maximized when events in autumn and spring differ in terms of whether they adhere to patterns of co- or countergradient variation.

Although the study of individual species and local species communities has revealed fine-tuning of species to local conditions (22), and a wealth of studies report shifts in phenology worldwide (23), we still lack a general understanding of how the two tie together: how strong is local adaptation in the timing of events, and how do they vary across the season? Here, a major hurdle to progress has been a skew in the focus of past studies: our current understanding of climatic effects on phenology has been colored by springtime events (24–26), whereas events with a mean occurrence later in the season have been disproportionately neglected (27). To achieve satisfactory insight into how climate and its change affect the timing of biological activity across the season, we should thus ask how strongly phenology is influenced by climatic variation, what part of this response reflects phenotypic plasticity and what part evolutionary differentiation, and how the relative imprint of the two varies across the season. Addressing these pertinent questions is logistically challenging (e.g., ref. 28). Therefore, few studies have tackled them outside of the laboratory (29).

Phillimore and coworkers (10, 30) proposed an elegant technique for identifying the relative roles of plasticity and local adaptation in generating spatiotemporal patterns of phenological variation. The rationale is to use a space versus time comparison (10, 30) (but see ref. 31 for criticism), drawing on the realization that at any one site, local conditions will vary between years. To be active at the right time, species will thus need to respond to temporal variation in climatic conditions. Let us assume that a focal species times some aspect of its annual activity (a species-specific “phenological event”) by reacting to a single environmental cue (e.g., the crossing of a given temperature sum). Now, if there were no differentiation between populations and all populations followed the same reaction norm, then with variation in the relative timing of the cue over time, all populations would react in the same way to the same cue regardless of spatial location (Fig. 1). At the level of population means across space (blue line in Fig. 1), we would then see a relationship between phenological event and cue timing identical to year-to-year variation within locations (red lines in Fig. 1). However, if populations differentiate across sites, then we will see an added component in the spatial slope, reflecting the contribution of local adaptation to the mean phenology of the populations. By subtracting the within-population temporal slope from the spatial slope, we will thus achieve a direct measure of local adaptation (10), henceforth called Δb (30).

Importantly, the temporal slope (i.e., the local phenological response to local year-to-year variation in the cue) can be either steeper or more shallow than the spatial slope (Fig. 1 vs. Fig. 1)—the former being a sign of countergradient local adaptation, the latter of cogradient local adaptation (20, 21, 32). For a worked-through example of how this methodology is applied to the current data, see .

Here, we adopt temperature sums as widely used predictors of phenological events (33–35) and treat the difference between the spatial and temporal slopes of phenological events on such sums as our estimates of local adaptation in reaction norms (). Pinpointing the relative roles of plasticity and microevolution from spatiotemporal observations in the absence of direct measures of fitness will, per necessity, rely on several assumptions (for a full discussion, see ref. 36). However, given the adequate precaution, such quantification allows a tractable way toward estimating local adaption on a large scale (8, 10, 30, 36–38).

A key requirement for the successful application of this approach to resolving patterns across events of different relative timing is the existence of abundant data covering a large geographic area (30, 36). The extensive phenological data-collection scheme implemented at hundreds of nature reserves and other monitoring sites within the area of the former Soviet Union offers unique opportunities for addressing community-level phenology across a large space and long time (39). From this comprehensive dataset spanning 472 monitoring sites, 510,165 events and a time series of up to 118 y (Fig. 2 and ref. 39), we selected those 178 phenological events for which we have at least 100 data points that represent at least 10 locations (). These events concerned 91 distinct taxa ().

Fig. 2.

Study sites and spatiotemporal patterns in climatic and phenological data. A shows the depth of the data and the spatial distribution of monitoring sites, with the size of the symbol proportional to the number of events scored locally. Since the selection of sites differed between events (39), in A, we have pooled sites located within 300 km from each other for illustration purposes. B shows the mean timing (day of year) of a phenological event: the onset of blooming in dandelion (Taraxacum officinale). C shows the mean timing (day of year) of a climatic event: the day of the year when the temperature sum providing the highest temporal slope for the onset of blooming in dandelion was first exceeded, computed as the mean over the years considered in B. For a worked-through example estimating reaction norms and metrics of local adaptation (Δb) for this species, see .

To express data on species phenology and abiotic conditions in the same currency, we related the dates of the phenological events (e.g., the first observation of an animal, or first flowering time of a plant species; ) to the dates when a given thermal sum (34, 35) was first exceeded. This choice of units has a convenient consequence in terms of the interpretation of slope values: if the date of phenology changes follows one-to-one the date of attaining a given temperature sum, then the slope will be one—an assumption frequently made but rarely tested in studies based on growth-degree days. The observed reaction norms can then be compared to this value. A value below 1 will signal undercompensation, i.e., that the earlier the cue, the larger the relative delay of the phenological event compared to its cue. By contrast, a value larger than 1 would signal overcompensation, i.e., that with an advancement of the cue, the timing of the phenological event will be advanced even more.

Since thermal sums can be formed using a variety of thresholds, we used a generic approach and considered dates for exceeding a wide range of both heating and chilling degree-day sums (34, 35) (see Material and Methods for more information). As there is also evidence that sensitivity to temperature arises after a certain time point (13, 36), we calculated each heating and chilling degree-days sum for a range of starting dates. For each of the resulting 2,926 events, we then picked the variable that offered the highest temporal slope estimate, i.e., the largest within-location change in the timing of the event with a change in the timing of the cue (see Material and Methods for more information). Following the rationale outline above, this will be the most appropriate optimization criterion, since it selects the cue to which the phenological event responds the strongest to over time.

Results

Overall, most metrics of heating and chilling degree-day sums proved highly correlated with each other, with strong negative correlations dominating among heating and chilling degree-day dates (). This pattern prevailed both in space () and in time (). Hence, warmer sites come with both early springs and late autumns (), and within years, if the local spring comes early, then autumn comes late—so summers tend to be either short or long ().

Due to their strong correlations (), multiple metrics offered similar slope estimates (), and our results are thus robust to the exact selection of climatic event. Overall, species proved finely attuned to the local thermal environment, with spatial slopes of phenology on temperature sums typically being close to (but still below) 1 (Fig. 3). In other words, a 1-d shift in attaining a given temperature sum was typically reflected in a 1-d shift in the commencement of (or end of) species’ activity, and species’ phenology was thus almost fully matched with temperature phenology from the geographic point of view. Importantly, early phenological events typically shifted earlier and late events shifted later at warmer sites (Fig. 3 and ). Within sites, phenological plasticity was too small to result in an equivalent match, with most temporal slopes being substantially shallower than 1 (Fig. 3 and ). Thus, the overall match in species phenology to the local thermal environment was not only shaped by phenotypic plasticity, but also bore a significant imprint of local adaptation (Fig. 3 and ). As a result, for the majority of species, the 95% credibility limits of Δb (the difference between the spatial and temporal slopes) did not overlap with 0 (Fig. 3 and ). This pattern varied little between different types of events (Fig. 3).

Fig. 3.

The relationship between the mean timing of an event (day of year) and the slope of phenology on dates of achieving specific degree-day sums. Shown are spatial (A) and temporal (B) slopes (i.e., temporal slopes within populations), with C showing the difference between them, Δb, as an estimate of local adaptation (see main text and Fig. 1 for details). Phenological events are shown by filled circles, with the trophic level in question identified by color: primary producers are shown in green, primary consumers in yellow, secondary consumers in black, and saprotrophs in orange. A quadrat around the circle identifies species for which the 95% HPD does not overlap with 0. For visual comparison, a black line has been added to A–C at a slope value of 0 (indicating no relationship), and a red line has been added at a slope value of 1 (indicating a perfect relationship, i.e., a shift of 1 d in the timing of the event with a shift of 1 d in the date of achieving the degree-day sum in question). Dashed curves refer to model estimates provided in . For the degree sum related to individual events, see .

In terms of the mean timing of individual phenological events within years, events occurring particularly early or late within the year proved insensitive to annual variation in any temperature sum (Fig. 3 and ). This lack of response occurred even though at a larger spatial scale the timing of these events varied widely with spatial variation in average climate (Fig. 3 and ) and despite substantial variation in the exact timing of these early and late events across years—with no apparent association between the timing of an event and its absolute variability among years (). Thus, early and late events show little plastic response to the heating and chilling degree-day sums. In contrast, events occurring closer to midseason were strongly affected by annual local temperature conditions (Fig. 3 and ), showing a strong plastic response but little imprint of local adaptation (Fig. 3 and ).

Our most important finding was a seasonal pattern in the local adaptation of responses to the local climate, i.e., in the local reactions norms. Local adaptation proved strongest for events occurring early and late in the season but weak for events occurring in midseason (Fig. 3 and ): the sign of the slope (Δb) crosses the x axis around midsummer, with early events dominated by cogradient adaptation and late events by countergradient adaptation. In other words, local adaptation tends to accentuate the effect of climatic variation in spring and attenuate it in autumn.

All of the above patterns prevailed in the face of increasingly stringent criteria, including the filtering of data by successively higher temporal slope values (to focus on cues of increasingly strong effects) and splitting the data into two subsets along latitude 57°N (≤57°N vs. >57°N). Within the southern half of the data, the slope and R2 were clearly lower than in the northern half () or the overall data. Yet, this added scatter may be related to substantially more complex environmental gradients in the southern (Fig. 1) than in the northern part. For detailed results of these added analyses, see .

Discussion

The salient insights emerging from our results are a consistent pattern of widespread responses to temperature and local adaptation in these responses—with the type of differentiation varying directionally from spring to autumn along a co- to countergradient continuum.

The current observation of apparent local adaptation is consistent with the notion that phenological traits may be strongly adaptive (5–7). As an example of how adaptive local evolution may shape spatial and temporal patterns in the mean flight dates of multiple species of butterflies, Roy et al. (8) found that all species examined showed a plastic response to temperature but that emergence dates were still synchronized among populations—suggesting pervasive countergradient local adaptation. As a demonstration of convergent evolution of local adaptation, Thuiller et al. (40) showed that species growing in similar regions had developed similar phenologies. In line with these results, we found that for 94 of 178 phenological events examined, the 95% credible interval of Δb did not overlap with 0, supporting a clear imprint of local adaptation on phenological reaction norms (Fig. 3). As we clearly lack any direct evidence on whether and to what extent the patterns are actually adaptive in nature (i.e., associated with fitness benefits), we stress that our results point to an explicit hypothesis amenable to testing by e.g., extensive translocation experiments: that over time local organisms should be better attuned to the local temperature regime than individuals from elsewhere, as reflected in higher fitness.

However, not all events seem to be equally differentiated in space. By including events spanning the whole season (27), we resolved strong differences in the added component of local adaptation between different parts of the year. This pattern appeared despite the lack of any detectable association between the timing of an event and its absolute variability among years (). Events occurring particularly early or late within the year showed little within-site response to annual variation in local temperature sums (Fig. 3) but responded distinctly differently to variation in average climatic conditions occurring at a larger spatial scale (Fig. 3). The patterns observed were not compatible with any separate impact of, e.g., photoperiod (), and what exact proximate cue they are triggered by is then an open question. Events occurring in mid-spring and autumn responded sensitively to heating degree-day sums and showed high levels of local adaptation. Most intriguingly, the pattern of variation was different in spring and autumn, with springtime events being dominated by cogradient variation and autumn events by countergradient variation.

Covariation between genetic and environmental influences across ecological gradients has been generally predicted (32) and sometimes observed (20, 21). Such patterns are essential for improving our understanding of the mechanisms that trigger evolutionary responses and shape the dynamics of populations (20, 21). Overall, environmental variation in time and space can be expected to generate genetic compensation, especially if these fluctuations are asymmetrical (i.e., when the environment varies differently in both space and time) (32). If the environment fluctuates equally in time and space, then phenotypic plasticity may be favored at the expense of local adaptation (32).

That the current patterns emerged so clearly may be attributed to the extent and depth of our data (39). Overall, our data stretched over an unusually large geographic area (Fig. 2 and ref. 39) and were acquired and analyzed in a consistent way (39) (in contrast to metaanalyses). These aspects will increase the signal to noise ratio. What may, on the other hand, reduce this ratio is the fact that our dataset is based on first dates instead of mean population event dates. Previous studies (41) have suggested that first phenological dates may be affected by changes in population density, observation effort, and observability. However, Phillimore et al. (10) offer explicit simulations showing that population density or recorder effort will only pose a bias if there is a latitudinal or temporal trend in them. Our dataset was systematically collected with sampling effort remaining nearly constant across sites—in strong contrast to studies based on volunteer observations. By including so many systematically collected events from so many sites, and by explicitly including the effect of timing, we may thus have revealed imprints perhaps hidden (yet extant) in previous data on events from shorter or inconsistent time periods.

In terms of the absolute values of the slope estimates, our choice of a coherent unit (dates) for both abiotic predictors and species’ response allows an intuitive interpretation of slopes. Here, we were able to show that phenological events oftentimes tracked temperature cues with high accuracy, including some shifting by exactly 1 d with a shift of 1 d in the date when a given temperature sum was achieved (Fig. 3). However, this match was based on neither the temporal reaction norm (Fig. 3) nor local adaptation per se (Fig. 3); either mechanism in isolation allowed species to slide somewhat earlier or later but not to an extent fully matching the realized variation (note that in Fig. 3 , average values are well below 1).

Clearly, the current approach is associated with limitations, too (10, 15, 30). Estimating the relative contribution of plasticity versus local adaptation from spatiotemporal data are not without complications (see ref. 14 for a full discussion), and the validity of this method has sometimes been challenged (15). Of the concerns raised in ref. 36, the fourth needs particular attention in the current case. First, the approach (10, 30) relies on the correct identification of phenological cue(s). Second, the reaction norms for phenological data to temperature sums might not always be strictly linear; third, phenotypic plasticity might not necessarily be constant among sites; and fourth, the time span considered by us may, in theory, be long enough to allow some of the within-site variation observed to be caused by evolutionary change through time (i.e., microevolution) (2). In terms of the first concern, we note that the current inference seems reasonably robust with respect to this assumption (), in terms of the second that the relations examined by us seem reasonably linear and homoscedastic (). In terms of the third item, we find that previous studies report little evidence for intraspecific geographic variation in mean population plasticity (8, 10, 30, 36, 42). When it comes to the fourth item, we note that most events emanated from the time period after the 1960s (39) and that change during this time would call for a change more rapid than compatible with our other findings.

Importantly, the current match between local climate and species phenology in space and time is the combined product of evolutionary differentiation and local phenotypic plasticity. The patterns uncovered suggest that in isolation, phenotypic plasticity as such is too weak to provide an accurate match between species’ phenology and variation in temperature (43): the slope of the timing of the event on the timing of the cue was consistently below 1 (Fig. 3). Today’s patterns of local adaptation will partly reflect historical selection pressures and processes, and these forces have apparently been strong enough to create ubiquitous differentiation across a wide set of species. However, the exact strength of this patterning varies substantially across the season, with local plasticity dominating phenological adjustments around midsummer and geographic differentiation governing species responses early and late in the year. In spring, adaptive differentiation added to geographic variation in phenology, whereas in autumn, it reduced it. How quickly current patterns will equilibrate with current changes in climate we do not know. In the best of all possible worlds, rapid, continuing evolutionary change will allow the efficient tracking of ambient conditions (44)—but the mechanisms for such tracking seem to vary substantially for events occurring during different parts of the season. Only time and intensive studies will tell how they play out next.

Materials and Methods

Phenological Records.

Our study is based on a dataset collected across 472 localities in Russia, Ukraine, Belarus, Latvia, Lithuania, and Estonia (Fig. 2) during a 115-y period (from 1899 to 2014). During this period, researchers recorded the first and/or last occurrence of a predefined list of phenological and weather-related events, with somewhat different time periods covered at different sites. The full dataset is described in ref. 39.

At some of the monitoring sites, researchers conducted observations at multiple locations. To perform comparable analysis among regions, we used only one locality per site. When more were available, we picked one at random. Of the data thus generated, we focused on the set for which we have at least 100 data points that represent at least 10 monitoring sites. The selection of monitoring sites and years varied among the events, since for each event, we used all available data. These criteria resulted in a final dataset of 178 phenological events from a total of 92 taxa. For the majority of events, most data emanated from the time period after the 1960s (39).

Climatic Descriptors.

As predictors of phenological events, we adopted widely used temperature sums (34, 35). To express data on species phenology and abiotic conditions in the same currency, we related the dates of the phenological events (e.g., the first observation of an animal or first flowering time of a plant species; ) to the dates when a given thermal sum (34, 35) was first exceeded. This choice of units has a convenient consequence in terms of the interpretation of slope values: if the date of activity changes is directly proportional to the date of attaining a given temperature sum, then the slope will be 1, and slopes of both plastic responses and local adaptation can then be compared to this value.

Since thermal sums can be formed using a variety of thresholds, we used a generic approach and considered dates for exceeding a wide range of both heating and chilling degree-day sums (34, 35). We define the heating degree-day sum at day as

where is the temperature of day The temperature threshold and the initial day are parameters of this model. Following ref. 35, we considered for the following range of values:

We define the heating-degree date as

If never exceeds we set

We applied different threshold values for each phenological event. To determine suitable threshold values we first computed the maximal value (i.e., the value that is obtained at of for each and for each year–site combination. For each we defined as the 50% quantile of the maximal values over the years and sites. defines a threshold value that is large in the sense that in only half of the cases, it is reached at all over the entire year. For each we applied the threshold values and where For 10 values for and 11 values for the threshold, this produces 110 different heating-degree dates (many of which can be expected to be highly correlated with each other) for each

We define the chilling degree-day sum at day as

Following ref. 34, we considered a range of values for :

We define the chilling-degree date as

If never exceeds we set

To determine suitable threshold values we first computed the maximal value (i.e., the value is obtained at of for each and for each year–site combination. For each and we defined as the 50% quantile of the maximal values over the years and sites. defines a threshold value that is large in the sense that in only half of the cases, it is reached at all over the entire year. Like with the heating sums, for each and we applied the threshold values and where With 4 values for and 11 values for the threshold, this produces 44 different chilling-degree dates (many of which can be expected to be highly correlated with each other; ) for each

With 110 heating-degree dates and 44 chilling-degree dates, we generate, in total, 154 temperature-related predictors for each all having the unit of date. In addition to the numeric values, the values may be missing in case of missing data or infinite in case of the threshold being never exceeded. Note that while the temperature predictors are calculated separately for each phenological event (due to variation of where and when they are sampled and how data from nearby sites are merged), we used the same numerical values for the maximal thresholds and for all phenological events. To do so, we determined these thresholds from pooled data over all phenological events.

Finally, as there is evidence that sensitivity to temperature arises after a certain time point (13, 36), we calculated each heating and chilling degree-day sum for a range of starting dates For heating-degree dates, we consider the range of =1,11,21,31,...181 (i.e., approximately the first 6 mo, in increments of 10 d), and for chilling degree-days, we considered dates from the beginning of June to the end of November (again 6 mo, starting from day 172 and adding increments of 10 d). Combined with the 154 heating and chilling-degree dates described above, this results in a total of 154 × 19 = 2,926 predictors considered.

Deriving Spatial and Temporal Slope Estimates.

We fitted models of 2,926 climatic variables to each of 178 events. Here, the proportion of variation in species phenology attributable to variation in degree-days, and the relevant slopes of species activity on degree-days, were quantified by fitting a bivariate analysis of variance using the R package MCMCglmm (10, 45). We applied the default settings of 13,000 iterations with burn-in of 3,000 and thinning of 10. This chain length was selected on the basis of pilot experiments, which indicated that the results were consistent with 10 and 100 times longer chains. We then selected the default level due to its computational feasibility. We assumed for both the R and the G matrices an inverse-Wishart prior with expected covariance set to the identity matrix and degree of belief parameter set to 0.002.

We estimated the posterior means as 95% highest posterior density (HPD) intervals for the spatial and temporal slopes following the approach of ref. 10. Thus, the spatial slopes were derived from the posterior distributions of the 2 × 2 variance–covariance matrices corresponding to the level of location, while also allowing for covariance at the residual level. The temporal slopes were then derived from the posterior distributions of the 2 × 2 variance–covariance matrices corresponding to the level of year.

Among candidate models, the ones including dates occurring after the phenological event occurred will seem unlikely (since a phenological event cannot be causally related to a climatic event occurring after it; ). Yet, since not every single combination of H0, C0 and t0 was sampled, there is some scope for estimation error. To avoid excluding the best predictor by applying overly stringent criteria, we thus added a small tolerance of 5 d for predictions after the event. Here, the number 5 is identified as a compromise reflecting the limited resolution of considering 110 different heating-degree dates and 44 different chilling-degree dates (and 19 values of t0), thereby accounting for a scenario where the metric identified as optimal were off by a few (but not very many) days.

Further, we also accounted for complications generated by the many metrics considered while varying t0. In essence, some combination of t0 and the other choices (see indices above) yielded no or very little variation in the climatic predictor among locations and years (e.g., day 1 for all cases, due to nonmeaningful combination of choices). Fitting the models with these predictors gave spurious results. Thus, we narrowed the set considered to metrics that offered not just a high mean estimate but also a reasonable confidence around it. In other words, we omitted metrics for which the length of the 95% credibility interval of Δb was more than 1 d. For five events, the interval length of 1 d was exceeded for all metrics, and the event thus exclude from final analyses.

Of the remaining set of candidate variables, for each event we picked the variable which offered the highest temporal slope estimate. We did so because this quantity, the temporal slope, reflects the covariance between the climatic variable in question and local phenology, and thus the sensitivity of realized phenology to climatic variability. A justified metric of the climate should thus be one to which species show a sensitive response.

Establishing Seasonal Patterns in Phenotypic Plasticity and Local Adaptation.

To establish the impact of timing on the temporal slope, the spatial slope and the local adaptation (Δb), we fitted univariate general linear models of the respective slopes as functions of the mean timing (day of year [DoY]) and its quadratic term—the latter to allow for nonlinearities in the pattern. Here, each data point was the Δb observed for an individual event (n = 178). In each case, we assumed normally distributed errors and used an explicit, likelihood-based rationale for model selection: starting from a constant-only model, we used type1 likelihood-ratio analysis to test for improved model fit when adding an independent variable. We first tested for the presence and shape of a relationship between the respective slope value and the timing of an event, by sequentially adding DoY and DoY2, retaining either the first or both if providing a significant improvement in deviance, tested against a χ2 distribution with the appropriate degrees of freedom. These models were fitted in R version 3.6.2 (46), function lm. To check for biases due to species- or event-specific impacts, we validated our results by running separate analyses using a single event per taxon, with the event selected to be the event with the highest number of independent records for each single taxon ().