Forest Phenology Dynamics to Climate Change and Topography in a Geographic and Climate Transition Zone: The Qinling Mountains in Central China

Forest ecosystems in an ecotone and their dynamics to climate change are growing ecological and environmental concerns. Phenology is one of the most critical biological indicators of climate change impacts on forest dynamics. In this study, we estimated and visualized the spatiotemporal patterns of forest phenology from 2001 to 2017 in the Qinling Mountains (QMs) based on the enhanced vegetation index (EVI) from MODerate-resolution Imaging Spectroradiometer (MODIS). We further analyzed this data to reveal the impacts of climate change and topography on the start of the growing season (SOS), end of the growing season (EOS), and the length of growing season (LOS). Our results showed that forest phenology metrics were very sensitive to changes in elevation, with a 2.4 days delayed SOS, 1.4 days advanced EOS, and 3.8 days shortened LOS for every 100 m increase in altitude. During the study period, on average, SOS advanced by 0.13 days year−1, EOS was delayed by 0.22 days year−1, and LOS increased by 0.35 day year−1. The phenological advanced and delayed speed across different elevation is not consistent. The speed of elevation-induced advanced SOS increased slightly with elevation, and the speed of elevation-induced delayed EOS shift reached a maximum value of 1500 m from 2001 to 2017. The sensitivity of SOS and EOS to preseason temperature displays that an increase of 1 ◦C in the regionally averaged preseason temperature would advance the average SOS by 1.23 days and delay the average EOS by 0.72 days, respectively. This study improved our understanding of the recent variability of forest phenology in mountain ecotones and explored the correlation between forest phenology and climate variables in the context of the ongoing climate warming.


Introduction
Research on phenology focuses on the phases of recurring biological life cycles, how these cycles are influenced by biotic and abiotic drivers and the interrelations among the same or different species [1]. Forest phenology is one of the most reliable biological indicators for evaluating the impacts of global warming in temperate and boreal zones because of its high sensitivity to climate change [2]. In turn, variations in forest phenology strongly affect the fitness and distribution of temperate and boreal forests [3] and, thus, the interspecific interactions among plants [4], ecosystem structure and function by varying the length of growing season [5,6], terrestrial ecosystem water and carbon balances [7,8], and climate system energy balances [9]. However, the processes underlying forest phenology are still unclear and need further research as they are essential to enhance the understanding of forest ecosystems and their responses to climate change. Global warming has become an undeniable fact in recent years, changing the timings of plant flowering, leaf-out, growth stages, and leaf senescence over the past few decades [10]. Extensive research primarily focused on the start of the growing season (SOS) and found an earlier SOS in high latitudes and high altitudes based on the ground observations, satellite remote sensing, and model simulations [11][12][13][14][15][16][17]. However, the SOS trend diverges significantly among previous studies [18][19][20]. Due to the complexity and perplexity of vegetation phenological processes during autumn, studies on the processes of the end of the growing season (EOS) and its associated control factors have not received sufficient attention [21,22]. As the final phase of the plants' seasonal cycle, EOS serves as an important occasion to recycle nutrients from senescing leaves to other plant tissues to sustain their growth in the following spring [23,24]. The EOS was recognized to be as important as the SOS in influencing the interannual variability of the carbon balance by controlling the length of the photosynthetically active period [10,25,26]. Therefore, combining the SOS, EOS, and length of the growing season (LOS) for vegetation and its associated environmental factors will improve our understanding of the response of vegetation to climate change.
Temperature is a primary climatic factor regulating the SOS and EOS in temperate and boreal forests. Recent studies have indicated that spring warming directly regulates the SOS, and the winter chilling requirement for dormancy release indirectly regulates the SOS [13]. Warming during summer and autumn can significantly delay the EOS [27], but the response of EOS to warming differed among species [28]. The warming trends also affect the time of phenology; for example, many studies show that phenology trends are significantly faster pre-1990 than post-2000 [19,20,29,30]. The EOS has a different temperature sensitivity compared with SOS; for example, the EOS response to temperature is larger than that of SOS [27]. Precipitation has also been found to be a primary driving factor affecting the SOS and EOS in a complex manner in temperate and boreal forests [11,12,[31][32][33][34]. For example, Fu et al. [12] found that winter precipitation has a positive correlation with the heat requirement for SOS at northern middle and high latitudes. Precipitation also has a considerable impact on the EOS, particularly for arid and semi-arid regions [31]. Therefore, understanding the response of forest phenology to climate change is valuable. However, the processes underlying the forest phenology, which may differ among regions and ecosystems due to their associated climatic factors and the complex response of spring and autumn phenology to climate change in the mountains, are still not fully understood.
Mountain forests typically cover a concentrated and abundant vertical band spectrum of biological climates, providing important ecological service functions in maintaining biological diversity, regulating regional climate, and conserving water [35]. The elevational gradients can be treated as space-for-time substitutions to evaluate the influences of climate change on forest ecosystems. The environmental gradient of mountainous systems is unique, making it a sensitive indicator of climate change, and has become a focus for current climate change research. Meanwhile, warming trends may vary across the bioclimatic gradient, with stronger trends at higher latitudes and altitudes [36,37]. Therefore, biomes at high elevations are experiencing a more rapid change in temperature than those at low elevations [36]. Recent studies have shown trends of earlier SOS and delayed EOS in many mountain zones, including the Greater Khingan Mountains [38], the Changbai Mountains [39], the Andes Mountains [40], and semi-arid mountain regions [41], indicating that rising temperature could influence the phenological metrics. However, there is a lack of answers to whether the phenological advancement along different elevation gradients is uniform or diverse and how is the dynamics of uniform or diverse.
Satellite-based information has been used for monitoring the spatial heterogeneity of land surface phenology over recent decades and for discerning the spatiotemporal patterns of land surface phenology responding to the effects of climatic change [13,42]. Recent studies have effectively described the vegetation phenology using fine-resolution (<30 m), which can be achieved from fusing spare fine-and frequency medium-resolution images [43], from multiple years into a single synthetic year [44], and from optical HJ-1 satellites [45] and Sentinel-2 [46]. Although finer-resolution sensors are more effective for capturing the heterogeneous land cover, such sensors usually need longer revisit times with too few cloud-free images to effectively capture the seasonal changes in greenness [47]. According to the spatiotemporal resolutions and the length of time series, three satellite sensors, namely, Moderate Resolution Imaging Spectroradiometer (MODIS), Système Pour I'Observation de la Terre (SPOT), and Advanced Very High Resolution Radiometer (AVHRR), have often been used to monitor vegetation phenology [16,34]. Vegetation indices (e.g., EVI, Enhanced Vegetation Index; NDVI, Normalized Difference Vegetation Index) based on the spectral reflectance of vegetation are well correlated with chlorophyll abundance, photosynthetically active biomass, and energy absorption [48][49][50]. The EVI, which could avoid the problem of NDVI saturation in high vegetation coverage areas and reduce the effects of atmospheric and soil background, has proven to be a robust indicator of vegetation growth [51]. Therefore, a MODIS EVI dataset with a 500 m spatial resolution and an 8-day temporal resolution was used in this study for retrieving the forest phenology in the Qinling Mountains (QMs), China, taking into account the heterogeneity of the mountains and NDVI saturation in areas of dense forest coverage.
The QMs are located in a typical and important ecological and geographical transition zone in China. The transitional nature of the physical conditions, which is sensitive to external interference, provides an ideal region for examining the climatic and altitudinal effects on forest phenology metrics that is, hence, ideal for understanding the response of the regional forest ecosystems to ongoing global warming. We hypothesize that (1) the recent climate change has altered the forest growing season in the QMs, but the changes in SOS, EOS, and LOS along different elevation gradients are uneven; (2) temperature and precipitation are the key driving factors for SOS and EOS in the QMs. In this study, we investigated the relationships between topography, climate factors, and forest phenological metrics. In particular, we aim to (i) quantify the spatial pattern and interannual variability of the forest SOS, EOS, and LOS in the QMs from 2001 to 2017; (ii) evaluate the influence of elevation gradients on SOS, EOS, and LOS on a spatiotemporal scale; and (iii) explore the relationships between climatic factors (e.g., temperature, precipitation) and the spatiotemporal variations of forest phenology metrics.

Study Area
The QMs is located in central China (Figure 1a), and it features important geographical transitions from a subtropical climate in the south to a warm-temperate climate in the north. The area contains a set of important geographical boundaries, including a boundary between the north and south, a 0 • C isotherm in January, an 800 mm isohyet, and an isoline of 2000 h of sunshine duration [52]. The annual mean temperature is −1.5~16.3 • C from high to low altitudes (Figure 1b), the cumulative precipitation is 452-1189 mm from north to south (Figure 1c deciduous broad-leaved forest belt grows, a deciduous oak forest belt grows above 1300 to about 2300 m, from 2300 to 2800 m a typical birch forest belt grows, a coniferous forest belt grows beyond 2800 to about 3200 m, and above 3200 m is the highest forest ecosystem of an alpine shrub-meadow belt. Because of the significant climate differences between the south and north of the QMs, there is a mildly different distributional spectrum of vegetation types [53]. ecosystem of an alpine shrub-meadow belt. Because of the significant climate differences between the south and north of the QMs, there is a mildly different distributional spectrum of vegetation types [53].

MODIS EVI Dataset
The study area falls in MODIS path 26/row 05 and path 27/row 05. At each path/row, MOD09A1 scenes acquired in an annual year, with a temporal resolution of 8 days and a spatial resolution of 500 m from 2001 to 2017, were downloaded from the publicly available website of the USGS (United States Geological Survey) EROS (Earth Resources Observation and Science) [54]. Based on the heterogeneity of the Mountains, we selected the higher spatial resolution of 500 m with the maximum temporal resolution of 8 days. The original EVI data were mosaicked and reprojected from a sinusoidal projection to the Lambert Azimuthal Equal Area Projection. Each composite was clipped for the QMs area using the MODIS Reprojection Tool (MRT), and then, sequentially, were stacked to generate the yearly time-series dataset.

Climate Data
The meteorological dataset used in this study included daily mean temperature and cumulative precipitation from 2001 to 2017, which were downloaded from the China Meteorological Data Sharing Service System (http://data.cma.cn/). To obtain reliable spatially distributed climate data, 57 stations were selected covering not only the QMs but also the neighboring counties distributed throughout the south of Shanxi province, the southeast of Gansu province, the north of Sichuan province, the northwest of Hubei province, and the west of Henan province ( Figure 1). The daily meteorological observations for all the sites from 2001 to 2017 were used to obtain the monthly mean temperature and cumulative precipitation. To analyze the correlation between the phenology metrics and the climate data, the ANUSPLIN method [55], which measures topographical effects on spatial climate interpolation [36], was used to interpolate and rasterize monthly climate data from the stations into 500 × 500 m grids to match the spatial resolution of EVI.

Auxiliary Data
To analyze the topographic effects of phenology metrics, we downloaded SRTM (Shuttle Radar Topography Mission) DEM (digital elevation model) data, with a 90 m spatial resolution, from the CGIAR (Consultative Group for International Agriculture Research) Consortium for Spatial Information website (http://srtm.csi.cgiar.org/). Using the mean algorithm, the SRTM (90 m) imagery was resampled into a new image with a cell size of 500 m to match the cell size of EVI. Elevation was also used as a co-variable in the ANUSPLIN method.
To overcome the mosaic effect of MODIS land cover (MOD12Q1), we used Globeland30, which is a global land cover at a 30 m resolution with good performance in classification accuracy [56,57]. Globeland30 was downloaded from a webpage (http://www.globallandcover.com/GLC30Download), and it was then resampled to 500 m to match the aforementioned EVI data, using the majority method.

Phenology Metrics Extraction
The EVI time series are often hindered by the poor quality of the original datasets due to various factors, for example, haze and cloud contamination and atmospheric absorption. To suppress atmospheric effects, EVI was composited by the Maximum Value Compositing (MVC) algorithm. Although MVC is surprisingly effective in reducing noise in EVI data, it cannot remove all the noise. However, the seasonal forest canopy EVI data tend to change gradually from a seasonal perspective. Thus, it is necessary to adopt a filter method to eliminate the negative influence, to interpolate the missing value, and finally to construct new EVI time-series datasets before extracting seasonal phenology metrics.
In general, the types of input data and desired results determine the choice of the smoothing method. The HANTS (Harmonic ANalysis of Time Series) method, which can fully consider the plant growth cycle and the dual characteristics of the data themselves, is an effective analysis approach. This approach has been widely used to reconstruct the EVI time series and NDVI time series for the past 20 years [58,59]. The HANTS algorithm is applied using an iterative way (for more detail on the curve fitting, see Roerink, Menenti, and Verhoef [58]).
In our analysis of EVI from 2001 to 2017, the HANTS parameters (the definition of the parameters can be found in Supplementary (Text1) were set as follows: (1) the number of frequency is 2, (2) the suppression flag is low, (3) the invalid data are 0 for low threshold and 1 for high threshold, (4) the fit error tolerance is 0.05, (5) the degree of over-determinedness is 5, and (6) the damping factor is 0.5. The 8-day EVI data were interpolated to daily timesteps in order to capture the rapid changes in land surface phenology. Then, the phenology metrics (SOS, EOS) were extracted based on the reconstructed EVI time series ( Figure S1). Regarding the extraction of phenology metrics, it is difficult to define a single set of smoothing, and phenology metrics extraction methods can work with all different ecosystems and across all different kinds of remotely sensed data [60,61]. In this study, we chose the smoothing method of HANTS and the relative threshold method, which was successfully applied to the QMs [52]. The dates of SOS or EOS were defined, using the relative threshold method, as the times at which EVI values increased (SOS) or decreased (EOS) to a certain proportion of the seasonal amplitude, in the left or the right edge of the annual curve of EVI [34,62]. The LOS is equal to the difference between the EOS and the SOS. Such threshold values were determined according to the lowest absolute error between the leaf-out and leaf senescence dates measured in situ by the National Geographic Resource Science Data Center, National Science and Technology Infrastructure of China, National Earth System Science Data Sharing Infrastructure (http://gre.geodata.cn), and the remotely sensed phenology metrics ( Figure S2). Ultimately, 30% of the amplitude was used to extract the dates of SOS and EOS, respectively, at both sides of the annual curve for each pixel.

Trend Analysis of the Phenology Metrics
We computed the mean annual EVI profile. Then, we analyzed the spatiotemporal variation trends of phenology metrics (SOS, EOS, and LOS), from 2001 to 2017, by means of a linear regression method. The linear trends (slope and p-value) are calculated based on pixel-by-pixel between the phenology metrics and years sequence. The slope is calculated as: where slope represents the temporal trend of the changes of land surface phenology, n = 17 is the total years of the monitoring period, and P j (j = 1, 2, and 3) is the phenology metric (SOS, EOS, or LOS).
Negative slopes indicate advanced SOS and EOS and shorter LOS, while positive slopes indicate delayed SOS and EOS and longer LOS. Thirdly, we further investigated how these spatially heterogeneous trends of the altitudinal gradient influenced the phenology metrics (SOS, EOS, and LOS) and tried to examine the differences in phenological responses along with elevation effects.

Determination of the Preseason of Forest Phenology
The temperatures in the preceding months are known to control the timing of spring leaf-out [15] and leaf senescence [21]. To define the length of the preseason whose mean temperature and accumulated precipitation had the most influence on SOS (or EOS), we analyzed the correlation coefficients between SOS (or EOS) and mean temperature and accumulated precipitation during 0, 1, 2, . . . n months before the multiyear averaged SOS (or EOS) using Spearman's rank correlation coefficients. The preseason for mean temperature was defined by the highest absolute correlation coefficient between mean temperature and the multiyear averaged SOS (or EOS) [13,15,21]. We obtained the length of the preseason for precipitation in the same way but replacing the mean temperature with accumulated precipitation. The length of the preseason for SOS (or EOS) was defined as Lsos or (Leos). After that, we applied a temporal partial correlation analysis between the SOS (or EOS) and preseason mean temperature (accumulated precipitation), controlling for corresponding accumulated precipitation (mean temperature) during 2001-2017. The preseason of SOS is generally from the preceding November to the multiyear mean date of SOS [15], and the preseason of EOS is generally from June to the multiyear mean date of EOS [21,63].

Analysis of the Climate Change Impact on Phenology
Partial correlation analyses were used to reveal the relationship between each dependent and independent variable, excluding the confounding effects of other variables [49,64]. In this study, the pixel-based approach was used to analyze the partial correlation of forest phenology metrics and climatic factors (e.g., mean temperature, cumulated precipitation) in the preseason before each phenology metric. Considering the cumulative effect of climatic factors on SOS (or EOS), we obtained the largest significant partial correlation coefficient between the SOS (or EOS) and the temperature in n-month, after controlling the remaining corresponding EOS (or SOS) and precipitation. For instance, if the n-month temperature had the highest correlation with SOS (or EOS), the cumulative effect of temperature on SOS (or EOS) for this pixel would be treated as n months. In other words, at each pixel, the continuous n months before SOS (or EOS) (comprising the SOS (or EOS) appearing month) are the most important hints to trigger SOS (or EOS). We obtained the cumulative effect of precipitation on SOS (or EOS) in the same way but replacing the mean temperature with accumulated precipitation. The partial correlations between SOS (or EOS) and the preseason mean temperatures and cumulative precipitations were assessed at a 95% confidence level during 2001-2017. The statistical significance of the partial correlation coefficient was evaluated by using the Student's t-test.

Spatial Patterns of Phenology Metrics (SOS, EOS, LOS)
The spatial patterns of the forest phenology metrics, which are averaged from multiyear phenology data from 2001 to 2017, in the QMs are shown in Figure 2. With increasing altitude, the SOS, EOS, and LOS were gradually postponed, advanced, and shortened, respectively, which was in conformance with the regional hydrothermal conditions. The SOS was generally between the 81st and 120th Julian days, using a 10% to 90% percentile ( Figure 2a). Owing to the high altitudes in the west of the QMs, the SOS began later, and in the eastern part, the SOS started earlier because of the lower altitude. The central QMs have the largest elevation differences, so the span of the SOS is larger. The EOS showed the opposite trend compared with the SOS and is shown in Figure 2b. Contrary to the SOS, the EOS ends earlier in the west of the QMs because of the higher altitudes. The eastern part of the QMs has relatively lower elevations, with later EOS. The central QMs have the largest elevation differences, so the span of EOS is also larger. The dates for the EOS mainly range from the 260th to 310th Julian days, using a 10% to 90% percentile. The LOS is determined by both the SOS and the EOS. Regions with earlier SOS averaged dates and later EOS averaged dates featured longer LOS results (Figure 2c), which are usually associated with low altitudes. The LOS ranged from 150 to 230 days across the study area, using a 10% to 90% percentile.
The SOS and EOS were getting later and earlier with increasing elevation, respectively; as a consequence, the LOS became shorter with elevation ( Figure 3). The SOS was delayed by 2.4 days per 100 m increase in elevation, whereas the EOS was advanced by 1.4 days, and the LOS was shortened by 3.8 days. The changes in phenology along elevation gradients were statistically significant. However, the EOS remained stable at altitudes above 2200 m above sea level (a.s.l)., and the shortening of LOS slowed down due to the stability of EOS above 2200 m a.s.l. the SOS, the EOS ends earlier in the west of the QMs because of the higher altitudes. The eastern part of the QMs has relatively lower elevations, with later EOS. The central QMs have the largest elevation differences, so the span of EOS is also larger. The dates for the EOS mainly range from the 260th to 310th Julian days, using a 10% to 90% percentile. The LOS is determined by both the SOS and the EOS. Regions with earlier SOS averaged dates and later EOS averaged dates featured longer LOS results (Figure 2c), which are usually associated with low altitudes. The LOS ranged from 150 to 230 days across the study area, using a 10% to 90% percentile. The SOS and EOS were getting later and earlier with increasing elevation, respectively; as a consequence, the LOS became shorter with elevation ( Figure 3). The SOS was delayed by 2.4 days per 100 m increase in elevation, whereas the EOS was advanced by 1.4 days, and the LOS was shortened by 3.8 days. The changes in phenology along elevation gradients were statistically significant. However, the EOS remained stable at altitudes above 2200 m above sea level (a.s.l)., and the shortening of LOS slowed down due to the stability of EOS above 2200 m a.s.l.

Interannual Variability and Trends of Phenology Metrics at Different Altitudes
At the regional scale, the average SOS, EOS, and LOS displayed advanced by 0.13 days year −1 , delayed by 0.22 days year −1 , or extended by 0.35 days year −1 from 2001 to 2017 ( Figure 4). Figure 5 further illustrates the spatial distributions of the interannual variability of phenology metrics in the QMs from 2001 to 2017, which depicts a large spatial difference. The variation trend of SOS is negative 66.9% for the whole forest land in the study area and is concentrated on −0.18 days each year, denoting that most of the regions underwent an earlier SOS (Figure 5a,d). The advanced SOS is mainly located at a middle-high altitude of the QMs (Figure 5b,e). The delayed EOS is focused on 0.17 days each year, and the distribution of EOS is a flat peak. Overall, 72.2% of the entire region presented a lengthened growing season, which is caused by both the earlier SOS and later EOS (Figure 5c,f). The flatness or sharpness between the LOS and EOS is similar, and the spatial pattern of LOS is also similar to that of EOS in most of the area. According to the slopes in Figure 4, it seems EOS contributed 64% in lengthening the growing season. Therefore, a later EOS was the more significant factor in lengthening the growing season. Furthermore, we further probed how the altitude gradient influenced the differences in phenological variation and found that the inter-annual changes of SOS, EOS, and LOS showed different patterns across altitudinal gradients ( Figure 6). The advanced trend in the SOS date showed a slight decrease with increasing elevation. The delayed trend in the EOS date was first increasing and then decreasing, and the peak area was about 1500 m. The changes in the trend of LOS corresponded well with the EOS.  Figure 5 further illustrates the spatial distributions of the interannual variability of phenology metrics in the QMs from 2001 to 2017, which depicts a large spatial difference. The variation trend of SOS is negative 66.9% for the whole forest land in the study area and is concentrated on −0.18 days each year, denoting that most of the regions underwent an earlier SOS (Figure 5a,d). The advanced SOS is mainly located at a middle-high altitude of the QMs (Figure 5b,e). The delayed EOS is focused on 0.17 days each year, and the distribution of EOS is a flat peak. Overall, 72.2% of the entire region presented a lengthened growing season, which is caused by both the earlier SOS and later EOS (Figure 5c,f). The flatness or sharpness between the LOS and EOS is similar, and the spatial pattern of LOS is also similar to that of EOS in most of the area. According to the slopes in Figure 4, it seems EOS contributed 64% in lengthening the growing season. Therefore, a later EOS was the more significant factor in lengthening the growing season. Furthermore, we further probed how the altitude gradient influenced the differences in phenological variation and found that the inter-annual changes of SOS, EOS, and LOS showed different patterns across altitudinal gradients ( Figure 6). The advanced trend in the SOS date showed a slight decrease with increasing elevation. The delayed trend in the EOS date was first increasing and then decreasing, and the peak area was about 1500 m. The changes in the trend of LOS corresponded well with the EOS.        Figure S3 shows that the spatial distribution of the lengths of the SOS and EOS preseason correspond to temperature and precipitation across altitudinal gradients. In terms of the total percentages at each month, we observed that the preseason length of SOS (EOS) to the preseason mean temperature was almost located in month 0 and 1, which is the mean temperature of the current month and the previous month for leaf-out (leaf-off) has the greatest impact of forest leaf-out and leaf-off. However, from the perspective of effects for precipitation, the preseason length of SOS and EOS are evenly distributed throughout the previous six months and have a strong relationship with the different lengths of accumulative precipitation.

Climate Impacts on Forest Phenology
The spatial patterns of the trends of preseason temperature and precipitation for SOS and EOS during 2001-2017 are shown in Figure S4. The preseason temperature for SOS shows an increasing trend in high latitudes in the middle and west of the QMs and in the west of the QMs. The preseason precipitation for SOS shows an increasing trend in the middle and west of the QMs and mainly displays a decreasing trend in the east of the QMs. The preseason temperature for EOS exhibits an increasing trend in the entire QMs, except the southwestern region. The preseason precipitation for EOS shows a non-uniform distribution and displays decreasing trends in the north and west areas.
The patterns of the partial correlation between SOS and the preseason mean temperature exhibited an inconsistent spatial distribution. Figure 7a and Table 1 show that spring SOS in the QMs were negatively correlated with temperature during Lsos across 73.2% of the QMs, and 13.2% areas, which mainly distributed in north of QMs, show significant (p < 0.05) negative correlation with the preseason mean temperature, i.e., increased spring temperature resulted in the advancement of SOS. The spatial distribution of the negative correlations between SOS and Leos-averaged temperatures mainly located in the north and east regions of QMs. Negative and positive partial correlations between preseason precipitation and SOS dominated about 60.5% and 39.5% of the QMs, and about 7.1% and 1.9% of them were significant (p < 0.05), respectively (Figure 7b, Table 1), indicating that increased precipitation resulted in the delayed of SOS, and increasing temperature and precipitation  Figure S3 shows that the spatial distribution of the lengths of the SOS and EOS preseason correspond to temperature and precipitation across altitudinal gradients. In terms of the total percentages at each month, we observed that the preseason length of SOS (EOS) to the preseason mean temperature was almost located in month 0 and 1, which is the mean temperature of the current month and the previous month for leaf-out (leaf-off) has the greatest impact of forest leaf-out and leaf-off. However, from the perspective of effects for precipitation, the preseason length of SOS and EOS are evenly distributed throughout the previous six months and have a strong relationship with the different lengths of accumulative precipitation.

Climate Impacts on Forest Phenology
The spatial patterns of the trends of preseason temperature and precipitation for SOS and EOS during 2001-2017 are shown in Figure S4. The preseason temperature for SOS shows an increasing trend in high latitudes in the middle and west of the QMs and in the west of the QMs. The preseason precipitation for SOS shows an increasing trend in the middle and west of the QMs and mainly displays a decreasing trend in the east of the QMs. The preseason temperature for EOS exhibits an increasing trend in the entire QMs, except the southwestern region. The preseason precipitation for EOS shows a non-uniform distribution and displays decreasing trends in the north and west areas.
The patterns of the partial correlation between SOS and the preseason mean temperature exhibited an inconsistent spatial distribution. Figure 7a and Table 1 show that spring SOS in the QMs were negatively correlated with temperature during Lsos across 73.2% of the QMs, and 13.2% areas, which mainly distributed in north of QMs, show significant (p < 0.05) negative correlation with the preseason mean temperature, i.e., increased spring temperature resulted in the advancement of SOS. The spatial distribution of the negative correlations between SOS and Leos-averaged temperatures mainly located in the north and east regions of QMs. Negative and positive partial correlations between preseason precipitation and SOS dominated about 60.5% and 39.5% of the QMs, and about 7.1% and 1.9% of them were significant (p < 0.05), respectively (Figure 7b, Table 1), indicating that increased precipitation resulted in the delayed of SOS, and increasing temperature and precipitation is the main factor to contribute to the advanced SOS. The positive correlations between preseason precipitation and SOS were largely sited in the middle-west high altitudes of the QMs. Table 1. Proportions of partial correlation coefficients between the SOS (EOS) and mean temperatures and accumulated precipitations during Lsos, i.e., length of the preseason for SOS (or Leos, i.e., length of the preseason for EOS) for significantly positive correlations (R < 0.4974), not significantly positive correlations (R > 0.4974), not significantly negative correlations (R < −0.4974), and significantly negative correlations (R < −0.4974). Values of R = ± 0.4974 correspond to p-values of 0.05, according to the two-tailed Student's t-tests.

Correlation
Lsos Leos is the main factor to contribute to the advanced SOS. The positive correlations between preseason precipitation and SOS were largely sited in the middle-west high altitudes of the QMs.   On the contrary, the partial correlations for EOS mainly exhibited positive correlations with the mean temperature and showed negative correlations with accumulative precipitation during Leos (Figure 7c,d, and Table 1). A positive partial correlation between EOS and the preseason temperature during Leos dominated about 56% of the study area, which was found in middle and west of the QMs, and about 6.6% of the correlations were significant (p < 0.05), denoting that an increased temperature during Leos delayed the EOS. Meanwhile, 44% area, which was positioned in the west of the QMs, show a negative correlation with preseason temperature during Leos, and only 2.4% of the correlations were significant (p < 0.05). We also found a negative partial correlation between accumulative precipitation and EOS in 55.6% of the area, which is located in the mid-west area during Leos, and only 7.2% of the correlations were significant (p < 0.05). The 44.4% area, which is situated in the west and east of the QMs, showed a positive correlation between preseason precipitation and On the contrary, the partial correlations for EOS mainly exhibited positive correlations with the mean temperature and showed negative correlations with accumulative precipitation during Leos (Figure 7c,d, and Table 1). A positive partial correlation between EOS and the preseason temperature during Leos dominated about 56% of the study area, which was found in middle and west of the QMs, and about 6.6% of the correlations were significant (p < 0.05), denoting that an increased temperature during Leos delayed the EOS. Meanwhile, 44% area, which was positioned in the west of the QMs, show a negative correlation with preseason temperature during Leos, and only 2.4% of the correlations were significant (p < 0.05). We also found a negative partial correlation between accumulative precipitation and EOS in 55.6% of the area, which is located in the mid-west area during Leos, and only 7.2% of the correlations were significant (p < 0.05). The 44.4% area, which is situated in the west and east of the QMs, showed a positive correlation between preseason precipitation and EOS, and only the 3.3% area was significant (p < 0.05). These results suggest that increasing temperature and precipitation have a countervailing effect on delayed EOS in the QMs. Figure 8 shows spatial patterns of the sensitivities of the SOS (EOS) to preseason mean temperature and accumulated precipitations in the QMs during 2001-2017. The sensitivity of SOS to preseason temperature displays that an increase of 1 • C in the regionally averaged preseason temperature would advance the average SOS by 1.23 days. The 28.9% and 23.6% region shows the SOS sensitivity to preseason temperature in the north and west areas ranged from −4 to −2 days • C −1 and −2 to −1 days • C −1 (Figure 8a). Across the QMs, the sensitivity of SOS to preseason precipitation is not obvious, and an increase of 1 mm will advance 0.03 days for SOS in the regionally level. There are 57.7% of pixels that show negative precipitation sensitivity values, and the sensitivity is ranged from −0.01 to −0.15 days mm −1 in 48% of the pixels. Meanwhile, the positive sensitivity is ranged from 0 to 0.1 days mm −1 in 35% of the pixels, which occurred mostly in the high altitudes of the QMs. The sensitivity of EOS to preseason temperature shows that an increase of 1 • C in regionally an averaged preseason temperature would delay the average EOS by 0.72 days. There are about 51.3% of the pixels showed positive temperature sensitivity values, the sensitivity is mainly concentrated from 1 to 2 days • C −1 in 15.1% of the pixels and from 0 to 1 days • C −1 in 13% of the pixels, which occurred mostly in the western and middle part of the QMs. There are also 32.5% of pixels that showed negative temperature sensitivity values, and the sensitivity ranged from −2 to −4 days • C −1 in 15.8% and from −4 to −6 days • C −1 in 16.7%, which occurred mostly in the west of the QMs. Across the QMs, the sensitivity of EOS to preseason precipitation is also not obvious, and an increase of 1 mm will advance 0.002 days for EOS at the regional level. Although 52.2% of the pixels showed negative precipitation sensitivity values, which ranged mostly from 0 to −0.1 days mm −1 , and the precipitation sensitivity value more than 0 days mm −1 in 47.8% of the pixels, which occurred mainly from 0 to 0.1 days mm −1 . EOS, and only the 3.3% area was significant (p < 0.05). These results suggest that increasing temperature and precipitation have a countervailing effect on delayed EOS in the QMs. Figure 8 shows spatial patterns of the sensitivities of the SOS (EOS) to preseason mean temperature and accumulated precipitations in the QMs during 2001-2017. The sensitivity of SOS to preseason temperature displays that an increase of 1 °C in the regionally averaged preseason temperature would advance the average SOS by 1.23 days. The 28.9% and 23.6% region shows the SOS sensitivity to preseason temperature in the north and west areas ranged from −4 to −2 days °C −1 and −2 to −1 days °C −1 (Figure 8a). Across the QMs, the sensitivity of SOS to preseason precipitation is not obvious, and an increase of 1 mm will advance 0.03 days for SOS in the regionally level. There are 57.7% of pixels that show negative precipitation sensitivity values, and the sensitivity is ranged from −0.01 to −0.15 days mm −1 in 48% of the pixels. Meanwhile, the positive sensitivity is ranged from 0 to 0.1 days mm −1 in 35% of the pixels, which occurred mostly in the high altitudes of the QMs. The sensitivity of EOS to preseason temperature shows that an increase of 1 °C in regionally an averaged preseason temperature would delay the average EOS by 0.72 days. There are about 51.3% of the pixels showed positive temperature sensitivity values, the sensitivity is mainly concentrated from 1 to 2 days °C −1 in 15.1% of the pixels and from 0 to 1 days °C −1 in 13% of the pixels, which occurred mostly in the western and middle part of the QMs. There are also 32.5% of pixels that showed negative temperature sensitivity values, and the sensitivity ranged from −2 to −4 days °C −1 in 15.8% and from −4 to −6 days °C −1 in 16.7%, which occurred mostly in the west of the QMs. Across the QMs, the sensitivity of EOS to preseason precipitation is also not obvious, and an increase of 1 mm will advance 0.002 days for EOS at the regional level. Although 52.2% of the pixels showed negative precipitation sensitivity values, which ranged mostly from 0 to −0.1 days mm −1 , and the precipitation sensitivity value more than 0 days mm −1 in 47.8% of the pixels, which occurred mainly from 0 to 0.1 days mm −1 .

Variations in Phenology
Understanding the forest phenological changes and variation trends is valuable for reducing the negative impacts of climate change and utilizing its advantages. Recent studies have shown advanced SOS, delayed EOS, and extended LOS at middle and high latitudes of the Northern Hemisphere [10,19,20,30,32,65,66]. Zeng and Jeong et al. [20,30] also found that phenological trends became slower during 2000-2008 compared to 1982-1999 and did not display strong significance. Comparing with the changes in phenology (SOS: −7.9, EOS: 3.7; LOS: 11.6 days/decade) in China's temperate vegetation [32] from 1982 to 1999, our results also show lower changes in phenology (SOS: −1.3, EOS: 2.2; LOS: 3.5 days/decade). Meanwhile, our results also show that the advanced magnitude of SOS is lower than the delayed magnitude of EOS, which is the opposite of what was found by Piao et al. [32]. This observation is similar to the results of Wang et al. [66], which indicates that the satelliteobserved phenology change rates become slowed down during a global warming hiatus between 1998 and 2012.
The change in elevation is a key factor influencing the phenology of vegetation in mountainous environments [67,68]. Although our results indicate that elevation change was a significant factor for

Variations in Phenology
Understanding the forest phenological changes and variation trends is valuable for reducing the negative impacts of climate change and utilizing its advantages. Recent studies have shown advanced SOS, delayed EOS, and extended LOS at middle and high latitudes of the Northern Hemisphere [10,19,20,30,32,65,66]. Zeng and Jeong et al. [20,30] also found that phenological trends became slower during 2000-2008 compared to 1982-1999 and did not display strong significance. Comparing with the changes in phenology (SOS: −7.9, EOS: 3.7; LOS: 11.6 days/decade) in China's temperate vegetation [32] from 1982 to 1999, our results also show lower changes in phenology (SOS: −1.3, EOS: 2.2; LOS: 3.5 days/decade). Meanwhile, our results also show that the advanced magnitude of SOS is lower than the delayed magnitude of EOS, which is the opposite of what was found by Piao et al. [32]. This observation is similar to the results of Wang et al. [66], which indicates that the satellite-observed phenology change rates become slowed down during a global warming hiatus between 1998 and 2012.
The change in elevation is a key factor influencing the phenology of vegetation in mountainous environments [67,68]. Although our results indicate that elevation change was a significant factor for variations in trends of the phenological metrics, the effect was not uniform for all elevations or phenology metrics. First, the SOS was delayed by 2.4 days per 100 m with increasing altitude; the EOS was integrally advanced by 1.4 days per 100 m but stayed stagnant above 2200 m. Meanwhile, the LOS was shortened by 3.8 days per 100 m below 2200 m, and the rate of shortening became slower above 2200 m. Shen et al. [16] found a similar trend that the green-up date stayed stagnant above a certain altitude in the Qinghai-Tibetan Plateau.
Second, as for spatial distribution, the amplitude of the advanced trend for SOS slightly increased, and the amplitude of the delayed trend of EOS reached a maximum at 1500 m. Therefore, the elevation-induced advanced trend of SOS declined slightly with increasing elevation, which was illustrated by Vitasse et al. [17] and Vandvik et al. [69]. Furthermore, while previous studies only reveal an elevation-induced SOS shift [17,69], our results also show the elevation-induced EOS shift and an elevation-induced LOS shift. The declining elevation-induced phenological shift will have important implications for the structure and function of mountain ecosystems, and the actual impact of such ongoing change needs to continue to be studied in the future. In addition, the advancement and delay of SOS and EOS coexisted at different altitudes. The elevation and aspect of mountains can shift over very short distances, altering the local temperature patterns, imposing very significant impacts on the timing of phenological events, and ultimately leading to heterogeneity [16]. Ding et al. [70] showed an obvious differentiation of trends in phenological changes along the gradient. The results of Ding et al. [70] indicate that as the altitude rises, the amplitude of advancement for SOS decreases from 2.4 days year −1 to 0.3 days year −1 , the amplitude of delay for EOS decreases gradually from 1.3 days year −1 to 0 days year −1 , and the amplitude of lengthening for LOS decreases from 3.4 days year −1 to 0.2 days year −1 . The different amplitudes of phenological change could be related to different data sources, temporal scales, and regional differences [38,59,62,70].

Relationships between Phenology Metrics and Climatic Factors
To reveal the variability of the phenology metrics (SOS, EOS, and LOS) in response to climatic factors, partial correlation analyses were carried out with preseason climatic factors (including mean temperature and cumulative precipitation). We found that a 0-2 month accumulated temperature (spring temperature), which accounted for 73.2% of the entire study area, had the greatest impact on the SOS, indicating that warming spring temperatures promoted the advancement of the SOS. This finding agrees with the results reported in the literature [32,38,62]. Zhao et al. [62] analyzed the relationship between temperature and SOS in the forest of northeast China and showed that the spring temperature had the greatest positive impact on the SOS. Tang et al. [38] also reported that the forest SOS is negatively correlated with mean spring temperature in the Greater Khingan Mountain area. Crabbe et al. [71] showed that warmer spring temperatures at high latitudes in Europe lead to a significantly earlier SOS in broad-leaf deciduous forest and mixed forest. Our study found that the duration, in months, of temperatures correlated with EOS increased with altitude. Growing degree days have a negative correlation with chilling accumulation [12]. Therefore, a longer duration of cooling, in days, is required at high altitudes to initiate a delayed EOS. This finding is consistent with other studies [71][72][73]. Crabbe et al. [71] showed that warmer temperatures in extreme warm autumn events lead to significantly delayed EOS for needle leaf evergreen forest, broad-leaf deciduous forests, and mixed forests. We also showed that the cumulative precipitation in different months affected the EOS for the forests. Some studies have indicated that LOS has been prolonged in recent years, a change that is more closely associated with temperature than with precipitation, and that the different patterns between the phenological metrics and climate factors may be associated with their growth environment [74]. Several studies have also proposed that a time-lagged and cumulative effect exists between vegetation growth and climate changes [32,[75][76][77]. Different lagged and cumulated effects between the phenological metrics and climate factors have been noticed and attributed to species, and locations, at diverse temporal and spatial scales [22]. Chmielewski and Rotzer [78] found a significant negative relationship between the SOS and the mean temperature of spring (February-April). Piao et al. [32] identified that the phenological metrics were correlated most significantly with the mean temperature during the previous 2-4 months. Our results indicated that the absolute value of the correlation coefficient between the SOS and temperature was higher in the preceding two months before SOS, implying that the cumulative temperature during the two months before the SOS contributed significantly to the growth of forests. A similar pattern can also be observed for the relationship between temperature and the EOS. The temperature of the preceding two months before the EOS contributed significantly to the growth of forests. Regional differences in temperature sensitivity may be derived from interactions among climatic factors [60]. Previous studies indicate that temperature sensitivity can vary in space and time [13,79]. Our result showed that the mean sensitivity of SOS to preseason temperature is −1.23 days • C −1 , which is lower than the temperate region (−2.5 days • C −1 ) [79]. Symmetric effects of daytime and nighttime warming also have different sensitivity to phenology metrics, for example, an increase of 1 • C Tmin and Tmax in spring would advance SOS 1.34 days and 0.64 days, respectively [80]. In contrast to the extensive research efforts on the sensitivity of SOS, the sensitivity of EOS to temperature is more challenging to understand [22,81]. Although there are about 0.72 days delay with 1 • C increasing, our result showed a different sensitivity of EOS preseason temperature between the eastern and mid-west area. Precipitation has no obvious effects on SOS or EOS at the regional level. A significant interaction between temperature and precipitation indicates the complex mechanisms of phenological responses to climate change [82].

Uncertainty and Future Research Needs
Our results showed that the average SOS was mainly distributed between the 81th and 120th Julian days, and the mean EOS was primarily between the 260th and 310th Julian days, which was in line with the results of previous studies [26,65]. According to the phenology statistics of eight representative woody plants, including Princess tree (Paulownia tomentosa), Babylon weeping willow (Salix babylonica), Dwarf elm (Ulmus pumila), Paper mulberry (Broussonetia papyrifera), Chinese wild peach (Prunus davidiana), and Black locust (Robinia pseudoacacia) in the northern edge of the Qingling Mountains, from 1998 to 2008, He [83] showed that the leaves were completely unfolded on the 96.5th Julian day of the year, and the total discoloration of the leaves was on the 308.8th Julian day of the same year. Li and Liu [84] conducted a statistical survey based on phenological observations from 1986 to 2010 at Xi'an Botanical Garden and Qujiang Scenic Spot. Their results indicated that 37.7% of the plants started leaf-out in March, 61.6% of the plants started leaf-out in April, 67.6% of the plants began to color and defoliate in September, while 32.4% of the plants discolored in October. These conclusions are consistent with the result of our study that 94.2% of trees' SOS in the QMs occurred between the 81th and 120th Julian days, and 98% of trees' EOS appeared between the 270th and 311th Julian days. The spatial extent and time span shown by Wang, Xue, Quan and He [53] is almost consistent; therefore, the spatial distribution of the SOS and EOS showed a similar pattern with our results. Our study demonstrated that using a MODIS EVI time-series dataset to obtain the phenological metrics offered acceptable reliability, and the results were comparable with the ground observation data and can be used to reveal the forest phenological characteristics.
Differences in the models for smoothing the vegetation index and the model parameters resulted in potential differences of uncertainty and bias among alternative phenology analyses [85]. Even with the same model, different thresholds often lead to changes in phenology metrics. The uncertainty of our research results is firstly associated with the reconstruction of remote sensing data by the HANTS algorithm. When the time-series dataset is reconstructed, the acquired images are not as time-critical as the Fast Fourier Transform (FFT) with equal time interval images. Differences in the time intervals of time-series images and in the calibrations of model parameters affect the amounts or degrees of uncertainty and bias in the results [38,86]. Therefore, ground-based phenological observations should be incorporated into the threshold of phenological estimates to improve the accuracy of phenological estimates in the future. This research expanded to an in-depth investigation of the variation in phenological response to climate change. Phenology changes are strongly influenced by a combination of climate factors such as temperature, precipitation, and solar radiation, which typically are significantly correlated with elevation.
The shortening of the growing season at low altitudes of the QMs may be related to mixed pixels containing both forests and cultivated vegetation. This is because the cultivated vegetation in the QMs is mainly planted in lowlands with gentle slopes and in reclamation zones of the low valley with altitudes between 984 m a.s.l. and 1375 m a.s.l. [87]. The increases in temperature and heat resources to the north of the QMs were found to be larger than those in the south [88]. It was also found that in recent years, the boundaries between the semi-arid and humid areas, and between the warm temperate zone and the subtropical zone, have shifted northward [89]. The question of how to quantify the differentiated impact of climatic factors such as temperature, precipitation, and sunshine duration on the phenology of the north side and south side of the QMs should be addressed independently in the follow-up studies.

Conclusions
Based on the time series of the MODIS EVI dataset with 8-day and 500 m resolutions and climate dataset, we estimated pixel-based SOS, EOS, and LOS and analyzed phenological change across altitudinal gradients and the influence of the preseason temperature and precipitation across the QMs in Central China, for the period 2001-2017. The results indicated that the SOS was more sensitive to changes in elevations than the EOS. The SOS, EOS, and LOS were delayed, advanced, and shortened by 2.4, 1.4, and 3.8 days per every 100 m increase in altitude, respectively. The EOS was found to be the most significant factor in prolonging the growing season. The elevation-induced leaf-out shift declined slightly from 2001 to 2017. The derived partial correlation coefficients between the phenological metrics and the climatic factors indicated that the current and previous-month mean temperature had strong impacts on SOS and EOS. Furthermore, the temperature was significantly correlated with variations in SOS and EOS, an increase of 1 • C in the regionally averaged preseason temperature would advance the average SOS by 1.23 days and delay the average EOS by 0.72 days. The precipitation was observed to have no significant correlation with SOS and EOS, the sensitivity of SOS and EOS to the preseason precipitation was not obvious across the QMs. This study provides valuable evidence to improve our understanding of the complex impacts of climatic and elevational changes on the spatiotemporal phenology variabilities and trends in the QMs, an important ecotone in China. Future studies are still needed to investigate the phenology response of different species and boundary effects in the QMs.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/10/11/1007/s1, Figure S1: An example of the determination of the start (end) of the growing season using MODIS EVI with SOS (EOS) extraction method. Figure S2: The threshold values were determined according to the lowest absolute error between the leaf-out and leaf senescence dates measured in situ. Figure S3: The length of the preseason period corresponds to temperature and precipitation. Figure S4: Change trend of the preseason temperature and precipitation for SOS and EOS.