13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Received: 11 April 2022 | Accepted: 26 September 2022 DOI: 10.1111/gcb.16502 RESEARCH ARTICLE Radiation-­constrained boundaries cause nonuniform responses of the carbon uptake phenology to climatic warming in the Northern Hemisphere Adrià Descals1,2 | Aleixandre Verger1,2,3 | Gaofei Yin1,2,4 | Iolanda Filella1,2 | Yongshuo H. Fu5 | Shilong Piao6 | Ivan A. Janssens7 | Josep Peñuelas1,2 1 CREAF, Cerdanyola del Vallès, Barcelona, Spain Abstract 2 CSIC, Global Ecology Unit CREAF-­C SIC-­ Climatic warming has lengthened the photosynthetically active season in recent dec- UAB, Barcelona, Spain 3 ades, thus affecting the functioning and biogeochemistry of ecosystems, the global CIDE, CSIC-­UV-­GV, València, Spain 4 carbon cycle and climate. Temperature response of carbon uptake phenology var- Faculty of Geosciences and Environmental Engineering, Southwest ies spatially and temporally, even within species, and daily total intensity of radiation Jiaotong University, Chengdu, China 5 may play a role. We empirically modelled the thresholds of temperature and radiation College of Water Sciences, Beijing Normal University, Beijing, China under which daily carbon uptake is constrained in the temperate and cold regions 6 College of Urban and Environmental of the Northern Hemisphere, which include temperate forests, boreal forests, alpine Sciences, Peking University, Beijing, China and tundra biomes. The two-­dimensionality of the temperature-­radiation constraint 7 Department of Biology, University of Antwerp, Wilrijk, Belgium was reduced to one single variable, θ, which represents the angle in a polar coordi- nate system for the temperature-­radiation observations during the start and end of Correspondence Adrià Descals, CREAF, Cerdanyola del the growing season. We found that radiation will constrain the trend towards longer Vallès, Barcelona 08193, Catalonia, Spain. growing seasons with future warming but differently during the start and end of sea- Email: a.descals@creaf.uab.cat son and depending on the biome type and region. We revealed that radiation is a Funding information major factor limiting photosynthetic activity that constrains the phenology response Catalan Government, Grant/Award Number: AGAUR-­2020PANDE00117 to temperature during the end-­of-­season. In contrast, the start of the carbon uptake and SGR 2017-­1005; CSIC Thematic is overall highly sensitive to temperature but not constrained by radiation at the hemi- Interdisciplinary Platform Teledetección, Grant/Award Number: PTI-­TELEDETECT; spheric scale. This study thus revealed that while at the end-­of-­season the phenology European Research Council Synergy response to warming is constrained at the hemispheric scale, at the start-­of-­season project, Grant/Award Number: SyG-­ 2013-­610028; Fundación Ramon Areces, the advance of spring onset may continue, even if it is at a slower pace. Grant/Award Number: ELEMENTAL-­ CLIMATE; National Funds of China for KEYWORDS Distinguished Young Scholar, Grant/ carbon uptake phenology, climatic warming, phenology modelling, photosynthetically growing Award Number: 42025101; Spanish season, radiation constraints, vegetation phenology Government, Grant/Award Number: PID2019-­110521GB-­I 00 1 | I NTRO D U C TI O N warming, the growing season has lengthened, with the start of the growing season advancing more than the end of the growing season Seasonal changes in weather conditions drive the timing of the start delaying (Menzel et al., 2006; Peñuelas & Filella, 2001). Vegetation and end of vegetation growth (Piao et al., 2019). With recent climatic phenology influences the carbon and water cycle (Keenan et al., 2014; This is an open access article under the terms of the Creative Commons Attribution-­NonCommercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes. © 2022 The Authors. Global Change Biology published by John Wiley & Sons Ltd. Glob Change Biol. 2023;29:719–730.  wileyonlinelibrary.com/journal/gcb | 719 | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 720 DESCALS et al. Myneni et al., 1997; Richardson et al., 2010), successional transition Parazoo, et al., 2020). The start and end of the season occur when (Chuine, 2010) and plays a role in the feedback between vegetation a limiting factor, such as air temperature and radiation intensity ex- and climate (Peñuelas & Filella, 2009; Richardson et al., 2013). The ceeds the threshold under which photosynthesis is constrained. warming-­induced lengthening of the growing season has increased One limiting factor of the phenology of photosynthesis is the the carbon uptake (Le Quéré et al., 2009), offsetting atmospheric low light availability (Zhang, Commane, et al., 2020). Photosynthetic carbon from human emissions. However, how phenology may re- activity as a function of radiation, with constant temperature and spond to future warming and whether vegetation will increase car- soil water content conditions, is equivalent to light-­response curves bon sequestration remains unclear (Piao et al., 2019), which adds (Herrmann et al., 2020), which represent the restrictions of photo- uncertainty to future atmospheric carbon concentration and, thus, synthetic activity by low levels of radiation, even if temperature, soil climate projections. water content, and other factors influencing photosynthesis are fa- Winter chilling, forcing requirements, and photoperiod trigger vourable. The restriction is due to the inoperability of chloroplasts the start of foliar phenology (Richardson et al., 2013). Previous re- and light-­dependent reactions in the absence of photosynthetically search suggests that chilling requirements will constrain the advance active radiation (Johnson et al., 2008). of leaf unfolding in deciduous forests by reducing temperature sen- Here, we aimed to clarify the role of temperature and radiation sitivity (Fu et al., 2015), whereas other factors such as precipitation constraints during the start and end of the photosynthetically ac- (Peaucelle et al., 2019) and photoperiod (Körner & Basler, 2010; tive season and how these constraints change among ecosystems Meng et al., 2021; Zohner et al., 2016) may increase the heat require- and over space. To achieve this, we proposed a phenology model ments during the ecodormancy stage and, as a result, slow down based on the law of limiting factors that predicts the start and end the warming-­induced advance of the leaf unfolding. With regards of the photosynthetically growing season. We evaluated the current to leaf senescence, a recent study found that increased productivity temperature and radiation constraints on the phenology of tem- during the growing season counteracts the warming-­induced delay perate and cold regions of the Northern Hemisphere, where tem- in leaf senescence (Zani et al., 2020). These findings are based on perature and radiation mostly regulated the vegetation phenology the study of leaf phenophases in deciduous forests obtained from in (Piao et al., 2019). Lastly, we estimated the potential lengthening of situ observations (Templ et al., 2018) or remotely-­sensed vegetation the growing season in a warming scenario in which radiation might greenness indices (Julien & Sobrino, 2009), such as the normalized take over the role as a limiting factor of photosynthetic activity vegetation difference index (NDVI). phenology. The key conceptual framework in these studies is that conditions preceding the phenological dates affect the timing of tree pheno- phases (Chuine et al., 2013). Phenology modelling has largely ac- 2 | M ATE R I A L S A N D M E TH O DS counted for chilling and forcing requirements in decisuous forests (Fu et al., 2020). However, the start and end of the photosyntheti- 2.1 | Data cally active season are directly influenced by current meteorological conditions, not by conditions preceding the phenological dates. This We used daily in situ records from 63 sites of the FLUXNET2015 is evidenced by a decoupling between remotely sensed vegetation Tier 1 data set (Pastorello et al., 2020) (see Figure S1 for the loca- greenness and proxies of photosynthetic activity (Jeong et al., 2017; tions of the sites). The biome types at the FLUXNET sites were Yin et al., 2020; Zhang, Commane, et al., 2020), which indicates that characterized using the MODIS MCD12Q1 V6 product (Friedl & vegetation might present leaves but these are not photosynthetically Sulla-­Menashe, 2015) for the mode of land cover during 2001–­2019 active because meteorological conditions are restricting photosyn- and the RESOLVE Ecoregions 2017 map (Dinerstein et al., 2017). We thetic activity at the moment. For instance, vegetation productivity selected the sites with at least 4 years of record and categorized as de- declines in agreement with the decrease in radiation intensity during ciduous broadleaved forests (DBF; 12 sites), evergreen needleleaved autumn (Zhang, Commane, et al., 2020), and a lack of available light forests (ENF; 25 sites), mixed forests (MX; 7 sites), and grasslands and may, thus, prompt the end of the growing season. Similarly, water wetlands (GRA and WET; 19 sites) in the tundra, boreal, and temper- availability determines the start and end of the photosynthetically ate biomes of the Northern Hemisphere according to the RESOLVE active season in tropical dryland ecosystems (Eamus & Prior, 2001; Ecoregions 2017 map. These FLUXNET sites represent 643 site-­years Jiao et al., 2021; Zhang, Parazoo, et al., 2020). and are located >30° N. The sites cover different range of years from Given the link between photosynthesis and current meteorolog- 1991 to 2014. We used daily gap-­filled GPP reference from GPP ver- ical conditions, the seasonality of carbon uptake can be modelled sions using model efficiency (MEF) and obtained from the day-­time with meteorological variables, and the start and end of the growing partitioning method (Pastorello et al., 2020). season can be determined by constraint functions (Jolly et al., 2005). For the climatic data, we used daily averages of air temperature The constraint functions define the photosynthesis-­inhibiting and incoming shortwave radiation at 0.1 arc degrees from the ERA5-­ thresholds that restrict the carbon uptake during the growing sea- Land hourly data (Muñoz-­Sabater et al., 2021). The phenology dates son. Using this phenology modelling framework, carbon uptake phe- were extracted from the 4-­day clear-­sky daily contiguous solar-­ nology can be, thus, determined by the most limiting factor (Zhang, induced chlorophyll fluorescence (SIF) estimates at 0.05° from the | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License DESCALS et al. 721 OCO-­2 (Orbiting Carbon Observatory-­2) (Zhang et al., 2018) in tem- model first generates a binary time series where daily observations perate and cold regions of the Northern Hemisphere for the 2001–­ are categorized as 0 (vegetation productivity is inhibited or con- 2020 period. SIF observations have been demonstrated to track the strained) and 1 (vegetation productivity is active or unconstrained seasonality of GPP globally (Sun et al., 2017). The temperate and by temperature and radiation). In the second step of the model, the cold regions include tundra, boreal forests, temperate broadleaf and start and end of season dates are extracted from the binary time mixed forests, and temperate coniferous forests in the RESOLVE series using the maximum-­separation method. We named the model Ecoregions 2017 map (Figure S1). as climate-­constraint model because climate time series are used to generate a binary time series of constrained and unconstrained daily observations. The binary time series is analogous to the growing-­ 2.2 | Phenology extraction season index (Jolly et al., 2005), which is generated from cut-­off functions from daily temperature, daylength, and vapor-­pressure We estimated the Start of Season (SoS) and End of Season (EoS) deficit (VPD), and ranges from 0 and 1 (0 = vegetation is constrained using the maximum-­separation method (Descals et al., 2020). The and 1 = vegetation is unconstrained). maximum-­separation method is a threshold-­based method that can The modelled SoS and EoS were compared with estimated SoS effectively estimate phenological metrics without the need of time and EoS extracted from in situ GPP time series at the FLUXNET series pre-­processing, which tends to distort the seasonality of the sites. The performance of the model was evaluated in terms of time series. The algorithm first runs a moving window that calculates root-­mean-­squared-­error (RMSE). As the growing-­season index, the proportion of observations that are above a given threshold th the climate-­constrained model assumes no forcing requirements before the central day of the moving window (pbefore) and after the for the SoS and EoS; the growing season starts as soon as cli- central day (pafter). A moving window greater than 90 days and lower mate conditions are favorable for vegetation growth, respectively. than the length of the growing season is advised for the method to To test this assumption, we compared the performance of the be effective (Descals et al., 2020). In this study, the length of the climate-­constrained model with the photo-­threshold model (Meng moving window was 120 days. Then, the algorithm computes the et al., 2021), a model based on the growing-­degree-­day model. We difference between the proportions (pbefore -­ pafter) for each day of also compared the performance of the SoS and EoS extracted from year in the time series. The SoS is calculated as the date when the the growing-­season index. difference in proportions (pbefore -­pafter) reaches a minimum, and the EoS is calculated as the date when the difference is maximal. The method was used to extract the SoS and EoS from the FLUXNET 2.4 | Climate-­constraint model GPP time series and in the climate-­constraint model and growing-­ season index. In this study, we used a dynamic threshold for the The approach for phenological modelling categorizes daily time extraction of SoS and EoS from the FLUXNET GPP and the growing-­ series of temperature T(t) (°C) and incoming shortwave radia- season index. The dynamic threshold th is computed as a percentage tion R(t) (W m−2) into a binary time series B(t) with values 0 (con- of the amplitude (Equation 1): strained) and 1 (unconstrained). The classification was done with a threshold function (Equation 2); vegetation was considered un- ( ) th = Imax − Imin × P + Imin (1) constrained when daily temperature and radiation were above the threshold function. The threshold function has the form of The Imax and Imin are the maximum and minimum values of the time a rational function with one asymptote for temperature (Tthresh) series (e.g. FLUXNET GPP or the growing-­season index). The param- and another for radiation (R thresh) (Figure 1a). The curvature of the eter P ranges from 0 to 1. We set P equal to 0.2, which represents a rational function is defined by parameter C (unitless). The param- percentage of 20% of the amplitude. This low percentage ensures the eter C in Equation (2) reflects the covariation of temperature and reliable detection of the first and last stages of the SoS and EoS. The radiation constraints. Instead of a constant threshold for tem- maximum-­separation method can be applied to binary time series, perature and radiation, the parameter C reflects the curvature of which is the case of the climate-­constraint model (See Section 2.4. the threshold function and provides a threshold that depends on Climate-­constraint model). In a binary time series in which values are 0 the levels of temperature and light availability (i.e., high levels of or 1, the SoS and EoS would be the same for any percentage P greater temperature reduce the radiation constraint). We used a rational than 0 and lower than 1. For this reason, we used a constant thresh- function owing to the relationship of FLUXNET GPP with air tem- old th equal to 0.5 for the climate-­constraint model. perature and shortwave radiation (Figure 1b). ⎧ 1 C 2.3 | Phenology modelling ⎪ T(t) − Tthresh + Rthresh − R(t) < 0 ⎪ C (2) B(t) = ⎨ 0 + Rthresh − R(t) ≥ 0 ⎪ T(t) − Tthresh We proposed a phenology model that incorporates the constraints ⎪ 0 T(t) < Tthresh of temperature and radiation on the vegetation productivity. The ⎩ | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 722 DESCALS et al. F I G U R E 1 Graphical representation of the calculation of parameters in the climate-­constraint model. (a) Graphical representation of the threshold function (black line) and the angle θ. The threshold function consists of a rational function defined by the asymptotes in Tthresh and Rthersh. The angle θpheno represents the radiation constraint of temperature and radiation during a phenological event. (b) Mean gross primary production (GPP) as a function of air temperature and incoming shortwave radiation for all observations of the FLUXNET sites in the cold and temperate regions of the Northern Hemisphere. The time series of FLUXNET GPP and ERA5-­L and air tempera- temperature (Equation 3; in °C), daylength (Equation 4; in hours), ture and shortwave radiation were used to calibrate the parame- and VPD (Equation 5, in kPa), which bound the climate time series ters Tthresh , R thresh , and C. To accomplish this, the FLUXNET GPP to values between 0 and 1. The cut-­off functions are defined by a was converted to a binary time series: 0 when GPP was < GPP thresh minimum and a maximum threshold for temperature (Tmin and Tmax), −2 −1 and 1 when GPP was > GPP thresh C m d . We used a GPP thresh daylength (DLmin and DLmax), and VPD (VPDmin and VPDmax). The equal to 2 C m−2 d−1, which is approximately the median GPP of growing-­season index is the multiplication of Tgsi, DL gsi, and VPDgsi the 192,709 observations and represents a low GPP threshold (Equation 6) and ranges from 0 (vegetation is fully constrained) to 1 that ensure the first and last stages of the SoS and EoS, respec- (vegetation is fully unconstrained). tively. The calibration was done using the simplex search method (Lagarias et al., 1998). ⎧ ⎪ 0 T(t) ≤ Tmin The second step of the model consists of the extraction of phe- ⎪ T(t) − Tmin Tgsi (t) = ⎨ Tmin < T(t) < Tmax (3) nological dates (SoS and EoS) from the binary time series B(t) using ⎪ Tmax − Tmin the maximum-­separation method (Descals et al., 2020). The binary ⎪ 1 T(t) ≥ Tmax ⎩ time series B(t) reflects that vegetation is unconstrained whenever temperature and radiation conditions are favorable for vegetation growth. Thus, the binary time series may reflect unconstrained ⎧ ⎪ 0 DL(t) ≤ DLmin observations when sporadic favorable conditions occur during the ⎪ DL(t) − DLmin DLgsi (t) = ⎨ DLmin < DL(t) < DLmax (4) dormant period, and the SoS could incorrectly identified in such ⎪ DLmax − DLmin scenario. To avoid this shortcoming, Jolly et al. (2005) proposed a ⎪ 1 DL(t) ≥ DLmax ⎩ 21-­day moving average of the growing-­season index to remove ex- treme events that can prematurely detect the start of the growing season. In this study, we used the maximum-­separation method be- ⎧ ⎪ 0 VPD(t) ≤ VPDmin cause the method detects the SoS only when there are many values ⎪ VPD(t) − VPDmin VPDgsi (t) = ⎨ VPDmin < VPD(t) < VPDmax (5) that are unconstrained after the dormant period. Similarly, the EoS ⎪ VPDmax − VPDmin is detected when there is a substantial decline in the number of un- ⎪ 1 VPD(t) ≥ VPDmax ⎩ constrained observations after the growing season. Thus, the daily observations of favorable temperature and radiation that may occur GSI(t) = Tgsi (t) × DLgsi (t) × VPDgsi (t) (6) sporadically during the dormant period are not detected as SoS or EoS by the maximum-­separation method. Equation (3) originally used daily minimum air temperature but, for consistency with the climate-­constrained function, we used daily mean 2.5 | Growing-­season index air temperature. The daylength was calculated as a function of time (day of the year) and latitude (°) using the model proposed in (Forsythe The growing-­season index was designed to reproduce the NDVI time et al., 1995). VPD was calculated using the formula in Yuan et al. (2019), series from air temperature, daylength, and VPD (Jolly et al., 2005). which uses air and dew point temperature, both variables included The growing-­season index is calculated from cut-­off functions for in the ERA5-­Land dataset. Jolly et al. (2005) defined the threshold | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License DESCALS et al. 723 parameters of the cut-­off functions. However, we calibrated these pa- and Tthresh) were scaled by dividing by 25°C, and radiation values rameters for our case study because (1) we used different climate data, (Rpheno and Rthresh) were scaled by dividing by 300 W m−2. These scal- (2) we used mean air temperature instead of minimum air temperature, ing values represent approximately the 95th percentile of all tem- and (3) we aimed to reproduce the FLUXNET GPP, representing vege- perature and radiation observations in the FLUXNET sites. tation productivity, instead of the NDVI, which is an indicator of green ( ( )) biomass. The calibration of the parameters was done using the simplex 25 Rpheno − Rthresh 𝜃 = tan −1 ( ) (9) search method. 350 Tpheno − Tthresh The growing-­season index can be used as a phenology model. As the climate-­constrained model, phenological dates can be extracted A high temperature and low amount of radiation at the time of the phe- from the growing-­season index. We used the maximum-­separation nological event lead to low values of θ proximal to 0°, and a low tem- method and a dynamic threshold of 20% of the amplitude to deter- perature and high amount of radiation at the time of the phenological mine the SoS and EoS dates. event lead to high values of θ proximal to 90°. Values of θ <0 or >90° were clamped to 0 and 90°, respectively. 2.6 | Photo-­threshold model 2.8 | Estimation of the potential lengthening of the The photo-­threshold works similarly as the growing-­degree-­day growing season model (Hänninen, 2016) but replaces the arbitrary date when forcing starts accumulating (generally January 1) with the date when day- We converted the current climate constraints on phenology into length reaches the threshold DL start (Meng et al., 2021). The state of potential lengthening of the growing season. The potential SoS and forcing Sf(t) (Equation 7) is calculated as the summation of the rate EoS (SoSpot and EoSpot) represent the dates under a warming sce- of forcing Rf(t) starting from DL start (Equation 8). Temperature T(t) nario in which SoS and EoS would be restricted only by radiation. is offset with a base temperature Tbase (0°C in this study). The SoS This gives a practical quantification of the climate constraints on the occurs the date when the state of forcing Sf(t) is higher than the forc- phenological dates. For instance, a small difference between EoS ing requirement F. If the forcing requirement is not met, the model and EoSpot indicates that future warming will not delay the EoS. For forces the SoS to the first day of year when daylength surpasses the calculation of SoSpot and EoSpot, we changed the threshold func- DLend. In the photo-­threshold model, DL start, DLend, and F are the tion in Equation 2 to a constant radiation threshold equal to Rthresh. parameters to calibrate. The simplex search method was used for The observation was considered inhibited when shortwave radiation the calibration. was Rthresh. Analytically, this would ∑ correspond to the climate-­constraint function (Equation 2) having a Sf(t) = Rf(T(t)) (7) high temperature T(t) that is nearly infinity. t0 ⎧ ⎪ (t) − Tbase T(t) > Tbase 2.9 | Estimation of temperature dependency Rf(t) = ⎨ (8) ⎪ 0 T(t) ≤ Tbase ⎩ The dependency of SoS and EoS on temperature was estimated for 2001–­2020 using ERA5-­L and and the phenological dates extracted 2.7 | Estimation of temperature and radiation from the SIF time series. The SoS and EoS were extracted using a constraints on phenological dates dynamic threshold of 20% the amplitude. The mean air tempera- ture during the phenological dates were calculated as the average To quantify the limitations of radiation and temperature on the start daily air temperature during ±7 days the SoS and EoS. Temperature and end of the growing season, we reduced the two-­dimensionality dependency was estimated using the coefficient of correlation be- of the temperature-­radiation constraint to one single variable; θ. tween temperature and phenological dates obtained from SIF time The angle θ was calculated geometrically (Equation 9) based on the series, which provides more spatial and temporal coverage (2001–­ threshold function (Equation 2) and given the points Pa, Pb, and Pc, 2020) than the FLUXNET GPP records. We used Pearson correla- where Pa and Pb have fixed coordinates in (Tthresh, Rthresh) and (25, tions, which is termed temperature dependency, instead of the Rthresh), respectively. The coordinates of Pc are given by the tem- conventional metric of temperature sensitivity (Fu et al., 2015) be- perature and shortwave radiation at the time of the phenological cause the latter overestimates the sensitivity when the temperature event (Tpheno and Rpheno). We calculated Tpheno and Rpheno as the mean time series is highly variable (Keenan et al., 2020), which could po- temperature and radiation 7 days before and after the phenological tentially weaken the spatial analysis in our study. Pixels were masked event (SoS and EoS) to avoid noisy values of θ due to variation on if minimum winter radiation was below Rthresh; latitudes where win- daily temperature and radiation. A graphical representation of the ter radiation is above Rthresh are potentially not limited by radiation calculation of θ is shown in Figure 1a. The temperature values (Tpheno and present SIF values greater than 0. The minimum winter radiation | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 724 DESCALS et al. was estimated as the average shortwave radiation from day of year function in ENF and MX (Figure S2), although on average the 345 to 365 and aggregated over the 2001–­2020 period. daily mean radiation differed between SoS (air temperature, 7.5°C; shortwave radiation, 212 W m−2) and EoS (air temperature, 7.7°C; shortwave radiation, 94 W m−2) in all vegetation types, 3 | R E S U LT S which suggests a contrasting limitation of temperature and ra- diation during the SoS and EoS in temperate and cold regions of The parameters Tthresh, Rthresh, and C, which define the threshold the Northern Hemisphere. Values of θ, extracted from SIF time function of temperature and radiation, were 2.56°C, 26.5 W m−2, and series, were on average high during the SoS for 2001–­2 020 in the 342 when considering all FLUXNET site-­years (Table 1). The param- cold and temperate regions of the Northern Hemisphere (median eters Tthresh and Rthresh, which define the asymptotes of the rational and 5–­95th percentile range was 87 [65, 98]), which contrasted function, did not change substantially among vegetation types; ENF with the low θ during the EoS (32 [14, 50]) (Figure 3a). Specifically, presented the lowest Tthresh (2.07°C) and highest Rthresh (34.8 W m−2), the regions where radiation limitation was highest during the SoS while DBF presented the highest Tthresh (2.71°C) and lowest Rthresh and EoS correspond to temperate Europe (Figure S3), while tem- −2 (25.9 W m ). Temperature and radiation during the FLUXNET GPP perature restricted the SoS and EoS in high-­latitude and high-­ SoS and EoS occurred closely to the threshold function in ENF and altitude regions. MX (Figure S2). The limitation of radiation to photosynthesis was linked to a low The SoS and EoS were modeled using climatic time series and then sensitivity of phenology to temperature in temperate and cold re- compared to the observed SoS and EoS in the FLUXNET sites. The gions of the Northern Hemisphere. SoS was on average more de- estimated and observed phenological metrics showed the highest pendent than EoS to temperature. The dependency of phenological agreement in ENF (RMSE was 11.78 d and 12.97 d for SoS and EoS, dates to temperature was significantly lower when the limitation respectively) (Figure 2). The model output agreed well with observed of radiation was high (Figure 4). Regions with the highest radiation SoS and EoS in ENF sites, despite the wide range of latitudes and var- limitation during the SoS, where values of θ were lower than the ious biome type; constraints of temperature and radiation explained 5th percentile (θ < 65), presented a significantly lower temperature the SoS and EoS in the FLUXNET site San Rossore (IT-­SRo), with a dependency (R = −0.50) than the rest of the regions (R = −0.63). Mediterranean climate, as well as the FLUXNET site Saskatchewan Similarly, regions with the lowest radiation limitation during the EoS, (Ca-­Obs), covering a boreal forest. The largest RMSE was found for defined as values greater than the 95th percentile (θ > 50), exhibited SoS in deciduous forests (34.38 d), which presented a large bias for the a significantly greater temperature dependency (R = 0.23) than the SoS (ME was 34.4 d in DBF compared to the 0.4 d in ENF). rest of the regions (R = 0.16). Model parameters in the growing-­season index and the photo-­ Lastly, we defined the dates when radiation could potentially threshold model were overall similar. However, the forcing re- become the climatic variable that fully restricted photosynthetic quirement F, in the photo-­threshold model, was a parameter that activity during spring onset and autumnal senescence as SoS pot differed substantially among vegetation types (Table 1). The forc- and EoS pot , representing the days of the year when daily short- ing requirement F ranged from 35.7°C in ENF to 258.4°C in DBF. wave radiation reaches threshold R thresh (27 W m−2) in spring The photo-­threshold method presented the highest RMSE in DBF and autumn. The 27 W m−2 threshold represents the radiation and GRA (Table 2) but, in ENF and MX, the RMSE was similar to threshold if temperature did not constrain photosynthesis. The the climate-­constraint function. The SoS and EoS modeled using the estimated EoS was generally similar to the EoS pot based on the growing-­season index showed a higher RMSE compared with the FLUXNET and gridded climatic data (Figures 3b and 5). The mean climate-­constraint model in all vegetation types. difference EoS pot –­EoS was 23.6 d in the FLUXNET sites. The ob- Temperature and radiation during SoS and EoS at the served SoS estimated using FLUXNET and gridded climatic data, FLUXNET sites were within the ranges defined by the threshold in contrast, was far from reaching SoS pot at the hemispheric scale TA B L E 1 Parameters of the climate-­constraint model, growing-­season index, and the photo-­threshold model Climate-­constraint model Growing-­season index Photo-­threshold model Tthresh Rthresh C Tmin Tmax DLmin DLmax VPDmin VPDmax DL start DLend F ALL 2.56 26.52 342.2 −1.57 4.33 13.1 13.7 0.88 4.18 12.6 22.1 73.9 ENF 2.07 34.83 112.4 −1.84 4.75 12.0 12.6 0.71 3.85 11.4 21.5 35.7 DBF 2.71 25.90 338.1 −1.99 4.92 10.2 11.2 0.91 4.30 12.6 22.1 258.4 MX 2.50 27.00 350.0 −1.49 4.92 12.1 13.0 0.67 3.75 11.3 14.6 118.8 GRA 2.56 27.30 367.8 −2.01 5.23 10.2 11.4 0.83 4.14 12.6 19.0 95.0 Note: The parameters were calibrated using an optimization algorithm that found the lowest RMSE between modeled and estimated phenological dates at the FLUXNET sites. | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License DESCALS et al. 725 F I G U R E 2 Comparison between observed and simulated start of the growing season (SoS) (in red) and end of the growing season (EoS) (in blue) at the FLUXNET towers in evergreen needleleaved forests (ENF), deciduous broadleaved forests (DBF), mixed forests (MX), and grasslands (GRA). Observed SoS and EoS were obtained from GPP time series from the FLUXNET towers. Simulated SoS and EoS were obtained from ERA5-­L and temperature and radiation time series using the climate-­ constrained model. The titles of the charts report the mean error (ME), the root-­mean-­squared-­error (RMSE), and the number of observations (n). TA B L E 2 Root-­mean-­squared error (RMSE) between modeled Zhang, Parazoo, et al., 2020) and, thus, determine the phenological and estimated phenological dates at the FLUXNET sites dates across latitudes. Our findings highlight the necessity of con- sidering radiation, temperature, and their covariance in modelling RMSE (start of season) RMSE (end of season) the photosynthetically active season in evergreen needleleaved Climate-­ Growing-­ Photo-­ Climate-­ Growing-­ forests and their responses to climatic warming. The results did constraint season threshold constraint season not prove that limiting factors determine the start and end of the model index model model index growing season in other vegetation types other than evergreen ALL 27.22 21.12 22.35 21.72 50.93 needleaved forests. The climate-­constraint model is not suited ENF 11.78 10.74 16.04 12.97 31.11 for modelling the spring onset in deciduous forests, which require DBF 34.38 9.16 58.69 16.71 28.54 heat accumulation (Richardson et al., 2013) and, in certain decidu- MX 21.74 17.87 23.26 13.97 43.12 ous species, chilling requirements (Hänninen, 2016) for bud devel- GRA 33.26 18.70 48.44 34.33 45.85 opment and leaf unfolding. Our results highlight that temperature and radiation constraints (mean difference SoS pot –­SoS was −79.3 d for FLUXNET). SoS was differ spatially during the SoS and EoS, indicating that vegetation only close to SoS pot in Europe. will respond to climate change differently depending on the region of the Northern Hemisphere. Radiation is already constraining the SoS and EoS in regions such as Europe, while the EoS is largely con- 4 | DISCUSSION strained by radiation in the Northern Hemisphere, as indicated pre- viously (Zhang, Commane, et al., 2020), but not in highly elevated Our results indicate that the photosynthetically growing season areas where temperature mainly limit photosynthesis in autumn. begins shortly after conditions for growth become favorable and Moreover, the advance of spring onset was overall not limited by ceases when these conditions decline in evergreen needleaved radiation, and future warming can potentially lead to an advance of forests. Climate constraint boundaries, which define the thresh- the spring onset of carbon uptake in evergreen needleaved forests, old under which photosynthesis is inhibited, can differentiate be- which respond promptly to favorable temperature and radiation tween active and restricted vegetation stages (Jolly et al., 2005; conditions. | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 726 DESCALS et al. F I G U R E 3 Potential lengthening of the growing season in temperate and cold regions of the Northern Hemisphere. (a) Magnitudes of the constraints on the start and end of season (SoS and EoS) by radiation and temperature. The temperature constraints increase with θ, and the radiation constraint decreases with θ. (b) Potential advance of the SoS (SoS − SoSpot) and potential delay of the EoS (EoS − EoSpot). SoSpot and EoSpot represent the start and end of the growing season if radiation was the only limiting factor. SoS and EoS were estimated using the OCO-­2 SIF data set and temperature and radiation were extracted from the ERA5-­L and for 2001–­2020. The climate-­constraint model is similar to the growing-­s eason index, which considers three separate cut-­off functions for tem- perature, daylength, and VPD (Jolly et al., 2005). A difference be- tween these models is that the climate-­constraint model considers the interaction between the constraints of radiation and tempera- ture, which might explain the better estimates of phenological dates. Our model used radiation instead of daylength. Daylength has been proposed as a factor controlling the spring onset (Körner & Basler, 2010) and leaf senescence (Way & Montgomery, 2015). Daylength also plays a role in regulating the forcing requirements in deciduous forests (Meng et al., 2021). However, radiation is more suited for modelling the start and end of the photosyn- thetically active season (Zhang, Commane, et al., 2020), primar- ily because photosynthesis is directly linked to light availability (Herrmann et al., 2020). The cold-­avoidance strategy of deciduous forests may explain their exception to the law of the limiting factors in the phenology of the carbon uptake, as revelead by the long difference between modelled SoS and flux-­based SoS in deciduous forests. Deciduous trees respond differently than evergreen conifers to cold winters. Most conifers in temperate and cold regions have a cold-­tolerant F I G U R E 4 Histograms of the coefficients (R) for the correlation strategy. Leaves in conifers resist cold, reject excessive radiation between temperature (ERA5-­L and) and phenological dates when temperatures are still low in the spring, and resume photosyn- (estimated from the OCO-­2 SIF time series) during 2001–­2020. The thetic activity as soon as climatic conditions are favorable (Chang asterisk represents that mean R was significantly different in the et al., 2021). In contrast, deciduous trees shed their leaves in autumn regions with the highest radiation constrain during the SoS (θ < 65) and leaf budbreak in spring (Chuine et al., 2013). This would account and in the regions with the highest temperature constrain during for our findings that deciduous trees remain dormant during spring, the EoS (θ > 50). The threshold of θ for the SoS represents the 5th percentile of all mean θ values and the threshold of θ for the EoS is even when photosynthesis conditions are favorable. Such strategy the 95th percentile of the mean θ values for the 2001–­2020 period. would also avoid leafing out prematurely to prevent frost damage | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License DESCALS et al. 727 sites provide key examples; daily mean air temperatures remained <2.5°C just before the summer solstice, when annual radiation was at its peak. Because spring is colder than autumn, the growing sea- son begins closer to the summer solstice than the senescence stage does. In contrast, the small difference betwee EoS and EoSpot in- dicates a strong constraint of radiation on GPP senescence in the Northern Hemisphere and suggests that future climatic warming cannot delay the end of the photosynthetically growing season due to low levels of radiation in autumn. During the EoS, the radiation constraint on photosynthetic activity is increased because autumn temperatures remain high. This suggests that, as a result of future warming, radiation constraints will increase further in autumn, limit- ing the delay of the EoS. The strongest radiation constraint was found in temperate forest biomes in Europe, where the maritime influence kept temperatures favorable for vegetation growth in autumn. In contrast, temperature mostly constrained carbon uptake phenology in the alpine biome, which suggests that warming could potentially advance the SoS and delay the EoS in high-­elevated regions more than radiation-­ restricted vegetation. This further support the claims of spatial ho- mogeneity in spring phenology (Liu et al., 2019) and the reduction of the altitudinal gradient in phenological dates as a result of global warming (Vitasse et al., 2018). As found in previous studies, high-­ elevated vegetation will lengthen the growing season more than low-­elevated vegetation with climatic warming, equalizing pheno- logical dates across elevations. Phenology research has had a major focus on the timing of phenophases in temperate deciduous forests (Piao et al., 2019). The PEP725, the most widely used in situ phenological database, mostly contains phenophase records of deciduous species in Central Europe (Templ et al., 2018). However, we showed that the climate of Europe has milder temperatures than other areas of the same lati- tude, exacerbating the limitation of radiation in both spring and au- tumn. Secondly, deciduous forests have high heat requirements for their release from dormancy that differ from other vegetation types F I G U R E 5 Potential lengthening of the growing season in (Richardson et al., 2013), which explains their different responses the FLUXNET sites. Y-­axis show the latitudinal gradients of the to favorable weather conditions during the spring onset. Lastly, observed start of season (SoS; red dots) and the End of Season the coverage of deciduous trees and mixed forests represents only (EoS; blue dots) obtained from the GPP time series at the FLUXNET towers. The mean potential start of season SoSpot and end of 18% of the temperate and cold regions of the Northern Hemisphere season EoSpot (black dots) represent the Day of Year (DoY) when (Figure S1). All these factors indicate that forecasted shifts in phe- shortwave radiation surpasses the 27 W m−2 threshold. The DoY nology under climatic warming and projections of carbon sequestra- was offset such that the summer solstice (DoYsolstice) is the origin of tion by vegetation, based on the sole modelling of phenophases of the axis. deciduous trees in Central Europe, cannot be scaled to other biomes of the Northern Hemisphere. Our findings highlight the importance (Meng et al., 2021). The ecodormancy stage is only broken when a of extensive carbon uptake observations for the study of the global certain amount of heat is accumulated (Richardson et al., 2013). carbon cycle. Carbon uptake phenology is more difficult to track Temperature was the main limiting factor during spring onset, than tree phenophases, and only in situ (FLUXNET) and satellite (SIF) but temperature and radiation have been limiting during autumn in observations have been made extensively over the last two decades temperate and cold regions of the Northern Hemisphere (Zhang, (Pastorello et al., 2020; Zhang et al., 2018). Commane, et al., 2020). Temperature is seasonally asymmetrical, The asymmetric temperature variability over the season explains with lower temperatures before the summer solstice than after. That the different radiation and temperature constraints during the SoS is because leafless temperate deciduous forests and snow-­covered and EoS. In turn, the different radiation and temperature constraints ecosystems have a high albedo in the springtime. Tundra FLUXNET explain the nonuniform seasonal temperature dependency during | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 728 DESCALS et al. the SoS and EoS and across the temperate and cold regions, suggest- C O N FL I C T O F I N T E R E S T ing divergent responses of vegetation phenology to future climatic The authors declare that there is no conflict of interest. warming. The dependency of the SoS and EoS to temperature was significantly lower in regions where radiation constraints were the DATA AVA I L A B I L I T Y S TAT E M E N T highest. These results indicate that the sensitivity of vegetation phe- The data that support the findings of the study are openly available from nology to temperature decreases with the increasing constraints of the FLUXNET2015 dataset at https://fluxn​et.org/data/fluxn​et201​ radiation. Our results show that the senescence stage has a low tem- 5-­datas​et/, the contiguous OCO-­2 SIF at https://figsh​are.com/artic​ perature dependency due to the constraints of radiation, and this les/datas​et/CSIF/6387494, the ERA5-­Land Hourly -­ ECMWF Climate might be a reason for the lower magnitude in the EoS delay than the Reanalysis at https://doi.org/10.24381/​cds.e2161bac, the MODIS SoS advance (Menzel et al., 2006). In deciduous forests, other alter- Land Cover Dynamics MCD12Q2v6 at https://doi.org/10.5067/ native mechanisms for leaf senescence include the growing-­season MODIS/​MCD12​Q2.006, and the RESOLVE Ecoregions 2017 at productivity (Zani et al., 2020). SoS could continue its advance with https://doi.org/10.1093/biosc​i/bix014 and in Google Earth Engine at future warming, even if it is at a lower pace (Fu et al., 2015), although https://devel​opers.google.com/earth​-­engin​e/datas​ets/catal​og/RESOL​ the advance of spring onset would be highly spatially variable de- VE_ECORE​GIONS_2017. Furthermore, a sample python code of the pending on the local radiation limitations. In needleaved forests, climate-­constraint model is available at https://github.com/adria​desca​ radiation should be taken into account when modeling how the phe- ls/clima​te_const​raint_model. The code for the Maximum Separation nology of the photosynthetically active season will change as the method, which was used for the phenology extraction, can be found at climate warms. https://github.com/adria​desca​ls/Maxim​umSep​arati​on_method. ORCID 5 | CO N C LU S I O N S Adrià Descals https://orcid.org/0000-0003-1644-3036 Aleixandre Verger https://orcid.org/0000-0001-9374-1745 This study showed that limiting factors determined the start and end Gaofei Yin https://orcid.org/0000-0002-9828-7139 of the photosynthetically active season in evergreen needleaved Iolanda Filella https://orcid.org/0000-0001-6262-5733 forests in temperate and cold regions of the Northern Hemisphere. Yongshuo H. Fu https://orcid.org/0000-0002-9761-5292 Our findings showed that temperature and radiation constraints on Shilong Piao https://orcid.org/0000-0001-8057-2292 vegetation productivity differed depending on the regions of the Ivan A. Janssens https://orcid.org/0000-0002-5705-1787 Northern Hemisphere. In addition, the radiation constraints were Josep Peñuelas https://orcid.org/0000-0002-7215-0150 different during the SoS and EoS, indicating that carbon uptake phenology will respond to future warming differently in spring than REFERENCES autumn. Temperature was overall limiting photosynthesis during the Chang, C. Y.-­Y., Bräutigam, K., Hüner, N. P., & Ensminger, I. (2021). spring onset, but temperature and radiation were both constraining Champions of winter survival: Cold acclimation and molecular reg- ulation of cold hardiness in evergreen conifers. New Phytologist, at the end of the growing season. The radiation constraint increases 229(2), 675–­691. in autumn because temperatures remain high at that stage of the Chuine, I. (2010). Why does phenology drive species distribution? season, suggesting that future warming will not delay the end of Philosophical Transactions of the Royal Society B: Biological Sciences, the photosynthetically active season. In contrast, radiation limita- 365(1555), 3149–­3160. Chuine, I., de Cortazar-­Atauri, I. G., Kramer, K., & Hänninen, H. (2013). tion was low during the spring onset and the start of the season Plant development models. In Phenology: An integrative environmen- was highly sensitive to temperature, which indicates that the photo- tal science (pp. 275–­293). Springer. synthetically active season could continue overall its advance with Descals, A., Verger, A., Yin, G., & Peñuelas, J. (2020). A threshold method future warming. Our results highlight the importance of taking radia- for robust and fast estimation of land-­surface phenology using Google Earth Engine. IEEE Journal of Selected Topics in Applied Earth tion, temperature, and their covariance into account when modeling Observations and Remote Sensing, 14, 601–­606. the photosynthetically active season in evergreen needleleaved for- Dinerstein, E., Olson, D., Joshi, A., Vynne, C., Burgess, N. D., Wikramanayake, ests and its responses to climatic warming. E., Hahn, N., Palminteri, S., Hedao, P., Noss, R., Hansen, M., Locke, H., Ellis, E. C., Jones, B., Barber, C. V., Hayes, R., Kormos, C., Martin, V., AC K N OW L E D G M E N T S Crist, E., … Saleem, M. (2017). An ecoregion-­based approach to pro- tecting half the terrestrial realm. Bioscience, 67(6), 534–­545. The authors acknowledge the financial support from the European Eamus, D., & Prior, L. (2001). Ecophysiology of trees of seasonally dry Research Council Synergy project SyG-­2013-­610028 IMBALANCE-­P, tropics: Comparisons among phenologies. In Advances in Ecological the Catalan Government grants AGAUR-­2020PANDE00117 and SGR Research (pp. 113–­197). Elsevier. 2017-­1005, the Spanish Government grant PID2019-­110521GB-­I 00, Forsythe, W. C., Rykiel, E. J., Jr., Stahl, R. S., Wu, H., & Schoolfield, R. M. (1995). A model comparison for daylength as a function of latitude the Fundación Ramon Areces grant ELEMENTAL-­CLIMATE, and the and day of year. Ecological Modelling, 80(1), 87–­95. National Funds of China for Distinguished Young Scholar grant No. Friedl, M., & Sulla-­Menashe, D. (2015). MCD12Q1 MODIS/Terra+ aqua 42025101. This work represents a contribution to CSIC Thematic land cover type yearly L3 global 500m SIN grid V006 [data set]. Interdisciplinary Platform Teledetección (PTI-­TELEDETECT). NASA EOSDIS Land Processes DAAC, 10, 200. | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License DESCALS et al. 729 Fu, Y., Li, X., Zhou, X., Geng, X., Guo, Y., & Zhang, Y. (2020). Progress Muñoz-­Sabater, J., Dutra, E., Agustí-­Panareda, A., Albergel, C., Arduini, in plant phenology modeling under global climate change. Science G., Balsamo, G., Boussetta, S., Choulga, M., Harrigan, S., Hersbach, China Earth Sciences, 63(9), 1237–­1247. H., Martens, B., Miralles, D. G., Piles, M., Rodríguez-­Fernández, N. Fu, Y. H., Zhao, H., Piao, S., Peaucelle, M., Peng, S., Zhou, G., Ciais, P., J., Zsoter, E., Buontempo, C., & Thépaut, J.-­N. (2021). ERA5-­L and: Huang, M., Menzel, A., Peñuelas, J., Song, Y., Vitasse, Y., Zeng, Z., A state-­of-­the-­art global reanalysis dataset for land applications. & Janssens, I. A. (2015). Declining global warming effects on the Earth System Science Data Discussions, 13(9), 4349–­4383. phenology of spring leaf unfolding. Nature, 526(7571), 104–­107. Myneni, R. B., Keeling, C., Tucker, C. J., Asrar, G., & Nemani, R. R. (1997). https://doi.org/10.1038/natur​e15402 Increased plant growth in the northern high latitudes from 1981 to Hänninen, H. (2016). The annual phenological cycle. In Boreal and tem- 1991. Nature, 386(6626), 698–­702. perate trees in a changing climate (pp. 35–­138). Springer. Pastorello, G., Trotta, C., Canfora, E., Chu, H., Christianson, D., Cheah, Herrmann, H. A., Schwartz, J.-­M., & Johnson, G. N. (2020). From empiri- Y.-­W., Poindexter, C., Chen, J., Elbashandy, A., Humphrey, M., Isaac, cal to theoretical models of light response curves-­linking photosyn- P., Polidori, D., Reichstein, M., Ribeca, A., van Ingen, C., Vuichard, thetic and metabolic acclimation. Photosynthesis Research, 145(1), N., Zhang, L., Amiro, B., Ammann, C., … Papale, D. (2020). The 5–­14. FLUXNET2015 dataset and the ONEFlux processing pipeline for Jeong, S.-­J., Schimel, D., Frankenberg, C., Drewry, D. T., Fisher, J. B., eddy covariance data. Scientific Data, 7(1), 1–­27. Verma, M., Berry, J. A., Lee, J.-­E., & Joiner, J. (2017). Application of Peaucelle, M., Janssens, I. A., Stocker, B. D., Ferrando, A. D., Fu, Y. H., satellite solar-­induced chlorophyll fluorescence to understanding Molowny-­Horas, R., Ciais, P., & Peñuelas, J. (2019). Spatial variance large-­scale variations in vegetation phenology and function over of spring phenology in temperate deciduous forests is constrained by northern high latitude forests. Remote Sensing of Environment, 190, background climatic conditions. Nature Communications, 10(1), 1–­10. 178–­187. Peñuelas, J., & Filella, I. (2001). Responses to a warming world. Science, Jiao, W., Wang, L., Smith, W. K., Chang, Q., Wang, H., & D'Odorico, P. 294(5543), 793–­795. (2021). Observed increasing water constraint on vegetation growth Peñuelas, J., & Filella, I. (2009). Phenology feedbacks on climate change. over the last three decades. Nature Communications, 12(1), 1–­9. Science, 324(5929), 887–­888. Johnson, A., Alberts, B., Lewis, J., Raff, M., Roberts, K., & Walter, P. Piao, S., Liu, Q., Chen, A., Janssens, I. A., Fu, Y., Dai, J., Liu, L., Lian, X., (2008). Molecular biology of the cell, garland science. Taylor & Francis Shen, M., & Zhu, X. (2019). Plant phenology and global climate Group. change: Current progresses and challenges. Global Change Biology, Jolly, W. M., Nemani, R., & Running, S. W. (2005). A generalized, bio- 25(6), 1922–­1940. https://doi.org/10.1111/gcb.14619 climatic index to predict foliar phenology in response to climate. Richardson, A. D., Andy Black, T., Ciais, P., Delbart, N., Friedl, M. A., Global Change Biology, 11(4), 619–­632. Gobron, N., Hollinger, D. Y., Kutsch, W. L., Longdoz, B., Luyssaert, S., Julien, Y., & Sobrino, J. (2009). Global land surface phenology trends Migliavacca, M., Montagnani, L., Munger, J. W., Moors, E., Piao, S., from GIMMS database. International Journal of Remote Sensing, Rebmann, C., Reichstein, M., Saigusa, N., Tomelleri, E., … Varlagin, 30(13), Art. 13. A. (2010). Influence of spring and autumn phenological transitions Keenan, T. F., Gray, J., Friedl, M. A., Toomey, M., Bohrer, G., Hollinger, D. on forest ecosystem productivity. Philosophical Transactions of the Y., Munger, J. W., O'Keefe, J., Schmid, H. P., Wing, I. S., Yang, B., & Royal Society B: Biological Sciences, 365(1555), 3227–­3246. Richardson, A. D. (2014). Net carbon uptake has increased through Richardson, A. D., Keenan, T. F., Migliavacca, M., Ryu, Y., Sonnentag, O., warming-­induced changes in temperate forest phenology. Nature & Toomey, M. (2013). Climate change, phenology, and phenological Climate Change, 4(7), 598–­604. control of vegetation feedbacks to the climate system. Agricultural Keenan, T. F., Richardson, A. D., & Hufkens, K. (2020). On quantifying and Forest Meteorology, 169, 156–­173. https://doi.org/10.1016/j. the apparent temperature sensitivity of plant phenology. New agrfo​rmet.2012.09.012 Phytologist, 225(2), 1033–­1040. Sun, Y., Frankenberg, C., Wood, J. D., Schimel, D., Jung, M., Guanter, Körner, C., & Basler, D. (2010). Phenology under global warming. Science, L., Drewry, D., Verma, M., Porcar-­C astell, A., Griffis, T. J., Gu, L., 327(5972), 1461–­1462. Magney, T. S., Köhler, P., Evans, B., & Yuen, K. (2017). OCO-­2 ad- Lagarias, J. C., Reeds, J. A., Wright, M. H., & Wright, P. E. (1998). vances photosynthesis observation from space via solar-­induced Convergence properties of the Nelder–­Mead simplex method in chlorophyll fluorescence. Science, 358(6360), eaam5747. low dimensions. SIAM Journal on Optimization, 9(1), 112–­147. Templ, B., Koch, E., Bolmgren, K., Ungersböck, M., Paul, A., Scheifinger, Le Quéré, C., Raupach, M. R., Canadell, J. G., Marland, G., Bopp, L., Ciais, H., Busto, M., Chmielewski, F.-­M., Hájková, L., Hodzić, S., Kaspar, P., Conway, T. J., Doney, S. C., Feely, R. A., Foster, P., Friedlingstein, F., Pietragalla, B., Romero-­Fresneda, R., Tolvanen, A., Vučetič, V., P., Gurney, K., Houghton, R. A., House, J. I., Huntingford, C., Zimmermann, K., & Zust, A. (2018). Pan European Phenological Levy, P. E., Lomas, M. R., Majkut, J., Metzl, N., … Woodward, F. I. database (PEP725): A single point of access for European data. (2009). Trends in the sources and sinks of carbon dioxide. Nature International Journal of Biometeorology, 62(6), 1109–­1113. Geoscience, 2(12), 831–­836. Vitasse, Y., Signarbieux, C., & Fu, Y. H. (2018). Global warming leads to Liu, Q., Piao, S., Fu, Y. H., Gao, M., Peñuelas, J., & Janssens, I. A. (2019). more uniform spring phenology across elevations. Proceedings of Climatic warming increases spatial synchrony in spring vegetation the National Academy of Sciences of the United States of America, phenology across the Northern Hemisphere. Geophysical Research 115(5), 1004–­1008. Letters, 46(3), 1641–­1650. Way, D. A., & Montgomery, R. A. (2015). Photoperiod constraints on tree Meng, L., Zhou, Y., Gu, L., Richardson, A. D., Peñuelas, J., Fu, Y., Wang, phenology, performance and migration in a warming world. Plant, Y., Asrar, G. R., De Boeck, H. J., Mao, J., Zhang, Y., & Wang, Z. Cell and Environment, 38(9), 1725–­1736. (2021). Photoperiod decelerates the advance of spring phenology Yin, G., Verger, A., Filella, I., Descals, A., & Peñuelas, J. (2020). Divergent of six deciduous tree species under climate warming. Global Change estimates of forest photosynthetic phenology using structural Biology, 27(12), 2914–­2927. and physiological vegetation indices. Geophysical Research Letters, Menzel, A., Sparks, T. H., Estrella, N., Koch, E., Aasa, A., Ahas, R., Alm-­ 47(18), e2020GL089167. Kübler, K., Bissolli, P., Braslavská, O., Briede, A., Chmielewski, F. M., Yuan, W., Zheng, Y., Piao, S., Ciais, P., Lombardozzi, D., Wang, Y., Ryu, Crepinsek, Z., Curnel, Y., Dahl, Å., Defila, C., Donnelly, A., Filella, Y., Chen, G., Dong, W., Hu, Z., Jain, A. K., Jiang, C., Kato, E., Li, S., Y., Jatczak, K., Måge, F., … Zust, A. (2006). European phenological Lienert, S., Liu, S., Nabel, J. E. M. S., Qin, Z., Quine, T., … Yang, S. response to climate change matches the warming pattern. Global (2019). Increased atmospheric vapor pressure deficit reduces Change Biology, 12(10), 1969–­1976. global vegetation growth. Science Advances, 5(8), eaax1396. | 13652486, 2023, 3, Downloaded from https://onlinelibrary.wiley.com/doi/10.1111/gcb.16502 by Readcube (Labtiva Inc.), Wiley Online Library on [13/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 730 DESCALS et al. Zani, D., Crowther, T. W., Mo, L., Renner, S. S., & Zohner, C. M. (2020). S U P P O R T I N G I N FO R M AT I O N Increased growing-­season productivity drives earlier autumn leaf Additional supporting information can be found online in the senescence in temperate trees. Science, 370(6520), 1066–­1071. Supporting Information section at the end of this article. Zhang, Y., Commane, R., Zhou, S., Williams, A. P., & Gentine, P. (2020). Light limitation regulates the response of autumn terrestrial carbon uptake to warming. Nature Climate Change, 10(8), 739–­743. Zhang, Y., Joiner, J., Alemohammad, S. H., Zhou, S., & Gentine, P. (2018). A global spatially contiguous solar-­induced fluorescence (CSIF) How to cite this article: Descals, A., Verger, A., Yin, G., dataset using neural networks. Biogeosciences, 15(19), 5779–­5800. Filella, I., Fu, Y. H., Piao, S., Janssens, I. A., & Peñuelas, J. Zhang, Y., Parazoo, N. C., Williams, A. P., Zhou, S., & Gentine, P. (2020). (2023). Radiation-­constrained boundaries cause nonuniform Large and projected strengthening moisture limitation on end-­ of-­season photosynthesis. Proceedings of the National Academy of responses of the carbon uptake phenology to climatic Sciences of the United States of America, 117(17), 9216–­9222. warming in the Northern Hemisphere. Global Change Biology, Zohner, C. M., Benito, B. M., Svenning, J.-­C ., & Renner, S. S. (2016). Day 29, 719–730. https://doi.org/10.1111/gcb.16502 length unlikely to constrain climate-­driven shifts in leaf-­out times of northern woody plants. Nature Climate Change, 6(12), 1120–­1123.