Light-intensity grazing improves alpine meadow productivity and adaption to climate change on the Tibetan Plateau

To explore grazing effects on carbon fluxes in alpine meadow ecosystems, we used a paired eddy-covariance (EC) system to measure carbon fluxes in adjacent fenced (FM) and grazed (GM) meadows on the Tibetan plateau. Gross primary productivity (GPP) and ecosystem respiration (Re) were greater at GM than FM for the first two years of fencing. In the third year, the productivity at FM increased to a level similar to the GM site. The higher productivity at GM was mainly caused by its higher photosynthetic capacity. Grazing exclusion did not increase carbon sequestration capacity for this alpine grassland system. The higher optimal photosynthetic temperature and the weakened ecosystem response to climatic factors at GM may help to facilitate the adaption of alpine meadow ecosystems to changing climate.

rapid plant growth 14,15 . In contrast, slowed nitrogen dynamics under long-term grazing can negatively influence plant growth and soil microbial activities 16 . All these changes have the potential to alter the magnitude of ecosystem carbon uptake and emission 4,17 .
The Tibetan Plateau ecosystem is a relatively fragile system and is sensitive to climate changes 18 . It acts as a critical "first response region" to climate change in China and East Asia 19 . Climatic factors and grassland management both have strong influences on the seasonal and inter-annual dynamics of carbon fluxes 3,7 .Grazing mediates the relationships between ecosystem function and carbon flux variability by means of altered canopy structure and plant physiology 4 . Changes in ecosystem function can dominate the inter-annual variability of carbon fluxes 20 and determine the ecosystem adaptability to climate change 4 . Thus, quantifying grazing effects on carbon dioxide uptake and emissions of alpine grassland systems is critical for understanding their role in global climate change.
To effectively monitor grazing effects on carbon flux, the eddy covariance technique is a reliable source of in situ measurements with its capability of automatic and high-frequency measurements of gas and energy fluxes. There have been few prior studies that have used the eddy covariance technique on paired treatments of grazing-nongrazing 12,21 , especially for actively managed grasslands 22,23 . These research gaps constrain our capacity to adequately address grazing disturbance effects. In particular, knowledge about grazing impacts on carbon fluxes is necessary to scale up carbon budget from measured sites to larger terrestrial entities 24 . To this end, in this study we focused on carbon fluxes of alpine meadow ecosystem under grazed and ungrazed treatments in the same environmental setting to examine how carbon fluxes respond to grazing in alpine meadows on the Tibetan Plateau. Our main research questions were: 1) Does light-intensity grazing could maintain higher ecosystem productivity compared with short-term fencing in these alpine meadows? and 2) How does grazing affect the responses of carbon fluxes to climatic variations?

Data and Methods
Site description. Our study site is located in the core zone of the Northern Tibetan Plateau, where a typical alpine Kobresia pygmaea meadow is distributed. The region belongs to the plateau subfrigid monsoon climate, with no absolute frost-free period throughout the year. The soil freezing period is from October to next May, and the annual mean air temperature is − 1.9 °C. The mean annual rainfall is 380 mm, with 80% falling between June and September. The annual average sunlight exceeds 2886 h. The soil is classified as meadow soil with sandy loam. The vegetation is dominated by Kobresia pygmaea, accompanied by Potentilla bifurca, Potentilla saundersiana, Leontopodium pusillum, and Carex moorcroftii.
The grazing exclusion treatment meadow (FM, 31.643686°N, 92.009683°E, 4598 m a.s.l.) has been fenced since October 2011 and lies in the same watershed with the grazed treatment meadow (GM, 31.648722°N, 92.007208°E, 4607 m a.s.l.), which ensures that the two treatment sites experience similar local climates. The topographic map around the paired eddy covariance towers is shown in Fig. 1. Prior to the grazing exclusion treatment, vegetation and climates were quite similar between fenced and grazed sites. They exhibited no significant net ecosystem exchange (NEE) differences (P = 0.997) and had relatively homogeneous vegetation, dominated by Kobresia pygmaea (> 65%). Before the experiment was established, both GM and FM was grazed using a common local grazing intensity of < 0.5 sheep and < 1.2 yak unit ha −1 half-y −1 for over 20 years, and the same grazing intensity maintained during the experiment at GM. Under the atmospheric neutral stability assumption, the footprint analysis 25 showed that approximately 95% and 97% of the measured scalar fluxes originated from the FM and GM towers, respectively.

Experimental measurements.
The open-path eddy covariance (OPEC) system was used to measure carbon and water vapor fluxes at 2.3 m above the ground. The OPEC system consists of a 3-D sonic anemometer (Model CSAT-3, Campbell Scientific Inc., Logan, UT, USA) to measure three-dimensional wind speed and temperature fluctuations, and an infrared gas analyzer (Model LI-7500A, Li-cor Inc., Lincoln, NE, USA) to measure CO 2 and water vapor densities. All signals were sampled at a frequency of 10 Hz. The CO 2 and H 2 O fluxes were calculated and recorded at 30 min intervals by a CR3000 datalogger (Model CR3000, Campbell Scientific).
The meteorological conditions, including solar radiation, net radiation (Rn) and photosynthetically active radiation (PAR), were observed at 1.5 m height above the ground using a four-component net radiometer (Model CNR-1, Kipp&Zonen, Netherlands) and a quantum sensor (LI190SB, Li-cor Inc.). Air temperature (Ta) and relative humidity (RH) were measured at 1.8 m height (Model HMP45C, Vaisala Inc., Helsinki, Finland). Soil temperatures (Ts) and soil water contents (SWC) were recorded at four depths (0.05, 0.1, 0.2 and 0.5 m) with thermometers (109-L, Campbell Scientific) and TDR probes (Model CS616-L, Campbell Scientific), respectively. The above meteorological data were summarized to half-hour intervals, except precipitation (PPT), which was measured continuously by the tipping bucket rain gauge (TE525MM-L, Campbell Scientific). All meteorological data were recorded using a CR1000 datalogger (Model CR1000, Campbell Scientific).
Replicate samples (n = 10) for aboveground biomass were collected during June-August from 2012 to 2014 once a month by clipping vegetation of 0.5 × 0.5 m 2 quadrats within a radius of 250 m around each observation tower.
Scientific RepoRts | 5:15949 | DOi: 10.1038/srep15949 Data processing. The analyses were based on half-hour means of CO 2 fluxes from 2012 to 2014. The flux data were first aligned with the coordinate system of mean wind direction using the three-dimensional rotation 26,27 , then corrected for bias caused by air density variation due to heat and water vapor fluctuations using the Webb, Pearman and Leuning density correction (WPL), a correction for the effects of air density fluctuations 28 . Anomalous or spurious values caused by interference from rain, extreme cloud cover, dew, hoarfrost, or birds were excluded from the analysis.
The CO 2 flux measured by the eddy covariance technique represents the net ecosystem CO 2 exchange (NEE), which is equivalent to the net ecosystem productivity (NEP) multiplied by − 1. The NEE was partitioned into gross primary productivity (GPP) and ecosystem respiration (Re) 29 . The daytime NEE was assumed to be valid for our analysis when photosynthetically active radiation (PAR) was greater than 1 μ mol m −2 s −1 . In calculating nighttime NEE, PAR was always less than 1 μ mol m −2 s −1 .To avoid possible flux underestimation under stable conditions during the night, effects of friction velocity u* were examined statistically. This analysis determined that the threshold values of u* were 0.16 m s −1 for fenced and 0.17 m s −1 for grazed meadows in 2012, 0.11 m s −1 for fenced and 0.13 m s −1 for grazed meadows in 2013 and 0.11 m s −1 for fenced and 0.14 m s −1 for grazed meadows in 2014, respectively. The negative nighttime NEE were excluded from our analysis. All the missing data for each 30 min were filled using the methods of linear fitting, nonlinear fitting and mean diurnal variations 27,30 . The number of missing observations, filtered values and the total observations were counted ( Table 1). The energy balance closure was evaluated based on half-hourly observations using the energy balance ratio (EBR).  Where LE represents the latent heat flux, H represents sensible heat flux, Rn represents net radiation, G represents soil heat flux, and S represents canopy heat storage. The closer EBR is to 1, the better the energy balance closure.
To investigate grazing affects, we used data from 2012 to 2013 when the grazing effect was strongest (The differences between FM and GM were significant, P < 0.01). We used the Michaelis-Menten equation 31 to calculate the ecosystem light response: where the unit of NEE daytime and PAR are mg CO 2 m −2 s −1 and μ mol m −2 s −1 , respectively; α (mgCO 2 μ molPhoton −1 ) is the apparent quantum yield or the initial slope of the light response curve; P max (mgCO 2 m −2 s −1 ) is the maximum apparent photosynthetic capacity of the canopy; and Re day (mgCO 2 m −2 s −1 ) is the ecosystem respiration during daytime.

Statistical analysis.
To determine the driving factors on carbon fluxes, stepwise multiple linear regression analysis was employed. First, we conducted a comprehensive analysis that used Grazing (value = 0 at FM; value = 1 at GM), Rn, Ta, Ts, SWC and VPD as independent variables and NEP, Re and GPP as separate dependent variable. Second, to explore the relative contribution of each factor to carbon fluxes under the grazing or non-grazing treatment, we conducted partial correlation analysis of carbon fluxes (NEP, Re and GPP) against climatic factors (Rn, Ta, Ts, SWC and VPD) for FM and GM sites separately. By comparing the two response patterns (with and without grazing being included), we can simultaneously obtain knowledge about how grazing affects carbon flux response to climate change.
The one-way analysis of variance (ANOVA) was used to examine differences of carbon fluxes and environmental variables between the fenced and grazed sites. We used OriginPro 8 SR0 software to perform all the statistical analysis and the level of significance was set at 0.05.

Characteristics of environmental variables and biotic environment. The net radiation (Rn),
photosynthetically active radiation (PAR), air temperature (Ta), soil temperature (Ts), and vapor pressure deficit (VPD) all exhibited unimodal dynamics at both FM and GM within a year. The annual total radiations were similar between FM and GM ( Fig. 2(a,b)). The average daily Ta was higher at GM than that at FM (Fig. 2(c), P = 0.048). The daily mean Ts was higher at FM than that at GM in growing season (P = 0.024), while the daily mean Ts was higher at GM in non-growing season ( Fig. 2(d), P < 0.01). Both Ta and Ts reached their peaks from early July through late August at the two sites. The VPD was always lower at FM than that at GM across the three years ( Fig. 2(e), P < 0.01). The mean SWC was lower at FM than at GM (P < 0.01) and the seasonal variations of SWC showed strong correspondence to rain events throughout each growing season ( Fig. 2(f), P < 0.01). A severe drought in August 2013 resulted in a continuous decrease of SWC from early August until September (Fig. 2(f)) and a VPD peak in August (Fig. 2(e)).
Both fenced and grazed grasslands achieved their maximum aboveground biomass in August 2012 and 2014 ( Fig. 3(a,c)). The maximum aboveground biomass was 139 and 110 g m −2 for FM and GM in 2012, and 188 and 160 g m −2 in 2014, respectively. In 2013, an earlier precipitation peak arrived in July and resulted in an earlier biomass peak ( Fig. 2(f)), with a maximum biomass values of 213 and 110 g m −2 for FM and GM, respectively. The biomass decreased sharply in August due to evident water stress in the alpine meadows.  (Fig. 3). The NEP and GPP were greater before noon than in the afternoon, especially in August 2013-a severe drought month, leading to an asymmetrical distribution of NEP and GPP around noon. The daily maximum of NEP and GPP occurred at 11:00, with values of 0.152 and 0.229 mgCO 2 m −2 s −1 , respectively. The daily maximum Re occurred around 16:00-17:00, and the values were higher at GM than at FM throughout each growing season. Both ecosystems showed net carbon gains in July and August and neutral values in June and September. However, the net carbon sink turned to a weak carbon source in August 2013, with a NEP value of − 1.64 gC m −2 in that month.
During the growing seasons from 2012 to 2014 (Fig. 5) Generally, carbon fluxes at GM were greater than those at FM in 2012. The NEP and GPP rebounded rapidly in the following two years, while recoveries of Re were not obvious. The carbon flux differences between FM and GM were obvious in 2012 and 2013 (P < 0.01), which provided an opportunity to investigate grazing effects on this alpine ecosystem.
Comparison of photosynthetic capacity between FM and GM. Solar radiations were very strong at both sites ( Fig. 2(a,b)). The NEE increased with light intensity and reached its maximum around 1500 μ mol m −2 s −1 of PAR, then leveled off (Fig. 6). The α was 0.000358 mgCO 2 μ molPhoton −1 for the FM grassland and 0.000531 mgCO 2 μ molPhoton −1 for the GM grassland, respectively. The P max was 0.21 mgCO 2 m −2 s −1 and 0.27 mgCO 2 m −2 s −1 for the FM and GM grasslands, respectively.     Under any environment group (Table 2), α and P max at the GM grasslands were greater than those at FM, indicating a greater carbon uptake capacity by the prior one. The α increased, but P max followed a flat trend along rising Ta at FM. At GM, both α and P max increased with rising Ta. The optimum temperatures were 10-13 °C and ~13 °C at FM and GM, respectively.

Figure 2. Monthly change in daily net radiation (Rn, (a)), photosynthetically active radiation (PAR, (b)), air temperature (Ta, (c)), soil temperature (Ts, (d)) at 5 cm depth, vapor pressure deficit (VPD, (e)), and soil volumetric water content (SWC, f) at 5 cm depth for the fenced and grazed alpine meadows
The value of α increased with enhanced SWC at both FM and GM. Across the two treatment grasslands, P max was greatest when SWC was in the range between 0.13 and 0.18 m 3 m −3 . The maximum α in GM and FM occurred in the VPD range between 0.3 and 0.6 and between 0.6 and 0.9, respectively, with a greater value for the prior treatment than for the latter one. Both ecosystems achieved their maximum P max when VPD fell in the range between 0.6 and 0.9 kPa, and P max of GM was 29% greater than that of FM.
Drivers of carbon fluxes. We analyzed the relationships between carbon fluxes (i.e., NEP, Re and GPP) and their potential drivers (Grazed, Rn, Ta, Ts, SWC and VPD) using multiple linear stepwise regression (Table 3). In this alpine ecosystem, grazing had the strongest influence on Re, and the Rn had the strongest influence on NEP and GPP. Grazing was the second most important driver of GPP.
The partial correlation analysis of carbon fluxes against climatic factors was conducted for the fenced and grazed meadow separately, to differentiate the role of grazing from climatic factors (Table 4). Rn was the main driver on NEP and GPP in both ecosystems. Re was mainly controlled by Ts at FM and by SWC at GM. In general, grazing weakened the correlations of GPP and NEP with climatic factors. Also, grazing weakened the correlations of Re with Rn, Ta and Ts, but enhanced the correlation of Re with SWC and VPD.    32,33 . A wider range of temporal variation has been identified at GM than at FM, though it has been fenced only for three years. This result indicates more carbon uptake in daytime and more carbon emission in nighttime at GM than at FM, which is in accord with similar studies conducted in tallgrass prairies in US 12 and desert steppe on the Mongolian Plateau 21 . The wider daily NEP amplitude at GM during growing seasons emphasizes the importance of including anthropogenic activities disturbances in future carbon flux modeling work. The grazed meadow acted as a stronger carbon sink than the fenced meadow in the first two years of fencing. This result shows that lightly-grazed alpine meadows can contribute more to atmospheric carbon uptake than similar meadows where grazing is excluded 12,21 . Similar findings have been reported in mixed prairie in Nevada, USA 34,35 . The higher GPP and Re at GM than that at FM were also in accord with previous studies 36 conducted for shortgrass steppe in northeastern Colorado, US.
The short-term fencing decreased ecosystem productivity. However, the ecosystem productivity of FM increased over time and reached a level almost equal to the GM site in the third year of fencing. Possible reasons might be that the ecosystem had adapted to the new environment 37 , and that grazing intensity in this alpine ecosystem is light. This experiment provided a way to judge the grazing intensity of a grazed grassland ecosystem. It also shows that in overgrazed areas, fencing would be a good method to facilitate ecosystem recovery.
Response of photosynthetic capacity to grazing in various environmental conditions. Ecosystem gross photosynthesis is driven by biotic and abiotic factors, as well as their interactions 30,30,38 . The higher ecosystem carbon sequestration potential at GM than at FM might be caused by a temporary negative response of the latter ecosystem to non-grazing. More suitable microclimate at GM might be another reason 21 . Temperature is a critical environmental factor for plant growth on the Tibetan Plateau 30 . The alpine plants have been acclimated to the long-term low temperature conditions 39 by developing photosynthetic and other related organisms optimal for low temperatures 30 . The higher optimum temperature for carbon uptake at GM than that at FM indicates that light-intensity grazing might improve alpine meadow ecosystem productivity under global warming.  Carbon assimilation was significantly reduced in both ecosystems when Ta < 7 °C. Both α and P max increased with Ta at GM, while P max almost leveled off along increased temperature at FM. This is the direct result of greater carbon uptake at GM than FM during the active growing season in the first two years of fencing. Grazing also indirectly affected carbon uptake by changing Ts 4 . Reduced biomass of litter and dead standing plants caused by grazing can weaken the greenhouse effect of canopy, and maintain suitable temperature and light conditions for plant growth.
Generally, water deficiency causes decreased net carbon uptake and internal leaf CO 2 concentrations through adjusting leaf stomata 40 . In this study, α increased with SWC in both treatments, which was consistent with other studies 41,42 . The GM had higher SWC than FM due to reduced evapotranspiration caused by grazing 13 . Therefore, the GM had a greater potential for carbon sequestration. The maximum P max appeared when SWC fell in the rang within 0.13-0.18, not at the point of maximum SWC, which indicates that moisture was not a limiting factor in this alpine meadow. The higher α and P max for GM than for FM might be related to the grazed meadow's more physiologically active leaves and higher single-leaf photosynthetic capacity 21 under the situation of reduced green net primary productivity (GNPP) and leaf area index (LAI) caused by grazing 14 .
Normally, there is a reduction in carbon uptake at higher VPD because of stomatal closure under drought conditions 43 . In this study, the lowest carbon uptake was not observed when VPD was greater than 0.9. On the contrary, both ecosystems had lower carbon sequestration under when VPD was lower than 0.3 because of low temperature. This indicated that Ta, not VPD was the major factor regulating ecosystem carbon uptake of this alpine meadow during active growing seasons.
Grazing effects on the alpine meadow ecosystem productivity. Grazing can increase ecosystem photosynthesis and respiration. Light-intensity grazing will stimulate development of new leaves, whose higher photosynthetic capacity can compensate for reduced photosynthesis caused by thinned leaves 12 . In addition, grazing can improve soil fertility by increasing available nitrogen, which generates another avenue for enhancing ecosystem productivity 44 and ecosystem respiration. The reason is that respiration is limited by supplies of carbohydrates fixed through photosynthesis 30 . Animal waste decomposition releases a large amount of CO 2 , which exerts a significant influence on ecosystem respiration. The improved microclimates at GM 45,46 and influences of animal saliva on grass growth 47 may also be responsible for increased carbon fluxes. For example, an open canopy allows for more adequate sunshine at GM, consequently prolonging carbon uptake time during thedaytime 21 . In contrast, the denser shading at FM can reduce carbon fluxes 48,49 .
The stepwise multiple linear regression analysis revealed the importance of grazing in controlling the carbon budget of this alpine meadow within three years since fencing. The absolute response magnitude of carbon fluxes to climate change should increase with plant biomass 4 . By reducing aboveground biomass, grazing lessened the strength of flux-climate correlations. It means that grazing weakened the response of carbon fluxes to climatic factors 3,4 , and light-intensity grazing might thus improve adaptation of this alpine ecosystem to changed climates 50 .