Precipitation Dominates the Relative Contributions of Climate Factors to Grasslands Spring Phenology on the Tibetan Plateau

: Temperature and precipitation are the primary regulators of vegetation phenology in temperate zones. However, the relative contributions of each factor and their underlying combined effect on vegetation phenology are much less clear, especially for the grassland of the Tibetan Plateau To quantify the contribution of each factor and the potential interactions, we conducted redundancy analysis for grasslands spring phenology on the Tibetan Plateau during 2000–2017. Generally, the individual contribution of temperature and precipitation to grasslands spring phenology (the start of growing season (SOS)) was lower, despite a higher correlation coefﬁcient, which further implied that these factors interact to affect the SOS. The contributions of temperature and precipitation to the grasslands spring phenology varied across space on the Tibetan Plateau, and these spatial heterogeneities can be mainly explained by the spatial gradient of long-term average precipitation during spring over 2000–2017. Speciﬁcally, the SOS for meadow was dominated by the mean temperature in spring (T spring ) in the eastern wetter ecoregion, with an individual contribution of 24.16% ( p < 0.05), while it was strongly negatively correlated with the accumulated precipitation in spring (P spring ) in the western drier ecoregion. Spatially, a 10 mm increase in long-term average precipitation in spring resulted in an increase in the contribution of T spring of 2.0% ( p < 0.1) for meadow, while it caused a decrease in the contribution of P spring of − 0.3% ( p < 0.05). Similarly, a higher contribution of P spring for steppe was found in drier ecoregions. A spatial decrease in precipitation of 10 mm increased the contribution of P spring of 1.4% ( p < 0.05). Considering these impacts of precipitation on the relative contribution of warming and precipitation to the SOS, projected climate change would have a stronger impact on advancing SOS in a relatively moist environment compared to that of drier areas. Hence, these quantitative interactions and contributions must be included in current ecosystem models, mostly driven by indicators with the direct and the overall effect in response to projected climate warming.


Introduction
Climate change was now documented in all major regions of the world and affected vegetation structure and function as well as public health [1][2][3][4]. As the dominant vegetation type on the Tibetan Plateau, grassland plays important roles in ecological balance and provides an important buffer against climate change [5,6], and the grassland phenology was identified as a key indicator of climate-vegetation interactions and reflects the response of living systems to climate change [7,8]. In particular, the spring phenology of grassland can provide insight into higher level dynamics of plant functioning and feedback to climate change through exchange of energy, the hydrological cycle, and carbon uptake, which is more relevant to global climate and regional natural environmental changes [9][10][11][12][13].
Temperature and precipitation are regarded as the primary regulators of plant phenology in the temperate zone. Highly credible evidence shows that global warming and regional changes in precipitation strongly affected vegetation phenology [12,[14][15][16][17]. Although there were many studies focused on understanding the mechanisms and influences of climate warming and regional precipitation on vegetation phenology, there are still a lot of uncertainties in the responses of the start of growing season (SOS) to these factors. As proof, Lucht et al. [18] and Xu et al. [19] showed that the growth of vegetation would be facilitated, because global warming eased climatic constraints. However, according to Jeong et al. [20], the vegetation spring phenology and preseason mean temperature are correlated across the Northern Hemisphere with a correlation coefficient varying from −0.3 to −0.7, suggesting that global warming would not necessarily induce greater advance in SOS. In addition, Yu et al. [13] showed that winter warming could delay vegetation spring phenology due to the adaptive ability of chilling requirements, while Suonan et al. [21] showed that winter warming increased the soil temperature, and then advanced spring leaf out and flowering phenology. Some studies also showed that the daylength could coregulate the SOS through its interaction with temperature or its influence on the stomata aperture and photosynthetic active radiation [22,23]. Furthermore, water is also necessary for plant growth, which indicates that changes in regional precipitation should be another main driving factor that regulate vegetation phenology, particularly in water-limited areas [9,12,16,24]. However, the SOS's responses to preseason precipitation were also diverse. Dai et al. [25] found that precipitation and the first leafing date do not show a significant negative correlation. Other studies, on the other hand, reported that, in arid or semiarid regions, preseason precipitation indeed advanced vegetation SOS because it determines spring water availability [5,9,26]. Moreover, some studies argued that preseason precipitation would increase the heat demand and then influence the sensitivity of SOS to precipitation and temperature [16,27]. That is, along with direct impact, precipitation may exert indirect impacts on SOS. These direct and indirect impacts indicated that temperature and precipitation would complicatedly interact to affect vegetation SOS. It is difficult to quantify the individual contributions of temperature and precipitation to the interannual variation of vegetation phenology because of the nonlinear impacts and underlying interactions of climate factors.
As the world's largest plateau, the Tibetan Plateau is considered as one of the most vulnerable ecosystems because of its typical alpine climate and its strong sensitivity to climate change [5,[9][10][11][12][13]. Although numerous recent studies reported the interannual variation of SOS and its responses to warming and precipitation [10][11][12][13]28], large uncertainties still exist about the advanced (i.e., earlier) or delayed (i.e., later) pattern of the vegetation spring phenology and their responses to different factors. For instance, Piao et al. [29] showed that the spring phenology of alpine meadow would be delayed by the increased accumulated preseason precipitation occurring in five months before the vegetation onset, while Shen et al. [30] found that declines in spring precipitation would delay SOS rather than increase preseason precipitation. Moreover, Yu et al. [13] showed that winter warming could delay grassland spring phenology on the Tibetan Plateau, while Suonan et al. [21] showed that asymmetric winter warming advanced plant phenology in alpine meadow. Furthermore, these explanations were later questioned by some studies that argued that the vegetation spring phenology was more strongly affected by the interaction of warming and precipitation than the unitary effect of driving factors [5,9,11,16]. Even so, the interaction effects of climate factors on vegetation phenology remain poorly understood, and the individual contribution of each factor to grassland SOS across the Tibetan Plateau is not yet well understood. Thus, quantitative estimation of the relative contribution of critical climate factors and their interaction effects on grassland SOS on the Tibetan Plateau is critical and valuable. In this study, we first analyze how SOS changed in the grassland in different ecoregions of the Tibetan Plateau. We untangle the individual contribution and interaction of the mean temperature in winter (T winter ), mean temperature in spring (T spring ), and the accumulated precipitation in spring (P spring ) for interannual variations in SOS based Remote Sens. 2022, 14, 517 3 of 18 on redundancy analysis and identify temporal-spatial pattern of SOS responses to critical climate factors on the Tibetan Plateau. This is an important step toward understanding the mechanistic effects and underlying causes of climate change on vegetation phenology.

Study Area
The Tibetan Plateau lies at an average altitude of over 4000 m in southwest China (Figure 1), with a typical alpine climate environment and unique composition and distribution of alpine grassland, along with a low level of human interference. Across the entire plateau, the annual mean temperature varies from −15 • C to 10 • C, and the annual cumulative precipitation decreases from more than 1000 mm to 50 mm. Recently, the Tibetan plateau experienced substantial climate change characterized by significant warming [9][10][11][12]. Meanwhile, numerous studies showed that vegetation was highly sensitive to climate changes on the Tibetan Plateau [9,12,15,27]. Hence, climate change and its increasingly pronounced effects on Tibetan Plateau grassland phenology are becoming an issue of global concern.
Remote Sens. 2022, 14, 517 3 of 19 contribution and interaction of the mean temperature in winter (Twinter), mean temperature in spring (Tspring), and the accumulated precipitation in spring (Pspring) for interannual variations in SOS based on redundancy analysis and identify temporal-spatial pattern of SOS responses to critical climate factors on the Tibetan Plateau. This is an important step toward understanding the mechanistic effects and underlying causes of climate change on vegetation phenology.

Study Area
The Tibetan Plateau lies at an average altitude of over 4000 m in southwest China (Figure 1), with a typical alpine climate environment and unique composition and distribution of alpine grassland, along with a low level of human interference. Across the entire plateau, the annual mean temperature varies from −15 °C to 10 °C, and the annual cumulative precipitation decreases from more than 1000 mm to 50 mm. Recently, the Tibetan plateau experienced substantial climate change characterized by significant warming [9][10][11][12]. Meanwhile, numerous studies showed that vegetation was highly sensitive to climate changes on the Tibetan Plateau [9,12,15,27]. Hence, climate change and its increasingly pronounced effects on Tibetan Plateau grassland phenology are becoming an issue of global concern. Figure 1. Location, ecoregions, and vegetation distribution for Tibetan Plateau. (I) Indicates Eastern Qinghai-Qilian montane steppe ecoregions; (II) indicates Golog-Nagqu high-cold shrub-meadow ecoregions; (III) indicates southern Qinghai high cold meadow steppe ecoregions; (IV) indicates Qiangtang high-cold steppe ecoregions; and (V) indicates Southern Xizang montane shrub-steppe ecoregions [31]. Pixels with annual maximum normalized difference vegetation index (NDVI) greater than 0.15 and mean NDVI of July to August larger than 1.35 × NDVI during November to February are shown.

Datasets
During the past two decades, the monitoring of vegetation phenology at large-scale regional units was widely achieved using satellite remote sensing. In this study, we first produced an NDVI (normalized difference vegetation index) dataset by using MODIS (moderate resolution imaging spectroradiometer) surface reflectance products with a spatial resolution of 500 m and a compositional period of 8 days from 2000 to 2017 [32], obtained from the LPDAAC (land processes distributed active archive center,  [31]. Pixels with annual maximum normalized difference vegetation index (NDVI) greater than 0.15 and mean NDVI of July to August larger than 1.35 × NDVI during November to February are shown.

Datasets
During the past two decades, the monitoring of vegetation phenology at large-scale regional units was widely achieved using satellite remote sensing. In this study, we first produced an NDVI (normalized difference vegetation index) dataset by using MODIS (moderate resolution imaging spectroradiometer) surface reflectance products with a spatial resolution of 500 m and a compositional period of 8 days from 2000 to 2017 [32], obtained from the LPDAAC (land processes distributed active archive center, https://lpdaacsvc. cr.usgs.gov/appeears/. Accessed 21 April 2021). To match the spatial resolution of the climatic factors, we then resampled these NDVI datasets to a resolution of 1 km and used them to determine vegetation SOS.
Information about the meadow and steppe distribution on the Plateau was obtained from the 1:1,000,000 vegetation map of China [33], created by the Chinese Academy of Sci-Remote Sens. 2022, 14, 517 4 of 18 ences, using datasets from field observations, aerial photos, TM (landsat thematic mapper) and ETM (enhanced thematic mapper) satellite images [34,35]. Climate factors at meteorological stations, including the mean air temperature ( • C), accumulative precipitation (mm) in spring (March-May), and the mean air temperature ( • C) in winter (December-February) were obtained from the China Meteorological Data Service Center. Based on the DEM (digital elevation model) data with 1 km spatial resolution produced by United States Geological Survey (USGS) and the 71 meteorological stations with complete records on the Tibetan Plateau (Figure 1), we calculated the spatial grid distribution of accumulated precipitation and mean temperature in spring using ANUSPLIN interpolation. Specifically, considering there were fewer than 2000 meteorological stations, we first employed the SPLNEA module to calculate the surface fitting parameters and their error covariance. Based on the latitude, longitude, and elevation information in a grid cell provided by DEM, the fitted climatic values and corresponding predicted standard error in the grid were then calculated by using the LAPGRD module [36][37][38].

Preprocessing of NDVI Data
Due to the inevitable impact of bare soils, sparsely vegetated grids, atmospheric contamination, etc., there were still many spurious changes in the NDVI time series [9,34,39]. Hence, to improve the quality of annual NDVI time series, the Savitzky-Golay filter for NDVI time series was conducted first. Additionally, we only selected pixels that met the following principles simultaneously to analyze the relative contributions of climate factors and their underlying combined effect on grassland spring phenology. First, the annual maximum NDVI must be greater than 0.15, and then the mean NDVI of July to August must be larger than 1.35 × NDVI during November to February. The spatial distribution of valid pixels is displayed in Figure 1.

Estimation of Grassland Spring Phenology
In this section, the 18-year averaged NDVI value with a compositional period of 8 days was calculated from the entire data set during 2000-2017. Based on this long-term averaged NDVI time series data, subsequent analysis and calculation were conducted. Specifically, two different methods were employed to detect the grassland spring phenology on the Tibetan Plateau.

Calculating the NDVI Thresholds Based on the Relative Rates of Changes
Based on the long-term averaged NDVI time series data, we first calculated the relative rates of changes of NDVI (NDVI_rate) for each valid pixel (Equation (1)). The dynamic NDVI threshold for each pixel was defined based on the highest positive relative change rate; that is, the maximum NDVI_rate corresponding to NDVI is the dynamic NDVI threshold for detecting SOS ( Figure 2).
where NDVI_rate(t) is the relative rate of change in NDVI, t is the time at 8-day intervals, and NDVI(t) is the averaged 8-day NDVI value at t time.

Polyfit Method
The polyfit method was initially applied to determine the vegetation spring phenology by Piao et al. [29]. Polyfit was widely shown to be one of the most efficient methods for the extraction of the vegetation phenology [8,29]. In general, the growing season of the grassland ranges from April to October in most areas of the Tibetan Plateau; that is, the NDVI value of grassland increased first, and then decreased with the increasing numbered day of the year during the growing season. Therefore, a 6-degree polynomial with a least-square analysis can be performed to reconstruct the daily NDVI time series from the original 8-day NDVI data during January to September. The polynomial function is shown in Equation (2): NDVI(t) = a0 + a1 × t + a2 × t 2 + a3 × t 3 + a4 × t 4 + a5 × t 5 + a6 × t 6 (2) where NDVI(t) is the NDVI at t time fitted by the polynomial equation. t is the numbered day of the year. a0, a1, a2, a3, a4, a5, a6 are the fitted coefficient derived from least-square regression analysis. The multiple linear regression function in MATLAB R2019b was used to estimate these fitted coefficients.

Harmonic Analysis of Time Series (HANTS) Method
The HANTS method was also used to reconstruct the daily NDVI time series from the original 8-day NDVI data. HANTS is regarded as an improved algorithm of the fast Fourier transform (FFT), which was proved to be one of the most efficient methods for the reconstruction of NDVI time series and the extraction of vegetation phenology [8,40,41]. The points with lower values than the neighbors are filtered during HANTS processing, and they are then reconstructed following Equation (3): where NDVI(t) is the NDVI fitted by the HANTS model, t is the numbered day of the year. a0 and ai are the fitted coefficient of the HANTS model. n is the frequency, which is set to 1 for the grassland on the Tibetan Plateau. φi is the phase of the NDVI time series. ωi is set

Polyfit Method
The polyfit method was initially applied to determine the vegetation spring phenology by Piao et al. [29]. Polyfit was widely shown to be one of the most efficient methods for the extraction of the vegetation phenology [8,29]. In general, the growing season of the grassland ranges from April to October in most areas of the Tibetan Plateau; that is, the NDVI value of grassland increased first, and then decreased with the increasing numbered day of the year during the growing season. Therefore, a 6-degree polynomial with a least-square analysis can be performed to reconstruct the daily NDVI time series from the original 8-day NDVI data during January to September. The polynomial function is shown in Equation (2): where NDVI(t) is the NDVI at t time fitted by the polynomial equation. t is the numbered day of the year. a0, a1, a2, a3, a4, a5, a6 are the fitted coefficient derived from least-square regression analysis. The multiple linear regression function in MATLAB R2019b was used to estimate these fitted coefficients.

Harmonic Analysis of Time Series (HANTS) Method
The HANTS method was also used to reconstruct the daily NDVI time series from the original 8-day NDVI data. HANTS is regarded as an improved algorithm of the fast Fourier transform (FFT), which was proved to be one of the most efficient methods for the reconstruction of NDVI time series and the extraction of vegetation phenology [8,40,41]. The points with lower values than the neighbors are filtered during HANTS processing, and they are then reconstructed following Equation (3): where NDVI(t) is the NDVI fitted by the HANTS model, t is the numbered day of the year. a 0 and a i are the fitted coefficient of the HANTS model. n is the frequency, which is set to 1 for the grassland on the Tibetan Plateau. ϕ i is the phase of the NDVI time series. ω i is set to 2π in our study because the grassland on the Tibetan Plateau has a unique seasonal cycle.
Based on the NDVI threshold and the fitted NDVI time series, the grassland SOS was determined for each pixel and each year. Specifically, the numbered day of the year when the pixel's NDVI value was first larger than the calculated NDVI threshold during the growing season was regarded as the grassland SOS.

Relationships between SOS and Climate
Generally, it is difficult to quantitatively estimate the individual contributions of temperature and precipitation for the interannual variation in vegetation phenology because of the nonlinear impacts and underlying interactions of climate factors. In this research, the redundancy analysis (RDA) was utilized to calculate the independent contribution of T winter , T spring and P spring to interannual variations in SOS [42,43]. The RDA is regarded as a direct extension of multiple linear regression and principal component analysis based on canonical multivariate analyses. Based on this method, multiple response variables could be regressed on multiple explanatory variables, and the explanation proportion of each explanatory variable to response variables can be estimated and summarized [44][45][46]. This processing of RDA was demonstrated by Eigen analysis in the following equation: where S yx is the covariance matrix among response variables and explanatory variables, S −1 xx is the standardized explanatory variables' inverse covariance matrix, λ k represents the eigenvalue of the corresponding axis k, I is unit matrix, and u k is normalized canonical eigenvectors.
In our study, a detrended correspondence analysis with the turnover units smaller than three was first conducted to confirm the suitability of an RDA. The adjusted coefficient of determination obtained from the Hierarchical Partitioning algorithm was regarded as the independent contribution of T winter , T spring and P spring to interannual variations in SOS (https://github.com/laijiangshan/rdacca.hp/, Accessed 27 July 2020). Specifically, we first conducted RDA for each ecoregion and the entire plateau to estimate the contribution of each factor at the ecoregional scale. Furthermore, the RDA was conducted for each valid pixel to analyze spatial patterns of relative contribution of each factor and determine the dominant factor and the proportion of corresponding dominating factors. The Bonferroni test and Monte Carlo permutation methods (Permutations = 499) were used to analyze the statistical significance [44,46,47]. In addition, to understand the spatial heterogeneity of the independent contribution of each factor, we performed linear regression in which the individual contribution of each climate factor was set to the dependent variable against the long-term average P spring and T spring across the Tibetan Plateau.

Interannual Variations in SOS
Across the Tibetan Plateau, most pixels' SOS showed an advanced trend (a negative linear trend) during 2000-2017, despite the spatial heterogeneity of trends in SOS revealed by the two different methods (Figure 3). That is, the timing of first and mean spring life history events advanced in time due to climate change. The average SOS of both methods advanced over about 72.42% of all the study areas, and most pixels advanced between 0-0.50 days per year −1 during the entire 18-year period (Figure 3c). Spatially, both methods showed widespread advancing trends over the central and eastern Plateau, while a slight delayed trend in spring phenology was found in a few areas of the central and western Plateau. More specifically, 19.80% of pixels advanced significantly (p < 0.05), which were mainly located in the eastern ecoregions. While the pixels with delayed trend were discretely distributed across the midwestern region of the Plateau, only 2.25% of pixels were statistically significant (p < 0.05). These similar proportions and spatial patterns were also revealed by the two different methods, and the results revealed an advanced trend of SOS from Polyfit and HANTS at 74.26% of pixels (pixels with a significant advanced trend Remote Sens. 2022, 14, 517 7 of 18 correspond to 24.58%, p < 0.05, Figure 3a) and 71.81% of pixels (22.03% with a significant advanced trend, p < 0.05, Figure 3b), respectively. mainly located in the eastern ecoregions. While the pixels with delayed trend were discretely distributed across the midwestern region of the Plateau, only 2.25% of pixels were statistically significant (p < 0.05). These similar proportions and spatial patterns were also revealed by the two different methods, and the results revealed an advanced trend of SOS from Polyfit and HANTS at 74.26% of pixels (pixels with a significant advanced trend correspond to 24.58%, p < 0.05, Figure 3a) and 71.81% of pixels (22.03% with a significant advanced trend, p < 0.05, Figure 3b), respectively. With respect to meadows and steppe in different ecoregions, all methods agreed on an advanced trend (negative) over the entire study period but without consistent significance (Table 1). Specifically, the trends in the multimethod averaged SOS were −0.30 (p < 0.01) and −0.20 (p < 0.01) days per year −1 for steppe and meadow in the entire plateau, respectively. Regarding the two different models' results, the trend of SOS for meadow across the entire plateau estimated from Polyfit was −0.21 (p < 0.01), while the result With respect to meadows and steppe in different ecoregions, all methods agreed on an advanced trend (negative) over the entire study period but without consistent significance (Table 1). Specifically, the trends in the multimethod averaged SOS were −0.30 (p < 0.01) and −0.20 (p < 0.01) days per year −1 for steppe and meadow in the entire plateau, respectively. Regarding the two different models' results, the trend of SOS for meadow across the entire plateau estimated from Polyfit was −0.21 (p < 0.01), while the result estimated from HANTS was slightly lower (Slope = −0.19, p < 0.01). Regarding steppe across the entire plateau, a similar pattern can be found: the advanced trend of SOS from Polyfit and HANTS was −0.32 (p < 0.01) and −0.27 (p < 0.01), respectively.  Concerning the different ecoregions, 42.88% of meadow and 77.53% of steppe showed a significant advanced trend (p < 0.05) (Figure 3d), with the mean linear trends of −0.30 (p < 0.01) and −0.82 (p < 0.01) days per year −1 in the I ecoregion in the eastern Tibetan Plateau. Moreover, the SOS for steppe in the III ecoregion also significantly advanced by −0.37 days year −1 (p < 0.01) according to the ensemble method, while no statistically significant trend for meadow in this ecoregion was observed at the p > 0.05 level. No statistically significant trend in SOS for meadow or steppe was found in IV and V ecoregions, which resulted from the substantial spatial heterogeneity of the trends in SOS, with a widespread delayed trend in these ecoregions.

Responses of SOS to Climate Factors
Across the Plateau, the interannual variations in SOS for meadow were dominated by T spring , with an individual contribution of 19.90% (p < 0.05) and partial correlation coefficient (PCC) of −0.48 (p < 0.05). By contrast, the steppe's SOS was more strongly affected by the interactions of precipitation and temperature, with an individual contribution of P spring and T spring of 15.90% (p < 0.1) and 13.30% (p < 0.1), respectively ( Table 2). Although higher PCC values were found between the SOS of steppe and P spring and T spring , with a PCC of −0.47 (p < 0.1) and −0.43 (p < 0.1), respectively, the individual contribution of P spring and T spring was still lower. Bold and symbols ** and * indicate significance levels at p < 0.05 and at p < 0.1, respectively. All indicates entire plateau. I, II, III, IV, and V indicate different ecoregions following Figure 1. T winter indicates mean temperature in winter; T spring indicates mean temperature in spring. P spring indicates accumulated precipitation in spring. Mean T winter , T spring , and P spring and SOS were first calculated for each ecoregion and entire plateau. RDA was then conducted at ecoregional scale.
Spatially, the individual contribution of T winter , T spring and P spring to interannual variations in SOS revealed a distinct east west disparity (Figure 4). The results showed that P spring most strongly dominated vegetation SOS over 18.44% of pixels, which were mainly located in the IV ecoregion (Qiangtang high-cold steppe ecoregion) in the western region of the Plateau. However, T spring dominated over 12.85%, mainly located in the I ecoregion (eastern Qinghai-Qilian montane steppe ecoregion) in the eastern part of the Plateau. In addition, the statistical analysis showed that T winter has the highest individual contribution in the central Plateau, with a proportion of 18.94% and 22.78% in the II (Golog-Nagqu high-cold shrub-meadow ecoregion) and III (southern Qinghai high-cold meadow steppe ecoregion) ecoregions in the central Plateau, respectively (Figure 4a).  Across the Tibetan Plateau, the responses of SOS to the interactions of temperature and precipitation were more prominent for steppe and meadow in the different ecoregions (Table 2 and Figure 5). Across the I ecoregion, 26.62% of the pixels for meadow were dominated by Tspring, with an individual contribution of 24.16% (p < 0.05) and PCC of −0.57 (p < 0.05). By contrast, the SOS of meadow in the IV ecoregion in the western part of the Plateau was dominated by Pspring, with an individual contribution of 18.47% (p < 0.05) and PCC of −0.58 (p < 0.05). With respect to the steppe on the Tibetan Plateau, 27.72% of the pixels for steppe in the I ecoregion were dominated by Pspring. The individual contribution of Pspring to the interannual variations in SOS in this ecoregion was 14.16% (p < 0.1), and the PCC was −0.41 (p < 0.1). While the interannual variations in SOS for steppe in the IV ecoregion in the western part of Tibetan Plateau were significantly affected by the underlying interactions of cumulative precipitation and mean temperature in spring. The individual contribution of Pspring and Tspring to the interannual variations in SOS was 20.98% (p < 0.05) and 12.78% (p < 0. 1), respectively. The PCC between SOS and Pspring and Tspring was −0.61 (p < 0.05) and −0.53 (p < 0.05), respectively. These significant individual contribution and Across the Tibetan Plateau, the responses of SOS to the interactions of temperature and precipitation were more prominent for steppe and meadow in the different ecoregions (Table 2 and Figure 5). Across the I ecoregion, 26.62% of the pixels for meadow were dominated by T spring , with an individual contribution of 24.16% (p < 0.05) and PCC of −0.57 (p < 0.05). By contrast, the SOS of meadow in the IV ecoregion in the western part of the Plateau was dominated by P spring , with an individual contribution of 18.47% (p < 0.05) and PCC of −0.58 (p < 0.05). With respect to the steppe on the Tibetan Plateau, 27.72% of the pixels for steppe in the I ecoregion were dominated by P spring . The individual contribution of P spring to the interannual variations in SOS in this ecoregion was 14.16% (p < 0.1), and the PCC was −0.41 (p < 0.1). While the interannual variations in SOS for steppe in the IV ecoregion in the western part of Tibetan Plateau were significantly Remote Sens. 2022, 14, 517 10 of 18 affected by the underlying interactions of cumulative precipitation and mean temperature in spring. The individual contribution of P spring and T spring to the interannual variations in SOS was 20.98% (p < 0.05) and 12.78% (p < 0. 1), respectively. The PCC between SOS and P spring and T spring was −0.61 (p < 0.05) and −0.53 (p < 0.05), respectively. These significant individual contribution and higher PCC indicated that interannual variations in SOS for steppe exhibited a close relationship with changes in climate conditions, illustrating the dependence of spring phenology on temperature and precipitation. In other words, an increase in preseason air temperature or cumulative precipitation would both correspond to a trend toward an earlier date of SOS. Our study also found that T winter was the most significant factor for steppe in the II and III ecoregions in the central Plateau, with an individual contribution of 31.32% (p < 0.05) and 11.40% (p < 0.1), respectively. temperature or cumulative precipitation would both correspond to a trend toward an earlier date of SOS. Our study also found that Twinter was the most significant factor for steppe in the II and III ecoregions in the central Plateau, with an individual contribution of 31.32% (p < 0.05) and 11.40% (p < 0.1), respectively.

Relationship between Individual Contribution and Climate Gradient
The spatial heterogeneities of responses of SOS to each factor can be fully explained by the long-term average precipitation gradient across the Tibetan Plateau ( Figure 6). Spatially, the SOS of grassland in the Tibetan Plateau was more strongly affected by mean temperature in the wetter areas (I and II ecoregions), while it was dominated by the accumulated precipitation during spring in drier areas (IV ecoregion). Considering the values averaged from the pixels with significant individual contribution at p < 0.05 level in different ecoregions in the Tibetan Plateau, a 10 mm increase in long-term average precipitation in spring responded to an increase in the individual contribution of Tspring of 1.70% (p < 0.05), while it caused a decrease in an individual contribution of Pspring of −0.70% (p < 0.1) (Figure 6a). In contrast, among different ecoregions on the Tibetan Plateau, this pattern linking the relative contribution of each climate factor to long-term average mean

Relationship between Individual Contribution and Climate Gradient
The spatial heterogeneities of responses of SOS to each factor can be fully explained by the long-term average precipitation gradient across the Tibetan Plateau ( Figure 6). Spatially, the SOS of grassland in the Tibetan Plateau was more strongly affected by mean temperature in the wetter areas (I and II ecoregions), while it was dominated by the accumulated precipitation during spring in drier areas (IV ecoregion). Considering the values averaged from the pixels with significant individual contribution at p < 0.05 level in different ecoregions in the Tibetan Plateau, a 10 mm increase in long-term average precipitation in spring responded to an increase in the individual contribution of T spring of 1.70% (p < 0.05), while it caused a decrease in an individual contribution of P spring of −0.70% (p < 0.1) (Figure 6a). In contrast, among different ecoregions on the Tibetan Plateau, this pattern linking the relative contribution of each climate factor to long-term average mean temperature in spring was not found in our study (Figure 6b). Furthermore, this similar pattern of the impacts of long-term average cumulative precipitation on the contribution of each climate factor to the SOS among the different vegetation types was found (Figure 7). Specifically, a 10 mm increase in long-term average precipitation in spring resulted in an increase in the individual contribution of Tspring of 2.0% (p < 0.05) to the SOS of meadow, while it caused a decrease in the individual contribution of Pspring of −0.30% (p < 0.1) (Figure 7a). Regarding steppe across the different ecoregions of the Tibetan plateau, these similar spatial variations were also found, with the sensitivity of individual contribution of Tspring and Pspring to the long-term average precipitation in spring of 0.07% mm −1 and −0.14% mm −1 , respectively (Figure 7c). Variations in individual contribution of each climate factor for grassland spring phenology to long-term average cumulative precipitation (a) and mean temperature (b) in spring in different ecoregions across Tibetan Plateau. Bars indicate mean long-term average precipitation or temperature in spring for each ecoregion. Point and dashed curves represent averaged value of pixels with significant individual contribution at p < 0.05 level for each ecoregion listed in right y-axis labels. Solid line represents linear regression of individual contribution of each climate factor to long-term average precipitation or temperature. p value denotes significance. NoSig indicates that slope is not significant at p ≥ 0.05 level. I, II, III, IV, and V indicate different ecoregions following Figure 1.
Furthermore, this similar pattern of the impacts of long-term average cumulative precipitation on the contribution of each climate factor to the SOS among the different vegetation types was found (Figure 7). Specifically, a 10 mm increase in long-term average precipitation in spring resulted in an increase in the individual contribution of T spring of 2.0% (p < 0.05) to the SOS of meadow, while it caused a decrease in the individual contribution of P spring of −0.30% (p < 0.1) (Figure 7a). Regarding steppe across the different ecoregions of the Tibetan plateau, these similar spatial variations were also found, with the sensitivity of individual contribution of T spring and P spring to the long-term average precipitation in spring of 0.07% mm −1 and −0.14% mm −1 , respectively (Figure 7c).  ) and steppe (c,d) to long-term average cumulative precipitation and mean temperature in spring in different ecoregions across Tibetan Plateau. Each point shows averaged values of pixels with significant individual contribution at p < 0.05 level for meadow and steppe in different ecoregion. Error bar is SEM (standard error of mean). Line represents linear regression, and shaded area represents 95% confidence interval. p value denotes significance. NoSig indicates that slope is not significant (p ≥ 0.05).

Interacting Effects of Temperature and Precipitation on Spring Phenology
Variability, particularly variability in temperature and precipitation, unequivocally affected the vegetation phenology on the Tibetan Plateau [5,9,12,15,27,28]. However, the relative contribution of each factor and the underlying combined effect on vegetation phenology received much less attention. In this study, we found that the individual contribution of temperature and precipitation to spring phenology was lower, despite the higher correlation coefficient consistently shown by most previous studies. More importantly, individual contribution of temperature and precipitation to interannual variations in SOS revealed a distinct east west disparity across the Tibetan Plateau, which can be fully explained by long-term average precipitation gradient over the entire study period. This spatial pattern of the quantitative contribution of temperature and precipitation to interannual variations in spring phenology further confirmed the stronger interactions of precipitation and temperature and indicated underlying mechanisms of the responses of ecological functions of vegetation to climate change.
With increasing temperature and precipitation [28,48], all methods in our study agreed on an advanced trend in SOS for meadow and steppe on the Tibetan Plateau during 2000-2017 despite the spatial heterogeneity, with an advanced SOS across most of the Plateau and a discretely distributed and delayed SOS in the midwestern region of the Tibetan Plateau. These complex trends in grassland spring phenology are generally supported by recent studies [28,30,[49][50][51]. Most previous studies consistently showed that the  ) and steppe (c,d) to long-term average cumulative precipitation and mean temperature in spring in different ecoregions across Tibetan Plateau. Each point shows averaged values of pixels with significant individual contribution at p < 0.05 level for meadow and steppe in different ecoregion. Error bar is SEM (standard error of mean). Line represents linear regression, and shaded area represents 95% confidence interval. p value denotes significance. NoSig indicates that slope is not significant (p ≥ 0.05).

Interacting Effects of Temperature and Precipitation on Spring Phenology
Variability, particularly variability in temperature and precipitation, unequivocally affected the vegetation phenology on the Tibetan Plateau [5,9,12,15,27,28]. However, the relative contribution of each factor and the underlying combined effect on vegetation phenology received much less attention. In this study, we found that the individual contribution of temperature and precipitation to spring phenology was lower, despite the higher correlation coefficient consistently shown by most previous studies. More importantly, individual contribution of temperature and precipitation to interannual variations in SOS revealed a distinct east west disparity across the Tibetan Plateau, which can be fully explained by long-term average precipitation gradient over the entire study period. This spatial pattern of the quantitative contribution of temperature and precipitation to interannual variations in spring phenology further confirmed the stronger interactions of precipitation and temperature and indicated underlying mechanisms of the responses of ecological functions of vegetation to climate change.
With increasing temperature and precipitation [28,48], all methods in our study agreed on an advanced trend in SOS for meadow and steppe on the Tibetan Plateau during 2000-2017 despite the spatial heterogeneity, with an advanced SOS across most of the Plateau and a discretely distributed and delayed SOS in the midwestern region of the Tibetan Plateau. These complex trends in grassland spring phenology are generally supported by recent studies [28,30,[49][50][51]. Most previous studies consistently showed that the increasing spring temperature was the major dominant driver of SOS advances [11,15,17,27,28,52]. However, the individual contribution of temperature and precipitation was lower, despite a higher correlation coefficient, which implied that these factors interact to affect the SOS across the Tibetan Plateau. The contributions of temperature and precipitation to the grassland spring phenology varied across space on the Tibetan Plateau, and these spatial heterogeneities can be mainly explained by the spatial gradient of long-term average precipitation over 2000-2017.
For meadow, the interannual variations in SOS were dominated by T spring in the eastern I ecoregion with relatively mild climate. In these wetter areas, the soil moisture constraints were released. Without limitation of water resources, vegetation tended to maximize the thermal benefit and showed higher sensitivity to temperature than to precipitation [12,28,53]. In other words, better thermal-hydraulic conditions would accelerate the rates of chemical reactions because of the effect on Rubisco enzymatic activity and then speed up development processes in photosynthetic organisms [54,55]. Consequently, the individual contribution of temperature to SOS would be significant and largest, that is, the increasing preseason temperature would significantly stimulate plant growth and advance the SOS under this better thermal-hydraulic conditions. Similarly, Shen et al. [12] showed that the sensitivity of SOS to temperature was negative and significant in eastern and northeastern part of the Plateau, and Piao et al. [28] showed that the SOS could advance by an extra 4.5 days with an increase in spring temperature of 1 • C. In contrast, the interannual variations in SOS were strongly negatively correlated with P spring in the western IV ecoregion with less precipitation. In those relatively dry ecoregions with spring precipitation of less than 100 mm [50], vegetation growth initiation after winters with low rainfall would be limited by water availability [56]. Limited water potential would inhibit plant growth and photosynthesis activities, increase the risk of chlorophyll degradation and plant mortality [57], and consequently delay the SOS. In addition, the soil moisture would be suboptimal because of the high evapotranspiration caused by increasing temperature [58]. Hence, the preseason precipitation would determine the water availability, and revealed a high individual contribution to the interannual variation in SOS.
Previous studies showed that meadow and steppe responded differently to climate change [50,53,59]. The different habitat conditions were primarily responsible for these various responsive characteristics. The steppe adapted to the long period of colder and drier weather with a shorter growth cycle. Thus, the steppe's ecological functions showed increased sensitivity to climate change [39,40,[60][61][62][63]. Moreover, although the warming accelerates plant growth, it can also cause the decline of soil moisture, and then increase the sensitivity of vegetation phenology to precipitation [5,27,53], which indicated that warming and precipitation would additionally interact to affect plant spring phenology in the drier area. Hence, changes and interaction between temperature and precipitation would affect the SOS for steppe significantly. Nonetheless, the long-term average precipitation gradient also dominated the spatial disparity of the individual contribution of precipitation and temperature; that is, the vegetation SOS would be less sensitive to temperature changes due to the lower long-term average precipitation, while the relative contribution of preseason precipitation would increase. These dynamic response patterns can be well explained by the vegetation adaptation strategy, which would maximize the benefit from the limiting climate factors and reduce the risk imposed by other factors [12].
Generally, precipitation and temperature were considered as two main regulators of vegetation activity on the Tibetan Plateau [24]. Their underlying combined effect on vegetation phenology suggested a vegetation adaptation strategy that achieves a balance between maximizing benefit and minimizing risk; that is, the vegetation would make the most of limiting climate factors and reduce the risks from other factors [12]. In the northeastern Tibetan Plateau, the vegetation spring phenology was strongly advanced by increasing temperature and slightly advanced by precipitation. The plant photosynthetic rate would be directly accelerated due to the effect of Rubisco enzymatic activity in these better thermal-hydraulic conditions [54]. In contrast, in the drier area, the higher spring temperature would cause the decline of soil moisture, and then increase the sensitivity of vegetation phenology to precipitation [5,27,53], which indicated that warming and precipitation would additionally interact to affect plant spring phenology. Alternatively, the higher individual contribution of preseason precipitation to SOS in more arid areas indicated that plant would intensify spring drought because bulk canopy water needs to be increased with an advanced SOS under future climate warming. Overall, our analysis highlighted stronger interactions between warming and precipitation and quantified the relative contribution of each factor, which would benefit to gain an accurate mechanism and a better understanding of terrestrial ecosystem processes and their responses to future climate change.

Uncertainties and Further Studies
To account for the potential interactions and the individual contribution of each factor, we employed two different methods to detect SOS and conducted redundancy analysis. Although a similar spatial pattern of the SOS trend was revealed by these two methods, distinct values in the interannual variations still existed. The proportion of delayed trend of SOS from Polyfit was 25.74% (Figure 3a), which was slightly less than that from HANTS, which was 28.19% (Figure 3b). This distinction maybe caused by the spurious oscillations of the NDVI time series conducted in the HANTS and Polyfit methods [8,[64][65][66][67]. Therefore, the results need to be further tested with more models to characterize more exact interannual variation in the SOS. In addition, previous studies showed that many other factors would also affect the vegetation phenology, for example, large diurnal temperature range [14,17,68,69], soil moisture [24,[70][71][72], and some vegetation functions [73][74][75]. In addition, some studies also showed that the photoperiod could coregulate the vegetation spring phenology through its interaction with temperature or its influence on the stomata aperture and photosynthetic active radiation [22,23,56]. However, Fan et al. [76] showed that sunshine duration in the mountain plateau zone displayed an insignificantly negative trend (p > 0.05) at the rates of −0.001 h year −1 during 1986-2015. Thus, the photoperiod was roughly assumed to be a valid constant without the effect of cloudiness during the study period. Our study also found that T winter was the most significant factor for the interannual variation of the steppe in the central Plateau, which could be attributed to warminginduced changes in soil moisture. Moreover, the vegetation would require adequate chilling conditions (vernalization) during endodormancy; therefore, grassland phenology is expected to be sensitive to winter warming [13,21]. Clearly, further analysis with more factors should be conducted to estimate more robust contribution of each factor, and then support these inferences and their role in the control over phenology. Nevertheless, our present work highlighted a stronger interaction between warming and precipitation, quantified the relative contribution of each factor, and identified temporal-spatial aspects of grassland SOS responses to critical climate factors on the Tibetan Plateau, which would provide a helpful reference and establish a better understanding for further studies on climate plant interactions.

Conclusions
This study quantified the individual contribution of warming and precipitation for the interannual variations in start of growing season (SOS) in the world's largest cold and arid/semiarid regions. Our results further confirmed the strong impacts of preseason precipitation on satellite-derived estimates of spring phenology of grassland across the Tibetan Plateau. First, the relative contribution of the accumulated precipitation in spring to interannual variations in SOS was higher in more arid than wetter areas. Alternatively, P spring (the accumulated precipitation in spring) strongly dominated vegetation SOS over most pixels located in the western part of the Plateau. Second, although the responses of SOS to each factor were complex and fragmented across the Tibetan Plateau, these spatial heterogeneities can be mainly explained by the long-term average precipitation during spring. Specifically, the increase in long-term average precipitation during spring resulted in an increase in the individual contribution of T spring to the SOS of grassland, while the individual contribution of the accumulated precipitation in spring would increase with the decrease in long-term average precipitation during spring. Alternatively, climate warming would have less impact on SOS than precipitation, which would lead to complex responses to climates change in these arid or semiarid regions. In addition, temperature in winter made a larger contribution, which would change soil moisture and strongly affect the SOS in the central Plateau. For the Tibetan plateau dominated by the arid/semiarid climate, considering the interactions of preseason temperature and precipitation and their stronger and more complex effect on SOS, vegetation would experience greater SOS advancement, with substantial warming in relatively wetter regions than in drier areas. Instead, substantial warming would slightly advance vegetation SOS and might delay SOS in the arid/semiarid regions because of the increase in evapotranspiration. Thus, the combined impacts and the quantitative contribution of warming and precipitation on SOS should be considered, while assessing the responses of vegetation to climates change rather than the unitary effect of one or more factors directly.