Climate Change and Livestock Management Drove Extensive Vegetation Recovery in the Qinghai-Tibet Plateau

: The vegetation of the Qinghai-Tibet Plateau (QTP), China, is diverse and sensitive to climate change. Because of extensive grassland degradation in the QTP, several ecological restoration projects, which affect the livestock population, have been implemented in the QTP. Although many studies have reported the impacts of climate change on vegetation in the QTP, our knowledge on the impacts of both climate change and livestock on vegetation remains very limited. Here, we investigated the impacts of climate change and livestock population on vegetation growth by using the annual maximum normalized difference vegetation index (NDVI max ) and growing-season climate data from 1981 to 2019. We analyzed the relationship between NDVI max and climate and livestock population using the modiﬁed Mann-Kendall trend Test and Pearson correlation analysis. For the entire QTP, NDVI max had a two-phase trend, with a slow rise during 1981–2000 and a rapid rise during 2000–2019. Overall, NDVI max in the QTP increased and decreased in 63.7% and 6.7% of the area in 2000–2019. In areas with signiﬁcant changes in NDVI max , it was strongly correlated with relative humidity and vapor pressure. The small positive trend in NDVI max during 1981–2000 was inﬂuenced by warmer and wetter climate, and the overgrazing by a large population of livestock slowed down the rate of increase in NDVI max . Livestock population for Qinghai and Tibet in recent years has been lower than in the 1980s.The warmer and wetter climate and substantial drops in the livestock population contributed to large recovery in vegetation during 2001–2019. Vegetation degradation in Qinghai during 1981–2000 and central-northern Tibet during 2000–2019 was driven mainly by drier and hotter climatic. Although 63.7% of the area in the QTP became greener, the vegetation degradation in central-northern Tibet should not be ignored and more measures should be taken to alleviate the impact of warming and drying climate. Our ﬁndings provide a better understanding of the factors that drove changes in vegetation in the QTP.


Introduction
The Qinghai-Tibet Plateau (QTP), China, has an average elevation of 4500 m, and is known as the roof and third pole of the world [1]. Its climate is characterized as arid, cold, windy, anoxic, sunny, and has large swings in daily air temperature [2]. The unique climate and geographical environment of the QTP have fostered diverse ecosystem types, including forests, shrubs, alpine grasslands, alpine meadows, and alpine desert. Studies have shown that air temperature in the QTP has increased twice the global average, and precipitation has significantly increased since 1982 [3,4]. The changes in climate on the QTP have altered the atmosphere and hydrological cycles and thus has reshaped the local environment [4]. As an indicator of climate change, vegetation plays a pivotal role in linking the biosphere, pedosphere, geosphere, hydrosphere, and atmosphere at the regional scale and across Asia [5,6]. Studying the impact of climate change on the vegetation of the QTP will help advance our understanding of changes in QTP vegetation under a warming climate and provide information that is vital to sustainable development [5,7].
NDVI is highly correlated with vegetation leaf area index, greenness, and productivity [8]. Spaceborne time series NDVI images are vital for monitoring the dynamics of vegetation at large spatial scales [9][10][11][12]. This index has been extensively used and applied at regional to global scales. Many scholars have used NDVI to study the response of vegetation to climate change in the QTP at different time scales and resolutions. Several studies have reported changes in vegetation in the QTP before 2013, and these studies have greatly contributed to monitoring vegetation in the QTP [3,5,6,[13][14][15]. However, these studies do not reflect the latest dynamics of vegetation change in the QTP. Chen [7] studied the characteristics of vegetation change during 2000-2019 and its correlation with climate factors, but they did not conduct in-depth study on the impact of the livestock population on vegetation changes.
The vastness of the QTP means that climate change and its impacts on vegetation growth vary from place to place. Therefore, it is necessary to study the difference in the spatial response of vegetation to climate change. Some previous studies used average NDVI (NDVI avg ) to represent the growth status of vegetation, such as mean annual NDVI, [5] and mean annual growing season NDVI [7], seasonal average NDVI [11], and some others used the annual maximum NDVI (NDVI max ) [16,17]. Both NDVI avg and NDVI max have been widely used. NDVI avg reflects the average state, while NDVI max is associated with peak greenness during the growing season [18]. Piao [19] reported that the grassland above-ground biomass was well correlated with the NDVI max .
Terrestrial ecosystems of the QTP are particularly sensitive to climate change during the growing season [20][21][22], and it is also an important time for biomass accumulation. Thus, climatic conditions during the growing season directly affect biomass accumulation and NDVI max . In previous research, air temperature and precipitation were the main climate variables studied for their effects on vegetation growth [6,13,15,23], while atmospheric water conditions have been rarely considered. According to Ding [24], atmospheric aridity affects grassland productivity in the QTP because plants close their stomata to prevent water loss when they are water stressed, which inhibits photosynthesis, carbon assimilation, and growth. Yuan [25] showed that increased atmospheric vapor pressure deficit reduced global vegetation growth.
In addition to climate change, the impact of livestock production (e.g., livestock population, types, grazing management, and policies) cannot be ignored. Overgrazing is widely considered to be a key factor affecting the structure, productivity, and nutritional value of natural grassland vegetation, as well as the main driver of grassland degradation [3,26,27]. The livestock population determines grazing intensity and whether it is overgrazed to some extent. Many researchers have used livestock population or the grazing pressure index to analyze the impact of anthropogenic activities on grassland productivity in the QTP and found that changes in livestock population affect vegetation productivity and grassland yield [3,28]. However, a better understanding of the impact of the latest livestock population on vegetation changes in the QTP is needed.
We found from our literature review that information on the spatial-temporal responses of vegetation to climate change during the growing season and livestock population in the QTP was very limited and outdated, which indicates that there is a knowledge gap. To address this gap, here, we aimed to analyze the spatiotemporal vegetation dynamics and its response to climate changes and livestock population in the QTP from 1981 to 2019. We used NDVI max as an index of vegetation greenness and we used relative humidity and vapor pressure as indicators of atmospheric water conditions to address the following questions: (1) What are the spatiotemporal changes of NDVI max in the QTP from 1981 to 2019? (2) How does NDVI max respond to climate change (temperature, precipitation, relative humidity, and vapor pressure) in the QTP during 1981-2019? (3) How does NDVI max respond to both climate change and livestock number changes in the QTP

Study Area
The QTP extends from the southern Himalayas in the south to the northern Qilian Mountains in the north, the Pamirs and Karakoram Mountains in the west, and the Hengduan Mountain and the Loess Plateau in the east [29], and covers all of China's Tibet, most of Qinghai, and some areas in Gansu, Xinjiang, and the Western Sichuan Plateau, with a total area of 2.58 million square kilometers ( Figure 1). The average elevation of Tibet is 4736 m, which is higher than Qinghai's average of 4034 m. The Tibet pastoral area and Qinghai pastoral area are two of the four largest pastoral areas in China. The central and eastern QTP is the source of many large rivers, such as the Yangtze River, Yellow River, Lancang River, Mekong River, Yarlung Zangbo River, and Indus River. The QTP has dense river networks, abundant water resources, and is known as Asia's Water Tower. The growing season for vegetation in the QTP is from May to September, and vegetation often reaches its maximum biomass in July or August.The precipitation in the QTP gradually decreases from southeast to northwest. Precipitation in summer (from June to August) accounts for as much as 40%-60% of the annual precipitation in most areas of the QTP [30].

NDVI Data from AVHRR and MODIS
Two data sources, the third generation GIMMS NDVI from the AVHRR sensors and the MOD13A2 from the Moderate Resolution Imaging Spectrometer (MODIS), were used to create time series NDVI from 1981-2019. According to Du [31] and Zhang [14], NDVI values of the two data sources were highly correlated in the QTP, and they both reflect vegetation dynamics. In addition, Zhao [32] found that the two data sources responded to all climatic variables and had the same trend. Therefore, it is credible to analyze the trend of changes in vegetation and its correlation with climatic variables using the time series NDVI images from the two data sources. We constructed time series NDVI from 1981 to 2000 using GIMMS 3 g NDVI with a spatial resolution of 8 km and constructed the NDVI

NDVI Data from AVHRR and MODIS
Two data sources, the third generation GIMMS NDVI from the AVHRR sensors and the MOD13A2 from the Moderate Resolution Imaging Spectrometer (MODIS), were used to create time series NDVI from 1981-2019. According to Du [31] and Zhang [14], NDVI values of the two data sources were highly correlated in the QTP, and they both reflect vegetation dynamics. In addition, Zhao [32] found that the two data sources responded to all climatic variables and had the same trend. Therefore, it is credible to analyze the trend of changes in vegetation and its correlation with climatic variables using the time series NDVI images from the two data sources. We constructed time series NDVI from 1981 to 2000 using GIMMS 3 g NDVI with a spatial resolution of 8 km and constructed the NDVI dataset from 2001 to 2019 using MOD13A2 with a resolution of 1 km. We resampled the MOD13A2 images to 8 km to match the GIMMS 3 g NDVI spatial resolution.
We used the maximum value compositing (MVC) method to calculate NDVI max because it can eliminate the influence of clouds on the vegetation index. Time series NDVI max Remote Sens. 2021, 13, 4808 4 of 18 images during 1981-2019 were created using the Google Earth Engine cloud computing platform, which includes widely used geospatial datasets, has massive computational capabilities, and provides numerous algorithms [33].

Meteorological Data
The monthly average air temperature, monthly precipitation, relative humidity, and monthly average vapor pressure were collected from 117 standard weather stations in the QTP and its surrounding area during 1981-2019, which were provided by the Meteorological Data Sharing Service System of China (http://data.cma.cn/, accessed on 17 October 2021). Those data were interpolated using the Thin Plate Spline (TPS) method to obtain the spatial distribution of climate variables over the QTP with a spatial resolution of 8 km.
The monthly average temperature, monthly average precipitation, monthly average relative humidity, and monthly average vapor pressure during the annual growing season was calculated using the interpolated meteorological data from May to September of each year.

Livestock Data
Tibet and Qinghai are the main part of the QTP. The livestock population in Tibet and Qinghai from 1980 to 2019 were provided by the National Bureau of Statistics of China (http://data.stats.gov.cn/easyquery.htm?cn=E0103, accessed on 17 October 2021), and included the number of large livestock (e.g., cattle, horses, camels) and small livestock (e.g., goats and sheep).
The GIMMS 3 g NDVI and MOD13A2 covered different time periods, so we divided the entire research period into two stages: 1981-2000 and 2001-2019. In the subsequent data processing and analysis of this paper, we divided the NDVI data, meteorological data, and livestock population according to these two periods.

Modified Mann-Kendall Test
The Modified Mann-Kendall (MK) test [34] was used to analyze the interannual trends of climate variables and NDVI max at the pixel level. In many situations, the observed time series data are autocorrelated and this autocorrelation will result in misinterpretation of trend tests results when using the MK trend test [35]. The modified MK trend Test method as described by Yue [34] was able to eliminate the influence of sequence correlation on the trend test, and it has been shown to be as powerful as the original MK trend test when applied to uncorrelated sequences [36].
When the rates of NDVI max and climate variables were significantly greater than 0 at a level of p < 0.05, they were considered a significant increase, and when the rates of NDVI max and climate variables were significantly less than 0 at a level of p < 0.05, they were considered a significant decrease. In other cases, the rate was insignificant.

Linear Regression
We used linear regression to calculate the trends in NDVI max , climate variables, and livestock population at the regional scale. The specific calculation formula can be found in the papers of Huang [6] and Chen [7].

Pearson Correlation Analysis
Pearson correlation analysis was used to calculate the relationships between NDVI max and monthly temperature, monthly precipitation, monthly relative humidity, and monthly vapor pressure. The Person correlation analysis was proposed by Karl Pearson, [37] and the Pearson correlation coefficient r is often used to measure the strength of correlation between two variables.

Spatial-Temporal Dynamics of NDVI max during 1981-2019
The spatial distribution of the average NDVI max from 1981-2000 and 2001-2019 was shown in Figure 2a,b, respectively. The trends in NDVI max and their significance during 1981-2000 and 2001-2019 were calculated using the modified MK Test method as described by Yue [34]. The results are shown in Figure 2c monthly vapor pressure. The Person correlation analysis was proposed by Karl Pearson, [37]and the Pearson correlation coefficient r is often used to measure the strength of correlation between two variables.

Spatial-Temporal Dynamics of NDVImax during 1981-2019
The spatial distribution of the average NDVImax from 1981-2000 and 2001-2019 was shown in Figure 2a,b, respectively. The trends in NDVImax and their significance during 1981-2000 and 2001-2019 were calculated using the modified MK Test method as described by Yue [34]. The results are shown in Figure 2c-f. The average NDVImax values had a negative trend from southeast to northwest (Figure 2a,b). The southeastern and middle eastern QTP had higher NDVI values, and the corresponding vegetation types there are forest, shrub, and meadow. In the west and northwest, especially in the Qaidam Basin and Kunlun Mountains, the NDVI values were The average NDVI max values had a negative trend from southeast to northwest (Figure 2a,b). The southeastern and middle eastern QTP had higher NDVI values, and the corresponding vegetation types there are forest, shrub, and meadow. In the west and northwest, especially in the Qaidam Basin and Kunlun Mountains, the NDVI values were relatively small, and the area is classified as alpine desert. In the subsequent analysis, we masked out the pixels with a perennial mean of NDVI max less than 0.1 to exclude sparse vegetation areas and the influence of water bodies and shadows.
The vegetation dynamics in the QTP during the two periods had different spatial characteristics. As shown in Figure 2c,e, the negative trends in NDVI max during the first period were in eastern QTP, central and southern Qinghai, the junction of Qinghai, Sichuan, and Tibet, and in the vicinity of Hengduan Mountain and northwestern Tibet. In 1981-2000, NDVI max increased mainly in central Tibet, at the junction of Qinghai, Tibet, and Xinjiang. In 2001-2019, NDVI max in most areas of the QTP had a positive trend, such as in central and northern Tibet, the western Sichuan plateau, and the southern and western parts of Tibet Figure 2d,f. North-central Qinghai had the highest positive trend in NDVI max , which was greater than 0.002 yr −1 .
Area with greenness increased significantly, from 36.3% in 1981-2000 to 63.7% in 2000-2019. The area with negative trends in NDVI max became smaller, from 11.3% in 1981-2000 to 6.7% in 2001-2019. During 2001-2019, there were still some areas with a negative trend in NDVI max , which were found in central-northern Tibet and southern Qinghai.
In eastern Qinghai and the Western Sichuan Plateau, degraded vegetation had a positive trend in NDVI max after 2000. Figures 1 and 2f showed that areas with decreased NDVI max during 2001-2019 were located at a higher altitude than in 1981-2000. According to Li [38], alpine grasslands at higher elevations in the QTP have been found to be more sensitive to climate variability than those at lower elevations.  Figure 3 shows that the four climatic variables have strong spatial differences. Among them, both precipitation and relative humidity gradually increased from northwest to southeast. We also found that the distribution of NDVI max was highly consistent with the distribution of precipitation and relative humidity in the QTP at regional scale ( Figure 3). Our finding agree with Fang [39] and Shi [40], in that precipitation and vegetation cover were positively correlated in the QTP.

The Relationships between NDVI max and Climate
Remote Sens. 2021, 13, x FOR PEER REVIEW 6 of 18 relatively small, and the area is classified as alpine desert. In the subsequent analysis, we masked out the pixels with a perennial mean of NDVImax less than 0.1 to exclude sparse vegetation areas and the influence of water bodies and shadows. The vegetation dynamics in the QTP during the two periods had different spatial characteristics. As shown in Figure 2c In eastern Qinghai and the Western Sichuan Plateau, degraded vegetation had a positive trend in NDVImax after 2000. Figures 1 and 2f showed that areas with decreased NDVImax during 2001-2019 were located at a higher altitude than in 1981-2000. According to Li [38], alpine grasslands at higher elevations in the QTP have been found to be more sensitive to climate variability than those at lower elevations.      Figure 3 shows that the four climatic variables have strong spatial differences. Among them, both precipitation and relative humidity gradually increased from northwest to southeast. We also found that the distribution of NDVImax was highly consistent with the distribution of precipitation and relative humidity in the QTP at regional scale ( Figure 3). Our finding agree with Fang [39] and Shi [40], in that precipitation and vegetation cover were positively correlated in the QTP.

The Relationships between NDVImax and Climate Changes
From 1981 to 2000, the average temperature in most areas of the QTP had a significant positive trend (Figures 4a and 5a). The temperature in the northern QTP rose rapidly, which ranged from 0.05 to 0.1 °C (Figure 4a). Both precipitation and relative humidity had obvious spatial differences during 1981-2000, with some areas showing an obvious positive trend, and some areas showing an obvious negative trend. Precipitation in Tibet had    Figure 3 shows that the four climatic variables have strong spatial differences. Among them, both precipitation and relative humidity gradually increased from northwest to southeast. We also found that the distribution of NDVImax was highly consistent with the distribution of precipitation and relative humidity in the QTP at regional scale ( Figure 3). Our finding agree with Fang [39] and Shi [40], in that precipitation and vegetation cover were positively correlated in the QTP.
From 1981 to 2000, the average temperature in most areas of the QTP had a significant positive trend (Figures 4a and 5a). The temperature in the northern QTP rose rapidly, which ranged from 0.05 to 0.1 °C (Figure 4a). Both precipitation and relative humidity had obvious spatial differences during 1981-2000, with some areas showing an obvious positive trend, and some areas showing an obvious negative trend. Precipitation in Tibet had From 1981 to 2000, the average temperature in most areas of the QTP had a significant positive trend (Figures 4a and 5a). The temperature in the northern QTP rose rapidly, which ranged from 0.05 to 0. Regionally, the temperature increased, and precipitation decreased in Qinghai during the first period, which indicated that there was a warming and drying trend that would not have been conducive to vegetation growth. Meanwhile, all the four climatic variables, temperature, precipitation, relative humidity, and vapor pressure, in Tibet increased, which signaled a trend in warming and humidification from 1981 to 2000. This trend of climate change is more beneficial to vegetation productivity, because a certain range of increased temperature and precipitation can increase photosynthesis and growth [6,41]. From 2000 to 2019, temperature, precipitation, and vapor pressure in most of Qinghai had a positive trend. Thus, it can be inferred that evaporation decreased because of the trend in relative humidity, and that the change in climate in Qinghai from 2000 to 2019 was beneficial to vegetation growth. During 2001-2019, the temperature in Tibet increased significantly, while the precipitation in most areas of Tibet did not significantly change. Meanwhile, considering that the relative humidity in this area decreased significantly, it can be inferred that the evaporation increased. The climate change in some areas of Tibet during the second stage may inhibit the productivity and growth of vegetation. Figures 2a,b and 3 show that the distribution of NDVI max was consistent with the water conditions of the QTP, especially precipitation and relative humidity. The southeastern QTP had a suitable temperature, sufficient precipitation, and high NDVI max , whereas the northwestern QTP had a relatively low temperature, less precipitation, and a low NDVI max .

Correlation between NDVI max and Climatic Variables
The spatial distributions of the correlation between NDVI max and climatic variables during the growing season were computed and are shown in Figure 6. The frequency of the correlations between NDVI max and the four climatic variables are shown in Figure 7.
As is shown in Figure 6, the correlation between NDVI max and the four climatic variables could vary among the different stages and places in the QTP. Overall, NDVI max and the four climatic variables mainly had a positive correlation (Figure 7). The correlation between NDVI max and climate variables in 2001-2019 was stronger than in 1981-2000 ( Figure 6), which may indicate that the vegetation change in 2001-2019 was more affected by climate change than in 1981-2000. The two data sources of NDVI may also lead to differences in the correlation between NDVI and climatic variables. From 2000 to 2019, the eastern QTP had a positive correlation between temperature and NDVI max , especially at the junction of Qinghai and Sichuan, where the correlation coefficient was greater than 0.8. But in most areas of Tibet, temperature and NDVI max were negatively correlated (Figure 6e). The reason for the negative correlation may be that when precipitation in those areas decreased (Figures 4b and 5b), increased temperature (Figures 4a and 5a) can exacerbate drought and inhibit vegetation activity. Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 18  As is shown in Figure 6, the correlation between NDVImax and the four climatic variables could vary among the different stages and places in the QTP. Overall, NDVImax and the four climatic variables mainly had a positive correlation (Figure 7). The correlation between NDVImax and climate variables in 2001-2019 was stronger than in 1981-2000 (Figure 6), which may indicate that the vegetation change in 2001-2019 was more affected by climate change than in 1981-2000. The two data sources of NDVI may also lead to differences in the correlation between NDVI and climatic variables. From 2000 to 2019, the eastern QTP had a positive correlation between temperature and NDVImax, especially at the junction of Qinghai and Sichuan, where the correlation coefficient was greater than 0.8. But in most areas of Tibet, temperature and NDVImax were negatively correlated ( Figure  6e). The reason for the negative correlation may be that when precipitation in those areas decreased (Figures 4b and 5b), increased temperature (Figures 4a and 5a) can exacerbate drought and inhibit vegetation activity.
To further analyze the relationship between NDVImax and climate in areas with significant changes in NDVImax, those areas where vegetation significantly increased or decreased were extracted, and annual average NDVImax of those areas and the corresponding monthly average temperature, monthly average precipitation, and monthly average  As is shown in Figure 6, the correlation between NDVImax and the four climatic variables could vary among the different stages and places in the QTP. Overall, NDVImax and the four climatic variables mainly had a positive correlation (Figure 7). The correlation between NDVImax and climate variables in 2001-2019 was stronger than in 1981-2000 (Figure 6), which may indicate that the vegetation change in 2001-2019 was more affected by climate change than in 1981-2000. The two data sources of NDVI may also lead to differences in the correlation between NDVI and climatic variables. From 2000 to 2019, the eastern QTP had a positive correlation between temperature and NDVImax, especially at the junction of Qinghai and Sichuan, where the correlation coefficient was greater than 0.8. But in most areas of Tibet, temperature and NDVImax were negatively correlated ( Figure  6e). The reason for the negative correlation may be that when precipitation in those areas decreased (Figures 4b and 5b), increased temperature (Figures 4a and 5a) can exacerbate drought and inhibit vegetation activity.
To further analyze the relationship between NDVImax and climate in areas with significant changes in NDVImax, those areas where vegetation significantly increased or decreased were extracted, and annual average NDVImax of those areas and the corresponding monthly average temperature, monthly average precipitation, and monthly average To further analyze the relationship between NDVI max and climate in areas with significant changes in NDVI max , those areas where vegetation significantly increased or decreased were extracted, and annual average NDVI max of those areas and the corresponding monthly average temperature, monthly average precipitation, and monthly average relative humidity, and monthly average vapor pressure were calculated and counted. The correlation between annual average NDVI max and the four climate variables in areas with significant changes in NDVI max was analyzed separately, and the results are shown in Figure 8.
As shown in Figure 8, the NDVI max in areas with increased NDVI max during 1981-2000 (I), NDVI max in areas with increased NDVI max area during 2001-2019 (III), and NDVI max in areas with decreased NDVI max during 2001-2019 (IV) show a strong correlation with the climatic variables. NDVI max in areas with decreased NDVI max during 1981-2000 (II) show a weak correlation with climatic variables. The correlation coefficients between the NDVI max and the temperature, precipitation, relative humidity, and vapor pressure in the areas with increased vegetation during 1981-2000 were 0.53 (p < 0.05), 0.38 (p > 0.05), 0.61, and 0.70 (p < 0.01) (Figure 8a-d). Thus, vegetation growth during the first period was significantly related to vapor pressure, followed by relative humidity and temperature. In areas with decreases in vegetation, the correlation coefficients between NDVI max and four climatic variables were small, −0.34, −0.19, −0.21, and −0.33, respectively. Although NDVI max showed a weak correlation with climatic variables in areas with decreases in vegetation during the first period, there were obvious spatial differences (Figures 2e and 6a-d). For example, in areas with vegetation degradation of northwestern Tibet, there was a strong positive correlation between NDVI max and precipitation. But in areas with vegetation degradation near Tanggula Mountains, NDVI max showed a weak correlation with the four climate variables.   (Figure 8a-d). Thus, vegetation growth during the first period was significantly related to vapor pressure, followed by relative humidity and temperature. In areas with decreases in vegetation, the correlation coefficients between NDVImax and four climatic variables were small, −0.34, −0.19, −0.21, and −0.33, respectively. Although NDVI max and temperature in the areas with increased NDVI max during the second period had the strongest correlation, with a correlation coefficient of 0.62 (p < 0.01) (Figure 8i), followed by precipitation, with a correlation coefficient of 0.57 (Figure 8j). In areas with decreased NDVI max during the second period, NDVI max was correlated with relative humidity and vapor pressure, with correlation coefficients of 0.58 and 0.60, respectively (Figure 8o,p).

The Relationships between NDVImax, Climate and Livestock Population at the Regional Scale
Due to the differences in livestock population, climate change, and vegetation, it is necessary to analyze the two provinces separately. In this section, we analyzed the impact of livestock population and climate change on NDVI max at the regional scale. Interannual variations of NDVI max were shown in Figure 9. According to Chen [3] and Fan [28], large livestock (such as cattle and horses) was converted to four sheep units. Then, we calculated the livestock population in Qinghai, Tibet, and total numbers, and then analyzed the trends in the livestock population. The results are shown in Figure 10a 8i), followed by precipitation, with a correlation coefficient of 0.57 (Figure 8j). In areas with decreased NDVImax during the second period, NDVImax was correlated with relative humidity and vapor pressure, with correlation coefficients of 0.58 and 0.60, respectively (Figure 8o,p).

The Relationships between NDVImax, Climate and Livestock Population at the Regional Scale
Due to the differences in livestock population, climate change, and vegetation, it is necessary to analyze the two provinces separately. In this section, we analyzed the impact of livestock population and climate change on NDVImax at the regional scale. Interannual variations of NDVImax were shown in Figure 9. According to Chen [3] and Fan [28], large livestock (such as cattle and horses) was converted to four sheep units. Then, we calculated the livestock population in Qinghai, Tibet, and total numbers, and then analyzed the trends in the livestock population. The results are shown in Figure 10a       To analyze the impact of livestock population on NDVI max , the interannual variations of NDVI max , trends in the livestock population, and interannual climate trends in Qinghai and Tibet were divided into three periods.
As shown in Figure 9a, there was a significant difference in the average NDVI max in the QTP during the two stages, which was caused by the different NDVI sources. NDVI max had a non-significant trend during the first stage, with an annual growth rate of 0.0004 yr −1 (p > 0.05), but NDVI max had a significant positive trend of about 0.0012 yr −1 (p < 0.01) during 2001-2019. Overall, the changes in livestock population in the QTP between 1981 and 2000 were not conducive to vegetation growth (Figure 10a). Although the livestock population decreased starting in 1991, the disturbed alpine vegetation can take many years to recover [42,43]. Therefore, the positive influence of warmer and wetter climate (Figure 11a-d) should be the main reason why the NDVI max did not decrease during 1981-2000 (Figure 10a). From 2001 to 2019, the livestock population in the QTP decreased at a rate of 0.07 million sheep units/year (p < 0.01) (Figure 11a). The decrease in the livestock population helped reduce the grazing pressure over the QTP. Simultaneously, the climate conditions became warmer and wetter (Figure 11a-d), and vegetation had a positive trend in NDVI max during this period (Figure 10a). In addition, NDVI max in the QTP peaked in 2018 (Figure 9a) with a peak in precipitation (Figure 11b), but other climatic factors did not significantly change. The simultaneous peak in NDVI max and precipitation indicated that under a constant temperature, higher precipitation can promote vegetation growth and be the driving force of changes in vegetation. Figure 9b, NDVI max in Qinghai changed little from 1981 to 1991, and it had a non-significant positive trend from 1992 to 2000 with a rate of 0.002 yr −1 (p > 0.05). From 2000 to 2019, NDVI max in Qinghai increased significantly at a rate of 0.002 yr −1 (p < 0.05). From 1981 to 1991, the livestock population in Qinghai increased significantly and peaked in 1991 (Figure 10b). During this period, temperature increased, and precipitation decreased (Figure 11e,f). Changes in weather conditions and livestock population are not beneficial to vegetation growth. However, there was no obvious change in NDVI max (Figure 9b), which may be due to the better elasticity of the vegetation ecosystem in this area, and this change was within its tolerance. From 1992 to 2000, the livestock population in Qinghai decreased significantly at a rate of 1.03 million sheep units/year (p < 0.01) (Figure 10b). During this period, the temperature and the precipitation increased (Figure 11e,f). Thus, the increase of NDVI max between 1992 and 2000 in Qinghai can be explained by the positive trend in temperature and precipitation and stark decreases in the livestock population.

As shown in
From 2000 to 2019, the livestock population in Qinghai maintained a downward trend (p < 0.05) (Figure 10a). During this period, both temperature and precipitation increased (Figure 11a It can be seen from Figure 10a-c that the livestock population for Qinghai and Tibet in recent years has been lower than in the 1980s. Overall, the control of the livestock population was beneficial to vegetation recovery.

Interannual Trends of NDVI in the Qinghai-Tibetan Plateau
As an indicator of vegetation greenness, NDVI has been widely used to assess vegetation dynamics [9-12] MOD13A2 and GIMMS-NDVI have been the most widely used source of NDVI data. In previous studies on vegetation dynamics in the QTP, nearly all used average NDVI (NDVI avg ) to represent the overall level of vegetation growth. For example, Zhang [5] used annual mean NDVI, and Chen [7] used averaged NDVI in the growing season.
Zhang [5] used annual NDVI avg derived from GIMMS and non-parametric rank-based MK test to analyze the interannual trends of vegetation change in the QTP and found that 26.4% of the pixels increased and 5.58% decreased significantly during 1981-2000. Using mean annual growing season NDVI from MOD13A2 and a univariate linear regression to determine the trend, Chen [7] found that NDVI avg in the growing season had a significant increase (0.0011 yr −1 ) during 2000-2019. They found that approximately 79.29% of the pixels had a greening trend and the greener area was mainly concentrated in northern QTP, especially in the northeast. Areas with a negative trend in greenness accounted for 20.71% of the area, which were mainly concentrated in the central plateau and north Tibet.
Here, we used NDVI max to represent the peak of vegetation greenness. Although our findings were basically consistent with previous studies, growth rate of NDVI and the areas where NDVI max changed were moderately different. Our results show that the annual NDVI max increased at a rate of 0.0012 yr −1 during 2000-2019, which is little higher than the growth rate of NDVI avg in the conclusion of Chen [7]. The difference in vegetation status represented by maximum NDVI and average NDVI may be the reason for the difference in the growth rate between them [16]. We found that annual NDVI max during 1981-2000 increased for 36.3% of the area in the QTP, and 11.3% of the area had a significant decrease in NDVI max (p < 0.05). Annual NDVI max during 2001-2019 showed that 63.7% of the vegetation became significantly greener and 6.7% of the QTP area had a significant decrease in NDVI max . The changed area of NDVI max during 1981-2000 in our results is larger than the results of Zhang [5]. The changed area of NDVI max during 2001-2019 is smaller than that of Chen [7]. The difference between our study and other studies is likely caused by different change detection methods (e.g., non-parametric rank-based MK test, modified Mann-Kendall test, and univariate linear regression). The linear regression method uses the least squares method to fit the linear model and determines the trend according to the slope of the fitted line [5]. The MK test is suitable for trend analysis of non-normally distributed hydrometeorological data and climatic time series [35]. Yue [34] proposed the modified MK test using effective sample size to eliminate the influence of serial correlation in time series. Therefore, the difference in these algorithms may lead to different results for the same climate data. In the future, perhaps we should compare the effectiveness of the different change detection methods in the QTP.

The Effects of Climate Change on the Interannual Trends of NDVI
Several studies evaluated the interannual trend of the climate variables [5,7]. Zhang [5] collected from meteorological stations and the non-parametric rank-based MK test to analyze trends of climate in QTP. They found that the temperature had a positive trend from 1981 to 2013, and that the annual rainfall in Tibet from 1981 to 2000 also had a positive trend, while the rainfall in Qinghai from 1981 to 2000 had a negative trend, which was consistent with our results. Chen [7] collected temperature data from the MOD11A2 product, calculated a moist coefficient to represent water availability, applied univariate linear regression to calculate the trends of climate change in the growing season, and calculated the correlation between vegetation and climatic factors in the QTP during 2000-2019. They found that the northeast, mid-east, and western edges of the plateau became cooler and wetter, while the southwest, mid-west, and southeast became warmer and drier. However, we found that from 2001 to 2019, nearly all of the QTP had increased temperature. The difference in these findings may be due to the use of different meteorological data and different trend analysis methods.
In terms of the correlation between NDVI and climate variables, Chen [7] showed that vegetation in northern QTP was mainly negatively correlated with temperature, while in the southern QTP there was a positive correlation. However, we found that NDVI in most areas of the QTP during 2000-2019, especially Qinghai and western Sichuan, was mainly positively correlated with temperature, while central-northern Tibet was negatively correlated with temperature. Our results were consistent with Zhou [44], who used GIMMS-NDVI data and temperature obtained from the China National Meteorological Information Center to investigate the relationship between alpine grassland vegetation and temperature from 1982 to 2006 over the QTP. They found that NDVI max was negatively correlated with temperature in areas with severe degradation, while NDVI max was mainly positively correlated with temperature in areas with greener vegetation. Both Zhou [44] and we used meteorological data from weather Station, while Chen [7] used MOD11A2 product. The different sources of meteorological data may be the reason why our results of the correlation analysis are similar with Zhou [44] but different with Chen [7].
In terms of the impact of climate on vegetation, our results agreed with Huang [6] who collected precipitation and temperature datasets from the Chinese Meteorological Data Service and used linear regression and the sequential MK method to detect trends. Their results showed that the climate in the arid regions across Tibet from 2000 to 2011 did not promote grassland growth and the effects of climate and human activities on vegetation varied with location and time period.
As shown in Figure 8, in areas where NDVI max changed significantly, relative humidity and vapor pressure had a strong correlation with NDVI max. Considering that relative humidity and vapor pressure both represent atmospheric moisture to some extent, our results indicated that atmospheric water content has a great influence on the vegetation growth in the QTP, which was consistent with the conclusions of Ding [24] and Yuan [45], whose research showed atmospheric water content affects vegetation growth.

The Effects of Livestock Management on the Interannual Trends of NDVI max
To strengthen the ecological protection of the QTP, the Chinese government adopted a series of measures in recent decades. Since 2003, the Chinese government has implemented the Grazing Withdrawal Program (GWP) in Tibet and other pastoral provinces to restore and improve vegetation conditions. Since then, many fencing projects have been implemented in the QTP to restore and protect the degraded grasslands. By 2012, the area of fenced land in northern QTP had reached 3.32 million hectares [46]. In these areas, grazing prohibition, no grazing, or grazing rotation policies were implemented. Since 2011, eight provinces and regions, including Tibet, Qinghai, Sichuan, Xinjiang, and Gansu, have implemented the Ecological Subsidy and Award System. These subsidies and rewards mainly include grazing prohibition subsidies, grass-livestock balance rewards, pasture seed subsidies, livestock seed subsidies, and production subsidies.
The size of the livestock population is the key to effective alpine grassland management in the QTP [3]. Many studies have found that these ecological protection measures have significantly reduced the livestock population, relieved the grazing pressure on grasslands, and contributed to the ecological restoration of pastures [3,46]. Here, we found that the livestock population in Tibet and Qinghai has decreased since 2004, which was after the implementation of the GWP. Li [47] conducted a short term field experiment (4-5 years) to quantify the impact of different grazing management regimes and found that fencing was the most suitable grazing management regime for QTP alpine meadows. Through meta-analyses and questionnaire-based surveys, Sun [48] reported that fences were effective in promoting aboveground vegetation growth for up to four years in degraded alpine meadows and for up to eight years in the alpine steppes of the QTP. Our results also show that the reduction in the livestock population had a positive impact on vegetation restoration. Since 2003, adverse effects of overgrazing on vegetation have been gradually curbed and the contribution of climate change to grassland degradation exceeds the impact of overgrazing in the QTP.
In terms of the impact of climate and livestock population on vegetation changes, our results were consistent with those of Lehnert [49]. Based on the precipitation data provided by the Tropical Rainfall Measurement Mission (TRMM, 3B42) from 2000 to 2013 and temperature from the ERA-Interim reanalysis data, Lehnert [49] used linear multitemporal correlation analysis methods to study the effects of climate change and livestock population on vegetation coverage. They found that the vegetation cover has increased since the year 2000 in north-eastern QTP due to an increase in precipitation, and that vegetation cover declined in central-western QTP due to rising temperature and declining precipitation. They concluded that climate variability, rather than overgrazing, was the primary cause for large-scale changes in vegetation in the QTP since the new millennium. However, vegetation change includes not only vegetation degradation, but also vegetation greenness. We found that both climate change and the decline in the livestock population were beneficial to vegetation greening. Future studies are needed to determine which factor contributed more to vegetation greening.

Conclusions
The annual growth rate of NDVI max during 2000-2019 had a significant positive trend with an annual growth rate of 0.0012 yr −1 (p < 0.01). From 1981 to 2000, 11.3% of the QTP area had a decrease in NDVI max , mainly in the central and eastern of QTP. From 2000 to 2019, 63.7% of the QTP area became greener and 6.7% of the QTP area had a decrease in NDVI max . These areas with a decrease in NDVI max were mainly located in central-northern Tibet and southern Qinghai.
We also found a strong correlation between NDVI max and vapor pressure in areas with significant changes in NDVI max , which indicated that atmospheric water content greatly impacted vegetation growth in the QTP. The reasons for NDVI max changing differed between the periods and for different places. Positive trends in average NDVI max in Tibet from 1981 to 2000 was mainly due to favorable climate change that vegetation found to be favorable. Vegetation degradation in large areas of Qinghai from 1981 to 2000 and central-northern Tibet during 2000-2019 was mainly affected by a drier and hotter climate, considering that the negative impact caused by the increase in livestock population was weakened or offset by later livestock management. Overall, significant reductions in the livestock population and favorable changes in climate during 2000-2019 contributed to vegetation restoration of the entire QTP.
Author Contributions: E.L. conceived of the idea, analyzed the data, wrote and revised the manuscript. X.X. contributed to the design of the research, commented and revised the manuscript. H.S. contributed to analysis of the results and revised the manuscript. X.Y. contributed to analysis of the results. Y.Z. and Y.Y. performed the computations. All authors have read and agreed to the published version of the manuscript.