Vegetation Phenology in the Qilian Mountains and Its Response to Temperature from 1982 to 2014

: The vulnerability of vegetation ecosystems and hydrological systems in high-altitude areas makes their phenology more sensitive and their response to climate change more intense. The Qilian Mountains, an important geographic unit located in the northeastern Tibetan Plateau (TP), has experienced the more signiﬁcant increases in temperature and precipitation in the past few decades than most areas of the TP. However, under such intense climate change, the temporal and spatial differences in phenology in the Qilian Mountains are not clear. This study explored the spatial and temporal heterogeneity of phenology in the Qilian Mountains from 1982 to 2014 and its response to three temperature indicators, including the mean daily temperature ( T mean ), mean daily daytime temperature ( T max ), and mean daily nighttime temperature ( T min ). The results showed that (1) as the altitude rose from southeast to northwest, the multiyear mean of the start of the growing season (SOS) was gradually delayed mainly from 120 to 190 days, the multiyear mean of the end of the growing season (EOS) as a whole was advanced (from 290 to 260 days), and the multiyear mean of the length of the growing season (LGS) was gradually shortened (from 150 to 80 days). (2) In general, there was an advanced trend in the annual average SOS (0.2 days per decade), a delayed trend in the annual average EOS (0.15 days per decade), and an extended trend in the annual average LGS (0.36 days per decade) over the study period. However, there has been no signiﬁcant phenological trend in recent years, especially for the SOS after 2000 and the EOS and LGS after 2003. (3) Higher preseason temperatures led to an advanced SOS and a delayed EOS at the regional scale. Moreover, the SOS and EOS were more triggered by T max than T min and T mean . The LGS was signiﬁcantly positively correlated with annual mean temperature ( r = − 0.82, p < 0.01).


Introduction
Vegetation phenology, the periodic life activity of plants [1,2], is a sensitive indicator of climate change. In recent decades, climate change has led to phenological changes in many plants, especially in temperate and northern regions [3,4], resulting in small-area changes in plant activity within the community and large-area influences of overall land surface processes, such as surface energy flux, carbon budget, and regional climate [5]. It is, therefore, critical to monitor vegetation phenology and explore its response to climate to elucidate the role of phenology in terrestrial ecosystem models and protect terrestrial ecosystems.
In high-altitude areas, the vulnerability of vegetation ecosystems and hydrological systems produces a more sensitive and intense response of phenology to climate change. Various studies have explored the changes in the phenology of the TP [8][9][10], but it is still unclear in the Qilian Mountains. The area of the Qilian Mountains is approximately one-tenth that of the TP, so regional studies on the TP fail to provide detailed information on the phenological changes in the Qilian Mountains. Moreover, a strong environmental gradient results in spatial heterogeneity in phenology and diverse phenological responses to climatic variables on the TP. Shen et al. [8] found a significantly delayed green-up date in the southwestern TP from 2000 to 2011, while other regions had an advanced trend. However, the preseason temperature increase in the southwestern TP region was more significant than that in other areas. Li et al. [11] explored the phenological trends in the Huangyuan, Menyuan, and Delhi areas of the Qilian Mountains using data from ground phenology observations. Zhou et al. [12] analyzed the spatial and temporal patterns of the Heihe River headwaters, which are located in the northeastern Qilian Mountains. However, our understanding of the spatial and temporal variation in phenology in the Qilian Mountains is still unclear. Hence, small-area phenological changes on the TP need further investigation.
Because of the unique geographical conditions and socioeconomic characteristics in the Qilian Mountains, it is important to explore the spatial and temporal heterogeneity of phenology and its response to climate. Piao et al. [13] found that phenology would affect the local net primary productivity (NPP) and forage yield. Meanwhile, animal husbandry, for which much foraging is needed, is a key economic industry in the Qilian Mountains. In addition, as an important ecological function zone, the Qilian Mountains have blocked the confluence of deserts effectively and maintained the ecological balance of desert and oasis. The Qilian Mountains are the source of the Shiyang, Heihe, Shule River in the northwest arid region and the runoff formation area. Meanwhile, it provides abundant water Various studies have explored the changes in the phenology of the TP [8][9][10], but it is still unclear in the Qilian Mountains. The area of the Qilian Mountains is approximately one-tenth that of the TP, so regional studies on the TP fail to provide detailed information on the phenological changes in the Qilian Mountains. Moreover, a strong environmental gradient results in spatial heterogeneity in phenology and diverse phenological responses to climatic variables on the TP. Shen et al. [8] found a significantly delayed green-up date in the southwestern TP from 2000 to 2011, while other regions had an advanced trend. However, the preseason temperature increase in the southwestern TP region was more significant than that in other areas. Li et al. [11] explored the phenological trends in the Huangyuan, Menyuan, and Delhi areas of the Qilian Mountains using data from ground phenology observations. Zhou et al. [12] analyzed the spatial and temporal patterns of the Heihe River headwaters, which are located in the northeastern Qilian Mountains. However, our understanding of the spatial and temporal variation in phenology in the Qilian Mountains is still unclear. Hence, small-area phenological changes on the TP need further investigation.
Because of the unique geographical conditions and socioeconomic characteristics in the Qilian Mountains, it is important to explore the spatial and temporal heterogeneity of phenology and its response to climate. Piao et al. [13] found that phenology would affect the local net primary productivity (NPP) and forage yield. Meanwhile, animal husbandry, for which much foraging is needed, is a key economic industry in the Qilian Mountains. In addition, as an important ecological function zone, the Qilian Mountains have blocked the confluence of deserts effectively and maintained the ecological balance of desert and oasis. The Qilian Mountains are the source of the Shiyang, Heihe, Shule River in the northwest arid region and the runoff formation area. Meanwhile, it provides abundant water resources to Hehuang Valley and Hexi corridor. The Qilian Mountains are rich in ecological diversity and have diverse vegetation types. Generally, with the elevation increasing from southeast to northwest, the vegetation types are cultivated vegetation, needleleaf forest, broadleaf forest, brush, meadow, grassland, and alpine vegetation ( Figure A1 in Appendix A).
There are three main methods to extract vegetation phenology: ground phenology observations, phenology models, and satellite-derived methods. Ground phenology observations can record specific plant phenology information, but only at the individual plant scale [14,15]. In addition, many phenology observation sites only focus on crops rather than natural vegetation, which cannot accurately capture phenology's response to climate change [16]. Phenology models use meteorological factors to simulate phenophases, but the plant species are limited, and the quantitative methods to evaluate the control of environmental factors need further study [17]. Satellite-derived methods, which depict the detailed vegetation growth process using high temporal resolution vegetation index time series derived from satellite remote sensing images, are the most widely used methods [18]. This method provides phenological metrics of various types of vegetation in a large area. Due to the limited number of available phenological sites and diverse vegetation types in the Qilian Mountains, in this study, the satellite-derived method was employed to detect annual seasonal phenology parameters.
The start of the growing season (SOS), the end of the growing season (EOS), and the length of the growing season (LGS) are the most widely used phenology parameters [19][20][21], which directly reflect the continuous decadal or interannual changes in phenology as a result of climate change. The SOS, extracted by the satellite-derived method, refers to the date when most plants in the pixel begin to develop leaves, after which the plants grow rapidly. The EOS refers to the date when the majority of plants in the pixel begin to turn yellow. After that, plants no longer continue to grow and gradually die or enter a dormant stage. The LGS refers to the length of the growing season, which is numerically equal to the difference between the SOS and EOS. Although the individual phenological events derived from satellite-derived methods attract particular attention when exploring the relationship between phenology and preseason temperature [22], a systematic study of the SOS, EOS, and LGS is essential.
Among the climate driving factors that affect phenology, temperature shows a definite correlation with phenology. Ge et al. [23] found that there was a significant correlation between the SOS and spring temperature (correlation coefficient > 0.7) in China, and the correlation was stronger than that with precipitation. Cong et al. [9] found that the interannual variation in the regional average EOS was driven mainly by the autumn temperature, and precipitation and insolation had a weak impact on the EOS. It has been widely acknowledged that an increase in temperature would advance spring phenology and delay autumn phenology in most areas. Ma et al. [24] found that spring phenology started 4.93 days earlier in response to a temperature increase of 1 • C over the last three decades.
Generally, the mean daily mean temperature (T mean ) is used as the temperature indicator when exploring the association between phenology and temperature, which ignores the different potential responses of phenology to daily daytime temperature (T max ) and daily nighttime mean temperature (T min ). Piao et al. [4] found that the satellite-derived spring phenology was triggered by T max rather than T min and T mean in the Northern Hemisphere (>30 • N). Shen et al. [25] studied the vegetation green-up date on the TP and found a stronger negative correlation between the green-up dates and T min than with T max at the regional scale. In addition, spatial heterogeneity in the phenological response to T max and T min was also found in this study. Ding et al. [26] found little difference in the EOS response to T max and T min on the TP. However, for the Qilian Mountains, the three types of temperature's various impacts on phenology are still unclear, which is an obstacle for precise phenological modeling in this area.
With the aim of exploring the spatial and temporal heterogeneity of phenology, including the SOS, EOS, and LGS, in the Qilian Mountains from 1982 to 2014 and its response to three temperature indicators, including T mean , T max , and T min , in this study, we extracted the phenology parameters using the midpoint threshold method, followed by the analysis of the spatial heterogeneity and dynamics of the phenology using linear regression. Finally, using the correlation coefficient, we compared the relationship between the phenology and three temperature indicators.

Data Source and Pretreatment
The time series of the Normalized Difference Vegetation Index (NDVI) is able to reflect the changes in vegetation biophysical with time and the seasonal changes [27,28], which has been widely used to quantify the vegetation growth and phenology research [10,29]. An advanced very high-resolution radiometer (AVHRR) NDVI dataset, developed by the global inventory modeling and mapping studies (GIMMS) group of the National Aeronautics and Space Administration (https://ecocast.arc.nasa.gov/data/pub/gimms/ 3g.v1/), has a spatial resolution of 1/12 • and covers 1982 to 2014 with a biweekly temporal resolution and has been corrected for various disturbances, such as terrain and volcanic ash. Owing to the consistent observation in a large area provided by the AVHRR, this dataset has been used in various phenology studies, including the alpine areas, which are similar to the Qilian Mountains [10,23]. We then selected the GIMMS NDVI dataset in the Qilian Mountains from 1982 to 2014 as the data source. Pixels without vegetation or with sparse vegetation were excluded. In particular, pixels with over 0.15 annual maximum NDVI and over 0.1 annual mean NDVI were used in this study [23]. The resultant NDVI dataset of the Qilian Mountains has 1827 pixels for 33 years (1982-2014), as shown in Figure 1b. To validate the satellite-derived vegetation phenology, field data of phenology sites in the Qilian Mountains were collected from the China Meteorological Administration (CMA) [30] (Figure 1a).
The China Meteorological Forcing (CMF) Dataset (http://westdc.westgis.ac.cn/zhhans/data/7a35329c-c53f-4267-aa07-e0037d913a21/), which covers 1979 to 2015 with a 0.1 • spatial resolution and a temporal resolution of 3 h [31], was used to explore the correlation between phenology and temperature. The forcing dataset was developed by the Data Assimilation and Modeling Center for Tibetan Multispheres, Institute of Tibetan Plateau Research, Chinese Academy of Sciences. In this study, near-surface air temperature at 2 m with a temporal resolution of 3 hours was downloaded, followed by transformation to T mean , T max , and T min . Phenology is significantly affected by the temperature of the preceding months before the mean phenophases [12,23], which is defined as the preseason. In this study, if the mean phenophases appear in the first half of one month, the preseason time refers to the two months before the current month; if the mean phenophases occur in the latter half of the month, the preseason time refers to the current month and previous month. Then, the mean daily temperature, maximum daily temperature, and minimum daily temperature during the preseason time were regarded as T mean , T max , and T min , respectively.

Extraction of Phenology
The phenology parameters were extracted via the following steps. The first step was to alleviate the negative impact of noise by using smoothing methods. In this study, asymmetric Gaussian functions [32] were used to address the noise and fit the GIMMS NDVI time series, since this method achieved the best fitting performance in the Qilian Mountains by comparison with two other smoothing methods, including the Savitzky-Golay filter [33] and double logistic smooth functions [34]. The core of the asymmetric Gaussian is to use a piecewise function to fit the NDVI time series as follows: where t refers to the time and f(t) refers to the NDVI value. The linear parameters c 1 and c 2 determine the base level and the amplitude, respectively. g(t; x 1 , x 2 , . . . , For Equation (2), x 1 determines the position of time variable t when the NDVI is at its maximum and minimum value. x 2 and x 3 determine the width and flatness of the right function. Similarly, x 4 and x 5 determine the width and flatness of the left function. The NDVI time series is fitted independently every year.
Based on the fitted NDVI time series, there are two main methods to detect phenology: (1) the threshold method and (2) the derivative method. The derivative method is sensitive to the selection of the fitting functions [8] as well as the noise in the time series. Meanwhile, the method may fail to determine the SOS and EOS when there are no abrupt or rapid changes during the growth stage [35]. The midpoint threshold method [21,23] is a typical threshold method, in which the SOS and EOS are defined as the days when the NDVI is equal to half of the range between the maximum and minimum of the NDVI time series (Figure 2). The derivative method is sensitive to different fitting functions. In contrast, the midpoint threshold method has a strong consistency for different fitting functions. We, thus, used the midpoint threshold method to extract the SOS and EOS, followed by the calculation of the LGS by the difference between the SOS and EOS. The phenology parameters are generally displayed as the day of the year (DOY) instead of the date. For example, January 25 can be instead of the 25 days of the year.
For Equation (2), x1 determines the position of time variable t when the NDVI is at its maximum and minimum value. x2 and x3 determine the width and flatness of the right function. Similarly, x4 and x5 determine the width and flatness of the left function. The NDVI time series is fitted independently every year.
Based on the fitted NDVI time series, there are two main methods to detect phenology: (1) the threshold method and (2) the derivative method. The derivative method is sensitive to the selection of the fitting functions [8] as well as the noise in the time series. Meanwhile, the method may fail to determine the SOS and EOS when there are no abrupt or rapid changes during the growth stage [35]. The midpoint threshold method [21,23] is a typical threshold method, in which the SOS and EOS are defined as the days when the NDVI is equal to half of the range between the maximum and minimum of the NDVI time series (Figure 2). The derivative method is sensitive to different fitting functions. In contrast, the midpoint threshold method has a strong consistency for different fitting functions. We, thus, used the midpoint threshold method to extract the SOS and EOS, followed by the calculation of the LGS by the difference between the SOS and EOS. The phenology parameters are generally displayed as the day of the year (DOY) instead of the date. For example, January 25 can be instead of the 25 days of the year.

Trend Analysis and Correlation Analysis
The linear regression method was used to analyze the phenology dynamics in this study. The variable t refers to one year, and p refers to the phenology date of this year. The linear regression between time and phenology was established as follows: When slope b is greater than 0, the phenology date increases with time, and the greater the slope is, the wider the change range. Conversely, the date of phenology decreases with time if the slope is lower than 0. The slope in Equation (3) was calculated by the least squares method as follows:

Trend Analysis and Correlation Analysis
The linear regression method was used to analyze the phenology dynamics in this study. The variable t refers to one year, and p refers to the phenology date of this year. The linear regression between time and phenology was established as follows: When slope b is greater than 0, the phenology date increases with time, and the greater the slope is, the wider the change range. Conversely, the date of phenology decreases with time if the slope is lower than 0. The slope in Equation (3) was calculated by the least squares method as follows: The correlation coefficient (r) and p-value (p) were calculated to analyze the association between phenology and temperature in this study. The correlation coefficient was calculated as follows: where x refers to the variable temperature, and y refers to the phenology parameters (SOS, EOS, and LGS). The significance of the correlation coefficient was also tested by the p-value. Additionally, the coefficients of determination (r 2 ) and p-values between the satellite-derived phenology and field data were used as indicators of the accuracy of extracted satellite-derived phenology in this study.

Spatial Variation in the Multiyear Mean Phenology in the Qilian Mountains
The satellite-derived vegetation phenology extracted in this study has displayed an agreement with field observations, as evidenced by r 2 > 0.5, p < 0.01 between the SOS and field data, r 2 = 0.44, p < 0.1 between the EOS and field data in Haibei phenology sites, and r 2 = 0.47, p < 0.01 between the EOS and corresponding field data in Huangyuan and Huzhu phenology sites, which ensures the reliability of our further analysis.
The spatial distribution of the multiyear mean phenology in the Qilian Mountains is shown in Figure 3. From northwest to southeast, the phenological parameters of the Qilian Mountains showed evident change. Overall, the multiyear mean SOS gradually advanced, the multiyear mean EOS gradually delayed, and the multiyear mean LGS gradually extended. Furthermore, there was a considerable difference between the phenology values of the northwest and southeast. Excluding the individual extreme values, the multiyear mean SOS in the northwestern Qilian Mountains was approximately 50 days later than that in the southeastern Qilian Mountains. Similarly, the multiyear mean EOS was 30 days earlier in the northwest than in the southeast. The multiyear mean LGS was shortened by 70 days in the northwest. The multiyear mean SOS in the Qilian Mountains mainly occurred between 130 days and 190 days, which accounted for 99.56% of the study region ( Figure 3a). From northwest to southeast, the multiyear mean SOS advanced gradually, ranging from 130 to 190 days, and the mean value was 159.24 days as a whole. In approximately 13.63% of the study region, the multiyear mean SOS was earlier than 150 days. These places were located in Ledu, Huangzhong, southern Huatong, and southern Datong areas on the southeast edge The multiyear mean SOS in the Qilian Mountains mainly occurred between 130 days and 190 days, which accounted for 99.56% of the study region ( Figure 3a). From northwest to southeast, the multiyear mean SOS advanced gradually, ranging from 130 to 190 days, and the mean value was 159.24 days as a whole. In approximately 13.63% of the study region, the multiyear mean SOS was earlier than 150 days. These places were located in Ledu, Huangzhong, southern Huatong, and southern Datong areas on the southeast edge of the Qilian Mountains, where they had lower elevations and the vegetation types are cultivated crops, cultivated crops, and grassland. The area where the multiyear mean SOS ranged from 150 to 160 days was concentrated in the low mountainous areas around Qinghai Lake, including Gonghe, Haiyan, Huangyuan, Huzhu, Menyuan, Tianzhu, and some northeastern areas. The vegetation of these areas is mainly brush and meadow. For the western Qilian Mountains with a high elevation, the multiyear mean SOS was between 160 and 170 days, including Delhi, Tianjun, Wulan, and southern Gangcha, Sunanyu, and Qilian Counties. The area with a multiyear mean SOS later than 170 days was mostly distributed in the high mountains on the northwestern edge, which accounted for 2.85% of the study region.
The multiyear mean EOS mainly occurred on days 260-290 of the year, and the average multiyear mean EOS was 276.26 days in the Qilian Mountains ( Figure 3b). From northwest to southeast, the multiyear mean EOS was gradually delayed as a whole, but the spatial differentiation of the EOS was not as significant as that of the SOS. In the high northwestern Qilian Mountains, the multiyear mean EOS was earlier than 270 days, including Tianjun, northern Gangcha, northern Qilian, northern Sunanyu, and Subei, where the grassland, meadow, and alpine vegetation are distributed there. In addition to these high mountainous areas, at the intersection of Huangzhong, Huzhu, and Datong, the multiyear mean EOS was also earlier than 270 days. This area is an important cultivated vegetation region. The multiyear mean EOS values ranging from 275 to 280 days were distributed in Tianzhu, Menyuan, southern Qilian, and the area in the southwest of the Qinghai Lake, which have lower elevations than those of the northwest. For northern Qinghai Lake, the multiyear mean EOS was between 280 and 285 days. Only in approximately 3.17% of the study region was the EOS later than 285 days.
The multiyear mean LGS was mainly concentrated in the period of 80-150 days, with a mean value of 117.02 days (Figure 3c). From northwest to southeast, the multiyear mean LGS extended gradually, showing a vertical differentiation law consistent with the SOS. The area where the multiyear mean LGS was longer than 120 days was around the southern and northeastern Qinghai Lake, with a lower elevation than the mean value of the study region. Other areas accounting for 55.06% of the study region showed a multiyear mean EOS shorter than 120 days, which was mainly distributed in the northwestern high mountains, where grassland, alpine vegetation, and large areas of meadow are distributed.
These results showed that the multiyear mean phenology parameters (SOS, EOS, and LGS) gradually changed with a gradient distribution from low to high elevation. The relationship between phenology and altitude was further analyzed in the discussion.

Spatial Variation of the Phenological Trend from 1982 to 2014 in the Qilian Mountains
As shown in Table 1, the Qilian Mountains have an advanced SOS, a delayed EOS, and an extended LGS during the period of 1982 to 2014 at the regional scale. The overall linear trend of the SOS was 0.2 days earlier per decade with a significant p-value < 0.01. The regional mean EOS was delayed by 0.15 days per decade. The regional mean LGS had a significantly extended trend of 0.36 days per decade. However, different phenological trends were extracted at the pixel scale, demonstrating an obvious spatial variation in phenological trends. More than 60% of the Qilian Mountains showed an advanced SOS trend, and the advanced trend was significant (p < 0.05) in 59.51% of the study region (Figure 4a). The areas with a delayed SOS from 1982 to 2014 accounted for 36% of the Qilian Mountains and were mainly distributed in Delhi, the junction of Menyuan and Tianzhu, and the northwestern area of Qinghai Lake, where the vegetation types in these areas are mainly grassland and brush.  Approximately 66.93% of the study region had a delayed EOS trend. The EOS was significantly delayed from 1982 to 2014 in 58.74% of the area and was mainly distributed in the central area of the Qilian Mountains (Figure 4b). The remaining pixels, accounting for approximately 30% of the study region, showed an advanced EOS, which was mainly distributed in the brush and alpine vegetation, including the Menyuan, Tianzhu areas, and the northwestern edges of the Qilian Mountains.
Similarly, the LGS of most pixels was extended during the period of 1982-2014 ( Figure 4c). In approximately 30% of the study region, the LGS showed a shortened trend from 1982 to 2014, which mainly included the area of the junction of Menyuan and Tianzhu, Delhi, and some individual pixels located in the Qilian Mountains.
Hence, it was crucial to analyze the phenology dynamics in the pixels with different phenology trends, considering the spatial heterogeneity of those trends. Approximately 66.93% of the study region had a delayed EOS trend. The EOS was significantly delayed from 1982 to 2014 in 58.74% of the area and was mainly distributed in the central area of the Qilian Mountains (Figure 4b). The remaining pixels, accounting for approximately 30% of the study region, showed an advanced EOS, which was mainly distributed in the brush and alpine vegetation, including the Menyuan, Tianzhu areas, and the northwestern edges of the Qilian Mountains.
Similarly, the LGS of most pixels was extended during the period of 1982-2014 (Figure 4c). In approximately 30% of the study region, the LGS showed a shortened trend from 1982 to 2014, which mainly included the area of the junction of Menyuan and Tianzhu, Delhi, and some individual pixels located in the Qilian Mountains.
Hence, it was crucial to analyze the phenology dynamics in the pixels with different phenology trends, considering the spatial heterogeneity of those trends. and 159.29 days. In summary, the SOS changed less than 1 day during the study period at the regional scale. In the Qilian Mountains, 64% of areas displayed an advanced trend of the SOS. The dynamics of the SOS in these areas were similar to those in the whole region: first fluctuating, advanced, and then fluctuating. From 1982 to 1989, the SOS fluctuated between 160.75 and 160.50 days. From 1989 to 2006, the SOS was significantly advanced by 1.7 days, much larger than the 0.65 days across the whole region. After 2006, the SOS fluctuated between 158.85 and 159.09 days, and from 2013 to 2014, the SOS was delayed. In summary, the SOS was advanced by approximately 2 days in these areas.
The SOS was delayed in 36% of the study region, and its dynamics could be divided into three phases: delayed, fluctuating, and then delayed. From 1982 to 1983, there was a significant delay in the SOS. From 1983 to 1996, the SOS fluctuated between 157.90 and 157.63 days. After 1996, the SOS tended to be delayed from 157.82 to 159.63 days, and the trend was slower from 2004 to 2014. In summary, the SOS was significantly delayed by 2.38 days in these areas during the study period.
In summary, there were similar SOS dynamics in different regions: first fluctuating, then significantly changing, and lastly fluctuating; however, the changed days and turning points were different. The SOS was significantly changed by approximately 2 days during the study period in the areas where the SOS showed an advanced trend and the areas with a delayed SOS, which was much larger than the 0.72 days of the whole region.  In the Qilian Mountains, 64% of areas displayed an advanced trend of the SOS. The dynamics of the SOS in these areas were similar to those in the whole region: first fluctuating, advanced, and then fluctuating. From 1982 to 1989, the SOS fluctuated between 160.75 and 160.50 days. From 1989 to 2006, the SOS was significantly advanced by 1.7 days, much larger than the 0.65 days across the whole region. After 2006, the SOS fluctuated between 158.85 and 159.09 days, and from 2013 to 2014, the SOS was delayed. In summary, the SOS was advanced by approximately 2 days in these areas.
The SOS was delayed in 36% of the study region, and its dynamics could be divided into three phases: delayed, fluctuating, and then delayed. From 1982 to 1983, there was a significant delay in the SOS. From 1983 to 1996, the SOS fluctuated between 157.90 and 157.63 days. After 1996, the SOS tended to be delayed from 157.82 to 159.63 days, and the trend was slower from 2004 to 2014. In summary, the SOS was significantly delayed by 2.38 days in these areas during the study period.
In summary, there were similar SOS dynamics in different regions: first fluctuating, then significantly changing, and lastly fluctuating; however, the changed days and turning points were different. The SOS was significantly changed by approximately 2 days during the study period in the areas where the SOS showed an advanced trend and the areas with a delayed SOS, which was much larger than the 0.72 days of the whole region. For the delayed and advanced regions, there was a turning point in approximately 2005. Before and after 2005, the SOS changed from a state of significant change to fluctuation. In the whole region of the Qilian Mountains, the corresponding turning point occurred in approximately 2000, which was different from 2005.

Dynamics of the EOS in Different Regions
The timelines in Figure 5b display the three dynamics of the EOS over the study period in different regions. In the whole region of the Qilian Mountains, the EOS had a delayed trend from 1982 to 2014, and its dynamic can be divided into two phases over time: In the Qilian Mountains, most areas had a delayed trend, which accounted for 67%. The dynamic of the EOS in these areas was similar to that in the whole region: fluctuating and then delayed. From 1982 to 2003, the EOS was significantly delayed by 2.04 days, which was much larger than the 1-day EOS across the whole region. After 2003, the EOS fluctuated drastically between 277.16 and 276.35 days. For the whole region and the delayed region, the same turning point occurred in 2003.
The EOS was advanced in 33% of the study region, and its dynamic could also be divided into two phases: advanced and then fluctuating. Before 2005, there was a significant advance in the EOS, and the trend was increasingly faster. From 1982 to 2005, the EOS was advanced by 2.08 days. From 2005 to 2014, the SOS fluctuated between 275.52 and 276.03 days.
In general, the EOS of different regions significantly changed during the first two decades of the study period, while there was no apparent trend in recent years. From 1982 to 2014, in the areas where the EOS showed an advanced trend and the areas with a delayed EOS, the EOS changed by 2 days, which was much larger than the 1-day EOS of the whole region. The weaker change in the whole region's EOS resulted from the opposing trends of different areas. For the delayed and whole regions, the same turning point occurred in 2003. For the advanced region of the Qilian Mountains, the turning point occurred in 2005.

Dynamics of the LGS in Different Regions
The dynamics of the LGS over the study period in different regions are displayed in Figure 5c In summary, the LGS was significantly extended from 1982 to 2003 and fluctuated from 2003 to 2014 in the whole region and the areas where the LGS showed an extended trend over the study period. In the remaining areas, the LGS was significantly shortened from 1982 to 2014. The trend during the period from 1997 to 2004 was faster than that of other periods. However, after 2004, the trend of the LGS was less significant.

The Relationship Between the SOS and Temperature
As shown in Table 2, the SOS had a significant negative correlation with T mean and T max and an insignificant negative correlation with T min at the regional scale of the Qilian Mountains. The correlation coefficient (r = −0.75) between the SOS and T max was slightly stronger than that (r = −0.60) between the SOS and T mean . For the relationship between the SOS and T mean , the SOS in approximately 64.28% of the study region was negatively correlated with T mean (Figure 6). Approximately 49.74% of the area had a significant correlation (p < 0.05) that was mainly distributed in the western and central Qilian Mountains. Other areas where the SOS and T mean were positively correlated accounted for 35.72% of the study region, which was distributed in the piedmont region, mainly including the junction area of Menyuan County and Tianzhu County, the northwestern part of the study region, and some areas around Qinghai Lake.
Occupying approximately 65.14% of the study region, the SOS showed a negative correlation with T max , and for other areas, a positive correlation was found ( Figure 6). The area with a very significant negative correlation (p < 0.01) between the SOS and T mean accounted for 49.86% of the study region, which was larger than 42.76% between the SOS and T mean . Similarly, the SOS and T max had a significant positive relationship in 20.38% of the study region. The area with a very significant negative correlation between the SOS and T mean only accounted for 16.37% of the Qilian Mountains. The spatial characteristics of the correlation coefficient between the SOS and T max were consistent with those between the SOS and T mean .
For the relationship between the SOS and T min , only approximately 18.48% of the study region was significantly negatively correlated with T min (p < 0.05), which was mainly distributed in Sunanyu, northern Gangcha, and Qilian. The area with a significant positive correlation (p < 0.05) between the SOS and T min accounted for 12.02% of the study region. Therefore, the SOS had a significant correlation with T min in 30% of the Qilian Mountains, which was obviously smaller than 71.38% between the SOS and T mean and 77.62% between the SOS and T max . Because of the weak correlation, T min was difficult to use to predict phenology in the phenology models. At the regional and pixel scales, different responses of the SOS to T mean , T max , and T min were found in this study. The correlation between the SOS and T max was significantly stronger than that of T min in the Qilian Mountains. T mean , eliminating extreme daily temperature, was also significantly correlated with the SOS. The spatial characteristics of the correlation coefficient between the SOS and T max and T mean were similar, while the general correlation between the SOS and T mean was slightly weaker than that with T max . Hence, T max was the most suitable temperature indicator for phenology models in the Qilian Mountains. For the relationship between the SOS and Tmin, only approximately 18.48% of the study region was significantly negatively correlated with Tmin (p < 0.05), which was mainly distributed in Sunanyu, northern Gangcha, and Qilian. The area with a significant positive correlation (p < 0.05) between the SOS and Tmin accounted for 12.02% of the study region. Therefore, the SOS had a significant correlation with Tmin in 30% of the Qilian Mountains, which was obviously smaller than 71.38% between the SOS and Tmean and 77.62% between the SOS and Tmax. Because of the weak correlation, Tmin was difficult to use to predict phenology in the phenology models. At the regional and pixel scales, different responses of the SOS to Tmean, Tmax, and Tmin were found in this study. The correlation between the SOS and Tmax was significantly stronger than that of Tmin in the Qilian Mountains. Tmean, eliminating extreme daily temperature, was also significantly correlated with the SOS. The spatial characteristics of the

The Relationship Between the EOS and Temperature
As shown in Table 2, a significant positive correlation between the EOS and T max and T mean exists at the regional scale in the Qilian Mountains, which was stronger than that between the EOS and T min . Different responses of the EOS to T max and T min have been found in the study region. The correlation between the EOS and T max was slightly stronger than that between the EOS and T mean (r = 0.53, p < 0.01).
In 68% of the Qilian Mountains, the EOS was positively related to T mean , which was mainly distributed in the central area of the study region ( Figure 6). The region with significant correlations (p < 0.05) accounted for 47% of the area. In 16.43% of the study region, the EOS was significantly negatively correlated with T mean , which mainly included Menyuan and Datong Counties and the northwestern edge of the Qilian Mountains, where the vegetation in these areas is mainly brush and alpine vegetation. Note that in the cultivated vegetation of Menyuan and Datong Counties, the EOS was advanced with increasing T mean , probably because human activities have changed the harvest time of local crops.
The correlation between the EOS and T max was slightly stronger than that between the EOS and T mean in the Qilian Mountains, and the spatial distribution of the correlation coefficient between them was similar. The EOS was significantly positive (p < 0.01), with T max in 47.68% of the study area, which was larger than the 40.93% between the EOS and T mean . Furthermore, these expanded areas were mainly distributed in southern Tianjun County. Accounting for approximately 29.36% of the study area, the EOS was negatively correlated with T max , which was mainly distributed in Menyuan County and the northwestern edge of the Qilian Mountains. The area with a significant negative correlation (p < 0.05) was approximately 15.06% of the study region.
The correlation between the EOS and T min was less significant than that between the EOS and T max . In 63.54% of the study region, the EOS was positively correlated with T min . Moreover, the area with a very significant correlation (p < 0.01) only accounted for 27.76% and was mainly distributed in northern Gangcha, Wulan, Qilian, and the southeastern part of Delhi. The area where the EOS had a negative correlation with T min expanded to Huangzhong and Tianzhu Counties. Accounting for 36.98% of the study region, the correlation between the EOS and T min was insignificant (p > 0.1) and was larger than the 25.53% between the EOS and T max . At the regional scale of the Qilian Mountains, the EOS was delayed by 0.1 day for a 1 • C increase in T min , and for a 1 • C increase in T max , the EOS was delayed by 0.18 days.
In summary, different EOS responses to T max and T min have been found in the Qilian Mountains. The response of the EOS to T min was weaker than that between the EOS and T max . T mean was also significantly correlated with the EOS. The spatial characteristics of the correlation coefficient between the EOS and the three temperature indicators were similar in the Qilian Mountains. However, the general correlation between the EOS and T max was the most significant. Hence, T max was a suitable temperature indicator for phenology and ecosystem models in the Qilian Mountains.

The Relationship Between the LGS and Annual Mean Temperature
At the regional scale of the Qilian Mountains, there was a significant positive correlation (r = 0.82, p < 0.01) between the LGS and annual mean temperature (Figure 7a). The LGS was extended by 0.56 days for a 1°C increase in mean temperature according to the regression analysis of the LGS and temperature. In particular, the LGS was positively correlated with the annual mean temperature in 74.41% of the study region. The positive correlation was significant (p < 0.05) in 63.48% of the areas. The areas with a negative correlation between the LGS and annual mean temperature were distributed in Menyuan and Tian Zhu Counties and some areas in Delhi, which accounted for 25.59% of the Qilian Mountains. The vegetation types in these areas are mainly meadow and brush.  In addition to the significant relationship, the interannual variation in the LGS was similar to that of annual mean temperature, as shown in Figure 7b. The annual mean temperature was lower from 1982 to 1996, and the LGS was shorter than the long-term trend. Similarly, the higher temperatures from 1996 to 2014 corresponded to the extended LGS In addition to the significant relationship, the interannual variation in the LGS was similar to that of annual mean temperature, as shown in Figure 7b. The annual mean temperature was lower from 1982 to 1996, and the LGS was shorter than the long-term trend. Similarly, the higher temperatures from 1996 to 2014 corresponded to the extended LGS in the Qilian Mountains. The phenology of the Qilian Mountains displayed significant spatial heterogeneity with increasing altitude (Figure 8). The multiyear mean phenology of pixels and their corresponding altitudes were extracted to explore the relationship between them. Generally, from low altitude to high altitude, the multiyear mean SOS was gradually delayed, the multiyear mean LGS was gradually shortened, and the multiyear mean EOS was first delayed and then advanced with the overall trend of advancement, which represented a significant altitude dependence. Figure 7. The spatial distribution of correlation coefficients and the interannual variation between the LGS and annual mean temperature in the Qilian Mountains. Correlation coefficient values of ±0.442, ±0.344, and ±0.291 are at the significance levels of 0.01, 0.05, and 0.1, respectively. (a) The spatial distribution of correlation coefficients between the LGS and annual mean temperature; (b) interannual variation in the LGS and annual mean temperature.
In addition to the significant relationship, the interannual variation in the LGS was similar to that of annual mean temperature, as shown in Figure 7b. The annual mean temperature was lower from 1982 to 1996, and the LGS was shorter than the long-term trend. Similarly, the higher temperatures from 1996 to 2014 corresponded to the extended LGS in the Qilian Mountains.

The Altitude Dependence of Phenology
The phenology of the Qilian Mountains displayed significant spatial heterogeneity with increasing altitude (Figure 8). The multiyear mean phenology of pixels and their corresponding altitudes were extracted to explore the relationship between them. Generally, from low altitude to high altitude, the multiyear mean SOS was gradually delayed, the multiyear mean LGS was gradually shortened, and the multiyear mean EOS was first delayed and then advanced with the overall trend of advancement, which represented a significant altitude dependence. There was a significant positive correlation between the multiyear mean SOS and altitude in the Qilian Mountains (Figure 8a). From 2500 to 4500 m, the SOS was delayed by approximately 30 days. The SOS value of pixels at the same altitude may also differ by There was a significant positive correlation between the multiyear mean SOS and altitude in the Qilian Mountains (Figure 8a). From 2500 to 4500 m, the SOS was delayed by approximately 30 days. The SOS value of pixels at the same altitude may also differ by approximately ten days, which may be due to local human activities and local microclimate influences.
The EOS had two different trends, with 3400 m as a turning point before and after it (Figure 8b). From 2500 to 3500 m, the EOS was gradually delayed with increasing altitude. Above 3400 m, the EOS showed a significant advanced trend from 280 to 265 days. The altitude of 3400 m was located in northern Qinghai Lake, with a significantly later EOS than the surrounding areas with higher altitudes and lower altitudes, which may be due to the abundant water resources enabling the growth of vegetation.
A negative correlation between the LGS and altitude was found in the Qilian Mountains (Figure 8c). From 2500 to 4500 m, the LGS was shortened by approximately 30 days. The shortened trend below 3200 m was not as significant as that above 3200 m, which was roughly consistent with the turning point of 3400 m in the EOS.
In summary, the SOS, EOS, and LGS of the Qilian Mountains were significantly related to altitude. The altitude dependence of phenology was mainly caused by temperature at different altitudes. Temperature is the key factor affecting phenological changes, and the relationship between the phenology and temperature of the Qilian Mountains was also explored in this study.

Phenology of Different Vegetation Types
We further analyzed the phenology of different vegetation types. In the Qilian Mountains, meadow, grassland, shrub, and alpine vegetation are the main vegetation types, occupying more than 80% of the study area (Table A1). Therefore, the phenological dif-ferences of these main vegetation types and their responses to increasing temperature were analyzed.
As shown in Table A1, the phenological differences of different vegetation types are quite significant. The differentiation of phenology in different vegetation types is consistent with the altitude distribution of phenology. The brush has the earliest SOS, the longest LGS, and the latest EOS among the main vegetation types and is distributed in the southeastern Qilian Mountains with a low elevation. Alpine vegetation has the latest SOS, the earliest EOS, and the shortest LGS, which is distributed in high-altitude areas.
A different response characteristics of different vegetation types to temperature can be seen in Table A2, especially the EOS. The EOS has a significant positive correlation with temperature in brush and meadow, while the EOS was significantly advanced in the alpine vegetation with the increasing temperature. During the period of 1982 to 2014, the EOS of alpine vegetation has been advanced about one day, partly due to the intersections of temperature and other environmental factors. Besides, note that the plant pieces in the same vegetation type of the same area may still have different response characteristics.
The AVHRR NDVI dataset used for phenology data extraction in this study has a low spatial resolution of 1/12 • , which may hide some spatial details of phenology parameters. The high-resolution NDVI dataset could better analyze the differentiation of different vegetation types, while the available years of these datasets are generally shorter than the AVHRR NDVI dataset.

The Turning Point of SOS Dynamics During the Study Period
This study showed an advanced SOS trend from 1982 to 2014 at the regional scale in the Qilian Mountains. However, the SOS fluctuated and showed no significant trend after 2000, and the advanced SOS mainly occurred before 2000. Other researchers have reported similar SOS dynamics in phenology on the Tibetan Plateau. Using GIMMS NDVI data of the TP, Piao et al. [36] found that the SOS had a significantly advanced trend in 1982-1999, while this advanced trend reversed from 1999 to 2006. Our results showed a similar SOS dynamic with a turning point in 2000, one year later than 1999 on the TP, which might be due to the smaller study region (Qilian Mountains in this study). In addition, Shen et al. [8] found no significant trend at the regional scale of the TP from 2000 to 2011 by combining NDVI data and multiple methods, which was consistent with the fluctuation in the SOS from 2000 to 2014 in the Qilian Mountains. Similarly, the SOS showed an advanced trend before 1998 in China and no evident trend thereafter.
The turning point in 2000 in the Qilian Mountains was consistent with that of larger regions such as the TP and China, which may be due to the similar insignificant climate change in recent years. The warming rate during the 1998-2012 period was significantly lower than the rate calculated since 1951 (IPCC, 2013), and this trend has been displayed in most areas of the globe. Phenology was significantly affected by temperature (found in this study), so the SOS trends of recent years were similarly smaller and insignificant. Note that in the areas where the SOS was advanced or was delayed in the Qilian Mountains, the turning points occurred in 2005, as opposed to 2000 in the whole region. In the research of large area, it is critical to explore the phenology dynamics in different areas.

Significantly Delayed Trend of the EOS in the Qilian Mountains
At the regional scale in the Qilian Mountains, the EOS mostly occurred in late September and October and was significantly delayed by 0.15 days per decade over the study period (p < 0.01). Compared with the TP, the delayed trend of EOS in the Qilian Mountains was more significant. Recent studies showed a delayed EOS trend on the TP, but the trend was insignificant. Che et al. [37], Ding et al. [26], and Cong et al. [9] reported that the EOS of the TP was delayed by 0.41 days per decade (p = 0.084), 0.15 days per decade (p = 0.75), and 0.7 days per decade (p = 0.184), respectively.
Because of the unique geographical conditions and the significant phenological change, the Qilian Mountains have become a typical research area for studying the phenology on the TP. The climate change in the Qilian Mountains is intense and the hydrological conditions are relatively uniform. On the contrary, different responses of phenology to climate change have been found on the TP, which result in the insignificant regional phenological change on the TP.

Different Responses of the SOS and EOS to T max and T min in the Qilian Mountains
This study revealed that the correlation between the SOS and T max was more significant than that between the SOS and T min in the Qilian Mountains. Plants need a certain amount of accumulative temperature above the specific threshold to trigger the SOS and the increasing temperature accelerated soil thawing and vegetation awakening so that the SOS advanced with increasing temperature [38]. Furthermore, T max played a more critical role in triggering the SOS. In the Qilian Mountains, T max was generally higher than 5 • C from 1982 to 2014, while T min was always below 0 • C, and its maximum value was −2.7 • C from 1982 to 2014. Because the temperature threshold mentioned before is generally 5 or 0 • C, the increase in T min cannot effectively prompt the accumulative temperature to advance the SOS from 1982 to 2014.
Some researchers have also verified the different SOS responses to T max and T min . However, the SOS is triggered by whichever one of the T max and T min needs to be explored, considering the local environmental conditions. Using the phenology data of in situ observations and satellite-derived data, Piao et al. [4] found that plant leaf onset was mainly driven by T max in the Northern Hemisphere (>30 • ), because in many areas, T min remained below 0 • C, which made it more difficult to fulfill the accumulative temperature than T max . In the warmer regions, the SOS was more driven by T min , because the increase in T max would increase evaporation, resulting in limited available water sources.
We also found that the EOS was more affected by T max in the Qilian Mountains. This result is similar to previous research [26]. Jeong et al. [39] found that plant senescence was triggered by cumulative cold temperatures below a specific threshold. Moreover, according to the Qinghai Meteorological Bureau, plant senescence is related to the pasture growth of grass and the end day of the local average temperature greater than or equal to 5°C. Hence, rising T max and T min can increase the average daily temperature to cause the EOS to be delayed in the Qilian Mountains, where most areas are located in Qinghai Province. An underlying explanation for the EOS being more affected by T max is that the annual variation in T max was greater than that in T min in the Qilian Mountains. This difference results in the overall daily mean temperature trend being more similar to T max .
Higher temperatures can alleviate frost damage by reducing freezing days and thus delaying vegetation senescence. In some arid areas, higher temperatures lead to intense evaporation, reducing the available water resources for vegetation and thus damaging ordinary growth. Yang et al. [40] showed that T min played a significant role in delaying the EOS by slowing chlorophyll degradation and alleviating frost injury, while increasing T max has little effect on delaying the EOS on the TP, which is due to the interacting effects of increasing T max and limited water resource caused by the high evaporation and less precipitation.
According to results of this study, the SOS and EOS were more related to T max in the Qilian Mountains than T min and T mean . However, T mean is widely used to project phenology in many phenology and ecosystem models. Therefore, these models based on T mean may not be able to extract phenology data effectively, and we should pay more attention to the important role of T max in phenology models.

Future Research
In addition to the temperature, the interaction of temperature and precipitation is also the influencing factor of phenological change [8,10,41]. For example, we found that in most parts of the Qilian Mountains, the SOS was advanced, and the EOS was delayed with the increasing temperature, while a significantly different response characteristic existed in the northeastern area of the Qinghai Lake, where the delayed SOS and advanced EOS can be seen. This phenomenon is probably due to the limited water resources, caused by the reduced precipitation and more intense evaporation under the high temperature in this area from 1961 to 2012 [42]. Besides, the influence mechanism of precipitation in different months on phenology is different [10,12]. It is, thus, our future work to reveal the influence mechanism of precipitation and the interaction of temperature and precipitation on phenology.

Conclusions
This study used the satellite-derived method to extract the phenology of the Qilian Mountains and further analyzed their spatial variation and temporal dynamics. Additionally, the response of phenology to temperature has been explored. First, we extracted the SOS, EOS, and LGS using the NDVI dataset. Then, we used correlation analysis to explore the different responses of the SOS and EOS to T mean , T max , and T mean . Furthermore, the relationship between the LGS and annual mean temperature was also researched in this study.
(1) The spatial distribution of plant phenology in the Qilian Mountains showed significant spatial heterogeneity and altitude dependence. From the northwestern high mountains to the lower elevation of the southeastern region, the multiyear mean SOS gradually advanced, the multiyear mean LGS gradually extended, and the multiyear mean EOS was delayed as a whole. In particular, the multiyear mean EOS was delayed in northern Qinghai Lake and was advanced around eastern Huangzhong County, an agricultural planting area. (2) At the regional scale in the Qilian Mountains, significant trends (p < 0.01) were detected in three phenology parameters over the study period: an earlier SOS (0. The phenology of the Qilian Mountains was significantly correlated with temperature over the study period at the regional scale. With the increase in the annual mean temperature, the LGS was significantly extended. Moreover, a higher preseason temperature led to an earlier SOS and a later EOS in most study areas. A more significant correlation between the SOS, EOS, and T max than between T min and T mean was found at the regional scale. Therefore, our results suggest that T max must be considered in the current phenology and ecosystem models.  Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
Publicly available datasets were analyzed in this study. The data supporting the results of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Publicly available datasets were analyzed in this study. The data supporting the results of this study are available from the corresponding author upon reasonable request.