Forest Phenology Shifts in Response to Climate Change over China–Mongolia–Russia International Economic Corridor

: Vegetation phenology is a sensitive indicator of climate change. With the intensiﬁcation of global warming, the changes in growing seasons of various vegetation types have been widely documented across the world. However, as one of the most vulnerable regions in response to the global climate change, the phenological responses and associated mechanisms in mid–high latitude forests are still not fully understood. In this study, long-term changes in forest phenology and the associated relationship with the temperature and snow water equivalent in the China–Mongolia–Russia International Economic Corridor were examined by analyzing the satellite-measured normalized di ﬀ erence vegetation index and the meteorological observation data during 1982 to 2015. The average start date of the growing season (SOS) of the forest ecosystem in our study area advanced at a rate of 2.5 days / decade, while the end date of the growing season (EOS) was delayed at a rate of 2.3 days / decade, contributing to a growing season that was approximately 15 days longer in the 2010s compared to that in 1980s. A higher April temperature is beneﬁcial to the advance in the SOS, and a higher summer temperature has the potential to extend the EOS in the forest ecosystem. However, our results also suggest that a single temperature cannot fully explain the advance of the SOS, as well as the delay in the EOS. The preseason Snow Water Equivalent (SWE) is also an essential factor in inﬂuencing the growing season. A higher SWE in February and March and lower SWE in April tend to advance the SOS, while higher SWE in pre-year December and lower SWE in current year October are beneﬁcial to the extension of the EOS.


Introduction
Vegetation phenology is a sensitive indicator of climate change [1][2][3][4][5]. The advance or delay of the growing season reflects the long-term biological impacts of climate on the terrestrial ecosystem and, as a result, influences the productivity of the ecosystem [6]. Temperature records, as well as model simulations, show that global temperatures have increased obviously and significantly, especially in mid and high latitudes since the beginning of the 20th century [7][8][9]. Over the past several decades, both field observations and satellite records showed an extension of the growing season at regional and global scales [10][11][12].
Remote sensing measurements are capable of providing continuous spatiotemporal information and have become the most commonly used tool to monitor the phenology variations at large scales in recent decades [13]. Previous studies have investigated phenology variations in the Northern

Study Area
The China-Mongolia-Russia International Economic Corridor is located in North and Northeast Asia and ranges from 56 • E to 141.5 • E in longitude and from 37.3 • N to 62 • N in latitude, as shown in Figure 1. The China-Mongolia-Russia International Economic Corridor mainly covers three countries, namely Northeast and Inner Mongolia of China, Mongolia and the southern part of Siberia in Russia. The elevation of this study area ranges from −272 m to 4390 m. The climate of this region was characterized by a long period of snow-covered winter and a relatively short, mild summer. The climate records show that the annual average temperature ranges from −11 • C to 13 • C, and the yearly total precipitation ranges from 35 mm to 990 mm in this region. Forests are the dominant land cover types in the area. Forests in this region include four different vegetation types, specifically broadleaved deciduous, needleleaved evergreen, needleleaved deciduous and mixed forest, which are classified as European Space Agency (ESA) Climate Change Initiative (CCI) land cover products used in this study [27].

Study Area
The China-Mongolia-Russia International Economic Corridor is located in North and Northeast Asia and ranges from 56° E to 141.5° E in longitude and from 37.3° N to 62° N in latitude, as shown in Figure 1. The China-Mongolia-Russia International Economic Corridor mainly covers three countries, namely Northeast and Inner Mongolia of China, Mongolia and the southern part of Siberia in Russia. The elevation of this study area ranges from −272 m to 4390 m. The climate of this region was characterized by a long period of snow-covered winter and a relatively short, mild summer. The climate records show that the annual average temperature ranges from −11 °C to 13 °C, and the yearly total precipitation ranges from 35 mm to 990 mm in this region. Forests are the dominant land cover types in the area. Forests in this region include four different vegetation types, specifically broadleaved deciduous, needleleaved evergreen, needleleaved deciduous and mixed forest, which are classified as European Space Agency (ESA) Climate Change Initiative (CCI) land cover products used in this study [27]. However, the mid-and high-latitude regions of the Northern Hemisphere have experienced extraordinary warming since the beginning of the 20th century-far higher than the global average level [9]. The warming here directly influences the snow cover, as well as the phenology stages, including the start or the end date of the growing season. Understanding the forest phenology changes, as well as the contributors in this area, is helpful in tracking the natural climate change and forecasting future forest shifts in the northern hemisphere.

Data
Global Inventory Modeling and Mapping Studies (GIMMS) normalized difference vegetation index (NDVI) time-series data, derived from the advanced very high-resolution radiometer (AVHRR) from 1982 to 2015, were used to detect the phenology changes in our study area [28]. Previous studies have reported that a 6-16 day temporal resolution of the satellite data can generate the optimal estimate of the phenological stages [29]. The temporal resolution of the GIMMS NDVI3g data is 15 days, and spatial resolution is 8 km, which is suitable for large-scale ground phenological monitoring. The dataset ranges from July 1981 to December 2015. To eliminate the uncertainty caused by the data quality, we used quality flags equal to 0 (best quality) and 1 (good quality) to select the good values for the phenology model. The fitting weight for best qualities was set to 1, and for the good quality pixels it was set to 0.5.
The 300 m resolution global land cover data series from 1992 to 2015 produced by the European Space Agency (ESA) Climate Change Initiative (CCI) were used to extract the forest distribution in our study area [27]. The CCI land cover data created for climate studies have been widely used to However, the mid-and high-latitude regions of the Northern Hemisphere have experienced extraordinary warming since the beginning of the 20th century-far higher than the global average level [9]. The warming here directly influences the snow cover, as well as the phenology stages, including the start or the end date of the growing season. Understanding the forest phenology changes, as well as the contributors in this area, is helpful in tracking the natural climate change and forecasting future forest shifts in the northern hemisphere.

Data
Global Inventory Modeling and Mapping Studies (GIMMS) normalized difference vegetation index (NDVI) time-series data, derived from the advanced very high-resolution radiometer (AVHRR) from 1982 to 2015, were used to detect the phenology changes in our study area [28]. Previous studies have reported that a 6-16 day temporal resolution of the satellite data can generate the optimal estimate of the phenological stages [29]. The temporal resolution of the GIMMS NDVI3g data is 15 days, and spatial resolution is 8 km, which is suitable for large-scale ground phenological monitoring. The dataset ranges from July 1981 to December 2015. To eliminate the uncertainty caused by the data quality, we used quality flags equal to 0 (best quality) and 1 (good quality) to select the good values for the phenology model. The fitting weight for best qualities was set to 1, and for the good quality pixels it was set to 0.5.
The 300 m resolution global land cover data series from 1992 to 2015 produced by the European Space Agency (ESA) Climate Change Initiative (CCI) were used to extract the forest distribution in our study area [27]. The CCI land cover data created for climate studies have been widely used to describe the vegetation patterns in previous studies [30][31][32]. To avoid the impact of land cover changes on the phenology variation, only the unchanged forest during 1992 to 2015 was analyzed in this paper. Due to fewer human activities, the land use and land cover changes were not significant from 1975 to 2010 in Siberia [33]. Although the land use changes were significant in Northeast China from the 1980s to 2000, the forests changes mainly concentrated on deforestation [34]. Based on the previous studies and the data availability, forest increase or conversion between different forest types from 1982 to 1991 were not taken into consideration in this study. To make sure the CCI land cover data could capture the spatial distribution for different forest types, we used the 1:4,000,000 vegetation map in Russia and the 1:1,000,000 vegetation map in China to qualitatively evaluate the CCI's accuracy. The two datasets showed good consistency in the forest and open woodlands, especially for dark and light coniferous forests and broad-leaved forests. It should be noted that the CCI data described the spatial distribution of different forest types with a higher resolution.
Mean air temperature, maximum air temperature, minimum air temperature and precipitation from the year 1982 to the year 2015 were obtained from Climate Research Unit (CRU) TS4.3 monthly data with a spatial resolution of 0.5 • [35]. The monthly SWE data in the same time period were obtained from GlobSnow SWE v2.0 data [36], which combine ground-based weather station data with satellite-based passive microwave measurements. We rescaled both the climate data and land cover data to 8 km using the bilinear spatial interpolation method to match the spatial resolution of the NDVI products.

Methods and Analyses
Logistic models have been commonly used to detect phenological events [15,37]. Logistic models have good conceptual coherency with biophysical inference in terms of being a "growth curve" and are the primary technique used to extract phenological information from satellite observation data [38]. The double logistic functions have proved to be an effective way to accurately monitor the ground phenological transition dates by optimally removing outliers and reducing noise in remote sensing time series data [39,40]. In this study, we used a double logistic function to fit the time series NDVI data for each year: where x 1 and x 3 determine the position of the left and right inflection point, respectively, while x 2 and x 4 represent the rate of change at the corresponding inflection points. Based on our previous studies [41], we selected the dynamic thresholds of 0.2 and 0.3 to extract the start of the season (SOS) and the end of the season (EOS), respectively. Then, the length of the season can be obtained by subtracting the SOS from the EOS. As the SOS and EOS extraction from the first (1982) and the last (2015) year was not robust, the results between 1983 and 2014 were analyzed. We used least square regression to detect the changing trend for the SOS and EOS from 1982 to 2015 at the pixel scale. Our previous study used this method to identify the change slope for the growing season in Northeast China, and detailed information about this method can be found in Yu et al. (2017) [41]. We evaluated the response of spring and autumn phenology to monthly temperatures and SWE based on the partial least-squares (PLS) regression (implemented in JMP 13, SAS Institute). PLS regression is demonstrated to be suitable for highly autocorrelated variables [18,20,42], such as monthly temperature or SWE, as used here. This analysis can provide both influence direction (model coefficients, MCs) and the importance of the independent variables (variable importance plot, VIP) in explaining the variation in the response variable [18,20,42]. We spatially averaged data for all forest pixels of the temperature and SWE, as well as phenology grids for this analysis. In this study, the MC and VIP values can reflect the direction and importance of monthly temperature and SWE in influencing the SOS and EOS.

Characterization of Forest Phenology
The start date (SOS) and the end date (EOS)of the growing season across the forest in the China-Mongolia-Russia International Economic Corridor was detected by remote sensing during our studying period, as shown in Figure 2. On average, the SOS and EOS occurred on 120.8 ± 7.8 and 290.9 ± 9.5 Julian days across the forest landscape in the study area. The start and the end date of the growing season varied significantly across our study area, with some areas in the northwest beginning to green up before the 110th Julian day (20 April), while some areas north of Baikal and the mountain areas started to green up after the 135th Julian day (15 May). The end date of the growing season also showed remarkable spatial heterogeneity. The EOS of the west part of the economic corridor occurred before the 275th Julian day (5 October), while the EOS in the southernmost of the study area was close to late October. It should be noted that the SOS and EOS varied across different forests. The detailed SOS and EOS days are shown in Table 1.

Characterization of Forest Phenology
The start date (SOS) and the end date (EOS)of the growing season across the forest in the China-Mongolia-Russia International Economic Corridor was detected by remote sensing during our studying period, as shown in Figure 2. On average, the SOS and EOS occurred on 120.8 ± 7.8 and 290.9 ± 9.5 Julian days across the forest landscape in the study area. The start and the end date of the growing season varied significantly across our study area, with some areas in the northwest beginning to green up before the 110th Julian day (20 April), while some areas north of Baikal and the mountain areas started to green up after the 135th Julian day (15 May). The end date of the growing season also showed remarkable spatial heterogeneity. The EOS of the west part of the economic corridor occurred before the 275th Julian day (5 October), while the EOS in the southernmost of the study area was close to late October. It should be noted that the SOS and EOS varied across different forests. The detailed SOS and EOS days are shown in Table 1.

Trends in Growing Season of Forest Phenology
To investigate the variations in forest phenology stages, we calculated the trend of the SOS and EOS from 1983 to 2014 at the pixel scale in our study area, as shown in Figure 3. The trend of the SOS, determined by the linear regression, was negative in approximately 95.7% of the forest, especially in

Trends in Growing Season of Forest Phenology
To investigate the variations in forest phenology stages, we calculated the trend of the SOS and EOS from 1983 to 2014 at the pixel scale in our study area, as shown in Figure 3. The trend of the SOS, determined by the linear regression, was negative in approximately 95.7% of the forest, especially in Forests 2020, 11, 757 6 of 13 the broadleaved deciduous forests and evergreen forests. The SOS trend was significantly negative (p < 0.05, t-test) in about 40.2% of the pixels, indicating significant advances in the SOS. The SOS advance exceeds 0.2 days/year, which is 6.4 days from 1983 to 2014, in nearly 41.7% of the area. Though the majority of pixels exhibited SOS advance, positive trends were also observed, and only 4.3% of the pixels showed an SOS delay during our study period. Significant spatial heterogeneity was also detected. Among the four forest types, the SOS advance in the broadleaved deciduous forest was the slowest, with 0.14 days/year, while the mixed forest showed the fastest SOS advances, which were 0.30 days/year on average. the broadleaved deciduous forests and evergreen forests. The SOS trend was significantly negative (p < 0.05, t-test) in about 40.2% of the pixels, indicating significant advances in the SOS. The SOS advance exceeds 0.2 days/year, which is 6.4 days from 1983 to 2014, in nearly 41.7% of the area. Though the majority of pixels exhibited SOS advance, positive trends were also observed, and only 4.3% of the pixels showed an SOS delay during our study period. Significant spatial heterogeneity was also detected. Among the four forest types, the SOS advance in the broadleaved deciduous forest was the slowest, with 0.14 days/year, while the mixed forest showed the fastest SOS advances, which were 0.30 days/year on average.
In contrast to the SOS, the trend for the EOS was positive in 64.9% of the forests, indicating a delayed end of growing season, and 75.6% of which was significantly positive (p < 0.05, t-test). The EOS delay exceeded 0.2 days/year in nearly 55.8% of the forest area. Similar to the SOS, the changing trend in the EOS showed significant spatial distribution differences. Areas, such as the Greater Khingan Mountains and central Siberia, showed a delayed trend higher than 0.6, indicating that the end of the growing season was extended by more than 19 days. A negative EOS trend was mainly found in the north part of our study area, and the pixels with an advance in EOS account for 4.8% of the forest area. To analyze the temporal trends in the SOS and EOS, we averaged the change slope over all pixels covered by forest, as shown in Figure 4c,d. The regional mean SOS over the forest in the China-Mongolia-Russia International Economic Corridor advanced over the past 32 years (R 2 = 0.264, p = 0.0026), with an annual advance of 0.25 days per year, as shown in Figure 4c. The regional mean of the EOS was delayed by 7.3 days during our study period (R 2 = 0.356, p = 0.0003), with a delay rate of 0.23 days per year, as shown in Figure 4d. To analyze the climate background change trend in this area, we calculated the annual average temperature and SWE for the forest grids. The regional mean temperature in this area increased at a rate of 0.019 °C/year from 1983 to 2014 (R 2 = 0.071, p = 0.1407). The mean SWE increased 2.74 mm (R 2 = 0.129, p = 0.0437), with an increased rate of 0.086 mm/year. In contrast to the SOS, the trend for the EOS was positive in 64.9% of the forests, indicating a delayed end of growing season, and 75.6% of which was significantly positive (p < 0.05, t-test). The EOS delay exceeded 0.2 days/year in nearly 55.8% of the forest area. Similar to the SOS, the changing trend in the EOS showed significant spatial distribution differences. Areas, such as the Greater Khingan Mountains and central Siberia, showed a delayed trend higher than 0.6, indicating that the end of the growing season was extended by more than 19 days. A negative EOS trend was mainly found in the north part of our study area, and the pixels with an advance in EOS account for 4.8% of the forest area.
To analyze the temporal trends in the SOS and EOS, we averaged the change slope over all pixels covered by forest, as shown in Figure 4c,d. The regional mean SOS over the forest in the China-Mongolia-Russia International Economic Corridor advanced over the past 32 years (R 2 = 0.264, p = 0.0026), with an annual advance of 0.25 days per year, as shown in Figure 4c. The regional mean of the EOS was delayed by 7.3 days during our study period (R 2 = 0.356, p = 0.0003), with a delay rate of 0.23 days per year, as shown in Figure 4d. To analyze the climate background change trend in this area, we calculated the annual average temperature and SWE for the forest grids. The regional mean temperature in this area increased at a rate of 0.019 • C/year from 1983 to 2014 (R 2 = 0.071, p = 0.1407). The mean SWE increased 2.74 mm (R 2 = 0.129, p = 0.0437), with an increased rate of 0.086 mm/year.  The regression analysis between the SOS and the climatic factors, including the annual average temperature and SWE, indicated negative sensitivity of SOS to temperature (R 2 = 0.18, p = 0.0165). An increase in annual average temperature would advance the SOS by 2.81 days across the forest in our study area, indicating that the temperature plays an essential role in regulating the SOS. The observed EOS changes are not linear responses to temperature, meaning that the changes in EOS do not result from a monotonic vegetation response to temperature. To further investigate the climatic driving factors for the phenology changes, we used PLS regression to analyze the dominant variables in determining the SOS and EOS.  The regression analysis between the SOS and the climatic factors, including the annual average temperature and SWE, indicated negative sensitivity of SOS to temperature (R 2 = 0.18, p = 0.0165). An increase in annual average temperature would advance the SOS by 2.81 days across the forest in our study area, indicating that the temperature plays an essential role in regulating the SOS. The observed EOS changes are not linear responses to temperature, meaning that the changes in EOS do not result from a monotonic vegetation response to temperature. To further investigate the climatic driving factors for the phenology changes, we used PLS regression to analyze the dominant variables in determining the SOS and EOS.  For the SOS, the variable importance plots (VIPs), as shown in Figure 5a,b indicate that the temperature in pre-year July, March, April and May, and the SWE in the winter half-year from November to the following April, play an essential role in determining the SOS in the forest ecosystem in our study area (VIP > 0.8). Temperatures in spring (March, April and May) have negative model coefficients, indicating that the warm temperature during these months advances the SOS. Temperatures in pre-year July have a positive model coefficient, indicating that the higher temperature may contribute to the delay of the SOS. It is clear that the temperature in April plays the dominant role in determining the SOS (MC = −0.86, VIP value = 2.45). We also calculated the relationship between SOS and Tmax and Tmin, which showed similar results, and no major difference was obtained. SWE in February and March has negative coefficients, suggesting that higher SWE in these months contributed to the advance of SOS, while the SWE in months such as November, December, January and April have positive model coefficients, indicating that the higher SWE in these months tends to delay the SOS. Figure 5c,d illustrates the EOS responses to monthly temperature and SWE. Combining the model coefficients and VIP values, the temperature in months such as June, July, October and the pre-year December and the SWE in pre-year December and October determined the EOS. The high positive model coefficients in June, July and October indicate that the warmer temperature tends to delay the EOS, while the high negative model coefficient in pre-year December indicates that the For the SOS, the variable importance plots (VIPs), as shown in Figure 5a,b indicate that the temperature in pre-year July, March, April and May, and the SWE in the winter half-year from November to the following April, play an essential role in determining the SOS in the forest ecosystem in our study area (VIP > 0.8). Temperatures in spring (March, April and May) have negative model coefficients, indicating that the warm temperature during these months advances the SOS. Temperatures in pre-year July have a positive model coefficient, indicating that the higher temperature may contribute to the delay of the SOS. It is clear that the temperature in April plays the dominant role in determining the SOS (MC = −0.86, VIP value = 2.45). We also calculated the relationship between SOS and Tmax and Tmin, which showed similar results, and no major difference was obtained. SWE in February and March has negative coefficients, suggesting that higher SWE in these months contributed to the advance of SOS, while the SWE in months such as November, December, January and April have positive model coefficients, indicating that the higher SWE in these months tends to delay the SOS. Figure 5c,d illustrates the EOS responses to monthly temperature and SWE. Combining the model coefficients and VIP values, the temperature in months such as June, July, October and the pre-year December and the SWE in pre-year December and October determined the EOS. The high positive model coefficients in June, July and October indicate that the warmer temperature tends to delay the EOS, while the high negative model coefficient in pre-year December indicates that the temperature chilling is beneficial to the delaying of the EOS. The positive coefficient in pre-year December, March and April suggests that the higher SWE in these months tends to delay the EOS, while the negative coefficients in pre-year November, January and October indicate that the higher SWE might contribute to advances in EOS. PLS regression for different forest types was also conducted, but the results did not show significant differences.

Trends and Spatiotemporal Variations in Forest Phenology
In this study, we investigated the start (SOS) and the end (EOS) date of the growing season in the forest ecosystem over the China-Mongolia-Russia International Economic Corridor by using long-term GIMMS NDVI records and the double logistic model from 1982 to 2015. The SOS and EOS of forests in our study area are 120.8 ± 7.8 and 290.9 ± 9.5 Julian days, respectively, and the spatial distribution pattern of these two phenological parameters is roughly consistent with two previous studies in midand high-latitudes in the Northern hemisphere [14,15].
Similar to numerous previous studies, our results showed a significant advance in the SOS with the background of climate warming [13,18]. We found that the average SOS of the forest ecosystem advanced at a rate of 2.5 days/decade in our study area. Myneni et al. (1997) reported that the SOS advanced by 8 ± 3 days from 1981 to 1991 in northern latitudes  • N) [11]. Deng et al. (2019) found that the SOS in boreal forests varied with an average advancement rate of 3.38 days/decade [17]. The difference in trends is mainly attributed to the mismatch of the spatial extension and temporal duration. Even utilizing the same time series data, Jeong et al. (2011) found that the SOS advanced by 5.2 days from 1982-1999 but advanced by only 0.2 days from 2000-2008 in the Northern Hemisphere's temperate vegetation [14]. It should also be noted that our study also detected a significant delay in the EOS, at a rate of 2.3 days/decade. Usually, the SOS is considered to be the most sensitive phenological parameter in response to global warming, while the EOS changes show more uncertainty [23]. Our results indicate that the EOS is as significant as the SOS in response to climate change and in determining the length of the growing season. Undeniably, however, the trends in the EOS appeared to be more spatially variable compared with those of the SOS, suggesting a more complicated feedback mechanism of the EOS to climate variation.

Relationship between SOS/EOS and Climatic Factors
The warming temperature triggering the advance of the SOS for diverse vegetation types has been widely reported [16,43]. Our results found that the SOS of the forest ecosystem in our study area was highly dependent on the temperature in April (MC = −0.86, VIP = 2.45). Though the PLS regression showed that the temperatures in other spring months, such as March and May, are also important variables in determining the advance of the SOS, their model coefficients are relatively small, indicating a relatively slight impact. Different to other studies, we chose the SWE instead of precipitation as a latent factor influencing forest phenology variations. Our results found that the winter SWE from pre-year November to current year April are essential factors in determining the SOS. The SOS is negative to SWE in February and March, while positive in other winter months. A greater snow depth can increase the amount of available water used for vegetation growth as well as maintain the soil temperature and nutrients [22]. Furthermore, the phenological changes in the boreal regions may be partially related to snowmelt timing variations [44,45]. Our results suggest that the higher SWE in February and March created adequate water, temperature and nutrient conditions for the forest green-up and the rapid increase in temperature accompanied by the decrease in snow water equivalent in April, which accelerated the advance of the SOS.
The autumn dormancy and the associated mechanisms in the forest are not well documented as compared with those of the SOS [23]. Our results found that the warmer summer is conducive to the extension of the growing season, which was consistent with most published studies [46,47]. Though the SOS showed different responses to monthly SWE, our results showed that the pre-year December and current year October SWE play a dominant role in determining the EOS. Previous studies showed that winter chilling has the potential to advance the SOS [48,49]. Our results suggest that the temperature chilling, as well as higher SWE in December, tend to have more potential to extend the EOS instead of the SOS. As the winter chilling requirements can always be sufficient in this mid-high region, the chilling environment and the snow cover protection may trigger more robust photosynthesis and a more extended photoperiod in the next growing season, as documented by Grippa et al. (2005) [24]. On the other hand, our study also suggested that using the annual values or preseason values to reveal the relationship between the SOS/EOS and its latent driving factors may be defective unless the response trends of these months are consistent between each other.

Uncertainties and Future Works
In this study, we estimated the forest phenology in the China-Mongolia-Russia International Economic Corridor, based on remote sensing observation records. By using the available in situ phenology observation data, we verified that the start and end dates of the growing season in Northeast China and the root mean square error of the simulated SOS and EOS are 3.14 and 3.90 days [41], indicating that the double logistic model and the threshold can detect the phenology effectively. Though our results showed a similar SOS and EOS pattern compared with previous studies [14,15], the in situ phenology observations in Mongolia and Russia are still further required to verify our results presented here.
It should be noted that, though the needleleaved evergreen forest shared similar average SOS and EOS values with other forest types, the standard deviations are much higher, suggesting that significant phenology uncertainties exist in the needleleaved evergreen forest. With evergreen leaves, the changes in the photosynthetic capacity and dormancy are hard to detect [50]. Wu et al. (2014) documented that widely used NDVI and EVI may exhibit limited potential in tracking growing season phenology, especially for the evergreen needleleaf forest ecosystems [39], which can explain why some needleleaved evergreen forest showed smaller SOS values (<110 days of Julian) in spring and larger values in autumn (>300 days in the autumn). A field-based survey, as well as a study of the temperature thresholds, should be further conducted to verify the accuracy of our results.
Furthermore, the boreal forests are vulnerable to global warming, which may contribute to the forest shift, and the increased fire or insect disturbances may change the forest pattern on the local scale [51]. Though we already eliminated the uncertainty caused by the land cover changes through using unchanged forest cover pixels based on the long-term land cover series , the lack of land use and land cover data before 1992 may result in some forest conversions due to human activities being ignored and, as a result, add to the uncertainty of our results.

Conclusions
In this study, the long-term changes in the forest phenology and the associated relationship with the temperatures and snow water equivalent over the China-Mongolia-Russia International Economic Corridor were examined by analyzing satellite-measured normalized difference vegetation index and meteorological observation data during 1982 to 2015. The average SOS of the forest ecosystem in our study area advanced at a rate of 2.5 days/decade, while the EOS was delayed at a rate of 2.3 days/decade, contributing to a growing season approximately 15 days longer in the 2010s compared to that in 1980s. A higher April temperature is beneficial to the advance in the SOS, and the higher summer temperature has the potential to extend the EOS in the forest ecosystem. However, our results also suggested that a single temperature cannot fully explain the advance of the SOS, as well as the delay in the EOS. The preseason SWE is also an essential factor in influencing the growing season. A higher SWE in February and March, and lower SWE in April, tend to advance the SOS, while higher SWE in pre-year December and lower SWE in current year October are beneficial to the extension of the EOS. Our results suggest both the temperature and snow water equivalent should be combined in predicting the forest phenology variations over China-Mongolia-Russia in the future. Furthermore, our findings would have further implications for reasonably estimating the productivity of the forest ecosystem and evaluating the regional ecosystem service value under the changing climate and economic development.
Author Contributions: Conceptualization, methodology and writing, L.Y.; climate data processing, Z.Y.; project administration and funding acquisition, S.Z. All authors have read and agreed to the published version of the manuscript.