Spring Phenological Sensitivity to Climate Change in the Northern Hemisphere: Comprehensive Evaluation and Driving Force Analysis

: Plant phenology depends largely on temperature, but temperature alone cannot explain the Northern Hemisphere shifts in the start of the growing season (SOS). The spatio–temporal distribution of SOS sensitivity to climate variability has also changed in recent years. We applied the partial least squares regression (PLSR) method to construct a standardized SOS sensitivity evaluation index and analyzed the combined effects of air temperature (Tem), water balance (Wbi), radiation (Srad), and previous year’s phenology on SOS. The spatial and temporal distributions of SOS sensitivity to Northern Hemisphere climate change from 1982 to 2014 were analyzed using time windows of 33 and 15 years; the dominant biological and environmental drivers were also assessed. The results showed that the combined sensitivity of SOS to climate change (S Com ) is most inﬂuenced by preseason temperature sensitivity. However, because of the asymmetric response of SOS to daytime/night temperature (Tmax/Tmin) and non-negligible moderating of Wbi and Srad on SOS, S Com was more effective in expressing the effect of climate change on SOS than any single climatic factor. Vegetation cover (or type) was the dominant factor inﬂuencing the spatial pattern of SOS sensitivity, followed by spring temperature (Tmin > Tmax), and the weakest was water balance. Forests had the highest S Com absolute values. A signiﬁcant decrease in the sensitivity of some vegetation (22.2%) led to a decreasing trend in sensitivity in the Northern Hemisphere. Although temperature remains the main climatic factor driving temporal changes in S Com , the temperature effects were asymmetric between spring and winter (Tems/Temw). More moisture might mitigate the asymmetric response of S Com to spring/winter warming. Vegetation adaptation has a greater inﬂuence on the temporal variability of SOS sensitivity relative to each climatic factor (Tems, Temw, Wbi, Srad). More moisture might mitigate the asymmetric response of S Com to spring/winter warming. This study provides a basis for vegetation phenology sensitivity assessment and prediction. This study showed that the combined sensitivity of SOS to climate change (S Com ) is most inﬂuenced by pre-season temperature sensitivity. However, because of the asymmetric response of SOS to Tmax/Tmin and the non-negligible moderating effect of Wbi and Srad on SOS, S Com was more effective in expressing the effect of climate change on SOS compared with the single climatic factor. Vegetation cover (or type) was the dominant factor inﬂuencing the spatial pattern of SOS sensitivity, followed by spring temperature (Tmin > Tmax); the weakest was Wbi. Forests had the highest S Com absolute values. In recent years, a signiﬁcant decrease in the sensitivity of some vegetation has led to a decreasing trend of sensitivity in the Northern Hemisphere, with vegetation in wetter/cooler climate zones and shrubs contributing more to the decrease in S Com . Biological factors (veg-etation adaptation) exerted a greater inﬂuence on the temporal variability of SOS sensitivity relative to each climatic factor (Tems, Temw, Wbi, Srad). Although temperature remains the main climatic factor driving temporal changes in S Com , the effects were asymmetric between spring and winter. More moisture might mitigate the asymmetric response of S Com to spring/winter warming. Despite the complex response of sensitivity to climate change, the dominant factors in different climatic regions still have different emphases. Identifying systematic differences in phenological climate sensitivity would facilitate the development of indicators and estimates of vulnerability for conservation and national adaptation programs The performance the current phenological models is unsatis-factory The quantitative of and driving force analysis of SOS in this study provides not only a new method for identifying the combined sensitivity of springtime phenology to climate change but also a basis for the construction of better phenological models to improve phenological prediction performance.


Introduction
Plant phenology, especially the start of the growing season (SOS) (or spring green-up date), is a highly sensitive indicator of the impacts of climate change on the biosphere [1][2][3]. The early onset of vegetation spring phenology has been widely reported as an indicator of the climate warming [1,3,4]. Vegetation phenology fine-tunes itself according to the environment; moreover, it is closely linked to climate variability [5,6]. Studies have shown that earlier SOS leads to higher photosynthetic rates and vegetation growth [7], while also increasing transpiration, resulting in increased drought stress [8]. Phenological changes also regulate vegetation feedback to the climate system by influencing the regional and global carbon budget, water flux, and energy balance [4,9]. Therefore, improving the knowledge on phenology change and its drivers is essential to better understand and model the relationship between ecosystems and the climate system.
Plant phenological sensitivity (phenological phase changes caused by unit changes in climatic conditions) determines the magnitude of phenological changes in response to future climate change and is key to understanding the relationship between vegetation phenology and climate change [3,10]. A number of studies have extensively assessed the impacts of global warming on the SOS over the past decades based on ground and satellite data [10][11][12][13][14][15][16][17][18][19][20][21]. Temperature is considered the dominant controlling factor of SOS [13,14,17,18], and studies on phenological sensitivity have also focused on the sensitivity of SOS to temperature (S Tem ) [22][23][24]. However, the relative importance of climate constraints on surface phenology has shifted since 1982 [13]. Recent studies have found that precipitation may play a key role in determining the spring [10,12]. Moreover, snow cover produced by preseason precipitation has also been shown to influence spring phenology at high latitudes and altitudes [15]. Radiation, a comprehensive measure of sunshine duration and light intensity, shows obvious interannual variability [16]. Moreover, the photoperiod plays an important role in controlling SOS in most Northern Hemisphere forests [17]. Radiation indirectly affects phenology, as it influences preseason temperature and water balance by controlling surface heat fluxes and evapotranspiration. Importantly, biological factors can also modulate the SOS [10,16]. The interaction between phenological events in spring and autumn may change the overall phenological response to climate warming [10]. Although SOS is influenced by multiple meteorological factors, most studies have examined only the temperature or precipitation sensitivity of SOS [22][23][24]; only a few studies have evaluated the combined effects of climate change on SOS [25]. This is partly because preseason meteorological factors, such as temperature, moisture utilization, and radiation, have different units of measurement. As a result, the sensitivity of each meteorological factor based on regression analysis cannot be quantitatively compared on the same scale. In addition, considering the S Tem as an example, the length of the phenologically relevant period (preseason) affected by temperature (including Tem, Tmax, and Tmin) varies in different climatic regions (or biomes), and changes in the temperature cumulative time and statistical characteristics may cause estimation biases in S Tem [22,25]. SOS responds differently to Tmax and Tmin for the same preseason duration in some areas [18]; this difference in response is controlled by a number of environmental factors and vegetation conditions [10]. Thus, the phenological response of vegetation to recent climate warming cannot be explained entirely by temperature [26,27]. Given the complex response of SOS to preseason meteorological elements and biological elements, most previous studies were single factor sensitivity studies based on linear regression analysis (univariate and multiple regression) [24] and correlation analysis (or partial correlation analysis) [28]. Therefore, it is necessary to establish a standardized sensitivity evaluation index considering the influence of biological factors to characterize the sensitivity of SOS to each meteorological factor, which would facilitate a comparative analysis of the SOS sensitivity to each meteorological factor, and to construct a comprehensive SOS sensitivity index for climate change.
The spatio-temporal differences in the sensitivity of SOS to climate change are not yet fully understood. The S Tem exhibits a large spatial variability [26] and has a downward trend with global warming [12,29,30]. Fu et al. [29] reported a declining S Tem in European temperate tree species. The causes for the spatio-temporal variation of S Tem are poorly understood, particularly on the global scale [10,24]. Studies have explored the causes of variation in phenological sensitivity from the aspects of climate, geography, and biology [12,[23][24][25]28,[31][32][33][34][35]. Interacting effects of temperature and precipitation have changed the climate sensitivity of spring vegetation green-up [25]. The SOS is more sensitive to pre-season precipitation in arid areas and more sensitive to preseason temperature in humid areas [12]. With warming-induced leaf-out advancement, a reduced photoperiod earlier in spring may decrease the S Tem [10]. The sensitivity of spring vegetation phenology is related to the rate of warming and phenological sensitivity with increasing local spring temperature variance at the species level from ground observations [23]. Spring phenological temperature sensitivity is higher at low latitudes than at cold or high latitudes [34,35]. The sensitivity of vegetation phenology to climate change varies by biome [32]; vegetation types with earlier SOS mean (mean SOS for 1982-2008) are more sensitive to spring temperature in the Northern Hemisphere [28]. Changes in phenological sensitivity are the result of a combination of multiple factors, but few studies have analyzed the dominant drivers of sensitivity variation. Understanding the large-scale spatio-temporal pattern of phenological sensitivity to climate change and its dominant controls will improve our ability to predict the timing of SOS.
In this study, we used four different methods to extract phenology dates from the normalized difference vegetation index (NDVI) time series of the Global Inventory Monitoring and Modeling System (GIMMS) NDVI3g.v1 dataset. We then applied the partial least squares regression (PLSR) to extract the standardized SOS sensitivity index to preseason climatic factors (air temperature, water balance, and radiation) and constructed a comprehensive SOS sensitivity index (S Com ) for climate change. The purpose of this study was to (1) quantify the spatial and temporal variability of the Northern Hemisphere (latitude ≥ 30 • ) SOS sensitivity to climate change (single factor and synthesis) from 1982 to 2014 in 33-year and 15-year time windows, (2) determine the influence of environmental and biological factors on the spatial distribution of SOS sensitivity, and (3) identify the dominant drivers of the temporal variability of SOS sensitivity.

Materials
The NDVI is widely used in global-scale phenological studies [10,20]. In this study, we used the latest NDVI dataset (GIMMS NDVI3g.v1 https://iridl.ldeo.columbia.edu/ SOURCES/.NASA/.ARC/.ECOCAST/.GIMMS/.NDVI3g/.v1p0/ (Accessed 15 May 2021) generated from NOAA/AVHRR series satellite images by the NASA GIMMS group, which covers the period from January 1982 to December 2015 and has a spatial resolution of 1/12 • and an interval of 15 days [36]. NDVI3g.v1 addresses various issues that were present in earlier datasets addressed in NDVI3g.v1 (e.g., updating satellite sensors, atmospheric disturbances, and non-vegetation dynamics), and NDVI3g.v1 further addresses the artifacts caused by snow cover and calibration changes [37]. In this study, we excluded areas of bare soil/sparse vegetation with an annual average NDVI of less than 0.1 [20].
TerraClimate is a high-resolution (1/24 • , 4 km) global monthly climate and climatic water balance dataset covering the period from 1958 to 2015 [38]. It combines high spatial resolution WorldClim climate data with coarse resolution data from other sources to generate monthly precipitation, maximum and minimum temperatures, wind speed, vapor pressure, and solar radiation datasets. TerraClimate also uses a water balance model to generate monthly surface water balance datasets. Compared with coarser datasets, Terra-Climate shows significant improvement with regard to spatial authenticity and the overall mean absolute error. We used the monthly average maximum temperature (Tmax), average minimum temperature (Tmin), accumulated precipitation (Pre), downward shortwave flux at the surface (Srad), soil moisture (Sm), actual evapotranspiration (Aet), and snow water equivalent (Swe) TerraClimate datasets covering 1981-2015.
The land cover (LC) data product from the European Space Agency's Climate Change Initiative (CCI) specifies different vegetation types. The project delivers consistent global land-cover maps at a spatial resolution of 300 m and an annual temporal resolution from 1992 to 2015. The land cover classification system used in this dataset was found to be compatible with the plant functional types used in climate models [39]. We used the new and improved Köppen-Geiger climate classification maps at 1-km resolution [40] for the present day  to specify the different climate regions by referring to the Intergovernmental Panel on Climate Change (IPCC) default climate regions recalculated by the Joint Research Centre (JRC) of the European Commission [41]. We excluded grid points where phenology could not be extracted and land cover types changed, especially for pixels with two or more growth cycles in one year. Detailed information and spatial distribution of climatic regions and plant types are shown in Table 1 and Figure 1. The spatial resolution of the land cover and climatic regions was resampled to match the resolution of the NDVI data.

Extracting Phenology Using NDVI Data
A number of methods have been developed to extract the SOS and EOS (end of the growing season) from seasonal NDVI variability. Phenology extraction is mainly divided into two steps: (1) reconstructing the NDVI time series, and (2) determining the phenological extraction from the change characteristics of the filtered NDVI curve. We used the double logistic function to fit the NDVI time series in this study [20]. The inflection point detection method or the threshold method can be used to extract phenological data from NDVI time series data. We used two inflection point detection methods (i.e., firstorder derivative [42] and second-order derivative [43] methods) and two fixed threshold methods (i.e., 0.2 [44] and 0.5 [45] fixed thresholds), to obtain a total of four methods to extract vegetation SOS and EOS in the Northern Hemisphere from 1982 to 2014. To reduce the uncertainty caused by the different methods, we used the average SOS and EOS of the four methods. Savizky-Golay filters were applied to the GIMMS NDVI3g.v1 datasets to minimize noise prior to phenological estimation [20].

Partial Least Squares Regression (PLSR)
PLSR is a multivariate statistical regression analysis that combines the advantages of multiple regression, typical correlation, and principal component analysis. It is advantageous in the analysis of multiple variables with multiple correlations, particularly when the number of variables is greater than the number of samples [46]. The standardized model regression coefficient (MC), an output of PLSR, indicates the extent to which the independent variables can explain the weight of the dependent variables. In this study, MC was defined as the normalization sensitivity of the dependent variable to the independent variable. Variable importance in the projection (VIP) is an important output of PLSR and represents the relative importance of independent variables in explaining the dependent variable in the model. When VIP ≥ 1, the independent variable is considered important and that it significantly affects the dependent variable [47]. The independent variable could be considered the dominant driver of the dependent variable when its VIP value is maximum and greater than 1. PLSR is widely used in the study of vegetation response to climate change and has played an important role in identifying influential factors in both lag and driving force analyses [36,47,48]. In this study, all variables were standardized before PLSR analysis.

Sensitivity of SOS to Climate Change
The process for calculating the sensitivity of the SOS to climate change ( Figure 2) is as follows: (1) Screening of indicators: Meteorological and biological factors can affect SOS. Therefore, we selected Tmax, Tmin, Pre, Aet, melting snow water equivalent (Ms), effective precipitation (Water = Pre − Aet), water balance index (Wbi = Pre − Aet − Ms), Sm, and Srad as meteorological variables in this study covering the six months of preseason; the lengths of the growing season (LOS), EOS, and SOS of the previous year were selected as biological variables. We performed a PLSR for each pixel using 57 (9 × 6 months + 3) indicators as independent variables and SOS as the dependent variable. Ms was calculated as the difference between the month's snow water equivalent (Swe) and that of the previous month. Based on our PLSR results in the Northern Hemisphere and the specific physical significance, we selected Tmax, Tmin, Wbi, and Srad as the final climatic indicators and the previous year's EOS and LOS as the final biological indicators affecting SOS in this study. (2) Optimal preseason determination: The end date of the preseason is usually fixed as the SOS date. Therefore, the start date should be optimized to strengthen the correlation between preseason climatic factors and phenology [29,33]. Based on the PLSR results, we defined preseason in this study as the period during which the preseason climatic factors (preseason average Tmax and Tmin, cumulative Wbi, and cumulative Srad at a time step of one month) showed the highest contribution to SOS dynamics in each pixel, that is, when the VIP value was at its maximum. The weather data for the month in which the SOS was located were used to calculate the preseason length when the SOS date occurred after the 15th [16]; this approach retains the influence of biological factors with a total of 26 (4 × 6 months + 2) independent variables. In this step, we obtained 20 optimal preseason datasets, each of which included the best preseason duration of Tmax, Tmin, Wbi, and Srad for each pixel in the Northern Hemisphere, one on a 33-year time window and the remaining 19 on a 15-year time window. (3) Standardized sensitivity index: To analyze the spatial distribution and temporal variability of SOS sensitivity to climatic factors in the Northern Hemisphere during 1982-2014, we conducted a PLSR for each pixel between SOS and the optimal preseason meteorological (mean Tmax, mean Tmin, cumulative Wbi, cumulative Srad) and biological (EOS and LOS of the previous year) factors with time windows of 33 and 15 years, respectively. MC is the standardized sensitivity index, and VIP is the relative importance of each climatic factor. S Com is the sum of the MC values of the four climatic factors. Finally, we obtained 20 sensitivity datasets, each of which included the sensitivity of SOS to preseason Tmax, Tmin, Tem, Wbi, Srad, and Com for each pixel in the Northern Hemisphere.

Effect Evaluation of the Sensitivity Index
A high standard deviation indicates dispersion across a wide range of values in a time series [49]. The standard deviation (SD) index, which only represents the stability of SOS under the influence of climate factors-was used with the standardized sensitivity index for the Spearman correlation analysis. Higher absolute values of the correlation coefficient indicate a stronger effect of the sensitivity index. The formula for calculating the SD index is as follows: where Std is the standard deviation of the SOS time series from 1982 to 2015, indicating the stability of SOS; R 2 is the determination coefficient of PLSR, which represents the contribution of climatic and biological factors to SOS dynamics; and W met is the weight of the meteorological factors. The SD index can therefore be used to measure the SOS fluctuations caused by climate change, with larger SD values indicating weaker SOS resistance to climate change.

Trend Analysis and Significance Testing
The Mann-Kendall (MK) trend test is a commonly used method for measuring hydrological, meteorological, and other non-normally distributed time-series trends. We used the Then-Sen (TS) estimator and the non-parametric MK test for trend analysis and significance testing of the sensitivity index. The TS-MK method minimizes the noise related to parametric uncertainty in the regression equation [50]. The climate tendency rate is a method to describe the long-term meteorological trends; the linear tendency method was used to calculate the climate tendency rate of each climatic factor for each pixel, where p < 0.05 indicates a significant change. The sensitivity intensity in this study was expressed as the absolute value of the sensitivity for trend analysis and driving force analysis.

Driving Force Analysis
The PLSR method was used to determine the influence of spring (March to May) climate (Tmin, Tmax, Wbi, Srad), spring climate tendency rate (Tmin', Tmax', Wbi', Srad'), biological (mean value of NDVI and SOS), and geographical (altitude and latitude) factors on the spatial distribution of S Com in the Northern Hemisphere during 1982-2014 with time windows of 33 and 15 years. The PLSR method was also used to explore the relative importance of climatic (average temperature in spring, average temperature in winter-December to February, accumulated Wbi in spring and winter-December to May, accumulated Srad in spring and winter) and biological (mean value of SOS with a 15-year time window) factors driving the temporal variation of S Com in each pixel. The values of MC and VIP obtained by the PLSR represent the respective driving direction and relative importance of each driving factor that affected the spatio-temporal variation of sensitivity.

Selection of Indicators and Preseason Duration for Sensitivity Assessment
Reasonable light, temperature, and moisture indicators can help us obtain more accurate SOS sensitivity to climate change. Statistical results on the influence of meteorological variables covering the six months of the preseason and previous year's phenology on SOS in the Northern Hemisphere (Figure 3a) showed that meteorological variables had the strongest influence on SOS during the month in which SOS occurred, and the intensity of the impact then decreased with increasing months. Tmax had the strongest impact on SOS followed by Tmin. The correlation between SOS and Tmax and Tmin may be the opposite (Figure 3b). Among the water-related indicators, Aet had the strongest effect on SOS followed by Ms. This may be due to the direct effect of temperature on evapotranspiration and snowmelt, indicating an SOS response to temperature. The water balance index (VIP Wbi = 1.26) considering snowmelt had a stronger influence on SOS dynamics than that of effective precipitation (VIP Water = 1.23). Sm was not considered because its VIP value was lower than that of Wbi and did not show a cumulative effect. Therefore, we selected Wbi as the moisture indicator. Srad had a significant effect on the SOS (VIP Srad = 1). To exclude the influence of biological factors, we included preseason EOS and LOS, which significantly (VIP ≥ 1) influenced SOS-as covariates in the sensitivity index calculation. Finally, we selected preseason Tmax, Tmin, Wbi, and Srad as climate indicators to construct the standard sensitivity index. When screening for meteorological variables, we found that the response of SOS to climatic factors (Tmax, Tmin, Srad) of neighboring and distant months was reversed (Figure 3b). Preseason duration implies the combined effects of spring and winter climate change on SOS, especially from winter warming. We calculated the preseason durations of each pixel for 19 time periods in 15-year time steps and found significant spatial variability in the preseason vegetation duration ( Figure S1 of Supplementary Materials), but fewer regions with significant temporal variability ( Figure S2 of Supplementary Materials). The application of the optimal preseason could reduce the sensitivity error caused by the change in the preseason over time.

Sensitivity of SOS to Preseason Meteorological Factors
The sensitivity statistics for the climatic region and plant type are shown in Figure 4. The S Com value for each biome occurred in the following order: forest > shrub > grassland > farmland > sparse vegetation, with the highest sensitivity in mixed forests (MF) (Figure 4a). S Com is also more sensitive in wet and cool zones than in dry, warm, and cold zones. The cool, temperate, and moist (CTM) zones, where most of the Northern Hemisphere deciduous forests are located, are the most sensitive, and warm, temperate, and dry (WTD) zones dominated by grass and sparse vegetation are the least sensitive (Figure 4b). Overall, SOS in the Northern Hemisphere ecosystem showed the highest sensitivity to Tmax (−0.   Table 1; N stands for the Northern Hemisphere (latitude ≥ 30 • ).
The spatial distribution of the SOS standardized sensitivity to preseason meteorological factors in the Northern Hemisphere from 1982 to 2014 is shown in Figure S3 of Supplementary Materials. The highest SOS sensitivities to preseason Tmax and Tmin (S Tmax and S Tmin ) were concentrated in the mid-latitude and non-farming areas. We also found that 19.3% of the pixels showed opposite SOS responses to preseason Tmax and Tmin ( Figure S4 of Supplementary Materials). The highest SOS sensitivities to preseason Wbi and Srad (S Wbi and S Srad ) were mainly distributed at the high latitudes. The proportions of pixels significantly affected by temperature, moisture, and radiation were 88%, 46%, and 29%, respectively; therefore, the spatial distribution of S Com was similar to that of S Tem ( Figure S3c-f of Supplementary Materials). The spatial distribution of each sensitivity's contribution to S Com (Figure 5d) also shows that S Tem limits S Com in most pixels of the mid-latitudes, whereas S Wbi and S Srad are the dominant drivers of S Com at high latitudes and high elevations. S Com in the Northern Hemisphere was dominated by S Tem , S Wbi , and S Srad from 1982 to 2014, contributing 74.9%, 15.7%, and 9.4%, respectively.

Temporal and Spatial Characteristics of SOS Sensitivity to Preseason Climatic Factors
We calculated the sensitivity of each pixel over 19 periods from 1982 to 2014 using a time window of 15 years. The spatial distribution of S Com varied significantly between the period of rapid warming (1984-1998) and interrupted warming (1998-2012), with the proportion of significantly negative sensitive pixels decreasing from 81% to 76%. The proportion of pixels in which S Tem , S Wbi , and S Srad dominated S Com was 66.8%, 19.0%, and 14.2% during the rapid warming period, respectively, which changed to 61.6%, 22.2%, and 16.1% during the warming interruption. The region in which sensitivity was dominated by Wbi and Srad expanded (Figure 5b,c,e,f). We observed obvious spatial heterogeneity in the time-varying sensitivity trends ( Figure 6 and Figure S5). A total of 39.5% of pixels showed significant changes in S Com (22.2% decrease, 17.3% increase), whereas 53.3% of pixels showed a decreasing trend in the Northern Hemisphere (Figure 6f). The trends in single-factor sensitivity and integrated sensitivity are not spatially consistent, even for the S Tem that contributes the most to S Com (Figure 6e Figure S6 of Supplementary Materials), the sensitivity of SOS to preseason Com, Tmax, Tmin, and Tem was significantly reduced. The sensitivity to preseason Wbi showed a tendency to first decrease and then increase, whereas the opposite was true for preseason Srad, although there were no significant changes. However, the proportion of pixels dominated by S Srad and S Wbi in S Com increased significantly with time.

Effectiveness Evaluation of the Standardized SOS Sensitivity Index
The correlation between the standardized sensitivity index and the SD index ( Figure S7 of Supplementary Materials) in most climate regions and plant types is expressed as Com > Tem > Tmax > Tmin > Srad > Wbi. That is, the comprehensive sensitivity index was superior to the single-factor sensitivity index and showed stronger applicability in boreal moist, polar alpine climate regions (BM = −0.60, ET = −0.61) and deciduous needleleaf forests (DNF = −0.67).

Drivers of the Spatial Pattern of S Com
Several factors affected the spatial pattern of S Com in the Northern Hemisphere from 1982 to 2014 (Figure 7). The spatial pattern was most strongly influenced by the NDVI (VIP = 1.17, MC = 0.05). Moreover, earlier SOS dates corresponded to higher sensitivities (MC = −0.02). SOS sensitivity increased with digital elevation model altitude (DEM, MC = 0.04) and decreased with latitude (LAT, MC = −0.02); the effects of LAT on the spatial distribution of sensitivity were significant (VIP = 1.11). Tmin was the dominant climatic control on the spatial pattern of S Com (VIP = 1.14, MC = −0.014), followed by Tmax (VIP = 1.12, MC = −0.004). Wbi had the weakest effect on the spatial pattern of sensitivity (VIP = 0.48); none of the rates of change of climatic factors had a significant effect on the sensitivity spatial pattern (VIP < 1). However, owing to the complexity of climate and land cover, the impacts of the different drivers on sensitivity were not consistent across different climate regions ( Figure S8

Drivers of the Temporal Characteristics of S Com
Biological factors had a stronger effect on the temporal variability of S Com than did climate change (Figure 8). SOS with the maximum VIP value was the dominant driver of the temporal variability of S Com in the Northern Hemisphere, followed by spring temperatures (Tems) and Wbi. Winter temperatures (Temw) and Srad had the weakest influence. The importance of Wbi on the temporal variability of S Com was greater in the arid zones (WTD and CTD) (Figure 8a). According to the spatial distribution of the drivers influencing the temporal variation of S Com ( Figure S10 of Supplementary Materials), regions significantly affected by SOS and climate factors were not concentrated. The proportion of pixels showing positive and negative influences on sensitivity for each driving factor in the different climate regions was relatively balanced at~50% and did not exceed 57.5% (Figure 8b). Considering climate alone, the pixel proportions in which Tem, Wbi, and Srad dominated the temporal variability of S Com in the Northern Hemisphere were 51.6% (Tems = 25.5%, Temw = 26.1%), 26.1%, and 22.3%, respectively. The pixel proportions dominated by Wbi in arid (WTD and CTD) and cold (BD, BM, ET) areas were >25%, and Srad performed best in wet and ET climates (Figure 8c). Our results suggest that temperature change is the dominant climatic control on the temporal variability of SOS sensitivity, with important contributions from water balance and radiation. Most notably, the effects of spring and winter warming on the temporal variability of sensitivity in different regions of the Northern Hemisphere were superimposed or offset. The proportions of pixels where spring/winter warming jointly promoted and inhibited sensitivity were 23.7% and 23.3%, respectively, whereas 53.1% of pixels showed the opposite effect from spring/winter warming ( Figure S11a of Supplementary Materials). The proportion of pixels where spring/winter warming had an opposite effect on sensitivity increased with increasing latitude and decreasing temperature in dry climate regions (WTD = 53.0%, CTD = 53.7%, BD = 58.1%). However, the opposite trend was observed in humid regions (WTM = 51.2%, CTM = 50.4%, BM = 48.8%). In addition, the aforementioned results also showed that the proportion of the asymmetric response of S Com to spring/winter temperature changes was smaller in the humid zone than that in the dry zone. This indicate that more moisture might mitigate the asymmetric response of S Com to spring/winter warming ( Figure S11b of Supplementary Materials).

Rationality of the Standardized SOS Sensitivity Index and Its Evaluation Results
Both existing studies [18,[51][52][53] and the present study showed an asymmetric effect of daytime and nighttime (spring and winter) warming on SOS. However, this asymmetry has rarely been considered in previous SOS sensitivity studies [23,24,33,37,54,55]. Biological and environmental factors combine to regulate the SOS. Interactions between phenological events in spring and autumn may change the phenological response to ongoing climate warming [10], and interspecies interactions may affect SOS [56]. Our results showed that preseason EOS and LOS significantly affect the SOS. Studies have shown that changes in SOS sensitivity to temperature are associated with preseason duration [16,33]. For the aforementioned facts, we constructed a composite sensitivity index to characterize the complex response of spring phenology to climate change (Tmax, Tmin, Wbi, Srad). The effects of the previous year's phenology (EOS, LOS) and optimal preseason duration were incorporated in the calculation process to exclude the effects of non-meteorological factors on the sensitivity index. The results demonstrated that the composite index was superior to the single-factor index in expressing the influence of SOS due to climate change. Temperature, especially the preseason Tmax, is considered to be the main controlling factor for SOS, as warmer temperature breaks the ecodormancy earlier with global warming [10,51,57]; our results also show that SOS sensitivity to Tmax is highest in the Northern Hemisphere ecosystems and most vegetation types and climate type ecosystems, followed by Tmin.
However, we also found that SOS changes were dominated by moisture and radiation in 15.7% and 9.4% of the grid points, respectively. The control of SOS by moisture and radiation can be explained by their indirect effect on the thermal demand of SOS [10]. The interaction between radiation, temperature, and water affects the sensitivity of SOS to climate change [17,25]. S Com can reflect the comprehensive impact of climate change (including changes in light, temperature, and water and their interactions) on SOS.

Driving Force Analysis of the Spatial Distribution of S Com
The combined interaction between climate characteristics, climate tendency, and biological and geographical factors influences the spatial pattern of sensitivity. However, few studies have examined these factors together to identify the dominant drivers for the spatial distribution of SOS sensitivity. S Tem was the highest in the forests [55]. SOS in warmer locations of the Northern Hemisphere was found to be less responsive to climate change, and S Tem increased with altitude (cold climate) [33]. Wang et al. [23] found that the temperature sensitivity decreased with increasing spring temperature. Moreover, vegetation at lower latitudes (or with earlier SOS) was found to be more sensitive to changes in spring temperature [37,54]. Climatic tendency is highly correlated with SOS [20]. The findings of this study are consistent with these observations. In addition, we found that the multi-year average NDVI (an index reflecting the status of surface vegetation cover) had the greatest influence on the spatial distribution of S Com . S Com decreased with increasing latitude (or decreasing temperature) in this study, despite warmer regions (lower latitudes) showing lower sensitivity. This may be due to the importance of vegetation coverage (NDVI) on the spatial distribution of sensitivity, with higher sensitivity at higher coverage. Sparse vegetation, grasslands, and tundra dominate the high latitudes, high altitudes, and cold regions; this prevents SOS from closely tracking changes in temperature to avoid frost risks [58].
The spatial distribution of S Com showed diametrically opposed responses to drivers in different climatic regions ( Figure S5 of Supplementary Materials) attributable to the regional characteristics of vegetation and climate types in the Northern Hemisphere, which might also be due to the interaction of hydrothermal conditions to influence the distribution of SOS sensitivity [25]. With the exception of LAT and Srad, the influence of biological factors (SOS, NDVI), climate characteristics (Tmax, Tmin, Wbi), climate trend rate (Tmax', Tmin', Wbi', Srad'), and geographical factors (DEM) on the spatial pattern of sensitivity has weakened in recent years ( Figure S6 of Supplementary Materials). The reasons for this phenomenon are as follows: First, the vegetation and the communities in which they are located have adapted to climate warming [31,56]. Furthermore, the rate of warming has decreased, resulting in no trends in SOS during the global warming hiatus [31,56]. Finally, a faster rate of warming at higher altitudes has led to SOS at different altitudes becoming more uniform [15,59].

Driving Force Analysis of the Temporal Variability of S Com
According to global average surface temperature data, a global warming interruption occurred between 1998 and 2012 [60,61]. To explore the effect of this warming slowdown on SOS sensitivity, we analyzed the temporal variability and dominant drivers of SOS sensitivity using a time window of 15 years. The simulated effect (R 2 = 0.57) of preseason climate change and phenology on SOS dynamics was significantly enhanced at the 15-year time step (Figure 4g,h,j) relative to the 33-year time window (R 2 = 0.40), which allowed us to obtain more accurate SOS sensitivity, as the SOS response to preseason climate change (especially temperature change) is nonlinear [29].
SOS sensitivity to climate warming generally shows a decreasing trend in the Northern Hemisphere [29]. The present study also showed a significant decreasing trend in the median of all pixels S Com (or S Tem ) in the Northern Hemisphere. However, we found that only 53.3% of the pixels showed a decrease in S Com ; 22.2% were statistically significant ( Figure 5 and Figure S2a). This may suggest that the decrease in the sensitivity trend exceeded the increase on the pixel scale. The significant reduction in SOS sensitivity in a small proportion of vegetation leads to a decreasing trend of S Com (or S Tem ) in the Northern Hemisphere, with vegetation in wetter/cooler climate zones and shrubs contributing more to the decrease in sensitivity. In all climate regions other than WTM, we found that biological factors (SOS) exert more influence on the temporal variability of SOS sensitivity relative to each climatic factor (Tems, Temw, Wbi, Srad), suggesting that vegetation autonomously regulates SOS to adapt to climate change [31]. Vegetation adaptation is the most important factor influencing temporal variation in S Com .
Temperature may be an important climate driver dominating temporal changes in sensitivity (5.16% of the pixels). Strong uncertainty remains in the main climatic factors that driver the temporal variability of sensitivity. In different climatic regions and the Northern Hemisphere, the proportion of positive and negative responses of sensitivity to various driving factors (biological and climatic) was relatively balanced (Figure 7b). This may be due to the spatial heterogeneity ( Figure 5 and Figure S2) of the sensitivity variation (only 53.3% of the pixels showed a decrease in sensitivity), or the asymmetric response of S Com to spring/winter warming (53.1% of pixels showed the opposite effect). More likely, it is due to the complex interaction of drivers [10]. For example, radiation showed clear timevarying characteristics [62], and the proportion of pixels dominated by radiation reached 22.3% in this study. Moreover, strong solar radiation in spring causes snow to melt, leading to higher temperatures in the soil relative to the atmosphere. The sensitivity of SOS may be affected by the influx of snowmelt water into partially frozen soils, as this promotes root activity before temperatures reach 0 • C [63]. The effect of Wbi on the temporal variability of S Com was greater than or equal to the effect of Tems in arid and cold climate regions. Water shortage in these regions is the dominant limiting factor in the development of spring phenology [64]. We also found that more moisture might mitigate the asymmetric response of S Com to spring/winter warming. Srad and Temw showed weaker influences on the temporal variability of sensitivity relative to Wbi and Tems; however, the proportion of pixels in which the temporal sensitivity was dominated by Temw still reached 26.1%. Studies have shown that warming in winter and early spring shortens the cold storage period, which reduces sensitivity, as the cold storage requirements of certain plant types are not met [65]. Although the response of sensitivity to climatic factors is complex, the main driving factors still have different emphases in different climatic regions.

Limitations
As the four methods to estimate SOS were based on satellite data, some uncertainties may arise due to resolution and data quality limitations. Moreover, although we selected pixels with no significant change in land cover type for statistical analysis, they may still include some impacts from human activity. We studied the spatial and temporal patterns of spring phenology sensitivity and analyzed the relative importance of each driving factor to the spatial and temporal distribution of sensitivity. However, we do not understand the underlying causes and ecological processes that lead to these patterns. Meteorological and phenological observations with higher spatio-temporal resolution and manipulatable physiological experiments can further our understanding of the underlying mechanisms of the SOS response to climate change.

Conclusions
We applied the PLSR method as the core algorithm and used meteorological and NDVI-extracted phenological data to screen the most suitable biological and climatic indicators affecting SOS in the Northern Hemisphere. We then constructed a standardized SOS sensitivity index for climate change based on the optimal preseason duration of the climatic indicators. The index can determine the weight and direction of individual climate factors as well as their combined effects on SOS. Based on the sensitivity evaluation results, we analyzed the spatial and temporal characteristics of SOS sensitivity and its dominant drivers in the Northern Hemisphere from 1982 to 2014. This study showed that the combined sensitivity of SOS to climate change (S Com ) is most influenced by pre-season temperature sensitivity. However, because of the asymmetric response of SOS to Tmax/Tmin and the non-negligible moderating effect of Wbi and Srad on SOS, S Com was more effective in expressing the effect of climate change on SOS compared with the single climatic factor. Vegetation cover (or type) was the dominant factor influencing the spatial pattern of SOS sensitivity, followed by spring temperature (Tmin > Tmax); the weakest was Wbi. Forests had the highest S Com absolute values. In recent years, a significant decrease in the sensitivity of some vegetation has led to a decreasing trend of sensitivity in the Northern Hemisphere, with vegetation in wetter/cooler climate zones and shrubs contributing more to the decrease in S Com . Biological factors (vegetation adaptation) exerted a greater influence on the temporal variability of SOS sensitivity relative to each climatic factor (Tems, Temw, Wbi, Srad). Although temperature remains the main climatic factor driving temporal changes in S Com , the effects were asymmetric between spring and winter. More moisture might mitigate the asymmetric response of S Com to spring/winter warming. Despite the complex response of sensitivity to climate change, the dominant factors in different climatic regions still have different emphases.
Identifying systematic differences in phenological climate sensitivity would facilitate the development of indicators and estimates of vulnerability for conservation and national adaptation programs [32]. The performance of the current phenological models is unsatisfactory [10]. The quantitative evaluation of sensitivity and driving force analysis of SOS in this study provides not only a new method for identifying the combined sensitivity of springtime phenology to climate change but also a basis for the construction of better phenological models to improve phenological prediction performance.  Figure S10: The spatial distribution of the driving forces influencing the temporal variability of S Com in the Northern Hemisphere; Figure S11: Effects of asymmetric warming in spring and winter on the temporal variation of S Com .
Author Contributions: All authors contributed meaningfully to this study. K.L., J.Z. designed the research; C.W., Q.S. performed the research, analyzed the data; G.R., Z.T., X.L. provided methodology support and revised the paper; K.L. drafted the manuscript, and all authors reviewed and edited the writing. All authors have read and agreed to the published version of the manuscript.