Factors Driving Changes in Vegetation in Mt. Qomolangma (Everest): Implications for the Management of Protected Areas

: The Mt. Qomolangma (Everest) National Nature Preserve (QNNP) is among the highest natural reserves in the world. Monitoring the spatiotemporal changes in the vegetation in this complex vertical ecosystem can provide references for decision makers to formulate and adapt strategies. Vegetation growth in the reserve and the factors driving it remains unclear, especially in the last decade. This study uses the normalized difference vegetation index (NDVI) in a linear regression model and the Breaks for Additive Seasonal and Trend (BFAST) algorithm to detect the spatiotemporal patterns of the variations in vegetation in the reserve since 2000. To identify the factors driving the variations in the NDVI, the partial correlation coefﬁcient and multiple linear regression were used to quantify the impact of climatic factors, and the effects of time lag and time accumulation were also considered. We then calculated the NDVI variations in different zones of the reserve to examine the impact of conservation on the vegetation. The results show that in the past 19 years, the NDVI in the QNNP has exhibited a greening trend (slope = 0.0008/yr, p < 0.05), where the points reﬂecting the transition from browning to greening (17.61%) had a much higher ratio than those reﬂecting the transition from greening to browning (1.72%). Shift points were detected in 2010, following which the NDVI tendencies of all the vegetation types and the entire preserve increased. Considering the effects of time lag and time accumulation, climatic factors can explain 44.04% of the variation in vegetation. No climatic variable recorded a change around 2010. Considering the human impact, we found that vegetation in the core zone and the buffer zone had generally grown better than the vegetation in the test zone in terms of the tendency of growth, the rate of change, and the proportions of different types of variations and shifts. A policy-induced reduction in livestock after 2010 might explain the changes in vegetation in the QNNP.


Introduction
Protected areas are designed to represent or sample regional biodiversity and are vital for safeguarding the biodiversity and maintaining crucial ecosystem services [1,2]. A total of 238,563 protected areas have been designated globally as of 2018, covering 14.9% of the Earth's land surface [3]. Research has shown that about 40-50% of global protected areas It is important to advance our understanding of the dynamics of the vegetation in the QNNP region. This work seeks to answer the following questions: How has vegetation in the reserve changed in recent years? What are the factors that have influenced this change? Has the establishment of the reserve protected the vegetation? In which area should the measures for vegetation protection be reinforced in the future? This study uses the NDVI to detect variations in the vegetation in the QNNP between 2000 and 2018. We examine the fluctuations induced by climate change and human activities, compare the changes in vegetation in different zones of the reserve, and discuss the factors influencing changes in vegetation in the QNNP.

Study Area
The QNNP is located close to the frontier of Nepal and China, and it includes Tingri, Gyirong, Nyalam, and part of the Dinggye counties, covering an area of 33,819 km 2 ( Figure 1). There are five mountains over 8000 m high in the QNNP, including Mt. Qomolangma (8844.43 m). The middle part of the terrain is a flat valley, with the Yarlung Zangbo River flowing across it, and the northern and western regions are much higher.  [46]. Land cover data obtained from the study by Nie [47]. Core zones (red color in (c)) are well-preserved natural ecosystems with concentrated distributions of rare and endangered animals and plants.
Weather and vegetation in the northern and southern slopes differ markedly. The southern slopes are more influenced by the humid and warm monsoon from the Indian Ocean. Precipitation over the course of a year peaks twice without a prominent wet and dry season. Warm moisture travels through deep valleys of the Himalayas, causing them to receive more rainfall; therefore, forests (2.38%) and shrublands (4.66%) are the main land cover types of these regions (Figure 1a). With much colder and dryer conditions, the northern slopes are characterized by grassland (44.02%) and sparse vegetation (22.09%) (Figure 1a), which are the main land cover types in the QNNP. On the northern slope, precipitation is mainly concentrated from June to September. Wetlands (1.84%) and croplands (0.35%) are mainly located near the main streams and tributaries of the Pumqu River (Figure 1a). The QNNP is sparsely populated. Urban and built-up areas are mainly located outside the preserve, near the boundaries. Within the reserve, their occupations are less than 0.1%. Industrial-scale activities are strictly controlled in the reserve and grazing and road construction are the major human activities [34].

Data Sources and Pre-Processing
The NDVI has been widely used to detect changes in vegetation and is an effective indicator of its health. To analyze variations in vegetation in the QNNP, we used the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices version 6.0 (from 2000 to present) from the United States Geological Survey (https://lpdaacsvc.cr.usgs. gov/appeears/ (accessed on 20 November 2019)), a 16-day maximum composite dataset with a spatial resolution of 250 m. The website provides online services on image mosaics, project transformation, format conversion, and image masking. A total of 434 original images (a single tile covered the study area, h25v06, with 23 images per year, and only 20 images in 2000) during 2000-2018 were acquired from the website. A Savitzky-Golay filter in the TIMESAT software was used to eliminate residual contaminations caused by clouds and snow. The TIMESAT software was developed by Jönsson and Eklundh [48] to analyze time series of sensor data and is an effective method to fix abnormal values in the MODIS preprocessing procedures [49,50]. The Savitzky-Golay filter is a filter method based on local simplified least-squares fit convolutions for smoothing and computing the time series [51]. After the filtering process, the NDVIs of the growing season (May to September, DOY 145 to 304; 10 images per year) were used as the inputs for trend detection (Section 2.3.1) and shift point detection (Section 2.3.2).
To analyze variations in the NDVIs of different land cover types, we used the QNNP landcover product generated by Nie [47], which reached a global accuracy of 93.87%. The spatial resolution of the product was 30 m. The classification system of the QNNP was designed based on the framework of the land cover classification system (LCCS) [52]. We divided the types of vegetation of the product into six main classes (sparse vegetation, grassland, forest, shrubland, cropland, and wetland), and masked the non-vegetation areas (construction land, glaciers, water, and bare lands). The specific definition of every vegetation type, and the dominant species/constructive species can be found in the study by Nie [47]. The product was interpolated to 250 m using nearest neighbor resampling to match the resolution of the NDVI product when analyzing the trends of NDVI of different vegetation types.
The China Meteorological Forcing Dataset (CMFD) was used to analyze variation in climate in the reserve [53]. The spatial resolution of this dataset was 0.1 degrees (approximately 11 km), and it was generated by combining data observed from the ground and several satellite-derived products. Ground-observed climate data from the National Meteorological Information Center (http://data.cma.cn/ (accessed on 20 November 2019)) were used to test the reliability of the product. There are only two stations in the QNNP: in Nyalam (station number 55655) on the southern slope of the reserve, and in Tingri (station number 55664) on its northern slope. The China Meteorological Forcing Dataset was in good agreement with the ground-observed data in terms of precipitation but overestimated the temperature of Tingri and underestimated that of Nyalam during the growing season ( Supplementary Information, Figure S1). When calculating partial correlation coefficient between climate and NDVI, the climate data was resampled to 250 m to match the resolution of NDVI and the resampled landcover product.
To analyze impact of human activities on vegetation in the QNNP, the statistical yearbook of Shigatse (2000-2016) was obtained from the local government. The yearbook covers information on husbandry, industry, transportation, construction, and business. We mainly used livestock numbers for our analysis.

Trend Analysis and Partial Correlation Analysis
A simple linear regression model was used to analyze variations in the NDVI during the past 19 years, and the slope of the NDVI was calculated using the least-squares method: where NDVI i is the annual mean NDVI in the growing season in the ith year. A positive number indicates a trend of greening while a negative number indicates a browning trend. We used the F test (p < 0.05) to test significance [54,55].

Break Point Detection
We used the BFAST algorithm [27] in R language (https://cran.r-project.org/web/ packages/bfast/index.html (accessed on 20 November 2019)) to explore shifts in the trend of the NDVI in the reserve. The algorithm is as follows: This algorithm decomposes time series Y t (i.e., hydrology, climatology, and economics) during period t into trend (T t ), seasonal (S t ), and remainder components (e t ); it then detects and characterizes abrupt changes in the time series. In our study, Y t denotes the 16-day NDVI time series in the growing season during 2000-2018.
We used BFAST01 implementation to detect either zero or one break point in the time series. Land use and land cover change are relative limited in our study area, and thus the detected break was likely to represent the most ecologically relevant shift in a time series [56]. We used the ordinary least-squares (OLS) residuals-based MOving-SUM (MOSUM) test [57] to evaluate whether break points occurred. We considered only points where both segments were significant (p < 0.05). According to de Jong et al. [58], six types of changes occur: monotonic greening, greening with setback, browning to greening, monotonic browning, browning with burst, and greening to browning.
We established buffer zones in ArcGIS software (version 10.3.1) to detect the influence of the natural reserve on the protection of vegetation. Considering the resolution of MODIS data and the previous research [12,59], the scope of the buffer zone was 25 km on both sides of the boundaries of the reserve, and we used the interval of 5 km when calculating variations in the NDVI. We did not conduct a buffer analysis of boundaries of the protected area that overlapped with national boundaries (Supplementary Information, Figure S2) because a large part of these buffer areas is in Nepal and have completely different vegetation types. Because we lacked a precise landcover product outside the QNNP, we masked pixels with an NDVI for the growing season of lower than 0.1 to minimize noise from ice, snow, water, sand, and stone when comparing variations in the NDVI between inside and outside the boundary of the reserve.

Factors Driving Changes in Vegetation
To analyze the impact of climate change on variations in the NDVI, we calculated the partial correlation coefficient (PCC) between the NDVI, and precipitation, temperature, and radiation. When calculating the association between the NDVI and a given factor, the two other factors were eliminated as control variables. The formula is as follows: where R 12,34 is the partial correlation coefficient of variables 1 and 2, and variables 3 and 4 are the controlled variables, R 12,3 refers to the partial correlation coefficient of variables 1 and 2, and variable 2 is the controlled variable, and R 12 is the Pearson correlation coefficient of variables 1 and 2, and the other variables have similar meanings as before. The absolute value of the PCC is used to determine the best-fitting time effects. When calculating the PCC, the effects of time lag and time accumulation were considered. Previous work has shown that the time effects are generally shorter than a quarter of a year [60,61]. Therefore, we considered the effects of time lag and time accumulation up to a maximum of 96 days (approximately 3 months). Figure 2 shows the 28 schemes considered in this study. In Figure 2, NDVI (t,t+16) is the maximum NDVI value during the period (t~t+16). Take the scheme marked with red color (L-16/A-16) as an example, the corresponding climate factors are the accumulate values during the time periods (t−32~t). To quantify the overall effect of climate change on variations in vegetation, a multiple linear regression model was developed: where A, B, and C are the regression coefficients, and ε is the error term. TMP, PRE, and SR are the adjusted time series of climatic factors (temperature, precipitation, and solar radiation, respectively) with the best effects of time lag and time accumulation identified through Equation (3). The absolute value of determination coefficient (R 2 ) of Equation (5) was used to quantify the overall explanation of climatic factors for variations in the NDVI.

Spatial Distributions of Tendency and Shift in NDVI in QNNP in 2000-2018
The annual average NDVI exhibited large spatial heterogeneity in the QNNP (Figure 3a). Large NDVI values were mainly concentrated in the southern region of the preserve, where forests and shrubland dominated. The NDVI values below 0.3 were widely distributed in the middle region of the QNNP. Elevations in these regions were higher than the southern and northern regions, and they were mainly covered with grassland and sparse vegetation. Alongside the rivers, where wetlands dominate, the NDVI values were much higher than in the surroundings. For different vegetation types, forests had the largest NDVI (over 0.7), followed by shrubland and cropland; sparse vegetation and grassland had relatively low NDVI values. The annual mean NDVI in the QNNP generally showed a tendency of growth (0.0008 per year, p < 0.05) during 2000-2018. For the entire QNNP, vegetation with a significant and an insignificant trend of greening accounted for 16.66% and 61.77% of the total land, respectively (blue color in Figure 3b). Vegetation with a significant trend of decrease accounted for only 1.27% and was mainly concentrated in the east part of the region (Dinggye county) and the land alongside roads. Using the BFAST model, we analyzed the shifts in the NDVI and their types (interruptions or reversals) in the QNPP during the past 19 years. Although the vegetation in the QNNP showed an overall tendency of growth, this overall period can be subdivided into two: 2000-2010, with an overall tendency of reduction (p < 0.05), and 2011-2018, with an overall tendency of increase (p < 0.05). Shift points in 2010 were also discovered for different vegetation types ( Figure 4): After 2010, the tendencies of growth of forest and shrubland increased, the trend of cropland and sparse vegetation changed from one of decrease to that of increase, and the trends of decrease of wetland and grassland decreased compared with those in the first 10 years of the study period.
In terms of spatial distribution, the NDVI underwent significant shifts, accounting for 38.34% of the entire QNNP ( Figure 5). A total of 29.24% of the shifts occurred during 2010-2011, and 15.57% in 2004. Among all the points with significant shifts, monotonic greening and greening with setback were widely distributed within the study area (green in Figure 5), accounting for 35.66% and 19.22%, respectively. The proportion of browning with burst (orange in Figure 5) was mainly distributed near rivers, roads, and the eastern region of Dinggye county. Notably, the reversal points, the points showing the transition from browning to greening (17.61%), accounted for much larger proportions of the QNNP than those representing the transition from greening to browning (1.72%). Points representing the transition from browning to greening were mainly distributed alongside rivers, lakes, and roads. Shifts from browning to greening mainly occurred before 2011 (90.38%), while a large number of shifts from greening to browning occurred after 2010 (85.11%).

Climate Change and Its Impact on Variations in NDVI
To analyze the factors that might have affected the vegetation in the QNNP, we examined the response of the vegetation to the changing climate by considering the effects of time lag and time accumulation. We first analyzed the climatic trend of the QNNP during the study period ( Figure 6). The entire QNNP exhibited a drying-warming trend in the growing season. During 2000-2018, most regions of the reserve showed a warming trend, especially the central part. Precipitation decreased significantly in the middle and southwest of the reserve and wetting in the northwestern and eastern parts of the reserve was insignificant. Radiation showed a trend of increase in the QNNP. Time effects were examined by analyzing the PCC between the NDVI and climate variables over different periods (Table 1). The NDVI time series had the strongest PCC with temperature and precipitation at L-16/A-16 (cumulative over 16 days with 16 days of lag), and the maximum PCC values were 0.82 and 0.70, respectively. The response of the NDVI to radiation showed no time lag but a strong time accumulation effect (L-0/A-96, accumulated over 96 days with no day lag), and the maximum PCC was 0.70. Table 1. Partial correlation coefficients between the NDVI and climatic factors when considering the time effect.  We used the data in Table 1 in the BFAST model to detect changes in the temperature and precipitation at L-16/A-16, and those in radiation at L-0/A-80 ( Figure 8). We found no shifts in 2010, which does not help to explain the changes in the NDVI in this year. There were other reasons for variations in the NDVI in the QNNP.

Impact of Human Activities on NDVI in QNNP
We first analyzed the variations in the NDVI in different zones of the QNNP to test the effectiveness of the protection of the reserve. The results showed that the tendencies of the growth of the vegetation in the core zone and the buffer zone were higher and more significant (p < 0.01) than those of the test zone ( Table 2). The core zone had the highest annual rate of change (0.42%), and the buffer zone had the highest total rate of change. With regard to the shifts in the NDVI, the core zone and the buffer zone showed a trend from browning to greening while the test zone showed a decrease in vegetation with a positive break (Figure 9).    We further made buffers outside and inside the boundary of the reserve to study the differences between the variation in vegetation near the boundaries of the QNNP ( Figure 10). Overall (in terms of all the indices), the vegetation growing inside the QNNP was generally better than that outside. The largest proportion of a significant decrease inside the QNNP was within 5 km of the border. With increasing distance, the proportion of significant/insignificant increases in the NDVI became larger closer to the core protected areas, as did the rate of change/annual rate of change in the NDVI. This characteristic was especially prominent for vegetation within 15 km of the border. In terms of shifts in the NDVI, the NDVI inside the reserve had larger proportions of monotonic greening areas and smaller proportions of monotonic browning and browning with burst. Compared with the area outside, a larger proportion of the vegetation inside the reserve had recovered, with trends of the NDVI shifting from that of decrease to one of increase. We determined the variation in the livestock in the reserve because it feeds on plants and thus directly influences the vegetation and, in turn, the variations in the NDVI. As shown in Figure 11, the livestock number in the QNNP has increased drastically during 1952-2018. It was only 379,000 in 1952, while in 2008 it increased to more than 950,000. During our study period, the number of livestock showed a monotonic decrease in Tingri, and an increasing-to-decreasing trend in Dinggye and Gyirong. In Nyalam, the livestock slowly decreased and then increased in 2012. For the entire reserve, the livestock numbers showed a significant trend of decrease during 2010-2018 ( Figure 11). This was consistent with the shift in the NDVI of the entire reserve in 2010.

Variation in Vegetation in the Reserve
In the past 19 years, the NDVI in the reserve has shown a tendency of growth (0.0008/yr), which is also the case in the Koshi River Basin [62], the Himalayan region [6,39], and the TP [63,64]. Compared with past work that has focused on vegetation in the QNNP [40][41][42][43], we have specified variations in the NDVI after 2010, and find that the NDVI of the QNNP has not undergone a linear change. It changed from a trend of decrease to one of increase around 2010, which has not been reported for the reserve before. The patterns of the NDVI detected by BFAST in the first period were consistent with previous studies [40,42,44,65]. For instance, we found that forests and shrublands showed a trend of monotonic greening while the other vegetation types showed a trend of decrease in the first period. Nie et al. [42] also found that, during 1998-2009, forests and shrublands had a larger proportion of significant increase than significant decrease, whereas the proportions of significant increases in grassland, sparse vegetation, and croplands were smaller than the proportion of significant decreases. Ma et al. [65] found that during 2000-2009, the vegetation in the northern slopes mainly showed a trend of decrease, and over 40% of the vegetation showed different degrees of degradation. These places were dominated by grassland and sparse vegetation.
The BFAST can help detect shifts beyond the linear regression model and can provide more detail on the variations in the NDVI. Although the linear regression showed that all the vegetation types had exhibited tendencies of growth in the past 19 years (Table 2), grasslands and wetlands still need to be attended to. They showed a slight trend of decrease during 2010-2018 (Figure 4). Figure 12 shows that among the points of insignificant increase and significant increase, a large number of points exhibited browning with burst, and a greening-to-browning trend. Grassland and sparse vegetation accounted for large proportions of these two shift types. If protections are not reinforced in the future, the vegetation in these places will decrease in the long run.

Response of NDVI to Climatic Variables
Temperature and radiation were dominant among the climatic factors influencing the variations in the NDVI in the QNNP. This is consistent with a previous study that focused on the response of the NDVI to climatic factors in different regions of the TP [54] and other energy-limited middle-to-high latitudes of the Northern Hemisphere [60]. Increasing temperature and radiation are important factors for an increase in the greenness of the vegetation, and can significantly affect the vegetation of the TP, increase its productivity, height, and the duration of its growing periods [66][67][68]. The impact of radiation on the vegetation was not always positive. In the regions with limited water, high solar radiation increased the soil temperature, thus accelerating the transpiration of moisture and inhibiting vegetation growth [69,70]. Li et al. [55] found that in the southwestern TP, increasing temperature together with a growing amount of solar radiation adversely affected vegetation growth in the case of an insufficient supply of water. The effects of time lag and time accumulation were noted in the relationship between the NDVI and climatic variables in the QNNP and have been reported for the Yumco Basin [54], Naqu [71], and TP [72]. Plants need an accumulation of climatic factors to initiate their lifecycles and to obtain a sufficient amount of nutrients from the soil, heat, and moisture from the environment [60].
Climatic factors could explain 44.04% of the variations in the NDVI in the QNNP. This is below its global average value (63%) in a previous study [60]. This region is fragile to disturbances (i.e., over-grazing, land use changes, drought, fire, and plant diseases) [61]. In the context of shift detection, Li et al. [73] found adverse trends in precipitation in the entire southwestern TP and attributed the shift in the NDVI of the region to changing precipitation. However, in the QNNP we did not find any shift in climatic factors in 2010, and this cannot explain the shift in the NDVI in this year. Other factors should thus be taken into consideration in examining the variations in the NDVI in this region.

The Effectiveness of Conservation
The effectiveness of the natural reserve in the QNNP can be verified from the following: First, the buffer zone and the core zone generally had better vegetation conditions. The results showed that, since 2000, these zones had first decreased and then increased. The test zone showed a decrease in both periods, but after 2010 its rate of reduction decreased. Second, near the boundaries of the reserve, the vegetation inside the reserve generally grew better than that outside the boundary. Third, the result of the BFAST showed that 17.61% of the area had turned from brown to green and was not evenly distributed over the entire reserve but was mainly located alongside roads and rivers, or areas featuring intense human activities.
We further analyzed the tendencies of the NDVI for different types of vegetation in the three zones to confirm that the higher tendencies of NDVI growth in the buffer and core zones were caused by the establishment of different levels of protected zones, and not simply because the large areas of these zones were covered by forests and shrublands. The results showed that, except for the croplands (mainly located in the test zone), all types of vegetation had a higher tendency of growth, annual rate of change, and total rate of change in the buffer zone and the core zone than in the other zones ( Table 2). In terms of the proportions of types of change (Figure 13), most vegetation had better growth conditions in the buffer and core zones. Land falling into categories of increase (significant/insignificant increase, monotonic greening, greening with setback, and browning to greening) generally constituted a higher proportion of land in these zones, with a higher proportion of land falling into categories of decrease (significant/insignificant decrease, monotonic browning, browning with burst, and greening to browning) in the test zone. Therefore, the core and buffer zones provided more effective protections than the test zone. Previous studies also indicated the positive impact of the establishment of reserves in the TP. For instance, Zhang et al. [10] found that the NPP in 76% of samples inside nature reserves and 82% of samples inside national nature reserves in the TP was higher than that of the corresponding samples outside the reserve. Through comparing two reserves in the TP, Li et al. [74] found that longer-protected reserves exhibit higher stability to climate change. The establishment of natural reserves on the TP has reduced human activity [34]. According to the Chinese government's regulations on nature reserves and the local policies, human activities (i.e., production facilities) were considerably restricted in the reserves, especially in the buffer zone and the core zone. Without permission, people are strictly prohibited to enter the core zone. While in the buffer zone, only non-destructive scientific research and observation are allowed. In the core zone and the buffer zone, no production facilities are allowed to be built. In the test zone, the production facilities cannot pollute the environment or destroy the landscapes. In the protected area, the core zone also constitutes areas featuring ecological engineering (i.e., the Natural Forest Protection Project, 1998, and the Green for Grain Project, 2002) [42].

The Impact of Livestock Decreasing
The policy-induced reduction of livestock could help to explain the NDVI increase in the reserve. Grazing has been practiced on the TP for three millennia [75]. However, during the past 50 years, due to the improved infrastructure, the growth of the population, and the market demand, grassland-livestock industry has quicky developed [76]. According to the study of Yu et al. [77], in 2010, the livestock number in the TP increased to nearly double the estimated capacity of all available pastures. Overgrazing is threatening the ecological security of the TP, becoming one of the major factors causing grassland degradation [78,79]. In our study area, during 1952-2008, the livestock number increased by 96.16%. This may cause a huge impact on the environment. Overgrazing can reduce the soil-water content, as well as nitrogen and phosphorus, causing a reduction in vegetation cover, vegetation height, and aboveground biomass [80], and changing the plant community structure (i.e., dominant species, richness, abundance, and life forms) [81]. Research on three regions of river sources has shown that increasing livestock numbers have eliminated the positive effect of climatic factors during 1984-1993 [82]. Field experiments have showed that it could influence the grassland to a greater extent than increasing temperature [83]. In the field experiment of Wang et al. [84], heavy grazing, rather than the warming climate, led to the degradation of alpine meadows.
We found that about 29.2% of the shifts occurred during 2010-2011, and the tendency of livestock growth in the reserve reversed around 2009 and has declined drastically since 2010. This variation is consistent with the implementation of engineering policies that have reduced livestock in Tibet. In 2011, the government reinforced the policy by introducing an award allowance payment policy, and has spent 1.6 million Yuan every year in pastoral areas of China to compensate for prohibitions on grazing, to reward efforts to achieve pasture-livestock equilibrium, and to provide production subsidies for herdsmen [85]. In addition, 15.57% of the shifts happened in 2004. The Chinese government enforced the ecological project, Grazing Withdraw Program (GWP), in 2003, bringing the implementation of many large-scale ecological projects, such as fencing degrading grassland, ecological compensation, and restoring the land cover vegetation [86]. The detection of the NDVI variation in the Qilian Mountain Region [21] also found positive abrupt change increased significantly since 2005, after the implementation of the GWP. Studies in the central TP found that the implement of ecological protection and restoration projects have mitigated grassland degradation, and even reversed the degradation in some areas [87]. Therefore, we thought that the policy-induced livestock decrease after 2010 might be one potential explanation for the NDVI shift in 2010.

Future Protection of Vegetation in the QNNP
The vegetation outside the reserve had a lower tendency of growth in the NDVI, and higher ratios of significant and insignificant trends of decrease. Human activity is one of the major factors leading to vegetation reduction near the boundaries. Li et al. [34] found that the human footprint has increased more significantly outside the QNNP reserve. Studies have also shown that the protected areas attract, rather than repel, human settlements by providing economic and occupational opportunities, and this threatens the effectiveness of conservation [88]. From Figure 11, we can see the population in the four counties showed a steady increase. Land cover changes may be another reason for a decrease in the NDVI [17]. The vegetation conditions outside the protected areas and the factors driving them should receive more attention because the environmental conditions and human activities in them can influence the natural resources within the reserves by affecting their biological and physical processes [8,89,90], and can create a ring of disturbance outside the reserve that isolates the protected area from its surroundings [81].
Apart from the overall tendencies of growth in the NDVI in the QNNP, reductions in it were observed in the wetlands at high altitudes and the forests alongside roads on the southern slope. We found a trend of decrease in wetlands during 2000-2010. A previous study has shown that during 2000-2008, degraded wetlands occupied 71.3 km 2 , or 1.8% of all wetlands in the QNNP [91]. The wetlands east of Dingye county underwent significant reductions in vegetation (red color in Figure 3). Such wetlands, especially flatlands along the riverbank, are perennially impacted by sand erosion [92]. In the valleys of the southern slopes, the destruction of forests, and the construction of roads and other infrastructure, together with frequent rainfall and unstable geological conditions, have triggered landslides, thus inflicting further damage on the surrounding vegetation. In addition, anthropogenic disturbances, i.e., fire, over-grazing, and lopping, hinder the natural succession of some species of trees by preventing the regeneration of seedlings and saplings from the understory of the forest [93].

Shortcomings of This Study
Our study gave a preliminary finding on the vegetation response and the impact from humans and the climate in the QNNP. The coarse resolution of the climate data could have affected the calculated variations in the vegetation due to climate change. In future work, weather stations in the surroundings could be gathered to generate more precise climate data.
With respect to the human influence, we considered only the livestock number in the region, while the surrounding industrial-scale activities (i.e., the construction of roads and built-up areas) could also bring impact on the environment. Although the impacted areas of those human activities are limited, mainly alongside the roads or surrounding the built-up areas, their influences are severe. However, the current study cannot differentiate the impacts from different human activities. In future studies, more human activities could be taken into consideration. Spatializing the statistical data can help to clarify the relevant influential factors in the different areas. Our work did not take the impact of the land use and the landcover change into consideration because it is difficult to capture by moderate resolution image. The applications of the satellite data with higher resolutions (i.e., Landsat, Sentinel) could help us to clarify the impact from the land use and the land cover change. In addition, more collaborative research involving the local stakeholders and researchers (i.e., a questionnaire survey in the local communities) could help us to better understand the influence of human activities in different scales, and make the study result more convincing [94].

Conclusions
We studied the variations in the NDVI of the QNNP during 2000-2018 in this paper. We applied the break point detection to monitor variations in the vegetation in the protected areas. This algorithm can help provide detailed information on the changes in the vegetation so to implement timely protection measures. The results showed that the vegetation in the reserve had grown during the past 19 years, but its growth was not monotonic. The vegetation declined in the first 10 years, and then grew after 2010.
We found that climatic factors could account for only 40.04% of the variation in the vegetation in the QNNP. Human activity has thus affected the vegetation in even the highest region of the world. Although the linear regression model showed that all vegetation types had an increasing tendency, the result of the BFAST algorithm showed that grasslands and wetlands had decreased after 2010. The vegetation in the test zone, especially near the boundary of the reserve, should be strictly protected. Our work here can provide information for the sustainable management of the reserve.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/rs13224725/s1, Figure S1: Comparison of the ground-observed climate data and the China Meteorological Forcing Dataset. Figure S2: NDVI trends around the reserve boundary. Insert figures is NDVI trends with p value higher than 0.05 being masked out. Data Availability Statement: Data used in this study will be available upon request from the first author.