Climatic Warming Increases Spatial Synchrony in Spring Vegetation Phenology Across the Northern Hemisphere

Climatic warming has advanced spring phenology across the Northern Hemisphere, but the spatial variability in temperature sensitivity of spring phenology is substantial. Whether spring phenology will continue to advance uniformly at latitudes has not yet been investigated. We used Bayesian model averaging and four spring phenology models, and demonstrated that the start of vegetation growing season across the Northern Hemisphere will become substantially more synchronous (up to 11.3%) under future climatic warming conditions. Larger start of growing season advances are expected at higher than lower latitudes (3.7–10.9 days earlier) due to both larger rate in spring warming at higher latitudes and larger decreases in the temperature sensitivity of start of growing season at low latitudes. The consequent impacts on the northern ecosystems due to this increased synchrony may be considerable and thus worth investigating.


Introduction
The rising spring temperatures in recent decades have led to a notable advance of spring phenology over the Northern Hemisphere (NH; Fu et al., 2016;Menzel et al., 2006;Peñuelas & Filella, 2001;Piao et al., 2015;Schwartz et al., 2006), which has substantially affected terrestrial water, carbon, and nutrient balances (Keenan et al., 2014;Peñuelas & Filella, 2009;Piao et al., 2007Piao et al., , 2017Richardson et al., 2010). A recent study reported a declining response of spring phenology to warming temperature, attributed to winter chilling loss  and questioning whether spring phenology across different regions would continue to advance as climate warms. Extrapolating current knowledge on the bioclimatic gradient in spring phenology (e.g., a progressive delay in spring phenology with increasing latitude, "Hopkin's bioclimatic law"; Hopkins, 1920Hopkins, , 1938 to the future is therefore challenging. Recent studies have reported an increase in synchrony of spring phenology (i.e., reduced differences between the start of growing season (SOS) among species) in recent decades from in situ observations at regional scales (C. Wang et al., 2016) and along altitudinal gradients (Vitasse et al., 2018), attributed to spatial differences in warming speed and phenological responses to rising temperature. Whether such increasing synchrony continues to occur across the entire NH under warmer scenarios, however, has not yet been investigated. We quantified the changes in the spatial variance of spring phenology across the NH under current and future conditions based on model simulations to explore the changes in the synchrony of SOS.
Models of spring phenology have been widely used to reproduce SOS in previous studies (Jeong et al., 2013;Liu et al., 2018;Melaas et al., 2015), but considerable uncertainties nonetheless remain within and among models due to limited understanding of mechanisms underlying spring phenology (Chuine & Régnière, 2017;Richardson et al., 2012). Multiple model predictions, combining the strengths of individual models and reducing the computed uncertainties, are therefore required, although comprehensive understanding of spring phenology still need well-designed experimental studies. Previous multimodel studies simply assigned equal weights to all models, regardless of their performance, and calculated their ensemble mean (i.e., simple model averaging (SMA)). We, however, applied Bayesian model averaging (BMA) to deal with individual model uncertainty by assigning higher weights to models with higher accuracies in the observational period. Four models of spring phenology (see section 2), were applied to simulate SOS over the NH under current (1990NH under current ( -2009 and future (2080-2099; four Representative Concentration Pathways (RCPs)) conditions.

Climatic Data Sets
Daily mean temperature (T) was obtained from the CRU-NCEP (Climatic Research Unit (CRU) and National Centers for Environmental Protection (NCEP)) v5 data set at a spatial resolution of 0.5 × 0.5°for 1982-2012. This data set combines the CRU-TS (time series) v3.1 0.5 × 0.5°monthly climatological data  and the NCEP reanalysis at 2.5 × 2.5°and a 6-hr time interval since 1948 (Mitchell & Jones, 2005;New et al., 2000). For future projections, we acquired data for daily T from predictions by Coupled Model Intercomparison Project phase 5 models (Table S1). The Coupled Model Intercomparison Project phase 5 model predictions for the historical , RCP2.6, RCP4.5, RCP6.0, and RCP8.5 (2006and 2070-2099 scenarios were processed and interpolated to a spatial resolution of 1 × 1°using Climate Data Operators (Schulzweida, 2013).

Satellite-Derived SOS
In this study, we extracted SOS from the latest Normalized Difference Vegetation Index (NDVI) data set by NASA's GIMMS group (NDVI 3g.v1 ) at a spatial resolution of 1/12 × 1/12°spanning the period 1982-2012. Apart from the issues which have already been addressed in NDVI 3g (e.g., update of satellite sensors, atmospheric interference, and nonvegetation dynamics; Pinzon & Tucker, 2014), artifacts due to snow coverage and changes in calibration were further processed in NDVI 3g.v1 . The SOS extraction methods include two main steps: (1) reconstructing NDVI time series to a daily basis using data filter function and (2) determining SOS from predefined thresholds or changing characteristics of the filtered NDVI curve. Following Liu, Fu, Zeng, et al. (2016) and Liu, Fu, Zhu, et al. (2016), we estimated the date of SOS using four commonly used methods (i.e., Hants-Mr (Jakubauskas et al., 2001), Polyfit-Mr (Piao et al., 2006), double logistic (Julien & Sobrino, 2009), and piecewise logistic (Zhang et al., 2003)). To reduce the uncertainty resulted from intermethod differences, we applied SOS averaged from four methods in our analysis.

Models of Spring Phenology
We used four models of spring phenology involving both endodormancy and ecodormancy phases. These models assume that chilling and forcing (subsequent to chilling) are required to break endodormancy and ecodormancy, respectively, and define the day when the accumulation of forcing becomes critical as the date of leaf unfolding (Cannell & Smith, 1983;Murray et al., 1989). The four models include three temperaturedriven models, that is, the Sequential Model (Hänninen, 1990;Kramer, 1994), Unified Model, and UniChill Model (I. Chuine, 2000). Besides, the DORMPHOT model additionally considered the effect of photoperiod and assumed that reduced temperature and photoperiod co-triggered dormancy, and a long day length would promote the accumulation of forcing . These models were calibrated using satellite-derived SOS and daily temperature at plant function type (MacBean et al., 2015) scale for the period 1982-2012 using Bayesian optimization technique (Martinez-Cantin, 2014). The detailed model implementation could be found in Table S2 and Text S1.

BMA
Previous studies have mostly been based on results from either a single model of spring phenology or the average of multiple models. Simple model averaging (SMA) is superior to a single model but does not identify the best model. BMA, however, assigns higher weights to models that have previously performed better and is therefore able to obtain the minimum root-mean-square error between the weighted model ensemble and observations (Raftery & Zheng, 2003). We applied BMA to determine the optimal weight for each model of spring phenology at a pixel level for 1982-2012. Supposing that SOS anomaly was forecasted by on the basis of training data SOS T anomaly (randomly selected as 80% from the entire study period), using models of spring phenology (M 1 , … , M k ), then the forecast probability density function, p(SOS), is written as where p(SOS|M k ) is the normal conditional probability density function based on model M k and p(M k |SOS T ) is the posterior probability of model M k , which is corrected by the training data and also indicates how well model M k fits the training data (i.e., weight of each model). This term can be calculated as where p(M k ) is the prior probability of model M k . We defined the probability of the training data (e.g., SOS T anomaly) as the likelihood function following Raftery et al. (2005), the optimal weights (e.g., p(M k |SOS T )) were determined when the likelihood function reached its maximum. More details about BMA are provided by Raftery et al. (2005) and Vrugt and Robinson (2007).

Sensitivity of Spring Phenology to Temperature
We calculated the sensitivity (S T ) of SOS to mean spring temperature (T spr ) as the ratio between changes in SOS (ΔSOS) and changes in T spr (ΔT spr ) during 2080-2099 and 1990-2009: We calculated BMA-based SOS for each year over these two periods. T spr for individual years was calculated as the average of daily mean temperature from 1 March to 31 May. We assessed the influence of future warming on SOS by calculating S T across the NH under four RCPs. We first calculated ΔSOS and ΔT spr for each model and then applied their ensemble mean to reduce the uncertainties within predictions by the Coupled Model Intercomparison Project phase 5 models (Table S1).

Analysis
We ran the models of spring phenology under five independent scenarios to investigate whether the response of SOS to temperature might be constrained by shorter photoperiod as SOS continues to advance or chilling loss due to winter warming. We chose the DORMPHOT model to test the effect of photoperiod.
We first ran this model with the actual photoperiod (S1), then increased the photoperiod during spring (March-May) by 1 and 3 hr, individually (S2). The remaining simulations were carried over all four models. We calculated winter temperature from 1 September of the previous year to 31 January of the current year, because the accumulation of chilling begins on 1 September of the previous year in these models. In the third scenario (S3), we randomly selected winter temperature from 1982 to 2009 and kept it constant during 2080-2099. Similarly, we considered varying winter temperature but kept spring temperature (March-May) constant during 2080-2099 (S4). To reduce the sampling bias, we additionally ran S3 and S4 for 9 times and applied the average of 10 simulations in our analysis. In the fifth scenario (S5), we considered varying winter and spring temperature for 2080-2099. We thereafter calculated BMA-based SOS projections for 2080-2099 and then the differences among these scenarios. We conducted this analysis based on temperature projections from four RCPs to test the influence of longer photoperiod (S2-S1), winter warming (S3-S5), and spring warming (S4-S5) on the date of SOS.

BMA
The weight assigned by BMA served as a good criterion for model evaluation for 1982-2012. The Unified Model performed best in~30% of the NH, that is, northeastern Siberia and western North America. The UniChill and Sequential Models produced the best agreement with satellite-derived SOS in~25 and 24% of the NH, respectively. The UniChill Model outperformed the other models in central Eurasia, and the Sequential Model performed the best in northeastern Europe and northern Siberia. The DORMPHOT model agreed best in~21% of the NH, albeit fragmentally distributed (Figure 1a). Plants differ across the regions, and species differ in their response to chilling, photoperiod, and spring warming, so different models are likely to be better for different regions. Benefiting from the strength of individual models, the BMA-based SOS simulations were in better agreement with satellite observations across the NH. As expected, rootmean-square errors were consistently lower compared to simulations based on either SMA or the single models, particularly over southwestern North America and central Eurasia (Figures 1b-1d and S1). The median root-mean-square error decreased substantially from 9.3-10.2 days for the individual models to 8.8 days for SMA and 5.8 days for BMA (Figure 1d). These differences were less apparent in the temporal correlations between the simulated and satellite-derived SOSs (Figures 1e-1g and S2).

Projected Changes in Spring Phenology and Its Synchrony
SOS across the NH was obtained for 2080-2099 by extrapolating the same weights of the phenological models to future simulations. Regionally different responses to the projected climate change were obtained (Figures 2 and S4; based on ensemble mean of model projections under each RCP, because the intermodel differences showed similar magnitude) ( Figure S3). SOS was generally significantly earlier in most of the NH (ranging from 75% in RCP2.6 to 91% in RCP8.5, p < 0.05) between 2080-2099 and 1990-2009. SOS, however, was sometimes also delayed, especially under warmer scenarios in the low latitudes (e.g., 30-50°N; Figures 2a, 2b, 3a, and 3b). SOS across the NH was 7.9 ± 6.8 days (RCP2.6) to 19.1 ± 12.3 days (RCP8.5) earlier (average ± SD) during 2080-2099 than the current period (1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009). SOS under RCP2.6 scenario advanced more (e.g., >12 days) at high latitudes (e.g., northwestern North America, northern Europe, northern Siberia) and the Tibetan Plateau, followed by low latitudes (e.g., southeastern Siberia and central Eurasia; Figure 2a). This spatial pattern of faster SOS advances at higher than lower latitudes was more apparent under warmer RCPs. For example, the difference in the average SOS advancement between higher and low latitudes increased from 4 days under RCP2.6 to 11 days under RCP8.5 (Figures 2b, S4a, S4b, 3a, 3b, S5a, and S5b). The weakening latitudinal gradient of SOS under future climatic warming indicated the increasing synchrony of spring phenology across the NH. The spatial standard deviation of SOS decreased by 4.0% from 32.7 days (~99% of 33.0 days based on satellite observations) for 1990-2009 to 31.4 days for 2090-2099 (RCP2.6) and continued to decrease under warmer RCPs, that is, 30.8 days (5.8%, RCP4.5), 30.4 days (7.0%, RCP6.0), and 29.0 days (11.3%, RCP8.5; Figure 4c). We found that the decrease in the spatial variation of SOS was related to the regional differences in the response of vegetation to the future climatic warming. Average SOS for 1990-2009 was significantly negatively correlated (p < 0.05) with its projected future change (Figures 4a, 4b, S6a, and S6b), implying that regions with later SOSs at current condition (usually at high latitudes) would likely have larger SOS advances under future climatic warming. These findings were also confirmed by each of the four spring phenology models ( Figure S7) and thus indicate that the vegetation across the NH will initiate spring growth more synchronously in the future, especially between the low and high latitudes.

Determinations of the Increasing Synchrony of Spring Phenology
The increasing synchrony of spring phenology between 2080-2099 and 1990-2009 is co-determined by spring temperature (T spr ) and the sensitivity of SOS to spring warming (S T ), which both differ between lower and higher latitudes. S T was higher at low latitudes (e.g., 30-50°N) than at high latitudes (e.g., >60°N,  2012), implying that vegetation at lower latitudes (or with earlier SOS) is more sensitive to changes in T spr .
The changes in SOS between these periods nonetheless did not always match the spatial pattern of S T , because of the much larger increase in T spr at higher latitudes (e.g., 0.2-0.7°C higher than at low latitudes; Figures 3c, 3d, S5c, and S5d). For example, S T was lower for northern Siberia than the Tibetan Plateau, but the advance in SOS was similar between these regions due to the faster increase in T spr across northern Siberia (e.g., RCP2.6; Figures 2a and 2e). This finding was consistently confirmed under the remaining RCPs (Figures 2 and S4). Although the S T was lower at higher latitudes, the faster increasing in T spr could also lead to larger advancement in SOS, thereby reducing the spatial differences in predicted SOS across the NH (e.g., RCP2.6; Figure 2). Moreover, S T decreased more at lower than higher latitudes across different RCPs. The decrease in S T between the results under RCP2.6 and RCP8.5 was largest at lower latitudes (30-50°N), averaging 4.1 days/°C (Figures 3e and 3f). The average decrease in S T at higher latitudes (>60°N), however, was only 2.7 days/°C (Figures 3e and 3f). The comparison between S T estimated from RCP2.6 and the two intermediate RCPs confirmed this finding (Figures 3e,  S5e, and S5f). Greater synchrony in spring phenology is thus projected across the NH (Figure 4c), accompanied with the faster spring warming at the high latitudes (Figures 3c, 3d, S5c, and S5d).
Several hypotheses have been proposed to explain the faster decrease in S T at low latitudes. A first hypothesis correlates the reduced S T with the variance of local spring temperature (T spr ). A recent ground-based study found that plants at locations with larger T spr variances were less sensitive to warming (T. Wang et al., 2014), because plants under such conditions tend to depend on more stable environmental cues, such as photoperiod (Körner & Basler, 2010) and/or higher chilling requirement (Schwartz & Hanes, 2010), to avoid premature exposure to damaging spring frost. Such reduction in S T due to higher T spr variances, however, was not directly implemented in the temperature-based models of phenology and could not be tested.
The hypothesized role for a negative impact of a shorter photoperiod on S T , however, could be tested across the NH using the DORMPHOT model (Table S2), which involves the effect of photoperiod . We compared SOS projections over 2080-2099 with actual (S1) and increased (S2) photoperiod (see section 2). The difference between S2 and S1 indicated that increasing photoperiod alone would exacerbate the advance of SOS across most of the NH (>98% (Figures S8a-S8d) and >99% ( Figures S8e-S8h), across all RCPs), especially at low latitudes (e.g., southern North America, the Tibetan Plateau) and in northern Europe. The effects of increased photoperiod on advancing SOS decreased significantly toward higher latitudes ( Figure S8). This finding not only confirmed the photoperiodic control of SOS (Caffarra, Donnelly, Chuine, & Jones, 2011;Körner & Basler, 2010;Zohner et al., 2016) but also suggested that the shorter days due to earlier SOS might restrict the advance of SOS in the future, particularly at lower latitudes, thus shaping the spatial pattern of a faster decrease in S T at low latitudes. Photoperiod does not exclusively account for our finding; the other models considering chilling and forcing effects also produced a similar spatial pattern of faster decrease in S T at lower latitudes ( Figure S9).
The control of SOS by cold temperatures during dormancy is a third hypothesis that could account for the regional difference in the future decrease in S T . We tested this hypothesis by applying four chilling-based models for 2080-2099 but kept the winter temperature (S3) or the spring temperature (S4) constant and then assessing their difference to simulations based on varying winter and spring temperatures (S5; Figures S10a-S10h; see section 2). The warming winter led to a delayed SOS, albeit primarily at lower latitudes (except for the Tibetan Plateau; Figures S10a-S10d). The predicted SOSs across most of the NH (mainly at higher latitudes) did not differ significantly between these two simulations (p > 0.05). These findings imply that increasing winter temperatures may indeed lead to insufficient chilling at lower latitudes (except for the Tibetan Plateau), which would subsequently minimize the response of SOS to warming temperatures because of the increasing requirement for forcing temperatures with decreased chilling Laube et al., 2014). In contrast, vegetation at high latitudes would still experience sufficient chilling due to the colder and longer winter, and the increasing temperature might even increase its exposure to chilling in regions where temperature was typically below the base temperature for chilling accumulation (Vandvik et al., 2018). Experimental studies comparing phenology in regions at lower and higher elevations have reported similar findings (Gusewell et al., 2017;Vitasse et al., 2018). Large-scale manipulative field or virtual numerical experiments are still required to further test these hypotheses, because the scenarios designed to explain the faster decrease in S T were less realistic. It should also be noted that our analysis assumes that model bias in the historical period is representative of model performance in the future, and model weights determined under current conditions are preserved for the future. This assumption needs to be tested a wide range of plant species and climatic conditions, because the underlying mechanisms for spring phenology are still not entirely elucidated. For example, photoperiod might influence the spring phenology of some species more strongly under a warmer climate (Körner & Basler, 2010;Way & Montgomery, 2015). In this study, we applied four empirical models, involving chilling, forcing, and photoperiod effects, to account for changes in SOS and observed an increased synchrony in spring phenology. Nonetheless, more efforts should also be placed on experimental studies focused on revealing mechanisms behind the phenological processes (e.g., temperature, photoperiod, and precipitation controls) and improving the predictability of spring phenology models (e.g., Figures 1e and 1f).

Conclusions
We report overall advancing, but heterogeneously distributed, trends in spring phenology across the NH between 1990-2009 and 2080-2099 (under four climatic scenarios). The synchrony of spring phenology is expected to increase (i.e., a decrease in spatial variation) across the NH in the future, because projected advances in SOS for 2080-2099 are larger for regions with later SOSs during 1990-2009 (e.g., higher latitudes). This phenomenon is due mainly to the larger increase in spring temperature at high latitudes, which even compensates for the influence of the generally lower S T at high versus low latitudes. Besides, S T under the warmer climatic scenarios decreased more at low latitudes than at high latitudes, probably due to the shorter day length and/or loss of winter chilling, thus amplifying the phenological synchrony as climatic warming proceeds. Changes in the synchrony of spring phenology may result in synchronous change of land-surface physical conditions (e.g., canopy conductance, albedo, and surface aerodynamic roughness) over large areas (Peñuelas & Filella, 2009;Richardson et al., 2013), and thereby surface water fluxes, energy balance, and atmospheric circulation during spring (Peñuelas & Filella, 2009;Stéfanon et al., 2012). To date, most phenological modules in dynamic global vegetation models relied on growing degree days (Sitch et al., 2003;Thornton et al., 2002) or carbon allocation (Ostle et al., 2009) and yielded poor predictions in spring phenology (Richardson et al., 2012). Therefore, updating the bioclimatic gradient of spring phenology (Hopkins, 1920(Hopkins, , 1938 across the latitude due to climate change and improving the representation of phenological synchrony in dynamic global vegetation models could enhance our understanding and predictive capacity of how northern ecosystems will respond to a future warmer climate.

Data Availability
All data sets and materials used in this study are freely available online, as described in Table S3.