Estimating the ability of plants to plastically track temperature‐mediated shifts in the spring phenological optimum

Abstract One consequence of rising spring temperatures is that the optimum timing of key life‐history events may advance. Where this is the case, a population's fate may depend on the degree to which it is able to track a change in the optimum timing either via plasticity or via adaptation. Estimating the effect that temperature change will have on optimum timing using standard approaches is logistically challenging, with the result that very few estimates of this important parameter exist. Here we adopt an alternative statistical method that substitutes space for time to estimate the temperature sensitivity of the optimum timing of 22 plant species based on >200 000 spatiotemporal phenological observations from across the United Kingdom. We find that first leafing and flowering dates are sensitive to forcing (spring) temperatures, with optimum timing advancing by an average of 3 days °C−1 and plastic responses to forcing between −3 and −8 days °C−1. Chilling (autumn/winter) temperatures and photoperiod tend to be important cues for species with early and late phenology, respectively. For most species, we find that plasticity is adaptive, and for seven species, plasticity is sufficient to track geographic variation in the optimum phenology. For four species, we find that plasticity is significantly steeper than the optimum slope that we estimate between forcing temperature and phenology, and we examine possible explanations for this countergradient pattern, including local adaptation. &NA; We use spatiotemporal spring phenology observations for 22 UK plant species to estimate the temperature‐mediated plasticity of each species and the degree to which optimum timing changes with temperature. We find that all species are highly plastic and that in most cases, this plasticity is adaptive (i.e. it partially tracks temperature‐mediated changes in the phenological optimum). Figure. No caption available.


Introduction
Shifts in phenology are among the most widely reported ecological responses to changing climate across different ecosystems and taxa (Walther et al., 2002;Parmesan & Yohe, 2003;Settele et al., 2014). For temperate plants, the timing of spring events, such as leafing and flowering, has been especially well recorded by both professional and citizen scientists. Analysis of the resultant longitudinal data sets reveals that spring phenology has advanced in many species over the past few decades, coincident with rising temperatures (Fitter & Fitter, 2002;Amano et al., 2010;Thackeray et al., 2016). Some of the advancement in phenology will be due to microevolutionary change (Franks et al., 2014), but plastic responses to temperature probably dominate (Nicotra et al., 2010).
Temperate plants often exhibit pronounced temperature-mediated plasticity in their spring phenology, as documented via longitudinal studies of individuals (Vitasse et al., 2010;www.trackatree.org.uk), geographic transplants (Kramer, 1995) and experimental approaches (Franks et al., 2014). Such plastic phenological responses are often assumed to be adaptive, allowing plants to track spatial and temporal variation in the onset of benign environmental conditions. Consequently, observed advances in spring phenology suggest that optimal phenologies have advanced as the climate has warmed. However, little information is available as to whether observed phenologies are advancing at the same rate as optimal phenologies, and what the demographic consequences of any shortfall are (see Wilczek et al., 2014 for an exception). Measuring how optimal phenologies change per unit change in temperature, termed the environmental sensitivity of selection (B) (Chevin et al., 2010), is logistically challenging. It requires that the phenology-fitness surface be estimated in several different environments, followed by characterization of the temperature dependence of the timing of the fitness peak . As a result, outside of the laboratory there have been very few estimates of this key parameter (Michel et al., 2014).
Usually B is used to map a temporal change in temperature to a temporal change in the optimum. However, temperatures also vary in space and we can map spatial changes in the temperature cue to spatial changes in the optimal phenology. In a simple model where temperature varies clinally (with latitude for example), and phenotypic optima depend linearly on temperature, spatial differences in the observed phenology quickly equilibrate to spatial differences in the optimal phenology if there is no spatial variation in population density and dispersal is symmetric (Felsenstein, 1977;Slatkin, 1978). Under these conditions, the spatial relationship between phenology and the driving temperature variable could be used as an estimator of B. However, when spatial changes in temperature also have a stochastic component, then this estimator is biased towards zero. Hadfield (2016) recently demonstrated that an unbiased estimator of B can be obtained from the ratio of the rate at which phenology changes with latitude to the rate at which the temperature cue changes with latitude. In a recent study of the nesting phenology of four passerine birds, Phillimore et al. (2016) extended the approach of Phillimore et al. (2010) for separating the contributions of plasticity and local adaptation from spatiotemporal phenological datato estimate B. An advantage of this spatial approach for estimating B is that there already exists an abundance of amenable data for a range of taxa. However, the theory and statistics required to estimate B from spatial observations in the absence of direct measures of fitness rests on several assumptions (for a full discussion, see Table 1).
Regression-based estimation of B from observational data relies on the correct phenological cue(s) having been identified (Chevin & Lande, 2015). In most temperate plants, the primary drivers of spring phenology are temperature and photoperiod (Rathcke & Lacey, 1985;Polgar & Primack, 2011), although the relative importance of these cues will vary among species (Laube et al., 2014) and is a source of debate (Chuine et al., 2010;K€ orner & Basler, 2010). Temperatures at different times of the year can have opposing effects, with the forcing effect of warm springs usually advancing phenology, while warm conditions in autumn/winter may delay phenology via effects on dormancy induction, breaking dormancy and stimulating growth (Murray et al., 1989;Polgar & Primack, 2011;Laube et al., 2014;Roberts et al., 2015). Experimental studies reveal that photoperiod plays a key role in some species (Caffarra & Donnelly, 2011), although the precise nature of any interactions between photoperiod and the response to forcing and/or chilling temperatures is not known for most species (Polgar & Primack, 2011;Vitasse & Basler, 2013). While longitudinal data from a single site are not informative about the role that photoperiod plays, spatial replication of longitudinal data across latitudes may be (Phillimore et al., 2013), and such data are often collected by citizen science schemes.
The UK Phenology Network (UKPN) citizen science scheme was set up in 1998 and now comprises millions of phenology observations in space and time. In this study, we apply a recently developed statistical framework (Hadfield, 2016;Phillimore et al., 2016) to the spring phenology (first leafing and first flowering) of 22 plant species that have been monitored by the UKPN (Table S1). Our main aims are to estimate (i) the temperature sensitivity of the phenological optimum (B) and (ii) the degree to which phenological plasticity is adaptive and tracks geographic variation in the optimal phenotype. As warmer spring conditions may herald an advance in favourable conditions for growth, we predict that B with respect to forcing will be negative. Theory predicts that the plastic slope will be in the same direction as B but will be shallower unless the environment of development and selection are the same (Tufto, 2015). Our analysis relies on us having identified the correct phenological cues, and a secondary aim of our study was to identify the sensitivity of species to forcing, chilling and photoperiod.

Spatiotemporal data
We used phenological records of first flowering and leafing collected by citizen scientists from the UK Phenology Network (UKPN; www.naturescalendar.org.uk) over the period 1998-2014 (see Table S1 for details). We selected species with >3500 records and excluded taxa for which there were known issues associated with data collection. This meant that we excluded species with common cultivars (e.g. snowdrop and primrose) or those that may be easily confused with other species. We selected first leafing events in preference to first budburst because we have found the first leaf phenophase to be more straightforward to observe, and preliminary analyses revealed that this measure was less subject to among recorder variance. We focused on spring events because the role of temperature as a cue is better understood for spring than autumn events (Gallinat et al., 2015). Prior to analysis, we visually inspected histograms and removed extreme outliers for each species. In order to minimize measurement errors introduced by novice recorders (Dickinson et al., 2010), for each species we excluded all data collected by participants who only contributed records for a single year. The number of filtered observations ranged between 2805 for sessile oak and 22 177 for lesser celandine (Table S1).
Each phenological observation was assigned to the 5 9 5 km grid cell (hereafter 5-km grid cell) in which it was reported, and matched to daily air temperature data interpolated between recording stations on the same grid for the appropriate year (Perry & Hollis, 2005;Perry et al., 2009); www.metoffice.gov.uk/climatechange/science/monitoring/ While the temperature sensitivity of plant phenology is often modelled using growing degree-day models (e.g. Chuine, 2000), here we adopt a reaction norm approach due to its amenability to linear statistical modelling and to facilitate comparisons with theoretical models of quantitative trait evolution (e.g. Chevin et al., 2010;Hadfield, 2016). Where growing degree-day and linear reaction norm approaches have been applied to the same data sets, they have been shown to provide similar insights into phenological cues and responses (Phillimore et al., 2013;Roberts et al., 2015). We use a sliding-window approach to identify the window(s) during which mean temperature best predicts phenology. It is possible that the window of thermal sensitivity varies geographically and our phototemp model allows for a latitudinal cline in the window of forcing temperature sensitivity. ii. The temperature cue(s) that determines plasticity (b) also determines the optimum (B).
Where the environment of selection (i.e. that which determines the optimum) and development (the cue) are the same, our estimate of B will correspond to the optimal slope of phenology on temperature. If the correlation between the two environments is <1, our estimate of B will correspond to the optimal phenological response to the environment of development, which is shallower than the optimal phenological response to the environment of selection (Tufto, 2015). To test whether the environment of selection is the same as the environment of development would require information on the fitness of individuals in different environments (e.g. . iii. The selected temperature variable(s) is/are the sole determinant of the optimum.
While temperature may have a direct effect on the optimum for a species, it is quite likely that some of its effect is indirect, via the phenology of interacting species (e.g. forest tree and understorey species competing for light in spring, or flowers competing to attract pollinators). If the identity of interacting species varies clinally, then this may cause B lat or B lon (see methods) to overestimate or underestimate B. Similarly, other environmental variables that vary geographically and affect phenology, such as precipitation, may cause our estimates of B lat or B lon to over-or underestimate B. Where chilling plays a role in addition to forcing, our ability to separate the effects of chilling and forcing on B depends on how correlated the two temperature variables are. In practice, we find that chilling and forcing temperatures are highly correlated in space (Table S4). This means that our estimate of B based only on forcing temperatures will be biased in the direction of any relationship between chilling temperatures and optimum timing that exists. iv. Population density is constant in space.
Violation of this assumption is anticipated to lead to underestimation of B (Garc ıa- Ramos & Kirkpatrick, 1997). Atlas data reveal little present-day geographic heterogeneity across Britain in the abundance of larch, rowan, silver birch, field maple (although this species is absent from the north of Britain) and alder. However, horse chestnut, beech, pedunculate oak and ash all appear about twice as frequently in plots in the south of Britain than they do further north, whereas sycamore has elevated abundance at mid-latitudes (San-Miguel-Ayanz et al., 2016). If spatial heterogeneity in abundance leads to a severe underestimation of B across species, we would expect to find a negative correlation between our estimates of |B| and the absolute change in abundance with latitude or longitude, which we do not (appendix S2). v. Migration is symmetric among populations.
At the range limits, migration will be from a single direction and migration load is expected to perturb such populations from the optimum (Hadfield, 2016). To assess whether this impacts on our estimation of B, we plot the residuals and 150-km grid cell best linear unbiased predictors (BLUPs) as a function of latitude and longitude, and visually inspected whether there was a tendency for values to depart from 0 at the latitudinal and longitudinal extremes. We observed such deviations in BLUPs at one or both latitudinal extremes for wood anemone, lesser celandine, sycamore, hazel and rowan, and for these species, B lat may be biased towards 0. Most species showed such departures in BLUPs over longitudes, implying that B lon will tend to be biased towards 0. vi. Populations are at migration-selection equilibrium.
Violation of this assumption would cause B lat and B lon to be biased away from B towards b. Introduced species, such as horse chestnut, larch and sycamore, will violate this assumption. We do not know whether the other remaining species obey this assumption, ukcp09/). To calculate the day length (time from sunrise to sunset) in minutes for each day, sunrise and sunset equations (Meeus, 1991) were applied to the centroid of each 5-km grid cell. Each 5-km grid was assigned to a 25-km grid cell and 150-km grid cell, with the latter treated as an arbitrary definition of a population in analyses, as in earlier studies using a similar approach (Phillimore et al., 2010(Phillimore et al., , 2016.

Statistical analyses
We fitted a series of linear mixed models in ASReml-R (Butler et al., 2009;Gilmore et al., 2009) that were designed to identify the environmental cues that best explain the spatiotemporal variation in phenology of each species. All models (parameters summarized in Table 2) treat the ordinal date of the phenological observation as a response variable and include 150-km grid cell, year, 25 km:year and a residual term as random effects. Our motivation for including the 25 km:year term was to account for pseudoreplication of interpolated temperatures within a 5-km grid cell and year.
Our null model included only the intercept as a fixed effect. We also considered geographic and temporal cline models in order to (i) identify broad spatial and temporal trends and (ii) act as an additional baseline against which the performance of cue-based models can be compared. Our simple clinal model included year (as a continuous variable), latitude and longitude as fixed effects (geo1). A more complex clinal model also included the interaction between latitude and longitude, as well as quadratic terms of latitude and longitude (geo2).
All subsequent models include environmental cues (Fig. 1). The first was consistent with a photoperiod threshold triggering phenology (photo). The ordinal date at which the specified minutes of daylight (we considered values between 486 and 980 minutes at intervals of 4 minutes) was first reached in each 5-km grid cell was used as an offset in the model, making the response the time lag (between a specified photoperiod being reached and the date of the phenological event. The only fixed effect in this model was the intercept. For models that incorporated an effect of temperature, we followed Phillimore et al. (2010Phillimore et al. ( , 2016 and fitted both phenology and temperature as a bivariate response. This approach allowed us to separately model the relationship between phenology and temperature over space (across locations) and time (across years). The temperature response in the temp model was the mean temperature during a predefined sliding window. The start and end dates for the sliding windows were the same for all locations, and we tested different window durations by varying the start date (from ordinal days À59 to 100 in 2-day intervals) and duration (from 4 to 120 days in 2-day intervals). Each time window was constrained, so it did not extend beyond ordinal day 150 (30th May). The end of the time window was included as an offset for the phenology response and this generated a model to that tested whether temperature within a time window predicts the lag time until the phenological response is observed.
To model the combined effects of temperature and photoperiod (phototemp), we allowed sliding windows of thermal sensitivity to be initiated once a specified day length (using the same range of values as the photo model) had been reached. This date then became the start of the local time window, and we considered the same range of window durations as in the temp model.
The final model included two sliding windows during which mean temperatures predict phenological response (doubletemp), with both temperature variables and the phenological lag (between the end of the second time window and the phenological event) fitted in a trivariate response. The time window immediately preceding the event (the forcing Assumption Comments but note that short-lived species (e.g. garlic mustard) will have had more generations over which to adapt. vii. The temporal slope of phenology regressed on temperature is attributable to mean population plasticity.
Based on average Central England Temperatures for February-May (Parker et al., 1992), there has been little directional trend in UK spring temperatures over the period 1998-2014 (slope = À0.06 AE 0.03). For long-lived species, such as the focal tree and shrub species, the contribution of microevolution to the temporal slope is likely to be negligible. For these species, our assumption that this slope is attributable to plasticity is also supported by similar estimates being obtained for individual trees (Vitasse et al., 2010). Several of the focal species are short-lived perennials (herbs, grasses), and for these species, we cannot discount the possibility that microevolution contributes to the temporal slope and biases our estimate of b towards B. viii. Populations share the same plastic response.
When we estimate the temporal slope separately for each 150-km grid cell, we find little evidence for intraspecific geographic variation in mean population plasticity (Fig. S3). Mean population plasticity has also been found to vary little between sites for a sample of European trees (Vitasse et al., 2009b). ix. Observations are random samples from a population.
UKPN observations are of population first dates, which means that the individuals we are sampling have more negative intercepts than the population they are drawn from. First dates are also sensitive to sampling effort and species abundance; if either covaries with spring temperatures over time and/or space, this can bias any of our slope estimates up or down (appendix A in Phillimore et al., 2012).
window) was identical to the best performing temp model for each species. We then explored mean temperatures over a preforcing time window during the autumn/winter preceding the phenological event. For simplicity, this window is referred to as the chilling window, although the timing of this window could reflect temperatures that impact on phenology through a mechanism other than chilling, such as dormancy induction (Heide, 2003). We varied start dates (from ordinal day À120 up to the beginning of the forcing window in 2-day intervals) and durations (from 4 to 120 days in 2-day intervals) in all combinations.
We used Akaike information criteria (AIC) and AIC weights (Burnham & Anderson, 2004) to compare among the best models of each class (see Table 2 for the parameters included in the AIC calculation for each model). We explain the calculation of model likelihood in Appendix S1. The sliding-window method involves multiple testing which will inflate type I errors (Bailey & Van De Pol, 2016), although the very high autocorrelation in daily temperatures will serve to reduce the severity of this problem (M. Morrissey pers comm).
We obtained an estimate of the variance-covariance between response variables for each random effect (r): v Pr v 0 Pr;Tr v Pr;Tr V Tr v Pr is the variance in phenology, v Pr,Tr is a vector of covariance (s) between phenology and the temperature cue(s) and V Tr is a matrix of (co)variances between the temperature cue(s). In the bivariate model, v Pr,Tr and V Tr are scalars. The slope estimate(s) of the phenological lag on the temperature cue(s) was obtained as V À1 Tr v Pr;Tr for each random term (Phillimore et al., 2012). When year was the random effect, we obtained a temporal slope (i.e. the change in phenology in response to year-to-year variation in temperature), and when 150-km grid is the random effect, we obtained a detrended spatial slope (i.e. the change in phenology in response to nonclinal spatial variation in temperature).
We assumed that temporal slopes are primarily due to the mean population-level temperature-mediated phenological plasticity (b) (for discussion of the validity of assumptions required by the theory and statistical models, see Table 1). Following the approach of Phillimore et al. (2016), we estimated the temperature sensitivity of selection over latitude and longitude, which we refer to as B lat and B lon . Assuming that the temperatures in the selected thermal window cue phenology, populations are at migration-selection equilibrium and that population density is constant in space, B can be estimated by dividing the slope of phenological lag on latitude (or longitude) by the slope of temperature on latitude (or longitude) (Hadfield, 2016). We expect that B lat = B lon and any significant difference between these slopes is potentially indicative of an unmeasured confounding variable or nonconstancy of B in space. Assuming that plasticity is constant among populations, we can use B-b to estimate the contribution made by clinal local adaptation. When |B-b| < |B|, then plasticity partially tracks the optimum and can be said to be adaptive, and B = b indicates perfect adaptive plasticity. In addition to clinal local adaptation, nonclinal local adaptation can be estimated as the difference between the detrended spatial slope and b (Phillimore et al., 2010(Phillimore et al., , 2016. When |B-b| > 0, migration is expected to reduce the efficiency of adaptation to track B across temperatures that vary stochastically across grid cells (Hadfield, 2016). Therefore, we predict that the detrended spatial slope will lie between B and b.
To get credible intervals for slopes and slope differences, we selected the lowest AIC model for each species and reestimated the parameters in a Bayesian setting using MCMCglmm (Hadfield, 2010). For species where the phenological response was best explained by the temp or phototemp models, we ran MCMCglmm (Hadfield, 2010) using forcing windows from the best performing model. For species where phenology was best predicted by the doubletemp model, we tested the correlation of mean temperatures across the two time windows over time and space and found that temperatures in the two time windows were highly correlated over space but not time (correlation over space ranges from 0.57 to 0.99, Table S4). We interpret these models as being effective at identifying the time windows during which temperature is most important as a phenological cue, and the plastic response to these cues. However, multicolinearity precludes interpretation of forcing and chilling slopes estimated across spatially varying temperatures. Therefore, for these species we focused solely on parameter estimates for the forcing window and reestimated parameters from the best performing temp or phototemp model.
We ran models for 60 000 iterations, discarding the first 10 000 as burn-in. We sampled every 10th iteration to get posterior sample sizes of 5000 for each species and visually inspected traces of the posterior distributions of focal parameters to check for model convergence. We used priors for the Fig. 1 A schematic depicting latitudinal variation in cues arising from models that include (a, d) photoperiod and (b-d) average temperature in a sliding window. Parameters that we optimize via iterative searches are in blue. Lag* indicates models where the lag duration (which can be positive or negative) is a linear response to spatial and temporal variation in the mean temperature during the time window.
(co)variance components which were drawn from the inverse Wishart distribution with V = I and m = 0.002. Phenological observations and derived temperature variables for the best cues are available from https://doi.org/10. 7488/ds/1688.

Spatiotemporal trends
We find that the phenology of first leafing and flowering dates varies spatially among 150 km 9 150 km (hereafter 150 km) grid cells for all focal species (Fig. 2). Species differ substantially in their spatial phenological variance, with levels highest in lesser celandine, wood anemone and meadow foxtail, and lowest in field maple and beech. Variance among years is of similar magnitude to the variance among 150-km grid cells, and tends to be higher for species with earlier phenology such as lesser celandine, hawthorn and blackthorn. For all species, residual variance of first dates within a single 25 km 9 25 km grid cell and year is considerably larger than other variance components.
All species show significant latitudinal and/or longitudinal trends in phenology. With the exception of field maple and garlic mustard, a model where the effects of latitude and longitude interact and are subject to quadratic relationships (geo2, Table 2) outperforms a model that considers only linear effects and no interaction (geo1). For most species, phenology is delayed as latitude increases (Fig. S1), although the magnitude of this gradient varies, being steepest in bluebell and pedunculate oak and shallow in hawthorn, horse chestnut and beech. For elder, sycamore, rowan, garlic mustard and field maple, the latitudinal slope is reversed and phenology advances as latitude increases. Longitudinal trends vary among species, with some species being most advanced in the west and others in the east. Several species show longitudinal clines that change sign as one goes north; in most cases, phenology is earliest in the east in the south, and west in the north. Directional temporal shifts in phenology are nonsignificant for all species, consistent with the weak temporal temperature trend over the focal time period.

Cues
We use AIC to compare the ability of four types of cuebased model to predict spatiotemporal variation in phenology (Fig. 1, Table 2). While all focal species are sensitive to spring forcing (Fig. 3), they vary in whether they are sensitive to chilling or photoperiod and in the parameters defining the sliding windows. The single sliding-window temp model is preferred for meadow foxtail, with the more complex doubletemp model performing best for thirteen species, most of which are typified by early phenology. For eight predominantly late spring species, the phototemp model performs best. For all species, the spring forcing windows precede and overlap with spatiotemporal variation in the event itself. Forcing windows are earlier for species with earlier phenology; however, there are no clear trends in the length of forcing window according to best model type or timing of phenological event (Fig. 4, Fig. S2). For species with phenology best predicted by the phototemp model, forcing windows start and end later in the north (Fig. 4). The time delay in when the photoperiod threshold is met at 50°N as compared with 56°N ranged from 3 to 11 days. The latitudinal trend in the photoperiod initiation of the sliding-window start date becomes shallower approaching the equinox. For species with phenology best predicted by the doubletemp model, the preforcing or chilling temperature sensitivity window is generally towards the end of the year preceding the phenological event (Fig. 4).

Plasticity and adaptation
Our estimates of population mean plasticity with respect to forcing temperature (b), obtained as the temporal (among-year) slope of phenology on temperature during the forcing window, are significantly negative for all focal species (Fig. 5a, Tables S2-S4), with posterior medians varying from À3 to À10 days°C À1 . We find little evidence of among 150-km grid cell variance in plasticity (Fig. S3), except in lesser celandine, for which plasticity is estimated to be shallower at the temperature extremes. Our estimates of population mean plasticity with respect to chilling temperature (b chill ), which we obtain as the temporal (among-year) slope of phenology on temperature during the chilling window, are close to 0 in most cases. We do, however, find a significant positive b chill for larch, horse chestnut and sessile oak (Table S4). The temporal chilling slopes are shallower than the forcing slopes, ranging from À2 to 2 days°C À1 . Spatial multicolinearity between forcing and chilling temperatures precludes us from interpreting spatial slopes obtained from the doubletemp model. Instead, when interpreting drivers of spatial variation in phenology, we consider only forcing temperatures obtained from the best performing temp or phototemp model. The temperature sensitivity of the optimum across latitudes (B lat ) and longitudes (B lon ) is significantly negative in most cases (Fig. S4a, b), with a median of~À3 days°C À1 . Sycamore, for which we estimate a positive B across both spatial clines, is an exception. In 16 of 22 cases, B lat and B lon are of the same sign (Fig. S4a, b, Tables S2 and S3).
Four species, larch, sycamore, bluebell and garlic mustard, show significant differences between B and b that are qualitatively consistent whether B is estimated as B lat or B lon (Fig. 5b, c). For each of these species, the gradient of the optimum (B) is shallower than the plastic slope b, consistent with countergradient local adaptation (i.e. temperature-mediated local adaptation acting in the opposite direction to plasticity), or the effect of a third variable on the optimum (Chevin & Lande, 2015). For lesser celandine alone, B is more steeply negative than b, consistent with cogradient local adaptation (i.e. temperature-mediated local adaptation acting in the same direction as plasticity), although the credible interval for B lonb includes zero.
For seven species, the two estimates of Bb do not depart significantly from zero, consistent with temperature-mediated plasticity tracking clinal variation in the phenological optimum (Fig. 5b, c). In a further three species (horse chestnut, pedunculate and sessile oak), while there is a significant difference between B b over either latitude or longitude, the point estimate does not depart greatly from 0, implying that plasticity is adaptive and partially tracks the optimum in these species. For the remaining seven species, Bb estimates over latitudes and longitudes are inconsistent. Across all 22 species, we find that the point estimates of forcing plasticity are adaptive (defined as |B-b| < |B|) for 12 and 16 species when B is estimated as B lat and B lon , respectively. Plasticity is consistently maladaptive for bluebell, garlic mustard, larch and sycamore. Our measure of nonclinal local adaptation, the detrended spatial slope of phenology on temperature, is in most cases negative and, consistent with theoretical expectations (Hadfield, 2016), intermediate between B and b (Fig. 5d, Fig. S4c). In fifteen cases, the difference between the detrended slope and b was Bars are coloured according to the ability of plasticity to track the optimum | B-b|; grey = inconsistent signal, gold = consistent in direction, green = consistent with the hypothesis that plasticity tracks the optimum | B-b| = 0. Circles = phototemp model, triangles = temp model. Slopes for species whose phenology was best predicted by the doubletemp model were plotted using results from the lowest AIC alternative model and are represented by unfilled symbols. Species are in ascending order of mean phenology from left to right.
nonsignificant, meaning that we cannot reject the null hypothesis that the detrended spatial slope of phenology on forcing temperature is due to plasticity alone.

Discussion
The absolute difference between the temperature sensitivity of the optimum B (across latitude and longitude) and plasticity (b) reveals the contribution that genetic adaptation must make for populations to track the optimum (Chevin et al., 2010). For wood anemone, silver birch, alder, cuckooflower, beech, ash and cocksfoot, the difference between these slopes is small and nonsignificant. Plasticity also closely tracks gradients of thermal optima for the leafing of horse chestnut (a non-native species), pedunculate and sessile oak, and dogrose, although a small but significant slope difference exists across latitude or longitude for each of these.
Assuming that parameters estimated across space can substitute for those that act over time, we project that these populations will be able to track temperaturemediated changes in the phenological optima, and that, all else being equal, climate change should pose least threat to such populations (Chevin et al., 2010). In more than half of species, we find that plasticity is adaptive and leads to populations tracking temperaturemediated variation in the optimum more closely than they would in the absence of plasticity (|B-b| < |B|).
For four species, we found consistent evidence that |B-b| departs significantly from zero; that is, plasticity does not track the temperature sensitivity of the optimum. In each case, the plastic response was significantly steeper and more negative than B, which is consistent with countergradient variation in the optimum timing (Conover & Schulz, 1995;Phillimore et al., 2012). We expect the spring phenology of temperate plants to be exposed to opposing selection pressures, for later phenology to reduce frost damage and early phenology to take advantage of the growing season (Lenz et al., 2013). In colder locations, growing seasons are expected to be shorter and this may result in the optimum phenology being earlier than it would be under equivalent spring temperatures at lower latitudes or elevations. In other words, the intercept of the temporal relationship between spring temperatures and the optimum timing may vary in space and may make it inappropriate to use our estimate of B lat as a substitute for B. Evidence for countergradient variation in spring phenology from common garden experiments on plants is quite limited, although examples do exist (Alberto et al., 2013;Kremer et al., 2014;Toftegaard et al., 2015). Unexpectedly, we found countergradient patterns for two non-native species, larch and sycamore. For these species, we suggest that a more likely explanation for the countergradient patterns is that another variable may have played a confounding role (Chevin & Lande, 2015). For instance, it is conceivable that the lower amounts of chilling received by plants in warmer regions biases our estimates of B downwards. Unfortunately, the spatial correlation between forcing and chilling temperatures means that we are unable to separate B chill and B force . The early flowering lesser celandine is the only species for which we identify a cogradient local adaptation pattern (although we only find significant evidence for local adaptation over latitude), meaning that local adaptation of phenology with respect to temperature acts in the same direction as plasticity.
For the remaining ten species, our estimate of the temperature sensitivity of the optimum differs in magnitude over latitude vs. longitude, which is inconsistent with the underlying theory (Hadfield, 2016). This may reflect the influence of a third variable, such as the frequency of late frosts or precipitation, which covaries with temperature and phenology differently over latitude vs. longitude (see Table 1). We note that for several species for which spatial slope estimates were not consistent, such as hawthorn, blackthorn and rowan, the timing of temperature sensitivity sliding windows was estimated with a higher degree of uncertainty (Fig. S2).
Our analyses revealed a broad trend in cue use; species with earlier mean phenology are better predicted by two temperature time windows, while photoperiod is an important cue for species with later phenology. Exposure to late frosts and the damage this incurs can impair new growth and reproductive success (Inouye, 2000). Therefore, the positive phenological response to winter temperatures during a chilling window identified for 11 of 13 species (for which the doubletemp model performed best) may be an adaptation to reduce the chances of initiating new growth during a warm winter spell. A reliance on temperature rather than photoperiod cues may also enable these early phenology species to respond more quickly to warm forcing temperatures early in the year (Polgar & Primack, 2011). Chilling requirements have been demonstrated for numerous woody species (Laube et al., 2014) and flowering annuals (Kim et al., 2009). Our finding that early spring species are generally more sensitive to a chilling window agrees with inferences drawn for an overlapping set of species by an analysis of the famous Marsham phenological time series (Roberts et al., 2015).
The precise timing of all forcing and chilling windows varied among species. Forcing windows were important in cuing the phenology of all focal species and directly preceded the mean phenology. In comparison, the timing of chilling windows varied more among species. This may reflect different processes that constitute a chilling signal in our analyses, such as autumn dormancy induction or winter chilling accumulation for dormancy release (H€ anninen & Tanino, 2011;Polgar & Primack, 2011). Species that show a positive phenological response to temperature during a chilling window (i.e. delay their phenology) may advance their phenology by less than chilling insensitive species if temperatures rise throughout the year (Murray et al., 1989;Fu et al., 2015;Roberts et al., 2015).
For species with phenological events later in the spring, photoperiod assumes a greater influence than chilling as a predictor of a species' spatiotemporal phenological variation. This is in broad agreement with findings from growth chamber experiments on tree cuttings (Basler & K€ orner, 2012). We included photoperiod as a threshold, but recent experimental and modelling work on trees finds that photoperiod and chilling interact, such that photoperiod assumes greater importance when chilling requirements are not fully met (Caffarra et al., 2011a,b;Laube et al., 2014). Our models do not capture these more subtle effects. It may be possible to extend our approach to incorporate such complexity, although we caution that expanding the parameter space would present a substantial computational challenge and multicolinearity of cues in space will remain an obstacle.
Three of the species included here, pedunculate oak, garlic mustard and cuckooflower, have been subject to earlier work in a simpler version of our framework (Phillimore et al., 2012(Phillimore et al., , 2013. Our estimates of species' plasticity are similar to those obtained in previous studies. However, in contrast to our finding that garlic mustard shows countergradient local adaptation, Phillimore et al. (2012) reported no evidence of local adaptation and that plasticity could account for the spatiotemporal covariation between temperature and phenology. We attribute this discrepancy to the earlier study relying on a measure of model fit that took only the temporal relationship between temperature and phenology into account. This serves to highlight the risks of relaying on a correlation-based approach to identify local adaptation and B. We advise particular caution when interpreting B for species best predicted by the doubletemp model, as we suspect that the spatial correlation between forcing and chilling temperatures will have biased our estimates towards zero. For these reasons, we recommend that our findings be viewed as hypotheses requiring validation via the classic approach of transplant experiments.
We are aware of two opportunities for assessing the validity of our inferences. The first is the Marsham record, which allows a comparison of estimates of plasticity (or more accurately the temporal slope of phenology on forcing temperature), for ten taxa at a single site within the same region estimated over a nonoverlapping time period (Table S1b in Roberts et al., 2015). We identify similar forcing sliding windows and estimate plasticity of the same sign and of similar magnitude, but all our estimates are shallower, with the average difference~1.65 days°C À1 . This discrepancy may reflect a true difference in the phenological response to temperatures during the two time periods, but is more likely due to methodological biases. For instance, enforcing a single sliding window across the United Kingdom might underestimate the true local responses to forcing temperatures or the spatially interpolated temperature data that we use may include more measurement error. Several tree species have been subject to extensive transplant experiments in the Pyrenees, providing an opportunity to test the validity of our inferences regarding plasticity and local adaptation. Vitasse et al.'s (2010) estimates of the plasticity of leaf unfolding with respect to spring temperatures ranged from~À4.9 to À5.8 days°C À1 in beech and À5.7 to À6.3 days°C À1 in sessile oak, similar to our temporal slope estimate of À4.27 and À5.30 days°C À1 (Table S3), respectively. A common garden study of tree provenances from different elevations revealed countergradient local adaptation of flushing in beech, cogradient variation in ash and no local adaptation in sessile oak (Vitasse et al., 2009a). Our results for beech (a weak countergradient tendency across latitudes) and sessile oak are in broad agreement with this (Fig. 5). However, we find no evidence of local adaptation in ash and weak/absent countergradient variation in sessile oak (Fig. 5). While differences between our inferences into local adaptation and those arising from reciprocal transplant experiments may reflect true biological differences in the nature of adaptation to elevation vs. mesoscale geographic clines, or differences between the United Kingdom and Pyrenees, they underscore the need to interpret our findings with caution. In summary, we have shown how spatiotemporal data can be used to infer the sensitivity of plant phenology to forcing, chilling and photoperiod and estimate the ability of plasticity to track temperature-mediated shifts in the optimum timing. We find that for many UK plant species, temperature-mediated phenological plasticity is adaptive and will allow populations to partially or completely track shifts in optimum timing arising from increases in spring temperatures. While the statistical approach we present here relies on a large number of assumptions, we propose that it can be a useful tool for estimating parameters that are key to projecting population responses to climate change.

Supporting Information
Additional Supporting Information may be found in the online version of this article:

Fig. S1
Latitudinal and longitudinal trends in phenology. Slope estimates obtained from the lowest AIC simple clinal model (geo1 = green or geo2 = blue). Filled grid cells represent locations with records available. Species are plotted in ascending order of mean phenology. Fig. S2 Time windows that received substantial support (AIC i -AIC min <2 for a model type), during which mean temperatures best predict the phenological events (median phenology shown as filled circle) for each species. Windows corresponding to AIC min appear as bold bars, with other windows that received substantial AIC support appearing as translucent bars. Species are plotted in ascending order of mean phenology from bottom to top. Event type is reported in parentheses, where F = flowering and L = leafing. Bars are coloured according to the AIC min model type: orange = temp, red = phototemp and blue = doubletemp. Time windows for the phototemp model covary with latitude; the thick bar depicts the time window at 50°N and the thin bar the time window at 56°N. Fig. S3 Temporal slopes across 150-km grid cells calculated for the forcing time window that was identified by the best temp or phototemp model for each species (a-v). Only grid cells with at least 20 records and spanning a minimum of 10 years were included. Fig. S4 Posterior medians and 95% credible intervals for slopes of phenology on forcing temperature that correspond to the temperature sensitivity of the optimum phenology across (a) latitude = B lat , (b) longitude = B lon and (c) the detrended spatial slope across 150-km grid cells. Table S1 Species records selected for analyses from the UKPN data set. Table S2 Parameters for the temp model estimated via MCMCglmm. Species listed in ascending order of mean phenology, * indicates species for which temp is the best performing alternative to the doubletemp model, CI = credible interval. Table S3 Parameters for the phototemp model estimated via MCMCglmm. Species listed in ascending order of mean phenology, * indicates species for which phototemp is the best performing alternative to the doubletemp model, CI = credible interval. Table S4 Parameters for the doubletemp model estimated via MCMCglmm for forcing window and chilling windows. Species listed in ascending order of mean phenology, CIs = credible intervals. Note that B lat , B lon and slope differences are not reported due to the issues of multicolinearity in the doubletemp model. Appendix S1 Supporting methods. Appendix S2 Assessing the impact of spatial variation in population heterogeneity on B lat and B lon.