Climate Variability Rather Than Livestock Grazing Dominates Changes in Alpine Grassland Productivity Across Tibet

Alpine grasslands on the Tibetan Plateau, being vulnerable to environmental and anthropogenic changes, have experienced dramatic climate change and intensive livestock grazing during the last half-century. Climate change, coupled with grazing activities, has profoundly altered alpine grassland function and structure and resulted in vast grassland degradation. To restore degraded grasslands, the Central Government of China has implemented the Ecological Security Barrier Protection and Construction Project since 2008 across the Tibetan Autonomous Region. However, the relative effect of climate change and grazing activities on the variation in alpine grassland productivity is still under debate. In this study, we quantified how aboveground net primary production (ANPP) varied before (2000–2008) and after (2009–2017) starting the project across different alpine grasslands and how much variance in ANPP could be attributed to climate change and grazing disturbance, in terms of temperature, precipitation, solar radiation, and grazing intensity. Our results revealed that Tibet’s climate got warmer and wetter, and grazing intensity decreased after starting the project. Mean ANPP increased at approximately 81% of the sites, on average from 27.0 g C m–2 during 2000–2008 to 28.4 g C m–2 during 2009–2017. The ANPP positively correlated with annual temperature and precipitation, but negatively with grazing intensity for both periods. Random forest modeling indicated that grazing intensity (14.5%) had a much lower influence in controlling the dynamics of grassland ANPP than precipitation (29.0%), suggesting that precipitation variability was the key factor for alpine grassland ANPP increase across Tibet.


INTRODUCTION
Climate change and human activities are primary drivers for changes in terrestrial ecosystems globally (Haberl et al., 2007;Chen et al., 2013;Tong et al., 2018), especially for unprecedented changes in ecosystem service and function (Knight and Harrison, 2012;Seiferling et al., 2014;Erb et al., 2018). A meta-analysis by Lin et al. (2010) found that climate warming can promote global vegetation productivity. Precipitation can also accelerate or decelerate vegetation growth in most terrestrial ecosystems, depending on its frequency and timing (Wu et al., 2011). Meanwhile, human activities can affect ecosystem response and feedback to climate change, via biomass utilization, biofuel consumption, land use, and land cover change (Haberl et al., 2007;Kröel-Dulay et al., 2015). For example, human disturbances might push an intact ecosystem away from its successional trajectory, alter species assembly, and turnover, and make it much more vulnerable to climate change than before (Kröel-Dulay et al., 2015).
Aboveground net primary production (ANPP) is a critical ecosystem service in maintaining global carbon balance (Ruimy et al., 1999;, often used to forecast ecological consequences under anthropogenic and climatic perturbations. A better understanding of how ecosystem ANPP responds to climate change and human activities can help mitigate environmental damages and optimize ecosystem management (Zhou et al., 2018). A considerable volume of literature tried to identify and quantify the relative influences of climate change and human activities on ecosystem productivity with various methods (Paudel and Andersen, 2010;Erb et al., 2018;Li L. et al., 2018), such as manipulative experiments, traditional statistical analysis, and residuals-trend modeling (Li L. et al., 2018). Recently, random forest modeling, with high accuracy and robust efficiency, is increasingly being used to quantify predictors' relative importance (Gill et al., 2017;Huang and Xia, 2019). The random forest modeling can also account for interactions and non-linear relationships between predictors, obstacles (the overfitting problems) and deal with data noises (Heung et al., 2014). Therefore, scientists often recommended random forest models for processing high-dimensional and -correlated datasets (Breiman, 2001).
Due to high elevations, alpine grasslands widely distributed on the Tibetan Plateau are susceptible to climate change and human activities, as vegetation in Arctic and Antarctic (Yao et al., 2012;Yang et al., 2017;Li et al., 2020). Climate change, combined with human activities, has significantly reshaped the structure and function of the Tibetan alpine grasslands (Wang et al., 2005), and resulted in about 0.4 × 10 6 km 2 grassland degradation in the1990s, which accounted for 33% of total grassland area on this plateau (Long et al., 2009). Therefore, rational utilization, restoration, and conservation of alpine grasslands under changing climate and anthropogenic disturbances have attracted increasing attention since 2000 (Harris, 2010;Yu et al., 2012;Dong et al., 2020). The central government of China initiated several ecological projects, like the "Ecological Security Barrier Protection and Construction Project" started since 2008, and the "Compensation and Rewards to Herders for Natural Grassland Conservation" since 2010, to promote the recovery and restoration of degraded alpine grasslands on the Tibetan Plateau (Wang et al., 2017). The implementation of these projects has eased the pressure of human disturbances (mainly referring to livestock grazing hereafter) (Fan et al., 2015) and also increased the alpine grassland productivity (Chen et al., 2014;Xu et al., 2016). However, the relative contributions of climate change and human activities to alpine grassland productivity are still unclear. Moreover, it is pertinent to examine the extent to which ecological projects and policies can affect alpine grassland productivity.
To fill these gaps in current research, we collected historical records of climatic variables and the Normalized difference vegetation index (NDVI) data before (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and after (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) the starting of the "Ecological Security Barrier Protection and Construction Project." These data were mainly used to simulate the changes in alpine grassland productivity during the two subperiods. We also conducted field measurements of aboveground biomass between fenced and grazed pastures in different alpine grassland types to validate the models' outputs. By comparing the difference in alpine grassland productivity, it is possible to clarify the effects of the "Ecological Security Barrier Protection and Construction Project." We also collected the livestock numbers from statistic yearbooks of 2000-2017 to quantify grazing intensity over time and across space. Thus, it is also possible to quantify how much the variance in alpine grassland productivity can be explained by the corresponding climatic (refers to temperature, precipitation, and solar radiation) and anthropogenic (refers to grazing intensity) factors.

Study Area
The Tibetan Autonomous Region (hereafter Tibet) is located in the southwest of China (Figure 1), covering about 1.22 million km 2 . With an average elevation of over 4,000 m above sea level, Tibet has a cold (mean annual temperature from 2000 to 2017 was -1.1 • C) and dry (mean annual precipitation from 2000 to 2017 was 307.8 mm) climate. This climate shaped various sensitive and vulnerable ecosystems. Alpine grasslands mainly distribute in northwestern places at high elevations, covary, and coevolve well with zonal climates . From southeast to northwest, grassland types vary from humid alpine meadow (AM) dominated by Kobresia pygmaea to semiarid alpine steppe Frontiers in Ecology and Evolution | www.frontiersin.org (AS) dominated by Stipa purpurea, and to arid alpine desertsteppe (ADS) co-dominated by S. purpurea and S. glareosa.

Data Collection
Normalized difference vegetation index (NDVI) is a remote sensing technique widely used in the regional ecosystem monitoring and evaluation. This study used the moderateresolution imaging spectroradiometer (MODIS) Version 6 NDVI (MOD13A3) from 2000 to 2017, with 1 km spatial resolution and 1 month time interval. NDIV data were downloaded from the National Aeronautics and Space Administration agency 1 . The monthly NDVI data were developed using the Maximum Value Composition method (MVC) and calibrated for geometrical and atmospheric effects and cloud contamination. Grassland distribution referred to the China Vegetation Atlas with a scale of 1:1,000,000 (Chinese Academy of Sciences, 2001), from the Resource and Environment Data Cloud Platform 2 .
Daily meteorological data from 2000 and 2017 were collected from the National Meteorological Information Center (NMIC), China Meteorological Administration (CMA) 3 . Daily meteorological records were aggregated at the monthly level. Then, we interpolated monthly precipitation, temperature, and sunshine duration into raster surfaces with a 1 km spatial resolution using the ANUSPLIN 4.3 (Hutchinson, 2004). According to Allen et al. (1998), solar radiation was calculated based on geographical position and sunshine duration. It has been examined that the grid climatic surfaces match well with field observations (Chen et al., 2014;Tao et al., 2015).
In 2009, a 1,200 km transect was established across the northern Tibetan Plateau. This transect crosses exceptionally arid, semiarid, to semi-humid alpine continental climates from west to east, encompassing three primary grassland type: alpine meadow, alpine steppe and alpine desert-steppe. From 2009 to 2015, we collected aboveground biomass (AGB, g/m 2 ) fieldmeasured during peak growing season (late July to early August) along with this transect, between fenced and open (grazed) alpine grasslands. In total, we measured AGB at 224 sites in open alpine grasslands and 138 sites in fenced ones (Figure 1). At each sampling site, five 0.5 m × 0.5 m quadrats were laid at 20 m intervals along a 100 m random transect line. Plant aboveground materials were dried at 65 • C for 48 h and then weighed for AGB measurements. For grassland in Tibet, the most sampled species sprout annually in early May and senesce in late September with peak biomass generally between late July and early August. Thus, the field-measured AGB can surrogate for yearly ANPP. Finally, the dry matter was converted to carbon, assuming a carbon content of 45% (Lieth and Whittaker, 1978).
The livestock number for each county was collected from the "Statistical Yearbooks" of 2000-2017. The numbers of sheep, goats, and large herbivore animals (mainly referring to yaks, donkeys, and horses) were finally calculated as standardized sheep units (SSU), by following the algorithms of Fan et al. (2010) that one sheep equals to one SSU and one large herbivore to four SSUs.

Calculation of Grassland ANPP
The Carnegie-Ames-Stanford Approach (CASA) model was driven by satellite-observed NDVI and climate (e.g., temperature, precipitation, and radiation) and other factors, such as land-use change. The NDVI changes can reflect human harvest from plant material. NPP is mainly determined by two variables, absorbed photosynthetically active radiation (APAR) and the light-use efficiency (LUE) in the CASA model, where ANPP(x,t) represents plant growth at spatial location x and time t (g C m −2 ), and PAR(x, t), and FPAR(x, t) are total solar radiation (MJ m −2 ) and the fraction of the incoming PAR intercepted by plant incident at spatial location x and time. ε * is the maximum possible light energy conversion efficiency and was set uniformly at 0.56 g C MJ −1 . T ε 1 (x, t), T ε 2 (x, t) are the effects of temperature to ε * , and W ε (x, t) accounts for effects of water stress, also at location x and time t. ANPP can be inferred by the ratio between BNPP and ANPP (R), around 0.587 for alpine grasslands on the Tibetan Plateau predicted by Wu et al. (2010).

Calculation of Grazing Intensity
Grazing intensity (GI) was quantified by the ratio between actual livestock carrying capacity (C a ) and theoretical livestock carrying capacity (Ct): The actual livestock carrying capacity for each county was determined as follows: where C n is the number of livestock inventory in a given year, C h is the number of livestock sold in a given year. A is the available grassland area (ha).
The theoretical livestock carrying capacity for each grassland pixel was calculated as follows: where Y is grassland yield; U is the utilization rate of herbage (%), estimated as 70%; C is the proportion of the area available for pasture (%), estimated as 0.84 (Fan et al., 2010;Cao et al., 2020); H is the proportion of edible forage (%), estimated as 0.76, 0.69, and 0.76 for alpine meadow, alpine steppe and alpine desert-steppe, respectively, according to observed data in North Tibet; S is the daily feed intake per SSU, set as 1.33 kg (Fan et al., 2010); G is grazing days (d), set as 365.
In our study, we regarded yield (Y) as equivalent to the potential aboveground biomass (AGB p ), which was not grazed by herbivores. Grassland ANPP p can be estimated according to the terrestrial ecosystem model (TEM).
TEM is one of the process-based ecosystem models driven by spatially referenced information on vegetation type, climate, elevation, soils, and water availability to calculate the monthly carbon and nitrogen fluxes and pool sizes of terrestrial ecosystems. TEM can only be applied in a mature and undisturbed ecosystem without considering the effects of land use. TEM NPP (NPP p ) was calculated by the difference between gross primary productivity (GPP) and autotrophic respiration (R a ) in the monthly time step. R a is considered as the sum of maintenance respiration (R m ) and growth respiration (R g ). Monthly GPP is driven by several factors and is calculated as: where C max represents the maximum rate of C assimilation by plants in 1 month (g C m −2 month −1 ). The function of f (PAR), f (T), and f (NA) accounts for the effects of photosynthetically active radiation, temperature and relative nutrient availability, respectively, on GPP. f (LEAF) is the leaf area relative to the maximum annual leaf area and depends on monthly estimated evapotranspiration. f (CO 2 , H 2 O) is the interactive effects of atmospheric CO 2 concentrations and moisture availability to GPP. In this study, the value of C max was set as 949.6 g C m −2 month −1 for alpine meadow, 617.9 g C m −2 month −1 for alpine steppe and 251.2 g C m −2 month −1 for desert steppe (Chen et al., 2014).

Statistical Analysis
The differences in climatic variables, grazing intensity, and alpine grassland ANPP before and after the implementation of the "Ecological Security Barrier Protection and Construction Project in Tibet" were calculated by the following formula: where V represents each of mean annual temperature (MAT), annual total precipitation (AP), annual total radiation (AR), GI, and ANPP, respectively. We first calculated the zonal mean of each variable at each county in Tibet. Then, we employed a random forest (RF) regression analysis to identify the relative importance of the four variables on the variation in grassland ANPP (Huang and Xia, 2019). In the RF model, the response variable was the differences in the mean ANPP between the two periods, and the predictors were the differences in the mean values of MAT, AP, AR, and GI. The importance of each predictor variable is defined as the percentage increase in the mean square error (%IncMSE) between observations and predictions, and the decrease is averaged over all the trees to produce the final estimation for importance (Delgado-Baquerizo et al., 2017;Huang and Xia, 2019). High%IncMSE values represent high importance for each predictor. The RF model was run using R 4.0.2 with randomForest_4.6-12 4 .

Validation of Grassland ANPP and Potential Aboveground Biomass
To assure the accuracy of the ANPP and AGB p , observed ANPP data in open grazed grasslands and observed AGB data in fenced grasslands were used to validate the simulated ANPP and AGB p , respectively. The results showed that both simulated ANPP and AGB p matched well with observed records. The simulated ANPP can explain 80% of the variance in observed ANPP of open alpine grasslands (Figure 2A), and the simulated AGB p can explain 78% of the variance in observed AGB p of fenced alpine grasslands ( Figure 2B).

Dynamics of Climatic Variables and Grazing Intensity Between Two Subperiods
Both MAT and AP showed a non-significant increasing trend in the period of before (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and after the implementation (2009-2017) of the "Ecological Security Barrier Protection and Construction Project." MAT trend was faster, but AP trend was slower in 2000-2008, compared to those in 2009-2017 (Figures 3A,B). AR trend was similar to the GI trend in the two periods, showing an initial non-significant increasing trend and then a significant decreasing trend (Figures 3C,D).
The changes in climatic variables and grazing intensity in Tibet were non-uniform across space (Figure 4). Warming was evident across the whole plateau, especially in northern Tibet, where MAT increased by 0.6 • C ( Figure 4A). Considerable increase and decrease in AP also occurred in northern and western Tibet ( AP > 20 mm) and in central and eastern regions, respectively ( Figure 4B). In contrast, AR was dimming (negative AR) in northern and western Tibet but lighting (positive AR) in eastern areas ( Figure 4C). Regarding grazing intensity, it increased in the middle Tibet after the implementation of the "Ecological Security Barrier Protection and Construction Project" and decreased in the edge counties of the Tibetan Autonomous Region (Figure 4D).

Dynamics of Grassland Productivity Between Two Subperiods
Grassland ANPP in Tibet decreased from southeast to northwest with grassland types, from alpine meadow, alpine steppe to alpine desert steppe (Figures 5A,B). Mean ANPP for all grassland types (28.4 g C m −2 ) during 2009-2017 increased   decreased in the remaining 19% of grassland pixels, mainly in central Tibet (Figure 5C).

Effects of Climate and Grazing Intensity on Grassland ANPP Changes
The correlation coefficients between ANPP and three climatic variables (MAT, AP, and AR) and grazing intensity (GI) were presented in  Table 1).
RF model was used to quantify the relative contributions of the four responsible variables ( Mean MAT , Mean AP , Mean AR , and Mean GI ) with respect to the differences in mean ANPP ( Mean ANPP ) between the two subperiods. Mean AP was the primary driver of Mean ANPP with% IncMSE of 29.0%, followed by Mean MAT (%IncMSE = 19.9%), Mean AR (%IncMSE = 16.2%) and Mean GI (%IncMSE = 14.5%), respectively (Figure 8).

DISCUSSION
Quantifying the key ecosystems' dynamics controlling factors is crucial for ecological management and adaptation (Li L. et al., 2018;Wei et al., 2020). Alpine grasslands in Tibet are one of the most vulnerable biomes to human activities and climate change in the world. However, it is still a question of debate how and at what extent different drivers contribute to grassland productivity change on this plateau. In this study, we first used the CASA model to quantify the alpine grassland ANPP in Tibet from 2000 to 2017. Then, we investigated the changes in grassland ANPP before (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008) and after (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) the implementation of the "Ecological Security Barrier Protection and Construction Project." Finally, an RF model was used to quantify the relative importance of temperature, precipitation, radiation, and grazing intensity to the dynamics of grassland ANPP.

Increased ANPP Due to More Favorable Climate
We found that ANPP after the implementation of the project increased to approximately 81% of alpine grasslands in Tibet. The mean values of grassland ANPP increased from 27.0 g C m −2 before to 28.4 g C m −2 after the  project. This finding was consistent with Wang et al. (2017) and Huang et al. (2018) who also examined that vegetation coverage and forage supply capacity in Tibet increased slightly after the start of the project. Correlation analysis revealed that this increase in ANPP might be attributed to more favorable climatic conditions and the decrease in the grazing intensity after the project. Further, the climate in Tibet became warmer and wetter after the project's start, which positively affected the grassland ANPP. We also found that the increased precipitation had a higher effect on the f grassland ANPP than temperature. Nonetheless, this finding is not in agreement with Zhu et al. (2016) who concluded that alpine plants are sensitive to temperature as the average growing season temperature is lower than the optimum air temperature of vegetation productivity in global alpine regions . However, the finding is in line with long-term in situ monitoring, manipulative experiments, satellite remote sensing, and model simulations on this plateau. For example, Wu et al. (2016) found that growing season precipitation explained the most variance of alpine grassland species richness and aboveground biomass across the northern Tibetan Plateau. Fu et al. (2018) reported that increased precipitation has more strong effects on plant production in alpine meadows than experimental warming in Tibet. A potential explanation is that colder habitats have long shaped alpine vegetation's functional traits to adapt to the low temperatures and large diurnal temperature ranges (Elmendorf et al., 2012;Shi et al., 2014). This might be due to the reason that arctic and alpine plants have excellent resistance to shortterm temperature fluctuations (Theurillat and Guisan, 2001;Elmendorf et al., 2012).

Increased ANPP Due to Weakening Grazing Intensity
Human activities were significantly intense across the Tibetan Plateau before the implementation of the project. For example, the population reached upto 3.12 million in the Tibetan Autonomous Region of China, in 2013, which is twice as in 1965; meanwhile, the livestock number also increased from 9.74 million in 1958 to 23 million in the early twenty-first century (Fan et al., 2015). This study found that grazing intensity was negatively correlated with ANPP, suggesting that grazing activities can mediate the interannual ANPP change. After the start of the project, grazing intensity decreased significantly (Figure 3D), owing to the increase in the potential grassland productivity and the reduction in the livestock (Supplementary Figure S1). This significant decrease in the grazing intensity had a positive effect on grassland ANPP. Grazing and human-induced land use/cover change are the two most significant human disturbances for alpine grassland ANPP in Tibet, and the former is supposed to be the dominant one (Arthur et al., 2008;Harris, 2010;Chen et al., 2014). The complicated process of grazing, such as forage selection, herbivores and their trampling (Parsons and Dumont, 2003;Paruelo et al., 2008), directly and indirectly, can modify grassland productivity (Chen et al., 2007;Xiong et al., 2014;Wang and Wesche, 2016). Meanwhile, grassland productivity is highly dependent on grazing and its intensity. The prevailing view is that moderate grazing promotes plant growth due to compensatory growth while overgrazing reduces vegetation productivity (de Mazancourt et al., 1998;Schuman et al., 1999;Luo et al., 2012). At the end of the twentieth century, the stocking rates of livestock in Tibet were high in most of its counties (Fan et al., 2015). Overgrazing affected the native species diversity and ecosystem stability in Tibet and altered the structure and function of grassland ecosystems, induced C losses and grassland degradation (Zhou et al., 2017;Zhang et al., 2019). The Chinese government continuously formulates policies to limit grazing throughout the project, such as grazing exclusion, conversion of grazing land to grass and grazing withdrawal programs. This decline in grazing intensity after the start of the project elevated grassland productivity. Thus, we concluded that the check on the livestock numbers was a useful tool for the restoration of grassland in Tibet.

Climate Variability Dominated the ANPP Dynamic
The RF model revealed that the grazing intensity is far less critical than climatic variables in controlling grassland ANPP. This signifies that climate change is the primary driver for the sustainability of alpine grassland ecosystems on the Tibetan Plateau Lehnert et al., 2016;Xu et al., 2016). Although the intensity of human activities on the Tibetan Plateau is increasing rapidly, however, the impact of human activities on ecosystems is less than other regions of the world (Venter et al., 2016;. Grazing activities have not altered the dominant role of climate change for grassland variation in Tibet. However, we must highlight that, as suggested by previous studies, the biochemical cycle and their feedbacks to climate change would get more complicated under human disturbance (Kröel-Dulay et al., 2015;Peters et al., 2019). For example, in our study, we found that the relationships between climatic variables and grassland ANPP were more robust before the project's start. Those changes in correlation coefficients might cause by the changes in grazing intensity. Thus, the mechanisms of how plants respond to the coupled effects of climatic and anthropogenic stresses should be explored in the future.

CONCLUSION
This study used 9 years of datasets before and after the implementation of Ecological Security Barrier Protection and Construction Project, including meteorological conditions, grazing intensity, to explore the combined effects of climate change and grazing activities on the dynamics of grassland ANPP in Tibet. We found an increase in grassland ANPP after the project's start. Precipitation was the dominant factor in controlling the observed alpine grassland changes during the study period. Furthermore, the weak grazing intensity after the project promoted grassland productivity. Thus, we suggested that the check on the livestock numbers has a positive effect on the restoration of degraded grasslands in Tibet.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
JW and ML conceptualized this study and led the writing. ML collected and analyzed the data. YF, BN, YH, and XZ interpreted the results and revised the text. All authors contributed to this work and approved the final manuscript before submission.

ACKNOWLEDGMENTS
We acknowledge the support of all co-authors for their constructive and helpful comments and organization of this study.

SUPPLEMENTARY MATERIAL
The