ARTICLE https://doi.org/10.1038/s41467-019-13365-1 OPEN Spatial variance of spring phenology in temperate deciduous forests is constrained by background climatic conditions Marc Peaucelle 1,2,6*, Ivan A. Janssens 3, Benjamin D. Stocker 1,2,7, Adrià Descals Ferrando 1,2, Yongshuo H. Fu4, Roberto Molowny-Horas1, Philippe Ciais 5 & Josep Peñuelas 1,2 1234567890():,; Leaf unfolding in temperate forests is driven by spring temperature, but little is known about the spatial variance of that temperature dependency. Here we use in situ leaf unfolding observations for eight deciduous tree species to show that the two factors that control chilling (number of cold days) and heat requirement (growing degree days at leaf unfolding, GDDreq) only explain 30% of the spatial variance of leaf unfolding. Radiation and aridity differences among sites together explain 10% of the spatial variance of leaf unfolding date, and 40% of the variation in GDDreq. Radiation intensity is positively correlated with GDDreq and aridity is negatively correlated with GDDreq spatial variance. These results suggest that leaf unfolding of temperate deciduous trees is adapted to local mean climate, including water and light availability, through altered sensitivity to spring temperature. Such adaptation of heat requirement to background climate would imply that models using constant temperature response are inherently inaccurate at local scale. 1 CREAF, Cerdanyola del Vallès, Barcelona 08193 CA, Spain. 2 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona 08193 CA, Spain. 3 Department of Biology, University of Antwerp, Wilrijk 2610, Belgium. 4 College of Water Sciences, Beijing Normal University Beijing, Beijing, China. 5 Laboratoire des Sciences du Climat et de l’Environnement, UMR 1572 CEA-CNRS UVSQ, 91191 Gif sur Yvette, France. 6Present address: Computational and Applied Vegetation Ecology - CAVElab, Department of Environment, Faculty of Bioscience Engineering, Ghent University, Gent, Belgium. 7Present address: Institute of Agricultural Sciences, Department for Environmental Systems Science, ETH Zürich, Universitätsstrasse 2, 8006 Zürich, Switzerland. *email: mpeau.pro@gmail.com NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 T he temporal variation in spring leaf unfolding (LU), spatial heterogeneity of spring LU and its GDDreq, reflecting including the effect of climatic warming, has been exten- evolutionary mechanisms through which plants have adjusted sively documented1–6. In temperate and boreal regions, their growth strategies in order to maximize their fitness under temporal variation of LU is dominated by the occurrence of warm those specific biogeographical conditions. We argue that bio- spring temperatures2, typically quantified by growing degree geographic constraints on plant phenology can be detected by days7–9 (GDD), with LU occurring when the supplied GDD analyzing the spatial response of stands long-term mean LU and equals the required GDD. The GDD requirement of LU (hereafter GDDreq, instead of their inter-annual variability. In this respect, GDDreq) is commonly defined as the accumulated air tempera- we hypothesize that not only spring, but also mean growing ture above a temperature threshold over the preseason. GDDreq season conditions are important controls of the spatial differences varies from year to year in response to differences in the occur- in leaf phenology among different locations. If so, these evolu- rence of low winter temperatures9,10 (chilling, hereafter NCD), tionary mechanisms control the sensitivity of LU to short-term and constraints imposed by optimal daylength11. Also, the spring temperature variations, and consequently are key com- availability of water and light have been shown to affect LU and ponents of the observed spatial variability in LU and its GDDreq. its GDD requirement11–13. Here, we relate spatial variations in LU and GDDreq in 8 While the temporal variation in LU and its GDDreq is exten- dominant European deciduous tree species (see Methods) to sively studied, the spatial heterogeneities of LU, and especially of long-term cross-sites differences in temperature, radiation and its controls, have been much less studied. Land surface models aridity (Supplementary Table 1). We combine long-term LU assume that the drivers of the temporal variation of LU are also observations at 27790 sites over 1970–2016 (Supplementary determining the spatial gradients in LU, but to our knowledge, Fig. 1) with climatic data from CRU-NCEP30 at a spatial reso- this assumption has not been thoroughly tested. The temperature lution of 0.5° (Supplementary Figs 2 and 3) to estimate the effects response of LU is often considered constant (spatially and tem- of long-term mean temperature, light (radiation and daylength) porally uniform), albeit species-dependent, in phenology mod- and aridity as determinants of the spatial variance of LU and els14. However, when applied at the regional scale, LU models GDDreq using statistical models. Our study provides evidence for were not able to accurately reproduce the observed spatial var- a significant control of leaf unfolding by long-term background iation of LU14–17. climatic conditions across sites, potentially representing long- Recent warming trends have advanced LU for virtually all term adaptation of species. For all species: (1) the spatial variance temperate-zone tree species, but there is large spatial hetero- of LU and GDDreq supports previous correlations with tem- geneity in these advancing LU trends at the global scale18,19. At perature and chilling requirement, but is also determined by the regional scale, different warming trends among sites have radiation intensity (W m−2) and site aridity, (2) the spatial var- been proposed as the main cause for these spatially differing iance of GDDreq is better explained by preseason radiation trends in mean LU date between high and mid-latitudes20,21, low intensity than day length at LU, (3) LU and GDDreq are more and high elevations22–24 or coastal and inland areas25,26. Pre- sensitive to water availability on drought-prone sites and (4) LU cipitation has also been identified as a key controlling factor of and GDDreq respond both to long-term preseason and growing- spatial differences of vegetation green-up in arid and semi-arid season climatic conditions. These findings suggest that at least regions12,27, while at high latitudes, photoperiodic controls were two mechanisms influence spring phenology: (i) the direct sen- proposed as an evolutionary safety mechanism that mitigates the sing of meteorological conditions during spring to optimize the risk of frost damage28,29 and causes LU to occur later at higher restart of plant activity and (ii) the long-term adjustment of bud latitudes. sensitivity to spring meteorological conditions in order to cope Studies assessing the effects of multiple environmental variables with growing season pressures at sites. It also suggests that have highlighted complex and species-specific behaviors12,27, common GDDreq and NCD metrics used for simulating the mainly due to the mixing of temporal and spatial aspects of temporal variability of LU are not suitable for spatial studies and phenology, rather than trying to deconvolute these two compo- should be used with caution. nents. Indeed, the spatial and temporal variances of phenology are expected to reflect two different mechanisms of control. First, short-term, fast responses to changing weather should drive Results & Discussion temporal variations in LU and its GDDreq, aiming to maximize Importance of background climate for LU spatial variability. fitness under inter-annually varying weather conditions. Second, Among-site spatial variability of long-term mean annual LU (day of an adaptive response to local biogeographical conditions may year -DoY-) and GDDreq (°C; see Methods) was of the same order maximize tree fitness under the local long-term mean climatic of magnitude than the typical interannual variation in LU and conditions and would select for spatially optimized LU and cli- GDDreq (Table 1). Figure 1 and 2 present the importance of each mate sensitivity, inducing spatial variation therein. Biogeographic climatic variable (taking into account co-linearity and spatial auto- constraints on LU include all environmental variables that impose correlation; see Methods) in determining median LU date and long-term adaptation of LU and its GDDreq to optimize fitness GDDreq. The full model, including long-term temperature, incident under local conditions. These include climatic variables, such as surface radiation and water availability predictors, explained 61 ± site-specific occurrence of late frost events, drought occurrence, 7% of the spatial variance of LU among species (regression coeffi- low or high light extremes, that may need to be avoided and cients in Supplementary Table 2). Chilling and GDDreq, which are therefore require shifts in growing season to enable maximum tree functions of temperature and are additionally modulated by pos- fitness. Also, site-specific interactions with neighboring competi- sible biotic adaptation to site-specific conditions, accounted for only tors, pathogens and herbivores may induce spatial differences in half of the spatial variance of LU (Fig. 2a). Average preseason LU and its weather dependency, in order to maximize tree fitness. temperature (TP) was positively correlated with site LU and was Taken together, this suggests a complex response of plant phe- selected as the best predictor instead of GDDreq for half of the nology to climate change, but also that models of LU that apply studied species. In addition to preseason temperature conditions, spatially uniform parameters may not capture the observed pat- growing season temperature (TG) was found to explain 29 ± 5% of terns of LU and its GDD requirement. the spatial variability in LU. Considering all temperature-related The key hypothesis that we explore in this study is that long- variables as predictors, including preseason and growing season, term mean background biogeographical conditions determine the models captured 52 ± 6% of total LU’s spatial variance. 2 NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 ARTICLE Table 1 Observed LU and GDDreq variability. Species Intra-site LU sd Inter-Site LU sd Intra-site GDDreq sd Inter-site GDDreq sd # sites # drought prone sites A. hippocastanum 8.9 8.4 91.0 80.0 4429 1052 A. glutinosa 9.7 8.7 105.4 90.2 3207 687 B. pendula 8.95 7.5 85.5 78.0 4577 1079 F. sylvatica 7.3 6.5 105.0 100.0 3974 785 F. excelsior 8.2 7.5 126.8 106.9 3003 518 Q. robur 7.7 7.9 110.9 96.5 3995 892 S. aucuparia 7.7 7.5 93.4 132.8 2327 472 T. cordata 8.1 7.6 97.2 81.4 2278 262 Description of the intra and inter-site variability in observed leaf unfolding (LU) dates from the PEP dataset and the corresponding estimated Growing Degree Days requirement (GDDreq), as well as the number of sites (#) per species. Source data are provided as a Source Data file Long-term mean LU across sites was also significantly spring phenology in the long-term. Based on the observed rela- positively correlated with long-term growing season mean tionships in Fig. 3, we therefore selected drought-prone sites (αE incident shortwave radiation (SWG, visible and near infrared) < 0.9) to assess if their responses to temperature, radiation and for all species (Fig. 2a), while no clear effect of long wave water differed (results for wet sites can be found in Supplemen- radiation (LWG, infrared) was observed. SWG contributed an tary Fig. 5 and Supplementary Tables 6 and 7). additional 9 ± 2% of explained variance in LU. The importance of At drought-prone sites, TP outcompeted GDDreq as the main SWG in the model determining LU was low compared to that of determinant of the spatial variance of LU. The contribution of temperature (Fig. 1a), but still statistically significant (Supple- temperature (NCD, GDDreq, TP and TG) did not change at mentary Table 2). We also examined the relations between LU drought-prone sites relative to all sites (Fig. 2b). However, the and precipitation (P), soil-moisture content (SM) and the ratio of relative importance of preseason precipitation (PP) was doubled actual to potential evapotranspiration31 (αE; see Methods). Long- at those drought-prone sites compared to the wet-sites analysis term mean LU was significantly and negatively correlated to long- (Supplementary Tables 4 and 6). Across the drought-prone sites, term mean growing-season P (PG) for A. hippocastanum, B. LU was negatively correlated with PP for A. hippocastanum, A. pendula, F. sylvatica, F. excelsior and Q. robur, for which it glutinosa, B.pendula, F. sylvatica and Q. robur. Instead, for A. captured 4.0 ± 0.6% of LU’s variance. Only T. cordata exhibited a glutinosa and S. aucuparia, LU was significantly and positively negative correlation with preseason P (PP) and SM (SMP). correlated with PG. LU also was positively correlated with αE for Growing season SM (SMG) was not correlated to LU, while αE F. excelsior, suggesting a differential effect of summer water showed a weak, but significant, positive correlation for F. sylvatica availability on LU for these three species when grown at drought- and F. excelsior, consistent with the relationship with PG (Fig. 2a; prone sites. We also observed a smaller and divergent effect of Supplementary Table 2). incoming radiation on LU at drought-prone sites (Fig. 2b) compared to the wet-sites analysis (Supplementary Fig. 5). Importance of background climate for GDDreq spatial varia- As for LU, the relative importance of temperature and bility. Long-term temperature, incident surface radiation and incoming radiation as determinants of the spatial variance of water availability explained 54 ± 12% of the spatial variance of GDDreq across drought-prone sites did not differ from the other GDDreq among species (Figs 1b and 2b; regression coefficients sites (Fig. 2d). However, we observed a significant contribution of can be found in Supplementary Table 3). Large differences in the water availability at drought-prone sites compared to other sites. impact of long-term conditions on GDDreq are observed among GDDreq was negatively correlated with PG and positively species, with background climate explaining between 33% of the correlated with PP, for all species except for T. cordata. GDDreq spatial variance in GDDreq for B. pendula and 73% for S. was also negatively correlated with αE for A.hippocastanum and aucuparia. A.glutinosa, while being positively correlated for T. cordata, For all species, GDDreq was negatively correlated with NCD which was consistent with the observed correlation with PG and positively correlated with TG, together explaining 52 ± 4% of (Fig. 2e, Supplementary Tables 5 and 7). the spatial variance in GDDreq (Fig. 2c). The proportion of spatial While we observed significant differences in explaining LU and variance in GDDreq explained by preseason incoming shortwave GDDreq spatial variance at drought-prone sites, the same radiation (SWP) was about 30%. Day length was selected as a approach showed no differences between warm/cold sites or predictor of GDDreq only for A. hippocastanum and S. aucuparia. low/high light sites. Growing season light conditions played a minor role, albeit still statistically significant, with incoming longwave radiation (LWG) GDD and NCD did not fully capture LU spatial variance. We capturing about 4.0 ± 1.5% of GDDreq’s explained variance focused on the spatial variance of spring LU and its GDDreq (Fig. 2c). GDDreq was significantly correlated with both PG and requirement. If the determinants of the temporal variation of LU PP, capturing together around 7.5 ± 2.3% of GDDreq variance. and GDDreq would not fully explain their spatial variation, this However, no significant correlations were found with soil would imply adaptation of spring phenology to long-term local moisture, nor αE. background environmental conditions. Results revealed that temperature and the interplay between chilling and heating Different response of phenology in drought-prone sites. While during winter and spring, the main determinants of the temporal precipitation was a weaker determinant of the spatial variance in variation of LU, also were important factors controlling the GDDreq than temperature and light, for all species we observed spatial variation of spring LU, in line with previous studies9,10. that the spatial variance in GDDreq decreased with increasing However, we also showed that these chilling and heat require- drought stress (αE) during the growing season (Fig. 3; Supple- ment metrics only captured 30% of the total spatial variance in mentary Fig. 4), suggesting that frequent and/or more intense LU for all species (Fig. 2a, Supplementary Table 2). The positive water stress during the growing-season results in adaptation of correlation observed between LU and TP was not expected. NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 a R 2 = 0.67 b R 2 = 0.48 800 1000 1200 150 50 40 140 Predicted GDDreq (°C) Predicted LU (DoY) 120 130 600 100 110 400 40 200 30 90 0 80 90 100 110 120 130 140 150 0 500 1000 1500 Fraction of explained R 2 Fraction of explained R 2 Observed LU (DoY) Observed GDDreq (°C) 30 20 20 10 10 0 0 L G G q TG PG PP G P D TP D TG G G P P PG PP G P αE αE re D SM D SW SM C LW SW LW SM C SW LW SM D N N G c d 1000 1200 R 2 = 0.68 R 2 = 0.56 50 40 140 Predicted GDDreq (°C) Predicted LU (DoY) 800 120 600 100 400 40 200 30 80 0 90 100 110 120 130 140 150 0 500 1000 1500 Fraction of explained R 2 Observed LU (DoY) Fraction of explained R 2 Observed GDDreq (°C) 30 20 20 10 10 0 0 G G q TG PG PP G P D TG G G PG PP G P P D P TP αE αE re SW SM SM LW D SW C SM C LW SW LW SM D N N G Fig. 1 Effect of background climate on LU and GDDreq. Relative importance (%) of each variable in estimating (a) observed LU date for all sites and species (n = 27790), (b) GDDreq for all sites and species, (c) LU date for drought-prone sites (n = 5747) and (d) GDDreq for drought-prone sites. See the last row of Fig. 2a–d (ALL) for the direction of each effect. NCD, number of chilling days estimated as the number of days between 1 November in the previous year and the LU date with temperatures between 0 and 5 °C; TG, mean growing-season temperature; TP, mean pre-season temperature; SWG, mean growing- season shortwave [visible and near infrared] radiation; LWG, mean growing-season longwave [infrared] radiation; SWP, mean pre-season shortwave [visible and near infrared] radiation; LWP, mean pre-season longwave [infrared] radiation; PG growing-season total precipitation; PP, pre-season total precipitation; DL, day length at LU date; SMG, growing season soil-moisture content; SMP, pre-season soil-moisture content αE, ratio of actual to potential evapotranspiration. drought-prone sites were defined as sites with long-term αE <0.9. Variables were selected based on penalized elastic net regression and corrected for spatial autocorrelation; the coefficient of determination (R²) of the selected models is given at the top of each panel. Error bars represent the 95% confidence interval. The sum of the relative importance equals 100% of R². The scatter plot in each panel represents modeled versus observed variables (LU or GDDreq) after correcting coefficients for collinearity and spatial autocorrelation. The red line represents the 1:1 line. Source data are provided as a Source Data file. 4 NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 ARTICLE a LU c GDDreq q re G G P G D G G D P G P D P SW SW D SW SM SM C PG LW LW SM SM TG αE C PG PP LW TG L αE PP TP G N D N A. hippocastanum 27.1 26.9 27.7 12.4 1.1 3.9 0.7 0.2 A. hippocastanum 26.6 25.6 7 1.7 3.9 27.1 4.5 3.7 A. glutinosa 34.9 31.1 24.4 7.6 1.6 0.3 0.2 A. glutinosa 27.3 24.2 3.2 3.1.7 5.3 5.6 2.6 B. pendula 26.7 27.3 29.9 9.1 1.6 4.3 0.2 0.5 0.4 B. pendula 25.7 31.2 4.4 29.2 1 8 0.2 0.2 F. sylvatica 29.8 28 26.2 7.8 0.4 4.5 0.9 0.2 2.2 F. sylvatica 29 30.5 4 22.3 3.5 5.2 4.9 0.5 ALL F. excelsior 29 24.2 35.4 5.7 3 0.1 2.7 F. excelsior 30.5 20.2 5.9 24.9 8 5.6 3.5 1.4 Q. robur 25.5 29.4 30.2 7.2 1.4 4.1 2 0.1 Q. robur 23.9 25.1 9.5 4.8 24.8 4.6 2 5 0.2 S. aucuparia 38 20 27.5 10.3 3.4 0.3 0.1 0.3 S. aucuparia 25.1 29.5 5 5.6 7.6 23 3.5 0.2 0.4 T. cordata 2.8 42.3 34.4 11.6 3.4 3.2 2.3 T. cordata 7.9 39.3 4.9 22.3 17.1 4 4.3 0.4 ALL 19.8 29.6 33.6 9.5 2 3.8 1.1 0.5 ALL 23 26.8 2 37.1 1.1 5.3 4.1 0.2 0.2 Summary Summary b d q re G P G G D P G P D G G D P SW SW SM SM C PG LW LW TG αE D SW PP L SM SM C PG LW TG αE PP TP N D G N A. hippocastanum 30.1 27.3 29.9 0.8 10.6 0.3 0.8 A. hippocastanum 37.1 8.7 22.6 15.6 3.3 1.3 2.6 8.8 A. glutinosa 25.8 26.6 11.1 6.3 9.3 9.4 5.9 5.6 A. glutinosa 28.3 9.7 29.8 9.1 11.6 4.2 7.3 B. pendula 24.3 33.8 27.8 0.3 1.4 9.4 2.7 0.2 B. pendula 24.3 30 1.9 21.8 7.2 2.8 7.8 3.2 1 DROUGHT-PRONE F. sylvatica 27.7 22.2 29.3 5 2.9 11.4 1.5 F. sylvatica 27.2 26.8 9.6 16.2 8.8 0.1 8.6 2.3 0.6 F. excelsior 31.4 40.7 13.8 2.4 2.9 8.8 F. excelsior 22.3 19.3 20.2 19.1 5.3 11.1 2.1 0.4 Q. robur 23.7 26.7 30.7 11 1 6.8 0.1 Q. robur 23.4 29.8 19.1 15.2 5.1 7.2 0.1 S. aucuparia 46.2 27.1 10.3 10.6 3.3 2.4 S. aucuparia 11.6 41.5 30.5 11.6 1 2 1.7 T. cordata 19.8 38.2 40.2 1.9 T. cordata 16.9 24.8 9.9 30.2 11.3 2.4 4.5 ALL 22.7 26.5 36.6 4.4 0.2 9.2 0.3 0.1 ALL 22.3 27.2 32.8 9.1 0.2 6.5 1.2 0.7 Summary Summary Fig. 2 Effect and direction of background climate on LU and GDDreq for each species. Relative importance of each variable (in percentage) in explaining leaf unfolding date (LU) and corresponding heat requirement (GDDreq, estimated as the sum of temperatures > 5 °C between 1 January and the LU date) for each species considering (a–c) all sites, (b–d) only drought prone sites. The sum of each row is equal to 100% of model R². Red colors indicate a positive correlation, while blue colors indicate a negative correlation. Color intensity reflects variable relative importance. Blanks represent variables that were discarded during the predictors selection. See caption of Fig. 1 for the description of each variable and Table 1 for samples size. The direction and importance of correlations among species is summarized by the size and color of arrows at the bottom of each panel. A double black arrow means that the direction of the response is species dependent. See Supplementary Tables 5–10 in appendix for a complete description of correlation coefficients. 200 Higher TP implies that GDDreq is reached earlier and should thus Species correlate negatively with LU. This positive correlation thus sug- A. hippocastanum gests an adaptation of LU to background temperature. 150 A. glutinosa Even when taking NCD, GDDreq and TP into account, which GDDreq iqr (°C) B. pendula all correlate with pre-season temperature, TG was still selected as F. sylvatica the major explanatory variable of LU’s spatial variance. Two 100 F. excelsior hypotheses could explain this result: (1) TG is a good integrative Q. robur 50 proxy of site biogeographical constraints on LU, and potentially S. aucuparia summarizes other variables not considered in our study, or (2) T. cordata trees growing at different locations have optimized the control of 0 their LU in response to biogeographical differences in preseason 0.7 0.8 0.9 1.9 and growing season conditions32. In both cases, this result αE indicates that adaptation to long-term mean site biogeographical Fig. 3 Effect of aridity on GDDreq spatial variance. Variation of the conditions, including growing season conditions but also a suite interquartile (iqr) of the long-term (1970–2016) heat requirement (GDDreq, of biological interactions that could not be included in this study, estimated as the sum of temperatures > 5 °C between 1 January and the LU constitutes an important evolutionary mechanism to optimize LU date) as a function of the ratio of actual to potential evapotranspiration at that location, and must be considered in addition to commonly (αE). Low αE values indicate drought conditions during the growing season. used preseason temperatures, which do not suffice to explain the The black line represents the iqr of GDDreq when considering all sites and spatial distribution of LU. species. See Supplementary Fig. 4 for the distribution of GDDreq values along αE. Source data are provided as a Source Data file. New models are needed for regional studies. In our study, long- term GDDreq was strongly correlated to long-term preseason NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 temperature (r = 0.95, Supplementary Fig. 6), indicating that water conditions on LU and its temperature sensitivity. Leaf GDD models are not independent of sites TP. With constant unfolding was also correlated with light and water conditions. GDDreq, higher temperature would imply earlier LU. LU indeed The common choice of using GDD models (even when making occurs earlier in warm areas and later in cold areas, but much less GDDreq dependent on NCD and day length) ignores the effects of than expected based on the temporal GDDreq. Trees have the availability of water and light. This partly explains why GDD adjusted their GDDreq to avoid too late LU in cold areas and risk models typically exhibit large uncertainties when used at regional damage with too early LU at warmer sites. In principle, GDDreq scales14,35. Despite their capacity to well explain the inter-annual should accounts for all warming effects, and therefore be constant variability of LU, the strong correlation between GDDreq and unless another driver is at play or adaptation has occurred. So long-term mean background climate observed in our study even if LU is partially taken as input to calculate GDDreq, they suggests that GDD models, in their current rigid parameteriza- would be uncorrelated. The fact that they are correlated highlights tion, are not suitable for studying phenology at the regional scale that acclimation/adaptation has occurred. The correlation if they do not include biogeographical constraints, especially for between GDDreq and TP and between LU and GDDreq are in line regions where precipitation or light are key controls of LU. with previous remote sensing studies33,34 and strongly suggest adaptation of spring phenology to long-term site temperature. It also highlights that, even if generally useful to describe the inter- Radiation intensity matters more than day length for GDDreq. annual variability, GDDreq in its current definition is a poor proxy Recent experimental spring phenology data have indicated that of biogeographical constraints of leaf unfolding. only 35% of northern hemisphere woody species relied on day To emphasize this point, we additionally looked at the length (DL) as a control of LU, and that the dependence on DL commonly used negative relationship between GDDreq and occurred predominantly at mid latitudes of the northern hemi- NCD. We observed that the relationship between NCD and sphere36. Other studies have described a greater impact of day- GDDreq for LU is modulated by SWP (Fig. 4). We argue that this time than nighttime temperatures on LU date37,38 as well as relationship is mainly an artifact induced by the fact that we differences in GDDreq and NCD requirement depending on DL11. applied uniform GDDreq and NCD definition for all sites. Because Day length and LU date are highly collinear (Supplementary trees can have different sensitivity to pre-season temperature in Fig. 6), making it difficult to separate the impacts of day length different regions (related to light and water availability), it implies (or preseason radiation) on LU. However, several previous studies that different GDDreq and NCD definitions have to be used at the did suggest a modulation of the temperature sensitivity of LU by spatial level for different sites if we want to be able to simulate LU light11,36–38. In our study, growing and preseason radiation at the regional scale. Compared to current regional models using intensity received by plants were clearly identified as an important constant GDDreq definition independently of the studied region, explanatory variable of GDDreq’s spatial variance (Figs 1 and 2) our results (Figs 2 and 4) suggested that northern sites need to and were more important than day length in explaining the have lower temperature threshold for temperature sum and/or spatial differences in GDDreq. lower GDDreq thresholds than southern sites. While the effect of light intensity was clear, contributions to the Since leaf unfolding determines the restart of the growing spatial heterogeneity of LU and GDDreq were very different season, we expect a requirement/threshold effect of light and between SWP and LWP. The ratio of visible to infrared light 1000 SWP (W m–2) GDDreq (°C) 200 175 150 125 100 500 75 0 0 25 50 75 100 NCD (days) Fig. 4 Adaptation of chilling-forcing requirement to preseason radiation. Scatterplot of average site long-term (1970–2016) GDDreq (°C) versus NCD (days) for all species and all sites. The color gradient represents the average site long-term incoming shortwave preseason radiation (SWP, W m−2). Gray color represents the median value of SWP (130 W m−2) of all sites, blue and red colors correspond to values lower and higher than the median SWP, respectively. Source data are provided as a Source Data file. 6 NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 ARTICLE determines several plant processes, such as the state of species responded differently to water availability at drought-prone phytochrome photoequilibrium, which controls growth rate, sites. A. glutinosa, T.cordata and F.excelsior showed opposite foliar and chloroplast development and even apical dominance39. correlation with water availability compared to the other species Some evidence also suggests that light modulates internal (Fig. 2). T. cordata and F. excelsior were proposed as species with a hormone-regulated growth40 and protein production in plants41 relatively high drought tolerance compared to other European by affecting signaling pathways of ethylene and abscisic acid, two species50,51, while A. glutinosa naturally occurs in wet sites. This phytohormones involved in bud set and leaf development42. A different behavior observed for A. glutinosa, T.cordata and F. recent experiment found a phytochrome-mediated photoperiodic excelsior could reflect either their drought tolerance, or the fact that control for Fagus sylvatica43, while a recent review highlighted the differences in edaphic heterogeneity were not captured by aggre- effects of the light spectrum on spring and autumn phenology44. gating data at the pixel level, resulting in a mismatch between We therefore hypothesize that trees adjusted their life cycle to the actual water availability and the estimates of αE used here. average light spectrum of the site that can be directly sensed by Differences were also observed for some species regarding their buds and plays a direct role in enhancing or inhibiting LU, but heat requirement. Biogeographical conditions of temperature, this remains to be experimentally verified also for other species. light and water only captured 33, 44 and 49% of GDDreq spatial Light plays a key role in plant activity during the growing variance for B.pendula, F. sylvatica and A. glutinosa, respectively. season by regulating photosynthesis and growth. The growing Previous studies showed a reduced sensitivity of leaf unfolding to season is also the period during which buds are created. In the climate warming in the last decades, mainly attributed to plant long-term, we expect that trees respond to growing season plasticity52. However, how trees acclimate or adapt to future meteorological conditions during which the formation of buds, climate change remains unclear and might be species-dependent, but also carbon reserves, can be affected, which in turn affects the due to differences in temperature, light and water sensitivities. sensitivity of buds to temperature in the subsequent winter and The remaining unexplained variance of LU and GDDreq can be spring as already shown with Populus tremula45. However, attributed to uncertainties in observations and climate data, with differences in stand canopy openness, leaf area or even plant a potential effect of data aggregation at the pixel level, but also of activity can lead to large uncertainties in the potential effect of unaccounted biotic and edaphic factors (e.g., stand age, growing season radiation on tree eco-physiology due to biogeography of pathogens or mycorrhizal associations, soil differences in local light regimes. As for TG, both SWG and structure and fertility, etc.). LWG may be good proxies of background biogeographical In the end, the co-limitation of spring phenology by light and constraints without having a direct effect on phenology. aridity may account for why LU does not keep pace with climate We argue that intensity of growing and preseason incoming change19,52, which may have vast and far implications on the radiation should be included in phenological studies, not only day carbon cycle53, with a possible alteration of the competitive length as a proxy of photoperiod. More research, especially balance among species54. experimental studies, is needed to clearly distinguish among the effects of the light spectrum, light intensity and day length on LU and its required GDDreq. Tree seasonality affects spring phenology. Temperature, avail- ability of light and water and their interactions thus affect the spatial heterogeneity in LU date and its associated GDDreq. On The role of aridity: The water-saving hypothesis. No theory has one hand, our results showed that LU date adjusted to spring yet been accepted that accounts for the effect of water on spring meteorological conditions, for which trees are co-limited by TP, foliar phenology in temperate forest ecosystems. Here, we propose SWP, LWP and PP. On the other hand, we observed different a hypothesis on the effect of drought stress on LU dates and responses of spring phenology to preseason and growing-season GDDreq. LU and GDDreq were more strongly correlated to water meteorological conditions, especially visible for drought-prone availability at drought-prone sites (Fig. 2) than at wetter sites, sites, highlighting that tree adjust their phenology to cope with potentially reflecting a long-term adaptation of trees to frequent seasonality on the long-term. drought stress. Temperature sensitivity tends to be lower in water- Since buds are formed during the growing season, we argue limited conditions, as indicated by the higher GDDreq with that any effect of growing season pressure might in turn have a increasing drought stress (Fig 2 and Supplementary Fig. 4). feedback on tree seasonality, with a potential long-term Organogenesis and primary growth in buds have been correlated adjustment of ecophysiological responses to these constraints. with hydraulic architecture46, and previous studies have high- In the long term, it may also alter the restart of plant activity in lighted a clear effect of growing-season water stress on bud pro- spring and thereby optimize the long-term growth, reproduction, duction and foliar development47,48. Here, we highlight an or survival of the trees, which influence the restart date of additional effect on the sensitivity of buds to temperature in spring. physiological processes. We speculate that bud acclimation to previous drought may The impact of elevated air temperature during the growing represent a water-saving strategy. By decreasing the bud season has been proposed to affect spring leaf unfolding by sensitivity to TP, trees delay LU and the associated start of modifying growth cessation and dormancy induction in tempe- evapotranspiration49. The identification of preseason precipita- rate and boreal trees55. Drought stress also affects the onset of tion as an important control of the spatial variance of LU at senescence and possibly dormancy56,57. Precipitation has indeed drought-prone site is in line with this hypothesis and might also been shown to play a large role in vegetation dynamics during the reflect a safety mechanism by which plant delay leaf unfolding senescence period of deciduous forests in the Northern Hemi- until water is available for the restart of plant activity. As a long- sphere58. A delayed dormancy induction can translate into a term strategy, a delay in evapotranspiration will lead to a slower delayed leaf unfolding during the next spring via the later start depletion of water resources at the beginning of the growing date of chilling and GDD accumulation59. Our results are fully in season and reduces the risk of water stress during late spring and line with these observations and provide additional evidence that summer when radiation is more favorable for photosynthesis. the timing of phenological events is impacted by previous phases in the annual cycle of trees. Species show different responses to long-term constraints. All We argue that tree seasonality and long-term biogeographical species exhibited the same response to temperature, however a few constraints are too often overlooked and should be taken into NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 account in phenological studies. Further research on the effect of ecosystems, and over which we are more likely able to gather information about tree seasonality on inter-annual variability of phenology is needed water and temperature pressures of each site. As for LU, extreme climatic years were excluded as we seek to estimate the to clearly identify the role played by biogeographical constraints. average response of the vegetation. We then analyzed the spatial relationships In conclusion, assessing the long-term spatial variance of LU between long-term climatic variables (averaged over the different periods of the and GDDreq is a step in developing a unified framework that will year), LU and GDD. We proceeded in four steps: allow an understanding of the multiple control of climate on 1. We first assessed potential collinearity issues between variables by examining plant phenology. Future research on the importance of plant pairwise correlation coefficients62 (Supplementary Fig. 6); 2. We then selected relevant predictors of observed LU and GDDreq using phenology on ecosystem functioning should focus on space-time penalized elastic net regressions (glmnet function from the glmnet63 package) in interactions with environmental conditions specifically to combination with collinearity information from step 1). Interaction terms were not address: 1) the effects of light and aridity on bud sensitivity to included in the analysis. In step 2–4) predictors were all standardized in order to temperature, and 2) the potential coordination between plant represent the relative contribution of each variable in explaining LU and GDD variability (Supplementary Tables 8 and 9). processes and phenology that could account for a co-limitation by 3. After selecting relevant variables in step 2, we assessed the remaining temperature and the availability of light and water. In line with collinearity between variables by estimating their respective variance inflation these recommendations, the use of current, constant, GDDreq and factor (VIF) as: NCD metrics for the study of spatio-temporal patterns in plant 1 VIFj ¼ ð3Þ phenology should be used with caution. 1  R2j where the VIF for variable j is the reciprocal of the inverse of R² from the Methods regression. VIF values increase with collinearity, and arbitrary threshold of 5–10 Datasets. Data for in situ leaf unfolding (LU) were obtained from the Pan Eur- are commonly used to define high VIF values. Here, when two predictors exhibited opean Phenology (PEP) network (www.pep725.com) and the GDR2968 database potential collinearity (i.e., high VIF values), we removed the one with a VIF value (http://www.gdr2968.cnrs.fr/) for France. Phenological observations followed the higher than 4 and the lowest correlation coefficient with LU or GDD using the Biologische Bundesanstalt, Bundessortenamt und Chemische Industrie (BBCH) stepVIF function from the VIF64 package. code, with LU corresponding to BBCH = 11. 4. We finally assessed the spatial structure of residuals of the reduced model Daily temperature, precipitation and incoming radiation (both shortwave from step 3) with semi-variograms (Supplementary Figs 7–9). We performed [visible and near infrared] and longwave [infrared]) were retrieved from the CRU- generalized least square regressions taking into account the spatial structure of NCEP climatic data set30 at a spatial resolution of 0.5°. Day length was calculated residuals to correct coefficients using the gls function from the nlme65 package. for each site using the ‘geosphere’ R package60. Monthly SM (10–40 cm) was Different spatial structures were tested (linear, exponential, spherical, Gaussian and retrieved from the land data assimilation systems (GLDAS) data set (https://ldas. rational quadratic) and the best model was selected using AIC criterion gsfc.nasa.gov/gldas/) aggregated at a spatial resolution of 0.5° to be consistent with (Supplementary Table 10). the CRU-NCEP data. All the above steps were applied to each species separately. Because results were We approximated drought conditions using the ratio of actual to potential consistent between species (Fig. 2) regarding aggregation uncertainties and because evapotranspiration as a proxy for drought periods (αE). This ratio accurately we aimed at exploring biogeographical constraints at the regional scale for all represents droughts61 and was calculated using the SPLASH model31. αE was deciduous tree species, we applied the same approach to the full dataset with all calculated daily, but only the growing-season αE was used to calculate long-term species pooled together. Analyses were performed with the R v.3.5 software66. drought stress. Reporting summary. Further information on research design is available in Analyses. We corrected site temperature for altitudinal differences between the the Nature Research Reporting Summary linked to this article. site and the mean elevation of each grid cell of the CRU-NCEP data set14 using a gradient of 6.4 °C km−1 and then estimated the heat requirement and associated chilling for each LU observation. Heat requirement (GDDreq) corresponds to the Data availability sum of mean daily temperatures above a threshold of 5 °C calculated from 1 All phenology data are available at http://www.pep725.eu/ and http://www.gdr2968.cnrs. January to the date of LU and was calculated for each observation as: fr. CRU-NCEP data can be downloaded at https://rda.ucar.edu/datasets/ds314.3/. The SPLASH model used to estimate evapotranspiration can be downloaded at https:// LU ¼ dðGDDd ¼ GDDreq Þ ð1Þ bitbucket.org/labprentice/splash/. Soil moisture data can be downloaded at https://ldas. gsfc.nasa.gov/gldas/. The source data underlying Figs 1, 3, 4 and Table 1 and Supplementary Figs 2 and 4 are provided as a Source Data file. X d GDDd ¼  d  Tth Þ; 0Þ maxððT ð2Þ t0 Code availability GDDreq defines the heat requirement of buds at the observed LU date, estimated No custom code nor mathematical algorithms were developed for this study. Only as the accumulated daily air temperature (GDDd), where T d is the mean daily air existing packages and software were used for the analysis, which can be found in the temperature, Tth the temperature threshold for GDD accumulation (5 °C) and to Methods. the starting date (1 January).The number of chilling days (NCD) was calculated as the number of days from 1 November to the LU date with mean daily temperatures between 0 and 5 °C. Received: 18 February 2019; Accepted: 4 November 2019; We were interested in the spatial variance of LU and the heat requirement, so we used the median LU date and the corresponding median GDD for 1970–2016 for each site, assuming that the medians represented the long-term “optimal” LU date and GDD for the background climatic and soil conditions. We only retained sites with more than 5 years of observations and removed years with a LU date outside two interquartiles around the median distribution (i.e., outside days 80–152), which potentially represent a response to extreme events. The same References 1. Peñuelas, J. & Filella, I. Responses to a warming world. Science 294, 793–795 analysis was performed using sites with more than 10 years of observations and led (2001). to similar results. We analyzed the correlations between biogeographical variables, 2. Schwartz, M. D. & others. Phenology: an integrative environmental science. LU dates and GDD for eight species of dominant European deciduous trees with (Springer 2003). many records in the PEP725 database: Aesculus hippocastanum, Alnus glutinosa, Betula pendula, Fagus sylvatica, Fraxinus excelsior, Quercus robur, Sorbus 3. Badeck, F.-W. et al. Responses of spring phenology to climate change. New aucuparia and Tilia cordata. Phytol. 162, 295–309 (2004). Long-term climatic variables (i.e., over 1970–2016) were averaged for two 4. Cleland, E. E., Chuine, I., Menzel, A., Mooney, H. A. & Schwartz, M. D. periods of the year: the growing season, defined as the period between days 180 and Shifting plant phenology in response to global change. Trends Ecol. Evol. 22, 250, and the preseason, defined as the three months before LU. Since we are 357–365 (2007). interested in biogeographical constraints, not the temporality of processes, and 5. Polgar, C. A. & Primack, R. B. Leaf-out phenology of temperate woody plants: because no information about the length of the growing season was available in the from trees to ecosystems. New Phytol. 191, 926–941 (2011). PEP data set, selecting a constant summer period ensures that we have a 6. Tang, J. et al. Emerging opportunities and challenges in phenology: a review. representative period for all trees and all years. It allowed a consistent comparison Ecosphere 7, (2016). between sites without introducing bias induced by different growing season 7. Hänninen, H. & Kramer, K. A framework for modelling the annual cycle of lengths. This period also corresponds to the peak of plant activity in temperate trees in boreal and temperate regions. Silva Fenn. 41, 167–205 (2007). 8 NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 ARTICLE 8. Richardson, A. D. et al. Climate change, phenology, and phenological control 36. Zohner, C. M., Benito, B. M., Svenning, J.-C. & Renner, S. S. Day length of vegetation feedbacks to the climate system. Agric. Meteorol. 169, 156–173 unlikely to constrain climate-driven shifts in leaf-out times of northern woody (2013). plants. Nat. Clim. Change 6, 1120 (2016). 9. Fu, Y. H. et al. Increased heat requirement for leaf flushing in temperate 37. Piao, S. et al. Leaf onset in the northern hemisphere triggered by daytime woody species over 1980–2012: effects of chilling, precipitation and insolation. temperature. Nat. Commun. 6, 6911 (2015). Glob. Change Biol. 21, 2687–2697 (2015). 38. Fu, Y. H. et al. Three times greater weight of daytime than of night-time 10. Zhang, X., Friedl, M. A., Schaaf, C. B. & Strahler, A. H. Climate controls on temperature on leaf unfolding phenology in temperate trees. New Phytol. 212, vegetation phenological patterns in northern mid-and high latitudes inferred 590–597 (2016). from MODIS data. Glob. Change Biol. 10, 1133–1145 (2004). 39. Leuchner, M., Menzel, A. & Werner, H. Quantifying the relationship between 11. Fu, Y. H. et al. Daylength helps temperate deciduous trees to leaf-out at the light quality and light availability at different phenological stages within a optimal time. Glob. Change Biol. (2019). mature mixed forest. Agric. Meteorol. 142, 35–44 (2007). 12. Peñuelas, J. et al. Complex spatiotemporal phenological shifts as a response to 40. Zhong, S. et al. A molecular framework of light-controlled phytohormone rainfall changes. New Phytol. 161, 837–846 (2004). action in Arabidopsis. Curr. Biol. 22, 1530–1535 (2012). 13. Fu, Y. H. et al. Unexpected role of winter precipitation in determining heat 41. Eckstein, A., Krzeszowiec, W., Banaś, A. K., Janowiak, F. & Gabryś, H. requirement for spring vegetation green-up at northern middle and high Abscisic acid and blue light signaling pathways in chloroplast movements in latitudes. Glob. Change Biol. 20, 3743–3755 (2014). Arabidopsis mesophyll. Acta Biochim. Pol. 63, 449–458 (2016). 14. Olsson, C. & Jönsson, A. M. Budburst model performance: The effect of the 42. Singh, R. K., Svystun, T., AlDahmash, B., Jönsson, A. M. & Bhalerao, R. P. spatial resolution of temperature data sets. Agric. Meteorol. 200, 302–312 Photoperiod-and temperature-mediated control of phenology in trees–a (2015). molecular perspective. New Phytol. 213, 511–524 (2017). 15. Botta, A., Viovy, N., Ciais, P., Friedlingstein, P. & Monfray, P. A global 43. Zohner, C. M. & Renner, S. S. Perception of photoperiod in individual buds of prognostic scheme of leaf onset using satellite data. Glob. Change Biol. 6, mature trees regulates leaf-out. New Phytol. 208, 1023–1030 (2015). 709–725 (2000). 44. Brelsford, C. C., Nybakken, L., Kotilainen, T. K. & Robson, T. M. The 16. Migliavacca, M. et al. On the uncertainty of phenological responses to climate influence of spectral composition on spring and autumn phenology in trees. change, and implications for a terrestrial biosphere model. Biogeosciences 9, Tree Physiol. 39, 925–950 (2019). 2063–2083 (2012). 45. Strømme, C. et al. UV-B and temperature enhancement affect spring and 17. Richardson, A. D. et al. Terrestrial biosphere models need better autumn phenology in P opulus tremula. Plant Cell Environ. 38, 867–877 (2015). representation of vegetation phenology: results from the North American 46. Cochard, H., Coste, S., Chanson, B., Guehl, J. M. & Nicolini, E. Hydraulic Carbon Program Site Synthesis. Glob. Change Biol. 18, 566–584 (2012). architecture correlates with bud organogenesis and primary shoot growth in 18. Garonna, I., de Jong, R. & Schaepman, M. E. Variability and evolution of beech (Fagus sylvatica). Tree Physiol. 25, 1545–1552 (2005). global land surface phenology over the past three decades (1982–2012). Glob. 47. Le Dantec, V., Dufrêne, E. & Saugier, B. Interannual and spatial variation in Change Biol. 22, 1456–1468 (2016). maximum leaf area index of temperate deciduous stands. Ecol. Manag. 134, 19. Huang, M. et al. Velocity of change in vegetation productivity over northern 71–81 (2000). high latitudes. Nat. Ecol. Evol. 1, 1649 (2017). 48. Bréda, N., Huc, R., Granier, A. & Dreyer, E. Temperate forest trees and stands 20. Zhang, X., Tarpley, D. & Sullivan, J. T. Diverse responses of vegetation under severe drought: a review of ecophysiological responses, adaptation phenology to a warming climate. Geophys. Res. Lett. 34, L19405 (2007). processes and long-term consequences. Ann. Sci. 63, 625–644 (2006). 21. Prieto, P. et al. Changes in the onset of spring growth in shrubland species in 49. Gonsamo, A., Ter-Mikaelian, M. T., Chen, J. M. & Chen, J. Does earlier and response to experimental warming along a north–south gradient in Europe. increased spring plant growth lead to reduced summer soil moisture and plant Glob. Ecol. Biogeogr. 18, 473–484 (2009). growth on landscapes typical of Tundra-Taiga interface? Remote Sens. 11, 22. Čufar, K., De Luis, M., Saz, M. A., Črepinšek, Z. & Kajfež-Bogataj, L. 1989 (2019). Temporal shifts in leaf phenology of beech (Fagus sylvatica) depend on 50. Radoglou, K., Dobrowolska, D., Spyroglou, G. & Nicolescu, V.-N. A review on elevation. Trees 26, 1091–1100 (2012). the ecology and silviculture of limes (Tilia cordata Mill., Tilia platyphyllos 23. Pellerin, M., Delestrade, A., Mathieu, G., Rigault, O. & Yoccoz, N. G. Spring Scop. and Tilia tomentosa Moench.) in Europe. Romania 15, 16 (2008). tree phenology in the Alps: effects of air temperature, altitude and local 51. Scherrer, D., Bader, M. K.-F. & Körner, C. Drought-sensitivity ranking of topography. Eur. J. Res. 131, 1957–1965 (2012). deciduous tree species based on thermal imaging of forest canopies. Agric. 24. Dantec, C. F. et al. Escape of spring frost and disease through phenological Meteorol. 151, 1632–1640 (2011). variations in oak populations along elevation gradients. J. Ecol. 103, 52. Fu, Y. H. et al. Declining global warming effects on the phenology of spring 1044–1056 (2015). leaf unfolding. Nature 526, 104 (2015). 25. Schleip, C., Sparks, T. H., Estrella, N. & Menzel, A. Spatial variation in 53. Lapenis, A., Henry, H., Vuille, M. & Mower, J. Climatic factors controlling onset dates and trends in phenology across Europe. Clim. Res. 39, 249–260 plant sensitivity to warming. Clim. Change 122, 723–734 (2014). (2009). 54. Kramer, K. Phenotypic plasticity of the phenology of seven European tree 26. Seyednasrollah, B., Swenson, J. J., Domec, J.-C. & Clark, J. S. Leaf phenology species in relation to climatic warming. Plant Cell Environ. 18, 93–104 (1995). paradox: Why warming matters most where it is already warm. Remote Sens. 55. Hänninen, H. & Tanino, K. Tree seasonality in a warming climate. Trends Environ. 209, 446–455 (2018). Plant. Sci. 16, 412–416 (2011). 27. Crimmins, T. M., Crimmins, M. A. & Bertelsen, C. D. Complex responses to 56. Xie, Y., Wang, X. & Silander, J. A. Deciduous forest responses to temperature, climate drivers in onset of spring flowering across a semi-arid elevation precipitation, and drought imply complex climate change impacts. Proc. Natl gradient. J. Ecol. 98, 1042–1051 (2010). Acad. Sci. USA 112, 13585–13590 (2015). 28. Vitasse, Y. & Basler, D. What role for photoperiod in the bud burst phenology 57. Marchin, R., Zeng, H. & Hoffmann, W. Drought-deciduous behavior reduces of European beech. Eur. J. Res. 132, 1–8 (2013). nutrient losses from temperate deciduous trees under severe drought. 29. Way, D. A. & Montgomery, R. A. Photoperiod constraints on tree phenology, Oecologia 163, 845–854 (2010). performance and migration in a warming world. Plant Cell Environ. 38, 58. Papagiannopoulou, C. et al. Vegetation anomalies caused by antecedent 1725–1736 (2015). precipitation in most of the world. Environ. Res. Lett. 12, 074016 (2017). 30. Viovy, N. CRUNCEP Version 7 - Atmospheric Forcing Data for the 59. Fu, Y. S. et al. Variation in leaf flushing date influences autumnal senescence Community Land Model (2018). and next year’s flushing date in two temperate tree species. Proc. Natl Acad. 31. Davis, T. et al. Simple process-led algorithms for simulating habitats (SPLASH Sci. USA 111, 7355–7360 (2014). v. 1.0): robust indices of radiation, evapotranspiration and plant-available 60. Hijmans, R. J., Williams, E. & Vennes, C. Geosphere: spherical trigonometry. R moisture. Geosci. Model Dev. 10, 689–708 (2017). package version 1.3-11 (2014). 32. Liu, G., Chen, X., Zhang, Q., Lang, W. & Delpierre, N. Antagonistic effects of 61. Stocker, B. D. et al. Quantifying soil moisture impacts on light use efficiency growing season and autumn temperatures on the timing of leaf coloration in across biomes. New Phytol. 218, 1430–1449 (2018). winter deciduous trees. Glob. Change Biol. 24, 3537–3545 (2018). 62. Dormann, C. F. et al. Collinearity: a review of methods to deal with it and a 33. Jenkins, J., Braswell, B., Frolking, S. & Aber, J. Detecting and predicting spatial simulation study evaluating their performance. Ecography 36, 27–46 (2013). and interannual patterns of temperate forest springtime phenology in the 63. Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for generalized eastern US. Geophys. Res. Lett. 29, 54–1 (2002). linear models via coordinate descent. J. Stat. Softw. 33, 1 (2010). 34. Maignan, F., Bréon, F., Vermote, E., Ciais, P. & Viovy, N. Mild winter and 64. Lin, D., Foster, D. P. & Ungar, L. H. VIF regression: a fast regression spring 2007 over western Europe led to a widespread early vegetation onset. algorithm for large data. J. Am. Stat. Assoc. 106, 232–247 (2011). Geophys. Res. Lett. 35 (2008). 65. Pinheiro, J. et al. nlme: Linear and nonlinear mixed effects models. R Package 35. Peaucelle, M. et al. Representing explicit budburst and senescence Version 3 (2012). processes for evergreen conifers in global models. Agric. Meteorol. 266, 97–108 66. R Core Team, D. R: A Language and Environment for Statistical Computing (2019). (2016). NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-13365-1 Acknowledgements Peer review information Nature Communications thanks Renee Marchin and the We acknowledge the financial support from the European Research Council Synergy anonymous reviewers for their contribution to the peer review of this work. Peer reviewer grant ERC-SyG-2013–610028 IMBALANCE-P. We acknowledge all members of reports are available. the PEP725 project for providing the phenological data. B.D.S. was funded by the Swiss National Science Foundation grant no.PCEFP2_181115. M.P. acknowledges Reprints and permission information is available at http://www.nature.com/reprints financial support from the Fonds WetenschappelijkOnderzoek (FWO) grant G018319N. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Author contributions M.P., I.J. and J.P. designed the study. M.P. performed the analyses and wrote the first Open Access This article is licensed under a Creative Commons Attri- version of the manuscript. B.D.S. provided evaporation data. R.M.H. contributed to statistical analysis. M.P., I.J., B.D.S., A.D.F., Y.H.F., R.M.H., P.C. and J.P. contributed to bution 4.0 International License, which permits use, sharing, adaptation, the interpretation of the results and to revisions of the manuscript. distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this Competing interests article are included in the article’s Creative Commons license, unless indicated otherwise in The authors declare no competing interests. a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the Additional information permitted use, you will need to obtain permission directly from the copyright holder. To Supplementary information is available for this paper at https://doi.org/10.1038/s41467- view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. 019-13365-1. Correspondence and requests for materials should be addressed to M.P. © The Author(s) 2019 10 NATURE COMMUNICATIONS | (2019)10:5388 | https://doi.org/10.1038/s41467-019-13365-1 | www.nature.com/naturecommunications