Temporal Trends and Spatial Variability of Vegetation Phenology over the Northern Hemisphere during 1982-2012

Satellite-derived vegetation phenology has been recognized as a key indicator for detecting changes in the terrestrial biosphere in response to global climate change. However, multi-decadal changes and spatial variation of vegetation phenology over the Northern Hemisphere and their relationship to climate change have not yet been fully investigated. In this article, we investigated the spatial variability and temporal trends of vegetation phenology over the Northern Hemisphere by calibrating and analyzing time series of the satellite-derived normalized difference vegetation index (NDVI) during 1982–2012, and then further examine how vegetation phenology responds to climate change within different ecological zones. We found that during the period from 1982 to 2012 most of the high latitude areas experienced an increase in growing period largely due to an earlier beginning of vegetation growing season (BGS), but there was no significant trend in the vegetation growing peaks. The spatial pattern of phenology within different eco-zones also experienced a large variation over the past three decades. Comparing the periods of 1982–1992, 1992–2002 with 2002–2012, the spatial pattern of change rate of phenology shift (RPS) shows a more significant trend in advancing of BGS, delaying of EGS (end of growing season) and prolonging of LGS (length of growing season) during 2002–2012, overall shows a trend of accelerating change. Temperature is a major determinant of phenological shifts, and the response of vegetation phenology to temperature varied across different eco-zones.


Introduction
Vegetation phenology is highly sensitive to climate change, so monitoring vegetation phenology is becoming an increasing important way to understand and quantify global environmental changes [1,2]. Since the 1950s, the world's mean surface temperature has increased by 0.6°C, and the warming trend has been faster in the high latitude areas, such as Eurasian continent [3]. As a result of the recent increase in surface temperature, earlier spring phenological events have been investigated in the northern hemisphere [1], [2], [4]. Shifts in vegetation phenology will affect ecosystem functions since they affect hydrological cycle, surface energy balance, and terrestrial carbon cycle [5][6]. Therefore, it is important to improve our understanding to quantify accurately the timing of phenological shifts for improving our understanding of terrestrial ecosystems response to global change.
In the past few decades, satellite remote sensing has been recognized as a valuable tool for acquiring and analyzing information about land surface development. In particular, time series of remote sensing data are an important data source for investigating vegetation phenological changes. The seasonal dynamics in a long-term vegetation index (VI) time series such as those from the AVHRR, SPOT, or MODIS sensors can be used for plant phenological detection. For example, time series of NDVI data derived from AVHRR or MODIS data have been used to estimate the inter-annual variations in vegetation phenology for the past decades [7][8]. Annual time series of Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) data have been used to identify key phenological events [9]. Recent, research efforts have been focused towards improving algorithms for calculating vegetation phenological shifts using other satellite-derived data (e.g. MODIS-enhanced vegetation index (EVI), leaf area index (LAI), FAPAR and Albedo) and modeling techniques.
Many previous studies investigated the timing of phenological events and their connection to climate change by using time-series of remote sensing data in the northern hemisphere [2], [10][11][12][13][14][15]. Due to the increase in land surface temperature, strong vegetation phenology shifts have occurred in the northern hemisphere. For example, in North America, a delayed dormancy onset date (0.551 days yr -1 ) and an extended growing season length (0.683 days yr -1 ) were observed over the mid and high latitudes during 1982-2006 [16]. In China, the beginning of the growing season (BGS) has advanced in spring by 0.79 days yr -1 and the end of the growing season (EGS) delayed in autumn by 0.37 days yr -1 in temperate regions during 1982-1999 [17]. In Europe, an earlier BGS (0.54 days yr -1 ) and a prolonged length of vegetation growing season (LGS) (0.96 days yr -1 ) were found during 1982-2001 [18]. For the Northern Hemisphere as a whole, BGS advanced by 5.2 days during 1982-1999 but advance by 0.2 days during 2000-2008, and EGS was delayed by 4.3 days in the period of 1982-1999 and it was further delayed by another 2.3 days in the period of 2000-2008 [2]. However, although it took a great deal of effort to examine phenological shifts over the northern hemisphere, great uncertainties still exist. Just as our literature review revealed that many studies present different magnitudes of phenological shifts, including earlier or later phenology events and longer or shorter vegetation growing seasons [13], [19][20][21]. Furthermore, some other studies only concentrated on earlier BGS or assessed later shifts of dormancy events [9], [22][23][24].
In this paper, we intend to investigate the spatial variability and temporal trends of vegetation phenology shifts in the northern hemisphere during 1982-2012, and further examine the mechanisms of phenology response to climate change within different ecological zones. Our objectives in this study are to better understand the following key scientific questions: (1) What is the vegetation phenology pattern on the pixel level in the northern hemisphere? (2) How did the patterns of phenology shift within different ecological zones over the past three decades? And (3) what are the key factors affecting vegetation growing season change in recent years? Answers to these questions are essential to improve our understanding of phenological shifts caused by climate warming and to formulate measures for climate warming and environmental problems at the global level. still some cloud contamination existing in some cases. To eliminate errors caused by clouds, snow and ice contamination, we applied a Savitzky-Golay filtering procedure to each annual NDVI cycle for smoothing and reconstrucing the NDVI time series, as described by Chen et al. [34] and Broich et al. [35]. For MODIS MOD13C1 product, a summary pixel reliability layer has been included in the MOD13C1, we first use the reliability layer information to locate the cloud or ice/snow pixels (where the pixel value equals 2 or 3 in the reliability layer) in the NDVI time-series. Then these cloud or ice/snow pixels can be replaced by a linear interpolation method using the adjacent NDVI data. The new NDVI values can be calculated as follows: ( Where PR (Pixel Reliability) is the value of pixel reliability and NDVI(i, t) is the value of NDVI of the ith pixel in the tth time.  To further minimize the influence of soil and non-vegetated signals in NDVI time series, we chose targeted study regions with multiyear average NDVI greater than 0.1 in the Northern Hemisphere based on GIMMS (1982GIMMS ( -2006 and MODIS (2002-2012) NDVI data.

Determination of Phenological Parameters
Several methods have been developed to detect the phenophases (greenup, maturity, senescence, and dormancy) of vegetation phenology by using NDVI data, such as NDVI thresholds [7,36], largest NDVI increase [37], backward-looking moving averages [38], and fitting logistic functions [38][39]. Comparative studies have shown that the Midpoint pixel threshold method based on variations of an NDVI ratio is the one of the most consistent methods for the estimation of the BGS, EGS, LGS and peak of growing season (PGS) over a variety of ecosystems and corroborates with ground-based phenology data [7]. Thus, for the determination of BGS, EGS based on NDVI, we used the Midpoint pixel threshold method, where the value midpoint was set at fifty percent of the maximum NDVI [7]. Firstly we calculated the annual NDVI ratio time series curve during 1982-2012. Subsequently, the phenology date was determined for each pixel by using the midpoint pixel method, and BGS or EGS was determined according to the slope of the fitted curve. Fig 2 shows the general scheme of the method for determining the phenological parameters. The NDVI ratio was calculated by using the following formula: Where NDVI ratio is the output ratio, ranging from 0-1, NDVI is the daily NDVI, NDVI max is the annual maximum NDVI, and NDVI min is the annual minimum NDVI. Thus the phenological parameters can be extracted with a Midpoint pixel thresholds method for each pixel.

Phenology Trend Analysis
The spatial distribution of the phenology shift trends in the Northern Hemisphere from 1982-2012 based on the GIMMS and MODIS NDVI datasets was characterized by the change rate of phenology shift (RPS), which was calculated by the following relationship: Where i is the year for 1 to n; n is the total number of years; and Day i is the seasonal phenology (BGS, EGS, LGS and PGS) of year i. A positive value of RPS indicates an increase of phenology in a time series. The timings of these parameters were calculated to the day because we took the day of the satellite observation into account. Correlation coefficients and significance tests were used to examine the ralationship betweeen climate and phenology shift over the period of analysis, 1982-2012. We analysed the effects of seasonal temperature and precipitation on the trends of phenology shift in different Holdridge life ecozones by Pearson correlation analysis. In addition, a t-test determined the significance of the correlation, with p values <0.1, 0.05 and 0.01 are considered to be low, middle and high significance levels, respectively. LGS and PGS over the northern hemisphere. The BGS dates ranged widely from approximately 15 days (15 January) at low latitudes to approximately 225 days (15 July) at high latitudes. The mean BGS date was progressively delayed with increasing latitude. Concurrently, the progressive patterns that vegetation turn green gradually along latitudinal gradient are mainly attributed to the climatological temperature [2], [40]. The BGS in Eurasia was much earlier than in North America at the same latitude, and the latest dates of BGS occured in northern Canada and northern Siberia owing to low temperatures. The EGS dates ranged from 15 May to approximately 30 November, the LGS value ranged from 64 days to 240 days, and the PGS value ranged from 64 days to approximately 240 days. There were opposite spatial patterns between BGS and EGS over the northern hemipshere. The dates of vegetation senescence occured in reverse order of the green-up onset. Compared with other regions in the hemisphere, east Asia, Europe and southern North America had earlier BGS and later EGS, and thus had longer growing seasons, with some regions lasting up to 10 months. In contrast, in northern Canada and northern Siberia the vegetation growing seasons lasted only 1 month.  America and Eurasia, with exceptions being parts of Russia, which showed a delayed BGS (Fig  4a). A slight pattern change occurred in EGS, with more than 65% of the study region experiencing an advancing trend in EGS (Fig 4b). A notable change occurred in LGS, which became about 80% longer (over 0.5 days yr -1 ) in northern America and Eurasia (Fig 4c), which was similar to BGS trend. Results indicated that a longer LGS trend is mainly caused by a significant advance in spring phenology during 1982-1992, which is consistend with previous studies [17], [41]. During 1992-2002, BGS advanced throughout most of the hemisphere, with exceptions being parts of central northern America and southern Asia, while EGS was delayed throughout almost the entire hemisphere, resulting in a longer LGS (Fig 4e). During 2002-2012, a marked advancing trend (over 0.8 days yr -1 ) in BGS occurred form 55°N to 75°N (Fig  4i), correspondingly, the LGS can also be found to be siginificantly prolonged in this region (Fig 4m). From 1982 to 2012, unlike BGS, EGS and LGS, PGS showed a more fragmented pattern, but overall there was not a significant trend in PGS. From the point of spatial patterns, PGS was more relvant with BGS than EGS or LGS (Fig 4d, 4h and 4m).  Table 1, changes in BGS, EGS, LGS and PGS showed high interannual variations. Long-term variations in BGS, EGS and LGS showed 'earlier green-up onset', 'delayed dormancy onset', and 'prolonged growing season' features, which are consistent with previous studies [2], [10],   Fig 5b), which reflects the vegetation phenology response to global warming for different latitude. Linear regression between phenological timing and latitude indicates that the rate of change is 1.24 days per degree of latitude (R 2 = 0.96, P < 0.001) for BGS, -0.51 days (R 2 = 0.57, P < 0.001) for EGS, -2.02 days (R 2 = 0.97, P < 0.001) for LGS and 0.54 days (R 2 = 0.76, P < 0.001) for PGS in Northern Hemisphere between 30°N and 75°N.

Spatial and Temporal Variations of Phenological Metrics over the Northern Hemisphere
The variation in mean BGS, EGS and PGS for specific vegetation classes were estimated by using AVHRR data and MODIS data.   BGS, EGS, LGS and PGS for specific vegetation classes during 1982-2012. It was found that deciduous broadleaf forest and woodland exhibited a higher rate of change for both BGS (−0.24 and −0.25 days yr −1 , respectively) and EGS (+0.48 and +0.33 days yr −1 , respectively) than the other classes. Shrubland showed the least variation in BGS (-0.08 days yr −1) and EGS (0.05 days yr −1 ) compared to the other classes. Interestingly, including deciduous broadleaf forest, woodland and deciduous needleleaf forest, the delay in EGS is consistently larger than the advancement in BGS,which leads to the inference that autumnal response of vegetation is greater and faster than that in spring, which is consistent with the findings of previous studies [42]. All these specific vegetation classes revealed a prolonging growth cycle with earlier BGS and delayed EGS.

Variation in Phenology within Different Holdridge Life Eco-Zones
In order to examine the phenological trends and their characteristic differences, the index RPS was calculated for the Holdridge life eco-zones in the study area. The index RPS was divided arbitrarily into four groups based on the value of the change rate of phenology shift. The first group represented the slight mode of phenological shift, in which the value of RPS was from 0.2 to 0.5 or from -0.5 to -0.2. The second group represented the moderate mode where the RPS was from 0.5 to 0.8 or from -0.8 to -0.5. The third group belonged to the fast mode with the RPS from 0.8 to 1.0 or from -1.0 to -0.8, and the group four was for the acute mode with the vulue of RPS above 1.0 or below -1.0. Fig 7 indicated that during 1982-1992, for BGS, most of the eco-zones, including boreal coniferous forest (Ba), boreal mountain system (BM), boreal tundra woodland (Bb), temperate continental forest (TeDc) and temperate steppe (TeBSK), were in the first RPS group (Fig 7a, a light blue area shows a slight advance trend in spring phenology). For EGS, no significant change occurred in many eco-zones (Fig 7b, yellow areas). Due to most of the eco-zones experienced an advancing trend in BGS, resulting in a longer growing season length. So for LGS, boreal coniferous forest belonged to the group one, boreal mountain system and temperate steppe belonged to the group two (Fig 7c, light red areas), all eco-zones experienced an extend trend in LGS. During 1992-2002, the same spatial pattern of RPS as the period of 1982-1992, but we should note that a remarkable extend occurred in the vegetation growing length, due to the 'earlier green-up onset', 'delayed dormancy onset' features in eco-zones. During 2002-2012, compared with the period of 1982-1992 and 1992-2002, the spatial pattern of RPS showed a more significant trend in advancing of BGS, delaying of EGS and prolonging of LGS during 2002-2012, overall shows a trend of accelerating change. Such as boreal coniferous forest ecozone, belonged to the group three, group two and group four for BGS, EGS and LGS, respectively (Fig 7g, 7h and 7i). It should be noted that in some eco-zones, such as polar, temperate desert and tropical desert, due to poor vegetation cover and snow (ice) and soil contamination, the spatial pattern of RPS shows a meaningless trend for the whole period of 1982-2012.

Analysis of the Relationship "Phenology-Climate" within Different Eco-Zones
In order to investigate the relationship between vegetation phenology and climate factors, correlation coefficients were calculated over Holdridge life eco-zones from 1982 to 2012 based on NOAA-derived phenology, MODIS-derived phenology and NCEP/NCAR reanalysis meteorological data. As we previously mentioned, in the Northern Hemisphere, phenology is primarily determined by temperature in the months period preceding the event (e.g., bud-burst, flowering, and leaf out), and higher pre-season temperatures may advance the BGS [2], [11], [17], [43], [44]. In this study, we determining the correlation analyses between the BGS and spring temperature (March, April and May), EGS and autumn temperature (August, September and October), BGS and precipitation, EGS and precipitation in the main ecozones of phenological changes. Fig 8 presents the statistical distributions of the correlation coefficients derived from Pearson correlation analyses between the interannual variations of BGS, EGS and of the climatic data (spring temperature, autumn temperature and mean annual precipitation). The presentation of the pixels were restricted to those eco-zones (Ba, Bb, BM, TeBSK, TeDc) only where the RPS exhibited a significant trend in the period 1982-2012. About 75% of all pixels showed a negative correlation with significance at p<0.05 level between the interannual variations of BGS and of the spring temperature, with Pearson's correlation coefficient ranging from -0.97 to -0.20, and histogram peak ranging from 0.4 to 0.6, which indicates that there is a strong negative correlation between BGS and spring temperature. About 45% of all pixels exhibited a positive correlation with significance at p<0.05 level between EGS and autumn temperature, with Pearson's correlation coefficient ranging from 0.2 to 0.95, and histogram peak ranging from 0.1 to 0.4, which indicates that there is a strong positive correlation between EGS and autumn temperature.This is consistent with previous studies that the earlier BGS and delayer EGS are the function of mean or accumulated temperatures over the Northern Hemisphere [2], [11], [45], [46]. Like the distribution pattern of the correlation coefficients between BGS, EGS and temperature, about 65% of all pixels showed a negative correlation between BGS Variability of Vegetation Phenology over the Northern Hemisphere and mean annual precipitation, with Pearson's correlation coefficient ranging from 0.2 to 0.90, and histogram peak ranging from 0.2 to 0.4. However, overall EGS and mean annual precipitation did not exhibit a significant correlation. It should be noted that comparing the correlation coefficient of vegetation phenology and temperature with that of phenology and precipitation, we can see the correlation between temperature and phenology is much higher than that between precipitation and phenology, displaying that temperature is a major determinant of changes in vegetation phenology. This is consistent with previous studies reporting the cycles of vegetation in temperate zones are primarily controlled by temperature [2], [17], [44], [47][48].
Even though temperature is considered as the major determinant for vegetation phenological change, little is known on vegetation phenological responses to temperature changes [44], [49][50]. Moreover, what is the phenological change per unit temperature. To further examine the temperature-sensitivity of the BGS and EGS, we analyzed the statistical distribution pattern of BGS and EGS changes in relation to mean annnual temperature within the different Holdridge life eco-zones (Fig 9). In those eco-zones whose RPS showed a significant trend during 1982-2012, the vegetation phenology was strongly dependent on temperature, in general, which indicates a marked advancing trend in BGS and a delaying trend in EGS along increasing temperature gradient. Interestingly, in all eco-zones, the BGS and EGS exhibited the largest variation in the temperature range of 5-12°C, which shows the vegetation phenological response to global warming is more sensitive in this temperature zone. Furthermore, in Ba, TeDC and TeBSK eco-zones, the advancement in BGS was larger than the delay in EGS, which leads to the inference that vernal vegetation response to temperature is greater and faster than that in autumn in these eco-zones. In the BM eco-zone, by contrast, the delay in EGS is Variability of Vegetation Phenology over the Northern Hemisphere consistently larger than the advancement in BGS, which also leads to the inference that autumnal vegetation response to temperature is greater and faster than that in spring in this zone. These inferences require detailed ground verification.
We next examined the vegetation phenological changes associated with temperature change within the different Holdridge life eco-zones (Fig 10). We found in all cases, higher mean temperature change associated with earlier mean BGS, such as during 2002-2012, the two largest advancement in BGS were found within the Ba (−0.84 days yr −1 ) and Bb (−0.55 days yr −1 ) zones, whose standard deviations of temperature change were 0.58 and 0.56, respectively. During 1992-2002, the largest advancement in BGS was found within the TeDc zone (−0.46 days yr −1 ), whose standard deviation of temperature change was 0.53.During 1982-1992, the largest advancement in BGS was found within the TeDc zone (−0.71 days yr −1 ), whose standard deviation of temperature change was 0.52. This is consistent with the results of previous studies showing that higher temperature-sensitivity associated with earlier onset of spring phenological events [44], [51][52]. Furthermore, in all eco-zones, overall vegetation phenology revealed a longer phenological cycle with earlier BGS and delayed EGS, which represents an extension of the growth cycle. Specifically, when comparing BGS with EGS, although the advanced BGS and delayed EGS have been reported to be the key factors that cause prolonged vegetation growing seasons, but in our study, the earlier BGS was found to be a more important factor that causes the LGS changes during 1982-2012, which differs from the results of Jeong et al. [2].

Uncertainties in Remote Sensing Estimates on Phenological Metrics
Satellite-based techniques have been widely used for phenological monitoring and eco-environmental evaluation. However, there were many uncertainties in the results due to the temporal, spatial, and ecological complexity of vegetation biochemical processes. In this study, the uncertainties in the phenological results are twofold: uncertainties related to the different datasets and uncertainties related to the methods for determing phenological metrics.
1. Uncertainties related to the different datasets. In this paper, the uncertainties related to datasets include five ways: satellite drift, sensor differences, calibration uncertainties, atmopsheric effects and sub-pixel clouds [10], [25][26][27]. For example, sensor differences, such as NOAA AVHRR sensor or MODIS sensor, with different sensor characteristics and spatial resolution (e.g., AVHRR and MODIS have different bandwidths and central wavelengths, different number of bands and radiometric resolution; and MODIS sensor has more better onboard geometric and radiometric calibrations than AVHRR sensor), can affect the quantification of phenology shifts. Some previous studies also show that the phenology paremeters (e.g., BGS, EGS) from the different datasets were significantly different over the nothern high lattitudes [20], [23][24], [53]. Furthermore, some studies have found that the newer generation sensors such as MODIS and SPOT performed better than AVHRR for investigating the phenology in alpine vegetation [20], [24], [54]. We found similar results with AVHRR and MODIS sensor, for instance, the BGS and EGS values appear big differences in both satellite sensors (see Figs 6 and 7), but the trends of BGS and EGS show generally similar spatio-temporal patterns and good agreement between the two sensors. To reduce the uncertainty of the results, in this paper the whole study period was divided into three periods (1982-1992, 1992-2002 and 2002-2012, respectively), in each period we only use the same dataset for trend analysis, thus these disagreements between AVHRR and MODIS sensors did not influence the results reported in this paper.
2. Uncertainties related to the methods for determing phenological metrics. On a hemispherical scale, satellite-based studies have been demonstrated to be useful for phenological monitoring. However, there are many uncertainties in the results when determing phenological metrics due to the temporal, spatial, and ecological complexity of phenological processes. Over the past few years, several models have been developed for detecting the phenological metrics by using satellite data, including NDVI thresholds, largest NDVI increase, backwardlooking moving averages, and fitting logistic functions, but each model has its uncertainty for estimating the phenological metrics. Overall, comparative studies have shown that the Midpoint pixel threshold model based on variations of an NDVI ratio is the one of the most consistent methods for determing phenological metrics over a variety of ecosystems and corroborates with field-based observational data [7], [55]. Thus, in order to determine the global phenological metrics based on AVHRR and MODIS data, we used the Midpoint pixel threshold method in this study. Furthermore, in response to concerns over the validity of satellite derived phenology, our research adopted a strict approach to the selection of pixels for analysis through statistical tests of significance at different levels.

Influencing Factors of Phenological Shifts
Many factors affect phenological shift patterns, such as topographic conditions, meteorological conditions, human disturbance and structural geology, among other factors. As far the global vegetation phenology, as reported by the IPCC's AR4 [56], many natural environmental factors strongly influence phenological shifts, such as temperature, precipitation/soil moisture, photoperiod, winter chilling. As we previously mentioned, in the Northern Hemisphere, phenology is primarily determined by temperature, specifically, alpine plant species are considered to be the most sensitive vegetation ecosystem in response to rising temperature. It should be noted that different plant species have different temperature requiements for development, also different plant species have different responses to temperature change. For instance, Jin et al. [57] reported that on the Tibetan Plateau, 1°C increase in soil temperature will advance the BGS by 4.6 to 9.9 days, and postpone the EGS by 7.3 to 10.5 days among the different plant species. Meanwhile, our study also indicates that there were different plant phenological responses to temperature within different eco-zones. In particular, our study further revealed that higher temperature warming related to earlier onset of BGS. Even though temperature clearly played the most important role in vegetation phenology at the high latitudes, precipitation/soil moisture was also an important factor for vegetation phenology. Especially in arid and semi-arid regions, precipitation/moisture was a more important factor in vegetation phenology. For example, Prieto et al. [58] found that flowering time was strongly influenced by precipitation patterns in Mediterranean shurbs. Lesica et al. [59], Broich et al. [60] and Krishnaswamy et al. [15] also reported that precipitation and moisture may play a more important role in phenological cycles in arid and semi-arid systems. In the high-latitude regions snowmelt is also an important factor affecting alpine vegetation phenology, especially in some cold regions, the time of snowmelt is often considered to be the observable indicators of phenological events in tundra plants [61]. Some recent studies have found that the spatial pattern of vegetation growing season mainly corresponded to the snowmelt gradient and plant species in the alpine ecosystem [62][63]. Jeganathan et al. [42] also found some link between the trend in the onset of snowmelt and inter-annual variation in BGS and EGS. Interestingly, in some snow-rich regions, such as the Kola Peninsula in Russia, and outer Troms County in northern Norway, mountain birch and herbaceous plants may have green leaves before snowmelt in spring when there is enough light through the snow [64][65]. Overall, in some high-latitude alpine areas, timing of snowmelt and temperature are the two most important factors controlling phenology.
In addition to temperature and precipitation, the photoperiod and winter chilling also have an influence on vegetation life cycle at the high latitudes. Although it is commonly assumed that winter chilling plays a important role in regulation of phenological phases, in fact, the interactions between chilling and photoperiod might play a more important role in controlling phenological phases at the high latitudes. For example, some recent experiments suggest that phtotoperiod might play an important role in controlling budburst and scenescence [66]. Some previous studies found that warmer temperatures during the winter will slow the dormancy breaking process by delaying the fulfillment of the chilling requirement. However, warm temperatures during the spring will advance the plant life growth by accelerating the accumulation of heat. [45], [52], [67]. Especially for humid extratropical areas, phenology is more controlled by winter chilling, photoperiod and temperature [68][69].

Conclusions
During the period from 1982 to 2012, most of the high-latitude regions experienced significant trends in advancing of BGS, delaying of EGS and prolonging of LGS, but no significant trends in PGS. Comparing the period of 2002-2012 with the periods of 1982-1992 and 1992-2002, the spatial pattern of RPS shows a more significant trend in advancing of BGS, delaying of EGS and prolonging of LGS during 2002-2012, overall shows a trend of accelerating change. The correlation coefficients between climate and phenological events for the Holdridge life ecozones from 1982 to 2012 showed that temperature was a major determinant of changes in vegetation phenology, with a significant advanced trend in BGS and a delayed trend in EGS along an increased temperature gradient. However, the phenological response of vegetation to temperature varied across different eco-zones. This study has shown that Midpointpixel threshold method is useful, and the trends in this study are generally in agreement with trends in previous studies. However, there are considerable uncertainties associated with the data and methods. To reduce uncertainies, it is clearly needed to design a biogeochemical model that integrates satallite-derived parameters, plant types, natural environmental factors and groundbased observations in future research.