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

: Land surface phenology (LSP), as a precise bio-indicator that responds to climate change, has received much attention in ﬁelds 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 signiﬁcant 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 quantiﬁed the degree of inﬂuence 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 signiﬁcantly in western QMs. Furthermore, we conﬁrmed 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 inﬂuenced by SWP. The soil condition has, overall, a minor inﬂuence the regional LSP. This study highlighted the speciﬁcity of different vegetation growth in the QMs under warming, which should be considered in the accurate prediction of vegetation growth in the future.


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 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.

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).

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).

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 where NDVI ratio(t) is the calculated relative changing rate of NDVI at time t and 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 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 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.

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): 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.

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. 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 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: 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.

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. 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.

SOS Drivers EOS Drivers
Meteorological factors 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: where X j is the jth variable, ntree is the number of trees, errOOB j t is the OOB error of each tree t, and err OOB j t 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 2 ) 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: where y i is the observed satellite-based SOS or EOS values, is y i the predicted RF-based SOS or EOS values, and 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. 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).

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).

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%).    Table S1. Note: The abbreviated variable names are the same as in Table 3.
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.  Table S1.
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.

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.

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.

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 (Tables 6 and 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.

Evaluation of RF Model
We validated the accuracy of our RF model ( Figure S2). The R 2 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.

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.
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. 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.