Specific Drivers and Responses to Land Surface Phenology of Different Vegetation Types in the Qinling Mountains, Central China

Abstract

Land surface phenology (LSP), as a precise bio-indicator that responds to climate change, has received much attention in fields concerned with climate change and ecology. Yet, the dynamics of LSP changes in the Qinling Mountains (QMs)—A transition zone between warm-temperate and north subtropical climates with complex vegetation structure—under significant climatic environmental evolution are unclear. Here, we analyzed the spatiotemporal dynamics of LSP for different vegetation types in the QMs from 2001 to 2019 and quantified the degree of influence of meteorological factors (temperature, precipitation, and shortwave radiation), and soil (temperature and moisture), and biological factors (maximum of NDVI and middle date during the growing season) on LSP changes using random forest models. The results show that there is an advanced trend (0.15 days/year) for the start of the growing season (SOS), a delayed trend (0.24 days/year) for the end of the growing season (EOS), and an overall extended trend (0.39 days/year) for the length of the growing season (LOS) in the QMs over the past two decades. Advanced SOS and delayed EOS were the dominant patterns leading to a lengthened vegetation growing season, followed by a joint delay of SOS and EOS, and the latter was particularly common in shrub and evergreen broadleaved forests. The growth season length increased significantly in western QMs. Furthermore, we confirmed that meteorological factors are the main factors affecting the interannual variations in SOS and EOS, especially the meteorological factor of preseason mean shortwave radiation (SWP). The grass and crop are most influenced by SWP. The soil condition has, overall, a minor influence the regional LSP. This study highlighted the specificity of different vegetation growth in the QMs under warming, which should be considered in the accurate prediction of vegetation growth in the future.

1. Introduction

Vegetation phenology is the seasonal timing of lifecycle events, such as leaf emergence, flowering, leaf coloration and fall, and it has become an important topic in the field of climate and ecology as a sensitive and precise indicator that is responsive to climate warming [1,2]. Shifts in spring and autumn vegetation phenology caused by climate warming can differentially alter the length of the growing season, which affects carbon, water, and energy exchange between terrestrial ecosystems and the atmosphere [3,4,5]. Recent studies have reported that in addition to climatic factors, soil and biological factors also influence shifts in vegetation phenology by affecting plant growth processes in the context of ongoing global climate change [6,7], due to the poor interpretation of phenology shifts among different vegetation types [8,9]. Hence, it is essential to study the dynamics and drivers of phenology among different vegetation types to improve phenology models and enrich our understanding of the carbon cycle of terrestrial ecosystem.

With the application of remote sensing in monitoring vegetation phenology, we traditionally use the term land surface phenology (LSP) to denote the dynamic variations in vegetation land surface as observed from satellite imagery [10]. Satellite-derived LSP metrics are usually focused on the start (SOS) and end (EOS) of growing seasons [11]. Satellite-based studies have shown that SOS was advanced by 10.6 days (i.e., 5.4 days per decade) throughout Europe and by 14 days (i.e., 7.9 days per decade) in temperate China before 2000 [12,13]. However, this trend of SOS advancement may have slowed or even reversed since the 2000s. For instance, in the entire northern hemisphere, it was advanced by only 0.2 days during 2000 to 2008, but a delayed SOS was revealed in the Tibetan Plateau [14,15]. Regarding satellite-derived EOS, the published results have not always been consistent. Across the entire Northern Hemisphere, EOS was delayed at a rate of 2.2 days per decade during 2000–2008 [14]. In the Yellow River Basin, EOS was delayed by 5.6 and 3.4 days in 1982–1999 and 2000–2015, respectively [16]. In the Qinghai-Tibet Plateau, however, Wang et al. [17] reported the opposite phenology change trends, in both the east and west zones. Overall, these varied results might be due to different study areas, periods, and methods of extracting phenology metrics. However, few studies have focused on the diversity of phenology across different vegetation types. In particular, the dynamic phenology characteristics of herbaceous or shrubs and evergreen forests in the Qinling Mountains have been little studied.

To date, the processes and drivers governing LSP remain poorly understood. Several studies reported that temperature is a major driver of early spring leaf development and delayed autumn leaf fall in plants and has less control over autumn phenology than spring phenology [18,19]. The impact of precipitation on LSP processes is mainly directed at plants in arid and semiarid regions, where water deficits limit the use of light and heat conditions by plants in arid and semiarid areas [20]. Some studies further considered solar radiation (i.e., shortwave radiation) and found that increasing photosynthetic active radiation can promote earlier leaf germination and delaying leaf senescence [21]. Besides meteorological factors, soil factors and biological factors have also been shown to be drivers of LSP change processes [22,23,24]. Soil temperature and moisture information, due to the prevalent freezing—Thawing process of soil in alpine and arctic regions, could more directly control vegetation growth; for example, soil wetting will, to some extent, reduce the effect of soil warming on LSP changes [22,25]. Slight fluctuations in the time interval between the middle date (MD) of the growing season and the autumn phenology have a strong effect on regulating EOS [7]. Peak growth in summer (i.e., maximum NDVI during the growing season, MN) can have an impact on vegetation greening and senescence, and its unique vegetation growth patterns may result in different allocations of green or carbon across the growing season [24]. To date, how these meteorological, soil, and biological factors affect regional-scale LSP variations has not been clearly and consistently studied, which seriously affects our ability to predict LSP periods.

The responses of vegetation growth processes to climate change are inherently nonlinear [26]. Random forest (RF), as a nonparametric multivariate method, can explain nonlinear processes to a large extent [27]. The advantage of the RF model is that it can consider many predictor variables and nonlinearly determine the relative importance of each predictor variable [28]. Due to its high efficiency in handling the potentially complex relationships between LSP periods and meteorological, soil, and biological factors, RF has been widely and successful applied in recent ecology studies [6,29].

The Qinling Mountains (QMs), rich in vegetation types, represent a demarcation line of climate in China and are also an area characterized by sensitivity or response to climate change, with a significant upward temperature trend in the past half century [30]. Here, we used the satellite derived normalized difference vegetation index (NDVI) records (2001–2019) from MOD13A2 to extract the LSP dates of QMs. The objectives were (a) to explore the temporal and spatial trends of LSP in the QMs, (b) to quantify the relative contribution of SOS and EOS of different vegetation types to the length of growing season (LOS) and to determine the dominant growth pattern during the growing season, and (c) to simulate the LSP dates and assess the relative importance of meteorological, soil, and biological factors on the interannual variations in LSP. This study focuses on the specificity of different vegetation growths in the QMs, and the results are helpful for future accurate prediction of vegetation growth and to develop scientific management strategies.

2. Materials and Methods

2.1. Study Area

Qinling is the highest mountain range in the central region of China and also a geographical boundary between the north and the south of China, with an elevation range of 51 to 5120 m and a spatial range from 30.8° to 35.5°N to 102.5° to 114.6°E (Figure 1a). Its climate differs significantly between north and south, with a humid northern subtropical climate in the south and a warm temperate semi-humid and semiarid climate in the north (Figure 1a). It has also been classified as one of the critical terrestrial biodiversity areas of world significance [31], with mixed coniferous and deciduous broadleaved forests widely distributed on its northern slopes, while the southern slopes are dominated by mixed evergreen and deciduous broadleaved forests (Figure 1b).

2.2. Datasets

NDVI is the most commonly applied vegetation index to characterize vegetation greenness and is strongly correlated with vegetation photosynthetic activity [32]. Climate change is the main factor affecting the change of vegetation greenness, and this change can be reflected by the spectral information of NDVI images. In this study, we used the NDVI datasets generated from NOAA/AVHRR series satellite images by the NASA MODIS13A2 group (Table 1). We used this dataset to extract the LSP dates for QMs from 2001 to 2019. We also excluded areas of bare soil/sparse vegetation with an annual average NDVI of less than 0.1 [33].

Meteorological data, including daily mean air temperature, daily total precipitation, and daily mean shortwave radiation and soil data, including daily mean soil temperature and moisture in 0–100 cm soil layer, from 2001 to 2019 were used in this study. These gridded data were derived from the ERA5-Land hourly data (Table 1). Moreover, we transformed the hourly climate data (24 hourly data were averaged for temperature, solar radiation, soil temperature and humidity, 24 hourly data were summed for precipitation) to daily-scale temporal resolution and resampled meteorological data to the same resolution as MODIS13A2 data. A time lag of 30 days before SOS and 60 days before EOS is defined as preseason.

The 300 m spatial resolution Climate Change Initiative Land Cover (CCI-LC) maps from 2001 to 2019 were available from the European Space Agency (ESA) (Table 1). CCI-LC discriminates 22 classes of land cover. In this study, we resampled these maps to 1 km and analyzed only pixels of unchanged vegetation types containing evergreen needle leaved forest (ENF), evergreen broadleaved forest (EBF), deciduous broadleaved forest (DBF), mixed forest (MF), shrubland (SL), grassland (GL), and cropland (CL).

2.3. Retrieval of Phenology Metrics from NDVI Time Series Data

The premise of quantitatively analyzing the phenology changes is to derive several key phenology metrics: SOS, EOS, LOS, MN, and MD (Figure 2). In this study, we firstly stacked the NDVI images from 2001 to 2019 in chronological order and smoothed the NDVI time series with a Savitzky-Golay (SG) filter for each pixel per year. The SG filter was chosen because it can best preserve the temporal vegetation dynamics and minimize atmospheric contamination and has also been integrated into the processing of the MODIS phenology product [34]. The smoothed data was used further for extracting phenology metrics of different vegetation types by detecting the inflection point (i.e., date) when the NDVI time series begins to ascend or descend for the specific year. This is the derivative method which the phenology metrics were extracted for each pixel per year, whereby the maximum value of NDVI_ratio corresponds to the greatest change of the smoothed NDVI time series [35]. Equation (1) is given as

(1)${\mathrm{NDVI}}_{ratio(t)}=\frac{{\mathrm{NDVI}}_{t+1}-{\mathrm{NDVI}}_t}{{\mathrm{NDVI}}_t}$

where ${\mathrm{NDVI}}_{ratio\left(t\right)}$ is the calculated relative changing rate of NDVI at time t and ${\mathrm{NDVI}}_t$ is the NDVI value at time t. Occurrence dates were obtained using these smoothed NDVI time series. SOS and EOS dates were determined as the day with the maximum and minimum ${\mathrm{NDVI}}_{ratio}$. LOS was determined to be the difference between EOS and SOS. MN was defined as the peak of vegetation growth, i.e., the NDVI value corresponding to ${\mathrm{NDVI}}_{ratio}$ closest to zero. MD date was the middle date between the EOS date and the SOS date. The description of phenology metrics correlations is shown in Figure 2.

2.4. Method and Statistical Analysis

2.4.1. Trend Analysis

The method used in this study is shown in Figure 3. The spatiotemporal trends of LSP during 2001–2019 were estimated using Sen’s slope method, also known as the Theil-Sen median method. The method is a robust nonparametric statistical method for trend calculation that is insensitive to measurement errors and is far more accurate than nonrobust simple linear regression [36]. Sen’s slope was calculated using Equation (2):

(2)$Sen’sslope=Median\left(\frac{x_j-x_i}{j-i}\right)$

where the median is the mean value of all the slopes, and x_i and x_j represent the LSP dates of years i and j. A negative Sen’s slope indicates an advancing trend, whereas a positive Sen’s slope indicates a delaying trend.

Then, we used the Mann-Kendall (MK) method to test the significance of time series trends, which is a nonparametric statistical test and is robust to outliers [37]. We used the normal cumulative distribution function to determine the p-value of the MK test statistic with a significant confidence level of p < 0.1. In this study, we used Sen’s slope and MK test to trend analysis and significance test the spatial distribution (average growth) and interannual variation (interannual growth) of LSP for different vegetation types from 2001 to 2019, all using MATLAB 2017a were completed.

2.4.2. Change Pattern and Relative Attribution Analysis

To further understand the seasonal changes of vegetation growth in the QMs in the past two decades, we divided the trend changes of LSP (SOS, EOS, and LOS) into six combinations, where each combination represents a pattern of vegetation growth (Table 2). We used spatial analysis to count the proportion and significance of different vegetation for each pattern and to derive the dominant pattern of seasonal changes in the growth of each vegetation. All analysis was accomplished using ArcGIS 10.4.

To evaluate the symmetry of SOS and EOS for LOS changes, we used the C-index proposed by Garonna et al. [38] to calculate the relative contribution of trends in SOS and EOS for LOS changes. It was calculated as follows:

(3)$C=-1+\frac{2\cdot \left|SO{S}_{slope}\right|}{\left|SO{S}_{slope}\right|+\left|EO{S}_{slope}\right|}$

where SOS_slope and EOS_slope are the Sen’s slope of SOS and EOS, respectively. A positive C value indicates that the trend in LOS is mainly attributable to changes in EOS, and a negative C value indicates that the trend in LOS is mainly attributable to changes in SOS. The variation of C value is from −1 to 1.

2.4.3. Analysis of the Relative Importance of Different Drivers

Based on previous studies on the drivers of interannual variation in vegetation phenology [18,19,20,21,22,23,24], we selected drivers such as Table 3 to simulate SOS and EOS. These drivers are divided into three main categories: meteorological, soil and biological factors, which are further divided into preseason cumulative values and cumulative values throughout the growing season.

In this study, the RF model was used to assess the relative importance of the drivers affecting interannual variations in LSP. First, for each pixel, we calculated the partial correlation coefficients between environmental factors (TP, TG, PP, PG, SWP, SWG, STP, STG, SMP, SMG, MN, and MD) and SOS and EOS during 0, 1, 2, … n months before SOS and EOS, and we separately derived the highest correlation with SOS and EOS for the time range of 30 days before SOS and 60 days before EOS for environmental data. Moreover, a subset of variables highly correlated with SOS and EOS was selected, and the values of the variables at selected years and locations (spatiotemporal models) were combined into a set of input feature vectors that are used as inputs for the RF algorithm. Then, these feature vectors were divided equally into two subsets, with 2/3 of the dataset used for model training (in bag) and the remaining 1/3 of the dataset used as an additional test of the RF internal computation (out of bag, i.e., OOB) to estimate the importance of each variable. Variable importance can also be measured by OOB, which compares the increases in OOB error with that variable randomly permuted and all others unchanged [39]. The importance score of a variable is as follows:

(4)$VI\left(X^j\right)=\frac{1}{ntree}{\displaystyle \sum _t\left(err’OO{B}_t^j-errOO{B}_t^j\right)}$

where $X^j$ is the jth variable, ntree is the number of trees, $errOO{B}_t^j$ is the OOB error of each tree t, and $err’OO{B}_t^j$ is the OOB error when $X^j$ is permuted, while all other variables remain unchanged among OOB data. For regression, the OOB error is the mean square error.

Finally, to optimize the model, the hyperparameter search was used to select the best tested hyperparameter set, and the optimal model was trained on the whole training set for accurate prediction of the LSP dates. Besides modeling multi-year-scale LSP variation for the entire region, subregional variation according to different vegetation types was also modeled to quantify the relative importance of different drivers. These models were constructed using the RF package in R statistical software.

To evaluate the predictiveness of the model and to further test the applicability of the model, we used a randomly selected subset (1/3 of the dataset) for model validation. Both the proportion of explained variance (R²) and the root mean square error (RMSE) were used to assess the performance of the model on the complete datasets. The stability of the model fit is explained by the mean absolute error (MAE) [40]. These statistics were calculated as follows:

(5)$R^2=1-\left[\frac{{\displaystyle \sum _{i=1}^n{(y_i-{\stackrel{⌢}{y}}_i)}^2}}{{\displaystyle \sum _{i=1}^n{(y_i-{\overline{y}}_i)}^2}}\right]$

(6)$RMSE=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{(y_i-{\stackrel{⌢}{y}}_i)}^2}}{n-1}}$

(7)$MAE={\displaystyle \sum _{i=1}^n\left|\frac{y_i-{\stackrel{⌢}{y}}_i}{n}\right|}$

where $y_i$ is the observed satellite-based SOS or EOS values, is ${\stackrel{⌢}{y}}_i$ the predicted RF-based SOS or EOS values, and ${\overline{y}}_i$ is the mean satellite-observed SOS or EOS values of all selected test pixels for 2001–2019. n is the sum of all selected test pixels.

3. Results

3.1. Spatiotemporal Variations of Phenology Metrics in the QMs

Figure 4 shows the spatial variation of the annual mean LSP and their corresponding standard deviations (Std) over the study period 2001–2019. Earlier (<90 days) sites of SOS (14.4%) were located at low elevations in the central QMs, and later (>130 days) sites (6.8%) were located at high elevations in the western QMs (Figure 4a). The earliest occurrence of SOS was for SL, with a mean SOS of 97 ± 14 days, and the latest occurrence was for GL, with a mean SOS of 109 ± 12 days (Figure 4a). The Std of SOS has significant spatial variation, with larger areas located in the southwestern QMs (57.2%) having higher Std (>9 days) (Figure 4b). The overall spatial variation in multi-year average EOS was not significant, and in the southeastern QMs, EOS was mainly concentrated in 290–300 days, accounting for 43.0% of the entire study area (Figure 4c). The earliest EOS occurred in GL, with a mean EOS of 287 ± 11 days, and the latest occurred in EBF, with a mean EOS of 295 ± 12 days (Figure 4c). There was higher Std (>13 days) in the southern QMs (25.2%) compared with the northern QMs (Figure 4d). LOS had clear spatial differences, with the central QMs (15.5%) having the longest LOS (>210 days) and the western high-altitude areas (8.5%) had the shortest LOS (<150 days). For SL in the southern QMs (27.9%), its mean LOS was 193 ± 20 days and the std (>18 days) was also the largest (Figure 4e,f).

We also characterized the spatial distribution of LSP trends for different vegetation types from 2001 to 2019 (Figure 5). For the whole QMs, SOS was advanced in 67.8%, the average rate of advance was 1.5 days/decade, and 27.5% of the area (mostly located in the northern QMs) was significant (Figure 5a,b). DBF advanced at a rate of 1.9 days/decade and was the fastest compared to other vegetation types (Figure 5a). EOS was delayed in 72.1% of the region and significant for 42.1% of the region (mostly located in the southern QMs), with an average delay rate of 2.4 days/decade across the region (Figure 5c,d). EBF had the fastest delay rate of 3.3 days/decade (Figure 5c). The average rate of LOS lengthening across the study area was 3.9 days/decade, and 74.6% of the areas (mostly in the southwestern QMs) were lengthened (Figure 5e). The rate of LOS lengthening was 4.7 days/decade for EBF, fastest among the seven vegetation types. Of these areas that were lengthened, the change was found to be significant in 40.3% (mostly in the western QMs) of cases (Figure 5f).

Data of interannual variation trends and the significance of LSP for different vegetation types are shown in Figure 6. Overall, the Sen’s slope of SOS is −0.09 days/year from 2001 to 2019, but this advance is insignificant. There was a trend of significant SOS advancement for DBF and GL, at 0.16 and 0.13 days/year, respectively. EOS shows a significant delay trend over the entire region of 0.29 days/year. The trend of EOS delay was more significant for both EBF and CL compared to other vegetation (p < 0.05), and the rate of EBF delay was the fastest (0.37 days/year). LOS is significantly lengthened at a rate of 0.48 days/year. ENF, EBF, DBF, MF, and CL show a more significant trend for lengthened LOS (p < 0.05), with the fastest being for EBF (0.68 days/year).

3.2. Change Pattern of LSP and Relative Attribution Analysis

For all vegetation types, the dominant change pattern of the growing season was Type I, with a proportion of 48.4%, which implies that most vegetation in the QMs had an advanced SOS, delayed EOS, and lengthened LOS (Table 4). Type V showed the second largest proportion (15.2%) which meant that there were also many plant species on the QMs having delayed SOS, delayed EOS, and lengthened LOS. Types III, IV, and VI had the smallest proportions (all lower than 10.0%), which indicated that the probabilities of shortened LOS were very low for all plants on the QMs. For five of these vegetation types (ENF, EBF, MF, SL, and CL), the dominant change pattern of the growing season was Type I, followed by Type V. This indicates that these types of plants on the QMs had delayed EOS and lengthened LOS. The main change pattern for DBF and GL was also Type I, with Type II being the second most prevalent. This implies that there is some DBF and GL showing advanced SOS, advanced EOS, and lengthened LOS.

Figure 7 shows the significance of each pattern. For all vegetation types, lengthened LOS, advanced SOS, and delayed EOS were significant in terms of Type I, II, and V change patterns, respectively. Types I, II, and V were the top three patterns in terms of percentage, as seen in Table 4, which were also the three patterns of LOS lengthening. In the Type I, lengthened LOS is significant for most vegetation types (EBF, MF, SL, GL, and CL). In terms of Type II change, the SOS of ENF, DBF, MF, GL, and CL were all significantly advanced, and advanced SOS resulted in lengthened LOS. Fewer changes in LSP trends were significant in terms of Type III, IV, and VI changes, only SL and GL showed significant Type IV changes, and SL (15.5%) and GL (12.2%) accounted for a large proportion of Type IV changes. In terms of Type V, the EOS of all six vegetation types (ENF, EBF, DBF, MF, GL, and CL) was significantly delayed. Delayed EOS resulted in lengthened LOS, a pattern which also happens to be the second largest in terms of proportion, and this pattern is also the one we should be concerned about.

For the entire study area, the calculated C-index values were negative for 106,385 pixels and positive for 112,514 pixels (Table 5 and Figure S1). Pixels with C-index values less than 0 are mainly distributed in the easternmost and southernmost parts of the QMs, which indicates that their LOS variations are mainly controlled by SOS shifts. All other regions have pixels with values over 0, which indicates that they are controlled by EOS shifts (Figure S1, in Supplementary Materials). The percentage of LOS changes controlled by SOS and EOS is 48.6% and 51.4%, respectively (Table 5). This also shows that LOS trends, except for DBF and GL, were mainly controlled by the shift in EOS for each vegetation type. The percentage of LOS changes controlled by SOS shifts was 53.4% and 52.0% for DBF and GL, respectively. The largest percentage of changes in LOS of all vegetation types controlled by EOS was 58.3% (SL), and the smallest was 46.6% (DBF).

3.3. Drivers of Interannual Variations in LSP

For the QMs, different drivers affect the interannual variability in SOS and EOS (Figure 8 and Table 6). The SWP, MD, and STP are the three most important factors influencing the interannual SOS variation, and the relative importance accounts for 54.4% of total (Figure 8a and Table 6). The total percentage of TG, PP, TP, STG, and PG was 39.2%, and the effect of TG and PP on the interannual SOS variation was almost the same. The remaining four variables (SWG, SMG, SMP, and MN) have a very small effect on interannual SOS variation, and their combined percentage was only 6.4%. Figure 8b and Table 6 also show that SWP, PP, and MD are the three most important factors influencing the interannual EOS variation, with a total relative importance of 54.0%, and the influence of SWP is much stronger than that of PP and MD. The effects of TP, SWG, STP, PG, TG, and STG on the interannual EOS variation totaled 41.9%, and the effects of TP and SWG on the interannual EOS variation were not very different, with a relative importance of 10.2% and 9.9%, respectively. There was also little difference in the relative importance of PG, TG, and STG. The effect of the single variable of STP on the interannual EOS variation (7.2%) is much larger than the sum of SMP, SMG, and MN (4.1%).

The drivers influencing interannual variations in LSP of different vegetation types were assessed (Figure 9 and Table 6). Figure 9a and Table 6 show that the main drivers affecting the interannual SOS variation of ENF, MF, and GL were SWP, MD, and STP, and the relative importance of SWP in these three vegetation types was ranked as GL (28.6%) > MF (22.1%) > ENF (21.2%). The interannual SOS variation of EBF, DBF, and CL was mainly influenced by MD and SWP. The main drivers influencing the interannual SOS variation of SL were STP, TP, and SWP, and STP (29.2%) was the most important factor influencing the interannual SOS variation of SL. As shown in Figure 9b and Table 6, in terms of the interannual EOS variation, SWP had the strongest effect on GL (36.4%) and the slightest effect on SL (20.5%). The main drivers of interannual EOS variation of ENF, DBF, and CL are SWP, TP, and PP. SWP, PP and MD are the main drivers of interannual variations in EOS for EBF, MF and SL. Besides SWP, which is the most important driver, PP is the second factor affecting MF and SL, and MD is the second factor affecting EBF with a relative importance of 14.5%. Moreover, the main factors affecting the interannual EOS variation of GL are SWP, SWG, and TP, and the relative importance of SWG and TP is 11.5% and 11.2%, respectively.

We used 19 years of data to construct an RF model for seven vegetation types (Table S2), and we randomly sampled 1/3 of the pixel dataset to assess their linear relationships (Figure S2). The actual and predicted SOS displayed good linear relationships, with their correlation coefficient R2 values ranging from 0.900 to 0.938, RMSE values ranging from 4.90 to 6.52, and MAE values ranging from 3.93 to 5.15 (Figure S2a). Both the actual and predicted EOS also show good linear relationships, with their correlation coefficients R2 values ranging from 0.911 to 0.942, RMSE values ranging from 4.23 to 5.63, and MAE values ranging from 3.26 to 3.74 (Figure S2b). These results indicate that it is appropriate to use RF models to analyze interannual variations in LSP in the QMs.

4. Discussion

4.1. Dynamics Changes in LSP in the QMs

Understanding the interannual variations in vegetation phenology and its trends is important for recognizing the patterns of vegetation growth dynamics as a response to climate warming. Our study showed that there is an advanced trend (1.5 days/decade) for SOS, a delayed trend (2.4 days/decade) for EOS, and an overall extended trend (3.9 days/decade) for LOS in the QMs during 2001–2019. In comparison with previous studies on phenology changes, there were different degrees of advanced SOS, delayed EOS, and lengthened LOS in different study areas and periods [32,41]. For example, during 1982–2006, the SOS was advanced by 0.56 days/decade, while the EOS delayed trend rate was 5.5 days/decade, and the growing season was significantly longer by 6.06 days/decade in North America [42]. In the Tibetan Plateau region, SOS was advanced at a rate of 0.17 days/decade, EOS was delayed at a rate of 5.29 days/decade, and LOS was lengthened at a rate of 5.46 days/decade for the period 1981–2017 [41]. These results are not entirely consistent with those of studies conducted for the QMs, and the differences may be mainly due to their different target periods and the different methods of phenology extraction. However, other investigations reported that, compared to 1982–1999, the phenology trend slowed down in 2000–2008 and the changes were not highly significant [14,43]. Meanwhile, our results also show a slower change in phenology trends over the last 20 years in the QMs, and the magnitude of SOS advance is also smaller than that of EOS delay. This observation is similar to the results of Wang et al. and Xia et al. [33,44], which indicates that the satellite-observed phenology change rates slowed down during a global warming hiatus between 1998 and 2012.

We also found that the trend in phenology changes of different vegetation is also highly variable, which is related to the microclimate of different regions or the geographical variation of plant origin [45]. Our results show a significantly greater trend of SOS changes in SL than other vegetation types, indicating that the earlier the SOS, the more significant the trend in SOS variations, and that this difference may be related to plant pollination type, life type, phylogenetic and wood type, etc. [46]. However, the trend of SOS changes in CL was weaker than in other vegetation types, and this trend was insignificant because farmers controlled the sowing time in each year, resulting in significantly smaller variability in crop phenology than in field observed plants [47]. The results show that the trend of delayed EOS is significantly stronger for EBF than other vegetation types, mainly because EBF is mainly located in the region south of the QMs, with a humid northern subtropical climate. Some researchers have shown that in subtropical mountainous and hilly areas, broadleaved forests can grow longer under the same similar climatic conditions compared to coniferous forests [48]. The study also found a significantly stronger trend in the lengthening of growing season for trees than for shrubs and herbs, which is the same as the findings of Zhu et al. [49] but in contrast to those of Ge et al. [50], who reported that the interannual variation trend for trees in China from the 1960s to the 2000s was significantly weaker than for herbaceous plants, and this difference in trend was due to differences in the study area and the species of the plants themselves.

4.2. Asymmetry in Contributions of SOS and EOS Trends to LOS

We found the asymmetry in contributions of the SOS and EOS trends to LOS variations by counting the percentage of pixels with positive and negative C-index values. The results show that SOS trends control 48.4% of LOS variations and EOS trends control 51.4% of LOS variations, which shows a stronger association between EOS trends and LOS variations compared to SOS (Figure S1). Previous studies illustrated that the lengthened growing season was mainly driven by delayed autumn phenology, which is consistent with our results [14,38,49]. However, other researchers found that it is the changes in SOS, and not EOS, that dominate the changes in growing season length [51,52]. It can be seen that there are differences in previous studies regarding the attribution of LOS variations. To investigate the reasons for such differences, we divided the trends of SOS, EOS, and LOS into six change patterns (Table 2 and Figure 7). Our results show that in addition to the main growth pattern of Type I (SOS advanced, EOS delayed, and LOS lengthened), 15.2% of the regions had Type V (SOS delayed, EOS delayed, and LOS lengthened), and the delayed EOS was significant in this pattern. Another 12% of the regions showed Type II (SOS advanced, EOS advanced, and LOS lengthened) growth pattern, and advanced SOS was significant. However, since the percentage of the region of the growth pattern Type II is smaller than that of Types I and V, it is still the trend of delayed EOS that dominates the variation in LOS for the whole study area, leading to asymmetry of the relative contribution of SOS and EOS to LOS.

In addition, we found that ENF, EBF, MF, SL, and CL were all controlled by EOS trends, while the variations in LOS for two vegetation types, DBF and GL, were controlled by SOS trends (Table 5). As Figure 4e shows, the growing season lengths of DBF and GL are short, and there are previous studies demonstrating that the effect of EOS shifts on vegetation with short growing season cycles is insignificant [53]. The percentage of growth pattern type II is higher than other vegetation types in DBF and GL, and the advance in SOS is also significant, resulting in SOS dominating the variation in LOS (Table 4). Therefore, we suggest that the asymmetry in SOS and EOS trends contributing to LOS is related to vegetation types, and that future studies should focus on vegetation types to accurately model and predict vegetation phenology periods.

4.3. Analysis of the Drivers of Interannual Variations in LSP

Previous studies showed that the interaction of meteorological, soil, and biological factors influenced the interannual variability of LSP [6,54]. Our results suggest that SWP is the most important driver of interannual variations in SOS and EOS across the QMs (Figure 8). This is mainly because shortwave radiation compensates for the lack of chilling demand during plant physiological dormancy through day length, i.e., longer daylight hours, and has an critical effect on SOS by delaying the accumulation of abscisic acid and slowing down the rate of leaf senescence, and also on EOS [21]. We also found that SWP and MD contributed a total of 42.3% to the interannual variations in SOS, and besides SWP, MD was also an important driver of the interannual SOS variation (Table 6 and Table S1). This is due to the fluctuation of the time interval between SOS and MD, which depends on the different developmental stages of the phenology events to a large extent and on the specific differences in the life history of the plant, and needs to be explained by the phenotypic plasticity of the individual and the adaptation to the environment [55]. Therefore, we suggest that the effect of MD on interannual SOS variation varies considerably among vegetation individuals. STP also contributed 12.0% importance in explaining the interannual SOS variation (Table S1). This is mainly due to the increased soil temperature, which accelerated the rate of leaf tip emergence and whole leaf expansion, thus promoting SOS [56]. Moreover, SWP and PP together explain the importance of 43.2% of the interannual EOS variation (Table S1). This is mainly because preseason shortwave radiation and precipitation control the availability of sunlight and water in vegetation, respectively, and reduced precipitation affects water transport capacity, which limits the photosynthetic rate of leaf, leading to lower utilization of light and water conditions by plants and affecting the interannual variation in EOS [20,57]. The combined contribution of MD, SWP, and PP to interannual EOS variation was also found to be as high as 54.0% (Table S1), suggesting that the lifecycle of vegetation is strongly regulated by its own rhythms under improved hydrothermal conditions, and that biological rhythms play a critical role in interannual EOS variation [7].

Furthermore, our study shows that the effect of each driver on interannual variations in LSP was varied for different vegetation types (Figure 9 and Table 6). For example, SWP is the most important driver for ENF, MF, GL, and CL; MD is the most important driver for EBF and DBF; and STP is the most important driver for SL (Table 6). This difference is mainly due to the diversity of plant physiological structures and the different adaptive strategies of plants to environmental changes [58]. SWP has the greatest effect on the interannual SOS variation in GL, which is mainly distributed in higher parts of the QMs and receives abundant solar radiation. The strong solar radiation promotes root activity and advances SOS [59]. The greatest contribution of MD to interannual SOS variation in EBF is related to the fact that EBF grows mainly in the southern part of the QMs, where its deeply rooted system and water conservation adaptations combine to reduce water stress under the influence of a humid northern subtropical climate. This adaptation to environmental changes is strongly regulated by its own rhythms, such that EBF is most affected by MD [60]. The effect of STP in SL is mainly due to the preseason accumulation of soil temperatures susceptible to specific thresholds that accelerate soil thaw and vegetation wake, triggering SOS [61]. SWP was the main driver of interannual EOS variation for all vegetation types, with GL being most influenced by SWP (Table 6). This is because abundant solar radiation increases surface evaporation and reduces water availability in grasslands, which subsequently inhibits vegetation growth, resulting in the EOS of GL being most influenced by SWP [62]. PP and MD have the strongest effect on interannual EOS variation in SL (Table 6). To our best knowledge, SL is mainly distributed in semi-humid and semiarid areas, and the control of plant metabolism by water stress affects its transpiration and photosynthesis, resulting in impaired ATPase synthesis and accelerated chlorophyll degradation. Meanwhile the adaptation of vegetation to such adversity changes also affects the interannual variations in EOS [63]. The relationship between regional climate and vegetation phenology growth will be further explored in future studies.

4.4. Evaluation of RF Model

We validated the accuracy of our RF model (Figure S2). The R² values of the linear regression formed by the predicted values and the observations inversions are both greater than 0.900 and 0.911, respectively. When predicting the dates of SOS or EOS, both RMSE and MAE are relatively small, showing that our RF model has good predictive performance. Machine learning, as a nonparametric multivariate approach, can integrate complex relationships between multiple spatial and temporal LSP dates and climate into a single model for predicting SOS and EOS. Then, the main drivers of SOS and EOS can be identified by estimating the importance of each variable [6]. However, it should be noted that although we have tried our best to adjust the hyperparameters of the algorithm to prevent overfitting in the model, some errors still appear in the test datasets (Table S2). Therefore, to obtain a better fit, it is necessary to further refine the study area and tree species in the future. Further study should compare different algorithms to better simulate the phenology period.

5. Conclusions

This study used the phenology metrics of vegetation in the QMs extracted from satellite NDVI data to analyze the spatiotemporal trends of LSP during 2001–2019, and to identify the dominant growth patterns of different vegetation types during the growing season. Furthermore, driving factors influencing interannual variations in LSP were emphatically investigated using the RF model. The main conclusions were as follows:

• (1). The average advance of SOS across QMs was 1.5 days/decade, with a significant advance in SOS observed for 27.5% of pixels. EOS was delayed by 2.4 days/decade, with a significant delay in EOS observed for 42.1% of pixels. LOS was lengthened by 3.9 days/decade, with a significant LOS lengthening observed for 40.3% of pixels.

• (2). The dominant pattern of change in the growing season for different vegetation types was advanced SOS, delayed EOS, and lengthened LOS, and this pattern had the highest percentage in evergreen broadleaved forests. The percentage of area shows that the patterns of delayed SOS and EOS and lengthened LOS were the highest percentage in shrubs.

• (3). For the whole QMs, LOS changes were mainly controlled by EOS, and the percentage was 51.4%. For deciduous broadleaved forests and grasses, LOS changes were attributed to SOS, while for other vegetation types, they were attributed to EOS.

• (4). SWP was found to be the most important factor influencing SOS and EOS, and grass and crop most influenced by SWP. Interannual variations in SOS were more influenced by biological factors (MD) than in EOS. The interannual variability of EOS is more influenced by preseason precipitation (PP) than SOS.

Author Contributions

J.G.: Formal analysis, Methodology, Writing—original draft. X.L.: Conceptualization, Funding acquisition, Methodology, Supervision, Writing—original draft. W.G.: Data curation, Visualization. X.N.: Data curation, Visualization. W.M.: Writing—review & editing. Q.L.: Writing—review & editing. X.X.: Resources. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the National Natural Science Foundation of China (Grants 41971104) and partially supported by the Fundamental Research Funds for the Central Universities (Projects GK202107009 and 2021TS015).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The datasets generated and/or analyzed during the study are available from the corresponding author upon reasonable request.

Acknowledgments

I am very grateful to the reviewers and editors for their suggestions on the revision of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

LSP	Land surface phenology
QMs	Qinling Mountains
NDVI	Normalized difference vegetation index
SOS	The start of the growing season
EOS	The end of the growing season
LOS	The length of the growing season
ENF	Evergreen needleleaved forest
EBF	Evergreen broadleaved forest
DBF	Deciduous broadleaved forest
MF	Mixed forest
SL	Shrubland
GL	Grassland
CL	Cropland
TP	Preseason average temperature
TG	Growing season average temperature
PP	Preseason total precipitation
PG	Growing season total precipitation
SWP	Preseason mean shortwave radiation
SWG	Growing season mean shortwave radiation
STP	Preseason soil temperature
STG	Growing season soil temperature
SMP	Preseason soil moisture
SMG	Growing season soil moisture
MD	The middle date of the growing season
MN	Maximum NDVI during growing season
RF	Random forest
OOB	Out of bag

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. An overview of (a) the study area and (b) spatial distribution of vegetation types. Only pixels with unchanged vegetation types were analyzed in this study period. The dataset with a total of 218,899 pixels was divided into seven vegetation types.

Figure 2. The description of phenology metrics correlations extracted using the NDVI time series datasets.

Figure 3. Flowchart of the method used in this study. MODIS 13A2 NDVI (2001–2019), in which the spatiotemporal distributions of NDVI for 2001, 2010, 2015 and 2019 are presented as examples of NDVI time series used as extract phenology metrics. Land cover (2001–2019), in which the spatiotemporal distribution of land cover for 2001, 2010, 2015 and 2019, used as an example of the changes in different vegetation types. Environmental factors, the first column indicates the drivers for SOS and the second column indicates the drivers for EOS, and shows the spatial distribution of each driver separately, which is used to example the predictors used for RF models.

Figure 4. (a,c,e) Spatial distribution of the average phenology metrics from 2001 to 2019 and (b,d,f) standard deviation (Std) of the phenology metrics. Insets at bottom left show the histogram of the average pixel values for different vegetation types.

Figure 5. (a,c,e) Spatial distribution of phenology metrics trends from 2001 to 2019 and (b,d,f) significant (p < 0.1) changes in phenology metrics trends for the study periods. Insets at bottom left show the average pixel values of the trends of the phenology metrics for the different vegetation types.

Figure 6. Interannual variations and significance of the mean phenology metrics for the entire study area and the areas covered by different vegetation types. The unit of Sen’s slope is days/year.

Figure 7. Trends and significance of each pattern during the growing season for different vegetation types. (a) growing season pattern I. (b) growing season pattern II. (c) growing season pattern III. (d) growing season pattern IV. (e) growing season pattern V. (f) growing season pattern VI. Note: “All” refers to all vegetation in the QMs.

Figure 8. Relative importance of drivers of affecting interannual variations in LSP for the entire study area, in decreasing order. (a) The ranking of the drivers of affecting interannual variations in SOS. (b) The ranking of the drivers of affecting interannual variations in EOS. The specific values for the relative importance of each driver are in Table S1. Note: The abbreviated variable names are the same as in Table 3.

Figure 9. The relative importance of drivers that affect interannual variations in LSP for different vegetation types. Relative importance was derived from the importance scores of variables (VI) based on RF models established for different vegetation types: (a) SOS. (b) EOS. Different colors indicate different factors. Note: SOS or EOS as response variable in the RF model. The specific values for the relative importance of each driver are in Table S1.

Table 1

Datasets and sources in the study area.

Dataset	Spatial Resolution	Temporal Resolution	Time Span	Source
MODIS13A2 NDVI	1 km	16 days	2001–2019	The Level-1 and Atmosphere Archive and Distribution System Distributed Active Archive Center (LAADS DAAC) (https://search.earthdata.nasa.gov/search/, accessed on 15 October 2020).
Land cover(CCI-LC)	300 m	Yearly	2001–2019	http://maps.elie.ucl.ac.be/CCI/viewer/index.php, accessed on 20 October 2020
Temperature	0.1°	hourly	2001–2019	The Reanalysis (ERA5) climatic datasets (https://cds.climate.copernicus.eu, accessed on 12 November 2020)
Precipitation	0.1°	hourly	2001–2019
Shortwave radiation	0.1°	hourly	2001–2019
Soil temperature	0.1°	hourly	2001–2019
Soil moisture	0.1°	hourly	2001–2019

Table 2

Six change patterns of plant growing seasons, i.e., six combination types of SOS, EOS, and LOS change trends. The plus and minus signs represent the trend direction corresponding to SOS, EOS, or LOS, respectively.

Change Pattern	Trend of SOS	Trend of EOS	Trend of LOS
I	Advanced (−)	Delayed (+)	Lengthened (+)
II	Advanced (−)	Advanced (−)	Lengthened (+)
III	Advanced (−)	Advanced (−)	Shortened (−)
IV	Delayed (+)	Advanced (−)	Shortened (−)
V	Delayed (+)	Delayed (+)	Lengthened (+)
VI	Delayed (+)	Delayed (+)	Shortened (−)

Table 3

Predictive variables used in the modeling of the LSP dates. The 12 predictive variables were classified into three categories: meteorological factors, soil factors, and biological factors. Meteorological factors include TP, TG, PP, PG, SWP, and SWG (6 in total). Soil factors include STP, STG, SMP, and SMG (4 in total). Biological factors include MN and MD (2 in total). Growing season is defined as the days between SOS and EOS.

Variables	SOS Drivers	EOS Drivers
Meteorological factors	Preseason average temperature * (TP)	Preseason average temperature ** (TP)
	Growing season average temperature (TG)	Growing season average temperature (TG)
	Preseason total precipitation * (PP)	Preseason total precipitation ** (PP)
	Growing season total precipitation (PG)	Growing season total precipitation (PG)
	Preseason mean shortwave radiation * (SWP)	Preseason mean shortwave radiation ** (SWP)
	Growing season mean shortwave radiation (SWG)	Growing season mean shortwave radiation (SWG)
Soil factors	Preseason soil temperature * (STP)	Preseason soil temperature ** (STP)
	Growing season soil temperature (STG)	Growing season soil temperature (STG)
	Preseason soil moisture * (SMP)	Preseason soil moisture ** (SMP)
	Growing season soil moisture (SMG)	Growing season soil moisture (SMG)
Biological factors	Maximum NDVI during growing season (MN)	Maximum NDVI during growing season (MN)
	Middle season date (MD)	Middle season date (MD)

* Predicted over a period of 30 days. ** Predicted over a period of 60 days.

Table 4

The percentage of phenology metrics datasets consisting of each pixel showing different change patterns in the growing seasons for all the vegetation types. Change patterns refer to the trend groupings in Table 2.

Change Patterns	All the Vegetation Types	ENF	EBF	DBF	MF	SL	GL	CL
I	48.4%	46.6%	50.1%	49.3%	47.2%	43.2%	44.8%	48.7%
II	12.0%	10.9%	9.5%	13.7%	11.8%	6.9%	15.2%	10.5%
III	7.3%	6.8%	6.2%	8.6%	7.3%	3.1%	7.6%	5.9%
IV	8.7%	9.0%	7.7%	7.9%	8.8%	15.5%	12.2%	9.3%
V	15.2%	16.8%	18.0%	12.4%	16.5%	23.5%	12.7%	17.2%
VI	8.5%	9.9%	8.5%	8.1%	8.4%	7.8%	7.5%	8.6%
Total	100.0%	100.0%	100.0%	100.0%	100.0%	100.0%	100.0%	100.0%

Table 5

The percentage of datasets in which LOS change was primarily attributable to the shift in SOS or EOS. For most vegetation types, the percentages show that the trend in LOS was mainly controlled by the shift in EOS.

All the Vegetation Types	SOS Controlled	EOS Controlled	Total
ENF	46.9%	53.1%	100%
EBF	43.4%	56.6%	100%
DBF	53.4%	46.6%	100%
MF	47.2%	52.8%	100%
SL	41.7%	58.3%	100%
GL	52.0%	48.0%	100%
CL	45.9%	54.1%	100%
The whole area	48.6%	51.4%	100%

Table 6

The top three dominant drivers affecting interannual variations in LSP. These three dominant factors are derived from the ranking of the importance scores of the variables (VI). Different vegetation types have different dominant drivers.

LSP	All the Vegetation Types	First Dominant Driver	Second Dominant Driver	Third Dominant Driver
SOS	ENF	SWP	MD	STP
	EBF	MD	SWP	PP
	DBF	MD	SWP	STP
	MF	SWP	MD	STP
	SL	STP	TP	SWP
	GL	SWP	MD	STP
	CL	SWP	MD	TG
	the whole area	SWP	MD	STP
EOS	ENF	SWP	TP	PP
	EBF	SWP	MD	PP
	DBF	SWP	PP	TP
	MF	SWP	PP	MD
	SL	SWP	PP	MD
	GL	SWP	SWG	TP
	CL	SWP	TP	PP
	the whole area	SWP	PP	MD

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/rs13224538/s1, Figure S1: The primary factor (SOS or EOS) controlling change in LOS across the entire study area calculated by the C-index, Figure S2: Red fitting curves between predicted LSP dates from RF model and observed LSP dates from MODIS13A2-NDVI, Table S1: The relative importance of each driver affecting interannual variation in LSP, Table S2: Correlation results of random forest models built for different vegetation type cover areas.

References

1. Güsewell, S.; Furrer, R.; Gehrig, R.; Pietragalla, B. Changes in temperature sensitivity of spring phenology with recent climate warming in Switzerland are related to shifts of the preseason. Glob. Chang. Biol.; 2017; 23, pp. 5189-5202. [DOI: https://dx.doi.org/10.1111/gcb.13781]

2. Fu, Y.H.; Piao, S.; Vitasse, Y.; Zhao, H.; De Boeck, H.J.; Liu, Q.; Yang, H.; Weber, U.; Hänninen, H.; Janssens, I.A. Increased heat requirement for leaf flushing in temperate woody species over 1980–2012: Effects of chilling, precipitation and insolation. Glob. Chang. Biol.; 2015; 21, pp. 2687-2697. [DOI: https://dx.doi.org/10.1111/gcb.12863]

3. Lang, W.; Chen, X.; Qian, S.; Liu, G.; Piao, S. A new process-based model for predicting autumn phenology: How is leaf senescence controlled by photoperiod and temperature coupling?. Agric. For. Meteorol.; 2019; 268, pp. 124-135. [DOI: https://dx.doi.org/10.1016/j.agrformet.2019.01.006]

4. Thackeray, S.J.; Henrys, P.A.; Hemming, D.; Bell, J.R.; Botham, M.S.; Burthe, S.; Helaouet, P.; Johns, D.G.; Jones, I.D.; Leech, D.I. et al. Phenological sensitivity to climate across taxa and trophic levels. Nature; 2016; 535, pp. 241-245. [DOI: https://dx.doi.org/10.1038/nature18608] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27362222]

5. Zeng, Z.; Piao, S.; Li, L.; Zhou, L.; Ciais, P.; Wang, T.; Li, Y.; Lian, X.; Wood, E.; Mao, J. Climate mitigation from vegetation biophysical feedbacks during the past three decades. Nat. Clim. Chang.; 2017; 7, pp. 432-436. [DOI: https://dx.doi.org/10.1038/nclimate3299]

6. Rodriguez-Galiano, V.F.; Sanchez-Castillo, M.; Dash, J.; Atkinson, P.M.; Ojeda-Zujar, J. Modelling interannual variation in the spring and autumn land surface phenology of the European forest. Biogeosciences; 2016; 13, pp. 3305-3317. [DOI: https://dx.doi.org/10.5194/bg-13-3305-2016]

7. Zu, J.; Zhang, Y.; Huang, K.; Liu, Y.; Chen, N.; Cong, N. Biological and climate factors co-regulated spatial-temporal dynamics of vegetation autumn phenology on the Tibetan Plateau. Int. J. Appl. Earth Obs. Geoinf.; 2018; 69, pp. 198-205. [DOI: https://dx.doi.org/10.1016/j.jag.2018.03.006]

8. Keenan, T.F.; Gray, J.; Friedl, M.A.; Toomey, M.; Bohrer, G.; Hollinger, D.Y.; Munger, J.W.; Keefe, J.O.; Schmid, H.P.; Wing, I.S. Net carbon uptake has increased through warming-induced changes in temperate forest phenology. Nat. Clim. Chang.; 2014; 4, pp. 598-604. [DOI: https://dx.doi.org/10.1038/nclimate2253]

9. Wu, C.; Hou, X.; Peng, D.; Gonsamo, A.; Xu, S. Land surface phenology of China’s temperate ecosystems over 1999–2013: Spatial–temporal patterns, interaction effects, covariation with climate and implications for productivity. Agric. For. Meteorol.; 2016; 216, pp. 177-187. [DOI: https://dx.doi.org/10.1016/j.agrformet.2015.10.015]

10. Beurs, K.M.D.; Henebry, G.M. Land surface phenology and temperature variation in the International Geosphere-Biosphere Program high-latitude transects. Glob. Chang. Biol.; 2010; 11, pp. 779-790. [DOI: https://dx.doi.org/10.1111/j.1365-2486.2005.00949.x]

11. Cleland, E.E.; Chuine, I.; Menzel, A.; Mooney, H.A.; Schwartz, M.D. Shifting plant phenology in response to global change. Trends Ecol. Evol.; 2007; 22, pp. 357-365. [DOI: https://dx.doi.org/10.1016/j.tree.2007.04.003]

12. Stoeckli, R.; Vidale, P.L. European plant phenology and climate as seen in a 20-year AVHRR land-surface parameter dataset. Int. J. Remote Sens.; 2004; 25, pp. 3303-3330. [DOI: https://dx.doi.org/10.1080/01431160310001618149]

13. Piao, S.; Fang, J.; Zhou, L.; Ciais, P.; Zhu, B. Variations in satellite-derived phenology in China’s temperate vegetation. Glob. Chang. Biol.; 2010; 12, pp. 672-685. [DOI: https://dx.doi.org/10.1111/j.1365-2486.2006.01123.x]

14. Jeong, S.J.; Chang-Hoi, H.O.; Gim, H.J.; Brown, M.E. Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982–2008. Glob. Chang. Biol.; 2011; 17, pp. 2385-2399. [DOI: https://dx.doi.org/10.1111/j.1365-2486.2011.02397.x]

15. Shen, M.; Tang, Y.; Chen, J.; Zhu, X.; Zheng, Y. Influences of temperature and precipitation before the growing season on spring phenology in grasslands of the central and eastern Qinghai-Tibetan Plateau. Agric. For. Meteorol.; 2011; 151, pp. 1711-1722. [DOI: https://dx.doi.org/10.1016/j.agrformet.2011.07.003]

16. Yuan, M.; Zhao, L.; Lin, A.; Li, Q.; Qu, S. How do climatic and non-climatic factors contribute to the dynamics of vegetation autumn phenology in the Yellow River Basin, China?. Ecol. Indic.; 2020; 112, 106112. [DOI: https://dx.doi.org/10.1016/j.ecolind.2020.106112]

17. Wang, C.; Guo, H.; Zhang, L.; Liu, S.; Qiu, Y.; Sun, Z. Assessing phenological change and climatic control of alpine grasslands in the Tibetan Plateau with MODIS time series. Int. J. Biometeorol.; 2015; 59, pp. 11-23. [DOI: https://dx.doi.org/10.1007/s00484-014-0817-5]

18. Flynn, D.; Wolkovich, E.M. Temperature and photoperiod drive spring phenology across all species in a temperate forest community. New Phytol.; 2018; 219, pp. 1353-1362. [DOI: https://dx.doi.org/10.1111/nph.15232]

19. Liu, Q.; Piao, S.; Fu, Y.H.; Gao, M.; Peñuelas, J.; Janssens, I.A. Climatic Warming Increases Spatial Synchrony in Spring Vegetation Phenology Across the Northern Hemisphere. Geophys. Res. Lett.; 2019; 46, pp. 1641-1650. [DOI: https://dx.doi.org/10.1029/2018GL081370]

20. Prevéy, J.S.; Seastedt, T.R. Seasonality of precipitation interacts with exotic species to alter composition and phenology of a semi-arid grassland. J. Ecol.; 2014; 102, pp. 1549-1561. [DOI: https://dx.doi.org/10.1111/1365-2745.12320]

21. Way, D.A.; Montgomery, R.A. Photoperiod constraints on tree phenology, performance and migration in a warming world. Plant Cell Environ.; 2015; 38, pp. 1725-1736. [DOI: https://dx.doi.org/10.1111/pce.12431] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25142260]

22. Ztürk, M.; Bolat, L.; Ergün, A. Influence of air-soil temperature on leaf expansion and LAI of Carpinus betulus trees in a temperate urban forest patch. Agric. For. Meteorol.; 2015; 200, pp. 185-191. [DOI: https://dx.doi.org/10.1016/j.agrformet.2014.09.014]

23. Jin, Z.; Zhuang, Q.; He, J.S.; Luo, T.; Shi, Y. Phenology shift from 1989 to 2008 on the Tibetan Plateau: An analysis with a process-based soil physical model and remote sensing data. Clim. Chang.; 2013; 119, pp. 435-449. [DOI: https://dx.doi.org/10.1007/s10584-013-0722-7]

24. Zhou, Y. Asymmetric Behavior of Vegetation Seasonal Growth and the Climatic Cause: Evidence from Long-Term NDVI Dataset in Northeast China. Remote Sens.; 2019; 11, 2107. [DOI: https://dx.doi.org/10.3390/rs11182107]

25. Yang, M.; Nelson, F.E.; Shiklomanov, N.I.; Guo, D.; Wan, G. Permafrost degradation and its environmental effects on the Tibetan Plateau: A review of recent research. Earth-Sci. Rev.; 2010; 103, pp. 31-44. [DOI: https://dx.doi.org/10.1016/j.earscirev.2010.07.002]

26. Park, H.; Jeong, S.J.; Ho, C.H.; Kim, J.; Brown, M.E.; Schaepman, M.E. Nonlinear response of vegetation green-up to local temperature variations in temperate and boreal forests in the Northern Hemisphere. Remote Sens. Environ.; 2015; 165, pp. 100-108. [DOI: https://dx.doi.org/10.1016/j.rse.2015.04.030]

27. Breiman, L. Random forests. Mach. Learn.; 2001; 45, pp. 5-32. [DOI: https://dx.doi.org/10.1023/A:1010933404324]

28. Li, Z.; Liang, M.; Li, Z.; Pierre, M.; Tong, X.; Zhang, J.; Dong, L.; Zheng, Y.; Ma, W.; Zhao, L. Plant functional groups mediate effects of climate and soil factors on species richness and community biomass in grasslands of Mongolian Plateau. J. Plant Ecol.; 2021; 14, pp. 679-691. [DOI: https://dx.doi.org/10.1093/jpe/rtab021]

29. Zhao, M.; Peng, C.; Xiang, W.; Deng, X.; Tian, D.; Zhou, X.; Yu, G.; He, H.; Zhao, Z. Plant phenological modeling and its application in global climate change research: Overview and future challenges. Environ. Rev.; 2013; 21, pp. 1-14. [DOI: https://dx.doi.org/10.1139/er-2012-0036]

30. Deng, C.; Bai, H.; Ma, X.; Zhao, T.; Huang, X. Spatiotemporal differences in the climatic growing season in the Qinling Mountains of China under the influence of global warming from 1964 to 2015. Theor. Appl. Climatol.; 2019; 138, pp. 1899-1911. [DOI: https://dx.doi.org/10.1007/s00704-019-02886-w]

31. Bai, W.G.; Zhang, L. Phytocoenological characteristics and community classification of Abies fargesii forests in Qinling Mountains. J. Beijing For. Univ.; 2007; S2, pp. 222-226.

32. Myneni, R.B.; Keeling, C.D.; Tucker, C.J.; Asrar, G.; Nemani, R.R. Increased plant growth in the northern high latitudes from 1981 to 1991. Nature; 1997; 386, pp. 698-702. [DOI: https://dx.doi.org/10.1038/386698a0]

33. Wang, X.; Xiao, J.; Li, X.; Cheng, G.; Ma, M.; Zhu, G.; Altaf Arain, M.; Andrew Black, T.; Jassal, R.S. No trends in spring and autumn phenology during the global warming hiatus. Nat. Commun.; 2019; 10, 2389. [DOI: https://dx.doi.org/10.1038/s41467-019-10235-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31160586]

34. Cao, R.; Yang, C.; Shen, M.; Chen, J.; Ji, Z.; Cong, W.; Wei, Y. A simple method to improve the quality of NDVI time-series data by integrating spatiotemporal information with the Savitzky-Golay filter. Remote Sens. Environ.; 2018; 217, pp. 244-257. [DOI: https://dx.doi.org/10.1016/j.rse.2018.08.022]

35. Kaur, R.; Kiran, G.S.; Shah, M.N.; Mistry, N.V.; Mohan, S. Applicability of Smoothing Techniques in Generation of Phenological Metrics of Tectona grandis L. Using NDVI Time Series Data. Remote Sens.; 2021; 13, 3343.

36. Gocic, M.; Trajkovic, S. Analysis of changes in meteorological variables using Mann-Kendall and Sen’s slope estimator statistical tests in Serbia. Glob. Planet. Chang.; 2013; 100, pp. 172-182. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2012.10.014]

37. Mann, H.B. Nonparametric test against trend. Econometrica; 1945; 13, pp. 245-259. [DOI: https://dx.doi.org/10.2307/1907187]

38. Garonna, I.; Jong, R.D.; Wit, A.D.; Mücher, C.A.; Schmid, B.; Schaepman, M.E. Strong contribution of autumn phenology to changes in satellite-derived growing season length estimates across Europe (1982–2011). Glob. Chang. Biol.; 2015; 20, pp. 3457-3470. [DOI: https://dx.doi.org/10.1111/gcb.12625]

39. Abdel-Rahman, E.M.; Ahmed, F.B.; Ismail, R. Random forest regression and spectral band selection for estimating sugarcane leaf nitrogen concentration using EO-1 Hyperion hyperspectral data. Int. J. Remote Sens.; 2013; 34, pp. 712-728. [DOI: https://dx.doi.org/10.1080/01431161.2012.713142]

40. Tian, F.; Cai, Z.; Jin, H.; Hufkens, K.; Scheifinger, H.; Tagesson, T.; Smets, B.; Van Hoolst, R.; Bonte, K.; Ivits, E. et al. Calibrating vegetation phenology from Sentinel-2 using eddy covariance, PhenoCam, and PEP725 networks across Europe. Remote Sens. Environ.; 2021; 260, 112456. [DOI: https://dx.doi.org/10.1016/j.rse.2021.112456]

41. Sun, Q.; Li, B.; Zhou, G.; Jiang, Y.; Yuan, Y. Delayed autumn leaf senescence date prolongs the growing season length of herbaceous plants on the Qinghai-Tibetan Plateau. Agric. For. Meteorol.; 2020; 284, 107896. [DOI: https://dx.doi.org/10.1016/j.agrformet.2019.107896]

42. Wang, X.; Piao, S.; Xu, X.; Ciais, P.; Macbean, N.; Myneni, R.B.; Li, L. Has the advancing onset of spring vegetation green-up slowed down or changed abruptly over the last three decades?. Glob. Ecol. Biogeogr.; 2015; 24, pp. 621-631. [DOI: https://dx.doi.org/10.1111/geb.12289]

43. Zeng, H.; Jia, G.; Epstein, H. Recent changes in phenology over the northern high latitudes detected from multi-satellite data. Environ. Res. Lett.; 2011; 6, pp. 45508-45518. [DOI: https://dx.doi.org/10.1088/1748-9326/6/4/045508]

44. Xia, H.; Qin, Y.; Feng, G.; Meng, Q.; Liu, G. Forest Phenology Dynamics to Climate Change and Topography in a Geographic and Climate Transition Zone: The Qinling Mountains in Central China. Forests; 2019; 10, 1007. [DOI: https://dx.doi.org/10.3390/f10111007]

45. Wang, H.; Ge, Q.; Dai, J.; Tao, Z. Geographical pattern in first bloom variability and its relation to temperature sensitivity in the USA and China. Int. J. Biometeorol.; 2015; 59, pp. 961-969. [DOI: https://dx.doi.org/10.1007/s00484-014-0909-2]

46. Panchen, Z.A.; Primack, R.B.; Nordt, B.; Ellwood, E.R.; Stevens, A.D.; Renner, S.S.; Willis, C.G.; Fahey, R.; Whittemore, A.; Du, Y. Leaf out times of temperate woody plants are related to phylogeny, deciduousness, growth habit and wood anatomy. New Phytol.; 2014; 203, pp. 1208-1219. [DOI: https://dx.doi.org/10.1111/nph.12892]

47. Kurukulasuriya, P.; Mendelsohn, R. Impact And Adaptation Of South-East Asian Farmers To Climate Change: Conclusions And Policy Recommendations. Clim. Chang. Econ.; 2017; 8, 1740007. [DOI: https://dx.doi.org/10.1142/S2010007817400073]

48. Qiu, B.W.; Zhong, M. Spatiotemporal variability of vegetation phenology with reference to altitude and climate in the subtropical mountain and hill region, China. Chin. Sci. Bull.; 2013; 58, pp. 2883-2892. [DOI: https://dx.doi.org/10.1007/s11434-013-5847-6]

49. Zheng, Z.; Zhu, W.; Chen, G.; Jiang, N.; Fan, D.; Zhang, D. Continuous but diverse advancement of spring-summer phenology in response to climate warming across the Qinghai-Tibetan Plateau. Agric. For. Meteorol.; 2016; 223, pp. 194-202. [DOI: https://dx.doi.org/10.1016/j.agrformet.2016.04.012]

50. Ge, Q.; Wang, H.; Rutishauser, T.; Dai, J. Phenological response to climate change in China: A meta-analysis. Glob. Chang. Biol.; 2015; 21, pp. 265-274. [DOI: https://dx.doi.org/10.1111/gcb.12648]

51. Park, T.; Ganguly, S.; T Mmervik, H.; Euskirchen, E.S.; H Gda, K.A.; Karlsen, S.R.; Brovkin, V.; Nemani, R.R.; Myneni, R.B. Changes in growing season duration and productivity of northern vegetation inferred from long-term remote sensing data. Environ. Res. Lett.; 2016; 11, 84001. [DOI: https://dx.doi.org/10.1088/1748-9326/11/8/084001]

52. Shen, M.; Piao, S.; Dorji, T.; Liu, Q.; Cong, N.; Chen, X.; An, S.; Wang, S.; Wang, T.; Zhang, G. Plant phenological responses to climate change on the Tibetan Plateau: Research status and challenges. Natl. Sci. Rev.; 2015; 2, pp. 454-467. [DOI: https://dx.doi.org/10.1093/nsr/nwv058]

53. Cong, N.; Huang, K.; Zhang, Y. Unsynchronized Driving Mechanisms of Spring and Autumn Phenology Over Northern Hemisphere Grasslands. Front. For. Glob. Chang.; 2021; 3, 610162. [DOI: https://dx.doi.org/10.3389/ffgc.2020.610162]

54. Wang, J.; Sun, H.; Xiong, J.; He, D.; Cheng, W.; Ye, C.; Yong, Z.; Huang, X. Dynamics and Drivers of Vegetation Phenology in Three-River Headwaters Region Based on the Google Earth Engine. Remote Sens.; 2021; 13, 2528. [DOI: https://dx.doi.org/10.3390/rs13132528]

55. Mcwatters, H.G.; Devlin, P.F. Timing in plants—A rhythmic arrangement. Febs. Lett.; 2011; 585, pp. 1474-1484. [DOI: https://dx.doi.org/10.1016/j.febslet.2011.03.051]

56. Stone, P.J.; Sorensen, I.B.; Jamieson, P.D. Effect of soil temperature on phenology, canopy development, biomass and yield of maize in a cool-temperate climate. Field Crop. Res.; 1999; 63, pp. 169-178. [DOI: https://dx.doi.org/10.1016/S0378-4290(99)00033-7]

57. Chen, X.; Ciais, P.; Maignan, F.; Zhang, Y.; Zhang, H. Vapor Pressure Deficit and Sunlight Explain Seasonality of Leaf Phenology and Photosynthesis Across Amazonian Evergreen Broadleaved Forest. Glob. Biogeochem. Cycles; 2021; 35, e2020GB006893. [DOI: https://dx.doi.org/10.1029/2020GB006893]

58. Diez, J.M.; Ibáñez, I.; Miller-Rushing, A.J.; Mazer, S.J.; Crimmins, T.M.; Crimmins, M.A.; Bertelsen, C.D.; Inouye, D.W. Forecasting phenology: From species variability to community patterns. Ecol. Lett.; 2012; 15, pp. 545-553. [DOI: https://dx.doi.org/10.1111/j.1461-0248.2012.01765.x]

59. Li, K.; Wang, C.; Sun, Q.; Rong, G.; Tong, Z.; Liu, X.; Zhang, J. Spring Phenological Sensitivity to Climate Change in the Northern Hemisphere: Comprehensive Evaluation and Driving Force Analysis. Remote Sens.; 2021; 10, 1972. [DOI: https://dx.doi.org/10.3390/rs13101972]

60. Zhao, A.; Zhang, A.; Cao, S.; Liu, X.; Liu, J.; Cheng, D. Responses of vegetation productivity to multi-scale drought in Loess Plateau, China. Catena; 2017; 163, pp. 165-171. [DOI: https://dx.doi.org/10.1016/j.catena.2017.12.016]

61. Chuine, I. A unified model for budburst of trees. J. Theor. Biol.; 2000; 207, pp. 337-347. [DOI: https://dx.doi.org/10.1006/jtbi.2000.2178] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11082304]

62. Chen, B.Z.; Zhang, H.F.; Yan, J.W.; Wang, G.Y.; Zhou, T.M.; Xu, G.; Innes, J.L.; Che, M.L.; Dou, X.M. Spatial and temporal variations in the end date of the vegetation growing season throughout the Qinghai-Tibetan Plateau from 1982 to 2011. Agric. For. Meteorol.; 2014; 189–190, pp. 81-90.

63. Gerten, D.; Schaphoff, S.; Lucht, W. Potential future changes in water limitations of the terrestrial biosphere. Clim. Chang.; 2007; 80, pp. 277-299. [DOI: https://dx.doi.org/10.1007/s10584-006-9104-8]