Discrepancies in vegetation phenology trends and shift patterns in different climatic zones in middle and eastern Eurasia between 1982 and 2015

Abstract Changes in vegetation phenology directly reflect the response of vegetation growth to climate change. In this study, using the Normalized Difference Vegetation Index dataset from 1982 to 2015, we extracted start date of vegetation growing season (SOS), end date of vegetation growing season (EOS), and length of vegetation growing season (LOS) in the middle and eastern Eurasia region and evaluated linear trends in SOS, EOS, and LOS for the entire study area, as well as for four climatic zones. The results show that the LOS has significantly increased by 0.27 days/year, mostly due to a significantly advanced SOS (−0.20 days/year) and a slightly delayed EOS (0.07 days/year) over the entire study area from 1982 to 2015. The vegetation phenology trends in the four climatic zones are not continuous throughout the 34‐year period. Furthermore, discrepancies in the shifting patterns of vegetation phenology trend existed among different climatic zones. Turning points (TP) of SOS trends in the Cold zone, Temperate zone, and Tibetan Plateau zone occurred in the mid‐ or late 1990s. The advanced trends of SOS in the Cold zone, Temperate zone, and Tibetan Plateau zone exhibited accelerated, stalled, and reversed patterns after the corresponding TP, respectively. The TP did not occurred in Cold‐Temperate zone, where the SOS showed a consistent and continuous advance. TPs of EOS trends in the Cold zone, Cold‐Temperate zone, Temperate zone, and Tibetan Plateau zone occurred in the late 1980s or mid‐1990s. The EOS in the Cold zone, Cold‐Temperate zone, Temperate zone, and Tibetan Plateau zone showed weak advanced or delayed trends after the corresponding TP, which were comparable with the delayed trends before the corresponding TP. The shift patterns of LOS trends were primarily influenced by the shift patterns of SOS trends and were also heterogeneous within climatic zones.

Traditional ground-based observations are widely used for the species level and for small-scale areas because detailed information can be documented at the observation; however, these observations are limited by the location and number of observation sites and spatial scopes (Cleland et al., 2007;Studer, Stöckli, Appenzeller, & Vidale, 2007;Wu & Liu, 2013;Zhu et al., 2012). Remote sensing, in another way, provides continuously spatial and temporal informations on a variety of surface vegetation at regional or global scales, and numerous studies have investigated the variations of vegetation phenology in different regions over different time periods based on satellite-measured Normalized Difference Vegetation Index (NDVI) datasets (Fu et al., 2014;Myneni, Keeling, Tucker, Asrar, & Nemani 1997;Shen et al., 2018;Wang et al., 2016;White et al., 2009;Wu & Liu, 2013;Zhang, Tarpley, & Sullivan 2007).
Most previous studies reported a consistent observation of a significant advanced trend in start date of vegetation growing season (SOS) at the middle and high latitudes of the Northern Hemisphere during the 1980s and 1990s (Chen, Hu, & Yu 2005;Myneni et al., 1997;Tucker et al., 2001;Zhou et al., 2001). With the accumulation of remote sensing data, some studies based on long-term GIMMS NDVI dataset found a weakened, even delayed trend of SOS occurring since the late 1990s in many temperate regions of the Northern Hemisphere, such as East Asia, the Tibetan Plateau, and temperate China (Fu et al., 2014;Jeong, Ho, Gim, & Brown 2011;Piao, Cui, et al., 2011b;Wu & Liu, 2013;Yu, Luedeling, & Xu 2010). However, Wang et al. (2016) revealed that the SOS was advanced throughout the study period spanning from 1982 to 2012 in high latitudes of the Northern Hemisphere, from 55°N to 75°N. As the third pole of the Earth, the Tibetan Plateau is an opener and regulator of climate change in the Northern Hemisphere, and its SOS changes have always been the focus of debate in phenology studies. Based on the GIMMS NDVI during the period spanning 1982, Yu et al. (2010 found that the advanced trend in SOS has significantly reversed after the mid-1990s under increasing winter and spring temperatures in Tibetan Plateau. Subsequently, by integrating three different NDVI datasets, Zhang, Zhang, Dong, and Xiao (2013) revealed that the SOS in Tibetan Plateau showed a continuous advanced trend between 1982 and 2011. However, this result was later questioned by the study from Shen et al. (2013), who posited that the continuously advanced trend in SOS estimated by  is due to the lack of sufficient evidence and that the result is probably caused by not correcting for the changed pregrowth NDVI. In comparison to the debates on SOS change, the studies about the end date of vegetation growing season (EOS) were less reported, particularly on regional scales. EOS is also a critical variable in determining the length of vegetation growing season (LOS; Garonna et al., 2014) and regulates the balance of energy flow and material cycling in global and regional ecosystems (Estiarte & Peñuelas, 2015;Piao, Friedlingstein, Ciais, Viovy, & Demarty 2007;Richardson et al., 2013). A delayed trend in EOS during the past three decades has been detected at the global and the Northern Hemisphere scales (Garonna et al., 2014;Liu, Fu, Zhu, et al., 2016a;Zeng, Jia, & Epstein 2011), and the trend has weakened in several regions of the Northern Hemisphere during the 1990s and 2000s (Jeong et al., 2011;Yang, Guan, Shen, Liang, & Jiang 2015;Zhao et al., 2015). However, the shifting pattern of EOS trend is still largely unclear among those studies.
In this paper, we monitored the variations of vegetation phenology (SOS, EOS, and LOS) in different climatic zones in the middle and eastern Eurasia region from 1982 to 2015. The middle and eastern Eurasia region has a vast territory, covering a wide range of ecosystems and climates, and including rich biodiversity. In the study area, the climate change on the Tibetan Plateau has a great impact on the climate of the Eurasia and even the Northern Hemisphere (Li et al., 2003;Piao, Cui, et al., 2011b;Shen, 2011), and vegetation phenology in the Tibetan Plateau is highly sensitive to climate change (Yu et al., 2010 (Fu et al., 2014;Jeong et al., 2011;Shen et al., 2013;Wang et al., 2015;Wu & Liu, 2013;Yu et al., 2010;Zhang et al., 2013). Therefore, the trend of vegetation phenology and its response to climate change should be further investigated with the accumulation of satellite observation data. At present, a new generation of GIMMS NDVI 3g dataset has been extended to 2015 Shen et al., 2018;Xia et al., 2018), and its quality is better than the previous version (Dai et al., 2018;Zhou et al., 2019). Longer-term data are more useful for us to understand the relationship between phenology and climate change. We hypothesized that there are divergent responses of vegetation phenology to climate change in different regions related to the different climate backgrounds. Therefore, we applied the Polyfit-Maximum (P-M) method to estimate the SOS, EOS, and LOS in the middle and eastern Eurasia region using long-term satellite NDVI data from 1982 to 2015. The aims of this study were (a) to quantify the long-term trends of vegetation phenology for the entire study area and its four climatic zones, and (b) to clarify the shift patterns of phenology trends in different climatic zones.

| Study area
We selected the middle and eastern Eurasia region (73°-170°E and 23.5°-85°N; Figure 1) as the study area, where seasonal changes of vegetation are clear and where the satellite-derived NDVI dataset is least affected by the solar zenith angle in the middle and high latitudes of the Northern Hemisphere (Cong et al., 2013;Jeong et al., 2011;Slayback, Pinzon, Los, & Tucker 2003). Our study excluded the pixels of croplands because their vegetation phenology is intensively influenced by anthropogenic activities (Liu, Fu, Zeng, et al., 2016b;White, Thornton, & Running 1997;Zhang, Friedl, Schaaf, & Strahler 2004). Areas covered by evergreen vegetation were also excluded from our study because of the lack of obvious seasonality in surface greenness (Yang et al., 2015). Moreover, this study removed the pixels from nonvegetation areas, water bodies, and areas with an annual mean NDVI value during 1982-2015 that is smaller than 0.1 in order to reduce the influence of bare land and sparse vegetation on the interannual curve of NDVI (Jeong et al., 2011;Zhou et al., 2001).
To analyze the spatial pattern of phenology trends, the study area was divided into four climatic zones, including Cold zone

| Data collection and generation
The study used the latest and longest release of GIMMS NDVI 3g version 1 dataset, which was provided by the Global Inventory  (Cong et al., 2013;Pinzon & Tucker, 2014). The new generation of the NDVI dataset added percentile data and recovered NDVI negative values of snow-covered F I G U R E 1 Location of the study area and four climatic zones as well as the spatial distribution of vegetation in the middle and eastern Eurasia region. The bottom right inset shows the location of the four climatic zones regions in the winter at high latitudes in the Northern Hemisphere.
Vegetation type data in our study were obtained from a land cover classification map, which was provided by the Joint Research Centre of the European Commission under the project of Global Land Cover 2000 (GLC2000; http://bioval.jrc.ec.europa.eu/produ cts/glc20 00/ glc20 00.php; Yang et al., 2015). The GLC2000 was made from the SPOT4\VEGETATION data by the Land Cover Classification System of Food and Agriculture Organization and included 22 land cover types (Bartholomé & Belward, 2005). The land classification map with a spatial resolution of 1 km was reprojected and resampled, to match NDVI 3g dataset.

| Determining the SOS and EOS
Before detecting vegetation phenology, we removed the cloud and snow contaminations for each annual time series of NDVI data.
First, we distinguished the pixels with potentially covered by cloud or snow using the flag data, which was calculated from the percentile data. These contaminated values were replaced by the temporally nearest valid NDVI values (Liu, Fu, Zhu, et al., 2016a). Second, the snow contamination was further minimized due to uncertainty of the flag data (Zhang et al., 2007). We screened out a time series of annual smallest NDVI values from 1982 to 2015 for each pixel. The median value of the smallest NDVI series in each pixel was extracted as surface background NDVI value, which was applied to replace the smaller values in annual NDVI time series (Zhang et al., 2007).
In our study, we adopted the P-M method, which has been widely applied for the extraction of vegetation phenology in middle and high latitudes of the Northern Hemisphere (Fu et al., 2014;Jeong et al., 2011;White et al., 2009)  where t is time (temporal resolution of 15 days). Second, the NDVI(t), corresponding to the maximum NDVI ratio at time t, is used as the threshold for detecting the SOS Wu & Liu, 2013). The NDVI(t + 1) at time (t + 1), corresponding to the minimum NDVI ratio at time t, is used as the threshold for detecting the EOS (Liu, Fu, Zeng, et al., 2016b;Yang et al., 2015). Finally, the NDVI time series data at 15 day intervals from January to September (for SOS) and from July to December (for EOS) for each year are fitted by a sixth-degree polynomial function Wu & Liu, 2013;Yang et al., 2015;Equation 2), to reconstruct the NDVI time series with a temporal resolution of 1 day. The SOS and EOS are determined from the reconstructed daily NDVI time series curves based on the corresponding threshold for each pixel.
where x is the day (Julian calendar). The regression coefficients (a 0 , a 1 , a 2 , a 3 …a n ) are determined by the least square regression. In addition, the EOS subtracts the SOS of the corresponding pixel, as well as the LOS for each pixel.

| Statistical analysis
The long-term trends in SOS, EOS, and LOS for the entire study to eliminate statistical uncertainties due to the first and last data points and individual abnormalities in the phenological time series (Wu & Liu, 2013).
where t is year; y is value of the SOS, EOS, and LOS time series in each climatic zone and the entire study area; α represents the TP in the SOS, EOS, and LOS trends; λ 0 , λ 1 , and λ 2 are regression coefficients; and ε is the residual random error. λ 1 and λ 1 + λ 2 represent the piecewise linear trends (SOS, EOS, and LOS) before and after the TP. Least square regression is applied to determine α and other coefficients. To assess the necessity of introducing TP, a t test was used to verify the following null hypothesis: "λ 2 is not significantly different from zero" Wu & Liu, 2013). A p value <0.05 was considered to be significant. Using simple linear regression, we also calculated linear trends in the SOS, EOS, and LOS before and after the TP.

| Trends in SOS and its discrepancies among different climatic zones
The mean SOS in the middle and eastern Eurasia region experienced a significant advance from 1982 to 2015 (R 2 = 0.42, p < 0.001),

| Trends in EOS and its discrepancies among different climatic zones
The mean EOS during 1982-2015 over the entire study area showed a slow delayed trend with an average rate of 0.07 days/ year (R 2 = 0.29, p = 0.001; Figure 4).

| Trends in LOS and its discrepancies among different climatic zones
The mean LOS over the middle and eastern Eurasia region significantly increased (R 2 = 0.54, p < 0.001) by a rate of 0.27 days/year over the past 34 years ( Figure 6). Furthermore, the longest LOS oc- TPs were also detected in the Cold-Temperate zone, Temperate zone, and Tibetan Plateau zone, and they occurred in 1989, 1997, and 1991 (Table 1)

| Changes in SOS trends
The  Barichivich et al. (2013). This finding is considered to be probably associated with spatial asymmetry of spring temperatures increasing in the Northern Hemisphere because spring phenology is a sensitive indicator of the response of terrestrial ecosystems to climate change (Cohen, Furtado, Barlow, Alexeev, & Cherry 2012;Menzel & Fabian, 1999;Menzel et al., 2006;Schwartz, 1998). Therefore, the regional differences of SOS trends should be further discussed on the condition that vegetation responses are associated with climate factors. Long-term trends of SOS in the four climatic zones illustrated that the higher latitude areas had a higher advanced rate, which is consistent with the finding of Post, Steinman, and Mann (2018), who also found rapidly advanced trend in SOS at higher latitudes of the Northern Hemisphere. Previous studies revealed that the advance in SOS was related to preseason temperature warming (Chen et al., 2005;Myneni et al., 1997; and that warming in higher latitudes of the Northern Hemisphere was stronger than that in lower latitudes (Hassol, 2004;Post et al., 2018). Warmer springs could easily thaw the topsoil and  (Chen et al., 2005;Cong et al., 2013;Menzel et al., 2006;. The weakened advance of SOS in the Temperate zone could be explained by the reversal or stalling of spring temperatures since the late 1990s (Jeong et al., 2011;Wang et al., 2011;Wu & Liu, 2013). Previous studies showed that global warming trends have decelerated over the past decade (Cane, 2010;Chen & Tung, 2014;Solomon et al., 2010). In contrast, warming in the spring at higher latitudes of the Northern Hemisphere has lasted over the recent decades (Kaufman, Schneider, McKay, Ammann, & Bradley, 2009;Li, Wang, Zhang, Li, & Zang, 2017;Serreze & Barry, 2011). Therefore, the SOS in Cold zone and Cold-Temperate zone still significantly advanced after the mid-1990s. In addition, the tundra vegetation in Cold zone is more sensitive to rises in temperature (Chapin, Eugster, Mcfadden, Lynch, & Walker, 2000;Hinzman et al., 2005).
Compared with other climatic zones, the reversed trend of SOS in the Tibetan Plateau zone is a unique phenomenon in Eurasia.
Although the Tibetan Plateau zone with its lower annual mean temperatures is similar to the Cold-Temperate zone, and while the vegetation in the Tibetan Plateau zone is similar to the Cold zone with the predominance of tundra vegetation, the Tibetan Plateau zone shows a very different shift pattern of spring phenology. Yu et al. (2010) concluded that rapid warming in the winter induced a lack of chilling in the Tibetan Plateau zone, which could lead to a delay in SOS. Zhang et al. (2007) also found that chilling accumulations in the winter decreased due to the reduction of chilling days and, as a result, that SOS was delayed starting at 40°N and moving southward in North America. However, this explanation of the delay in SOS was later questioned by the following studies. Chen, Zhu, Wu, Wang, and Peng (2011) believed that the combined effect of the thawing-freezing process, grassland degradation, and climate changes, but not the unitary effect of spring and winter warming, delayed the SOS.
However, a recent study from Shen, Piao, Cong, Zhang, and Jassens (2015) found that the delay in SOS mainly occurred in the southwestern region of the Tibetan Plateau zone and was due to the decrease in spring precipitation. The Tibetan Plateau zone has a special plateau climate and is strongly influenced by the East Asia monsoon events. A strong 1997/98 El Nino event significantly advanced the SOS in the Tibetan Plateau zone in 1998 (Klein et al., 2014;, and the SOS began to delay after 1998 despite the continued warming in winter and spring. Therefore, climate change and other factors (e.g., human activity) in the Tibetan Plateau zone can affect the changes in SOS. All of these issues require further study from different perspectives, such as dividing the Tibetan Plateau zone into multiple subunits, paying more attention to precipitation in the western part and strengthening the combination of ground observations and remote sensing observations.

| Changes in EOS trends
Our study estimated that the EOS over the middle and eastern Eurasia region was delayed on average by a rate of 0.07 days/year from 1982 to 2015. Previous studies also documented that the EOS experienced a delayed trend in regions of the Northern Hemisphere, such as North America (Zhu et al., 2012), Eurasia (Li et al., 2017), Europe (Garonna et al., 2014;Stöckli & Vidale, 2004), temperate China (Liu, Fu, Zeng, et al., 2016b;Yang et al., 2015), and the Tibetan Plateau (Cong, Shen, & Piao 2017;Jin, Zhuang, He, Luo, & Shi 2013 Hemisphere. It may be related to asynchronous changes in summer and autumn temperatures between different regions (Cohen et al., 2012). Previous studies have proven that autumn phenology is severely affected by summer and autumn temperature (Che et al., 2014;Yang et al., 2015). On climatic zone scales,  (Yang et al., 2015).
Studies suggested that warmer preseason (summer and autumn) temperatures delayed the EOS (Liu, Fu, Zhu, et al., 2016a;Yang et al., 2015). The preseason warming trends decelerated or even reversed in recent decades (Cane, 2010 Zeng et al. (2011) showed that Arctic, Eurasian, and American regions had different rates of LOS increase.

| Changes in LOS trends
Advanced SOS and delayed EOS lengthened LOS in the middle and eastern Eurasia region over the past 34 years. However, the advanced rates of SOS in the entire study area and over the four climatic zones were higher than the delayed rates of EOS, except for the Tibetan Plateau zone. This result is consistent with the conclusion of Zhao et al. (2015), who reported that the advanced rate of SOS was higher than that of EOS in Asia. In addition, similar to the variations of SOS trends, LOS trends also varied among different climatic zones, and the extended rates increased with increasing latitude. This revealed that the SOS is a key factor that regulates the LOS changes in middle and eastern Eurasia over the past 34 years.
However, Zhu et al. (2012) concluded that the LOS was prolonged mainly due to the delay in EOS in North America during the period spanning from 1982 to 2006.
The increased rate of LOS slowed down in the Cold-Temperate zone during the period after the corresponding TP, whereas it was reversed in the Tibetan Plateau zone. The increased trend of LOS in the Temperate zone disappeared after the mid-1990s. The shift patterns of LOS trends in some regions of the Northern Hemisphere have been observed in previous studies (Jeong et al., 2011). This result was primarily due to the decrease and reversal of SOS trends or to the decrease and advance of the EOS trends during the period after the corresponding TP .

| CON CLUS IONS
In this study, we investigated the trend characteristics of SOS, EOS, and LOS in the middle and eastern Eurasia.

CO N FLI C T O F I NTE R E S T
None declared.