Spatiotemporal Patterns of Vegetation Greenness Change and Associated Climatic and Anthropogenic Drivers on the Tibetan Plateau during 2000–2015

Alpine vegetation on the Tibetan Plateau (TP) is known to be sensitive to both climate change and anthropogenic disturbance. However, the magnitude and patterns of alpine vegetation dynamics and the driving mechanisms behind their variation on the TP remains under debate. In this study, we used updated MODIS Collection 6 Normalized Difference Vegetation Index (NDVI) from the Terra satellite combined with linear regression and the Break for Additive Season and Trend model to reanalyze the spatiotemporal patterns of vegetation change on the TP during 2000–2015. We then quantified the responses of vegetation variation to climatic and anthropogenic factors by coupling climatic and human footprint datasets. Results show that growing season NDVI (GNDVI) values increased significantly overall (0.0011 year−1, p < 0.01) during 2000–2015 and that 70.37% of vegetated area on the TP (23.47% significantly with p < 0.05) exhibited greening trends with the exception of the southwest TP. However, vegetation greenness experienced trend shifts from greening to browning in half of the ecosystem zones occurred around 2010, likely induced by spatially heterogeneous temporal trends of climate variables. The vegetation changes in the northeastern and southwestern TP were water limited, the mid-eastern TP exhibited strong temperature responses, and the south of TP was driven by a combination of temperature and solar radiation. Furthermore, we found that, to some extent, anthropogenic disturbances offset climate-driven vegetation greening and aggravated vegetation browning induced by water deficit. These findings suggest that the impact of anthropogenic activities on vegetation change might not overwhelm that of climate change at the region scale.


Introduction
Climate-induced changes in vegetation can feedback into the climate-vegetation-soil continuum via various biotic and abiotic interactions [1,2]. Additionally, anthropogenic activities can also drive Remote Sens. 2018, 10 vegetation change significantly at the global scale [3]. Therefore, the spatiotemporal variability in vegetation growth can be attributed to the heterogeneity of either climate change or anthropogenic influences, or of both [1,[3][4][5]. Thus, monitoring and attributing such spatiotemporal changes in vegetation growth is not only critical to understanding the role that this land cover type plays in the Earth's climatic system [6], but is also a requirement for developing more sustainable strategies and policies for ecosystem management [7,8].
Long-term records of vegetation indices derived from satellite remote sensing provide unparalleled information regarding the vegetation response to climatic and anthropogenic factors at regional-to-global scales over the last two decades, such as the Normalized Difference Vegetation Index (NDVI) derived from Generation of Global Inventory Modeling and Mapping Studies 3rd Version (GIMMS 3g ), Système Pour l'Observation de la Terre (SPOT), and MODerate resolution Imaging Spectroradiometer (MODIS) [9][10][11]. A range of evaluation studies have, however, noted the fact that these remote sensing products did not yield consistent patterns of vegetation change [12][13][14][15], likely due to the impacts of sensor shifts or degradation [14,16,17]. In some studies, the widespread areas of browning identified by MODIS Collection 5 (C5) vegetation indices from the Terra satellite could not be detected by Collection 6 (C6) vegetation indices, which is mainly because the algorithm used in the former did not take the effects of sensor degradation into account [14,16].
The Tibetan Plateau (TP), referred to as the "Earth's third pole" [18], is an ideal location for investigating how fragile and sensitive alpine ecosystems are in response to global change [19]. Numerous studies have addressed vegetation variation and its interactions with the atmosphere on the TP, generally using NDVI as a proxy of vegetation activity. Nevertheless, the magnitude and patterns of vegetation change in this plateau are still debated, to a large extent due to data uncertainties [20][21][22]. For example, growing season NDVI (GNDVI) derived from GIMMS 3g exhibited a decreasing trend from 2000 to 2012, however, the GNDVI from C5 products of MODIS Terra remained showing increasing trends for the same period on the TP [23][24][25]. Moreover, these inconsistencies of vegetation changing trends might further intensify the ongoing disputes over drivers for vegetation dynamics on the TP [24][25][26][27]. Studies based on GIMMS 3g NDVI, for example, have found a positive correlation between vegetation greenness and temperature [28,29], whereas studies based on MODIS Terra-C5 NDVI argued that climate warming adversely affected vegetation growth [24]. In addition, the TP has experienced pronounced warming since the 1950s [30], while a temporary warming slowdown since the late 1990s has been detected in this region [31], especially after the mid-2000s [32]. However, how recent warming fluctuation influenced vegetation growth on the TP remains poorly understood.
Furthermore, although vegetation change is a consequence of both climate change and anthropogenic activities, relatively less attention has been paid to the latter in studies of the TP [22,33]. The intensity of anthropogenic activities is low overall on the TP but increases rapidly [34,35]. The increase of human population has promoted the development of agriculture and animal husbandry, transportation, urbanization, and tourism, as well as the increase of other forms of human activities [36]. The human population, the number of livestock, and the accessibility factors were known as important driving factors in vegetation degradation at local scales [37][38][39]. Such studies were mostly conducted within an individual factor perspective, e.g., road network and settlement locations. The role of human activities is becoming more and more important in influencing vegetation change [27]. Several studies even argued that anthropogenic activities were the primary contributors to changes in one-third alpine vegetation on the TP since 2000 [33,40]. However, little is known about the relationship of vegetation change with comprehensive anthropogenic disturbances [22].
In light of MODIS NDVI datasets updated from C5 to C6 [16], it is reasonable to reanalyze the changes in vegetation greenness and its major underlying drivers across the TP since the start of the new millennium. Specifically, our aim is: (1) to clarify spatiotemporal patterns of changes in vegetation greenness on the TP; and (2) to quantify the attribution of vegetation change from the perspective of both climatic and anthropogenic factors.

Study Area
The TP encompasses an area of 2.56 × 10 6 km 2 at an average elevation of higher than 4000 m above sea level (a.s.l.) and is, thus, the highest and largest plateau on Earth [18]. This region is therefore characterized by low temperatures, especially at high altitudes. Annual precipitation in the southeast TP is more than 1000 mm, while it is less than 50 mm in the northwest [41], and this results in a corresponding humidity gradient from humid in the southeast to arid in the northwest. The interactions between these diverse climatic regimes and terrestrial ecosystems contribute to the development of diverse biomes, ranging from subtropical rainforest in the southeast to alpine desert in the northwest of the TP. Ecological regionalization divides the plateau into 11 distinct ecosystem zones ( Figure 1, Table 1) [41]. The intensity of anthropogenic activities exhibits significant spatial heterogeneity on the TP [34,35]. Specifically, the intensity of anthropogenic activities is mainly higher in the population-concentrated in the southern and eastern TP than in the sparsely-populated northwestern regions. The northwestern TP was identified to be one of the ten largest wilderness areas [42].

Study Area
The TP encompasses an area of 2.56 × 10 6 km 2 at an average elevation of higher than 4000 m above sea level (a.s.l.) and is, thus, the highest and largest plateau on Earth [18]. This region is therefore characterized by low temperatures, especially at high altitudes. Annual precipitation in the southeast TP is more than 1000 mm, while it is less than 50 mm in the northwest [41], and this results in a corresponding humidity gradient from humid in the southeast to arid in the northwest. The interactions between these diverse climatic regimes and terrestrial ecosystems contribute to the development of diverse biomes, ranging from subtropical rainforest in the southeast to alpine desert in the northwest of the TP. Ecological regionalization divides the plateau into 11 distinct ecosystem zones ( Figure 1, Table 1) [41]. The intensity of anthropogenic activities exhibits significant spatial heterogeneity on the TP [34,35]. Specifically, the intensity of anthropogenic activities is mainly higher in the population-concentrated in the southern and eastern TP than in the sparsely-populated northwestern regions. The northwestern TP was identified to be one of the ten largest wilderness areas [42].    NDVI is one of the most sensitive available indicators of large-scale vegetation growth and has therefore been utilized quite widely in investigations of vegetation greenness and underlying drivers [4,11]. C6 NDVI dataset from MODIS-Terra (MOD13A2) during 2000-2015 used in this study was acquired from National Aeronautics and Space Administration (NASA) (https://ladsweb.modaps.eosdis.nasa.gov/). The comprehensive MODIS products (e.g., NDVI) updated from C5 to C6 was initialized by NASA in February 2015 using a more advanced calibration approach, in order to address this sensor problem as well as other issues. Several evaluation studies have shown that the major sensor degradation impacts have been removed for MODIS-C6 NDVI as compared to MODIS-C5 NDVI [14,16]. This dataset comprises 16-day maximum composite with a spatial resolution of 1 km. We applied an adaptive Savitzky-Golay approach to smooth annual NDVI cycle and to reconstruct a time series in order to mitigate contamination caused by clouds, snow and ice cover; this was done using the TIMESAT software package in the MATLAB environment [43]. We then calculated growing season (May-September) averaged NDVI (GNDVI) in order to analyze change in vegetation greenness and its responses to climate change and anthropogenic disturbance. We masked out the pixels with a GNDVI less than 0.1 to exclude non-vegetated areas [23].

Climatic, Human Footprint Pressure, and Auxiliary Datasets
Climatic datasets, including mean monthly air temperature, precipitation, and solar radiation, were obtained from the China Meteorological Forcing Dataset (http://dam.itpcas.ac.cn/) [44]. The temperature dataset was the combination of observations from 740 meteorological stations across China and Princeton forcing data [45]. The precipitation dataset within south of 40 • N was generated by merging meteorological station records, GLDAS precipitation for the period between 1979 and 1998, and TRMM precipitation data for the period between 1998 and present. The solar radiation dataset was produced by combining the meteorological station data, GEWEX-SRB and GLDAS datasets. These climatic datasets have relatively high accuracy and a spatial resolution of 0.1 • × 0.1 • from the period between 1979 and the present [44], and have previously been applied to investigate the responses of vegetation growth to climate change [46,47]. In this study, we resampled all climatic datasets to a spatial resolution of 1 km using the nearest neighbor algorithm, to match with NDVI data.
We used a recently updated human footprint pressure dataset (https://www.nature.com/articles/ sdata201667) to quantify the relationship between vegetation change and anthropogenic activities [48]. This grid dataset provides a globally-standardized measure for cumulative human pressure on the terrestrial environment and has been used to investigate the impact of anthropogenic activities on environment change such as biodiversity loss [49] and vegetation change [50]. This dataset has a spatial resolution of 1 km covering both 1993 and 2009. In this study, we used human footprint data for 2009 because this is the year that falls within our study period. The values of global human footprint pressure range from 0 to 64 [48], and a high value corresponds with a high level of human pressure. In this study region, the maximum value of human footprint pressure is 46, and over 38% of the area has values of human footprint pressure below 1, while just 5% of the region has values of human footprint pressure above 12.
Auxiliary data, encompassing vegetation cover type, drought, topography, and livestock numbers, were also collected to support the results of this analysis and to enable a discussion of vegetation variation and major driving forces on the TP. The vegetation cover type data was obtained from Institute of Botany, Chinese Academy of Sciences (CAS) (Figure 1) [51]. The Standardized Precipitation-Evapotranspiration Index (SPEI) database was obtained from the Climatic Research Unit of the University of East Anglia (http://sac.csic.es/spei/) [52], which was based on monthly precipitation and potential evapotranspiration at a spatial resolution of 0.5 degrees for 1901-2015. A lower the SPEI value indicates the climate to be drier. The topography data were acquired from the United States Geological Survey global digital elevation model (DEM) with a 1 km spatial resolution grid (https://lta.cr.usgs.gov/GTOPO30). Livestock numbers in the study were obtained from the statistical yearbooks for Qinghai Province and Tibet Autonomous Region.

Linear Regression Analysis
We performed a linear regression analysis to analyze monotonic trends in GNDVI. Thus, the trend slope in a multi-year regression equation represents inter-annual change, and can be solved using the ordinary least squares method (OLS), as follows: where Slope refers to the inter-annual trend in GNDVI, n is the number of years simulated, and GNDVI t is the value of this index in the tth year. A positive slope therefore indicates vegetation greening, while a negative one indicates vegetation browning [14]. We determined the significance of variation via F-tests to calculate confidence levels. Additionally, the change percentages of GNDVI in each case correspond to the ratio between Slope and averaged GNDVI.

The Break for Additive Season and Trend Model (BFAST)
The BFAST model was developed to identify gradual and abrupt changes in land surface [53]. This approach is essentially an iterative procedure to decompose observed time series data into three additive components of trend, seasonality, and remainder components (e.g., noise), and has been widely applied to time series data derived from satellite images (e.g., NDVI) to detect climate-driven biophysical indicators [53,54]. The method is described mathematically, as follows: where Y t denotes the NDVI time series at time t, while T t , S t , and e t refer to the trend, seasonal, and remainder components of corresponding records, respectively. We used a model with either zero or one breakpoint (namely bfast01classify in the R project, available at https://cran.r-project.org/web/packages/bfast/index.html) as modified by de Jong et al. [54] to detect the most influential trend shifts in vegetation greenness. This model was always fitted to pre-processed monthly NDVI time series. However, as the plant growing season in most regions of the TP is very short (May-September) [55], we fitted pre-processed 16-day NDVI time series in this model. The existence of a breakpoint was estimated though moving sums (MOSUMs) of OLS residuals (OLS-MOSUM), and if OLS-MOSUM test signaled significant instability (at the 5% significance level), then the breakpoint was captured and its date was also extracted. Finally, the presence of no break or a single 'most important' event was confined to within the period 2004-2012 in order to avoid regression in one period with too few data points ( Figure S1).

Partial Correlation Analysis
Partial correlation analysis has been widely used to detect the major climatic driving factors on vegetation growth via exploring the links between vegetation growth and a single climate factor while eliminating the effects of other climatic factors [47,56]. Ongoing work has led to an increasing recognition of time lag effects on vegetation responses to climate change at the regional and global scale [5,56]. The time lag of the vegetation response to climate at the monthly scale is generally shorter than three months, but those time lags differ among different vegetation type [56]. Since the length of the growing season on the TP is generally less than five months, we considered time lags of up to two months in this study. In consideration of the time-lag effects of the vegetation greenness responses to different climatic factors, including temperature, precipitation and solar radiation, we selected growing season (May to September) without monthly lag, one-month lag of growing season (April to August), and two-month lag of growing season (March to July) as the three types of time lag for analyses. We then calculated partial correlations between GNDVI and each climatic factor, setting the other two variables for control in each ecosystem zone. For each climatic factor, the lag months (e.g., 0, 1, and 2) of the growing season exhibiting the maximum partial correlation coefficient were considered as the best periods to evaluate the responses of vegetation greenness to that climatic factor. We used F-tests to calculate the significance level for each partial correlation. Figure 2 shows inter-annual change of area-averaged GNDVI and its spatial distribution on the TP for the period of 2000-2015. The area-averaged GNDVI significantly increased at a rate of 0.0011 year −1 (p < 0.01, Figure 2a). GNDVI trend patterns from 2000 to 2015 were spatially heterogeneous, but the main GNDVI trend was increasing ( Figure 2b). GNDVI exhibited greening trends at 70.37% of vegetated area (23.47% significantly with p < 0.05) on the TP, primarily located in eastern and northern areas of the TP; about 16.56% of vegetated areas characterized by the greening trend over 0.003 year −1 , which mainly located in the northeast part of the TP. On the contrary, 29.63% of the vegetated area displayed browning trends, primarily in the south-central and southwest TP, but most of these trends were not significant; only 3.31% of vegetated areas characterized by the browning trend below −0.003 year −1 , which sparsely distributed within local areas of the southwest TP. However, such greening trends were not completely detected on the basis of GNDVI from Terra-C5 and GIMMS 3g , especially the latter. The details on differences within inter-annual variations of NDVI data derived from MODIS Terra-C5, Terra-C6, and GIMMS 3g can be found in Text S1 in the Supporting Information.

Inter-Annual Variations in Growing Season Vegetation Greenness
Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 16 setting the other two variables for control in each ecosystem zone. For each climatic factor, the lag months (e.g., 0, 1, and 2) of the growing season exhibiting the maximum partial correlation coefficient were considered as the best periods to evaluate the responses of vegetation greenness to that climatic factor. We used F-tests to calculate the significance level for each partial correlation. Figure 2 shows inter-annual change of area-averaged GNDVI and its spatial distribution on the TP for the period of 2000-2015. The area-averaged GNDVI significantly increased at a rate of 0.0011 year −1 (p < 0.01, Figure 2a). GNDVI trend patterns from 2000 to 2015 were spatially heterogeneous, but the main GNDVI trend was increasing (Figure 2b). GNDVI exhibited greening trends at 70.37% of vegetated area (23.47% significantly with p < 0.05) on the TP, primarily located in eastern and northern areas of the TP; about 16.56% of vegetated areas characterized by the greening trend over 0.003 year −1 , which mainly located in the northeast part of the TP. On the contrary, 29.63% of the vegetated area displayed browning trends, primarily in the south-central and southwest TP, but most of these trends were not significant; only 3.31% of vegetated areas characterized by the browning trend below −0.003 year −1 , which sparsely distributed within local areas of the southwest TP. However, such greening trends were not completely detected on the basis of GNDVI from Terra-C5 and GIMMS3g, especially the latter. The details on differences within inter-annual variations of NDVI data derived from MODIS Terra-C5, Terra-C6, and GIMMS3g can be found in Text S1 in the supporting information. We further calculated the change percentages of GNDVI for different elevation bins and ecological zones across the plateau over the study period. As shown in Figure 3, the trends of GNDVI were positive at all elevation bins (Figure 3a). The change percentages of GNDVI increased from about 0.40% year −1 at 2600 m to 0.49% year −1 at 3300 m and then decreases to 0.19% year −1 at 5500 m. The greening trends at the elevation bins below 4700 m a.s.l. were significant, while that above 4700 m a.s.l. were not significant. Moreover, the change percentages of GNDVI varied among the ecological zones (Figure 3b). The high value of greening percentages among the ecosystem zones occurred in the northern TP, including IID2, IID3, ID1, IC1, and IIC2 zones. On the contrary, the value of trends in IIC1 and IC2 zones were negative, but these browning trends were not significant. We further calculated the change percentages of GNDVI for different elevation bins and ecological zones across the plateau over the study period. As shown in Figure 3, the trends of GNDVI were positive at all elevation bins (Figure 3a). The change percentages of GNDVI increased from about 0.40% year −1 at 2600 m to 0.49% year −1 at 3300 m and then decreases to 0.19% year −1 at 5500 m. The greening trends at the elevation bins below 4700 m a.s.l. were significant, while that above 4700 m a.s.l. were not significant. Moreover, the change percentages of GNDVI varied among the ecological zones ( Figure 3b). The high value of greening percentages among the ecosystem zones occurred in the northern TP, including IID2, IID3, ID1, IC1, and IIC2 zones. On the contrary, the value of trends in IIC1 and IC2 zones were negative, but these browning trends were not significant.

Major Shifts in Vegetation Greenness Trends
Based on the BFAST model, we examined the major trend shifts in vegetation greenness across the entire TP and different ecosystem zones (Figure 4 and Figure S1). Monotonic greening trends in vegetation greenness were detected across the vegetated areas of entire TP. The trend shifts in vegetation greening and browning regimes, however, varied depending on ecosystem zones. All ecosystem zones have experienced abrupt changes with the exception of IIC2 and IIAB1 zones. Vegetation greenness in five ecosystem zones (i.e., IID2, IID3, ID1, IC1, and IB1 zones), located in the northern and middle regions of TP, shifted from greening to browning in either 2009 or 2010. Conversely, the vegetation greenness of IIC1 and IC2 zones in the southwestern part of TP was characterized by a shift regime of browning with a burst in 2011.  IID2  IID3  ID1  IC1  IIC2  IID1  OA1  IIAB1  IB1  IC2 TP  IID2  IID3  ID1  IC1  IIC2  IID1  OA1  IIAB1  IB1  IC2  IIC1 Year of trend shifts

Major Shifts in Vegetation Greenness Trends
Based on the BFAST model, we examined the major trend shifts in vegetation greenness across the entire TP and different ecosystem zones (Figure 4 and Figure S1). Monotonic greening trends in vegetation greenness were detected across the vegetated areas of entire TP. The trend shifts in vegetation greening and browning regimes, however, varied depending on ecosystem zones. All ecosystem zones have experienced abrupt changes with the exception of IIC2 and IIAB1 zones. Vegetation greenness in five ecosystem zones (i.e., IID2, IID3, ID1, IC1, and IB1 zones), located in the northern and middle regions of TP, shifted from greening to browning in either 2009 or 2010. Conversely, the vegetation greenness of IIC1 and IC2 zones in the southwestern part of TP was characterized by a shift regime of browning with a burst in 2011.
vegetation greening and browning regimes, however, varied depending on ecosystem zones. All ecosystem zones have experienced abrupt changes with the exception of IIC2 and IIAB1 zones. Vegetation greenness in five ecosystem zones (i.e., IID2, IID3, ID1, IC1, and IB1 zones), located in the northern and middle regions of TP, shifted from greening to browning in either 2009 or 2010. Conversely, the vegetation greenness of IIC1 and IC2 zones in the southwestern part of TP was characterized by a shift regime of browning with a burst in 2011.  TP  IID2  IID3  ID1  IC1  IIC2  IID1  OA1  IIAB1  IB1  IC2  IIC1 Year of trend shifts

Vegetation Greenness Responses to Climate Change
Based on partial correlation analyses, we demonstrated that the presence of clear disparities in vegetation greenness responses to air temperature, precipitation, and radiation in different ecological zones of the TP between 2000 and 2015. Time lag effects in vegetation responses to major climatic factors were also seen in some ecological zones (Figures 5-7).
Based on partial correlation analyses, we demonstrated that the presence of clear disparities in vegetation greenness responses to air temperature, precipitation, and radiation in different ecological zones of the TP between 2000 and 2015. Time lag effects in vegetation responses to major climatic factors were also seen in some ecological zones (Figures 5-7).
Increases in GNDVI were closely related to increasing precipitation on the northeastern TP, including within the IIC2 and IC1 zones (Figures 6 and 7), and a one-month lag effect response of GNDVI to precipitation existed in the IIC2 zone ( Figure 5). Although the positive partial correlation between GNDVI and temperature was significant within the IC1 zone, this was not the case within the IIC2 zone; and despite the fact that temperature declined in the IC1 zone before the GNDVI breakpoint, both GNDVI and precipitation increased (Figure 8a,b). We also demonstrated that both GNDVI and precipitation declined after the GNDVI breakpoint within IC1 zone. These results indicate vegetation change in the northeastern TP was sensitive to precipitation variation.
By contrast, GNDVI declined in the southwestern TP under the comprehensive influence of warming and increased solar radiation, combined with inter-annual fluctuations in precipitation (Figures 6 and 7). We also found a significant decrease of SPEI in this region during 2000-2015 ( Figure  S2), indicating that this region seems to be warming and drying. Besides, there was a two-month lag effect response of GNDVI to temperature, precipitation, and radiation within the IIC1 zone, and a one-month lag effect response of GNDVI to radiation within the IC2 zone ( Figure 5). Prior to the breakpoint of GNDVI, temperature and solar radiation increased, but precipitation decreased; all of these changes led to decreases in GNDVI, especially within the IIC1 zone (Figure 8c,d). Subsequent to the breakpoint of GNDVI, however, this index continued to decrease as warming and solar radiation increased despite a concomitant increase in precipitation. These findings might suggest that under conditions of insufficient water in the southwestern TP, increases in temperature and solar radiation both adversely affected vegetation growth as water availability decreased.   Figure 5). (a) The partial correlation between GNDVI and temperature (Tem) after controlling for total Based on partial correlation analyses, we demonstrated that the presence of clear disparities in vegetation greenness responses to air temperature, precipitation, and radiation in different ecological zones of the TP between 2000 and 2015. Time lag effects in vegetation responses to major climatic factors were also seen in some ecological zones (Figures 5-7).
Increases in GNDVI were closely related to increasing precipitation on the northeastern TP, including within the IIC2 and IC1 zones (Figures 6 and 7), and a one-month lag effect response of GNDVI to precipitation existed in the IIC2 zone ( Figure 5). Although the positive partial correlation between GNDVI and temperature was significant within the IC1 zone, this was not the case within the IIC2 zone; and despite the fact that temperature declined in the IC1 zone before the GNDVI breakpoint, both GNDVI and precipitation increased (Figure 8a,b). We also demonstrated that both GNDVI and precipitation declined after the GNDVI breakpoint within IC1 zone. These results indicate vegetation change in the northeastern TP was sensitive to precipitation variation.
By contrast, GNDVI declined in the southwestern TP under the comprehensive influence of warming and increased solar radiation, combined with inter-annual fluctuations in precipitation (Figures 6 and 7). We also found a significant decrease of SPEI in this region during 2000-2015 ( Figure  S2), indicating that this region seems to be warming and drying. Besides, there was a two-month lag effect response of GNDVI to temperature, precipitation, and radiation within the IIC1 zone, and a one-month lag effect response of GNDVI to radiation within the IC2 zone ( Figure 5). Prior to the breakpoint of GNDVI, temperature and solar radiation increased, but precipitation decreased; all of these changes led to decreases in GNDVI, especially within the IIC1 zone (Figure 8c,d). Subsequent to the breakpoint of GNDVI, however, this index continued to decrease as warming and solar radiation increased despite a concomitant increase in precipitation. These findings might suggest that under conditions of insufficient water in the southwestern TP, increases in temperature and solar radiation both adversely affected vegetation growth as water availability decreased.   Figure 5). (a) The partial correlation between GNDVI and temperature (Tem) after controlling for total Figure 6. The responses of the GNDVI to temperature (Tem), total precipitation (Pre), and radiation (Rad) within different ecosystem zones during 2000-2015 considering time-lag effects (according to Figure 5). (a) The partial correlation between GNDVI and temperature (Tem) after controlling for total precipitation (Pre) and radiation (Rad); and (b) the partial correlation between GNDVI and total precipitation (Pre) after controlling for temperature (Tem) and radiation (Rad); (c) the partial correlation between GNDVI and radiation (Rad) after controlling for temperature (Tem) and total precipitation (Pre). The significance level for partial correlation was calculated by F-tests.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 16 precipitation (Pre) and radiation (Rad); and (b) the partial correlation between GNDVI and total precipitation (Pre) after controlling for temperature (Tem) and radiation (Rad); (c) the partial correlation between GNDVI and radiation (Rad) after controlling for temperature (Tem) and total precipitation (Pre). The significance level for partial correlation was calculated by F-tests.  In the more humid areas, the impact of temperature and solar radiation became more prominent. GNDVI was mainly affected by climate warming across the IB1 zone, located in the mid-eastern TP (Figures 6 and 7). This increase of GNDVI was closely related to rising temperature prior to the breakpoint, and vice versa. In the southeastern region, including the IIAB1 and OA1 zones, GNDVI was mainly affected by a combination of temperature and solar radiation. Additionally, a one-month lag effect was found in the response of GNDVI to radiation in the IIAB2 zone ( Figure 5). In other Increases in GNDVI were closely related to increasing precipitation on the northeastern TP, including within the IIC2 and IC1 zones (Figures 6 and 7), and a one-month lag effect response of GNDVI to precipitation existed in the IIC2 zone ( Figure 5). Although the positive partial correlation between GNDVI and temperature was significant within the IC1 zone, this was not the case within the IIC2 zone; and despite the fact that temperature declined in the IC1 zone before the GNDVI breakpoint, both GNDVI and precipitation increased (Figure 8a,b). We also demonstrated that both GNDVI and precipitation declined after the GNDVI breakpoint within IC1 zone. These results indicate vegetation change in the northeastern TP was sensitive to precipitation variation.
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 16 precipitation (Pre) and radiation (Rad); and (b) the partial correlation between GNDVI and total precipitation (Pre) after controlling for temperature (Tem) and radiation (Rad); (c) the partial correlation between GNDVI and radiation (Rad) after controlling for temperature (Tem) and total precipitation (Pre). The significance level for partial correlation was calculated by F-tests.  In the more humid areas, the impact of temperature and solar radiation became more prominent. GNDVI was mainly affected by climate warming across the IB1 zone, located in the mid-eastern TP (Figures 6 and 7). This increase of GNDVI was closely related to rising temperature prior to the breakpoint, and vice versa. In the southeastern region, including the IIAB1 and OA1 zones, GNDVI was mainly affected by a combination of temperature and solar radiation. Additionally, a one-month lag effect was found in the response of GNDVI to radiation in the IIAB2 zone ( Figure 5). In other By contrast, GNDVI declined in the southwestern TP under the comprehensive influence of warming and increased solar radiation, combined with inter-annual fluctuations in precipitation ( Figures 6 and 7). We also found a significant decrease of SPEI in this region during 2000-2015 ( Figure S2), indicating that this region seems to be warming and drying. Besides, there was a two-month lag effect response of GNDVI to temperature, precipitation, and radiation within the IIC1 zone, and a one-month lag effect response of GNDVI to radiation within the IC2 zone ( Figure 5). Prior to the breakpoint of GNDVI, temperature and solar radiation increased, but precipitation decreased; all of these changes led to decreases in GNDVI, especially within the IIC1 zone (Figure 8c,d). Subsequent to the breakpoint of GNDVI, however, this index continued to decrease as warming and solar radiation increased despite a concomitant increase in precipitation. These findings might suggest that under conditions of insufficient water in the southwestern TP, increases in temperature and solar radiation both adversely affected vegetation growth as water availability decreased.
In the more humid areas, the impact of temperature and solar radiation became more prominent. GNDVI was mainly affected by climate warming across the IB1 zone, located in the mid-eastern TP (Figures 6 and 7). This increase of GNDVI was closely related to rising temperature prior to the breakpoint, and vice versa. In the southeastern region, including the IIAB1 and OA1 zones, GNDVI was mainly affected by a combination of temperature and solar radiation. Additionally, a one-month lag effect was found in the response of GNDVI to radiation in the IIAB2 zone ( Figure 5). In other ecosystem zones in the northwestern TP such as IID2 and IID3 zone, the partial correlations between GNDVI and temperature, precipitation, and solar radiation were insignificant.

Relationship between Vegetation Greenness and Anthropogenic Disturbance
Over the entire TP, an enhanced degree of vegetation greening is seen in areas where anthropogenic disturbance is lower (Figure 9). In order to evaluate this trend in more detail we, therefore, selected zones in two regions characterized by significant differences in both vegetation change, anthropogenic disturbance, and climate change-northeastern ecosystem zones (IIC2 and IC1 zones) and southwestern ecosystem zones (IC2 and IIC1 zones)-as the typical regions for further comparisons (Figures 2b, 7, and 8).
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 16 ecosystem zones in the northwestern TP such as IID2 and IID3 zone, the partial correlations between GNDVI and temperature, precipitation, and solar radiation were insignificant.

Relationship between Vegetation Greenness and Anthropogenic Disturbance
Over the entire TP, an enhanced degree of vegetation greening is seen in areas where anthropogenic disturbance is lower (Figure 9). In order to evaluate this trend in more detail we, therefore, selected zones in two regions characterized by significant differences in both vegetation change, anthropogenic disturbance, and climate change-northeastern ecosystem zones (IIC2 and IC1 zones) and southwestern ecosystem zones (IC2 and IIC1 zones)-as the typical regions for further comparisons (Figures 2b, 7, and 8). Along with increasing human footprint pressure, the percentage of greening trend declined gradually within the northeastern ecosystem zones (IIC2 and IC1 zones) while the percentage of browning trend increased in the southwestern ecosystem zones (IC2 and IIC1 zones) (Figure 9). This result implies an obvious negative impact of anthropogenic disturbance on vegetation greenness. Specifically, the trends in all levels of human footprint pressure in the northeastern ecosystem zones (IIC2 and IC1 zones) were positive, which indicates the fact that although the positive impact of wetting was weakened along with anthropogenic disturbance increase, the former was not overwhelmed by the latter. By contrast, in the southwestern ecosystem zones (IC2 and IIC1 zones), vegetation greenness in areas with lower human footprint pressure showed a greening trend, while a browning trend was clear in areas characterized by higher disturbance, suggesting that a larger degree of anthropogenic disturbance lead to a more negative impact of the water deficit on vegetation growth. Along with increasing human footprint pressure, the percentage of greening trend declined gradually within the northeastern ecosystem zones (IIC2 and IC1 zones) while the percentage of browning trend increased in the southwestern ecosystem zones (IC2 and IIC1 zones) (Figure 9). This result implies an obvious negative impact of anthropogenic disturbance on vegetation greenness. Specifically, the trends in all levels of human footprint pressure in the northeastern ecosystem zones (IIC2 and IC1 zones) were positive, which indicates the fact that although the positive impact of wetting was weakened along with anthropogenic disturbance increase, the former was not overwhelmed by the latter. By contrast, in the southwestern ecosystem zones (IC2 and IIC1 zones), vegetation greenness in areas with lower human footprint pressure showed a greening trend, while a browning trend was clear in areas characterized by higher disturbance, suggesting that a larger degree of anthropogenic disturbance lead to a more negative impact of the water deficit on vegetation growth.

Discussion
In this study, based on the updated NDVI form MODIS Terra-C6, we found that 70.37% of vegetated areas of the TP (23.47% significantly with p < 0.05) exhibited greening trends and the change percentages of GNDVI were positive across all elevation bins and most ecosystem zones over the TP during the past 16 years (Figures 2 and 3). These results further corroborate the finding that vegetation greenness on the TP has increased overall, albeit with the exception of local areas that have been degraded [57,58]. A similar pattern of greening has also been reported for the neighboring Himalayan region based on Terra-C6 NDVI [7]. However, GNDVI derived from Terra-C5 and GIMMS 3g underestimate the greening trend on the TP to different extents (the details can be found in Text S1 in the Supporting Information), largely due to the impacts of sensor shifts or degradation [14,16,17,20].
GNDVI increased in the eastern and northern TP but decreased in the south-central and southwest TP during 2000-2015 (Figure 2), and most ecosystem zones experienced major shifts in GNDVI between 2009 and 2011 ( Figure 4). The main explanation for this is that spatially heterogeneous temporal trends in major climatic factors induced complex responses in vegetation growth on the TP [59][60][61]. This study shows that the role of temperature, precipitation, and solar radiation played in vegetation dynamics varied depending on ecosystem zones, where the vegetation composition and humidity are different.
Vegetation growth in the southwestern and northeastern TP requires more water as these are mostly arid and semi-arid areas of alpine meadow and alpine steppe. These areas are consequently greatly affected by water availability. We show that vegetation greening was mainly the result of increasing precipitation on the northeastern TP, which is in accordance with the findings of other similar studies [24]. Geographically, the northeastern and southwestern TP are characterized by distinct water vapor sources and different thermal heating systems [62]. In contrast, within the context of the weakening Indian monsoon [62], vegetation browning in the southwestern TP was caused by a water deficit ( Figure S2), which induced by interactions among changes in temperature, solar radiation, and precipitation [60,63]. In this respect, the results of this study contradict the previous study that increasing temperature drove vegetation greening within the Yarlung Zangbo River Basin [26], perhaps because the greening trend was over-estimated on the basis of SPOT-VGT-NDVI ( Figure S3). Furthermore, a two-month lag effect response of vegetation to temperature, precipitation and radiation appeared in IIC2 zone ( Figure 5), which was consistent with the previous result that the browning trend in this region was associated with a delayed vegetation green-up date, likely affected by pre-monsoonal droughts [23,64,65].
In more humid areas on the mid-eastern and southeastern TP, the impact of temperature became prominent [66][67][68][69]. Furthermore, in forested areas in the southern TP, solar radiation appeared to be the main climatic driver of vegetation change, which is in good agreement with previous findings [70]. Interestingly, a recent study also reported that the solar radiation became an increasingly important factor for vegetation growth in the wet region of South Asia [71]. However, in the northwest TP, the partial correlations between GNDVI and major climatic factors were insignificant; one potential explanation for this is that other factors in this region play an important role in changing greenness such as ephemeral plants [72,73], frozen soil [61] and winter snow [74]. Apart from the direct precipitation from the atmosphere, permafrost soils, winter snow, and glacier melting might also provide an indirect supply of soil moisture in this high elevation region [66], but more underlying mechanism for vegetation change in this region need further research.
Although vegetation growth at high altitudes and latitudes is very sensitive to temperature [5,47], the monotonic greening trend of vegetation greenness across the entire TP and in most elevation bins over the past 16 years could not mirror the warming hiatus [32]. Spatially, even though vegetation greenness shifted from greening to browning in the half of zones (Figure 4), vegetation greenness in these zones were also sensitive to other factors (Figures 6 and 7). These results suggest that the recent temporary warming hiatus might not be useful as a direct indicator of a major trend shift in vegetation greenness across the TP.
Furthermore, the positive impact of increasing precipitation on vegetation growth has been obviously offset by an increase in anthropogenic disturbance in the northeastern TP, while a higher degree of anthropogenic disturbance has led to a more negative impact of water deficit in the southwestern TP (Figure 9). This result is not only in line with the previous finding that climate change is the key factor responsible for the vegetation change at the region scale [57,61], but also indicates that local degradation is mainly caused by the dual effects of climate change and anthropogenic disturbance, particularly in areas characterized by higher degree of anthropogenic disturbance.
Human footprint pressure increased rapidly overall on the TP [34]. At the same time, a series of projects have been implemented to facilitate ecosystem adaptation and restoration on the TP, especially after 2004, which reduced human footprint pressure to some extent. Several assessments of ecological projects also pointed out that a certain efficacy has been achieved [33,[75][76][77]. The reduction of livestock was considered as an important effective ecosystem management strategy in these projects over the TP. However, vegetation greenness dynamics mismatched to the change in total number of livestock on the TP [33]. In this study, we demonstrated that both vegetation greenness and total number of livestock increased between 2000 and 2005, whereas the total number of livestock decreased after 2005, but vegetation greenness continuously increased until 2010 in most zones of the TP, except for southwestern ecosystem zones (IC2 and IIC1 zones) (Figures 2 and 4 and Figure S4). These further indicate that increases in livestock number weakened the vegetation greenness greening or exacerbated browning but might not overwhelm the impact of climate change, which corresponds with previous research [61]. In addition, within the context of water deficit, even though both the total number of livestock and vegetation greenness decreased after implementation of these ecological projects in the southwestern ecosystem zones (Figures 2 and 8) [40,61], a larger degree of human footprint pressure led to a more negative impact of the climate condition on vegetation growth (Figure 9). Thus, we suggest that future research should pay more attention to evaluate the ecosystem natural capacity and improve animal husbandry management strategies on the areas with moderately and high human footprint pressure, in order to mitigate the negative effect caused by changing climate to this extremely fragile ecosystem.

Conclusions
We revisited the spatiotemporal patterns of vegetation change and its major underlying drivers on the TP during 2000-2015 based on the updated NDVI form MODIS Terra-C6. Overall, the area-averaged GNDVI exhibited monotonic greening trends on the TP during the past 16 years, with a rate of 0.0011 year −1 (p < 0.01). Over 70% of vegetated areas across the TP displayed greening trends apart from the southwest TP. However, GNDVI in most ecosystem zones underwent major shifts around 2010. Spatially heterogeneous temporal trends in major climate factors have induced different responses of vegetation dynamics on the TP. We demonstrated that water availability was the main factor of vegetation growth in the northeastern and southwestern TP, while the impact of temperature became prominent in the mid-eastern TP, and solar radiation and temperature were the two primary climatic drivers in the southern TP. Additionally, along with anthropogenic disturbance increase, the greening trend declined in the northeastern TP while the browning trend increased in the southwestern TP, suggesting that impacts of anthropogenic disturbance on vegetation change cannot be ignored and more attention should be paid to the areas with moderate and high human footprint pressure. This study provides an objective assessment of vegetation greenness and its climatic and anthropogenic drivers, which are of scientific significance for better environment management on the TP.