Received: 8 December 2021 | Accepted: 15 July 2022 DOI: 10.1111/gcb.16358 RESEARCH ARTICLE Winter warming offsets one half of the spring warming effects on leaf unfolding Huanjiong Wang1 | Junhu Dai1,2,3 | Josep Peñuelas4,5 | Quansheng Ge1,2 | Yongshuo H. Fu6 | Chaoyang Wu1 1 Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographic Abstract Sciences and Natural Resources Research, Winter temperature-­related chilling and spring temperature-­related forcing are two Chinese Academy of Sciences, Beijing, China major environmental cues shaping the leaf-­out date of temperate species. To what 2 University of Chinese Academy of degree insufficient chilling caused by winter warming would slow phenological re- Sciences, Beijing, China 3 sponses to spring warming remains unclear. Using 27,071 time series of leaf-­out dates China–­Pakistan Joint Research Center on Earth Sciences, Chinese Academy of for 16 tree species in Europe, we constructed a phenological model based on the Sciences-­Higher Education Commission of linear or exponential function between the chilling accumulation (CA) and forcing Pakistan, Islamabad, Pakistan 4 requirements (FR) of leaf-­out. We further used the phenological model to quantify CSIC, Global Ecology Unit CREAF-­C SIC-­ UAB, Bellaterra, Barcelona, Spain the relative contributions of chilling and forcing on past and future spring phenologi- 5 CREAF, Cerdanyola del Valles, Barcelona, cal change. The results showed that the delaying effect of decreased chilling on the Spain 6 leaf-­out date was prevalent in natural conditions, as more than 99% of time series College of Water Sciences, Beijing Normal University, Beijing, China 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 phe- Correspondence Huanjiong Wang and Chaoyang Wu, nological advance caused by the increase in forcing. In future warming scenarios, if Key Laboratory of Land Surface Pattern the same model is used and a linear, stable correlation between CA and FR is assumed, and Simulation, Institute of Geographic Sciences and Natural Resources Research, declining chilling will continuously offset the advance of leaf-­out to a similar degree. Chinese Academy of Sciences, Beijing Our study stresses the importance of assessing the antagonistic effects of winter and 100101, China. Email: wanghj@igsnrr.ac.cn and wucy@ spring warming on leaf-­out phenology. igsnrr.ac.cn KEYWORDS Funding information chilling, climate warming, forcing, future change, leaf unfolding, phenology National Key R&D Program of China, Grant/Award Number: 2018YFA0606102; National Natural Science Foundation of China, Grant/Award Number: 41871032 and 42125101; Chinese Academy of Sciences, Grant/Award Number: 2018070 and 2018RC101 1 | I NTRO D U C TI O N Species that have advanced their spring phenology due to warm- ing have increased their biomass and percent cover (Cleland The spring phenophases (e.g., leaf-­out date, flowering date) of tem- et al., 2012). Changes in phenology also have major implications for perate tree species are important functional traits that influence the ecosystem processes, including carbon uptake (Keenan et al., 2014; performance and fitness of species (Wolkovich & Ettinger, 2014). Richardson et al., 2010), water flux (Zha et al., 2010), surface energy This is an open access article under the terms of the Creative Commons Attribution-­NonCommercial-­NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-­commercial and no modifications or adaptations are made. © 2022 The Authors. Global Change Biology published by John Wiley & Sons Ltd. Glob Change Biol. 2022;28:6033–6049.  wileyonlinelibrary.com/journal/gcb | 6033 6034 | WANG et al. budget (Hiestand & Carleton, 2019), as well as the mediation of feed- the advancing effect of increased forcing during the past decades backs of terrestrial vegetation on climate (Peñuelas et al., 2009; Xu (1951–­2014) and in future warming scenarios (2022–­2099). This et al., 2020). study aimed to provide empirical evidence for the delaying effect of The spring phenology of temperate species is mainly affected by insufficient chilling on the leaf-­out date and to quantify to what de- three factors: forcing, chilling, and photoperiod (Ettinger et al., 2020). gree chilling would slow phenological responses of plants to global The forcing temperature is positively related to the development warming. 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 re- quirement (FR; Peaucelle et al., 2019). Spring warming could advance 2 | M ATE R I A L S A N D M E TH O DS the leaf-­out date because the species satisfies the FR more quickly in a warmer spring. However, experimental results showed that insuf- 2.1 | Phenological data ficient chilling in autumn and winter would increase the FR of leaf-­ out (Wang, Wu, et al., 2020). Thus, winter warming has a delaying The phenological data (1951–­2014) in Europe was derived from effect on the leaf-­out date due to the reduction in chilling related to the Pan European Phenology Project (PEP725). The PEP725 data- warming itself. The effect of photoperiod on the leaf-­out date is hard base presently holds almost 12 million records of 265 plant species to be detected in observational data, because forcing and photope- (Templ et al., 2018), which is available at http://www.pep725.eu/. riod are inherently correlated in spring (Elmendorf & Ettinger, 2020; We focused on the leaf-­out dates, defined as the date when the first Wolkovich et al., 2021). A controlled experiment showed that only leaf unfolds (BBCH code: 11; Meier, 2001). To have enough records 2% of the 173 species are photoperiod-­sensitive in high chilling con- to analyze the relationship between CA and FR, we selected the spe- ditions (Zohner et al., 2016). Thus, this study mainly focused on the cies that had been observed in at least 15 stations and with each sta- effects of chilling and forcing on the leaf-­out date, but only quanti- tion having at least 20 years of leaf-­out dates. The stations with less fied the effect of photoperiod for Fagus sylvatica, a species known to than 20-­year data were discarded and removed from the analysis. be highly sensitive to photoperiod (Basler & Körner, 2012; Fu, Piao, We thus ended up with 16 dominant tree species in Europe, includ- et al., 2019; Zohner & Renner, 2015). ing 13 deciduous broad-­leaved trees, 2 evergreen coniferous trees, To date, most of the evidence of the chilling effect on the and 1 deciduous coniferous tree (Table 1). leaf-­ out date was based on controlled experiments (Ettinger Quality controls of phenological data were conducted to remove et al., 2020; Hänninen et al., 2019). The twigs or seedlings of tem- stations having large variations and unreliable records showing ex- perate species need more FR to unfold their leaves if they experi- treme early or late leaf-­out dates, possibly due to changes in the ence insufficient chilling in natural conditions or low-­temperature observers. We excluded individual time series, for which the stan- treatments (Baumgarten et al., 2021; Flynn & Wolkovich, 2018; dard deviation of leaf-­out dates across years was higher than 25 days Laube et al., 2014; Nanninga et al., 2017). A global meta-­analysis of (Zohner et al., 2018). Within each time series, leaf-­out records cor- all published experiments found that species express the strongest responding to more than three times the median absolute deviation response to chilling as opposed to forcing and photoperiod (Ettinger were removed (Zohner et al., 2018). The selection and cleaning of et al., 2020). Although recent studies have shown the delaying ef- the data, following the criterion described above, resulted in 3754 fect of reduced chilling on the leaf-­out date in natural conditions stations (Figure 1). A total of 873,752 records and 27,071 time se- (Fu, Piao, et al., 2015; Fu, Zhao, et al., 2015), the effect of chilling ries (site–­species combinations) of leaf-­out dates were analyzed. still needs to be explicitly quantified. If the dominant role of chilling, The number of leaf-­out records, in particular, varied among species, as demonstrated in controlled experiments, also holds true in natu- ranging from 589 to 115,003 (Table 1). ral 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 2.2 | Climatic data if the dominant signal of climate change is increased forcing. Thus, quantifying the relative effects of chilling and forcing has critical The past climatic data of Europe we used were the E-­OBS data forecasting implications. set (Cornes et al., 2018) from the European Climate Assessment & Here, we used long-­term observations (≥20 years) of leaf-­out Dataset (ECA&D) project (https://www.ecad.eu/), where the daily dates derived from 16 dominant tree species in Europe to test the mean temperature during the period 1950–­2020 was available at relative effects of chilling and forcing on spring phenology. We first 1/10° resolution (~10 km). We matched each phenological station in analyzed the relationship between chilling accumulation (CA) and FR a specific grid based on its latitude and longitude. for each time series to examine the effect of chilling on the leaf-­out For the future climatic data, we obtained the daily resolu- date. Subsequently, we developed a phenological model based on tion and bias-­ corrected data simulated by the global climatic the functions between CA and FR to separate the chilling and forc- model HadGEM2-­ES, covering the period 2011–­2099, with a ing contributions to phenological change. Finally, we compared the ­resolution of 0.5 × 0.5° from the Fast Track input-­data catalog of relative magnitudes of the delaying effect of reduced chilling and the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP, WANG et al. | 6035 TA B L E 1 Summary of the species Life Number of stations Number investigated in this study No. Species Common name form (time series) of records 1 Acer platanoides Norway maple DBT 25 589 2 Acer Sycamore maple DBT 33 821 pseudoplatanus 3 Aesculus Horse chestnut DBT 3277 115,003 hippocastanum 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. 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). 2.3 | Relationship between CA and FR For each leaf-­out record, CA was calculated as the accumulated chill- ing 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). ∑ tmedian CA = CU(t), (1) t=t0 ⎧ ⎪ 1 if T(t) < Tc CU(t) = ⎨ , (2) F I G U R E 1 Distribution of the 3754 stations in Europe. These ⎪ 0 else ⎩ 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. where CA is the chilling accumulation for each leaf-­out record; tmedian is the median value of the current leaf-­out time series. t 0 is https://www.isimip.org/). We used future climatic data under two the starting date of chilling accumulation, which could be set as 1 representative concentration pathways (RCP 4.5 and 8.5). Although November of the previous year following previous studies (Laube the ISIMIP data set has corrected the original data (retrieved from et al., 2014; Wang, Wang, et al., 2020). However, a recent warming the CMIP5 archive) for systematic deviations of the simulated his- experiment found a much earlier start date of chilling since October torical data from real observations (Warszawski et al., 2014), the res- (Beil et al., 2021). Thus, we also calculated CA with t 0 = 1 October olution of 0.5° is too coarse for comparisons with the past climatic of the previous year. CU(t) is the chilling unit on date t. T(t) is the 6036 | WANG et al. daily mean temperature (°C) on date t; Tc is the upper-­temperature Renner, 2015). Thus, we fitted the following linear function for each threshold for effective chilling. Many studies have reported that time series of leaf-­out dates. 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 re- FR = f + g × CA, (5) cent study found that the effective chilling temperature could reach up to 10°C (Baumgarten et al., 2021): thus, we also calculated the CA where FR is the forcing requirement for spring events, CA is the chilling with Tc = 10°C. Overall, six chilling algorithms (2 t 0 × 3 Tc) were used accumulation, and f and g are the parameters. to calculate the CA. Many other chilling models also assume that However, several studies described the relationship between temperature below 0°C is ineffective for breaking dormancy, but CA and FR as an exponential function (Cannell & Smith, 1983; Lin these models could not reproduce the negative correlation between et al., 2022; Man et al., 2017). Thus, we further fitted the following CA and FR in the observation data (Wang, Wu, et al., 2020). A recent function: experiment proved that freezing temperature (e.g., −2°C) is also ef- fective in breaking dormancy (Baumgarten et al., 2021). Therefore, FR = a + b × e−c×CA , (6) we did not test chilling models assuming that temperature below 0°C is inefficient in breaking dormancy. where a, b, and c are the parameters. The calculation of FR was based on growing degree-­days For each time series, the parameters of the linear or exponen- (Equation 3), which integrated the forcing units [daily mean tem- tial function were fitted using the least square method. At several peratures above a specific threshold (Tb) from a fixed starting date stations, the FR or CA was very small (FR <10°C or CA <10 days), (t1) to the observed leaf-­out date (Equation 4)]. Tb is often set as 5°C suggesting the species may have adopted a different threshold tem- (Fu et al., 2013; Peaucelle et al., 2019; Wang, Wang, et al., 2020) or perature to start growth or break dormancy due to the adaptation 0°C (Basler & Körner, 2012; Heide, 1993; Piao et al., 2015). t1 is com- to the local climate. We removed these stations for further analysis monly set as 1 January (Fu, Piao, et al., 2015; Peaucelle et al., 2019; because we could not obtain a stable estimate of the parameters for Wang, Wang, et al., 2020), but 1 February has also been reported time series with small CA or FR. by others (Cannell & Smith, 1983; Fu et al., 2013). To test whether For each time series (species−station combination), we simulated the choice of the forcing algorithms would affect the result, we used the leaf-­out date on a daily basis since 1 January. The CA was cal- three forcing algorithms (Tb = 5°C, t1 = 1 January; Tb = 0°C, t1 = 1 culated based on the past (1951–­2014) or future (2022–­2099) cli- January; Tb = 5°C, t1 = 1 February) based on different threshold mate data. The FR for leaf-­out was estimated using the predefined temperatures and starting dates. function between FR and CA (Equations 5 or 6). The leaf-­out date was determined as the day when the accumulated degree-­day began ∑ tLOD FR = FU(t), (3) to exceed the FR. We assessed the performance of the phenolog- t=t1 ical models with two indicators: the coefficient of determination (R 2) and the root-­mean-­square error (RMSE) between the simulated ⎧ ⎪ T(t) − Tb if T(t) < Tb and observed dates. Although the phenological models based on FU(t) = ⎨ , (4) different chilling and forcing algorithms generate interrelated leaf-­ ⎪ 0 else ⎩ out dates (e.g., Figure S1), the RMSE and R 2 between the simulated and observed dates varied greatly among combinations of chilling where FR is the forcing requirement for leaf-­out and FU(t) is the forcing and forcing algorithms (e.g., Figure S2). For all species, the phe- unit at day of year t. tLOD is the observed leaf-­out date, and T(t) is the nological model based on the chilling algorithm with Tc = 5°C and daily mean temperature (°C) on date t. t 0 = 1 November and the forcing algorithm with Tb = 0°C and t1 = 1 To examine how chilling affects the leaf-­out date, we calculated January was most accurate since it showed the lowest RMSE and the Pearson’ r between CA and FR (rCA–­FR) for each time series of highest R 2 between the simulated and the observed leaf-­out date leaf-­out. Because six chilling algorithms and three forcing algorithms among 18 combinations. Thus, we only reported the results based were adopted, 18 rCA–­FR values were calculated for different combi- on this optimal chilling and forcing algorithm combination. nations of chilling and forcing algorithms. 2.5 | Quantifying chilling and forcing contributions 2.4 | Phenological models based on the CA–­FR relationship We developed a novel method to quantify the contribution of chill- ing and forcing on phenological change. This method separated the We calibrated the specific function to establish the relationship be- simulated leaf-­out date into three parts: the leaf-­out date during the tween CA and FR for each time series of leaf-­out dates. According reference period (1961–­1990), the chilling contribution, and the forc- to previous studies, the relationship between CA and FR may be ing contribution relative to the reference period, as expressed in the linear (Nanninga et al., 2017; Wang, Wu, et al., 2020; Zohner & following equation: WANG et al. | 6037 ( ) ( ) Pi FRi = Ci + Fi + Pref FRref , (7) during this period (Equation 8), and the reference leaf-­out date [Pref(FRref)] is calculated as the mean Pi(FR i) during the reference where Pi(FRi) is the simulated leaf-­out date (the date when FRi is ful- period (Equation 9). filled) at a certain year i. Pref(FRref) is the mean leaf-­out date during the ∑1990 reference period (1961–­1990). FRref represents the mean forcing re- i=1961 FRi (8) FRref = , 30 quirement during the reference period. Ci is the chilling contribution in a certain year i, representing the deviation in the leaf-­out date from ∑1990 � � � � i=1961 Pi FRi the reference period induced by changes in the chilling condition; Fi is Pref FRref = . (9) 30 the forcing contribution in a certain year i, representing the deviation in the leaf-­out date from the reference period induced by changes in 3. The forcing contribution (Fi) is calculated as the difference in the the forcing condition. date when FRref is fulfilled in year i and the reference leaf-­out The above quantitative splitting of the predicted phenological date [Pref(FRref )], because the year-­to-­year variation of forcing change into the effects caused by changes in forcing and those temperature alters the date when FRref is fulfilled (Equation 10; caused by changes in chilling is realized by the following steps: Figure 2). ( ) ( ) 1. For each time series of leaf-­out dates, the forcing requirement Fi = Pi FRref − Pref FRref , (10) in a year i (FR i) is calculated based on the amount of chilling (CAi, Equation 1) and the CA–­FR curve (Equations 5 or 6). where Fi and Pref(FRref) describe the same parameters as Subsequently, Pi(FR i), that is, the simulated leaf-­ out date at a Equation (7). Pi(FRref) is calculated as the date when FRref is fulfilled certain year i, is calculated as the date when FRi is fulfilled. in year i. 2. We defined 1961–­1990 as the reference period, and the refer- 4. The chilling contribution (Ci) is calculated as the difference between ence forcing requirement (FRref) is calculated as the mean FR i Pi(FRi) and Pi(FRref) because the year-­to-­year variation in FRi is caused F I G U R E 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 i is smaller than CAref, and then FR i is larger than FRref. (b) The methods to separate the effect of chilling and forcing. Pref(FRref ) is the mean leaf-­out date during the reference period; Pi(FR i) is the simulated leaf-­out date at year i, that is, the date when the accumulated growing degree days (with a threshold Tb, starting date t1) is larger than FR i. Pi(FRref ) is the date in year i when FRref is fulfilled. The chilling contribution (Ci) is calculated as the difference between Pi(FR i) and Pi(FRref ) because the year-­to-­year variation in FR i is caused by the year-­to-­year variation of chilling. The forcing contribution (Fi) is calculated as the difference in Pi(FRref ) and Pref(FRref ) because the year-­to-­year variation of forcing temperature alters the date when FRref is fulfilled. 6038 | WANG et al. by the year-­ to-­ year variation of chilling (Equation 11; Figure 2) if where DL is the day length, which is calculated as a function of lati- any potential contribution of photoperiod is neglected. tude 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 Ci = Pi FRi − Pi FRref , (11) more significant in low chilling conditions than that in high chilling con- where C i and Pi(FR i) describe the same parameters as Equation (7). ditions. Therefore, h and i should vary with CA. Pi(FRref ) describes the same parameters as Equation (10). In order to quantify h and i, we defined a variable w as the ratio Following Equations (10) and (11), Equation (7) is a mathematical of the additional effect of a 16 h photoperiod to the additional effect necessity. of an 8 h photoperiod: Following the methods above, we estimated the interannual changes in chilling and forcing contributions during the past period AE(16) 16h + i w= = . (15) (1951–­2014) and future period (2022–­2099, under RCP 4.5 and 8.5) AE(8) 8h + i for each species−station combination. At the species level, we cal- culated the time series of the chilling contribution, forcing contribu- According to the experimental results of Zohner and Renner (2015), tion, and leaf-­out date averaged from all stations. To obtain a regional under the same chilling conditions, the FR of leaf-­out under 16 h pho- time series, we averaged the time series of the chilling contribution, toperiod was significantly smaller than that under 8 h photoperiod. forcing contribution, and leaf-­out date from all species. The slopes Meanwhile, under the same photoperiod conditions, CA was linearly of the linear regression of the chilling contribution, forcing contribu- and negatively correlated with FR (Figure S3). Thus, w could be es- tion, and leaf-­out date against year were used to measure the trends timated as the ratio of FR under 8 h photoperiod to FR under 16 h of each variable. photoperiod. ⎧ FR8 f + g × CA ⎪ = 8 8 if CA > 123 2.6 | Estimating photoperiod contributions on leaf-­ w = ⎨ FR16 f16 + g16 × CA , (16) ⎪ 1 out date of F. sylvatica ⎩ if CA > 123 In previous experiments, F. sylvatica is identified as a photoperiod-­ where FR8 and FR16 correspond to the FR under an 8 or 16 h photo- sensitive species whose leaf-­o ut date is significantly delayed under period in the experiment, respectively, CA is the chilling accumulation, short photoperiod than long photoperiod (Basler & Körner, 2012; and f8, g8, f16, and g16 are parameters of the linear function. According Fu, Piao, et al., 2019; Zohner & Renner, 2015). In order to quan- to the experimental results of Zohner and Renner (2015), f8 = 2397.0, tify the effect of photoperiod on the leaf-­o ut date of F. sylvatica g8 = −17.14, f16 = 755.9, and g16 = −3.97 (Figure S3). in natural conditions, we further developed a CA-­p hotoperiod The solution of the system of linear equations (Equations 14 and model. In Section 2.4, the phenological model is based on the 15) was linear function between CA and FR (hereinafter called the CA-­ w−1 only model). However, the CA-­p hotoperiod model is based on the h= 4(w + 1) , (17) linear function between CA and FR incorporating the effect of photoperiod (FRP). i=1−3× w−1 . (18) w+1 ∑ tLOD Overall, FRP at each day of the year could be calculated as fol- FRP = FU(t) × AE(t), (12) t=t1 lows: (i) by calculating the CA based on the climate data (Equations 1 and 2 with Tc = 5°C, t 0 = 1 November); (ii) by calculating w based where FRP is the forcing requirement incorporating the effect of pho- on the CA and Equation (16); (iii) by calculating h and i based on w, toperiod and tLOD is the observed leaf-­out date. t1 is the starting date Equations (17) and (18); (iv) by calculating the day length (DL) based of degree-­day accumulation (set as 1 January). FU(t) is the forcing unit on equations from Forsythe et al. (1995); (v) by calculating AE(DL) at the day of year t (Equation 4, Tb is set as 0°C). AE represents the based on h, i, DL and Equation (13); (vi) by calculating FRP based on additional effect of photoperiod on leaf-­out. Equation (12). We assumed that the AE is linearly correlated with day length. At each station of F. sylvatica, we fitted the linear function be- A long photoperiod accelerates the accumulation of the FR (higher tween the CA and FRP. Likewise the CA-­only model, we simulated AE), while, on the contrary, a short photoperiod suppresses it. Thus, the leaf-­out date on a daily basis since 1 January. The CA was calcu- AE >1 if day length >12 h, AE = 1 if day length = 12 h, and AE <1 if lated based on the past (1951–­2014) or future climate data (2022–­ day length <12 h. 2099). FRP for the leaf-­out was estimated using the predefined function between FRP and CA. The day when the accumulated AE(DL) = h × DL + i, (13) photoperiod-­associated degree-­day began to exceed FRP was es- timated as the leaf-­out date. The contribution of photoperiod on AE(12) = h × 12 + i = 1, (14) the leaf-­out date of F. sylvatica is estimated as the difference in the WANG et al. | 6039 simulated leaf-­out date between the CA-­photoperiod model and February, Figure 3b). The CA thus decreased significantly by the CA-­only model. 3.86 days/decade (p < .01) in Europe (Figure 3a). We further calculated the FR [based on Equations (3) and (4) with Tb = 5°C and t1 = 1 January] and CA [based on Equations (1) and (2) 3 | R E S U LT S with Tc = 5°C and t0 = 1 November] for each record of the leaf-­out date. When analyzing the Pearson's r between CA and FR (rCA–­FR) for each 3.1 | Delaying effect of decreased chilling on the time series, we found that the FR was negatively correlated with the leaf-­out date CA in most cases (see the examples in Figure 3c). The proportion of the time series with negative rCA–­FR varied among species, ranging from The winter temperature (previous November to February) 93.8% to 100% (37.5%–­95.8% at p < .05, Figure 3d). We further used averaged from all phenological stations increased by 0.24°C/ other five chilling and two forcing algorithms to calculate the CA and decade (p < .01) from 1951 to 2014, respectively (Figure 3a). FR. The results showed that the prevalent negative rCA–­FR was indepen- At the interannual scale, the winter temperature correlated dent of the algorithms applied to calculate chilling or forcing (Figure S4). negatively with the CA (calculated as the number of days with Thus, over the past six decades, the signal of the delaying effect of de- a daily mean temperature <5°C during previous November to creased chilling on the leaf-­out date was shown to be strong. F I G U R E 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 (rCA-­HR) for each species. The number of time series (i.e., stations) for each species is shown in parentheses. 6040 | WANG et al. 3.2 | Contribution of chilling and forcing on past among species (Figure 5c; Figure S6), ranging from 0.09 to 0.18 days/ phenological change year. For each specific species, the trends of the chilling contribution also differed among stations with an interquartile range of 0.06–­ We developed a phenological model for each time series of leaf-­out 0.09 days/year (Figure 5b). dates using the linear function between CA and FR (Equation 5) based In order to examine to what degree the delaying effect of chill- on the optimal chilling and forcing algorithms. The leaf-­out dates could ing could offset the advancing effect of forcing, we estimated the be simulated as the dates when the CA-­associated FR was satisfied. forcing contribution as the deviation in the leaf-­out date from the The simulated and observed leaf-­out dates corresponded well given reference period (1961–­1990) induced by changes in the forcing the general uncertainty in phenological modeling (Figure 4). The RMSE temperature. The forcing contribution, averaged from all the spe- 2 and the coefficient of determination (R ) of all the records ranged from cies, enhanced during 1951–­2014 with a trend of −0.27 days/year 6.27 to 9.81 days and 0.40 to 0.81, respectively. Furthermore, our (p < .01). The forcing contributions were negative in almost all years model was able to reproduce the past spatial pattern of leaf-­out dates after 1988 (except for 1996, 2006, and 2013), indicating that the when comparing the simulated and observed leaf-­out dates at the sta- recent spring warming had an advancing effect on the leaf-­out date tion level (Figure S5). Thus, this model could be used to separate the of temperate species (Figure 5a). contribution of chilling and forcing on phenological change. The magnitude of the increasing trends of the chilling contribu- We estimated the chilling contribution to the past phenological tion was weaker than the decreasing trends of the forcing contribu- change (1951–­2014) as the deviation in the leaf-­out date from the tion (Figure 5a), and this result was also robust at the species level reference period (1961–­1990) induced by changes in the chilling. (Figure 5c) and station level (Figure 5b). Therefore, the simulated The chilling contribution, averaged from all the species, showed a leaf-­out date became earlier by 0.12 days/year (p < .01, Figure 5a) positive trend of 0.15 days/year (p < .01), suggesting an enhanced from 1951 to 2014, which matched previous results based on long-­ delaying effect of chilling caused by winter warming (Figure 5a). The term observations (Cayan et al., 2001; Ge et al., 2015; Menzel chilling contributions were positive in almost all years after 1988 et al., 2006). Furthermore, the chilling contribution correlated sig- (except in 1996 and 2006). Especially in 2007, the insufficient chill- nificantly with the forcing contribution at the interannual scale ing caused a delaying effect of 20.7 days relative to the reference (Figure 5d). The slope of the forcing contribution against the chilling period (1961–­1990). The trends of the chilling contribution varied contribution indicated that 1 day delay in the leaf-­out date induced F I G U R E 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. WANG et al. | 6041 F I G U R E 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). by chilling corresponded to 1.91 days advance induced by forcing the exponential function-­based model were similar to the linear (Figure 5d). Thus, the winter warming-­related reductions in chilling function-­based model (Figure 7). 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). 3.3 | Future changes in leaf-­out dates We further developed a phenological model based on the ex- ponential function between CA and FR (Equation 6). Similar to the We first predicted the future changes in the leaf-­out dates of temper- linear function-­based model, the RMSE and R 2 of the exponential ate species from 2022 to 2099 under RCP 4.5 and RCP 8.5 scenarios function-­based model ranged from 6.26 to 9.84 days and 0.39 to using the linear function-­based model. In both scenarios, the leaf-­out 0.95, respectively (Figure S7). Moreover, the exponential function-­ date would continue to become earlier. The trends of the leaf-­out based model also could accurately reproduce the past spatial pat- dates were −0.016 days/year (not significant) under RCP 4.5 averaged tern of leaf-­out dates when comparing the simulated and observed from 16 tree species (Figure 8a). Under RCP 8.5, the trends increased leaf-­out dates at the station level (Figure S5). Thus, the exponen- to −0.166 days/year (p < .01; Figure 8b). The earlier leaf-­ out date tial function-­ based model had a similar accuracy to the linear could be attributed to the weaker trends of the chilling contribution function-­based model. Furthermore, the trend and amplitude of the than the forcing contribution (Figure 8a,b). The future trends in the chilling and forcing contributions from 1951 to 2014 derived from leaf-­out date were stronger under RCP 8.5 with intensified warming 6042 | WANG et al. F I G U R E 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. than under RCP 4.5, but the trends varied among species (Figures S8 the future. The reason is that the delaying effect of decreased chill- and S9). For future scenarios (2022–­2099), the interannual chilling ing simulated by the exponential function is stronger than that sim- contribution (averaged from all the species) correlated significantly ulated by the linear function, and the advancing effect of increased with the forcing contribution (Figure 8c). Under RCP 4.5, a 1 day de- forcing is the same between the two types of function (Figure 9c). laying 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 3.4 | Contribution of photoperiod to leaf-­out used and a linear, stable correlation between CA and FR is assumed dates of European beech in the future warming scenarios, it is likely that the delaying effect of chilling will substantially offset the advance of leaf unfolding to a We specially developed a CA-­photoperiod model (incorporating similar degree in comparison to the past period. the effect of photoperiod) for simulating the leaf-­out date of F. syl- When using the exponential function between CA and FR to sim- vatica. The photoperiod contribution to the leaf-­out date of F. syl- ulate the leaf-­out date, the trends of the leaf-­out dates averaged from vatica could be estimated as the difference in the simulated leaf-­out 16 tree species were 0.026 and 0.016 days/year (both not significant) date between the CA-­photoperiod model and the CA-­only model. under RCP 4.5 and RCP 8.5, respectively (Figure 9a,b). Thus, if the The CA-­photoperiod model could slightly improve the simulation of exponential function between CA and FR is assumed in the future leaf-­out date compared to the original model (Figure S10), and pre- warming scenarios, the spring phenological advance may reverse in dicted similar trends in the leaf-­out date with the CA-­only model WANG et al. | 6043 F I G U R E 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). (Figure 10a,c,e). At the interannual scale, the photoperiod contribu- leaf-­out in more than 99% of the time series. A global meta-­analysis tion was positive when the leaf-­out date was earlier, but was nega- of all published experiments found that the effect of chilling on the tive when the leaf-­out date was later (Figure 10b,d,f). However, the spring phenology of species was stronger than the effect of forcing amplitude of the photoperiod contribution among years was only (Ettinger et al., 2020). However, this study suggested that the delaying ±6 days. 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 4 | DISCUSSION the wider range of chilling treatments in the experiments. For example, the difference in the CA under low chilling and high chilling As results from the literature show, the main evidence that a treatments reached 77 days in an experiment focused on 36 woody reduction in chilling could delay the leaf-­out date has come from the species in Europe (Laube et al., 2014). However, the CA in the field controlled experiment on saplings or twigs of temperate species (Du only reduced by 24.3 days during the period 1951–­2014 in Europe et al., 2019; Flynn & Wolkovich, 2018; Man et al., 2017; Nanninga (Figure 3a). Therefore, a stronger response of spring phenology to et al., 2017). This study provided empirical evidence in natural chilling is found if experiments expand observed chilling beyond the conditions since the CA was negatively correlated to the FR of current conditions. In current climate conditions, the counteracting 6044 | WANG et al. F I G U R E 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). effects of chilling are insufficient to stop the advance of leaf increased exponentially with the decrease in CA (e.g., Figure 3c), unfolding in response to warmer forcing temperatures. and the effect of chilling would be amplified in an extremely This study used two hypothesised CA–­FR functions to sim- warm winter. In the past period (1951–­2 014), the trend and am- ulate the leaf-­o ut date of temperate species. Compared to the plitude of the chilling and forcing contributions derived from linear function, the exponential function assumes that the FR the exponential function-­b ased model were similar to the linear WANG et al. | 6045 F I G U R E 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). function-­b ased model. However, in future warming scenarios, hypothesis of an exponential function between CA and FR, the the delaying effect of decreased chilling become stronger and future advance in the leaf-­o ut date may stop for tree species even exceed the advancing effect of increased forcing in sev- with the relatively high chilling requirement. We need to further eral extremely warm winters (Figure 9c). Therefore, under the validate the forms of relationship between CA and FR since the 6046 | WANG et al. F I G U R E 1 0 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. linear and exponential functions were shown to both accurately Zhang, et al., 2019). In warm springs, when leaf-­out is early, the short describe the phenological change that happened in the past, but day length mitigates the warming-­induced advancement of leaf-­out. to differ in performance regarding future warming scenarios. In cold springs, when leaf-­out is late, the long day length promotes Besides forcing and chilling, photoperiod has been reported to be a leaf-­out, ensuring that trees do not leaf-­out too late. However, the am- notable factor affecting spring phenology (Way & Montgomery, 2015), plitude of the photoperiod contribution was significantly smaller than although only 35% of species relied on spring photoperiod as a leaf-­ that of the chilling contribution (±40 days). Therefore, decreased chill- out signal in low chilling conditions (Zohner et al., 2016). We selected ing plays a dominant role in mitigating the advance of leaf-­out dates. F. sylvatica, a well-­known photoperiod-­sensitive species (Basler & The phenological models developed in this study, which describe Körner, 2012; Fu, Piao, et al., 2019; Zohner & Renner, 2015), to inves- the response of spring phenology of temperate tree species to warm- tigate the relative contribution of chilling and photoperiod on pheno- ing, are valid within the temperature range we have observed during logical change. The results confirmed the hypothesis that day length the past several decades (1951–­2014). The CA–­FR functions we es- helps temperate deciduous trees to leaf-­out at the optimal time (Fu, tablished might not be stable if chilling declines beyond the range of WANG et al. | 6047 current temperature conditions. The plants might be adapted only to REFERENCES a certain amount of chilling, and might respond unexpectedly (likely Basler, D., & Körner, C. (2012). Photoperiod sensitivity of bud burst stronger) if chilling becomes much more incomplete. Therefore, pro- in 14 temperate forest tree species. Agricultural and Forest Meteorology, 165, 73–­81. https://doi.org/10.1016/j.agrfo​r met.​ jected changes in future spring phenology in this study are based on 2012.​0 6.​0 01 the assumption that the CA–­FR function observed in the past is stable Baumgarten, F., Zohner, C. M., Gessler, A., & Vitasse, Y. (2021). Chilled to in the future and may overestimate the spring phenological advance. be forced: The best dose to wake up buds from winter dormancy. Overall, using 27,071 time series of leaf-­out dates for 16 temperate New Phytologist, 230(4), 1366–­1377. https://doi.org/10.1111/ nph.17270 tree species in Europe, we showed that the delaying effect of chilling Beil, I., Kreyling, J., Meyer, C., Lemcke, N., & Malyshev, A. V. (2021). Late on the leaf-­out date found in experiments was prevalent in natural to bed, late to rise—­Warmer autumn temperatures delay spring conditions, as more than 99% of time series exhibited a negative re- phenology by delaying dormancy. Global Change Biology, 27(22), lationship between CA and FR. Based on the linear function between 5806–­5817. https://doi.org/10.1111/gcb.15858 Cannell, M. G. R., & Smith, R. I. (1983). Thermal time, chill days and pre- CA and FR, we first quantified the relative contribution of chilling and diction of budburst in Picea sitchensis. Journal of Applied Ecology, forcing on phenological change. The winter warming-­related reduc- 20(1), 951–­963. https://doi.org/10.2307/2403139 tions in chilling during the time span 1951–­2014 could offset about Cayan, D. R., Kammerdiener, S. A., Dettinger, M. D., Caprio, J. M., & 52% of spring phenological advances caused by increased forcing. Peterson, D. H. (2001). Changes in the onset of spring in the Assuming that the linear CA–­FR function is stable in the future warm- western United States. Bulletin of the American Meteorological Society, 82(3), 399–­416. https://doi.org/10.1175/1520-­0 477(2001)​ ing scenarios, the delaying effect of chilling would be enhanced in the 082<0399:CITOO​S>2.3.CO;2 future and slow down the spring phenological advance to a similar de- Cleland, E. E., Allen, J. M., Crimmins, T. M., Dunne, J. A., Pau, S., Travers, gree. However, it is possible that the relationship between FR and CA S. E., Zavaleta, E. S., & Wolkovich, E. M. (2012). Phenological track- is or would become exponential (maybe in particular if temperature ing enables positive species responses to climate change. Ecology, 93(8), 1765–­1771. https://doi.org/10.1890/11-­1912.1 rises), and then the contribution of chilling could fully compensate for Cornes, R. C., van der Schrier, G., van den Besselaar, E. J. M., & Jones, forcing in extremely warm years. More controlled experiments and P. D. (2018). An ensemble version of the E-­O BS temperature continuous monitoring of spring phenology are needed to investigate and precipitation data sets. Journal of Geophysical Research: whether a linear CA–­FR relationship holds in a broader range of tem- Atmospheres, 123(17), 9391–­9409. https://doi.org/10.1029/​ 2017J​D 028200 peratures than the one we can observe by now, or whether the rela- Du, Y., Pan, Y., & Ma, K. (2019). Moderate chilling requirement con- tionship appears to be exponential if temperatures rise. trols budburst for subtropical species in China. Agricultural and Forest Meteorology, 278, 107693. https://doi.org/10.1016/j.agrfo​ AC K N OW L E D G M E N T S rmet.2019.107693 Elmendorf, S. C., & Ettinger, A. K. (2020). Is photoperiod a domi- This work was funded by the National Key R&D Program of China nant driver of secondary growth resumption? Proceedings of the (Grant No. 2018YFA0606102), National Natural Science Foundation National Academy of Sciences, 117(52), 32861–­32864. https://doi. of China (Grant No. 41871032, 42125101), Youth Innovation org/10.1073/pnas.20194​11117 Promotion Association, CAS (Grant No. 2018070), and Program for Ettinger, A. K., Chamberlain, C. J., Morales-­C astilla, I., Buonaiuto, D. M., Flynn, D. F. B., Savas, T., Samaha, J. A., & Wolkovich, E. M. “Kezhen” Excellent Talents in IGSNRR, CAS (Grant No. 2018RC101). (2020). Winter temperatures predominate in spring phenological responses to warming. Nature Climate Change, 10(12), 1119–­1137. C O N FL I C T S O F I N T E R E S T https://doi.org/10.1038/s4155​8-­020-­0 0917​-­3 The authors declare no known conflicts of interest. Flynn, D. F. B., & Wolkovich, E. M. (2018). Temperature and photope- riod drive spring phenology across all species in a temperate for- est community. New Phytologist, 219(4), 1353–­1362. https://doi. DATA AVA I L A B I L I T Y S TAT E M E N T org/10.1111/nph.15232 The derived data that support the results and figures in this study Forsythe, W. C., Rykiel, E. J., Stahl, R. S., Wu, H. I., & Schoolfield, R. M. are openly available in Dryad (https://doi.org/10.5061/dryad.83bk3​ (1995). A model comparison for daylength as a function of latitude j9v6). Phenological data in Europe can be accessed from http:// and day of year. Ecological Modelling, 80(1), 87–­95. https://doi. org/10.1016/0304-­3800(94)00034​-­F www.pep725.eu/. E-­OBS gridded climatic dataset (24.0e) from 1950 Fu, Y. H., Campioli, M., Deckmyn, G., & Janssens, I. A. (2013). to 2014 in Europe can be accessed from https://www.ecad.eu/. Sensitivity of leaf unfolding to experimental warming in Future climate data (ISIMIP2b) were available on the ISIMIP ESGF three temperate tree species. Agricultural and Forest server (https://data.isimip.org/). Meteorology, 181, 125–­132. https://doi.org/10.1016/j.agrfo​ rmet.​2 013.​07.016 Fu, Y. H., Piao, S., Vitasse, Y., Zhao, H., De Boeck, H. J., Liu, Q., Yang, ORCID H., Weber, U., Hanninen, H., & Janssens, I. A. (2015). Increased Huanjiong Wang https://orcid.org/0000-0002-2325-0120 heat requirement for leaf flushing in temperate woody species Junhu Dai https://orcid.org/0000-0001-6401-1735 over 1980–­2012: Effects of chilling, precipitation and insolation. Global Change Biology, 21(7), 2687–­2697. https://doi.org/10.1111/ Josep Peñuelas https://orcid.org/0000-0002-7215-0150 gcb.12863 Quansheng Ge https://orcid.org/0000-0001-8712-8565 Fu, Y. H., Piao, S., Zhou, X., Geng, X., Hao, F., Vitasse, Y., & Janssens, I. Yongshuo H. Fu https://orcid.org/0000-0002-9761-5292 A. (2019). Short photoperiod reduces the temperature sensitivity Chaoyang Wu https://orcid.org/0000-0001-6163-8209 of leaf-­out in saplings of Fagus sylvatica but not in horse chestnut. 6048 | WANG et al. Global Change Biology, 25(5), 1696–­1703. https://doi.org/10.1111/ Peñuelas, J., Rutishauser, T., & Filella, I. (2009). Phenology feedbacks gcb.14599 on climate change. Science, 324(5929), 887–­888. https://doi. Fu, Y. H., Zhang, X., Piao, S., Hao, F., Geng, X., Vitasse, Y., Zohner, org/10.1126/scien​ce.1173004 C., Penuelas, J., & Janssens, I. A. (2019). Daylength helps tem- Piao, S., Tan, J., Chen, A., Fu, Y. H., Ciais, P., Liu, Q., Janssens, I. A., Vicca, perate deciduous trees to leaf-­o ut at the optimal time. Global S., Zeng, Z., Jeong, S., Li, Y., Myneni, R. B., Peng, S., Shen, M., & Change Biology, 25(7), 2410–­2418. https://doi.org/10.1111/ Penuelas, J. (2015). Leaf onset in the northern hemisphere trig- gcb.14633 gered by daytime temperature. Nature Communications, 6, 6911. Fu, Y. H., Zhao, H., Piao, S., Peaucelle, M., Peng, S., Zhou, G., Ciais, P., https://doi.org/10.1038/ncomm​s7911 Huang, M., Menzel, A., Penuelas, J., Song, Y., Vitasse, Y., Zeng, Z., Richardson, A. D., Andy Black, T., Ciais, P., Delbart, N., Friedl, M. & Janssens, I. A. (2015). Declining global warming effects on the A., Gobron, N., Hollinger, D. Y., Kutsch, W. L., Longdoz, B., phenology of spring leaf unfolding. Nature, 526(7571), 104–­107. Luyssaert, S., Migliavacca, M., Montagnani, L., William Munger, https://doi.org/10.1038/natur​e15402 J., Moors, E., Piao, S., Rebmann, C., Reichstein, M., Saigusa, N., Ge, Q., Wang, H., Rutishauser, T., & Dai, J. (2015). Phenological response Tomelleri, E., … Varlagin, A. (2010). Influence of spring and au- to climate change in China: A meta-­analysis. Global Change Biology, tumn phenological transitions on forest ecosystem productiv- 21(1), 265–­274. https://doi.org/10.1111/gcb.12648 ity. Philosophical Transactions of the Royal Society B: Biological Hänninen, H., Kramer, K., Tanino, K., Zhang, R., Wu, J., & Fu, Y. H. Sciences, 365(1555), 3227–­3246. https://doi.org/10.1098/ (2019). Experiments are necessary in process-­based tree phenol- rstb.2010.0102 ogy modelling. Trends in Plant Science, 24(3), 199–­209. https://doi. Templ, B., Koch, E., Bolmgren, K., Ungersboeck, M., Paul, A., org/10.1016/j.tplan​t s.2018.11.006 Scheifinger, H., Rutishauser, T., Busto, M., Chmielewski, F., Heide, O. M. (1993). Daylength and thermal time responses of bud- Hajkova, L., Hodzic, S., Kaspar, F., Pietragalla, B., Romero-­ burst during dormancy release in some northern deciduous trees. Fresneda, R., Tolvanen, A., Vucetic, V., Zimmermann, K., & Zust, Physiologia Plantarum, 88(4), 531–­540. https://doi.org/10.1111/ A. (2018). Pan European Phenological database (PEP725): A j.1399-­3 054.1993.tb013​68.x single point of access for European data. International Journal Hiestand, M. P., & Carleton, A. M. (2019). Growing-­season synoptic and of Biometeorology, 62(6), 1109–­1113. https://doi.org/10.1007/ phenological controls on heat fluxes over forest and cropland sites s0048​4 -­018-­1512-­8 in the Midwest U.S. Corn Belt. Journal of Applied Meteorology and Wang, H., Wang, H., Ge, Q., & Dai, J. (2020). The interactive effects of Climatology, 59(3), 381–­400. https://doi.org/10.1175/JAMC-D-​19- chilling, photoperiod, and forcing temperature on flowering phe- 0019.1 nology of temperate woody plants. Frontiers in Plant Science, 11, Keenan, T. F., Gray, J., Friedl, M. A., Toomey, M., Bohrer, G., Hollinger, D. 443. https://doi.org/10.3389/fpls.2020.00443 Y., Munger, J. W., O Keefe, J., Schmid, H. P., Wing, I. S., Yang, B., & Wang, H., Wu, C., Ciais, P., Peñuelas, J., Dai, J., Fu, Y., & Ge, Q. (2020). Richardson, A. D. (2014). Net carbon uptake has increased through Overestimation of the effect of climatic warming on spring phenol- warming-­induced changes in temperate forest phenology. Nature ogy due to misrepresentation of chilling. Nature Communications, Climate Change, 4(7), 598–­604. https://doi.org/10.1038/NCLIM​ 11(1), 4945. https://doi.org/10.1038/s4146​7-­020-­18743​-­8 ATE2253 Warszawski, L., Frieler, K., Huber, V., Piontek, F., Serdeczny, O., & Schewe, Laube, J., Sparks, T. H., Estrella, N., Hoefler, J., Ankerst, D. P., & Menzel, J. (2014). The inter-­s ectoral impact model Intercomparison proj- A. (2014). Chilling outweighs photoperiod in preventing precocious ect (ISI-­MIP): Project framework. Proceedings of the National spring development. Global Change Biology, 20(1), 170–­182. https:// Academy of Sciences of the United States of America, 111(9), 3228–­ doi.org/10.1111/gcb.12360 3232. https://doi.org/10.1073/pnas.13123​3 0110 Lin, S., Wang, H., Ge, Q., & Hu, Z. (2022). Effects of chilling on heat re- Way, D. A., & Montgomery, R. A. (2015). Photoperiod constraints on tree quirement of spring phenology vary between years. Agricultural and phenology, performance and migration in a warming world. Plant Forest Meteorology, 312, 108718. https://doi.org/10.1016/j.agrfo​ Cell and Environment, 38(9), 1725–­1736. https://doi.org/10.1111/ rmet.2021.108718 pce.12431 Man, R., Lu, P., & Dang, Q. (2017). Insufficient chilling effects vary among Wolkovich, E. M., Auerbach, J., Chamberlain, C. J., Buonaiuto, D. M., boreal tree species and chilling duration. Frontiers in Plant Science, Ettinger, A. K., Morales-­C astilla, I., & Gelman, A. (2021). A simple 8, 1354. https://doi.org/10.3389/fpls.2017.01354 explanation for declining temperature sensitivity with warming. Meier, U. (2001). Growth stages of mono-­ and dicotyledonous plants: Global Change Biology, 27(20), 4947–­4949. https://doi.org/10.1111/ BBCH monograph (2nd ed.). Federal Biological Research Centre for gcb.15746 Agriculture and Forestry. Wolkovich, E. M., & Ettinger, A. K. (2014). Back to the future for plant Menzel, A., Sparks, T. H., Estrella, N., Koch, E., Aasa, A., Ahas, R., Alm-­ phenology research. New Phytologist, 203(4), 1021–­1024. https:// Kübler, K., Bissolli, P., Braslavská, O., Briede, A., Chmielewski, doi.org/10.1111/nph.12957 F. M., Crepinsek, Z., Curnel, Y., Dahl, Å., Defila, C., Donnelly, A., Xu, X., Riley, W. J., Koven, C. D., Jia, G., & Zhang, X. (2020). Earlier leaf-­ Filella, Y., Jatczak, K., Måge, F., … Zust, A. (2006). European phe- out warms air in the north. Nature Climate Change, 10(4), 370–­375. nological response to climate change matches the warming https://doi.org/10.1038/s4155​8-­020-­0713-­4 pattern. Global Change Biology, 12(10), 1969–­1976. https://doi. Zha, T., Barr, A. G., van der Kamp, G., Black, T. A., McCaughey, J. H., & org/10.1111/j.1365-­2486.2006.01193.x Flanagan, L. B. (2010). Interannual variation of evapotranspiration Nanninga, C., Buyarski, C. R., Pretorius, A. M., & Montgomery, R. A. from forest and grassland ecosystems in western Canada in rela- (2017). Increased exposure to chilling advances the time to bud- tion to drought. Agricultural and Forest Meteorology, 150(11), 1476–­ burst in north American tree species. Tree Physiology, 37(12), 1727–­ 1484. https://doi.org/10.1016/j.agrfo​rmet.2010.08.003 1738. https://doi.org/10.1093/treep​hys/tpx136 Zohner, C. M., Benito, B. M., Svenning, J., & Renner, S. S. (2016). Day Peaucelle, M., Janssens, I. A., Stocker, B. D., Descals Ferrando, A., length unlikely to constrain climate-­driven shifts in leaf-­out times Fu, Y. H., Molowny-­H oras, R., Ciais, P., & Penuelas, J. (2019). of northern woody plants. Nature Climate Change, 6(12), 1120–­ Spatial variance of spring phenology in temperate deciduous 1123. https://doi.org/10.1038/NCLIM​ATE3138 forests is constrained by background climatic conditions. Nature Zohner, C. M., Mo, L., & Renner, S. S. (2018). Global warming reduces Communications, 10, 5388. https://doi.org/10.1038/s4146​7-­019-­ leaf-­out and flowering synchrony among individuals. eLife, 7, 13365​-­1 e40214. https://doi.org/10.7554/eLife.40214 WANG et al. | 6049 Zohner, C. M., & Renner, S. S. (2015). Perception of photoperiod in in- dividual buds of mature trees regulates leaf-­out. New Phytologist, How to cite this article: Wang, H., Dai, J., Peñuelas, J., Ge, Q., 208(4), 1023–­1030. https://doi.org/10.1111/nph.13510 Fu, Y. H., & Wu, C. (2022). Winter warming offsets one half of the spring warming effects on leaf unfolding. Global Change Biology, 28, 6033–6049. https://doi.org/10.1111/gcb.16358 S U P P O R T I N G I N FO R M AT I O N Additional supporting information can be found online in the Supporting Information section at the end of this article.