Climate and Spring Phenology Effects on Autumn Phenology in the Greater Khingan Mountains , Northeastern China

Vegetation phenology plays a key role in terrestrial ecosystem nutrient and carbon cycles and is sensitive to global climate change. Compared with spring phenology, which has been well studied, autumn phenology is still poorly understood. In this study, we estimated the date of the end of the growing season (EOS) across the Greater Khingan Mountains, China, from 1982 to 2015 based on the Global Inventory Modeling and Mapping Studies (GIMMS) normalized difference vegetation index third-generation (NDVI3g) dataset. The temporal correlations between EOS and climatic factors (e.g., preseason temperature, preseason precipitation), as well as the correlation between autumn and spring phenology, were investigated using partial correlation analysis. Results showed that more than 94% of the pixels in the Greater Khingan Mountains exhibited a delayed EOS trend, with an average rate of 0.23 days/y. Increased preseason temperature resulted in earlier EOS in most of our study area, except for the semi-arid grassland region in the south, where preseason warming generally delayed EOS. Similarly, EOS in most of the mountain deciduous coniferous forest, forest grassland, and mountain grassland forest regions was earlier associated with increased preseason precipitation, but for the semi-arid grassland region, increased precipitation during the preseason mainly led to delayed EOS. However, the effect of preseason precipitation on EOS in most of the Greater Khingan Mountains was stronger than that of preseason temperature. In addition to the climatic effects on EOS, we also found an influence of spring phenology on EOS. An earlier SOS led to a delayed EOS in most of the study area, while in the southern of mountain deciduous coniferous forest region and northern of semi-arid grassland region, an earlier SOS was often followed by an earlier EOS. These findings suggest that both climatic factors and spring phenology should be incorporated into autumn phenology models in order to improve prediction accuracy under present and future climate change scenarios.


Introduction
Vegetation plays an important role in the land-atmosphere interface [1].Vegetation phenology, which refers to the growth cycle of flora in a particular region, can be seriously affected by climatic changes and is an important component of land surface process models and terrestrial carbon cycle models [2].Since it can serve as an indicator of global climate change, studies on vegetation phenology have become increasingly important.
The third-generation Normalized Difference Vegetation Index (NDVI3g) dataset released by the NASA's Global Inventory Modeling and Mapping Studies (GIMMS) group was the longest NDVI time series to date and has been proved to have good ability for the long-term monitoring of vegetation phenology variability.Additionally, the dataset has been normalized to account for issues such as sensor calibration loss, orbital drift, and atmospheric effects such as volcanic eruptions, so it can provide higher quality data for regions in mid to high latitudes [22].Although the spatial resolution (8 km) of the GIMMS NDVI3g dataset is relatively coarse, for our study area (3.35 × 10 7 ha), it is sufficient to detect long-term changes in vegetation phenology.
Most previously conducted studies of autumn phenology in China have been conducted over short time periods (10~20 years).For example, Chen et al. [23] calculated leaf senescence dates for deciduous species in temperate eastern China from 1982 to 1993, and Piao et al. [24] provided a more detailed analysis on autumn phenology over the entire temperate zone of China from 1982 to 1999.Recent studies, however, have shown that warming and greening trends may have slowed during the first decade of the 21st century compared to those observed during the 1980s and 1990s [25].Therefore, an evaluation of autumn phenology covering the entire period of the past three decades is necessary for providing more detailed insights into the relationship between climatic changes and vegetation phenological variation.
The Greater Khingan Mountains, China, span the mid-temperate and cold temperate zones of mid-high latitudes, the climatic change of which may precede that in low latitudes.Additionally, it is located in a climatic and topographic transition zone between the humid forest biome and the semiarid grassland biome.The transitional nature of the physical conditions makes the vegetation of this region sensitive to climatic change.Therefore, study on the phenology variability in this region can be used as an indicator for the phenology variability of vegetation distributed in regions at low latitudes and is necessary for understanding the response of regional ecosystem to the climatic change in the future.The primary objectives of this study were to (i) explore the EOS trends across different ecogeographical regions in the Greater Khingan Mountains; (ii) investigate the effects of climatic factors (e.g., temperature, precipitation) on the spatial and interannual variation of the EOS; and (iii) reveal the linkage between EOS and SOS.

Study Area
The Greater Khingan Mountains (40 • 59 -53 • 33 N, 115 • 05 -125 • 16 E), China, are located in the eastern portion of the Inner Mongolia Autonomous Region, north of Heilongjiang Province, with a northeast-southwest orientation (Figure 1).The eastern and western slopes of this region are asymmetric; the terrain of the east facing slopes is steep, quickly incising to the Songliao Plain, whereas the terrain of the west facing slopes is relatively gentle, gradually transiting into the Mongolian Plateau.The Greater Khingan Mountains are an important climate demarcation line in China, with a decreasing temperature gradient (mid-temperate to cold temperate) from north to south and a decreasing humidity gradient (humid to semi-arid) from east to west.The summer marine monsoon is blocked by the east slope of the mountains, thus forming a precipitation contrast between the wetter eastern slopes and the arid western slopes.The variety of vegetation types in this region is abundant, with significant differences in vegetation composition among different regions.Cold-temperate coniferous forest cover is prevalent in the high-latitude regions, gradually transitioning to mid-temperate grassland landscapes with decreasing latitude.
The Greater Khingan Mountains are divided into four ecogeographical regions [26] according to the ecogeographical regionalization system of China, which reflects the regional differences of natural elements such as temperature, moisture, and vegetation cover (Figures 1 and 2, Table 1): (i) the mountain deciduous coniferous forest region in the north (MFN); (ii) the forest grassland region in the northwest (FGN); (iii) the centrally located mountain grassland forest (CMF); (iv) and the grassland region in the south (GLS).

GIMMS NDVI3g Dataset
Since it can accurately represent vegetation growth status from different remote sensing information sources, NDVI is a commonly used indicator for monitoring vegetation phenology [28].The third-generation NDVI (NDVI3g) used in this study is the latest and longest production of satellite NDVI records released by the Global Inventory Modeling and Mapping Studies (GIMMS) group [21].It contains over 30 years of observations of 15-day maximum value compositions (MVC) at a spatial resolution of 8 km, and can be acquired from the NASA Ames Ecological Forecasting Lab (http://ecocast.arc.nasa.gov/data/pub/gimms/3g.v0/).Multiple corrections have been applied to eliminate errors and noise caused by satellite sensors, radiation, geometry, atmosphere, and other factors [22].So, compared to the older versions, this dataset can provide higher quality data for regions in mid to high latitudes.In this study, we extracted the GIMMS NDVI3g pixels covering complete 34 years (January 1982-December 2015) to produce the NDVI time-series curve and derive phenological metrics.

Climate Data
Daily mean temperature and cumulative precipitation at 839 meteorological stations distributed throughout China from 1982 to 2015 were acquired from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn).The daily climate data were compiled into monthly data, and then the stations based monthly climate data were interpolated into 8 km × 8 km grids to match the GIMMS NDVI3g dataset using ANUSPLIN 4.3.Finally, the gridded monthly temperature and precipitation data with spatial resolution 8 km × 8 km of the Greater Khingan Mountains from 1982 to 2015 were extracted by mask analysis.

GIMMS NDVI3g Dataset
Since it can accurately represent vegetation growth status from different remote sensing information sources, NDVI is a commonly used indicator for monitoring vegetation phenology [28].The third-generation NDVI (NDVI3g) used in this study is the latest and longest production of satellite NDVI records released by the Global Inventory Modeling and Mapping Studies (GIMMS) group [21].It contains over 30 years of observations of 15-day maximum value compositions (MVC) at a spatial resolution of 8 km, and can be acquired from the NASA Ames Ecological Forecasting Lab (http://ecocast.arc.nasa.gov/data/pub/gimms/3g.v0/).Multiple corrections have been applied to eliminate errors and noise caused by satellite sensors, radiation, geometry, atmosphere, and other factors [22].So, compared to the older versions, this dataset can provide higher quality data for regions in mid to high latitudes.In this study, we extracted the GIMMS NDVI3g pixels covering complete 34 years (January 1982-December 2015) to produce the NDVI time-series curve and derive phenological metrics.

Climate Data
Daily mean temperature and cumulative precipitation at 839 meteorological stations distributed throughout China from 1982 to 2015 were acquired from the China Meteorological Data Sharing Service System (http://cdc.cma.gov.cn).The daily climate data were compiled into monthly data, and then the stations based monthly climate data were interpolated into 8 km × 8 km grids to match the GIMMS NDVI3g dataset using ANUSPLIN 4.3.Finally, the gridded monthly temperature and precipitation data with spatial resolution 8 km × 8 km of the Greater Khingan Mountains from 1982 to 2015 were extracted by mask analysis.
The ANUSPLIN 4.3 algorithm applies a thin-plate spline function for interpolating meteorological data.Using this approach, a number of impact factors are used as covariates for spatially interpolating meteorological elements, greatly improving interpolation precision, and multiple surfaces can be interpolated simultaneously [29].The program is particularly suited to interpolating time series meteorological data and has been widely applied in recent phenological research [30].The measures of interpolation accuracy generated in the log file include the root of the generalized cross-validation (RTGCV), the root of the mean square error (RTMSE), the root of the mean square residual (RTMSR), and signal value.The signal value gives an indication of the degrees of freedom of the fitted spline.The best model judgment criteria [29] are the values of RTGCV, RTMSE, RTMSR, which are minimal, and the signal value, which is less than half the number of meteorological stations.In this study, we used the latitude, longitude, and monthly temperature or precipitation of each meteorological station as independent variables.Elevation, which is known to have a strong influence on climate [31], was a covariate incorporated into the interpolation process.Then, by comparing three test schemes (the number of splines was two, three, and four, respectively) based on the best model judgment criteria, we finally determined the second scheme as optimal to interpolate climate data.Elevations for this study were based on a resampled 8 km grid derived from a 90 m resolution digital elevation model (DEM) from the NASA and National Imagery and Mapping Agency (NIMA) Shuttle Radar Topography Mission (SRTM) (http://glcf.umd.edu/data/srtm/).

SOS and EOS Extraction
There is usually snow cover on the Greater Khingan Mountains from late autumn to the following late spring [32].Snow cover would reduce the NDVI value of the land surface and misrepresent the photosynthetic capacity of vegetation during the non-growing season, which can lead to errors in extraction of phenological parameters [33].Therefore, snow-covered NDVI during non-growing season should be eliminated first.According to the preprocessing method that has been validated in precious studies [20], we replaced the NDVI values of the pixels covered by snow with the median value of the snow-free winter NDVI values between November and the following March.Subsequently, we smoothed and reconstructed the NDVI time series by removing abnormal (extreme high/low) values, which were present due to atmospheric interference and sensor errors.Three time series reconstruction algorithms in the TIMESAT software [34] were tested to smooth and reconstruct the NDVI time series.These were Savitzky-Golay filtering (S-G), the double logistic function (DL), and the asymmetrical Gaussian function (AG).Since each of these aforementioned methods presents unique advantages and disadvantages, each was tested for the study area in order to choose the most appropriate algorithm.We randomly selected the samples in the study area, and the NDVI time series were fitted by the S-G, DL, and AG methods.
A comparison of the original and fitted curves of the GIMMS NDVI3g data at 15-day intervals (Figure 3) indicated that the S-G filtering approach yielded the closest match to the original data.However, the fitted curve of S-G had great volatility and was not smooth enough compared to the other two methods.A smoothed curve is more desirable as it prevents NDVI values from varying directionally over short intervals, and is thus more representative of actual annual vegetation growth.The DL and AG methods both produced smoothed curves with considerable noise reduction, while the performance of the two methods was similar and it was difficult to judge which was better.Therefore, we used three statistical indicators to evaluate the performance of DL and AG methods: (1) the root mean square error (RMSE); (2) Akaike's Information Criterion (AIC); and (3) Bayesian Information Criterion (BIC).The number of free parameters of each method was the same as Atkinson et al. [35].For the DL and AG methods, six and seven parameters were needed, respectively.The RMSE and AIC of AG method were both smaller than the values of DL method (Table 2).Thus, the AG method performed better than DL method.In addition, the AG function can generate a reconstructed NDVI curve that better describes subtle changes in the NDVI sequence data [36], more accurately reflects the characteristics of vegetation growth, and reduces the influence of outside interference factors (e.g., fire and pests).Therefore, the AG method was chosen as the time series reconstruction function to be used in this study.After reconstructing the NDVI time series, the next step is to extract specific phenological parameters such as SOS and EOS.The most common methods for this step include the maximum slope method, the inflection-based method, and the dynamic threshold method.The maximum slope method defined SOS or EOS as the day of year (DOY) when NDVI begins to either rapidly increase (SOS) or decrease (EOS) [37], identified based on the maximum absolute slope of the fitted NDVI curve during growth or senescence [38].The inflection-based method is an algorithm detecting points where maximum curvature occurs in a fitted NDVI time series [39].SOS was defined as a valley point at the beginning of a growing cycle, and EOS was defined as a valley point at the decaying end of a phenology cycle [39].The dynamic threshold method defines SOS and EOS as the DOY when the NDVIratio reaches a predefined threshold of the vegetation growth amplitude during the NDVI rising stage and decline stage, respectively [38].However, the setting of thresholds is based on the characteristics of the NDVI curve, and these characteristics vary among vegetation types [40].In addition, the extent of the study area is large and the type of vegetation is abundant.Therefore, in this study, different NDVI thresholds (e.g., 20%, 30%, 35%, 40%, 45%, or 50%) of multiple sampling  After reconstructing the NDVI time series, the next step is to extract specific phenological parameters such as SOS and EOS.The most common methods for this step include the maximum slope method, the inflection-based method, and the dynamic threshold method.The maximum slope method defined SOS or EOS as the day of year (DOY) when NDVI begins to either rapidly increase (SOS) or decrease (EOS) [37], identified based on the maximum absolute slope of the fitted NDVI curve during growth or senescence [38].The inflection-based method is an algorithm detecting points where maximum curvature occurs in a fitted NDVI time series [39].SOS was defined as a valley point at the beginning of a growing cycle, and EOS was defined as a valley point at the decaying end of a phenology cycle [39].The dynamic threshold method defines SOS and EOS as the DOY when the NDVI ratio reaches a predefined threshold of the vegetation growth amplitude during the NDVI rising stage and decline stage, respectively [38].However, the setting of thresholds is based on the characteristics of the NDVI curve, and these characteristics vary among vegetation types [40].In addition, the extent of the study area is large and the type of vegetation is abundant.Therefore, in this study, different NDVI thresholds (e.g., 20%, 30%, 35%, 40%, 45%, or 50%) of multiple sampling sites in the study area were repeatedly tested.By comparison, we found that the extracted phenological parameters at the threshold of 30% were most consistent with the results of previous studies in similar areas [41,42].Finally, we applied 30% as the dynamic threshold, which was also widely used for satellite-based phenology detection in previous studies [41,43] to extract EOS.The NDVI ratio is defined as: where NDVI t is the NDVI value at time t; NDVI max represents the annual maximum NDVI value; for SOS, NDVI min represents the annual minimum NDVI value during the growth period; and for EOS, NDVI min represents the annual minimum NDVI value during the senescence period.We applied the above three methods to extract the multiyear averaged EOS in the Greater Khingan Mountains.By comparison, we found the spatial patterns of the average EOS varied with different methods (Figure 4).The average EOS based on maximum slope method was earliest, ranging from Julian day 252 to 291.The average EOS based on dynamic threshold method occurred between Julian days 284 and 331, while the average EOS based on the inflection-based method was the latest, ranging from Julian day 288 to 354.The above results were similar to those of previous studies.You et al. [40] founded that the maximum slope method tended to estimate the EOS date earlier by comparing with the threshold and curvature change rate methods; Shang et al. [44] revealed that the inflection-based approach would overestimate the EOS when the vegetation was affected by disturbances, whereas the dynamic threshold approach can estimate it correctly.Additionally, the slope of vegetation NDVI curve was easily affected by external conditions, so the vegetation growth period could not be determined correctly based on the slope, which rendered the maximum slope method unsuitable for the extraction of phenological parameters of long time series [43].Studies also showed that if some vegetation growth trajectories cannot be fitted well with a logistic function, it will prevent an inflection point from being determined correctly [45].The vegetation growth curve of this study was fitted well by an AG function rather than a DL function, so the inflection-based method is not suitable for extracting the phenological parameters in our study area.Comparative studies [46] concluded that the dynamic threshold method is one of the simplest and most effective methods to extract phenological parameters, as it generally keeps dates within a reasonable range based on the threshold conditions and could achieve relatively high accuracy [40].Moreover, compared with the other two methods, the range of the multiyear averaged EOS extracted from the threshold method was the most reasonable referring to the results of previous studies in similar areas.Thus, we finally used the dynamic threshold method to estimate the SOS and EOS over the Greater Khingan Mountains.

Investigating Trends in EOS, SOS, and Climatic Factors
We used Sen's slope estimator to investigate temporal trends in EOS, SOS, and preseason climatic factors across the Greater Khingan Mountains from 1982 to 2015 at the pixel level, along with the Mann-Kendall trend test to determine their significance at each pixel.Compared with the linear trend based on the least squares method, Sen's slope estimator can help diminish the influence of missing time series observations and non-normally distributed data on the analysis results, as well as reduce the interference of abnormal values in a time series [47].Thus, this method is often used to

Investigating Trends in EOS, SOS, and Climatic Factors
We used Sen's slope estimator to investigate temporal trends in EOS, SOS, and preseason climatic factors across the Greater Khingan Mountains from 1982 to 2015 at the pixel level, along with the Mann-Kendall trend test to determine their significance at each pixel.Compared with the linear trend based on the least squares method, Sen's slope estimator can help diminish the influence of missing time series observations and non-normally distributed data on the analysis results, as well as reduce the interference of abnormal values in a time series [47].Thus, this method is often used to analyze long-term sequence datasets to detect the magnitude of the trend: where Q is the Sen's slope; x j and x i represent the sequence value at time j and i, respectively; and median is calculated from all pairs of observations in the time series.The Mann-Kendall trend test is one of the most common methods used for testing time series trends [48].It does not require samples to follow a certain distribution, is not subject to the interference of a few extreme values, is suitable for hydrological, meteorological, and other non-normally distributed data, and has been widely used in the analysis of remotely-sensed data [47].The combined use of Sen's slope estimator and the Mann-Kendall trend test has become an important method for analyzing time series vegetation data [49].

Detecting the Correlation between EOS and SOS as Well as Climatic Factors
During multivariate correlation analysis, the influence of each variable is mutual.Therefore, the correlation between variables will not be truly reflected when correlation analysis is conducted using only two variables [15].So, when analyzing the correlation between the two variables, the influence of other variables must be considered.To address this problem, we applied a temporal partial correlation analysis between EOS and SOS, and the preseason climatic factors.This approach was aimed at statistically investigating the relationship between EOS and a single driving factor while controlling the influence of the two remaining factors.This method has been commonly applied in climate change and ecological studies [3].The partial correlation coefficients were also calculated for each ecogeographical region as well as the entire study area.

Determination of the Preseason
The preseason was defined as the period (with a one-month time-step) before the EOS date for which the correlation between EOS and climatic factors (e.g., temperature and precipitation) was highest during 1982-2015.Following previous studies, the maximum range of this period was set from May to the multiyear average date of EOS [50].The average onset date of vegetation dormancy in the Greater Khingan Mountains over all pixels and years is around late October (as shown in Section 3), so the preseason in this study was defined as May through October.

Spatial Pattern of Autumn Phenology
The average EOS dates during the study period  ranged from Julian day 284 to 331, with an average date of Julian day 292 (late October) for the study area (Figure 5).EOS dates at more than 90% of pixels occurred in October, from Julian day 285 through 305.Both the earliest and latest EOS dates were observed in the MFN region, with a difference of 47 days between them.The CMF region had the earliest EOS date, which occurred on Julian day 290 on average.The average EOS dates for the GLS and FGN regions were Julian days 291 and 292, respectively, slightly later than that observed in the CMF region.
with an average date of Julian day 292 (late October) for the study area (Figure 5).EOS dates at more than 90% of pixels occurred in October, from Julian day 285 through 305.Both the earliest and latest EOS dates were observed in the MFN region, with a difference of 47 days between them.The CMF region had the earliest EOS date, which occurred on Julian day 290 on average.The average EOS dates for the GLS and FGN regions were Julian days 291 and 292, respectively, slightly later than that observed in the CMF region.

Trends in Autumn Phenology and Climatic Factors
The trends of the EOS, SOS, preseason temperature, and preseason precipitation were calculated at the 95% confidence level (Figure 6).The mean EOS date in the whole Greater Khingan Mountains had been delayed by 7.82 days over the past 34 years (an average rate of 0.23 days/y).More than 94% of the study area exhibited delayed trends of EOS, with roughly 72% showing significant delays (p < 0.05) (Figure 6a).The other ~6% of pixels showed earlier EOS, mainly in the southeastern of the study area, but only 12% were significant (p < 0.05).The mean EOS date of all the ecogeographical regions uniformly exhibited delayed trends.In the semi-arid GLS region, the mean EOS date exhibited a minimum delay of 3.06 days (0.09 days/y), while the mean EOS date of FGN region exhibited a maximum delay of 11.56 days (0.34 days/y) (Figure 7a).
of the study area exhibited delayed trends of EOS, with roughly 72% showing significant delays (p < 0.05) (Figure 6a).The other ~6% of pixels showed earlier EOS, mainly in the southeastern of the study area, but only 12% were significant (p < 0.05).The mean EOS date of all the ecogeographical regions uniformly exhibited delayed trends.In the semi-arid GLS region, the mean EOS date exhibited a minimum delay of 3.06 days (0.09 days/y), while the mean EOS date of FGN region exhibited a maximum delay of 11.56 days (0.34 days/y) (Figure 7a).During 1982-2015, a preseason temperature increase was observed in more than 88% of the study area (around 61% was statistically significant at p < 0.05), and occurred mainly in the FGN and CMF regions (Figure 6c).The preseason temperature in all of the ecogeographical regions (MFN, FGN, CMF, and GLS) exhibited a clear increasing trend (+0.02 • C/y, +0.05 • C/y, +0.04 • C/y, and +0.03 • C/y, respectively) (Figure 7c).For the spatial distribution of preseason precipitation variation, neither negative nor positive trends (both nearly 50%) dominated the entire study area (Figure 6d).In addition, the variation trends of preseason precipitation in all the ecogeographical regions were different.In the MFN and FGN regions, the preseason precipitation exhibited decreasing trends (−0.33mm/y and −2.6 mm/y, respectively), while in the CMF and GLS regions, the preseason precipitation exhibited increasing trends (+0.29 mm/y and +0.15 mm/y, respectively) (Figure 7d).
°C/y, respectively) (Figure 7c).For the spatial distribution of preseason precipitation variation, neither negative nor positive trends (both nearly 50%) dominated the entire study area (Figure 6d).In addition, the variation trends of preseason precipitation in all the ecogeographical regions were different.In the MFN and FGN regions, the preseason precipitation exhibited decreasing trends (−0.33mm/y and −2.6 mm/y, respectively), while in the CMF and GLS regions, the preseason precipitation exhibited increasing trends (+0.29 mm/y and +0.15 mm/y, respectively) (Figure 7d).

Effects of Climatic Factors on Autumn Phenology
The partial correlation coefficients between EOS and preseason temperature, preseason precipitation, SOS were also calculated at 95% confidence level (Figure 8).Over 61% of the entire study area showed negative partial correlations between preseason temperature and EOS (about 6% of the relationships were significant at p < 0.05), the remaining ~38% of the study area showed positive partial correlations (Figure 8a).In the humid and semi-humid regions (MFN, FGN, and CMF), more than 60% of each presented negative partial correlations between preseason temperature and EOS (Figure 9), suggesting that higher temperature during the preseason led to an earlier EOS in most of these areas.In most of the semi-arid GLS region, however, increased preseason temperature played a role in delaying the EOS date (about 72% of this region exhibited positive partial correlations) (Figure 9).

Effects of Climatic Factors on Autumn Phenology
The partial correlation coefficients between EOS and preseason temperature, preseason precipitation, SOS were also calculated at 95% confidence level (Figure 8).Over 61% of the entire study area showed negative partial correlations between preseason temperature and EOS (about 6% of the relationships were significant at p < 0.05), the remaining ~38% of the study area showed positive partial correlations (Figure 8a).In the humid and semi-humid regions (MFN, FGN, and CMF), more than 60% of each presented negative partial correlations between preseason temperature and EOS (Figure 9), suggesting that higher temperature during the preseason led to an earlier EOS in most of these areas.In most of the semi-arid GLS region, however, increased preseason temperature played a role in delaying the EOS date (about 72% of this region exhibited positive partial correlations) (Figure 9).Negative partial correlations between preseason precipitation and EOS dominated about 71% of the study area, and about 34% of them were significant (p < 0.05), indicating that increased precipitation during preseason would advance the EOS date (Figure 8b).Interestingly, in each of all the four ecogeographical regions, the relationships between preseason precipitation and EOS were same as that between preseason temperature and EOS (Figure 9).In the humid and semi-humid regions (MFN, FGN, and CMF), most of each region (94%, 98%, and 78%, respectively) experienced negative partial correlations between preseason precipitation and EOS.While in 74% of semi-arid GLS region, the partial correlations were positive, suggesting that increasing precipitation during the preseason could contribute to the extension of EOS.

Relationship between Spring Phenology and Autumn Phenology
Contrary to the delayed trend observed for EOS change, the mean date of the SOS in the Greater Khingan Mountains advanced an average of 0.20 days/y from 1982 to 2015.More than 94% of the pixels exhibited earlier trends of SOS, with more than 65% significant (p < 0.05) (Figure 6b).In addition, in each of the four ecogeographical regions, the mean SOS date all exhibited an earlier trend over the past 34 years (Figure 7b).
Negative partial correlations between EOS and SOS dominated the whole study area (more than 85%), with more than 22% of pixels having statistically significant relationships between the two variables (p < 0.05).Positive partial correlations were observed primarily in the southern MFN and northern GLS regions (Figure 8c).In each of the four ecogeographical regions, SOS was negatively associated with EOS.This was especially true in the FGN region, over 93% of which was absolutely dominated by negative relationships, with more than 45% statistically significant (p < 0.05) (Figure 9).Negative partial correlations between preseason precipitation and EOS dominated about 71% of the study area, and about 34% of them were significant (p < 0.05), indicating that increased precipitation during preseason would advance the EOS date (Figure 8b).Interestingly, in each of all the four ecogeographical regions, the relationships between preseason precipitation and EOS were same as that between preseason temperature and EOS (Figure 9).In the humid and semi-humid regions (MFN, FGN, and CMF), most of each region (94%, 98%, and 78%, respectively) experienced negative partial correlations between preseason precipitation and EOS.While in 74% of semi-arid GLS region, the partial correlations were positive, suggesting that increasing precipitation during the preseason could contribute to the extension of EOS.

Relationship between Spring Phenology and Autumn Phenology
Contrary to the delayed trend observed for EOS change, the mean date of the SOS in the Greater Khingan Mountains advanced an average of 0.20 days/y from 1982 to 2015.More than 94% of the pixels exhibited earlier trends of SOS, with more than 65% significant (p < 0.05) (Figure 6b).In addition, in each of the four ecogeographical regions, the mean SOS date all exhibited an earlier trend over the past 34 years (Figure 7b).
Negative partial correlations between EOS and SOS dominated the whole study area (more than 85%), with more than 22% of pixels having statistically significant relationships between the two variables (p < 0.05).Positive partial correlations were observed primarily in the southern MFN and northern GLS regions (Figure 8c).In each of the four ecogeographical regions, SOS was negatively associated with EOS.This was especially true in the FGN region, over 93% of which was absolutely dominated by negative relationships, with more than 45% statistically significant (p < 0.05) (Figure 9).

Spatial Pattern of the Dominant Factors Affecting Autumn Phenology
Preseason precipitation was the main driver of interannual EOS variations for 26% of all pixels, followed by SOS (17%), and preseason temperature (3%) (Figure 8d).The remaining pixels were not statistically significant at the p < 0.05 level.Preseason precipitation was the primary determinant of the EOS changes in the southwestern MFN, northeastern FGN, and western CMF regions.Preseason temperature sporadically controlled the EOS changes of a small number of pixels, however, indicating that precipitation was the main climatic driver of EOS changes in most of the Greater Khingan Mountains during the preseason.Additionally, SOS showed a certain driving effect on the EOS changes.

Relationship between Autumn Phenology and Climatic Factors
Our results indicated that the average EOS date in the Greater Khingan Mountains, northeastern China during 1982-2015 had a delayed trend of 0.23 days/y, generally consistent with the results of previous studies also showing a delayed EOS trend in the Northern Hemisphere [9], Europe [51], North America [52], and temperate China [50].Zhao et al. [41] estimated the EOS across the entire Northeast China from 1982 to 2013 and found a delay of 0.13 days/y.Tang et al. [42] reported EOS changes in the Hulunber of northeastern China during 1982-2012, with an average delay of 0.29 days/y.Compared to the two above studies investigating the similar study area (northeastern China) and using the same data source (GIMMS NDVI3g), the delayed rate of EOS observed in this study was slightly different from the results of these studies, possibly because the specific location of study area and the phenology extraction method were different.However, the comparison showed that the results of our study could explain well the change of autumn phenology in the Greater Khingan Mountains.
The terrain of the Greater Khingan Mountains is complex and climatic conditions vary across the region.Thus, the response of the EOS to climatic changes varied by latitude.We found that preseason temperature was generally negatively correlated with the EOS date in the MFN, FGN, and CMF regions.For the MFN and CMF regions, which are mainly covered in deciduous trees, this finding may be due to warming-induced enhancement of tree canopy transpiration, in which the roots of trees lacked moisture, thus reducing water utilization [53] and subsequently leading to earlier EOS.In addition, because permafrost is present in portions of these two regions, trees growth, especially the formation and development of Larix gmelini forest, was closely related to the alpine and permafrost environment [54].Permafrost could provide the needed moisture for tree growth, whereas high temperature would decrease the soil freezing depth and the protection of the aquifers, weaken the soil water retention, and thus in the dry season prevent adequate water for trees growth, accelerating leaf senescence [55].Another reason may be that leaves of deciduous trees need a certain accumulated temperature from growth to senescence, and preseason temperature increase makes the tree leaves obtain the required accumulated temperature in a short period of time, thus leading to earlier wilting dates [56].For the FGN region, which is dominated by herbaceous plants (Stipa baicalensis meadow steppe, Filifolium sibiricum meadow steppe), warmer temperatures during the preseason would also accelerate the end of the growing season.This may be because the preseason temperature in this region exhibited the most significant increasing trend (Figure 7c).Thus, an excessive increase in preseason temperature will limit vegetation growth in this region [57].Similarly, preseason precipitation was also negatively correlated with the EOS in most of the MFN, FGN, and CMF regions.The reason may be that the MFN, FGN, and CMF regions are located in humid and semi-humid areas with rain and heat over the same period and abundant rainfall.Soil moisture will rapidly increase with the increasing precipitation during the preseason, and then enhance the vegetation photosynthesis by affecting the vegetation carboxylation.The vegetation will accelerate growth and complete the entire growing season ahead of schedule [58].Additionally, due to the impact of permafrost, high soil moisture in these colder areas may limit the availability of nutrients for vegetation growth [59], thus prematurely ending vegetation growth.
In most of the grassland semi-arid region (GLS), which is mainly covered in herbaceous plants, a warming preseason temperature delayed the EOS date, which was consistent with findings in previous studies based on satellite data and field experiments [60].This result is probably because warming in summer and autumn can enhance the activities of photosynthetic enzymes [61], reduce the rate of chlorophyll degradation during leaf senescence in autumn, and thus delay the time of leaf senescence [62].Furthermore, herbaceous plants are more sensitive to chilling injury compared to deciduous trees, so the warmer preseason may increase the number of days available for vegetation growth and photosynthesis [63], reduce the probability of frost damage [64], and thereby postpone the autumn phenology of herbaceous plants in the GLS region.Positive effects of preseason precipitation on the EOS date of GLS region were also revealed in this study.Possibly due to shallow roots, herbaceous plants in semi-arid region are often subjected to water stress, and their growth is bound to respond rapidly to water.Increasing precipitation during the preseason can effectively alleviate water stress, thus delaying the EOS date and prolonging the growing season [65].In summary, preseason precipitation had a greater impact than temperature on the EOS in most of the Greater Khingan Mountains.Climate change has a direct impact on vegetation phenology and also indirectly affects ecological processes such as species competitive balance, which may alter forest composition (e.g., proportions of coniferous vs. deciduous), which ultimately affects phenological change [66].Previous studies [67] showed that temperature increase may result in the area of coniferous species contracted and broadleaf species expanded, with potential to replace coniferous forest completely.However, tree species migration usually lags behind warming for decades or even a century [68].Therefore, we believe over the past 34 years the impact of climate change on ecological processes may not greatly change forest composition and consequently may not alter our results on the phenological change in the Greater Khingan Mountains.Overall, the relationship between autumn phenology of different vegetation and climatic factors is complex and the impact mechanisms affecting delayed autumn phenology need to be further explored.

The Influence of Spring Phenology on Autumn Phenology
Our results showed that there was a strong positive correlation between SOS and EOS at the pixel level in the southwestern MFN and the northeastern GLS regions, suggesting that an earlier spring phenology would lead to an earlier autumn phenology.Many mechanistic studies have been conducted to explain this phenomenon.First, as the influencing factor of leaf traits, programmed cell death [69] and leaf longevity [70] were reported to constrain the timing of leaf senescence.Second, earlier springs may result in soil water losses via increases in snow sublimation and evapotranspiration in the early part of the growing season; thus, summer and autumn drought duration may increase, leading to earlier EOS dates [71].Third, earlier spring budburst can increase the risk of vegetation suffering damage such as spring frost [72] and insect disease [73], thereby inducing earlier leaf senescence.Finally, the size of the vegetation carbon sink was also found to restrict the correlation between spring and autumn phenology because the accumulation of unstructured carbohydrates in earlier spring may contribute to an earlier peak in autumn carbon content [15].Negative partial correlations between the SOS and EOS also occurred in the Greater Khingan Mountains; specifically, an earlier SOS was followed by a delayed EOS.This special situation may not be explained by the internal mechanism of vegetation alone.The SOS is also significantly affected by climatic factors, showing an earlier trend in other regions under the context of climate warming [4,28], which is consistent with our results.

Uncertainty
In general, there are six methods to reconstruct NDVI time series [18,35].These are (1) Savitzky-Golay filter (S-G); (2) double logistic function (DL); (3) the asymmetric Gaussian function (AG); (4) the harmonic analysis of time series (HANTS); (5) the Fourier transform (FT); and (6) the innovative Whittaker smoother (WS).Each method has advantages and disadvantages, and there is no universal method suitable for all vegetation types of specific study areas [35].The choice of the method depends on the purpose of the study.TIMESAT is a free software package for processing satellite time series data.The software can adapt to the upper envelope of the data, taking missing data and quality flags into account, and has been widely used in a large number of applied studies for data smoothing and phenology parameter extraction [34].TIMESAT only has three models (S-G, DL, AG), so we only tested the S-G, DL, and AG models in this study.The parameter settings of S-G might be off the norm, and we adopted the settings from Zhao et al. [41] (an adaptation strength of 2.0, no spike filtering, a seasonal parameter of 0.5, a seasonal parameter of 0.5, a window size of 2).In this study, we did not take the other common smooth methods (e.g., HANTS, FT, WS) into account, especially the innovative Whittaker smoother, which fits a discrete series to discrete data and penalizes the roughness of the smooth curve, balances reliability of the data, and roughness of the fitted data and may produce a better resulting curve [74].Thus, the finally selected AG method in this study may not be optimal, and a comprehensive comparison of these methods needs to be conducted.
The spatial resolution of GIMMS NDVI3g dataset used in our study was 8 km.Coarse spatial resolution (8 km) that caused mixed vegetation types was the main source of uncertainties.Each pixel value is the mean reflectance of several land cover types over an area of 64 km 2 .Consequently, the processed value approximated the phenology of the dominant vegetation type and overlooked the minor cover types, leading to uncertainties in capturing the phenology of all vegetation types [41].However, our study was based on a regional scale, and researchers have used the NDVI3g dataset with 8 km spatial resolution for phenology study at regional to continental scales [9,42,50], which provided effective information on seasonal variation of vegetation phenology over a long time period and with continuous spatial coverage.Moreover, higher resolutions may not necessarily lead to more accurate results [75].Studies showed that fusing multi-source remote sensing data with different spatial resolutions is considered a feasible way to reduce this uncertainty [76].
The curve of long time series remote sensing data (e.g., NDVI) reflected vegetation dynamics.When the vegetation is affected by disturbances (e.g., wind, pests, fire, and human activities) [77], aberrations would be found in the curve.Disturbances affect phenology in a stochastic manner, leading to another source of uncertainties.When large-scale phenological parameters are extracted, it is difficult to remove the influences of all disturbances during all the study years, thus the results of this study could also potentially be affected.
In addition, solar radiation [9], atmospheric CO 2 concentrations [78], nutrient availability [6], and other factors could also affect autumn phenology.Especially solar radiation, as another important climate factor and a combination of both photoperiod and solar intensity, would trigger autumn phenology in mid-high latitudes [9].Specifically, increased insolation can retard the accumulation of abscisic acid and subsequently slow down leaf senescence [79].Furthermore, the chlorophyll levels of leaves will increase due to the enhanced photosynthetic capacity and CO 2 sequestration caused by the increased insolation, and result in a delayed EOS date [80].However, it is difficult to take all of these factors into consideration in a large-scale remote sensing application.The phenological remote sensing models used in the present studies only make use of time series remote sensing data (e.g., NDVI, EVI) to extract phenological parameters without incorporating other potential phenological restriction factors [9,41].In addition, we only analyzed the effects of temperature and precipitation on autumn phenology in this study.Therefore, we should consider these factors as much as possible in follow-up studies to reduce uncertainties in the results.

Conclusions
In this study, we extracted phenological parameters for the Greater Khingan Mountains from 1982 to 2015 based on the GIMMS NDVI3g dataset and an asymmetric Gaussian fitting function.We analyzed the interannual variation and the influence of major environmental controls on the EOS.Results showed that the EOS in the entire Greater Khingan Mountains and each of its ecogeographical regions all exhibited a delayed trend over the past 34 years.Furthermore, the driving factor of the EOS varied among the different ecogeographical regions (MFN, FGN, CMF, and GLS regions).Both increasing temperature and precipitation during the preseason (May through October) accelerated leaf senescence in autumn at mid-high latitudes (MFN, FGN, and CMF regions), with the exception of the semi-arid region (GLS), where the warming temperature and increased precipitation in summer/autumn delayed the EOS.Importantly, the preseason precipitation was a stronger driver of EOS changes relative to preseason temperature.Apart from climatic factors, we found that the SOS also had effects on the EOS, with mainly negative correlations across most of the study area.However, in the southwestern MFN and northeastern GLS regions, an earlier SOS resulted in an earlier EOS.In addition, compared with climatic factors, SOS played a more significant role than preseason temperature.The results of our study provide a useful reference for understanding the interannual variation of autumn phenology in mid-high latitudes and its response to climate change.In addition, given the important role of autumn phenology in regulating the global carbon budget, our study suggests that both climatic factors and SOS should be assembled into the EOS models to more accurately simulate changes in autumn phenology.

21 Figure 1 .
Figure 1.Location of the study area.

Figure 1 .
Figure 1.Location of the study area.

Figure 3 .
Figure 3.Time series reconstruction for NDVI using three functions.Displayed SOS and EOS estimates were calculated by combining the asymmetrical Gaussian function with dynamic threshold method.

Figure 3 .
Figure 3.Time series reconstruction for NDVI using three functions.Displayed SOS and EOS estimates were calculated by combining the asymmetrical Gaussian function with dynamic threshold method.

21 Figure 4 .
Figure 4. Spatial patterns of multiyear averaged EOS derived from three methods: maximum slope method (a); dynamic threshold method (b); and inflection=based method (c).

Figure 4 .
Figure 4. Spatial patterns of multiyear averaged EOS derived from three methods: maximum slope method (a); dynamic threshold method (b); and inflection=based method (c).

Figure 5 .
Figure 5. Spatial distribution of EOS in the Greater Khingan Mountains during the period 1982-2015.Figure 5. Spatial distribution of EOS in the Greater Khingan Mountains during the period 1982-2015.

Figure 5 .
Figure 5. Spatial distribution of EOS in the Greater Khingan Mountains during the period 1982-2015.Figure 5. Spatial distribution of EOS in the Greater Khingan Mountains during the period 1982-2015.

Figure 6 .
Figure 6.Spatial pattern of interannual temporal trends in EOS (a); SOS (b); preseason temperature (c); and preseason precipitation (d) in the Greater Khingan Mountains over 1982-2015.The pink pixels in the top-left inset indicate that the detected trends were significant at p < 0.05.

Figure 6 .
Figure 6.Spatial pattern of interannual temporal trends in EOS (a); SOS (b); preseason temperature (c); and preseason precipitation (d) in the Greater Khingan Mountains over 1982-2015.The pink pixels in the top-left inset indicate that the detected trends were significant at p < 0.05.

Figure 7 .
Figure 7. Average trend and standard deviation of EOS (a); SOS (b); preseason temperature (c); and preseason precipitation (d) across four ecogeographical regions in the Greater Khingan Mountains from 1982 to 2015.

Figure 7 .
Figure 7. Average trend and standard deviation of EOS (a); SOS (b); preseason temperature (c); and preseason precipitation (d) across four ecogeographical regions in the Greater Khingan Mountains from 1982 to 2015.

Figure 8 .
Figure 8. Spatial distribution of partial correlation coefficients between EOS and preseason temperature (a); preseason precipitation (b); SOS (c) in the Greater Khingan Mountains during 1982-2015.The pink pixels in the top-left inset indicate a significant correlations at p < 0.05.(d) The spatial distribution of major controls on the EOS over the Greater Khingan Mountains during 1982-2015.The grey pixels indicated none of the three factors was significantly (p < 0.05 level) correlated with EOS.

Figure 8 .
Figure 8. Spatial distribution of partial correlation coefficients between EOS and preseason temperature (a); preseason precipitation (b); SOS (c) in the Greater Khingan Mountains during 1982-2015.The pink pixels in the top-left inset indicate a significant correlations at p < 0.05.(d) The spatial distribution of major controls on the EOS over the Greater Khingan Mountains during 1982-2015.The grey pixels indicated none of the three factors was significantly (p < 0.05 level) correlated with EOS.

Figure 9 .
Figure 9. Percentages of the partial correlation coefficients between EOS and preseason temperature, precipitation, and SOS of each ecogeographical region in the Greater Khingan Mountains from 1982 to 2015.Bars above 0 and below 0 represent the percentage of positive and negative correlations, respectively.Colored sections show the percentage of significant correlations at p < 0.05.

Figure 9 .
Figure 9. Percentages of the partial correlation coefficients between EOS and preseason temperature, precipitation, and SOS of each ecogeographical region in the Greater Khingan Mountains from 1982 to 2015.Bars above 0 and below 0 represent the percentage of positive and negative correlations, respectively.Colored sections show the percentage of significant correlations at p < 0.05.

Table 1 .
Natural [27]itions and major vegetation composition of each ecogeographical region in the Greater Khingan Mountains.Stable snow indicates the seasonal snow that has a relatively continuous spatial distribution and snow cover time (more than 60 days); unstable snow indicates seasonal snow that has discontinuous spatial distribution (mostly speckled distribution), with intermittent and relative short snow cover time (10-60 days)[27].

Table 1 .
Natural [27]itions and major vegetation composition of each ecogeographical region in the Greater Khingan Mountains.Stable snow indicates the seasonal snow that has a relatively continuous spatial distribution and snow cover time (more than 60 days); unstable snow indicates seasonal snow that has discontinuous spatial distribution (mostly speckled distribution), with intermittent and relative short snow cover time (10-60 days)[27].

Table 2 .
Root mean square error (RMSE), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) based assessment of the DL and AG fitting methods (best-fitting method shown in bold).

Table 2 .
Root mean square error (RMSE), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC) based assessment of the DL and AG fitting methods (best-fitting method shown in bold).