Vegetation Response to Goats Grazing Intensity in Semiarid Hilly Grassland of the Loess Plateau, Lanzhou, China

: Quantitatively estimating the grazing intensity (GI) effects on vegetation in semiarid hilly grassland of the Loess Plateau can help to develop safe utilization levels for natural grasslands, which is a necessity of maintaining livestock production and sustainable development of grasslands. Normalized difference vegetation index (NDVI), ﬁeld vegetation data, and 181 days (one goat per day) of GPS tracking were combined to quantify the spatial pattern of GI, and its effects on the vegetation community structure. The spatial distribution of GI was uneven, with a mean value of 0.50 goats/ha, and 95% of the study area had less than 1.30 goats/ha. The areas with utilization rates of rangeland (July) lower than 45% and 20% made up about 95% and 60% of the study area, respectively. Grazing signiﬁcantly reduced monthly aboveground biomass, but the grazing effects on plant growth rate were complex across the different plant growth stages. Grazing impaired plant growth in general, but the intermediate GI appeared to facilitate plant growth rate at the end of the growing seasons. Grazing had minimal relationship with vegetation community structure characteristics, though Importance Value of forbs increased with increasing GI. Flexibility in the number of goats and conservatively deﬁning utilization rate, according to the inter-annual variation of utilization biomass, would be beneﬁcial to achieve ecologically healthy and economically sustainable GI.


Introduction
Managed grazing covers more than 25% of the land surface of the world [1], and approximately 40% of China [2], making up the largest form of land use. Grazing has resulted in widespread degradation or destruction of lands and wetlands [3][4][5]. The semiarid hilly grassland in the Loess Plateau has been grazed for over a thousand years [6], and grazing activities are the main type of land use. The "Forage Livestock Balance" policy, balancing the relationship between forage productivity and grazing capacity of grasslands, has been implemented in this region since the year 2000 [7]. Quantitatively estimating the grazing effect on grasslands can help the government make a policy for proper grazing intensity (GI) of natural grasslands [4]. This is required to maintain livestock production and sustainable development of grasslands.
Calculating and assessing GI are basic to study grazing effects and make adequate management decisions in natural grasslands [31,32]. Most studies about the effect of grazing on vegetation are based on grazed and ungrazed controlled or fenced experiments [12,16,21,33]. Controlled experiments usually lack the long evolutionary history of plant-livestock interactions, which changes the physiognomy of the grasslands [34] and leads to the grassland community adapted to and tolerant to grazing [18]. GI is difficult to obtain by direct observation on a landscape scale, because the spatial distribution of GI is uneven, affected by forage resource, topography, and thermal conditions [35][36][37][38][39]. GPS tracking of free-grazing animals [31,32], and remote sensing of vegetation [40], can help us to obtain multi-scale spatiotemporal data, and determine the spatial interaction between vegetation and grazing animals, and the influence long-term grazing on natural grasslands.
In this study, GI calculated from 181 days (one goat per day) of GPS tracking, aboveground biomass (AGB), based on time-series normalized difference vegetation index (NDVI), and field vegetation data were combined to quantitatively evaluate the spatial pattern of GI, the effects of grazing on the dynamic vegetation production, and community structure in semiarid grassland of the Loess Plateau. The novelty of this study is in using an uncontrolled experiment to assess the impacts of GI on vegetation in a location that has been grazed for over a thousand years [6]. The objectives of this study were to (1) describe and assess the spatial distribution of GI in the study area, and explore its spatial relationship with AGB and plant growth rate; (2) analyze the response of vegetation community structure to continuous GI level; and (3) explore the proper grazing management strategies in the semiarid hilly rangeland of the western Loess Plateau. We hypothesized that (1) the spatial distribution of GI in the mountain is uneven, and relatively high GIs are distributed in the areas that must be passed by goats every day; (2) grazing reduces AGB, but intermediate GI will facilitate plant growth rate; and (3) grazing has minimal relationship with vegetation community structure characteristics over the time period of study.

Study Area
This study was conducted on the semiarid hilly rangeland of the western Loess Plateau, near Lanzhou, China ( The study area has been grazed for over a thousand years [6]. The herd used in this study is composed of 200-300 Zhongwei goats. Every shepherd has their own specific pasture. The shepherd drives his herd of goats into the fixed pasture in the morning and picks them up in the evening. The herd of goats are free-ranging with rare interference from the shepherd and other herds in the mountains, where no enclosure is set. No predation events have been reported in recent years. Goats are primarily grazed in the study area, and only kids are provided with supplemental feed in dens.

GI Obtained from GPS Tracking of Goats
Each day, one GPS unit (HOLUX GPSport 245 weighting less than 0.16% goats' body weight, with the accuracy of GPS about 10 m) was fitted to the stomach of a randomly selected adult goat in the morning at the den, and unloaded when the goat came back to the den in the evening. From 2013 to 2017, 181 days of GPS tracking data of goats from one herd were collected ( Table 1). The battery for the GPS units can support up to 24 h of continuous data collection and was recharged every day (one fix every second). The study area was approximately 714.81 ha (Figure 1) based on the minimum convex polygon derived from the 181 GPS tracking trajectories. The partial study area with GI of zero was defined as ungrazed area. The ungrazed area accounted for 7.12% of the study area, and the grazed area accounted for 85.86% of the study area ( Figure 2a). GI, as defined by Holechek, et al. [42], is the cumulative effect of grazing animals on rangelands during a particular time period. The GI distribution ( Figure 2) was calculated on the basis of the 181 sets of GPS tracking data using a grid cell method conducted in ArcGIS 10.2 software [43]. The study area was divided into 30 × 30 m grid cells, and the number of GPS points in each grid cell was counted. The GI during a specific period of time was calculated as follows [32]: where GI herd,i (goats/ha) is the GI in grid cell i; n i is the counts of GPS points within grid cell i; N herd is the number of goats in the herd (300); R is the GPS recording interval (1 s); D is the number of grazing seconds in a day (8 h × 3600 s); GP is the grazing period of total days recorded (181 days); S is the area of the grid cell (0.09 ha).

Field Vegetation Sampling Data
Eight (August 2016) and twenty-four (August 2017) field vegetation sampling sites ( Figure 1) were randomly selected to characterize the vegetation community structure characteristics during the peak of biomass period. At each sampling site, three 1 × 1 m plots were randomly laid out along the ridge line (approximately 30~60 m). For each plot, the number of individuals and heights were recorded for each species. The aboveground vegetation was collected, separated by species, and oven dried at 60 • C until constant weight to determine the AGB for each species. The Importance Value (IV) of grasses and forbs within a plot was calculated as follows: IV (%) = (relative density + relative AGB + relative frequency)/3. ( Species richness was expressed as the number of species within a plot, species diversity was expressed by Shannon-Weiner index (H), species evenness was expressed by Pielou index (J), and dominance was expressed by Simpson index (I).
where P i is the relative density of species i in a plot; S is the number of species in a plot. N i is the number of individuals of species i in a plot, and N is the total number of individuals of all species in a plot.

NDVI
Landsat 8 OLI/TIRS C1 Level-1 cloud-free datasets in August 2016 and 2017, corresponding to the field vegetation sampling date and covering the study area were downloaded (https://earthexplorer.usgs.gov/, accessed in 10 April 2020). The monthly Landsat 8 cloud-free datasets from January to December 2016 were also downloaded, but the Landsat 8 datasets in 2017 was not downloaded due to clouds in January, March, June, and July. Atmospheric correction of Landsat 8 images was conducted in Environment for Visualizing Images (ENVI) software [44]. The NDVI values were calculated in ArcMap 10.2 as follows: where NIR is the reflection value of near-infrared band; RED is the reflection value of red band.

Estimated AGB Using NDVI
A 30 m circular buffer was created for each vegetation sampling site in ArcMap 10.2, due to the vegetation sampling location error and spatial resolution of GI. The unary regression equation between the mean NDVI of the 30 m circular buffer and its corresponding measured AGB was mapped by the package basicTrendline in R [45]. The summarized result of Akaike's information criterion (AIC) [46] was used to choose the linear, quadratic, power, exponential, or log fitted model. Outliers were deleted by the standardized residuals of the predicted value greater than absolute value of 2. Based on the regression model, the monthly spatial distribution of AGB within the study area was calculated using NDVI in ArcMap 10.2.

Response of AGB and Growth Rate to Grazing
The spatial distribution of plant growth rate was calculated in terms of the spatial distribution of AGB in ArcMap 10.2 as follows: where AGB i and AGB j are the spatial distributions of AGB in months i and j, respectively, in the growing season from March to August; d is the time interval of days. A total of 40,000 random points were created in the study area. The GI, AGB in each month, and plant growth rate were extracted to the random points in ArcMap 10.2. The number of random points was 35,034 in the grazed area and 2801 in the ungrazed area. Least Significant Difference multiple comparison was conducted in R to test the differences in monthly AGB using a significance level of p = 0.05. We sorted the random points from smallest to largest in terms of GI (larger than 0), averaged AGB, and averaged growth rate every 0.05 unit increment of GI. The unary regression functions were applied to analyze the relationships between GI and AGB in each month, GI and growth rate at different plant growth stages in the package basicTrendline in R, respectively. AIC was used to choose the best models of linear, quadratic, power, exponential, or log in R. The random points in the grazed area were randomly subsampled and reduced to 2801 points. AGB and growth rate of the grazed area were compared with that of ungrazed area using t.test in R, respectively.
The effects of grazing on grassland, interannual variability of utilization of biomass, and utilization rate [47] by goats were calculated as follows: where AGB i,j , utilization biomass i,j and utilization rate i,j are the AGB, utilization biomass, and utilization rate by goats in month i and GI j, respectively.

Response of Vegetation Community Structure to Grazing
We analyzed the unary relationship between the GI of 30 m circular buffer of 32 vegetation sampling sites and AGB, density, height, richness, diversity, evenness, dominance, IV of grasses and forbs in the package basicTrendline in R, respectively. AIC was used to choose the best models of linear, quadratic, power, exponential, or log in R.

Grazing Intensity (GI)
The mean GI was 0.50 ± 0.44 goats/ha (means + std), ranging from 0 goats/ha to 3 goats/ha (Figure 2a). The frequency distribution of GI showed a decreasing trend ( Figure 2b); 95% of the study area had GI less than 1.30 goats/ha and 99% less than 2.00 goats/ha (Figure 2b). The spatial distribution of GI in the mountain was uneven, and relatively heavy GI was mainly distributed in the entrance pass, exit pass, and valley ( Figure 2a).

Estimated AGB Using NDVI
The relationship between AGB and NDVI had a significantly positive linear correlation for the vegetation sampling sites (R 2 = 0.35, p < 0.01) (Figure 3a). The spatial distribution of AGB in each month in 2016 was estimated using the time-series NDVI data based on their regression function. The mean AGB in each month, ranging from 11.05 to 58.92 g/m 2 , showed a trend of increasing from February to August and decreasing from August to December (Figure 3b).

Response of AGB and Growth Rate to Grazing
The relationship between GI and AGB showed a significantly negative linear relation from February to December 2016, respectively (p < 0.05) (Figure 4b-l). The relationship between monthly AGB and GI was insignificant in January (p > 0.05) (Figure 4a), maybe due to snow cover. The AGB of ungrazed areas was not significantly different to grazed area in January (p > 0.05), and higher than grazed areas from February to December (p < 0.01) (Figure 4).
The relationships between plant growth rate and GI level varied, showing unimodal curve, linear decrease, opposite unimodal curve, and logarithmic decrease at different plant growth stages ( Figure 5). During the whole growing season (March-August) and May-July, the growth rate significantly decreased with the increase in GI (p < 0.0001) (Figure 5a,d,e). In March-April and July-August, the unimodal curves best described the relationships between the growth rate and GI level (p < 0.01), and the growth rate peaked at a GI of 0.73 and 0.59 goats/ha, respectively (Figure 5b,f). In April-May, opposite unimodal curves best described the relationships between the growth rate and GI level (p < 0.0001), and the growth rate was lowest at a GI of 0.43 goats/ha (Figure 5c). The growth rate of grazed areas had no significant differences with the ungrazed areas in March-April (p < 0.01) (Figure 5b). The growth rate of grazed areas exceeded the ungrazed areas in July-August (p < 0.01) (Figure 5f). In other times, the growth rate of grazed areas were lower than that of ungrazed areas (p < 0.01) (Figure 5a,c-e).  The monthly utilization of biomass and rate changed over different months, and generally peaked at the periods of high biomass ( Figure 6). The monthly utilization biomass ranged from 1.37 to 9.73, 2.54 to 13.38, and 4.87 to 21.16 g/m 2 , with mean monthly values of 4.70, 7.01, and 11.93 g/m 2 at a GI of 0.50, 1.30, and 3.00 goats/ha, respectively ( Figure 6a). The monthly utilization rates ranged from 3.23% to 22.78%, 5.94% to 31.34%, and 11.39% to 49.56%, with mean monthly values of 11.00%, 16.41%, and 27.93%, at a GI of 0.50, 1.30, and 3.00 goats/ha, respectively (Figure 6b).

Response of Vegetation Community Structure to Grazing
The responses of vegetation AGB, density, height, species richness, diversity, evenness, and dominance to GI of 0-1.6 goats/ha were insignificant (p > 0.05) on the basis of vegetation sampling data (Figure 7a-g). There were eight species of grasses and 27 species of forbs in the study area (Supplementary File S1). The IV of species of grasses and forbs made up 22.31% and 77.69% of the total species, respectively (Supplementary File S1). The relationship between GI and the IV of grasses at each site was insignificant (p > 0.05), but the IV of forbs at each site increased with the increase in GI (p < 0.001) (Figure 7h,i). , evenness (f), dominance (g), and Importance Value (IV) of grasses (h) and forbs (i). The relationships are best described by dark blue lines. Blue shades refer to 95% confidence interval of the fitted curves. Significant regression at p < 0.05 is represented.

Relationship between AGB and NDVI
The relationship between AGB and NDVI had a significantly positive linear correlation for the vegetation sampling sites (R 2 = 0.35, p < 0.01) (Figure 3a). NDVI is consistent correlation with vegetation AGB [25,48,49] and NDVI acquired from TM data was the most important predictive factor for AGB in the Loess Plateau [49]. The explained coefficient (R 2 ) between the NDVI and AGB was 0.45 in semiarid grassland [31] and 0.57 in typical steppe region [50] of China. The explained coefficient (R 2 = 0.35) between NDVI and AGB found in this study was lower than other studies in grasslands. The explained coefficients of this study need further improvement by increasing vegetation sampling points, improving GPS accuracy, extending study region, or adding other vegetation indices. The monthly spatial distribution of AGB within the study area was calculated using NDVI, and then the utilization measurement and growth rate were calculated based on the estimated AGB.

Assessing the GI
The spatial distribution of GI in the mountain was uneven (Figure 2). The relatively high GIs were mainly distributed in the entrance pass, exit pass, and valley, and GI decreased with greater distance from these areas ( Figure 2). GI varied in accordance with the movements of livestock and the spatial distribution of available forage resource, topography, and thermal conditions [35][36][37][38][39]. High GIs are often found around the water points or surrounding of settlements, resulting in overgrazing and degeneration of grasslands [32,51]. Bare soil patches appeared around the entrance pass, exit pass, and valley (relatively flat) due to heavy grazing pressure and trampling. The grazing distributions of these goats and cattle in semiarid rangelands of North America both showed preference for flat and lowland areas [39].
In this study, the mean GI was approximately 0.50 goats/ha, ranging from 0 to 3 goats/ha, and 95% of the study area had less than 1.30 goats/ha (Figure 2). The mean GI in our study is lower than other regions with similar annual precipitation grasslands [33,52]. The carrying capacity in a mountainous area ranged from 0.63 to 0.92 goats/ha in southeastern Spain with mean annual precipitation of 324 mm [52]. A GI of 2-3 sheep/ha was sustainable in the Horqin sandy grassland in Inner Mongolia with annual mean precipitation of 370 mm [33]. However, only relying on the GI of animals to evaluate the effect on vegetation is difficult, because the effect varies with grazing animal, climatic condition, and grassland type.
The utilization rates of rangeland can be used as a management tool to achieve a GI that is ecologically healthy and economically sustainable [28][29][30]. The monthly utilization rates ranged from 3.23% to 22.78%, 5.94% to 31.34%, and 11.39% to 49.56%, with mean monthly values of 11.00%, 16.41%, and 27.93%, at GI of 0.50, 1.30, and 3.00 goats/ha, respectively (Figure 6b). Thus, the highest utilization rate (July) in at least 95% of the study area was lower than the agricultural standard of proper utilization rate of temperate typical steppe of 45% [47]. The utilization rate of approximately 20% was taken as moderate grazing pressure in the typical steppe of China, and moderate grazing pressure was used for ecosystem services and livestock production [53]. Johnston, et al. [30] proposed that the long-term safe utilization rate by domestic livestock is 15-20% in south-west Queensland, Australia. In approximately 60% of the study area, the highest utilization rate (July) was lower than 20%. Conclusively, the GI of most study areas were light or moderate in terms of the utilization rate.
The monthly utilization of biomass and rate generally peaked at the periods of high biomass (Figure 6), when conservative livestock numbers should be adjusted to maintain a safe utilization rate and grassland condition over a long-term period. Hunt [28] suggested that managing utilization is the priority during the growing seasons, when grazing has detrimental effects on plant survival, regrowth, and productivity [54]. The goats cannot obtain sufficient forage in low biomass periods (for example, the utilization biomass by goats in February was only 14.18% of that in July at a GI of 0.5 goats/ha). This effect was accentuated by high GI. Supplementary forage grass or reducing the number of goats at the end of the growing season can achieve the best compromise between livestock production and rangeland sustainability [28].

Response of AGB to Grazing
The monthly AGB (from February to December) exhibited a significant decreasing trend with an increase in GI from the spatial analysis of the study area (p < 0.05) (Figure 4), though an insignificant relationship was found from the vegetation sampling data (p > 0.05) (Figure 7a). Yan, et al. [16] found that light GI reduces AGB, but the effects were insignificant from a meta-analysis. AGB exhibited a significant decreasing trend with an increase in GI in most grassland regions of northern China [16,23,[55][56][57] and the world [10]. The percentage reduction of AGB in our study area (13.24% in July at a GI of 0.5 goats/ha) was lower than the proportional reduction 47.86% from other site of Loess Plateau [25,58], 42.77% from 251 data sets in China [16], and 23% from 236 data sets around the world [10].
The relationship between GI and AGB showed a significantly negative linear rela-tion ( Figure 4), indicating that the AGB did not exhibit the phenomenon of plant compensatory growth, and plant growth rates failed to compensate for the tissue lost to grazers. During most of the growing periods, the plant growth rates of grazed areas were lower than ungrazed areas ( Figure 5), indicating that grazing impairs the plant growth rate. At the beginning (March-April) and end (July-August) of growing seasons, the plant growth rates of grazed areas were higher than ungrazed areas (Figure 5b,f). During the above two periods, the plant growth rates showed unimodal curves correlation with GI, peaking at intermediate GI (Figure 5b,f). Liu, et al. [59] found that intermediate GI level can promote plant growth rate in meadow steppe. The fast transition of the grazing effect on plant growth rate and the complex and various relationships between plant growth rate and GI level at different plant growth stages ( Figure 5) were possibly related to local climatic factors [11] and species regrowth ability at different growth stages [60]. Exploring the influence factors about the grazing effect on plant growth rate can help us to explain the controversy of compensatory growth.

Response of Species Diversity and Composition to Grazing
Our results indicated that grazing had minimal relationship with vegetation community structure characteristics on the basis of field sampling data (Figure 7). The responses of plant AGB, density, height, species richness, diversity, evenness, and dominance were insignificant (p > 0.05) (Figure 7a-g). Mackey, et al. [61] concluded from a meta-analysis that nonsignificant relationships between grazing disturbance and diversity, richness, and evenness are common. From the meta-analyses of Wang, et al. [62], the plant diversity significantly increases after grazing, but the change of species richness and diversity affected by grazing is smaller from a world-wide meta-analysis [21].
Species diversity did not peak at intermediate GI ( Figure 7e). Yuan, et al. [17] believed that many studies miss the species diversity peak at intermediate disturbance, because the disturbance levels are either extremely small or extremely large. The spatial scale in this study area was less than Yuan, et al. [17]. Many other studies found that grazing has a negative impact on plant diversity and height, especially for degraded, nutrient poor, or low productivity ecosystems [20][21][22]24], where the ability to recover after grazing is low [22,27]. The long evolutionary history of plant-livestock interactions changed the physiognomy of grasslands [34] and led to the grassland community adapted to and tolerant to grazing [18], which may partially explain why no significant correlation was found between GI and vegetation diversity.
There was no significant response in the IV of high forage quality grass species in 2016 and 2017 to GI (p > 0.05) (Figure 7h). The dominant grass species of S. capillata and E. borealis exhibited relatively high grazing tolerance and regrowth ability [60,63]. Desirable species composition of grasses can be maintained at proper GI [28][29][30]. However, the IV of forb species in 2016 and 2017, relatively low forage quality, increased with the increase in GI from the vegetation sampling analysis in 2016 and 2017 (p < 0.01) (Figure 7i). Extending multiple years of analysis on the effects of grazing to species composition would be helpful to verify these results. Grazing increased the forbs component of vegetation in the temperate grassland on the Mongolian Plateau, and grazing removal increased the dominance of grasses [3]. Livestock prefer species of high forage quality [64]. Changes in vegetation composition from palatable grass species to less preferred forb species resulting from heavy grazing have been reported in many ecosystems [20,[23][24][25]. In the fire-grazing system, the preference of herbivores for recently burned areas rather than greater time, since fire resulted from herbivores selecting for high forage quality [65]. The effect of grazing on species composition was determined by livestock preference and the different strategies of plants to compete and compensate for herbage removal under grazing [60,66,67].

Conclusions
The spatial distribution of GI in the mountain was uneven, with a mean GI of 0.50 goats/ha, and 95% of the study area had less than 1.30 goats/ha. The monthly utilization of biomass and rate changed over different months, and generally peaked at the periods of high biomass. Utilization rates (July) lower than 45% and 20% made up approximately 95% and 60% of the study area, respectively. This finding indicated that most study areas belong to light or moderate GI. The spatial relationship between GI and monthly AGB showed a significantly negative linear relation. However, the grazing effect on the plant growth rate was complex, and various relationships were developed at different plant growth stages. At the end (July-August) of growing season, the plant growth rates of grazed areas were higher than ungrazed areas. During this period, the plant growth rates showed unimodal curve correlation with GI, peaking at intermediate GI. Further exploring the grazing effect on plant growth rate can help us to explain the phenomenon of compensatory growth. On the basis of field sampling vegetation data, the GI had minimal relationship with vegetation community structure characteristics, such as AGB, density and height, diversity, evenness, dominance, and IV of grasses, though IV of forbs increased with the increase in GI.
In order to achieve ecologically healthy and economically sustainable GI, utilization rate should be adjusted in the high biomass period and forage grass should be supplemented or the number of goats should be reduced at the end of the growing season. In addition, we should prohibit over exploitation of the areas of entrance pass and exit pass by periodically changing these places though human intervention. Balancing of GPS tracking data and field sampling should be addressed, and integrating vegetation dynamics with spatial movements of livestock and abiotic environmental factors would be helpful to make adequate management decisions for grasslands.  Institutional Review Board Statement: Ethical review and approval were waived for this study. All the procedures involving goats capturing and GPS installing were performed by the local shepherd under the current ethical guidelines. No drug injections are involved in this study.
Informed Consent Statement: Not applicable.