Front. Plant Sci. Frontiers in Plant Science Front. Plant Sci. 1664-462X Frontiers Media S.A. 10.3389/fpls.2019.00398 Plant Science Original Research Daily Maximum Temperatures Induce Lagged Effects on Leaf Unfolding in Temperate Woody Species Across Large Elevational Gradients Bigler Christof 1 2 * Vitasse Yann 2 3 1Forest Ecology, Institute of Terrestrial Ecosystems, Department of Environmental Systems Science, ETH Zürich, Zurich, Switzerland 2SwissForestLab, Birmensdorf, Switzerland 3Disturbance Ecology, Forest Dynamics, Swiss Federal Research Institute WSL, Birmensdorf, Switzerland

Edited by: Heikki Hänninen, Zhejiang Agriculture and Forestry University, China

Reviewed by: Qianqian Ma, South China Botanical Garden (CAS), China; Peter Prislan, Slovenian Forestry Institute, Slovenia

 *Correspondence: Christof Bigler, christof.bigler@env.ethz.ch

This article was submitted to Functional Plant Ecology, a section of the journal Frontiers in Plant Science

 28 03 2019 2019 10 398 12 12 2018 15 03 2019 Copyright © 2019 Bigler and Vitasse. 2019 Bigler and Vitasse

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The timing of leaf unfolding in temperate woody species is predominantly controlled by the seasonal course of temperature in late winter and early spring. However, quantifying lagged temperature effects on spring phenology is still challenging. Here, we aimed at investigating lagged and potentially non-linear effects of daily maximum temperatures on the probability of leaf unfolding in temperate woody species growing across large elevational gradients. We analyzed 5280 observations of leaf-out time of four tree species (European beech, horse chestnut, European larch, Norway spruce) and one shrub species (common hazel) that were recorded by volunteers over 40 years at 42 locations in Switzerland. We used a case-crossover sampling design to match leaf-out dates with control dates (i.e., dates before or after leaf-out), and analyzed these data with conditional logistic regression accounting for lagged temperature effects over 60 days. Multivariate meta-analyses were used to synthesize lagged temperature and elevational effects on leaf unfolding across multiple phenological stations. Temperature effects on the probability of leaf unfolding were largest at relatively short lags (i.e., within ca. 10 days) and decreased with increasing lags. Short- to mid-term effects (i.e., within ca. 10 to 20 days) were larger for late-leafing species known to be photoperiod-sensitive (beech, Norway spruce). Temperature effects increased for the broadleaved species (horse chestnut, hazel, beech) with decreasing elevation, particularly within ca. 10 to 40 days, i.e., leaf unfolding occurs more rapidly at low elevations for a given daily maximum temperature. Our novel findings provide evidence of cumulative and long-term temperature effects on leaf unfolding, whereby the efficiency of relatively high temperatures to trigger leaf-out becomes higher shortly before bud burst. These lagged associations between temperature and leaf unfolding improve our understanding of phenological responses across temperate woody species with differing ecological requirements that occur along elevational gradients.

 phenology conifers broadleaved species maximum temperature distributed lag models lag effects multivariate meta-analysis elevation Introduction

In temperate ecosystems, recurring biological phenomena are predominantly controlled by the cyclic, seasonal course of weather conditions (Forrest and Miller-Rushing, 2010; Schwartz, 2013). For example, temperature is a key environmental driver of the timing of trees’ spring phenophases such as leaf unfolding or flowering (Chuine, 2000; Polgar and Primack, 2011), with photoperiod playing an additional role for some species (Basler and Körner, 2012; Way and Montgomery, 2015). The strong dependency of phenological processes on temperature is further reflected in the phenological shift of trees’ leaf-out time along elevational gradients where temperature varies strongly over short distances (Vitasse et al., 2009; Bigler and Bugmann, 2018; Vitasse et al., 2018).

Leaf-out is typically observed when relatively warm spring temperatures prevail for a certain time. While one or a few warm days in late winter or early spring do not immediately result in leaf unfolding, a series of warm days in late spring may lead to premature leaf-out. Thus, the processes that ultimately induce leaf unfolding have already started several weeks or months before. Prior to leaf-out, vegetative buds of temperate woody species undergo a period of dormancy during which visible growth of buds is temporary suspended (Lang et al., 1987). Bud dormancy has been commonly divided into three physiological phases (Lang et al., 1987; Horvath et al., 2003): First, paradormancy from summer to early fall, which is imposed by distal hormonal control from outside the buds but within the plant (Cline and Deppong, 1999). Second, endodormancy from early fall to mid-winter, which is induced by chilling temperatures and shorter days and further regulated by physiological mechanisms occurring within the buds (Horvath et al., 2003). Endodormancy is released by exposure to cool temperatures in mid-winter (“chilling” temperatures between ca. 0 to 10°C; Polgar and Primack, 2011; Hänninen, 2016), and for some species photoperiod may also play a role (Caffarra et al., 2011a,b). Third, unfavorable environmental conditions (e.g., low temperatures, short days or drought) during ecodormancy prevent the buds from leafing out, until warm temperatures (“forcing” temperatures above ca. 5°C; Chuine, 2000; Polgar and Primack, 2011) and increasing day length in early spring induce bud burst and leaf unfolding. Chilling and forcing temperatures are assumed to be negatively related, i.e., a longer period of chilling requires a shorter period of forcing until bud break occurs (Murray et al., 1989; Laube et al., 2014). A longer photoperiod further increases the sensitivity of buds to forcing temperatures for some species such as European beech (Basler and Körner, 2014). Formulating and testing hypotheses about potential temperature effects on leaf unfolding remains challenging due to (1) the contrasting effects of chilling and forcing on bud development (Cannell and Smith, 1983; Murray et al., 1989), (2) the loosely defined ranges of chilling and forcing temperatures (Chuine, 2000; Luedeling et al., 2013; Delpierre et al., 2016), (3) the gradual transitions between or the overlap of the three dormancy phases (Cooke et al., 2012), (4) the relative importance of photoperiod versus temperature (Chuine et al., 2010; Körner and Basler, 2010; Vitasse and Basler, 2013), and (5) the correlation between temperature and photoperiod due to the similar seasonal course.

Because weather conditions in spring can drastically change from year to year (Körner and Basler, 2010; Vitasse et al., 2014b), dates of leaf unfolding may differ by one or even 2 months at the same site (Sparks and Carey, 1995; Lenz et al., 2016). Leaf unfolding is further known to be delayed at higher elevations (Vitasse et al., 2018), where it occurs after less accumulated warmth compared to lower elevations (Güsewell et al., 2017). To improve our understanding of leaf phenological processes, we would actually require to link continuous measurements of the development, growth or physiological activity of vegetative buds from bud set in summer to bud burst in spring with high-resolution weather data (Cooke et al., 2012). However, molecular or physiological changes that would reflect transient phases during bud dormancy remain to be found. Furthermore, long-term measurements of leaf phenology are typically restricted to observed dates of bud burst or leaf unfolding only. Thus, relating these phenological events to the driving forces such as lagged temperature effects is notoriously challenging.

A large range of statistical phenology models have been developed to predict the leaf-out time of forest trees in response to weather variability (Hänninen, 2016). The simplest statistical approach uses a linear model to relate the time of leaf unfolding (e.g., day of year) to some aggregated measure of temperature (e.g., growing degree days, mean temperature in late winter and spring, or monthly temperature variables; Sparks and Carey, 1995; Menzel et al., 2006; Ibáñez et al., 2010; Roberts, 2012). A drawback of this approach is that the period needs to be fixed over which monthly or seasonal weather variables are considered. Because the leaf-out time varies in response to the year-to-year changes of temperature and along elevation (Güsewell et al., 2017), fixed periods are likely to increase uncertainty in these models or may induce artifacts, e.g., if weather conditions following mean leaf unfolding dates are considered that turn out to have significant influences in specific years only. Unlike these linear models that include relatively coarse weather variables as predictors, more complex statistical approaches have been devised that rely on daily weather data such as (i) partial least squares regression (Roberts, 2008; Luedeling and Gassner, 2012); (ii) penalized signal regression (Roberts, 2008; Hudson, 2010; Roberts et al., 2015); or (iii) hierarchical state-space models (Clark et al., 2014a,b). Although these statistical models consider temperature data prior to leaf unfolding, they do not explicitly model lagged temperature effects, which may yet improve our understanding of how temperature regulates leaf unfolding (Hudson, 2010).

In addition to statistical phenology models, many process-based models have been proposed to predict spring phenology based on the hypothetical processes occurring during the endo- and ecodormancy (Chuine et al., 2013; Basler, 2016; Hänninen, 2016). These models may be classified as whether they account for (i) forcing temperature in spring only; (ii) chilling temperature during winter and forcing temperature in spring; or (iii) additionally photoperiod as an interacting factor during both the endodormancy and ecodormancy phases. These process-based models do not explicitly account for lagged daily weather effects, but rather integrate daily or hourly temperatures over a defined period. Besides, simpler process-based models often perform as well as or even outperform more complex models that require more parameters and are therefore often overparameterized, reflecting the poor understanding of the processes regulating spring phenology (Linkosalo et al., 2008; Vitasse et al., 2011; Clark et al., 2014b; Olsson and Jönsson, 2014; Basler, 2016).

The goal of this study was to investigate lagged and potentially non-linear effects of daily maximum temperature on the probability of leaf unfolding in temperate woody species along large elevational gradients. Although most of the phenological models have used daily mean temperature, here we used maximum temperature, because (i) daytime temperature seems to exert stronger effects on leaf unfolding than nighttime temperature (Fu et al., 2016) and (ii) we expect stronger lag effects of temperature beyond specific thresholds. The investigated tree and shrub species included early- to late-leafing species as well as photoperiod-sensitive and -insensitive species. We addressed the following research questions:

For how long do daily maximum temperatures of varying intensity affect leaf unfolding? We expect longer-term effects of cooler temperatures due to lower development rates of buds, but more immediate effects of warmer temperatures.

Do species differ with respect to the lagged associations between temperature and leaf unfolding? We expect stronger effects of temperature in late-leafing species, because leaf-out tends to occur at higher temperatures.

Does the lagged association between temperature and leaf unfolding change along elevation? We expect stronger effects of higher temperatures at lower elevations, as chilling requirements to break endodormancy may not always be fulfilled due to the warmer conditions.

 Materials and Methods Phenological Data

We used phenological observations of five tree and shrub species that were available from the phenological network of MeteoSwiss, the Swiss Federal Office of Meteorology and Climatology (Defila and Clot, 2001). We considered both broadleaved species (European beech, Fagus sylvatica L.; horse chestnut, Aesculus hippocastanum L.; common hazel, Corylus avellana L.) and conifer species (European larch, Larix decidua Mill.; Norway spruce, Picea abies (L.) H. Karst). The set of species ranges from early-leafing species with high frost resistance during leaf unfolding (larch, horse chestnut, hazel), to beech as an intermediate-leafing species with intermediate frost resistance and Norway spruce as a late-leafing species with low frost resistance (Lenz et al., 2013; Vitasse et al., 2014a; Bigler and Bugmann, 2018). While larch, horse chestnut and hazel rely mainly on temperature as a trigger of leaf unfolding, beech and Norway spruce rely on both temperature and photoperiod (Heide, 1993; Basler and Körner, 2012).

We included only phenological series without missing observations of leaf unfolding from 1972 to 2011 (i.e., 40 years of data; Supplementary Table S1 and Supplementary Figure S1). Phenological observations were conducted weekly by one voluntary observer per station, applying the same protocol for phenology monitoring across stations. The date of leaf unfolding was recorded when approximately 50% of the leaves of one or several trees or shrubs were unfolded, i.e., the leaf surface and leaf base is visible in broadleaved species, or the young needle bundles start to open and spread in conifer species. The dataset included 5280 observed dates of leaf unfolding from 42 phenological stations (beech: 26 stations; hazel: 23 stations; horse chestnut: 29 stations; larch: 35 stations; Norway spruce: 19 stations; Supplementary Table S1). The data had been checked for plausibility and consistency in a previous study (Bigler and Bugmann, 2018). All five species were observed jointly at five stations, four species at fourteen stations, three species at nine stations, two species at ten stations, and at four stations only one species was observed (Supplementary Table S1). The stations were distributed across Switzerland covering a perimeter of 12896 km2 and an elevational range from 200 to 1800 m a.s.l. (beech: 200 to 1240 m a.s.l.; hazel: 350 to 1120 m a.s.l.; horse chestnut: 200 to 1120 m a.s.l.; larch: 200 to 1800 m a.s.l.; Norway spruce: 200 to 1800 m a.s.l.; Supplementary Table S1). The climate regimes included mean annual temperature ranges from 2.9 to 12.4°C, and annual precipitation ranges from 751 to 1960 mm.

 Weather Data

We used spatially interpolated daily temperature data from 1 January 1972 to 31 December 2011 (i.e., 14610 values) at each of the 42 phenological stations. We selected daily maximum temperature Tmax as predictor variable (Supplementary Figure S1), because leaf unfolding responds more to daytime temperature than to nighttime temperature (Hanes, 2014; Piao et al., 2015; Fu et al., 2016). However, because daily mean temperature generally provides more accurate phenological predictions when using process-based models, we have also tested daily mean temperature (Tmean) as predictor variable to check the consistency of the resulting patterns. The daily temperature data were available for a grid of 100 m resolution across Switzerland, and were derived using the DAYMET interpolation algorithm (Thornton et al., 1997) applied to measured weather data from MeteoSwiss climate stations and a digital elevation model. The interpolated temperature data were provided by the Landscape Dynamics group at the Swiss Federal Research Institute WSL (Birmensdorf, Switzerland). For each phenological station, we averaged Tmax (Tmean) of the nearest grid cell and the surrounding eight neighboring cells to get robust estimates of Tmax (Tmean). Although the locations of the observed trees and shrubs may differ from the coordinates of the phenological stations by several hundreds of meters or in some cases even by several kilometers, we assumed that the interpolated temperature data still reflect relatively accurately the weather conditions of the trees and shrubs.

 Data Analysis

Rather than predicting leaf-out dates, we quantified the probability of leaf unfolding to occur under specific weather conditions. To quantify lagged temperature effects on leaf unfolding, we combined several statistical methods that are common in environmental epidemiology (Schwartz, 2000; Nitta et al., 2010; Gasparrini and Armstrong, 2013) but that have not been used in phenological studies (Hudson, 2010). We used the following approaches (see below for further details): (1) a case-crossover design to compare Tmax during case days (i.e., dates when leaf unfolding was observed) with control days (i.e., dates before or after leaf unfolding), accounting for lagged temperature effects up to 60 days before each case and control (Figure 1); (2) a conditional logistic regression to predict for each species and phenological station the effect of Tmax on the probability of leaf unfolding, which extends the logistic regression by accounting for stratified data (i.e., sets of cases and controls); and (3) multivariate meta-analyses to synthesize associations between Tmax and leaf unfolding across multiple phenological stations. The same approaches were used to analyze the data also based on Tmean, which provided qualitatively similar results as with Tmax as predictor variable (see Supplementary Figures S2, S3).

Visualization of the time-stratified case-crossover design. The example shows the observed leaf-out time of beech at the phenological station Liestal (350 m a.s.l.; Supplementary Table S1) in Switzerland and Tmax (maximum temperature) from July 2006 until end of June 2007. The red dot indicates the case (date of leaf unfolding: 17 April 2007), the black dots indicate the controls (dates before or after leaf unfolding: 3 April 2007, 10 April 2007, and 24 April 2007). The arrows indicate the 60-day lags that are considered in the distributed lag models. April 2007 is delimited by a white shaded box.

 Time-Stratified Case-Crossover Design

Case-crossover designs allow to compare environmental conditions during case days (i.e., dates when the phenological event occurred) that are matched to one or several control days (Maclure, 1991). This sampling design constitutes a self-matching of each phenological station and thus controls for time-invariant or only slowly time-varying confounders such as topography, vegetation etc. Several strategies have been proposed to match cases to controls. We selected the time-stratified case-crossover design (Lumley and Levy, 2000), i.e., each case (date of leaf unfolding) was matched to controls of the same year, month and weekday (Figure 1). Each set of 1 case and 3 to 4 controls is considered a stratum. The gaps of 6 days between case and control days (see Figure 1) reduce serial autocorrelation (Janes et al., 2005a). By restricting the controls to the same month as the case, only periods are considered when leaf unfolding of vegetative buds may actually be observed. The time-stratified case-crossover design further controls for potential bias due to seasonal and long-term changes in Tmax (Janes et al., 2005a), which ensures unbiased parameter estimates in conditional logistic regression (see below). Sampling of control days only before rare or non-recurrent events such as leaf unfolding would result in biased parameter estimates (Janes et al., 2005b). The dataset across all species and the entire 40-year period included 5280 cases and 17931 controls (beech: 1040 cases, 3563 controls; hazel: 920 cases, 3113 controls; horse chestnut: 1160 cases, 3921 controls; larch: 1400 cases, 4726 controls; Norway spruce: 760 cases, 2608 controls).

 Conditional Logistic Regression

Matched data from time-stratified case-crossover designs are commonly analyzed using conditional logistic regression (Breslow et al., 1978), which accounts for the non-independence of the stratified data. The probability of observing a case in stratum k (πk) may be defined as a stratum-specific logistic regression (Breslow et al., 1978; Hosmer et al., 2013):

πk(y=1 | x)=eαk+βx1+eαk+βx

where y is the dependent variable (y = 1: case, i.e., leaf unfolding occurs; y = 0: control, i.e., leaf unfolding does not occur), x is the predictor variable (e.g., Tmax or a function of Tmax), αk is a stratum-specific constant, and β is the coefficient for variable x, which is estimated across all strata. In this study, the strata k correspond to the years during the observation period (i.e., k = 1972, …, 2011). Equation 1 may be rewritten as log-odds:

logπk(y=1 | x)1−πk(y=1 | x)=αk+βx

Thus, a unit increase in variable x changes the log-odds by β. Equation 2 contains no intercept, but a stratum-specific constant αk, which cancels out in the conditional likelihood.

The effect of x on the probability of observing a case is typically expressed as odds ratio, i.e., the ratio between the odds of a value x1 and the odds of a reference value x0:

odds ratio= πk(y=1 | x1)/(1−  πk(y=1 | x1))πk(y=1 | x0)/(1−  πk(y=1 | x0))= eβ(x1−x0)

We used the statistical computing software R (R Core Team, 2018), version 3.5.1, to fit conditional logistic regression models (function “clogit” in the package “survival,” version 2.42-4). We selected an exact calculation method to maximize the conditional likelihood.

 Distributed Lag Models

We expanded the conditional logistic regression (eq. 2) to account for lagged effects of Tmax on leaf unfolding. In particular, we fitted distributed lag non-linear models (DLNMs) to the data, which allow for flexible and delayed associations between a response variable and a predictor variable along time (Gasparrini et al., 2010). The DLNM was defined as follows:

logπk(y=1 | xt, …, xt−l, …, xt−L)1−πk(y=1 | xt, …, xt−l, …, xt−L)=αk+s(xt, …, xt−l, …, xt−L;  η)

where the smoothing function s(…) was represented by a bi-dimensional natural spline for the predictor variable xt (here Tmax) at day t and lag days l (l = 0, …, L) and a vector of coefficients η (Gasparrini et al., 2010). We used natural cubic spline functions with two interior knots for both the predictor variable Tmax and the lag dimension. Tmax was restricted to the interval [0°C, 25°C] to prevent artifacts that may have occurred at particularly warm or cool stations due to extrapolations of temperature beyond the range of observed values. The maximum lag L was set at 60 days (Figure 1) based on findings in previous studies (Fu et al., 2015; Vitasse et al., 2018). The interior knots were selected to evenly split the range of Tmax (knots at 8.3 and 16.7°C) and lag dimension (knots at 20 and 40 days) into three intervals. We calculated predictor-specific summaries at Tmax of 5, 10, 15, and 20°C. Odds ratios (eq. 3) were calculated with respect to the reference value 0°C. Thus, we predicted odds ratios for lags from 0 to 60 days with Tmax set to 5, 10, 15 and 20°C. We refer to Figure 2 for an illustrative example. Because odds ratios not only depend on the estimated model coefficients but also on the difference between Tmax and the reference value of 0°C (eqs. 3 and 4), comparisons of the odds ratios are only valid across species, lags and elevations, but not across different values of Tmax. However, significance of the odds ratios may be compared at different values of Tmax. We used the R packages “dlnm” (Gasparrini, 2011), version 2.3.4, and “splines,” version 3.5.1, to calculate the DLNMs.

Visualization of a DLNM (distributed lag non-linear model). The example shows the odds ratio of beech at the phenological station Liestal (350 m a.s.l.; Supplementary Table S1) in Switzerland along Tmax (maximum temperature) and lag dimension. A conditional logistic regression model was fitted to time-stratified case-crossover data (Figure 1) from 1972 to 2011. The odds ratios are interpreted as the ratio between the odds of a specific Tmax at a specific lag compared to the odds of the reference Tmax (0°C) at the same lag (eq. 3). The red line indicates the odds ratio along Tmax at a lag of 0 days, the orange lines indicate the odds ratios along the lag dimension at Tmax of 0°C (odds ratio = 1), 5°C, 10°C, 15°C, 20°C, and 25°C.

 Multivariate Meta-Analyses

A multivariate meta-analysis combines estimates of several parameters (e.g., coefficients from a regression model) across models from several sites, which allows to estimate average parameters. We used multivariate meta-analyses to synthesize lagged effects of Tmax on leaf unfolding across multiple phenological stations (Jackson et al., 2011; Gasparrini et al., 2012). Multivariate meta-analysis allows to (i) pool estimates of multi-parameter associations and quantify between-station variability, and (ii) consider meta-variables that vary across stations. Based on the estimated coefficients ^ η i (eq. 4) and the corresponding variance-covariance matrix ^ S i from each phenological station i, random-effects multivariate meta-analyses were conducted (Gasparrini and Armstrong, 2013). We further included elevation as a meta-variable in random-effects multivariate meta-regression models (Gasparrini and Armstrong, 2013). The models were fitted to the data using restricted maximum likelihood from the R package “mvmeta,” version 0.4.11 (Gasparrini et al., 2012).

 Results

Leaf unfolding was observed when relatively high daily maximum temperatures (Tmax) prevailed, which increased from early- to late-leafing species (Figure 3 and Supplementary Figure S1): larch, 14.6 ± 4.7°C (mean ± standard deviation); horse chestnut, 15.0 ± 4.9°C; hazel, 15.1 ± 5.0°C; beech, 16.3 ± 4.9°C; Norway spruce, 17.0 ± 4.6°C (species ordered from early- to late-leafing; Bigler and Bugmann, 2018). However, in some cases, leaf unfolding was recorded at low temperatures with Tmax ≤ 5°C (Figure 3): larch, 2.43% of all observations; horse chestnut, 2.50%; hazel, 2.50%; beech, 1.63%; Norway spruce, 0.26%. The expected Tmax during the day of observed leaf unfolding decreased with increasing elevation for all species (Figure 3): larch, -0.22°C/100 m; horse chestnut, -0.41°C/100 m; hazel, -0.27°C/100 m; beech, -0.38°C/100 m; Norway spruce, -0.25°C/100 m. Because the distribution of daily Tmax during the day of leaf unfolding does not provide a proper account of temperature effects on leaf unfolding, we need to investigate the temperature profiles before leaf unfolding, which will be considered in the following analyses.

Change of Tmax (maximum temperature) during the day of leaf unfolding with increasing elevation for larch, horse chestnut, hazel, beech, and Norway spruce. The species are ordered from early-leafing (top) to late-leafing (bottom). The blue lines indicate the regression lines based on the linear models.

For each species and station, we fitted a conditional logistic regression model with distributed lags (eq. 4) to the time-stratified data over the entire 40-year period (see Figure 2 for an example). The species-specific summaries of the odds ratios (eq. 3) for specific values of Tmax (5, 10, 15, and 20°C) along the lag dimension are shown in Figure 4. At relatively short lags (0 to ca. 10 days), the odds ratios of most phenological stations were larger than 1, which tended to increase with decreasing elevation for some species (e.g., beech). The odds ratios of many phenological stations decreased with increasing lag (Figure 4), thus, the short-term (0 to ca. 10 days) effects of Tmax ≥ 5°C on the probability of leaf unfolding were higher compared to the reference Tmax of 0°C than the longer-term effects. At longer lags (ca. 50 to 60 days), the odds ratios were distributed around 1 when Tmax was 5 or 10°C. Particularly for warmer temperatures (Tmax of 15 and 20°C), the odds ratios dropped below 1 at longer lags for quite many stations, i.e., the probability of observing leaf unfolding 50 to 60 days following a warm day was lower than at 0°C. Such effects may occur, when leaf unfolding has already occurred (Figure 1 and Supplementary Figure S1). Overall, there was a relatively large variability among stations regarding the lagged association between Tmax and leaf unfolding (Figure 4).

Summaries of DLNMs (distributed lag non-linear models) based on conditional logistic regression for larch, horse chestnut, hazel, beech, and Norway spruce. The odds ratios along the lag dimension are shown for Tmax of 5°C, 10°C, 15°C, and 20°C. The species are ordered from early- to late-leafing species (left to right). The station-specific odds ratios are represented by the dark-blue lines (low-elevation stations) to light-blue lines (high-elevation stations). For sake of clarity, no confidence intervals are shown for the station-specific odds ratios. The red lines are the pooled estimates of the odds ratios (including 95% confidence intervals) from the multivariate meta-analysis. An inverse hyperbolic sine transformation has been applied to the y-axis to increase the visibility of smaller values.

Based on the multivariate meta-analyses, the pooled estimates and 95% confidence intervals of the odds ratios provide a general overview on the change of the odds ratios and their significance with increasing lags (Figure 4). The early-leafing species (larch, horse chestnut, hazel) showed quite a similar pattern, i.e., short- (0 to ca. 10 days) to long-term effects (ca. 50 to 60 days) with significant odds ratios >1 were detected at Tmax of 5 and 10°C, whereas short-term effects dominated at 20°C. For the later-leafing species beech and Norway spruce, the pooled odds ratios at short lags (0 to ca. 10 days) were larger than for the early-leafing species, but also suggest a higher variability among stations as indicated by the relatively wide confidence intervals (Figure 4). Larger odds ratios translate into more rapid leaf unfolding under identical temperatures. For example, the pooled odds ratio at Tmax of 15°C and at a lag of 0 days (Figure 4) increased from hazel (1.6), larch (1.7), horse chestnut (2.1), beech (3.8) to Norway spruce (4.2), i.e., leaf unfolding of beech and Norway spruce occurs more rapidly than for hazel, larch and horse chestnut at this temperature. The pooled estimates based on the multivariate meta-analyses using Tmean resulted in qualitatively similar results (Supplementary Figure S2), though the odds ratios at short lags (0 to ca. 10 days) were a bit lower than for Tmax.

The pooled odds ratios based on the multivariate meta-regression models showed distinct and significant effects along elevation (Figure 5). Particularly the broadleaved species (horse chestnut, hazel, beech) showed consistent elevational effects within ca. 10 to 60 days, i.e., odds ratios generally decreased from low elevation (300 m a.s.l.) to mid elevation (700 m a.s.l.) and high elevation (1100 m a.s.l.), which translates into more rapid leaf unfolding at lower elevations under identical temperatures. For example, the pooled odds ratio for beech at Tmax of 15°C and at a lag of 0 days is 5.5 at 300 m a.s.l. and 2.8 at 1100 m a.s.l. (Figure 5), i.e., leaf unfolding of beech occurs more rapidly at lower elevations than at higher elevations at this temperature. The pooled estimates based on the multivariate meta-regressions using Tmean also resulted in similar results (Supplementary Figure S3), but with slightly lower odds ratios at short lags and generally wider 95% confidence intervals compared to Tmax.

Summaries of DLNMs (distributed lag non-linear models) based on conditional logistic regression for larch, horse chestnut, hazel, beech, and Norway spruce. The odds ratios along the lag dimension are shown for Tmax of 5°C, 10°C, 15°C, and 20°C. The species are ordered from early- to late-leafing species (left to right). The odds ratios for three elevations (300 m, 700 m, 1100 m a.s.l.) and 95% confidence intervals were estimated from the multivariate meta-regression. An inverse hyperbolic sine transformation has been applied to the y-axis to increase the visibility of smaller values.

 Discussion

Long-term series of observed leaf-out timings from conifers and broadleaved species provide evidence of species- and elevation-specific responses to temperature. The time-stratified case-crossover design combined with the statistical framework proved to be a suitable approach to (i) describe simultaneously the non-linear and delayed dependencies of leaf unfolding on daily temperature, and (ii) synthesize lagged associations between temperature and leaf unfolding across multiple phenological stations. We detected a decline in the lagged effects of varying maximum temperature (Tmax) on leaf unfolding. Thus, the strongest effects of Tmax on the development of vegetative buds are expected within a few days, while lagged effects after ca. 30 days tend to fade out. Beech and Norway spruce, which are both sensitive to photoperiod and tend to leaf-out late, showed stronger short-term effects of Tmax on leaf unfolding than larch, horse chestnut and hazel. Broadleaved species (horse chestnut, hazel, beech) showed a stronger short-term response to Tmax at low elevations than at higher elevations.

Lagged Effects of Temperature on Leaf Unfolding

We established lagged and non-linear effects of Tmax on the probability of leaf unfolding. Temperature effects decrease with increasing lags as demonstrated by the pooled odds ratios of the multivariate meta-analyses (Figure 4). Thus, immediate temperature effects within ca. 10 days are more likely to induce leaf unfolding than lagged effects after ca. 30 days or more. Lower Tmax (5 and 10°C) tends to stimulate leaf unfolding even at longer lags (i.e., after more than ca. 30 days). In contrast at Tmax of 20°C this is the case only at relatively short lags and after more than ca. 10 days even turns into a lower odds compared to the reference Tmax of 0°C. The latter effect is known as “harvesting” (Gasparrini et al., 2010). In other words, odds ratios decrease for events such as leaf-out dates after the event has occurred, because they occur only once per year.

Despite these relatively distinct effects of pooled odds ratios, we observed quite a high variability of lagged associations between Tmax and leaf unfolding among phenological stations (Figure 4). Several reasons are likely to account for this heterogeneity: (1) Environmental conditions at the phenological stations do not only differ in terms of Tmax, but also with respect to further influences such as precipitation, radiation, topography or species composition, which may affect microclimate and site conditions as perceived by the trees. (2) The observed trees of the phenological network differ in terms of tree height, diameter, age, genetic status, social position within the forest stand, and the surrounding forest structure (e.g., forest edge vs. interior stand). For example ontogenetic differences between adult trees and seedlings may substantially affect the timing of leaf unfolding in broadleaved trees (Vitasse, 2013). (3) The true positions of the observed trees may differ from the coordinates of the phenological stations, i.e., the interpolated weather data do not fully represent the conditions experienced by the trees. (4) There might be an observer bias, because phenological observers may slightly differ in their assessment of the timing of leaf unfolding (Bison et al., 2018; Güsewell et al., 2018). Overall, however, there seems to be a pattern of lagged and non-linear associations between Tmax and leaf unfolding that is common across stations.

The preseason, which is related to the concept of lagged temperature effects on leaf unfolding, is defined as the period before leaf-out time when specific weather conditions are associated with leaf unfolding (Güsewell et al., 2017). The length of the preseason has been estimated in several studies during which temperatures exert a statistically significant influence on the timing of leaf unfolding. Correlation coefficients and partial correlation coefficients, respectively, between mean dates of leaf unfolding and temperatures averaged over varying periods before leaf-out have been used to determine the length of the preseason. Fu et al. (2015) found for seven broadleaved tree species across Europe, including horse chestnut and beech among other species, optimal preseason lengths of 45 to 60 days based on daily Tmean. However, there was substantial within-species variability across sites ranging from 15 days to 4 months. The preseason for 24 plant species across Europe was estimated based on monthly Tmax to range between 0 and 3 months for most combinations of species, sites and years (Piao et al., 2015). Using data from the phenological network of MeteoSwiss, Güsewell et al. (2017) estimated the optimal preseason to range between less than 60 days for beech and up to more than 70 days for Norway spruce, while Vitasse et al. (2018) found the highest correlations for the preseason to range from 40 to 50 days (beech) to ca. 65 days (Norway spruce). Thus, these estimates of preseason length are generally consistent with the findings in our study, even if methodologically quite different approaches were used. However, correlations between mean leaf-out time and mean temperature do not account for (i) variability of leaf-out time across sites or years; (ii) daily resolution of temperature data; and (iii) the time-ordered structure of the predictors (e.g., sequence of daily temperature).

Lagged temperature effects on leaf unfolding such as based on distributed lag models in our study or based on partial correlation coefficients (Fu et al., 2015; Piao et al., 2015; Güsewell et al., 2017) may be interpreted as indirect evidence of instantaneous thermal effects on the continuous development, growth or physiological activity of vegetative buds during dormancy and bud burst (Clark et al., 2014b). Increasing temperatures in spring induce deacclimation processes, whereby frost resistance decreases for example due to decreased levels of soluble carbohydrates, that finally result in bud burst (Pagter and Arora, 2013). Deacclimation occurs with some lag following exposure to warm temperatures, typically within a few days to weeks (Pagter and Arora, 2013). While it would be desirable to take continuous and more process-related measurements during dormancy and bud burst, observational and experimental studies are restricted in this regard. Assessing the state of dormancy in vegetative buds tends to be difficult, because molecular or physiological markers are lacking (Cooke et al., 2012) and because changes of bud morphology are hardly visible prior to bud break (but see Basler and Körner, 2014).

 Associations Between Temperature and Leaf Unfolding Differ Across Species

The lagged associations between temperature and leaf unfolding revealed differences between early-leafing, photoperiod-insensitive species (larch, horse chestnut, hazel) and late-leafing, photoperiod-sensitive species (beech, Norway spruce). Both beech and Norway spruce showed stronger short- to mid-term responses (i.e., within ca. 10 to 20 days) to any of the specified temperatures (Figure 4). Thus, leaf unfolding of beech and Norway spruce is expected to occur more rapidly than for larch, horse chestnut and hazel when periods with the same temperatures prevail. These patterns may reflect different leafing strategies to reduce frost risk following unseasonably warm spells, i.e., either by slow deacclimation in the early-leafing species or by more rapid deacclimation in the late-leafing species (Heide, 1974; Pagter and Arora, 2013).

Beech and Norway spruce are sensitive to photoperiod (Heide, 1974; Basler and Körner, 2012, 2014; Vitasse and Basler, 2013), and both species tend to leaf out rather late because of the lower frost resistance at leaf-out compared to the other species (Repo, 1992; Taschler et al., 2004; Lenz et al., 2013; Vitasse et al., 2014a; Bigler and Bugmann, 2018). Leaf unfolding of beech lags the early-leafing species in this study by ca. 1 to 2 weeks, Norway spruce by ca. 2 to 3 weeks or more (Bigler and Bugmann, 2018). Sensitivity to photoperiod or high chilling requirement to break dormancy may be considered as an additional safeguard allowing frost sensitive species to further minimize the risk of frost damage (Körner and Basler, 2010; Way and Montgomery, 2015), e.g., if a warm early spring or late winter is followed by a series of severe frost events. Species-specific rates of deacclimation in spring result in different development rates of vegetative buds (Pagter and Arora, 2013). Thus, when the photoperiodic and chilling requirements are fulfilled, the vegetative buds of the trees respond more rapidly to increasing temperatures (Körner and Basler, 2010).

The photoperiod-insensitive species and species having low chilling requirements (larch, horse chestnut, hazel), respectively, tend to leaf out rather early (Bigler and Bugmann, 2018). These species mainly rely on forcing temperatures as a trigger of bud burst and leaf unfolding (Basler and Körner, 2012). To compensate for the higher frost risk, they maintain a higher frost resistance during leaf-out. The weaker leaf phenological response to warmth compared to the late-leafing photoperiod-sensitive species may be interpreted as a further strategy to reduce the risk of frost damages (Lenz et al., 2016), because an overly rapid bud development may result in early leaf unfolding occurring simultaneously with or shortly after severe frost. Because the weather conditions may strongly change during early spring, and because there is a particularly high uncertainty regarding the occurrence of frost, plants have to take precautionary measures to reduce the risk of frost damages following early warm spells (Pagter and Arora, 2013). Advancing the timing of leaf unfolding by just 1–3 weeks strongly increases the risk of frost damages due to the shorter return intervals of frost events during the leaf-out period (Lenz et al., 2016).

Differences among tree species in terms of leaf-out times, frost resistance at the time of leaf-out and lagged temperature effects on leaf unfolding reflect varying strategies to cope with uncertain environmental conditions. Particularly late frost events in spring are considered ecologically and evolutionarily important events that may have severe consequences for the dynamics of tree populations, because physical and physiological damages occur that affect parts of trees or entire tree individuals (Inouye, 2000; Gu et al., 2008). Thus, frost is a strong selective pressure (Lenz et al., 2013) and is assumed to ultimately act as a key factor that determines tree species ranges (Chuine and Beaubien, 2001; Chuine, 2010; Körner et al., 2016).

 Effects of Elevation on the Association Between Temperature and Leaf Unfolding

The lagged association between temperature and leaf unfolding changes along elevation, particularly for broadleaved species (horse chestnut, hazel, beech), less so for conifers (larch, Norway spruce; Figure 5). The temperature effects increased for the broadleaved species with decreasing elevation, particularly within the first ca. 10 to 40 days, depending on species and specified temperature. Thus, under identical temperature conditions, leaf unfolding of these broadleaved species is expected to occur more rapidly at low elevation than at high elevation, assuming that chilling and photoperiod requirements are fulfilled. Some of these elevational effects may be explained with the decreasing Tmax with increasing elevation (Figure 3 and Supplementary Figure S1), i.e., trees at higher elevations tend to leaf out at lower temperatures than trees at lower elevations.

The stronger response of horse chestnut, hazel and beech to temperature at lower elevations may indicate that a lack of chilling could be compensated by increased forcing temperatures. Beech as a photoperiod-sensitive species is further known to be more limited by daylength and chilling at lower elevations than at higher elevations (Vitasse and Basler, 2013). In this context, our findings may be interpreted that leaf unfolding of low-elevation beech trees responds strongly and relatively rapidly as soon as the photoperiodic and chilling requirements are fulfilled.

Leaf unfolding is delayed with increasing elevation as a consequence of the cooler temperatures (Lenz et al., 2016; Bigler and Bugmann, 2018; Vitasse et al., 2018). While temperatures generally decrease with increasing elevation, temperature variability tends to increase, thus uncertainty regarding the occurrence of frost is relatively high. Return intervals of frost events during the leaf-out period typically decrease with increasing elevation (Lenz et al., 2016), therefore, even relatively warm temperatures in spring or early summer at high elevation do not guarantee absence of frost. We interpret the lower odds ratios of the lagged association between temperature and leaf unfolding at higher elevations in our study as an adaption to reduce frost risk.

 Conclusion

The findings from our study provide evidence of lagged, non-linear temperature effects on leaf unfolding of temperate broadleaved species and conifers with contrasting leaf-out timing, frost resistance and photoperiodic requirements. Maximum temperatures induce relatively strong effects on leaf unfolding within a few days, while lagged effects over more than 1 month are weaker. Our findings demonstrate cumulative and long-term temperature effects on leaf unfolding, whereby relatively high temperatures before leaf unfolding accelerate bud development inducing rapid bud burst. We assume that such effects are not observed earlier in spring, because the buds are in a premature developmental state. Tree species respond differently to the impact of temperature: the immediate reaction of late-leafing, photoperiod-sensitive species is stronger, i.e., leaf unfolding occurs more rapidly than for early-leafing, photoperiod-insensitive species. Broadleaved species further show stronger temperature effects at low elevation than at high elevation, which translates into more rapid leaf unfolding with decreasing elevation under identical temperature conditions. The lagged associations between temperature and leaf unfolding improve our understanding on how the phenology of tree species with differing ecological requirements that occur along large elevational gradients respond to the cyclic, seasonal course of weather.

 Data Availability

Publicly available datasets were analyzed in this study. This data can be found here: https://gate.meteoswiss.ch/idaweb.

 Author Contributions

CB conceived and designed the study, processed and analyzed the data, and wrote the draft manuscript. CB and YV discussed the results. YV contributed to the writing of the manuscript. Both authors contributed to manuscript revision, read and approved the submitted version.

 Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding. YV was partly funded by the SwissForestLab.

We acknowledge discussions with Regula Gehrig and Barbara Pietragalla (MeteoSwiss, Swiss Federal Office of Meteorology and Climatology), the many voluntary observers of the Swiss phenology network for recording the phenological data, and the Swiss Federal Institute for Forest, Snow and Landscape Research WSL for providing the interpolated DAYMET climate data. We would like to thank Antonio Gasparrini for developing and maintaining the R packages “dlnm” and “mvmeta.”

 Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2019.00398/full#supplementary-material

 References Basler D. (2016). Evaluating phenological models for the prediction of leaf-out dates in six temperate tree species across central Europe. Agric. For. Meteorol. 217 10–21. 10.1016/j.agrformet.2015.11.007 Basler D. Körner C. (2012). Photoperiod sensitivity of bud burst in 14 temperate forest tree species. Agric. For. Meteorol. 165 73–81. 10.1016/j.agrformet.2012.06.001 Basler D. Körner C. (2014). Photoperiod and temperature responses of bud swelling and bud burst in four temperate forest tree species. Tree Physiol. 34 377–388. 10.1093/treephys/tpu021 24713858 Bigler C. Bugmann H. (2018). Climate-induced shifts in leaf unfolding and frost risk of European trees and shrubs. Sci. Rep. 8:9865. 10.1038/s41598-018-27893-1 29959342 Bison M. Yoccoz N. G. Carlson B. Z. Delestrade A. (2018). Comparison of budburst phenology trends and precision among participants in a citizen science program. Int. J. Biometeorol. 63 61–72. 10.1007/s00484-018-1636-x 30382351 Breslow N. E. Day N. E. Halvorsen K. T. Prentice R. L. Sabai C. (1978). Estimation of multiple relative risk functions in matched case-control studies. Am. J. Epidemiol. 108 299–307. 10.1093/oxfordjournals.aje.a112623 Caffarra A. Donnelly A. Chuine I. (2011a). Modelling the timing of Betula pubescens budburst. II. Integrating complex effects of photoperiod into process-based models. Clim. Res. 46 159–170. 10.3354/cr00983 Caffarra A. Donnelly A. Chuine I. Jones M. B. (2011b). Modelling the timing of Betula pubescens budburst. I. Temperature and photoperiod: a conceptual model. Clim. Res. 46 147–157. 10.3354/cr00980 Cannell M. G. R. Smith R. I. (1983). Thermal time, chill days and prediction of budburst in Picea sitchensis. J. Appl. Ecol. 20 951–963. 10.2307/2403139 Chuine I. (2000). A unified model for budburst of trees. J. Theor. Biol. 207 337–347. 10.1006/jtbi.2000.2178 11082304 Chuine I. (2010). Why does phenology drive species distribution? Philos. Trans. R. Soc. B-Biol. Sci. 365 3149–3160. 10.1098/rstb.2010.0142 20819809 Chuine I. Beaubien E. G. (2001). Phenology is a major determinant of tree species range. Ecol. Lett. 4 500–510. 10.1046/j.1461-0248.2001.00261.x Chuine I. Garcia de Cortazar Atauri I. Kramer K. Hänninen H. (2013). “Plant Development Models,” in Phenology: An Integrative Environmental Science ed. Schwartz M. D. (Dordrecht: Springer) 275–293. 10.1007/978-94-007-6925-0_15 Chuine I. Morin X. Bugmann H. (2010). Warming, photoperiods, and tree phenology. Science 329 277–278. 10.1126/science.329.5989.277-e 20647448 Clark J. S. Melillo J. Mohan J. Salk C. (2014a). The seasonal timing of warming that controls onset of the growing season. Glob. Chang. Biol. 20 1136–1145. 10.1111/gcb.12420 24115181 Clark J. S. Salk C. Melillo J. Mohan J. (2014b). Tree phenology responses to winter chilling, spring warming, at north and south range limits. Funct. Ecol. 28 1344–1355. 10.1111/1365-2435.12309 Cline M. G. Deppong D. O. (1999). The role of apical dominance in paradormancy of temperate woody plants: a reappraisal. J. Plant Physiol. 155 350–356. 10.1016/S0176-1617(99)80116-3 Cooke J. E. K. Eriksson M. E. Junttila O. (2012). The dynamic nature of bud dormancy in trees: environmental control and molecular mechanisms. Plant Cell Environ. 35 1707–1728. 10.1111/j.1365-3040.2012.02552.x 22670814 Defila C. Clot B. (2001). Phytophenological trends in Switzerland. Int. J. Biometeorol. 45 203–207. 10.1007/s004840100101 11769321 Delpierre N. Vitasse Y. Chuine I. Guillemot J. Bazot S. Rutishauser T. (2016). Temperate and boreal forest tree phenology: from organ-scale processes to terrestrial ecosystem models. Ann. For. Sci. 73 5–25. 10.1007/s13595-015-0477-6 Forrest J. Miller-Rushing A. J. (2010). Toward a synthetic understanding of the role of phenology in ecology and evolution. Philos. Trans. R. Soc. B-Biol. Sci. 365 3101–3112. 10.1098/rstb.2010.0145 20819806 Fu Y. S. H. Liu Y. J. De Boeck H. J. Menzel A. Nijs I. Peaucelle M. (2016). Three times greater weight of daytime than of night-time temperature on leaf unfolding phenology in temperate trees. New Phytol. 212 590–597. 10.1111/nph.14073 27376563 Fu Y. S. H. Zhao H. F. Piao S. L. Peaucelle M. Peng S. S. Zhou G. Y. (2015). Declining global warming effects on the phenology of spring leaf unfolding. Nature 526 104–107. 10.1038/nature15402 26416746 Gasparrini A. (2011). Distributed lag linear and non-linear models in R: the package dlnm. J. Stat. Softw. 43 1–20. 10.18637/jss.v043.i08 22003319 Gasparrini A. Armstrong B. (2013). Reducing and meta-analysing estimates from distributed lag non-linear models. BMC Med. Res. Methodol. 13:1. 10.1186/1471-2288-13-1 23297754 Gasparrini A. Armstrong B. Kenward M. G. (2010). Distributed lag non-linear models. Stat. Med. 29 2224–2234. 10.1002/sim.3940 20812303 Gasparrini A. Armstrong B. Kenward M. G. (2012). Multivariate meta-analysis for non-linear and other multi-parameter associations. Stat. Med. 31 3821–3839. 10.1002/sim.5471 22807043 Gu L. Hanson P. J. Mac Post W. Kaiser D. P. Yang B. Nemani R. (2008). The 2007 eastern US spring freeze: Increased cold damage in a warming world? BioScience 58 253–262. 10.1641/B580311 Güsewell S. Furrer R. Gehrig R. Pietragalla B. (2017). Changes in temperature sensitivity of spring phenology with recent climate warming in Switzerland are related to shifts of the preseason. Glob. Chang. Biol. 23 5189–5202. 10.1111/gcb.13781 28586135 Güsewell S. Pietragalla B. Gehrig R. Furrer R. (2018). Representativeness of stations and reliability of data in the Swiss phenology network. Tech. Rep. MeteoSwiss 267:100. Hanes J. M. (2014). Spring leaf phenology and the diurnal temperature range in a temperate maple forest. Int. J. Biometeorol. 58 103–108. 10.1007/s00484-012-0603-1 23150087 Hänninen H. (2016). Boreal and Temperate Trees in a Changing Climate: Modelling the Ecophysiology of Seasonality. Dordrecht: Springer. 10.1007/978-94-017-7549-6 Heide O. M. (1974). Growth and dormancy in Norway spruce ecotypes (Picea abies): I. Interaction of photoperiod and temperature. Physiol. Plant. 30 1–12. 10.1111/j.1399-3054.1974.tb04983.x Heide O. M. (1993). Dormancy release in beech buds (Fagus sylvatica) requires both chilling and long days. Physiol. Plant. 89 187–191. 10.1111/j.1399-3054.1993.tb01804.x Horvath D. P. Anderson J. V. Chao W. S. Foley M. E. (2003). Knowing when to grow: signals regulating bud dormancy. Trends Plant Sci. 8 534–540. 10.1016/j.tplants.2003.09.013 14607098 Hosmer D. W. Lemeshow S. Sturdivant R. X. (2013). Applied Logistic Regression. Hoboken, NJ: John Wiley & Sons, Inc. 10.1002/9781118548387 Hudson I. L. (2010). Interdisciplinary approaches: towards new statistical methods for phenological studies. Clim. Chang. 100 143–171. 10.1007/s10584-010-9859-9 Ibáñez I. Primack R. B. Miller-Rushing A. J. Ellwood E. Higuchi H. Lee S. D. (2010). Forecasting phenology under global warming. Philos. Trans. R. Soc. B-Biol. Sci. 365 3247–3260. 10.1098/rstb.2010.0120 20819816 Inouye D. W. (2000). The ecological and evolutionary significance of frost in the context of climate change. Ecol. Lett. 3 457–463. 10.1046/j.1461-0248.2000.00165.x Jackson D. Riley R. White I. R. (2011). Multivariate meta-analysis: potential and promise. Stat. Med. 30 2481–2498. 10.1002/sim.4172 21268052 Janes H. Sheppard L. Lumley T. (2005a). Case-crossover analyses of air pollution exposure data: referent selection strategies and their implications for bias. Epidemiology 16 717–726. 16222160 Janes H. Sheppard L. Lumley T. (2005b). Overlap bias in the case-crossover design, with application to air pollution exposures. Stat. Med. 24 285–300. 10.1002/sim.1889 15546133 Körner C. Basler D. (2010). Phenology under global warming. Science 327 1461–1462. 10.1126/science.1186473 20299580 Körner C. Basler D. Hoch G. Kollas C. Lenz A. Randin C. F. (2016). Where, why and how? Explaining the low-temperature range limits of temperate tree species. J. Ecol. 104 1076–1088. 10.1111/1365-2745.12574 Lang G. A. Early J. D. Martin G. C. Darnell R. L. (1987). Endo-, para-, and ecodormancy: physiological terminology and classification for dormancy research. Hortscience 22 371–377. Laube J. Sparks T. H. Estrella N. Höfler J. Ankerst D. P. Menzel A. (2014). Chilling outweighs photoperiod in preventing precocious spring development. Glob. Chang. Biol. 20 170–182. 10.1111/gcb.12360 24323535 Lenz A. Hoch G. Körner C. Vitasse Y. (2016). Convergence of leaf-out towards minimum risk of freezing damage in temperate trees. Funct. Ecol. 30 1480–1490. 10.1111/1365-2435.12623 Lenz A. Hoch G. Vitasse Y. Körner C. (2013). European deciduous trees exhibit similar safety margins against damage by spring freeze events along elevational gradients. New Phytol. 200 1166–1175. 10.1111/nph.12452 23952607 Linkosalo T. Lappalainen H. K. Hari P. (2008). A comparison of phenological models of leaf bud burst and flowering of boreal trees using independent observations. Tree Physiol. 28 1873–1882. 10.1093/treephys/28.12.1873 19193570 Luedeling E. Gassner A. (2012). Partial least squares regression for analyzing walnut phenology in California. Agric. For. Meteorol. 158 43–52. 10.1016/j.agrformet.2011.10.020 Luedeling E. Kunz A. Blanke M. M. (2013). Identification of chilling and heat requirements of cherry trees-a statistical approach. Int. J. Biometeorol. 57 679–689. 10.1007/s00484-012-0594-y 23053065 Lumley T. Levy D. (2000). Bias in the case-crossover design: implications for studies of air pollution. Environmetrics 11 689–704. 10.1002/1099-095X(200011/12)11:6<689::AID-ENV439>3.0.CO;2-N Maclure M. (1991). The case-crossover design: a method for studying transient effects on the risk of acute events. Am. J. Epidemiol. 133 144–153. 10.1093/oxfordjournals.aje.a115853 Menzel A. Sparks T. H. Estrella N. Koch E. Aasa A. Ahas R. (2006). European phenological response to climate change matches the warming pattern. Glob. Chang. Biol. 12 1969–1976. 10.1111/j.1365-2486.2006.01193.x Murray M. B. Cannell M. G. R. Smith R. I. (1989). Date of budburst of fifteen tree species in Britain following climatic warming. J. Appl. Ecol. 26 693–700. 10.2307/2404093 Nitta H. Yamazaki S. Omori T. Sato T. (2010). An introduction to epidemiologic and statistical methods useful in environmental epidemiology. J. Epidemiol. 20 177–184. 10.2188/jea.JE20100010 Olsson C. Jönsson A. M. (2014). Process-based models not always better than empirical models for simulating budburst of Norway spruce and birch in Europe. Glob. Chang. Biol. 20 3492–3507. 10.1111/gcb.12593 24700779 Pagter M. Arora R. (2013). Winter survival and deacclimation of perennials under warming climate: physiological perspectives. Physiol. Plant. 147 75–87. 10.1111/j.1399-3054.2012.01650.x 22583023 Piao S. L. Tan J. G. Chen A. P. Fu Y. H. Ciais P. Liu Q. (2015). Leaf onset in the northern hemisphere triggered by daytime temperature. Nat. Commun. 6:6911. 10.1038/ncomms7911 25903224 Polgar C. A. Primack R. B. (2011). Leaf-out phenology of temperate woody plants: from trees to ecosystems. New Phytol. 191 926–941. 10.1111/j.1469-8137.2011.03803.x 21762163 R Core Team (2018). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Repo T. (1992). Seasonal changes of frost hardiness in Picea abies and Pinus sylvestris in Finland. Can. J. For. Res. 22 1949–1957. Roberts A. M. I. (2008). Exploring relationships between phenological and weather data using smoothing. Int. J. Biometeorol. 52 463–470. 10.1007/s00484-007-0141-4 18193297 Roberts A. M. I. (2012). Comparison of regression methods for phenology. Int. J. Biometeorol. 56 707–717. 10.1007/s00484-011-0472-z 21789728 Roberts A. M. I. Tansey C. Smithers R. J. Phillimore A. B. (2015). Predicting a change in the order of spring phenology in temperate forests. Glob. Chang. Biol. 21 2603–2611. 10.1111/gcb.12896 25731862 Schwartz J. (2000). The distributed lag between air pollution and daily deaths. Epidemiology 11 320–326. 10.1097/00001648-200005000-00016 Schwartz M. D. (ed.) (2013). Phenology: an Integrative Environmental Science. Dordrecht: Springer. 10.1007/978-94-007-6925-0 Sparks T. H. Carey P. D. (1995). The responses of species to climate over two centuries: an analysis of the Marsham phenological record, 1736-1947. J. Ecol. 83 321–329. 10.2307/2261570 Taschler D. Beikircher B. Neuner G. (2004). Frost resistance and ice nucleation in leaves of five woody timberline species measured in situ during shoot expansion. Tree Physiol. 24 331–337. 10.1093/treephys/24.3.331 14704142 Thornton P. E. Running S. W. White M. A. (1997). Generating surfaces of daily meteorological variables over large regions of complex terrain. J. Hydrol. 190 214–251. 10.1016/S0022-1694(96)03128-9 Vitasse Y. (2013). Ontogenic changes rather than difference in temperature cause understory trees to leaf out earlier. New Phytol. 198 149–155. 10.1111/nph.12130 23347086 Vitasse Y. Basler D. (2013). What role for photoperiod in the bud burst phenology of European beech. Eur. J. For. Res. 132 1–8. 10.1007/s10342-012-0661-2 Vitasse Y. Delzon S. Dufrêne E. Pontailler J. Y. Louvet J. M. Kremer A. (2009). Leaf phenology sensitivity to temperature in European trees: do within-species populations exhibit similar responses? Agric. For. Meteorol. 149 735–744. 10.1016/j.agrformet.2008.10.019 Vitasse Y. Francois C. Delpierre N. Dufrene E. Kremer A. Chuine I. (2011). Assessing the effects of climate change on the phenology of European temperate trees. Agric. For. Meteorol. 151 969–980. 10.1016/j.agrformet.2011.03.003 Vitasse Y. Lenz A. Hoch G. Körner C. (2014a). Earlier leaf-out rather than difference in freezing resistance puts juvenile trees at greater risk of damage than adult trees. J. Ecol. 102 981–988. 10.1111/1365-2745.12251 Vitasse Y. Lenz A. Körner C. (2014b). The interaction between freezing tolerance and phenology in temperate deciduous trees. Front. Plant Sci. 5:541. 10.3389/fpls.2014.00541 25346748 Vitasse Y. Signarbieux C. Fu Y. S. H. (2018). Global warming leads to more uniform spring phenology across elevations. Proc. Natl. Acad. Sci. U.S.A. 115 1004–1008. 10.1073/pnas.1717342115 29279381 Way D. A. Montgomery R. A. (2015). Photoperiod constraints on tree phenology, performance and migration in a warming world. Plant Cell Environ. 38 1725–1736. 10.1111/pce.12431 25142260