Phenological sensitivity to climate across taxa and trophic levels

Differences in phenological responses to climate change among species can desynchronise ecological interactions and thereby threaten ecosystem function. To assess these threats, we must quantify the relative impact of climate change on species at different trophic levels. Here, we apply a Climate Sensitivity Profile approach to 10,003 terrestrial and aquatic phenological data sets, spatially matched to temperature and precipitation data, to quantify variation in climate sensitivity. The direction, magnitude and timing of climate sensitivity varied markedly among organisms within taxonomic and trophic groups. Despite this variability, we detected systematic variation in the direction and magnitude of phenological climate sensitivity. Secondary consumers showed consistently lower climate sensitivity than other groups. We used mid-century climate change projections to estimate that the timing of phenological events could change more for primary consumers than for species in other trophic levels (6.2 versus 2.5–2.9 days earlier on average), with substantial taxonomic variation (1.1–14.8 days earlier on average). An ambitious study has used more than 10,000 datasets to examine how the phenological characteristics—such as the timing of reproduction—of various taxa alter in response to climate change, and suggests that differing levels of climate sensitivity could lead to the desynchronization of seasonal events over time. Variations in the phenological responses of different species to climate change have fuelled concerns that key species interactions may desynchronize over time, with consequences for ecosystem functioning. Stephen Thackeray et al. examine the climate sensitivity of 812 terrestrial and aquatic taxa across the United Kingdom, using more than 10,000 phenological data sets spanning 1960 to 2012, together with temperature and precipitation data. There was a systematic difference in the magnitude and direction of phenological climate sensitivity across trophic levels, despite marked heterogeneity among organisms sharing taxonomic affinities and trophic position. In particular, secondary consumers showed lower levels of climate sensitivity than primary producers and consumers. The authors suggest that the differential sensitivity of phenology to climate across trophic levels could result in the desynchronization of seasonal events in the future.

or resource availability 14 . Thus, even if species are most sensitive to climate change during the same seasonal period (window), they show different phenological responses to a given climatic change. Second, co-occurring species may vary in their seasonal periods of climate sensitivity, each being typified by different levels of directional climate change [15][16][17] . We conceptualize these two aspects of phenological responses as species-specific (or population-specific) climate sensitivity profiles (CSPs; Fig. 1). The CSP approach differs fundamentally from attempts to identify single 'critical' seasonal periods within which climatic change most strongly affects seasonal events 17 , by quantifying the full range of phenological responses to seasonal climatic change. We ask, "How sensitive are phenological events to temperature and precipitation change at different times of year?". By applying this approach to a large, taxonomically diverse, national-scale data set, we discern coherent patterns within a multitude of idiosyncratic biological climate responses. We assess whether systematic differences in climate sensitivity underpin differences in phenological change among taxonomic and trophic groups in the UK 8 .
We elected against using published climate responses that may be biased in favour of species showing an effect. Instead, we analysed 10,003 long-term (≥20-year) phenological time series for 812 marine, freshwater and terrestrial taxa over the period 1960-2012. Our data set aggregates many of the UK's foremost long-term biological monitoring schemes (Supplementary Table 1), including phenological information on amphibians (spawning), birds (egg laying, migration), planktonic Differences in phenological responses to climate change among species can desynchronise ecological interactions and thereby threaten ecosystem function. To assess these threats, we must quantify the relative impact of climate change on species at different trophic levels. Here, we apply a Climate Sensitivity Profile approach to 10,003 terrestrial and aquatic phenological data sets, spatially matched to temperature and precipitation data, to quantify variation in climate sensitivity. The direction, magnitude and timing of climate sensitivity varied markedly among organisms within taxonomic and trophic groups. Despite this variability, we detected systematic variation in the direction and magnitude of phenological climate sensitivity. Secondary consumers showed consistently lower climate sensitivity than other groups. We used mid-century climate change projections to estimate that the timing of phenological events could change more for primary consumers than for species in other trophic levels (6.2 versus 2.5-2.9 days earlier on average), with substantial taxonomic variation (1.1-14.8 days earlier on average).
Article reSeArcH crustaceans (population peaks), fish (spawning, migration), insects (flight periods), mammals (birth dates), phytoplankton (population peaks) and plants (flowering, fruiting, leafing). These taxa represent three broad trophic levels: primary producers (phytoplankton and plants), primary consumers (granivorous birds, herbivorous insects, mammals and planktonic crustaceans) and secondary consumers (predatory amphibians, birds, fish, insects, mammals and planktonic crustaceans). We spatially matched all 10,003 phenological time series with local temperature and precipitation data from a 5×5-km resolution gridded data set, before statistically modelling the relationship between seasonal timing and climatic variables. Between 1960 and 2012, mean UK air temperature increased in all months, and mean precipitation increased in most months (Fig. 2a).
Spatial variability in climatic change (Fig. 2b, c) necessitates local matching of phenological and climatic data sets rather than the use of regionally averaged climate data (for example, Central England Temperatures) or large-scale climatic indicators (for example, North Atlantic Oscillation). We did not make the restrictive assumption that biological events would be related to annual mean climatic conditions, or to conditions within periods based upon calendar months. Our CSP approach identified seasonal periods within which climatic change had the most positive and negative correlations with phenology (hereafter referred to as upper and lower limits of climate sensitivity, respectively). We could identify, for each phenological series, up to two seasonal periods within which climatic variation had a marked correlation with seasonal timing. The method was flexible enough to allow situations in which climatic variation within only a single period had a marked correlation with seasonal timing, and to identify seasonal windows ranging from a few days to a whole year in length. Our analysis captured the idiosyncrasies of phenological responses, allowed us to categorize them into generic types of climate response, and is consistent with current biological understanding of climate-phenology relationships 15,16 .

Climate response types in the UK
CSPs fall into three categories. The qualitative type of climatephenology correlation (positive or negative) may remain consistent, irrespective of when in the year climatic change occurs. In this case only the magnitude of the phenological response differs with the time of year at which climatic variables change. The climate-phenology correlation may be consistently negative (CSP type I; Fig. 1, red curve) or positive (CSP type III; Fig. 1, blue curve). Alternatively, opposing correlations between seasonal climatic change and the timing of biological events may exist; that is, the direction and magnitude of the phenological response may vary (CSP type II; Fig. 1, orange curve). We determined CSPs for responses to temperature (CSP temp ) and precipitation (CSP precip ).
Among responses to temperature changes, CSP type II was most common (Extended Data Table 1; 69.7% of phenological series). Seasonal events were advanced by (that is, negatively correlated with) warming during one period of the year, and delayed by (that is, positively correlated with) warming in another period. After multiple testing correction, 44.8% of the observed phenological advances (but only 1.0% of delays) with warming were statistically significant (P<0.05). CSP type I was the next most common response type: warming in different seasonal windows was consistently correlated with earlier seasonal events (that is, negative correlations; 24.7% of series). In this case the lower and upper limits of CSPs represent the 'strongest' and 'weakest' phenological advances with warming, respectively, and 58.1% of the 'strongest' responses were statistically significant (P<0.05, correcting for multiple testing).
Phenological events most commonly demonstrated opposing ( Fig. 1; CSP type II, 53.0% of series) or consistently positive ( Fig. 1; CSP type III, 28.0% of phenological series) correlations with increasing seasonal precipitation. Although delayed phenological events may commonly be associated with higher precipitation (81.0% of events show this type of response), few of these associations were statistically significant (Extended Data Table 1).

Climate sensitivity at the UK-wide scale
We matched each phenological series with four climate variables: mean temperature during the seasonal windows at the upper and lower limits of CSP temp , and similarly averaged precipitation data for the seasonal windows at the upper and lower limits of CSP precip . We then combined all 10,003 phenological series and their matched climate data, and modelled the relationships between seasonal timing and climate variables using linear mixed effects (LME) models. Initially we fitted a 'global' model to quantify the upper and lower limits of temperature and precipitation sensitivity, averaged across all phenological events. Marine plankton data were excluded at this stage, owing to a lack of precipitation data.
Most phenological events occurred earlier with seasonal warming (average rate −2.6 days per °C; Fig. 3a and Extended Data Table 2). Variation in the strength of this correlation was similar among sites and species (random-effects variances in site and species level seasonal timing-temperature slopes were 2.1 and 1.9, respectively). Some phenological events occurred later with seasonal warming (Fig. 3a) although, in other cases, the upper limit of CSP temp was in fact a 'weak' advance with warming. The upper limit of temperature sensitivity was more variable among species than among sites (random effects variances in species and site level seasonal timing-temperature slopes were 2.3 and 0.4, respectively). Averaged across species and populations, temperature responses were most consistent with CSP type II.
Most phenological events showed opposing responses to increasing seasonal precipitation ( Fig. 1; CSP    Article reSeArcH delays with rising precipitation was greatest: the average upper limit of CSP precip exceeded the lower limit (1.4 days per mm and −0.4 days per mm, respectively; Fig. 3b and Extended Data Table 2). The upper limit of CSP precip was more variable among species than among sites (species and site level random-effects variances in the seasonal timingprecipitation slopes were 1.9 and 1.2, respectively). The fitted climatephenology model was better supported by the data than a year-only model with the same random effects structure (ΔAIC (Akaike's information criterion) 293,516). This indicates the presence of real associations between climate and seasonality, rather than purely spurious correlations resulting from shared temporal trends. Average sensitivity to temperature was very similar in the model that included marine plankton data, but excluded precipitation effects (see Supplementary Discussion and Extended Data Fig. 1).

Taxonomic and trophic group sensitivity
We tested the hypothesis that the limits of seasonal climate sensitivity differ coherently among taxonomic groups by including a fixedeffect interaction between taxonomic group and each climatic variable ( Fig. 4 and Extended Data Table 2). The lower limit of CSP temp was negative for all groups ('earliness' with warming), the strongest responses being found for plants, freshwater phytoplankton, insects and amphibians (4.3, 4.1, 3.7 and 3.4 days earlier per °C, respectively). The upper limits of CSP temp indicated that freshwater phytoplankton and mammals experienced the greatest phenological delays with seasonal warming (2.9 and 2.0 days later per °C, respectively) but that plants showed little evidence of such delays. The strongest phenological delays with rising seasonal precipitation were found for freshwater phytoplankton and insects (2.5 and 2.2 days later per mm, respectively), while freshwater phytoplankton also exhibited the strongest phenological advances with rising precipitation during other seasonal windows (1.1 days earlier per mm). Average temperature and precipitation responses were consistent with CSP type II in most cases. There was considerable within-group variability in sensitivity. We examined trophic-level differences in climate sensitivity by including trophic level in interaction with each climate variable in the global model. The lower limit of CSP temp showed greater systematic variation than the upper limit among trophic levels (Fig. 3c, e). The tendency towards 'earliness' with seasonal warming was strongest at lower trophic levels (−4.1, −3.7 and −1.9 days per °C for primary producers, primary consumers and secondary consumers, respectively; Extended Data Table 2), consistent with observations of more rapid phenological changes at lower trophic levels, in the UK 8 . Conversely, the lower limit of CSP precip varied less than the upper limit among trophic levels ( Fig. 3d, f). The tendency for seasonal events to be later with higher seasonal precipitation was greater for primary producers and primary consumers (1.8 and 2.2 days per mm on average, respectively) than for secondary consumers (1.0 days per mm). Variations in climate Article reSeArcH sensitivity were described more parsimoniously by taxonomic groups than by trophic levels (AICs of taxonomic and trophic-level models 3,237,611 and 3,238,061, respectively).
The results were affected little when we analysed only pre-and post-1980 data, to minimize among-group variation in time series length, and after Monte Carlo re-sampling to assess the potential effects of taxonomic bias (see Supplementary Discussion and Extended Data Figs 2-4). The same qualitative trophic-level differences in climate sensitivity were apparent when we included marine plankton data in a temperature-only LME model (see Supplementary Discussion and Extended Data Fig. 1). In contrast to trophic-level differences in the magnitude of sensitivity, there was little evidence of similar variation in the seasonal timing of climate sensitivity (see Supplementary Discussion and Extended Data Figs 5-7).

Estimating future change
Overall, net phenological responses to climatic change combine potentially opposing responses to conditions in different seasonal periods. We estimated net responses by the 2050s by applying our fitted models to UKCP09 probabilistic projections (bias-corrected relative to a 1961-1990 baseline) of temperature and precipitation change under low-, medium-and high-emissions scenarios. Rather than predicting the absolute timing of future phenological events, we contrasted possible changes in seasonal timing among organism groups based upon established climate scenarios and contemporary patterns of climate sensitivity. Estimated average phenological changes were less for primary producers and secondary consumers than for primary consumers (Fig. 5a). This occurred because, averaged across species, the opposing climate responses of primary producers and secondary consumers are more similar in magnitude than are those of primary consumers (Fig. 3), effectively cancelling each other out. Our models suggest greater average advances for crustacea, fish and insects than for other groups, such as freshwater phytoplankton, birds and mammals (Fig. 5b). However, response variation is high for crustacea (Fig. 5b).

Discussion
In the UK, phenological climate sensitivity varies greatly, suggesting that it is influenced by locally varying, non-climatic drivers such as population structure 18 , resource availability 19 and adaptation 20 . This is relevant to the use of phenological change as a tangible climate change indicator 1,21 . Mediators of phenological climate sensitivity are known only locally for some of the groups in our data set (for example, nutrient availability for freshwater phytoplankton) 22 . However, for others, the climate sensitivity of different biological traits is known to be mediated by alternative drivers 23,24 . High climate response variability necessitates wide site and species coverage in long-term monitoring schemes designed to develop robust aggregate indicators of change 21 . As climatic conditions are more spatially variable across broader geographic domains, site-level replication of phenological monitoring is particularly important when interpreting phenology at continental-to-global scales. In the UK, average responses for fish and insects appear to be sensitive indicators of climate effects. These groups show consistently strong phenological advances with seasonal warming, and only weak opposing responses, resulting in relatively large (net) changes in seasonal timing. Interpretation of phenological changes for other groups is more complex. For example, freshwater phytoplankton show strong evidence of opposing phenological responses to climatic variation at different times of year and these responses are near-equivalent in magnitude, such that estimated net changes are negligible. This finding emphasizes that long-term observations represent the net effect of potentially opposing biological responses 25 . To fully capitalize on the indicator potential of phenological change, we must advance our mechanistic understanding of responses to potentially opposing climatic and non-climatic drivers.

Article reSeArcH
Despite this variability, we identified coherent patterns in climate sensitivity among the idiosyncratic responses of many wild plant and animal populations. We have shown that, on average, species in different trophic levels differ in the magnitude of seasonal climate sensitivity, but not in the time of year within which climatic change has its most pronounced effects. This may be a key mechanism underpinning observations of trophic level differences in phenological change in the UK 8 . Lower trophic levels demonstrated more pronounced variation in their sensitivity to changing temperature and precipitation at different times of year, and stronger phenological responses to climatic change during defined (taxon-and population-specific) seasonal periods.
In response to climatic changes projected for the 2050s, relative changes in seasonal timing are likely to be greatest for primary consumers, particularly in the terrestrial environment. The difference in magnitude between opposing climate responses is greatest for primary consumers, resulting in greater net change. Our approach makes the simplifying assumption that climatic change has an overriding influence upon seasonality. Nevertheless, our results suggest that systematic differences in climate sensitivity could result in widespread phenological desynchronization. However, factors that shape phenological climate responses introduce uncertainty into projections of future phenological change. These results should catalyse research to improve predictive capacity in the face of multiple environmental and demographic drivers that not only mediate rates of change, but might also confer resilience to desynchronization (for example, population density dependence 26 , compensatory range shifts 27 , and the formation of novel inter-specific interactions 28,29 ). These findings also underscore the importance of developing our capacity to manage ecosystems within a 'safe operating space' with respect to the likely impacts of projected climate change 30 .
Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper.

MethOdS
Data sets. We integrated data from many major UK biological monitoring schemes (Supplementary Table 1), resulting in 10,003 long-term (at least 20 years between 1960 and 2012) phenological series for 812 marine, freshwater and terrestrial taxa. The amassed data sets included records for plants, phytoplankton, zooplankton, insects, amphibians, fish, mammals and birds (379,081 individual phenological observations). For each study we used a single population-level phenological measure per year (Supplementary Table 1). Because the sampling resolution for the marine plankton data was monthly, before analysis we re-scaled these data into units of days. Trophic level, taxonomic class and environmental affinity were assigned to each taxon, to permit analyses of correlations between these attributes and climate sensitivity.
Daily air temperature and precipitation data were extracted from the Met Office National Climate Information Centre (NCIC) 5-km-resolution gridded data set 31 for the spatial locations of all biological monitoring sites across the UK land surface. If available, recorded water temperatures from the same sites were used in place of air temperatures for phenological time series representing obligate aquatic taxa (freshwater plankton and fish). Water temperatures were interpolated onto a daily time-step before analysis 32 . If these data were not available, daily water temperature data were estimated from air temperatures using a fitted empirical site-specific relationship between air and water temperature. For the sea trout (Salmo trutta) data, an existing linear relationship 33 was used, while for the Atlantic salmon (Salmo salar) data, a nonlinear relationship 34 was calculated for a nearby river, the Tarland Burn, and applied to air temperatures from the sampling site. For the marine plankton, mean monthly sea surface temperatures were extracted from the Met Office Hadley Centre Sea Ice and Sea Surface Temperature (HadISST) data set 35 for each of the Standard Areas 36 in which phenological data were available. Precipitation data were not available for marine Standard Areas. Statistics. Our analysis was conducted in two distinct phases (see Supplementary  Notes). First, the CSP for each phenological series was estimated using generalized linear models to quantify associations between the timing of seasonal events and mean temperature and precipitation (within defined seasonal time windows) at the same location. Second, the phenological time series were aggregated and a single LME model was run, capturing upper and lower limits of climate sensitivity across many species. CSPs for precipitation were not estimated for marine plankton data (see above), so the second-phase LME models were run twice: once to examine correlations with temperature and precipitation for all but the marine plankton phenological series (9,800 series), and once to examine only correlations with temperature for the whole data set (10,003 series). Phase 1: Estimating CSPs for each time series. We used consistent methods to 'screen' all phenological events with respect to their climate sensitivity, finding periods of the year in which temperature and precipitation had the most positive and negative correlations with seasonal timing (the upper and lower limits of climate sensitivity). This approach was flexible enough to detect when these limits represented opposing correlations between temperature or precipitation and seasonality, depending upon the seasonal timing of climatic change; for example, spring warming may advance budburst, but winter warming may delay it 37 ( Fig. 1; CSP type II). This approach could also detect when the direction of the correlation between climatic variables and seasonal timing was consistent irrespective of the seasonal timing of climatic change, with only the magnitude of the correlation varying between the limits of the CSP ( Fig. 1; CSP types I and III).
For each phenological time series, we calculated the day of the year by which 95% of the recorded seasonal events had occurred (doy 95 ). Inter-annual variations in seasonal timing were statistically modelled as a function of daily mean temperatures on doy 95 each year. Then, a series of 365 statistical models was run that instead used daily mean temperatures on doy 95 -1 to doy 95 -365 as predictors. Slope coefficients and R 2 values for the temperature terms in these models were collated, capturing seasonal variations in the sign and magnitude of the phenology-temperature relationship (that is, the CSP; Fig. 1). Generalized linear models (GLMs) were used.
For two data sets (BTO Nest Record Scheme and PTES National Dormouse Monitoring Scheme; Supplementary Table 1) we modified the above analytical framework. In both of these schemes, the precise location of the biological observations changed among years (compared with other schemes in which monitoring sites were static over time). We extracted matching climatic data for each specific location in each year, as for all other schemes, but then grouped the phenological and climatic data at county level (mean area 3,440 km 2 ). Then, for each taxon in each county we used the fixed-effect slope parameters and R 2 values from a series of LME models, instead of GLMs, as a basis for estimating CSPs. In these models, we included fixed effects of temperature on doy 95 to doy 95 -365 as before, and included a year random effect to account for replicate phenological records for each taxon in each county in each year. For the SAHFOS marine plankton data set, we modified our iterative approach to analyse seasonal timing-temperature relationships at monthly, instead of daily, time steps (the temporal resolution of the sea surface temperature data).
As a final step in estimating the CSP for each series, temporal variation in the sign and magnitude of the seasonal timing-temperature correlation was itself modelled (Extended Data Fig. 8). This was done by fitting generalized additive models (GAMs, gamma error distribution) to the time series of slope coefficients and R 2 values from the models described above. By smoothing these time series, the GAMs identified periods of the year in which slope coefficients were consistently negative (that is, warming advanced seasonal timing), or consistently positive (that is, warming delayed seasonal timing), and during which the climate-phenology models generating the slope estimates had a their highest goodness-of-fit.
Seasonal windows in which the upper and lower limits of temperature sensitivity occurred were identified as periods during which: 1) the 95% confidence interval for the GAM fitted to the slope coefficients surpassed the limits of the 2.5th and 97.5th percentiles of the original slope coefficients; and 2) the 95% confidence interval for the GAM fitted to the R 2 values surpassed the 97.5th percentile of the original R 2 values. This ensured that seasonal windows were defined by periods combining the greatest climate effect size and relatively strong predictive power (determined by R 2 ). Using this framework, we identified the lower limit of CSP temp : the period of the year in which an advancing effect of increasing temperature upon seasonal timing was most likely. This was estimated by determining when the 95% confidence interval of the GAM intersected the lower percentile of the seasonal timing-temperature slope coefficients, by tracking the most negative coefficients (Extended Data Fig. 8). In addition, we identified the upper limit of CSP temp by determining when the 95% confidence interval of the GAM intersected the upper percentile of the seasonal timing-temperature slope coefficients, by tracking the most positive (or least negative) coefficients. Excluding the marine plankton data, the whole modelling process was repeated with precipitation as a predictor instead of air temperature, culminating in the estimation of seasonal periods capturing the limits of phenological responses to changing precipitation.
After this process, temperature and precipitation were each averaged within the two seasonal windows in which the limits of phenological sensitivity occurred. With the exception of the marine plankton data, the final seasonal timing-climate model for each series was then fitted using a GLM with gamma error distribution including four predictors: inter-annual variations in 1) mean temperature during the period at the lower limit of CSP temp , 2) mean temperature during the period at the upper limit of CSP temp , 3) mean precipitation during the period at the lower limit of CSP precip , 4) mean precipitation during the period at the upper limit of CSP precip . For the marine plankton data, only the first two terms were fitted. For the BTO Nest Record and PTES National Dormouse Monitoring Scheme data sets we implemented these final models in a mixed effects framework with a random effect of year, as before. Therefore, although we modelled changes in statistical parameters (which are not estimated without error) to identify seasonal periods, this step was used only to find the original climatic data to be used in subsequent modelling. Inferences were not, therefore, directly based upon statistical modelling of uncertain parameter estimates. We categorized the results of all 10,003 CSPs according to three broad response types (CSP types I-III; Fig. 1), and retained P values for each fitted model term to infer which of the modelled climatic effects were statistically significant. We examined the evidence for trophic-level differences in the mean seasonal timing of climate sensitivity by modelling the relationship between the start date, end date and duration of the seasonal windows capturing the upper and lower limits of phenological sensitivity to temperature and rainfall as a function of trophic level (fixed effect), with random effects of phenological metric, within species, within site. Analyses were conducted using the base, mgcv and lme4 packages in R (refs 38-40). Phase 2: Global models of phenological climate sensitivity. We estimated the upper and lower limits of phenological climate sensitivity at a multi-species scale by matching each phenological series with data on mean temperature and precipitation, during the seasonal windows characterizing the CSP for that series (Phase 1, above). We aggregated all 10,003 of these matched phenology-climate data sets. To quantify the average, multi-species, upper and lower limits of climate sensitivity we constructed an LME model in which phenology (day of year) was modelled as a function of mean temperature and precipitation within the seasonal windows of the amassed CSPs (fixed effects) with random effects of phenological metric, within species, within site. These random effects were necessary because our data could not be considered independent. The timings of events are more likely to be similar for the same species than for different species; the same is true for different sites and the phenological metric-types used to describe the events (for example, first flight time or seasonal peak abundance). Random slopes and intercepts were allowed to ensure that each phenological event, for a species at a site, was allowed a different rate of climate response.
For some species, more than one phenological event was recorded in the same year, at the same site. For example, butterflies may have more than one flight period in the same year, and plankton populations may be characterized by more than one seasonal abundance peak. As climate responses are unlikely to be the same for the first event of the year and subsequent events, we introduced a voltinism factor in the analysis. This allowed us to distinguish between data representing the first or only events of each year (for example, a spring plankton bloom or butterfly generation) and second events in each year (for example, the subsequent summer plankton bloom or butterfly generation). This distinction captured all possibilities within our data set.
For site i, species j, voltmetric k (where voltmetric is a unique combination of voltinism class and the metric-type used to identify the event), the corresponding day of year (DOY) of a particular seasonal event is modelled as: where ε ijk ~ N(0, σ 2 ) and the model includes temperature at the upper limit of each CSP (TU), temperature at the lower limit of each CSP (TL), precipitation at the upper limit of each CSP (PU) and precipitation at the lower limit of each CSP (PL). Owing to the non-independence within the data, we allow the intercepts and coefficients corresponding to all four covariates to vary by site, species and voltmetric. Preserving the natural nesting of a metric for a species at a particular site, this gives: where each of the μ terms is a normally distributed random effect. This nesting of random effects is most conservative in terms of inference at the global level and is as flexible as possible, allowing each time series to have its own set of model parameters. This permits a high degree of biological realism because each distinct phenological event, for a given species, at a given site, is permitted to have a different slope for the effects of temperature and precipitation (that is, a different climate sensitivity).
In this model framework we are specifically testing the null hypotheses that each of the climate variables shows no relation to the seasonal timing of biological events. Because of this, and the fact that each parameter is estimated directly, without distributional form assumed a priori or as the target distribution, we follow a frequentist approach to analysis. However, because the exact degrees of freedom cannot be evaluated when using restricted maximum likelihood, and hence no exact P values can be obtained, we present full summaries of all the parameters estimated at species level (as given by γ + μ ij,k + μ i,jk , above). Approximate P values could be presented by taking conservative estimates of the degrees of freedom although, given the volume of data available, this will typically lead to the detection of many statistically significant results that may not be biologically significant. Examining the full range of estimated coefficients across the random effects levels ensures that we present the full range of variation around global parameters and can make more informed inferences. In this way we encourage the reader to interpret our results by using biological insight, not by depending upon P values alone.
To examine high-level differences in climate sensitivity among trophic levels and taxonomic groups, we re-fitted the LME model with these attributes as fixedeffect factors, interacting with the fixed-effect climate variables. The fixedeffect slopes from the resulting models allowed us to compare differences in phenological climate sensitivity among these broad organism groups, averaged across all taxa within each group. Supplementary Table 2 shows the number of phenological series, sites and distinct taxa that contributed data to each of these groups. All models were run twice: once to examine correlations with both temperature and precipitation excluding marine plankton data (9,800 time series), and once to examine only temperature-phenology correlations for the whole data set (10,003 time series). Potential biases. Data availability differed among taxonomic groups. To assess the extent to which mean responses were biased by data inequality we conducted Monte Carlo re-sampling, iteratively selecting 5, 20, 50 and 100 phenological series from each taxonomic group and re-fitting climate-phenology models with these sampled data sets. For taxonomic groups with fewer data than the larger sample sizes, we retained all available data (see Supplementary Discussion). This allowed us to compare taxonomic group and trophic level responses based on sampled and all data, to fully investigate potential bias.
Another potential bias in our analysis is that phenological time series length is variable, affecting the length of time over which climate-phenology correlations are assessed. In order to assess the extent to which differences in mean trophic level and taxonomic group responses are biased by variable time series length, we also re-fitted our models but based only on pre-and post-1980 data. All models were run in the lme4 package in R (refs 38, 40). Estimating future change. To estimate the potential future net effects of temperature and precipitation change, we compared predictions of seasonal timing under baseline conditions, and under established climate change scenarios. First, estimates of seasonal timing (day of year) were obtained for the same baseline period used in the UKCP09 projections (long term average 1961-1990), using modelled correlations between phenology, temperature and precipitation (from Phase 1). Having obtained these baseline estimates, we applied our models to projected changes in monthly temperature and precipitation for the 2050s (UK Climate Projections, UKCP09; http://ukclimateprojections.metoffice.gov.uk/). We used 10th, 50th and 90th percentile changes under low-, medium-and high-emissions scenarios (relative to the 1961-1990 baseline). The spatial location of each phenological series was matched to climate projection data for the 25 × 25-km grid square in which it occurred, and temporally matched to climatic data from the months of year in which its respective climate sensitivity windows occurred. Relative changes in timing, in response to climatic change of the magnitude projected to occur by the 2050s, were summarized by trophic levels and taxonomic groups. Code availability. R code to run the described analyses can be found on GitHub (https://github.com/NERC-CEH/Phenology_Climate). and precipitation (d, f) sensitivity are shown by trophic level. Inverted triangles: average sensitivity for all species (a, b) or trophic levels (c-f). Curves, kernel density plots: probability density distributions of specieslevel climate sensitivity (that is, the relative likelihood of different climate sensitivities within each species group) (n = 370,725).