ARTICLE https://doi.org/10.1038/s41467-020-18743-8 OPEN Overestimation of the effect of climatic warming on spring phenology due to misrepresentation of chilling Huanjiong Wang 1,2 ✉, Chaoyang Wu 1,2 ✉, Philippe Ciais 3, Josep Peñuelas 4,5, Junhu Dai1,2, Yongshuo Fu6 & Quansheng Ge 1,2 ✉ 1234567890():,; Spring warming substantially advances leaf unfolding and flowering time for perennials. Winter warming, however, decreases chilling accumulation (CA), which increases the heat requirement (HR) and acts to delay spring phenology. Whether or not this negative CA-HR relationship is correctly interpreted in ecosystem models remains unknown. Using leaf unfolding and flowering data for 30 perennials in Europe, here we show that more than half (7 of 12) of current chilling models are invalid since they show a positive CA-HR relationship. The possible reason is that they overlook the effect of freezing temperature on dormancy release. Overestimation of the advance in spring phenology by the end of this century by these invalid chilling models could be as large as 7.6 and 20.0 days under RCPs 4.5 and 8.5, respectively. Our results highlight the need for a better representation of chilling for the correct understanding of spring phenological responses to future climate change. 1 Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China. 2 University of Chinese Academy of Sciences, Beijing 100049, China. 3 Laboratoire des Sciences du Climat et de l’Environnement, IPSL- LSCE CEA CNRS UVSQ, 91191 Gif sur Yvette, France. 4 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona 08193 Catalonia, Spain. 5 CREAF, Cerdanyola del Valles, Barcelona 08193 Catalonia, Spain. 6 College of Water Sciences, Beijing Normal University, Beijing 100875, China. ✉email: wanghj@igsnrr.ac.cn; wucy@igsnrr.ac.cn; geqs@igsnrr.ac.cn NATURE COMMUNICATIONS | (2020)11:4945 | https://doi.org/10.1038/s41467-020-18743-8 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18743-8 S pring phenology (e.g., budburst, leaf-out, and flowering) of Results perennials has advanced in recent decades around the Linear trends in CA. We assessed 12 chilling models coded as middle and high latitudes (north of 30° N) of the Northern C1–C12 (see “Methods” for a full description). The linear trends in Hemisphere1–3. Shifts in spring phenology have essential impli- CA during winter to early spring (previous November to April) cations for ecosystems and the climatic system, ranging from from 1951 to 2018 were calculated for each station based on all interactions among plants and their herbivores4,5 to the balance models and E-OBS gridded temperature data set. We expected among surface carbon, energy, and water6–8. Projecting future that the accumulated chilling in Central Europe has progressively phenological changes are thus increasingly important for asses- decreased as the winter temperature (Previous November to sing the impact of climate change on ecosystems. Phenological February) increased significantly by 0.25 °C/decade from 1951 to prediction, however, still has large uncertainties, because climatic 2018. The trends in CA and its spatial pattern, however, varied warming has dual effects on spring phenology9. Spring warming among the models (Fig. 2). CA for Models C1, C2, C4, C5, and C12 (forcing temperature) advances spring phenological events10, but decreased significantly at >88% of the stations (p < 0.05) because this effect is counteracted by the reduced chilling due to winter the winter temperature exceeded the maximum effective tem- warming11,12. Several forcing and chilling models have been perature under climate warming (e.g., 5 °C for model C1). The developed based on controlled experiments or inverse modeling proportion of stations reporting a significant decrease in CA for (calibrated the model with the observation data) to quantitatively the other models only ranged from 1.0% (Model C8) to 61.5% describe the effects of temperature on spring phenology9. A (Model C3). The increasing trends even dominated across Central widely used forcing model based on growing degree days Europe for models C7–C9, possibly because the freezing tem- (GDDs)13 calculates the accumulated number of degrees above a perature became effective under climate warming in these models, particular temperature threshold after a specific date. The gen- and this effect of the increase was higher than the effect of erality of chilling models, however, limits operational applic- decrease causing by the loss of maximum effective temperature. ability, because they were developed for specific species and Pearson’s r between each pair of chilling models for each station climatic zones14. (Supplementary Fig. 4) indicated that not all chilling models were The dual role of temperature in regulating the spring phenol- positively correlated with each other. Some models were even ogy of perennials can be described by a negative relationship negatively interrelated (e.g., Model C1 and C8), suggesting that the between the heat requirement (HR, the accumulated forcing current chilling models could not reflect a consistent trend when temperature required for a phenological event) and the chilling measuring the change in the amount of chilling. accumulation (CA, the amount of chilling received by plants during endodormancy) (see the references in Supplementary Table 1). Experimental evidence suggests that saplings and twigs Relationship between CA and HR. We calculated CA (from 1 need less time and fewer degree days to budburst if pre-burst November in the previous year to the date of onset of spring chilling in natural or controlled environments is longer15,16. events) for each species, station, and year using the 12 chilling Phenological observations also support this phenomenon when models to determine if the models met the physiological comparing HR with CA (based on several chilling models) in assumption of Fig. 1. Also, HR was calculated based on a com- time and space12,17,18. From a detailed literature survey, we found monly used forcing model (integrating daily mean temperatures that among 95 perennials investigated, 91 (95.8%) belonging to 46 >0 °C from 1 January to the date of onset of spring events). The genera followed this relationship (Supplementary Fig. 1, Supple- correlations between CA and HR varied among the chilling mentary Table 2). These species involve many life forms (herbs, models for all species. Figure 3 shows an example of leaf-out for deciduous or evergreen coniferous trees, and deciduous broad- Betula pendula. Five chilling models (C1, C2, C4, C5, and C12) that leaved trees or shrubs), and include the most dominant forest tree predicted decreasing trends in CA exhibited significant negative species in Europe (Supplementary Fig. 2). The negative correla- correlations between CA and HR, and the other models exhibited tion between CA and HR is therefore common among perennials positive correlations. Twenty-nine of the 30 species had sig- and stable across studies, species, and locations (Fig. 1). The nificantly negative Pearson’s r between HR and CA based on the molecular basis of the control of budburst mediated by chilling above five models (Supplementary Fig. 5). We therefore classified was recently identified for aspen19. Chilling temperatures down- models C1, C2, C4, C5, and C12 as valid models, because they met regulat the SVL and TCP18 genes, which are negative regulators the physiological assumption in Fig. 1. No or only a few species of budburst, and promote budburst under subsequent forcing had negative CA–HR relationships in the other seven models conditions19. (Supplementary Fig. 5), which were thus classified as invalid Many chilling models have been used to measure the amount models. We also analyzed the station-level relationship between of chilling by quantifying the rate of chilling in daily (or hourly) CA and HR for all stations with at least 15-years of records. Few temperatures20,21, but whether or not current chilling models can stations (<10% for all species) had significantly negative Pearson’s describe this physiological process is unclear. We use long-term r between CA and HR when using the invalid models, but >30% in situ observations of leaf unfolding or flowering for 30 per- of the stations (for most of the species) had significantly negative ennials at 15,533 phenological stations across central Europe (see r between CA and HR when using the valid models (Supple- the distribution of the stations in Supplementary Fig. 3 and the mentary Fig. 6). Since we only used one forcing model (Model F1, list of species in Supplementary Table 3) to compare long-term see the equations in “Methods”) to calculate HR, we further trends in CA and the relationship between CA and HR based on validated the chilling models by correlating them with HR based different chilling models. We aim to assess the ability of current on the other seven forcing models (Models F2–F8). The percen- chilling models to accurately simulate physiological processes and tage of stations with significantly negative Pearson’s r between CA estimate uncertainties for predicting future changes in spring and HR was larger for the valid than the invalid models, phenology. We show that 7 of 12 current chilling models fail to regardless of the forcing model chosen (Supplementary Fig. 7). account for the correct relationship between CA and HR, leading Thus, the performance of the chilling models was not affected by to substantial overestimates of the advance of spring phenology the choice of forcing models. Although some of the invalid under climate change. There is an urgent need to address the models were calibrated specifically for fruit and nut trees (e.g., representation of phenological modeling, in this case specifically Model C7 and C9), they still did not perform better than the valid chilling models. models for Prunus avium (a fruit tree supplying edible cherry) 2 NATURE COMMUNICATIONS | (2020)11:4945 | https://doi.org/10.1038/s41467-020-18743-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18743-8 ARTICLE a b Deciduous coniferous tree Herbaceous perennial Evergreen coniferous tree Deciduous broadleaved shrub Deciduous broadleaved tree Heat requirement 0 10 20 30 40 50 c Larix Tilia Viburnum Syringa Salix Fagus Cornus Juglans Alnus Fraxinus Picea Betula Pinus Prunus Populus Quercus Acer Chilling accumulation 0 1 2 3 4 5 6 7 Number of species Fig. 1 Conceptual scheme of the relationship between HR and CA. Previous studies (Supplementary Table 1) found that saplings and twigs of 91 perennials need less heat for budburst (or leaf-out and flowering) if they receive more or longer chilling in natural or controlled environments. a Chilling accumulation (CA) is linearly or nonlinearly (the figure only shows a nonlinear relationship) negatively correlated with heat requirement (HR). The subplot shows the number of species exhibiting such a relationship for each life form (b) and genus (c). The CA–HR relationships for each species are shown in Supplementary Fig. 1. and Corylus avellana (widely grown hazelnut plants for com- Fig. 10). Across different locations, the simulated mean date mercial nut production) (Supplementary Fig. 7). correlated significantly with the observed mean dates with R2 In order to assess the scale effect caused by large grid size of E- ranging from 0.65 and 0.72 (p < 0.01) for all the models. However, OBS data set (about 10 km) and the non-uniform cover of plants when comparing the root-mean-square error (RMSE), the valid and temperature in the grid cells, we compared the HR and CA models (RMSE = 4.2–6.9 days) were overall more accurate than based on E-OBS data set and a station-based Global Historical the invalid model (RMSE = 6.5–10.7 days), because the invalid Climatology Network (GHCN) data set. The results showed that models produce an earlier leaf-out date than observation in cold CA or HR calculated by these two data sets was significantly regions but later leaf-out date in warm regions. Second, we correlated, and the CA–HR relationship was stable among two assessed whether different chilling models could reproduce the data sets (see the example of B. pendula in Supplementary Fig. 8). observed temporal trends of spring events in Europe. Across all Thus, the scale effect did not affect the relationship between CA species and locations, the different models reproduced similar and HR. We also evaluate the difference caused by the different trends with the observed data from 1980 to 2018, but the invalid starting dates of temperature accumulation. Compared to models simulated an earlier spring event compared to the temperature accumulated since 1 January, the later starting date observed data from the 1990s to the 2010s (Fig. 4a). Also, we (15 January and 1 February) would result in less HR (based on compare the simulated and observed trends from 1980 to 2018 at Model F1), especially for areas with a warm winter (see the different locations for the leaf out of B. pendula (Supplementary example of B. pendula in Supplementary Fig. 9). Although the Fig. 11). For all models, the simulated and observed data exhib- starting dates of heat accumulation altered the slope of HR ited consistently earlier trends at most locations. However, the against CA (based on Model C1), the significantly negative valid model (R2 = 0.15–0.17, RMSE = 0.18–0.20 days year−1) correlation between HR and CA remained. performed better than the invalid models (R2 = 0.02–0.09, RMSE = 0.22–0.34 days year−1) when comparing the R2 and RMSE between the simulated and observed trends. The invalid model Phenological trends based on chilling models. We tested the seemed to overestimate the past earlier trends in the leaf-out date impact of the choice of the chilling model on past and predicted (Supplementary Fig. 11). future trends in spring phenology by developing a process-based Based on the daily future climate data, the predicted advance in model incorporating the relationship between CA and HR. First, spring phenology during the 2020s to the 2090s was generally we assessed whether different chilling models could reproduce stronger for the invalid than the valid chilling models for all spatial gradients of spring events across warm and cold regions in species in scenarios for Representative Concentration Pathways Europe. Most of the species showed a similar result, so here we (RCPs) 4.5 and 8.5 (Supplementary Figs. 12 and 13). Averaged only show an example of leaf out for B. pendula (Supplementary from all species, the date of onset of spring events gradually NATURE COMMUNICATIONS | (2020)11:4945 | https://doi.org/10.1038/s41467-020-18743-8 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18743-8 a Model C1 b Model C2 c Model C3 d Model C4 65°N 60°N 55°N 50°N 45°N 40°N 35°N e Model C5 f Model C6 g Model C7 h Model C8 65°N 60°N 55°N 50°N 45°N 40°N 35°N i Model C9 j Model C10 k Model C11 l Model C12 65°N 60°N 55°N 50°N 45°N 40°N 35°N 0° 10°E 20°E 0° 10°E 20°E 0° 10°E 20°E 0° 10°E 20°E Linear trends in chilling accumulation (1951 to 2018) 0 600 km Decrease (p < 0.05) Decrease (p > 0.05) Increase (p > 0.05) Increase (p < 0.05) Fig. 2 Spatial patterns of the linear trends in chilling accumulation based on 12 chilling models. The chilling accumulation from 1 November to 30 April was calculated using the 12 chilling models. a–l Models C1–C12. Dark blue circles, significant and negative slopes (two-sided F-test, p < 0.01, n = 68 years). Light blue circles, insignificant and negative slopes (two-sided F-test, p > 0.05, n = 68 years). Dark red circles, significant and positive slopes (two-sided F- test, p < 0.05, n = 68 years). Light red circles, insignificant and positive slopes (two-sided F-test, p > 0.05, n = 68 years). advanced until the 2060s and then remains stable under RCP 4.5 Discussion but continuously advanced under RCP 8.5 (Fig. 4). The An earlier leaf-out date can increase the photosynthetic produc- phenological change matched the warming trends of two climate tion of forests22,23, so using an invalid chilling model to predict scenarios (Supplementary Fig. 14). The advance in spring events the start of the growing season would likely overestimate terres- under RCP 4.5 from the 2010s to the 2090s averaged for all trial photosynthesis and carbon uptake in spring. In addition to species was 11.2–13.7 days for the valid models but 16.8–22.5 days predicting future phenological change, the chilling models were for the invalid models (Fig. 4c). The advances in spring events used to estimate the chilling requirement of commercially under RCP 8.5 by the end of this century averaged 20.0 days important fruit and nut trees, such as apples, pears, and larger for the invalid than the valid models (Fig. 4d). cherries24,25. The cultivars selected based on the amount of winter 4 NATURE COMMUNICATIONS | (2020)11:4945 | https://doi.org/10.1038/s41467-020-18743-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18743-8 ARTICLE 800 a 800 b 800 c 800 d Model C1 Model C2 Model C3 Model C4 700 700 700 700 600 600 600 600 Heat requirement 500 500 500 500 400 400 400 400 300 300 300 300 200 200 200 200 100 100 100 100 0 0 0 0 0 50 100 150 200 250 0 50 100 150 200 250 0 20 40 60 80 100 120 140 0 50 100 150 200 250 800 e 800 f 800 g 800 h Model C5 Model C6 Model C7 Model C8 700 700 700 700 600 600 600 600 Heat requirement 500 500 500 500 400 400 400 400 300 300 300 300 200 200 200 200 100 100 100 100 0 0 0 0 0 50 100 150 200 250 0 50 100 150 0 50 100 150 0 50 100 150 800 i Model C9 800 j Model C10 800 k Model C11 800 l Model C12 700 700 700 700 600 600 600 600 Heat requirement 500 500 500 500 400 400 400 400 300 300 300 300 200 200 200 200 100 100 100 100 0 0 0 0 0 2 4 6 8 0 50 100 150 0 50 100 150 0 50 100 150 Chilling accumulation Chilling accumulation Chilling accumulation Chilling accumulation Fig. 3 Relationship between CA and HR for leaf-out of Betula pendula. Chilling accumulation (CA) from the previous 1 November to the leaf-out date was calculated using the 12 chilling models. Heat requirement (HR) was an integration of daily mean temperatures >0 °C from 1 January to the leaf-out date. a–l Models C1–C12. Data are divided into 20 groups with the same number of records. Error bars, mean values ± SD of HR or CA within each group (n = 6301 records). The linear regression lines (blue for negative correlations and red for positive correlations) are shown for each model. The shaded area shows the uncertainties, which were estimated by using the 95% confidence interval for slope and interception of regression lines. chilling at their production sites would likely fail to grow nor- effective) to measure the amount of chilling, winter warming mally in another location if an invalid model was used to estimate would lead to more number of days with the daily mean tem- the adaptative capacity of CA. In this study, the amount of perature higher than 5 °C, and thus decrease the CA. However, chilling based on invalid models even increase at mid- and high for Model C3 (temperature between 0 and 5 °C is effective), this latitudes under climatic warming (Fig. 2), possibly because in effect would be counteracted by the increase in the number of these areas the daily mean temperature during the coldest period days with temperature >0 °C. Considering the increased HR of of winter was below 0 °C, which did not contribute to CA for spring phenology in Europe12, Model C1 would produce a most of the invalid chilling models (e.g., Model C3 regarded that negative relationship between CA and HR, but Model C3 may only the temperature between 0 and 5 °C was effective). Due to produce an opposite result. Thus, our study suggested that the winter warming, the number of days with temperature below 0 °C valid models, which considered the effect of freezing temperatures reduced, causing an increase in the amount of chilling. If using on breaking dormancy, could better explain the CA–HR the invalid models to estimate the amount of chilling, the mechanism. Further experiments, however, are necessary to increased trends in the amount of chilling would continue at mid- confirm this assumption, such as giving saplings or twigs a and high latitudes in scenarios of future climatic warming26, so treatment of different chilling temperatures of the same duration growers in these regions may miss the opportunity to adapt to and then a test of regrowth under the same forcing conditions9. climate change, e.g., by breeding tree cultivars for lower chilling However, the effectiveness of temperatures <0 °C is rarely tested requirements. before, because plants cannot survive if they are exposed to The effect of freezing temperatures is the most important freezing temperatures immediately at the end of the growing structural difference between valid and invalid chilling models season27. Tolerance to freezing under natural conditions gradu- (Supplementary Fig. 15). The rate of chilling (the effectiveness of ally increases with exposure to low temperatures (a process different temperatures on dormancy break) for the invalid models known as cold acclimation), which allows plants to survive winter is zero or low for temperatures <0 °C. The rate of chilling for the conditions28. Sudden exposure to freezing temperatures could valid models, however, is still effective even for temperatures < extensively damage plant tissues due to the lack of cold accli- −5 °C. For example, if using Model C1 (temperature <5 °C is mation (by frost cavitation). More complex experiments are NATURE COMMUNICATIONS | (2020)11:4945 | https://doi.org/10.1038/s41467-020-18743-8 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18743-8 a b 140 RCP 4.5 (11-year moving average) 140 RCP 8.5 (11-year moving average) Mean date spring events (day of the year) 130 130 120 120 110 110 100 100 90 90 80 80 70 70 60 60 50 50 1980 1990 2000 2010 2020 2030 2040 2050 2060 2070 2080 2090 2100 1980 1990 2000 2010 2020 2030 2040 2050 2060 2070 2080 2090 2100 60 c RCP 4.5 (Advance by the 2090s compared to the 2010s) 60 d RCP 8.5 (Advance by the 2090s compared to the 2010s) 50 50 40 40 Difference (days) 30 30 20 20 10 10 0 0 1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12 Chilling model Chilling model Fig. 4 Changes in spring phenology averaged for all species and stations for 2019–2099. Spring phenological events were simulated using linear regression between chilling accumulations based on the 12 chilling models (from 1 November to the onset of spring events) and the heat requirement (accumulated temperature >0 °C from 1 January to the onset of spring events) for all records. The observed (thick gray line) and simulated trends across spring events of all species over 1980–2018 are shown. Red colors represent valid chilling models (C1, C2, C4, C5, and C12), and blue colors represent invalid models. a An 11-year moving average of spring phenological events under Representative Concentration Pathways (RCP) 4.5. b An 11-year moving average of spring phenological events under RCP 8.5. c Advance in spring events by the 2090s compared to the 2010s under RCP 4.5. d Advance in spring events by the 2090s compared to the 2010s under RCP 8.5. Error bars in (c) and (d), mean values ± SD (n = 30 species). therefore necessary to confirm the effect of freezing temperatures models are used to predict future phenological change, the on the rate of chilling, e.g., by using gradually decreasing tem- advance in spring phenology at the end of this century would peratures to simulate cold acclimation in situ. Furthermore, be twice as early as predicted using the valid models, and future acclimation in plants may alter the HR of plants. The acclimation spring terrestrial net carbon uptake would consequently be could be accounted for by using different slopes and interceptions overestimated. to describe the CA–HR relationship for individuals growing at different climatic conditions. However, when comparing the Methods CA–HR relationship between low latitude (lower than 50.65° N) Phenological and climatic data. We used data from the Pan European Phenology and high latitude (higher than 50.65° N), we did not found a Project (PEP725)29, an open-access database with long-term plant phenological significant difference in CA–HR relationship (Supplementary observations across 25 European countries (http://www.pep725.eu/). The regional/ Fig. 16). Thus, we did not consider the impact of acclimation in national network partners of PEP725 are following a consistent guideline for phenological observations30 and prepare the data for submission to the PEP725 this study. database curators29. We selected 30 species for which sufficient observational data Overall, the advance in budburst and leaf formation of per- were available: 21 deciduous broadleaved trees or shrubs, 6 herbaceous perennials, ennials in recent decades were associated with spring warming. 2 evergreen coniferous trees, and 1 deciduous coniferous tree (Supplementary Due to the negative relationship between CA and HR, the Table 3). Particularly, our data set included one fruit tree (Prunus avium) and one nut tree (Corylus avellana) since some of the chilling models are specifically reduced winter chilling slowed the spring phenological advance. developed for fruit and nut trees. A total of 2,493,644 individual records from We report that the majority of the current chilling models failed 15,533 phenological stations were used. The stations were mainly distributed in to represent the negative relationship between CA and HR, using moderate climates in Central Europe (Supplementary Fig. 3). Four spring events 2,493,644 phenological records for 30 species from 1951 to 2018, based on the BBCH code were investigated: BBCH 10, 11, 60, and 69, representing first leaves separated, first leaves unfolded, first flowers open, and end of flowering, potentially because these models were developed for limited respectively31. species in specific geographic regions and did not consider the We used the E-OBS v19.0eHOM data set32 with a spatial resolution of 0.1 × 0.1° effectiveness of freezing temperature. If these invalid chilling for 1950–2018 for calculating CA and HR of the in situ phenological records. This 6 NATURE COMMUNICATIONS | (2020)11:4945 | https://doi.org/10.1038/s41467-020-18743-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18743-8 ARTICLE data set is provided by the European Climate Assessment & Data set project and Hänninen44 using previous experimental results for Finnish birch seedlings includes homogenized series of daily mean, minimum, and maximum (Eq. (10)). Zhang et al.45 recently fitted observational data to the triangular model temperatures. We also use the daily maximum and minimum temperature data for 24 plant species and found that a mean optimal chilling temperature of 0.2 °C from the GHCN data set33 to assess the scale effect. The GHCN data set contains and an upper limit of the chilling temperature of 6.9 °C were most effective. Model station-based measurements from over 90,000 land-based stations worldwide, but C12, therefore, uses the triangular form with parameters of 0.2 and 6.9 °C (Eq. (11)). only parts of PEP725 stations match with the GHCN stations. 8 For future climatic data (2019–2099), we used daily minimum and maximum > <1 2:516 ; ð9Þ 0.5 × 0.5°) under two climatic scenarios (RCP 4.5 and RCP 8.5). These data have > :  2:10 ðTþ4:66Þ3:10 3:13 Tþ4:66 10:93 e 10:93 else been bias-corrected by applying the method used in the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP)34, which were available on the ISIMIP 8 server (https://esg.pik-potsdam.de/projects/isimip2b/). > <0 T ≤  3:4 or T ≥ 10:4 CU11 ¼ Tþ3:4 3:4 5þ3:4 : T10:4 Chilling models. We used 12 chilling models to measure the amount of chilling. 5 <0 T ≤  6:5 or T ≥ 6:9 commonly used upper-temperature threshold for chilling is 5 °C11,18. Some studies, CU12 ¼ Tþ6:5 6:5 6:90:2 : 6:9T however, have suggested that subfreezing temperatures were effective36,37, so we 0:20 °C for calculating CA11,38. Models C1–C6 were developed temperature in °C. based on various combinations of the upper and lower temperature limits. The rate of chilling was 1 for daily temperatures <5 °C for Model C1 (Eq. (1)), between −10 Forcing models. Forcing models were used to measure HR for the spring events of and 5 °C for Model C2 (Eq. (2)), between 0 and 5 °C for Model C3 (Eq. (3)), <7 °C plants. The GDD model is the most commonly used forcing model, which assumes for Model C4 (Eq. (4)), between −10 and 7 °C for Model C5 (Eq. (5)), and between that the rate of forcing is linearly correlated with temperature if the temperature is 0 and 7 °C for Model C6 (Eq. (6)). The equations for Models C1–C6 are as follows: above a particular threshold. We mainly used Model F1 (Eq. (12)), which adopts a  1 T ≤5 temperature threshold of 0 °C46–48, for examining the relationship between CA and CU1 ¼ ; ð1Þ HR: 0 T>5  FU1 ¼ maxðT; 0Þ; ð12Þ 1 10 ≤ T ≤ 5 CU2 ¼ ; ð2Þ 0 T>5 or T<  10 where FU1 is the rate of forcing for Model F1, and T is the daily mean tempera- ture (°C).  1 0≤T ≤5 We also validated the chilling models by correlating them with HR based on CU3 ¼ ; ð3Þ seven other forcing models to assess the impact of the choice of forcing model on 0 T>5 or T<0 the results. Model F2 (Eq. (13)) is also a GDD model but has a temperature  threshold of 5 °C12,18: 1 T ≤7 CU4 ¼ ; ð4Þ 0 T>7 FU2 ¼ maxðT  5; 0Þ; ð13Þ  1 10 ≤ T ≤ 7 where FU2 is the rate of forcing for Model F2, and T is the daily mean CU5 ¼ ; ð5Þ temperature (°C). 0 T>7 or T<  10 Piao et al.48 found that leaf onset in the Northern Hemisphere was triggered  more by daytime than nighttime temperature. They thus proposed a GDD model 1 0≤T ≤7 CU6 ¼ ; ð6Þ using maximum instead of mean temperature. Models F3 (Eq. (14)) and F4 (Eq. 0 T>7 or T<0 (15)) are based on maximum temperature with thresholds of 0 and 5 °C, where CUi is the rate of chilling for Model Ci, and T is the daily mean tempera- respectively. ture (°C). FU3 ¼ maxðT max ; 0Þ; ð14Þ Model C7 is also known as the Utah Model39, which assigned different weights to different ranges of temperatures and was first used to measure the chilling requirements of peach (Eq. (7)). The Utah Model was modified to produce Model FU4 ¼ maxðT max  5; 0Þ; ð15Þ C8 (Eq. (8)), which removed the negative contributions of warm temperatures to accumulated chilling20. where FUi is the rate of forcing for Model Fi, and Tmax is the daily maximum 8 temperature (°C). > > 0 T ≤ 1:4 A recent experiment demonstrated that the impact of daytime temperature on > > > > 0:5 1:4 > > > impact of nighttime temperature49. Model F5 (Eq. 16) thus uses two parameters > < 1 2:4 > > >0 12:4 > > > > > 0:5 15:9 : 1 T>18 maximum and minimum temperatures (°C), respectively. Many studies have suggested that the rate of forcing followed a logistic function 8 > 0 T ≤ 1:4 of temperature44,50. Model F6 (Eq. (17)) uses a logistic function proposed by > > > > Hänninen44, and Model F7 (Eq. (18)) uses another logistic function proposed by > < 0:5 1:4 > 28:4 > > 0:5 9:10 > > FU6 ¼ 1þe ; ð17Þ : 0 else 0 T>12:4 where CUi is the rate of chilling for Model Ci, and T is the daily mean 1 temperature (°C). FU7 ¼ ; ð18Þ 1 þ e0:47Tþ6:49 Model C9 is a dynamic model developed for peach in Israel and South Africa40,41 and now adopted for apricot cultivars42. The most important where FUi is the rate of forcing for Model Fi, and T is the daily mean characteristic of Model C9 was that a previous intermediate product affected the temperature (°C). rate of chilling in the current hour or day. We did not provide equations for Model Model F8 is a growing degree hour (GDH) model, where species have an C9 for simplicity (see the equation in Luedeling et al.20). optimum temperature for growth and where temperatures above or below that Harrington et al.43 summarized published results for chilling units and optimum have a smaller impact51. Model F8 (Eq. (19)) was first designed for constructed a chilling function based on a three-parameter Weibull distribution, calculating HR at hourly intervals, but we applied it at a daily interval. The stress coded as Model C10 (Eq. (9)). Model C11 has a triangular form, which was fitted by factor in the original GDH model was ignored, because we assumed that the plants NATURE COMMUNICATIONS | (2020)11:4945 | https://doi.org/10.1038/s41467-020-18743-8 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-18743-8 were not under other stresses. land-based-datasets/global-historical-climatology-network-ghcn. Future climate data 8 (2019–2099) were available on https://esg.pik-potsdam.de/projects/isimip2b/. > > 0 TTc < Tu TL  >   FU8 ¼ 1 þ cos πþπ TTT L TL ≥ T ≥ Tu ; ð19Þ Code availability 2 u TL > >    > : ðT  T Þ 1 þ cos π þ π TTu All Matlab scripts, from the initial processing of data sets to final analyses, are archived u L 2 2 T T Tu