Winter warming offsets one half of the spring warming effects on leaf unfolding.

Winter temperature-related chilling and spring temperature-related forcing are two major environmental cues shaping the leaf-out date of temperate species. To what degree insufficient chilling caused by winter warming would slow phenological responses to spring warming remains unclear. Using 27,071 time series of leaf-out dates for 16 tree species in Europe, we constructed a phenological model based on the linear or exponential function between the chilling accumulation (CA) and forcing requirements (FR) of leaf-out. We further used the phenological model to quantify the relative contributions of chilling and forcing on past and future spring phenological change. The results showed that the delaying effect of decreased chilling on the leaf-out date was prevalent in natural conditions, as more than 99% of time series exhibited a negative relationship between CA and FR. The reduction in chilling linked to winter warming from 1951 to 2014 could offset about one half of the spring phenological advance caused by the increase in forcing. In future warming scenarios, if the same model is used and a linear, stable correlation between CA and FR is assumed, declining chilling will continuously offset the advance of leaf-out to a similar degree. Our study stresses the importance of assessing the antagonistic effects of winter and spring warming on leaf-out phenology.

INTRODUCTION

The spring phenophases (e.g., leaf‐out date, flowering date) of temperate tree species are important functional traits that influence the performance and fitness of species (Wolkovich & Ettinger, 2014). Species that have advanced their spring phenology due to warming have increased their biomass and percent cover (Cleland et al., 2012). Changes in phenology also have major implications for ecosystem processes, including carbon uptake (Keenan et al., 2014; Richardson et al., 2010), water flux (Zha et al., 2010), surface energy budget (Hiestand & Carleton, 2019), as well as the mediation of feedbacks of terrestrial vegetation on climate (Peñuelas et al., 2009; Xu et al., 2020).

The spring phenology of temperate species is mainly affected by three factors: forcing, chilling, and photoperiod (Ettinger et al., 2020). The forcing temperature is positively related to the development rate of buds in spring so that the leaf‐out or flowering date may be expressed as a function of accumulated degree‐day or forcing requirement (FR; Peaucelle et al., 2019). Spring warming could advance the leaf‐out date because the species satisfies the FR more quickly in a warmer spring. However, experimental results showed that insufficient chilling in autumn and winter would increase the FR of leaf‐out (Wang, Wu, et al., 2020). Thus, winter warming has a delaying effect on the leaf‐out date due to the reduction in chilling related to warming itself. The effect of photoperiod on the leaf‐out date is hard to be detected in observational data, because forcing and photoperiod are inherently correlated in spring (Elmendorf & Ettinger, 2020; Wolkovich et al., 2021). A controlled experiment showed that only 2% of the 173 species are photoperiod‐sensitive in high chilling conditions (Zohner et al., 2016). Thus, this study mainly focused on the effects of chilling and forcing on the leaf‐out date, but only quantified the effect of photoperiod for Fagus sylvatica, a species known to be highly sensitive to photoperiod (Basler & Körner, 2012; Fu, Piao, et al., 2019; Zohner & Renner, 2015).

To date, most of the evidence of the chilling effect on the leaf‐out date was based on controlled experiments (Ettinger et al., 2020; Hänninen et al., 2019). The twigs or seedlings of temperate species need more FR to unfold their leaves if they experience insufficient chilling in natural conditions or low‐temperature treatments (Baumgarten et al., 2021; Flynn & Wolkovich, 2018; Laube et al., 2014; Nanninga et al., 2017). A global meta‐analysis of all published experiments found that species express the strongest response to chilling as opposed to forcing and photoperiod (Ettinger et al., 2020). Although recent studies have shown the delaying effect of reduced chilling on the leaf‐out date in natural conditions (Fu, Piao, et al., 2015; Fu, Zhao, et al., 2015), the effect of chilling still needs to be explicitly quantified. If the dominant role of chilling, as demonstrated in controlled experiments, also holds true in natural conditions, the advance of spring phenology would slow down, and even stop or reverse (showing a later trend) in future warming scenarios. In contrast, spring phenology would continue to advance if the dominant signal of climate change is increased forcing. Thus, quantifying the relative effects of chilling and forcing has critical forecasting implications.

Here, we used long‐term observations (≥20 years) of leaf‐out dates derived from 16 dominant tree species in Europe to test the relative effects of chilling and forcing on spring phenology. We first analyzed the relationship between chilling accumulation (CA) and FR for each time series to examine the effect of chilling on the leaf‐out date. Subsequently, we developed a phenological model based on the functions between CA and FR to separate the chilling and forcing contributions to phenological change. Finally, we compared the relative magnitudes of the delaying effect of reduced chilling and the advancing effect of increased forcing during the past decades (1951–2014) and in future warming scenarios (2022–2099). This study aimed to provide empirical evidence for the delaying effect of insufficient chilling on the leaf‐out date and to quantify to what degree chilling would slow phenological responses of plants to global warming.

MATERIALS AND METHODS

Phenological data

The phenological data (1951–2014) in Europe was derived from the Pan European Phenology Project (PEP725). The PEP725 database presently holds almost 12 million records of 265 plant species (Templ et al., 2018), which is available at http://www.pep725.eu/. We focused on the leaf‐out dates, defined as the date when the first leaf unfolds (BBCH code: 11; Meier, 2001). To have enough records to analyze the relationship between CA and FR, we selected the species that had been observed in at least 15 stations and with each station having at least 20 years of leaf‐out dates. The stations with less than 20‐year data were discarded and removed from the analysis. We thus ended up with 16 dominant tree species in Europe, including 13 deciduous broad‐leaved trees, 2 evergreen coniferous trees, and 1 deciduous coniferous tree (Table 1).

TABLE 1

Summary of the species investigated in this study

No.	Species	Common name	Life form	Number of stations (time series)	Number of records
1	Acer platanoides	Norway maple	DBT	25	589
2	Acer pseudoplatanus	Sycamore maple	DBT	33	821
3	Aesculus hippocastanum	Horse chestnut	DBT	3277	115,003
4	Alnus glutinosa	Common alder	DBT	1932	60,084
5	Betula pendula	Silver birch	DBT	3281	112,507
6	Betula pubescens	Downy birch	DBT	36	930
7	Fagus sylvatica	Common beech	DBT	2666	89,402
8	Fraxinus excelsior	Common ash	DBT	1930	62,065
9	Larix decidua	European larch	DCT	2445	81,305
10	Picea abies	Norway spruce	ECT	2641	89,068
11	Pinus sylvestris	Scots pine	ECT	2104	66,558
12	Populus tremula	European aspen	DBT	33	806
13	Quercus robur	Common oak	DBT	2795	93,572
14	Sorbus aucuparia	Mountain ash	DBT	863	19,720
15	Tilia cordata	small‐leaved lime	DBT	1182	31,926
16	Tilia platyphyllos	large‐leaved lime	DBT	1828	49,396

Abbreviations: DBS, deciduous broad‐leaved shrub; DBT, deciduous broad‐leaved tree; DCT, deciduous coniferous tree; ECT, evergreen coniferous tree.

Quality controls of phenological data were conducted to remove stations having large variations and unreliable records showing extreme early or late leaf‐out dates, possibly due to changes in the observers. We excluded individual time series, for which the standard deviation of leaf‐out dates across years was higher than 25 days (Zohner et al., 2018). Within each time series, leaf‐out records corresponding to more than three times the median absolute deviation were removed (Zohner et al., 2018). The selection and cleaning of the data, following the criterion described above, resulted in 3754 stations (Figure 1). A total of 873,752 records and 27,071 time series (site–species combinations) of leaf‐out dates were analyzed. The number of leaf‐out records, in particular, varied among species, ranging from 589 to 115,003 (Table 1).

FIGURE 1

Distribution of the 3754 stations in Europe. These stations have at least one species with time series of leaf‐out dates ≥20 years. Map lines delineate study areas and do not necessarily depict accepted national boundaries.

Climatic data

The past climatic data of Europe we used were the E‐OBS data set (Cornes et al., 2018) from the European Climate Assessment & Dataset (ECA&D) project (https://www.ecad.eu/), where the daily mean temperature during the period 1950–2020 was available at 1/10° resolution (~10 km). We matched each phenological station in a specific grid based on its latitude and longitude.

For the future climatic data, we obtained the daily resolution and bias‐corrected data simulated by the global climatic model HadGEM2‐ES, covering the period 2011–2099, with a resolution of 0.5 × 0.5° from the Fast Track input‐data catalog of the Inter‐Sectoral Impact Model Intercomparison Project (ISIMIP, https://www.isimip.org/). We used future climatic data under two representative concentration pathways (RCP 4.5 and 8.5). Although the ISIMIP data set has corrected the original data (retrieved from the CMIP5 archive) for systematic deviations of the simulated historical data from real observations (Warszawski et al., 2014), the resolution of 0.5° is too coarse for comparisons with the past climatic data. Therefore, we corrected the future temperature data for each phenological station by adding the difference in mean temperature between the past temperature data and future temperature data during the period in common (2011–2020).

Relationship between CA and FR

For each leaf‐out record, CA was calculated as the accumulated chilling units from a constant date of the previous year to the median value of the multi‐year leaf‐out date (Equation 1). The chilling unit measures the rate of chilling in response to temperature (Equation 2). where CA is the chilling accumulation for each leaf‐out record; t median is the median value of the current leaf‐out time series. t 0 is the starting date of chilling accumulation, which could be set as 1 November of the previous year following previous studies (Laube et al., 2014; Wang, Wang, et al., 2020). However, a recent warming experiment found a much earlier start date of chilling since October (Beil et al., 2021). Thus, we also calculated CA with t 0 = 1 October of the previous year. CU(t) is the chilling unit on date t. T(t) is the daily mean temperature (°C) on date t; T c is the upper‐temperature threshold for effective chilling. Many studies have reported that temperatures below 5 or 7°C are effective in breaking the dormancy of temperate species (Fu, Zhao, et al., 2015; Laube et al., 2014). A recent study found that the effective chilling temperature could reach up to 10°C (Baumgarten et al., 2021): thus, we also calculated the CA with T c = 10°C. Overall, six chilling algorithms (2 t 0 × 3 T c) were used to calculate the CA. Many other chilling models also assume that temperature below 0°C is ineffective for breaking dormancy, but these models could not reproduce the negative correlation between CA and FR in the observation data (Wang, Wu, et al., 2020). A recent experiment proved that freezing temperature (e.g., −2°C) is also effective in breaking dormancy (Baumgarten et al., 2021). Therefore, we did not test chilling models assuming that temperature below 0°C is inefficient in breaking dormancy.

The calculation of FR was based on growing degree‐days (Equation 3), which integrated the forcing units [daily mean temperatures above a specific threshold (T b) from a fixed starting date (t 1) to the observed leaf‐out date (Equation 4)]. T b is often set as 5°C (Fu et al., 2013; Peaucelle et al., 2019; Wang, Wang, et al., 2020) or 0°C (Basler & Körner, 2012; Heide, 1993; Piao et al., 2015). t 1 is commonly set as 1 January (Fu, Piao, et al., 2015; Peaucelle et al., 2019; Wang, Wang, et al., 2020), but 1 February has also been reported by others (Cannell & Smith, 1983; Fu et al., 2013). To test whether the choice of the forcing algorithms would affect the result, we used three forcing algorithms (T b = 5°C, t 1 = 1 January; T b = 0°C, t 1 = 1 January; T b = 5°C, t 1 = 1 February) based on different threshold temperatures and starting dates. where FR is the forcing requirement for leaf‐out and FU(t) is the forcing unit at day of year t. t LOD is the observed leaf‐out date, and T(t) is the daily mean temperature (°C) on date t.

To examine how chilling affects the leaf‐out date, we calculated the Pearson’ r between CA and FR (r CA–FR) for each time series of leaf‐out. Because six chilling algorithms and three forcing algorithms were adopted, 18 r CA–FR values were calculated for different combinations of chilling and forcing algorithms.

Phenological models based on the CA–FR relationship

We calibrated the specific function to establish the relationship between CA and FR for each time series of leaf‐out dates. According to previous studies, the relationship between CA and FR may be linear (Nanninga et al., 2017; Wang, Wu, et al., 2020; Zohner & Renner, 2015). Thus, we fitted the following linear function for each time series of leaf‐out dates. where FR is the forcing requirement for spring events, CA is the chilling accumulation, and f and g are the parameters.

However, several studies described the relationship between CA and FR as an exponential function (Cannell & Smith, 1983; Lin et al., 2022; Man et al., 2017). Thus, we further fitted the following function: where a, b, and c are the parameters.

For each time series, the parameters of the linear or exponential function were fitted using the least square method. At several stations, the FR or CA was very small (FR <10°C or CA <10 days), suggesting the species may have adopted a different threshold temperature to start growth or break dormancy due to the adaptation to the local climate. We removed these stations for further analysis because we could not obtain a stable estimate of the parameters for time series with small CA or FR.

For each time series (species−station combination), we simulated the leaf‐out date on a daily basis since 1 January. The CA was calculated based on the past (1951–2014) or future (2022–2099) climate data. The FR for leaf‐out was estimated using the predefined function between FR and CA (Equations 5 or 6). The leaf‐out date was determined as the day when the accumulated degree‐day began to exceed the FR. We assessed the performance of the phenological models with two indicators: the coefficient of determination (R 2) and the root‐mean‐square error (RMSE) between the simulated and observed dates. Although the phenological models based on different chilling and forcing algorithms generate interrelated leaf‐out dates (e.g., Figure S1), the RMSE and R 2 between the simulated and observed dates varied greatly among combinations of chilling and forcing algorithms (e.g., Figure S2). For all species, the phenological model based on the chilling algorithm with T c = 5°C and t 0 = 1 November and the forcing algorithm with T b = 0°C and t 1 = 1 January was most accurate since it showed the lowest RMSE and highest R 2 between the simulated and the observed leaf‐out date among 18 combinations. Thus, we only reported the results based on this optimal chilling and forcing algorithm combination.

Quantifying chilling and forcing contributions

We developed a novel method to quantify the contribution of chilling and forcing on phenological change. This method separated the simulated leaf‐out date into three parts: the leaf‐out date during the reference period (1961–1990), the chilling contribution, and the forcing contribution relative to the reference period, as expressed in the following equation: where P (FR ) is the simulated leaf‐out date (the date when FR is fulfilled) at a certain year i. P ref(FRref) is the mean leaf‐out date during the reference period (1961–1990). FRref represents the mean forcing requirement during the reference period. C is the chilling contribution in a certain year i, representing the deviation in the leaf‐out date from the reference period induced by changes in the chilling condition; F is the forcing contribution in a certain year i, representing the deviation in the leaf‐out date from the reference period induced by changes in the forcing condition.

The above quantitative splitting of the predicted phenological change into the effects caused by changes in forcing and those caused by changes in chilling is realized by the following steps:

For each time series of leaf‐out dates, the forcing requirement in a year i (FR ) is calculated based on the amount of chilling (CA , Equation 1) and the CA–FR curve (Equations 5 or 6). Subsequently, P (FR ), that is, the simulated leaf‐out date at a certain year i, is calculated as the date when FR is fulfilled.

We defined 1961–1990 as the reference period, and the reference forcing requirement (FRref) is calculated as the mean FR during this period (Equation 8), and the reference leaf‐out date [P ref(FRref)] is calculated as the mean P (FR ) during the reference period (Equation 9).

where F and P ref(FRref) describe the same parameters as Equation (7). P (FRref) is calculated as the date when FRref is fulfilled in year i.

The forcing contribution (F ) is calculated as the difference in the date when FRref is fulfilled in year i and the reference leaf‐out date [P ref(FRref)], because the year‐to‐year variation of forcing temperature alters the date when FRref is fulfilled (Equation 10; Figure 2).

FIGURE 2

Schematic diagram showing how to separate the contribution of chilling and forcing on spring phenological change. (a) Chilling accumulation (CA)‐forcing requirement (FR) curve for the leaf‐out date. CAref and FRref: mean CA and HR during the reference period (1961–1990). In a warm year i, CA is smaller than CAref, and then FR is larger than FRref. (b) The methods to separate the effect of chilling and forcing. P ref(FRref) is the mean leaf‐out date during the reference period; P (FR ) is the simulated leaf‐out date at year i, that is, the date when the accumulated growing degree days (with a threshold T b, starting date t 1) is larger than FR . P (FRref) is the date in year i when FRref is fulfilled. The chilling contribution (C ) is calculated as the difference between P (FR ) and P (FRref) because the year‐to‐year variation in FR is caused by the year‐to‐year variation of chilling. The forcing contribution (F ) is calculated as the difference in P (FRref) and Pref(FRref) because the year‐to‐year variation of forcing temperature alters the date when FRref is fulfilled.

The chilling contribution (C ) is calculated as the difference between P (FR ) and P (FRref) because the year‐to‐year variation in FR is caused by the year‐to‐year variation of chilling (Equation 11; Figure 2) if any potential contribution of photoperiod is neglected.

where C and P (FR ) describe the same parameters as Equation (7). P (FRref) describes the same parameters as Equation (10). Following Equations (10) and (11), Equation (7) is a mathematical necessity.

Following the methods above, we estimated the interannual changes in chilling and forcing contributions during the past period (1951–2014) and future period (2022–2099, under RCP 4.5 and 8.5) for each species−station combination. At the species level, we calculated the time series of the chilling contribution, forcing contribution, and leaf‐out date averaged from all stations. To obtain a regional time series, we averaged the time series of the chilling contribution, forcing contribution, and leaf‐out date from all species. The slopes of the linear regression of the chilling contribution, forcing contribution, and leaf‐out date against year were used to measure the trends of each variable.

Estimating photoperiod contributions on leaf‐out date of F. sylvatica

In previous experiments, F. sylvatica is identified as a photoperiod‐sensitive species whose leaf‐out date is significantly delayed under short photoperiod than long photoperiod (Basler & Körner, 2012; Fu, Piao, et al., 2019; Zohner & Renner, 2015). In order to quantify the effect of photoperiod on the leaf‐out date of F. sylvatica in natural conditions, we further developed a CA‐photoperiod model. In Section 2.4, the phenological model is based on the linear function between CA and FR (hereinafter called the CA‐only model). However, the CA‐photoperiod model is based on the linear function between CA and FR incorporating the effect of photoperiod (FRP). where FRP is the forcing requirement incorporating the effect of photoperiod and t LOD is the observed leaf‐out date. t 1 is the starting date of degree‐day accumulation (set as 1 January). FU(t) is the forcing unit at the day of year t (Equation 4, T b is set as 0°C). AE represents the additional effect of photoperiod on leaf‐out.

We assumed that the AE is linearly correlated with day length. A long photoperiod accelerates the accumulation of the FR (higher AE), while, on the contrary, a short photoperiod suppresses it. Thus, AE >1 if day length >12 h, AE = 1 if day length = 12 h, and AE <1 if day length <12 h. where DL is the day length, which is calculated as a function of latitude and day of the year (Forsythe et al., 1995). h and i are parameters. According to Zohner and Renner (2015), the effect of photoperiod is more significant in low chilling conditions than that in high chilling conditions. Therefore, h and i should vary with CA.

In order to quantify h and i, we defined a variable w as the ratio of the additional effect of a 16 h photoperiod to the additional effect of an 8 h photoperiod: According to the experimental results of Zohner and Renner (2015), under the same chilling conditions, the FR of leaf‐out under 16 h photoperiod was significantly smaller than that under 8 h photoperiod. Meanwhile, under the same photoperiod conditions, CA was linearly and negatively correlated with FR (Figure S3). Thus, w could be estimated as the ratio of FR under 8 h photoperiod to FR under 16 h photoperiod. where FR8 and FR16 correspond to the FR under an 8 or 16 h photoperiod in the experiment, respectively, CA is the chilling accumulation, and f 8, g 8, f 16, and g 16 are parameters of the linear function. According to the experimental results of Zohner and Renner (2015), f 8 = 2397.0, g 8 = −17.14, f 16 = 755.9, and g 16 = −3.97 (Figure S3).

The solution of the system of linear equations (Equations 14 and 15) was Overall, FRP at each day of the year could be calculated as follows: (i) by calculating the CA based on the climate data (Equations 1 and 2 with T c = 5°C, t 0 = 1 November); (ii) by calculating w based on the CA and Equation (16); (iii) by calculating h and i based on w, Equations (17) and (18); (iv) by calculating the day length (DL) based on equations from Forsythe et al. (1995); (v) by calculating AE(DL) based on h, i, DL and Equation (13); (vi) by calculating FRP based on Equation (12).

At each station of F. sylvatica, we fitted the linear function between the CA and FRP. Likewise the CA‐only model, we simulated the leaf‐out date on a daily basis since 1 January. The CA was calculated based on the past (1951–2014) or future climate data (2022–2099). FRP for the leaf‐out was estimated using the predefined function between FRP and CA. The day when the accumulated photoperiod‐associated degree‐day began to exceed FRP was estimated as the leaf‐out date. The contribution of photoperiod on the leaf‐out date of F. sylvatica is estimated as the difference in the simulated leaf‐out date between the CA‐photoperiod model and the CA‐only model.

RESULTS

Delaying effect of decreased chilling on the leaf‐out date

The winter temperature (previous November to February) averaged from all phenological stations increased by 0.24°C/decade (p < .01) from 1951 to 2014, respectively (Figure 3a). At the interannual scale, the winter temperature correlated negatively with the CA (calculated as the number of days with a daily mean temperature <5°C during previous November to February, Figure 3b). The CA thus decreased significantly by 3.86 days/decade (p < .01) in Europe (Figure 3a).

FIGURE 3

Warming‐related reductions in chilling accumulation (CA) delay leaf‐out dates of temperate species. (a) Annual mean winter (previous November to February) temperature and winter CA (number of days with a daily mean temperature <5°C during previous November to February) from 1951 to 2014 averaged from 3754 European stations. The dotted lines are linear regressions (the slopes are shown in the parentheses. **p < .01). (b) CA correlates negatively and significantly with winter temperature at the interannual scale (**p < .01). (c) The increase in CA decreases the forcing requirement (FR) of leaf‐out date. The examples of Acer platanoides at a European station (48.15°N, 15.15°E) are shown. The CA–FR relationship is fitted with linear and exponential functions. (d) The proportion of the time series with negative and positive Pearson's r between CA and FR (r CA‐HR) for each species. The number of time series (i.e., stations) for each species is shown in parentheses.

We further calculated the FR [based on Equations (3) and (4) with T b = 5°C and t 1 = 1 January] and CA [based on Equations (1) and (2) with T c = 5°C and t 0 = 1 November] for each record of the leaf‐out date. When analyzing the Pearson's r between CA and FR (r CA–FR) for each time series, we found that the FR was negatively correlated with the CA in most cases (see the examples in Figure 3c). The proportion of the time series with negative r CA–FR varied among species, ranging from 93.8% to 100% (37.5%–95.8% at p < .05, Figure 3d). We further used other five chilling and two forcing algorithms to calculate the CA and FR. The results showed that the prevalent negative r CA–FR was independent of the algorithms applied to calculate chilling or forcing (Figure S4). Thus, over the past six decades, the signal of the delaying effect of decreased chilling on the leaf‐out date was shown to be strong.

Contribution of chilling and forcing on past phenological change

We developed a phenological model for each time series of leaf‐out dates using the linear function between CA and FR (Equation 5) based on the optimal chilling and forcing algorithms. The leaf‐out dates could be simulated as the dates when the CA‐associated FR was satisfied. The simulated and observed leaf‐out dates corresponded well given the general uncertainty in phenological modeling (Figure 4). The RMSE and the coefficient of determination (R 2) of all the records ranged from 6.27 to 9.81 days and 0.40 to 0.81, respectively. Furthermore, our model was able to reproduce the past spatial pattern of leaf‐out dates when comparing the simulated and observed leaf‐out dates at the station level (Figure S5). Thus, this model could be used to separate the contribution of chilling and forcing on phenological change.

FIGURE 4

Comparisons between observed and simulated leaf‐out dates of 16 temperate species. For each time series, leaf‐out dates were simulated by a process‐based phenological model based on the linear regression between chilling accumulation (CA) and forcing requirement (FR). For each species, the root‐mean‐square error (RMSE) and the coefficient of determination (R 2) are shown at the bottom right, and the number of records is shown in parentheses.

We estimated the chilling contribution to the past phenological change (1951–2014) as the deviation in the leaf‐out date from the reference period (1961–1990) induced by changes in the chilling. The chilling contribution, averaged from all the species, showed a positive trend of 0.15 days/year (p < .01), suggesting an enhanced delaying effect of chilling caused by winter warming (Figure 5a). The chilling contributions were positive in almost all years after 1988 (except in 1996 and 2006). Especially in 2007, the insufficient chilling caused a delaying effect of 20.7 days relative to the reference period (1961–1990). The trends of the chilling contribution varied among species (Figure 5c; Figure S6), ranging from 0.09 to 0.18 days/year. For each specific species, the trends of the chilling contribution also differed among stations with an interquartile range of 0.06–0.09 days/year (Figure 5b).

FIGURE 5

Interannual changes in the leaf‐out date of temperate species during the past decades (1951–2014) simulated by a phenological model based on the linear function between FR and CA. DOY, day of the year. The chilling (forcing) contribution represents the impact of interannual changes in chilling (forcing) on the leaf‐out date relative to the reference period (1961–1990). (a) The mean time series of the chilling contribution, forcing contribution and leaf‐out date from 1951 to 2014. Error bars, mean values ± SD (n = 16 species). The trendlines are shown in the dotted line. The regression slope of the chilling contribution, forcing contribution, and leaf‐out date against year are shown in parentheses (**p < .01). (b) Boxplot of trends in chilling and forcing contributions for all time series of each species. Bottom and top of boxes: 25th and 75th percentiles; bands within boxes: medians; whiskers: 25th and 75th percentiles. (c) The trend (1951–2014) of chilling and forcing contributions for each species. (d) The relationship between chilling and forcing contributions across years. The linear regression lines between chilling and forcing contributions are shown (**p < .01).

In order to examine to what degree the delaying effect of chilling could offset the advancing effect of forcing, we estimated the forcing contribution as the deviation in the leaf‐out date from the reference period (1961–1990) induced by changes in the forcing temperature. The forcing contribution, averaged from all the species, enhanced during 1951–2014 with a trend of −0.27 days/year (p < .01). The forcing contributions were negative in almost all years after 1988 (except for 1996, 2006, and 2013), indicating that the recent spring warming had an advancing effect on the leaf‐out date of temperate species (Figure 5a).

The magnitude of the increasing trends of the chilling contribution was weaker than the decreasing trends of the forcing contribution (Figure 5a), and this result was also robust at the species level (Figure 5c) and station level (Figure 5b). Therefore, the simulated leaf‐out date became earlier by 0.12 days/year (p < .01, Figure 5a) from 1951 to 2014, which matched previous results based on long‐term observations (Cayan et al., 2001; Ge et al., 2015; Menzel et al., 2006). Furthermore, the chilling contribution correlated significantly with the forcing contribution at the interannual scale (Figure 5d). The slope of the forcing contribution against the chilling contribution indicated that 1 day delay in the leaf‐out date induced by chilling corresponded to 1.91 days advance induced by forcing (Figure 5d). Thus, the winter warming‐related reductions in chilling could, on average, offset about 52% of spring phenological advances caused by spring warming. However, this proportion varied among species, ranging from 38% to 62% (Figure 6).

FIGURE 6

Relationship between the chilling and forcing contributions to the changes in leaf‐out date across years (1951–2014) for 16 species in Europe. The chilling and forcing contributions were simulated by a process‐based phenological model based on the linear function between chilling accumulation (CA) and forcing requirement (FR). The linear regression lines (dotted line) and equations for each species are shown. **p < .01.

We further developed a phenological model based on the exponential function between CA and FR (Equation 6). Similar to the linear function‐based model, the RMSE and R 2 of the exponential function‐based model ranged from 6.26 to 9.84 days and 0.39 to 0.95, respectively (Figure S7). Moreover, the exponential function‐based model also could accurately reproduce the past spatial pattern of leaf‐out dates when comparing the simulated and observed leaf‐out dates at the station level (Figure S5). Thus, the exponential function‐based model had a similar accuracy to the linear function‐based model. Furthermore, the trend and amplitude of the chilling and forcing contributions from 1951 to 2014 derived from the exponential function‐based model were similar to the linear function‐based model (Figure 7).

FIGURE 7

Interannual changes in the leaf‐out date of temperate species (1951–2014) simulated by a phenological model based on the exponential function between FR and CA. DOY, day of the year. The chilling (forcing) contributions represent the impact of interannual changes in chilling (forcing) on the leaf‐out date relative to the reference period (1961–1990). (a) The mean time series of the chilling contribution, forcing contribution and leaf‐out date from 1951 to 2014. Error bars, mean values ± SD (n = 16 species). The trendlines are shown in the dotted line. The regression slope of the chilling contribution, forcing contribution, and leaf‐out date against year are shown in parentheses (**p < .01). (b) The relationship between the chilling and forcing contributions across years based on the linear function and exponential function between FR and CA. The linear regression lines between the chilling and forcing contributions based on the exponential function are shown (**p < .01).

Future changes in leaf‐out dates

We first predicted the future changes in the leaf‐out dates of temperate species from 2022 to 2099 under RCP 4.5 and RCP 8.5 scenarios using the linear function‐based model. In both scenarios, the leaf‐out date would continue to become earlier. The trends of the leaf‐out dates were −0.016 days/year (not significant) under RCP 4.5 averaged from 16 tree species (Figure 8a). Under RCP 8.5, the trends increased to −0.166 days/year (p < .01; Figure 8b). The earlier leaf‐out date could be attributed to the weaker trends of the chilling contribution than the forcing contribution (Figure 8a,b). The future trends in the leaf‐out date were stronger under RCP 8.5 with intensified warming than under RCP 4.5, but the trends varied among species (Figures S8 and S9). For future scenarios (2022–2099), the interannual chilling contribution (averaged from all the species) correlated significantly with the forcing contribution (Figure 8c). Under RCP 4.5, a 1 day delaying effect of chilling corresponded to a 1.37 days advancing effect of forcing, and under RCP 8.5 the same corresponded to a 1.65 days advancing effect of forcing (Figure 8c). Thus, if the same model is used and a linear, stable correlation between CA and FR is assumed in the future warming scenarios, it is likely that the delaying effect of chilling will substantially offset the advance of leaf unfolding to a similar degree in comparison to the past period.

FIGURE 8

Future changes in the leaf‐out date of temperate species simulated by a phenological model based on the linear function between FR and CA for 16 tree species in Europe during the period 2022–2099 and under the climate scenarios of RCP 4.5 and 8.5. DOY, day of the year. (a) The mean time series of the chilling contribution, forcing contribution (relative to the reference period 1961–1990) and leaf‐out date under RCP 4.5. (b) The mean time series of the chilling contribution, forcing contribution and leaf‐out date under RCP 8.5. Error bars and mean values ± SD (n = 16 species) are shown. The trendlines are shown in the dotted line. The regression slope of the chilling contribution, forcing contribution, and leaf‐out date against year are shown in parentheses (**p < .01). (c) The relationship between the chilling and forcing contributions across years (2022–2099) based on the linear function between FR and CA. The linear regression lines between the chilling and forcing contributions based on the linear function‐based phenological model are shown (**p < .01).

When using the exponential function between CA and FR to simulate the leaf‐out date, the trends of the leaf‐out dates averaged from 16 tree species were 0.026 and 0.016 days/year (both not significant) under RCP 4.5 and RCP 8.5, respectively (Figure 9a,b). Thus, if the exponential function between CA and FR is assumed in the future warming scenarios, the spring phenological advance may reverse in the future. The reason is that the delaying effect of decreased chilling simulated by the exponential function is stronger than that simulated by the linear function, and the advancing effect of increased forcing is the same between the two types of function (Figure 9c).

FIGURE 9

Future changes in the leaf‐out date of temperate species simulated by a phenological model based on the exponential function between FR and CA in Europe during the period 2022–2099 and under the climate scenarios of RCP 4.5 and 8.5. DOY, day of the year. (a) The mean time series of the chilling contribution, forcing contribution and leaf‐out date under RCP 4.5. (b) The mean time series of the chilling contribution, forcing contribution and leaf‐out date under RCP 8.5. Error bars, mean values ± SD (n = 16 species). The trendlines are shown in the dotted line. The regression slope of the chilling contribution, forcing contribution, and leaf‐out date against year are shown in parentheses (**p < .01). (c) The relationship between the chilling and forcing contributions across years (2022–2099) based on the linear function and exponential function between FR and CA. The linear regression lines between the chilling and forcing contributions calculated from the exponential function‐based phenological model are shown (**p < .01).

Contribution of photoperiod to leaf‐out dates of European beech

We specially developed a CA‐photoperiod model (incorporating the effect of photoperiod) for simulating the leaf‐out date of F. sylvatica. The photoperiod contribution to the leaf‐out date of F. sylvatica could be estimated as the difference in the simulated leaf‐out date between the CA‐photoperiod model and the CA‐only model. The CA‐photoperiod model could slightly improve the simulation of leaf‐out date compared to the original model (Figure S10), and predicted similar trends in the leaf‐out date with the CA‐only model (Figure 10a,c,e). At the interannual scale, the photoperiod contribution was positive when the leaf‐out date was earlier, but was negative when the leaf‐out date was later (Figure 10b,d,f). However, the amplitude of the photoperiod contribution among years was only ±6 days.

FIGURE 10

Comparisons between the change in multi‐station mean leaf‐out date of Fagus sylvatica simulated by the CA‐only model and CA‐photoperiod model. (a, b) Annual leaf‐out date (a) and the relationship between photoperiod contribution and leaf‐out date during the period 1951–2014 (b). The photoperiod contribution represents the difference in the number of days between the leaf‐out dates simulated by the CA‐photoperiod model and simulated by the CA‐only model. (c, d) Annual leaf‐out date (c) and the relationship between photoperiod contribution and leaf‐out date during the period 2022–2099 under the RCP 4.5 scenario (d). (e, f) Annual leaf‐out date (e) and the relationship between photoperiod contribution and leaf‐out date during the period 2022–2099 under the RCP 8.5 scenario (f). The linear slopes of each time series are shown. **p < .01. *p < .05.

DISCUSSION

As results from the literature show, the main evidence that a reduction in chilling could delay the leaf‐out date has come from the controlled experiment on saplings or twigs of temperate species (Du et al., 2019; Flynn & Wolkovich, 2018; Man et al., 2017; Nanninga et al., 2017). This study provided empirical evidence in natural conditions since the CA was negatively correlated to the FR of leaf‐out in more than 99% of the time series. A global meta‐analysis of all published experiments found that the effect of chilling on the spring phenology of species was stronger than the effect of forcing (Ettinger et al., 2020). However, this study suggested that the delaying effect of decreased chilling on the leaf‐out date was always weaker than the advancing effect of increased forcing. Such differences between experiments and observed data may be attributed to the wider range of chilling treatments in the experiments. For example, the difference in the CA under low chilling and high chilling treatments reached 77 days in an experiment focused on 36 woody species in Europe (Laube et al., 2014). However, the CA in the field only reduced by 24.3 days during the period 1951–2014 in Europe (Figure 3a). Therefore, a stronger response of spring phenology to chilling is found if experiments expand observed chilling beyond the current conditions. In current climate conditions, the counteracting effects of chilling are insufficient to stop the advance of leaf unfolding in response to warmer forcing temperatures.

This study used two hypothesised CA–FR functions to simulate the leaf‐out date of temperate species. Compared to the linear function, the exponential function assumes that the FR increased exponentially with the decrease in CA (e.g., Figure 3c), and the effect of chilling would be amplified in an extremely warm winter. In the past period (1951–2014), the trend and amplitude of the chilling and forcing contributions derived from the exponential function‐based model were similar to the linear function‐based model. However, in future warming scenarios, the delaying effect of decreased chilling become stronger and even exceed the advancing effect of increased forcing in several extremely warm winters (Figure 9c). Therefore, under the hypothesis of an exponential function between CA and FR, the future advance in the leaf‐out date may stop for tree species with the relatively high chilling requirement. We need to further validate the forms of relationship between CA and FR since the linear and exponential functions were shown to both accurately describe the phenological change that happened in the past, but to differ in performance regarding future warming scenarios.

Besides forcing and chilling, photoperiod has been reported to be a notable factor affecting spring phenology (Way & Montgomery, 2015), although only 35% of species relied on spring photoperiod as a leaf‐out signal in low chilling conditions (Zohner et al., 2016). We selected F. sylvatica, a well‐known photoperiod‐sensitive species (Basler & Körner, 2012; Fu, Piao, et al., 2019; Zohner & Renner, 2015), to investigate the relative contribution of chilling and photoperiod on phenological change. The results confirmed the hypothesis that day length helps temperate deciduous trees to leaf‐out at the optimal time (Fu, Zhang, et al., 2019). In warm springs, when leaf‐out is early, the short day length mitigates the warming‐induced advancement of leaf‐out. In cold springs, when leaf‐out is late, the long day length promotes leaf‐out, ensuring that trees do not leaf‐out too late. However, the amplitude of the photoperiod contribution was significantly smaller than that of the chilling contribution (±40 days). Therefore, decreased chilling plays a dominant role in mitigating the advance of leaf‐out dates.

The phenological models developed in this study, which describe the response of spring phenology of temperate tree species to warming, are valid within the temperature range we have observed during the past several decades (1951–2014). The CA–FR functions we established might not be stable if chilling declines beyond the range of current temperature conditions. The plants might be adapted only to a certain amount of chilling, and might respond unexpectedly (likely stronger) if chilling becomes much more incomplete. Therefore, projected changes in future spring phenology in this study are based on the assumption that the CA–FR function observed in the past is stable in the future and may overestimate the spring phenological advance.

Overall, using 27,071 time series of leaf‐out dates for 16 temperate tree species in Europe, we showed that the delaying effect of chilling on the leaf‐out date found in experiments was prevalent in natural conditions, as more than 99% of time series exhibited a negative relationship between CA and FR. Based on the linear function between CA and FR, we first quantified the relative contribution of chilling and forcing on phenological change. The winter warming‐related reductions in chilling during the time span 1951–2014 could offset about 52% of spring phenological advances caused by increased forcing. Assuming that the linear CA–FR function is stable in the future warming scenarios, the delaying effect of chilling would be enhanced in the future and slow down the spring phenological advance to a similar degree. However, it is possible that the relationship between FR and CA is or would become exponential (maybe in particular if temperature rises), and then the contribution of chilling could fully compensate for forcing in extremely warm years. More controlled experiments and continuous monitoring of spring phenology are needed to investigate whether a linear CA–FR relationship holds in a broader range of temperatures than the one we can observe by now, or whether the relationship appears to be exponential if temperatures rise.

CONFLICTS OF INTEREST

The authors declare no known conflicts of interest.

Appendix S1

Click here for additional data file.