Grazing intensity enhances spatial aggregation of dominant species in a desert steppe

Abstract Understanding how grazing activity drives plant community structure or the distribution of specific species in a community remains a major challenge in community ecology. The patchiness or spatial aggregation of specific species can be quantified by analyzing their relative coordinates in the community. Using variance and geostatistical analysis methods, we examined the quantitative characteristics and spatial distribution of Stipa breviflora in a desert steppe in northern China under four different grazing intensities (no grazing, NG, light grazing, LG, moderate grazing, MG, and heavy grazing, HG) at three small spatial scales (10 × 10 cm, 20 × 20 cm, 25 × 25 cm). We found that grazing significantly increased cover, density, and proportion in standing crop of S. breviflora, but decreased height. The spatial distribution of S. breviflora was strongly dependent upon the sampling unit and grazing intensity. The patchiness of S. breviflora reduced with sampling scale, and spatial distribution of S. breviflora was mainly determined by structural factors. The intact clusters of S. breviflora were more fragmented with increasing grazing intensity and offspring clusters spread out from the center of the parent plant. These findings suggest that spatial aggregation can enhance the ability of S. breviflora to tolerate grazing and that smaller isolated clusters are beneficial to the survival of this dominant species under heavy grazing.

found that grazing had limited the increase of shrub cover and had redistributed soil carbon in the profile. Wan, Bai, Schönbach, Gierus, and Taube (2011) tested the responses of aboveground biomass to two grazing management systems across different levels of organization (i.e., species, plant functional group, community), and found that under continuous grazing or haymaking aboveground biomass production at all organizational levels was reduced, whereas annual alternation of grazing and haymaking had no pronounced effects on aboveground biomass.
Many studies have also examined the impact of grazing on population spatial patterns. Grazing may affect individual plants and plant populations through several mechanisms, such as removal of plant shoot issues, dung and urine return, and trampling (Chen, Christensen, Nan, & Hou, 2017). Grazing alters life history strategy and resource utilization responses of individuals and plant populations, alters competition for limited resources and changes adaptive strategies, resulting in different spatial distributions of plant population or communities. In general, spatial distribution and population dynamics reflect the impacts of environment on individuals' survival and population growth as well as the responses through ecological adaptive strategies of indicator plants (Dieckmann, Herben, & Law, 1997). Spatial scale is an important dimension to assess when investigating mechanisms underpinning the observed changes in a population (Bai et al., 2012), and spatial distribution is strongly dependent on spatial scale. Few studies, however, have examined the mechanisms of spatial distribution of plant individual or population responses to grazing at different small scales (Kleijn & Steinger, 2002).
Strategies to cope with grazing by herbivores differ among plant species, and spatial aggregation-in which an intact plant cluster divides into several smaller isolated clusters (Figure 1)-can enhance the tolerance of clone plants to stresses (Wang et al., 2017). There are two forms of clonal growth: phalanx and guerrilla (Lovett-Doust, 1981). Liu, Yu, Ye, and Dong (2007) suggest that Cleistogenes squarrosa, a phalanx clonal grass, has an advantage in growing ramets under stressful and low nutrient conditions. Wang et al. (2017) (Ding et al., 2016). The vegetation in this area is dominated by Stipa breviflora, Artemisia frigida, and Cleistogenes songorica, with average vegetation height of 5 cm and a canopy cover ranging from 17% to 20%.

| Experimental design
To quantitatively test the effects of grazing intensity on the desert steppe ecosystem, twelve adjacent plots (each ca. 4.4 ha) were established at a grazing experimental site in 2004. The plots were F I G U R E 1 Stipa breviflora in the desert steppe. (a) Sampling quadrat (4 m 2 ); (b) Divided into several smaller isolated clusters; (c) Intact cluster. The square indicates the sampling quadrat, and the circles indicate clusters of S. breviflora arranged in a randomized complete block design, which included four stocking rate treatments with three repeats for each stocking rate ( Figure 2). Seasonal grazing started on 1 June and ended on 31 November in each year since 2004. The initial sheep from the same cohort were 2-year-old Mongolian wethers and individuals were replaced after 3 years. The daily grazing schedule was from 6:00 a.m.

| Vegetation sampling
In August 2016, for sampling of herbaceous species, ten 50 cm × 50 cm quadrats were randomly positioned to record species name, height, density, and cover in each plot. As shown in Figure 1 Figure 1c). Then, the standing biomass was clipped at the soil surface, oven-dried at 65°C, and weighed. Four plots were selected to represent the LG, MG, HG plots in block I and the NG plot was selected in block II in order to avoid edge effects because the NG plot in block I is located at the edge of the experimental site ( Figure 2). In each of the four plots, one 2 × 2 m quadrat was selected. In each quadrat, a 1 × 1 m quadrat frame with 10 × 10 cm grids was placed four times from left to right and from top to bottom of the quadrat sequentially, and a tape measure used to identify the precise spatial locations of S. breviflora in the quadrat. The origin of the coordinates was defined as the upper left corner of the quadrat.

| Data analysis
The basic quantitative characteristics of S. breviflora in the NG, LG, MG, and HG treatments were indicated by cover (%), density (clusters/m 2 ), height (cm), and proportion of standing crop in the community (%). Since the data on cover, density, height, and proportion of standing crop in the community did not conform to the normal distribution, square root transformation was carried out for cover, density, and height, and extreme values of cover were removed, so that the cover, density, and height data conformed to normality.
Since the proportion of standing crop is limited to the range 0-1, the arcsine transformation was carried out to transform the data to a normal distribution. A generalized linear model (GLM) was used to test the effects of grazing intensity on the basic quantitative characteristics of S. breviflora in the treatments, and we used Duncan's test (Levene's test for homogeneity) to compare grand means among the grazing treatments. Then, subquadrats of 10 × 10 cm, 20 × 20 cm, and 25 × 25 cm sampling scales were defined by dividing the 2 × 2 m quadrat into 400, 100, and 64 subquadrats, respectively. 10 × 10 cm was the size of the majority of intact clusters of S. breviflora in the NG treatment (as measured by the authors), while 20 cm × 20 cm was double the minimum scale and 25 × 25 cm was defined considering sample size. The design of the spatial scale gradient can reflect whether the density of S. breviflora increased with increasing scale and whether the scale of spatial autocorrelation increased. The density of S. breviflora was calculated in each subquadrat. Two-way GLM was used to test the effects of grazing intensity and spatial scale on the density of S. breviflora in the plots and generalized linear model was used to test the effects of spatial scale on the density of S. breviflora in the plots. Variance analysis was undertaken using SAS 9.4 (SAS Institute Inc.) at the p < 0.05 level of significance.
Before geostatistical semivariogram analysis, the skewness (S), kurtosis (K), and confidence intervals of the sample data distribution were calculated. If all skewness and kurtosis were contained within the intervals, the sample data were considered to have a normal distribution.
The calculations were performed in Excel 2010 (Microsoft Inc.), and the sample data fitted within the skewness and kurtosis confidence intervals. The density of S. breviflora in each square (10 × 10 cm, 20 × 20 cm, and 25 × 25 cm) was analyzed by geostatistics. Using a kriging method for spatial interpolation, the pattern F I G U R E 2 Schematic diagram for the grazing experiment plots. Dark color plots indicate sampling plots. The grazing experiment plots (each ca. 4.4 ha) were arranged in a randomized complete block design, which included four stocking rate treatments with three repeats at each stocking rate. The stocking rates were 0, 0.91, 1.82, 2.71 sheep·hm −2 ·half year −1 representing no grazing (NG), light grazing (LG), moderate grazing (MG), and heavy grazing (HG), respectively F I G U R E 3 Semivariogram model maps of S. breviflora were graphed according to the semivariogram (Matheron, 1963): where r(h) is the semivariogram when the sample spacing is h; N(h) refers to the number of samples with interval h; and Z(x i ) and Z(x i + h) represent the measured values in the corresponding locations, respectively.
In order to quantitatively study the spatial autocorrelation of S. breviflora density, spatial interpolation optimal theoretical models (linear, exponential, Gaussian, or spherical model, Figure 3) were used for semivariance optimum fitting. Analysis was performed using GS + 9.0 (Gamma Design Software, LLC). The number of lags was 1, which stands for the 10, 20, and 25 cm interval scales, respectively, and is related to the coordinates of the response scale. The lag interval used was the default setting and was generally equal to 1/2 of the maximum sample scale. Spatial autocorrelation is used to reflect spatial heterogeneity. Thus, the stronger the spatial autocorrelation, the lower the heterogeneity, the higher the uniformity, and the higher the degree of spatial aggregation. The optimal fit model was used to determine the preliminary type of curve according to the scatter diagram (Clark, 1981), and then the principle of least squares was used to estimate the parameters of the initial type curve and to determine the optimal curve (Cressie, 1991). In this model, Nugget Variance (C 0 ), Sill (C 0 + C), Structural ratio (C/(C 0 + C)), and Range parameter (A 0 ) are important parameters (Table 1, Figure 3).
When C/(C 0 + C) is less than 25%, this represents weak spatial autocorrelation, while between 25% and 75% represents moderate spatial autocorrelation, and greater than 75% represents strong spatial autocorrelation. A (Range) indicates the maximum distance of spatial correlation. Sometimes this is called the effective range in order to distinguish the range (A) from a model's range parameter (A 0 ). In GS+, the Range (A) is calculated from A 0 as described in the formulae for the different models, where the range of spatial autocorrelation of the linear, exponential, and spherical models was A 0 , 3 A 0, and A 0 , respectively.

| The quantitative characteristics of S. breviflora
Grazing intensity significantly affected the basic quantitative characteristics of S. breviflora ( Table 2). Cover of S. breviflora was significantly lower under NG than under LG and HG, but did not differ between NG and MG or between LG and HG (

| Change in the density of S. breviflora population at different scales
Spatial scale significantly affected the density of S. breviflora (Table 2). Density increased markedly with increasing spatial scale, and it dramatically increased by approximately 300% as the spatial scale shifted from 10 × 10 cm to 20 × 20 cm, while density increased by approximately 56% from 20 × 20 cm to 25 × 25 cm ( Figure 5). The

Parameters
Abbreviation Interpreting Nugget or Nugget variance The y-intercept of the model; spatial variation caused by random factors.
Sill C 0 + C The model asymptote; The proportion of structural spatial distribution factors in the maximum spatial variation.

| The effects of grazing intensity on S. breviflora spatial heterogeneity
In the NG treatment, spherical, linear, and exponential models were best fitted at 10 × 10 cm, 20 × 20 cm, and 25 × 25 cm scales, respectively (Table 3). We found the largest spatial variation was caused In the LG treatment, the best fitted models of the semivariogram were linear and exponential models ( Table 3)  In the MG treatment, the best fitted models of the semivariogram were spherical and exponential models ( Table 3)  In the HG treatment, the best fitted models of the semivariogram were spherical and exponential models ( Table 3). The largest spatial variation caused by random factors was 0.035 at 10 × 10 cm scale (Figure 6a4), and the largest spatial variation caused by structural factors was 3.394 at 25 × 25 cm scale (Figure 6c4). The maximum spatial variation was 3.414 at 25 × 25 cm scale, and the largest structure ratio was 99.99% at 20 × 20 cm scale. Structural factors increased with increasing scale. The range of spatial autocorrelation of S. breviflora was 7.8, 26.8, and 36.75 cm at 10 × 10 cm, 20 × 20 cm, and 25 × 25 cm scale, respectively. Thus, the range of spatial autocorrelation of S. breviflora increased with spatial scale. Combined with the structural ratio of the semivariogram, we found that the range of spatial autocorrelation and structural factors increased with increasing scale, which indicates that the heterogeneity of S. breviflora decreased with spatial scale and that patterns were mainly caused by structural factors.

| The spatial distribution of S. breviflora population under grazing
The 2-d spatial pattern map directly presents the heterogeneity and complexity of S. breviflora spatial distribution and shows the patchiness, hierarchy, and the mosaic distribution of S. breviflora.
At 10

| Response of S. breviflora population to grazing intensity
Plants have different response strategies to defoliation by herbivores. We found that grazing decreased the height of S. breviflora, but increased its cover, density, and proportion in the standing crop.
These results suggest that S. breviflora has a relatively strong ability to survive grazing. With increased grazing intensity, the height of S. breviflora was reduced. However, the proportion in the standing crop was increased. This suggests that domestic animals prefer to search and select palatable species though selective foraging (Jamieson & Hodgson, 1979). In this study, we measured the effects of long-term grazing across a gradient of grazing intensity on the S. breviflora community. We found that density increased with grazing intensity. Long-term overgrazing (e.g., HG) leads to a decrease in the palatability of plant species, thus increasing livestock selective foraging time and frequency of trampling. Frequent trampling results in intact clusters dividing into several smaller clusters, thereby increasing the density of S. breviflora. The reason is that S. breviflora is a perennial dense cluster grass, and within one growing season every basic tillering node of the mother culm usually produces a tiller. Tillers can also undertake tillering that produces secondary tillers, finally forming tillering clusters in an exponential progression.
In addition, the growing point of S. breviflora is above ground, and it is therefore sensitive to livestock trampling (Branson, 1953), which promotes spatial aggregation.
Furthermore, the phalanx structure of S. breviflora forms an effective barrier for fixing and accumulating sand (Liu, Lv, Wang, Yan, & Wei, 2018), thereby forming phytogenic hillocks that further prevent land erosion and induce sand deposition (Wang, Wang, Dong, Liu, & Qian, 2006). When a cluster of S. breviflora is buried by sand, the tiller node undergoes displacement (Chen, Zhang, Wang, Zhan, & Zhao, 2001). This is another reason why intact clusters of S. breviflora break into several smaller isolated clusters.
Compensatory growth of plants may well explain the response of S. breviflora to defoliation by grazing. Compensatory growth associated with grazing may result from the stimulation of photosynthesis in remaining green tissues (Anten & Ackerly, 2001), reallocation of resources (Zhao, Chen, & Lin, 2008), and/or activation of additional meristems because of release of apical dominance (Liu, Yu, He, Chu, & Dong, 2009).

| Propagation of offspring around the center of the mother plant
The 2-d spatial pattern maps clearly show that the contour lines became closer and the value of S. breviflora individual density decreased gradually from the center to the periphery. These findings indicate the regeneration capacity and diffusivity into surrounding space around the center of the mother plant. At 10 × 10 cm scale, the clusters were clearly independent of each other under NG (shown by independent closed contour lines, and mean basal diameter of S. breviflora branching generally around 10 cm, personal observation), whereas the clusters presented larger dense patches with increasing grazing intensity. Simultaneously, the range of spatial autocorrelation of S. breviflora showed that spatial heterogeneity of S. breviflora distribution reduced with increasing scale (from 10 × 10 cm to 25 × 25 cm) under MG and HG. These results suggest that heavy grazing intensity reduced spatial heterogeneity and promoted the mother plants to spread into the surrounding area. Phalanx plants have an advantage in acquiring local resources and therefore may have a competitive advantage in a homogeneous environment with higher spatial aggregation of clusters (Saiz, Bittebiere, Benot, Jung, & Mony, 2016). The pattern of spatial aggregation may result from limited seed or clonal dispersal, environmental heterogeneity (Xue, Huang, Yu, & Bezemer, 2018) and positive herbivore-plant feedback, which can enhance the capacity of plants to colonize different microhabitats, to tolerate resource heterogeneity, to compete, and to recover from herbivory predation (Schmid, Puttick, Burgess, & Bazzaz, 1988).
Offspring surrounding the mother plant may be advantageous because older ramets can share resources to support the development of younger ramets (Herben, 2004). Studies have shown that it is more common for resources to be shared from developmentally older to developmentally younger ramets than for resources to be shared from developmentally younger to developmentally older ramets (Song et al., 2013). For example, Alpert (1996) reported that nitrogen moves mainly from the older to the younger ramets. Thus, it appears that physiological integration in many plants primarily provides support to the establishment of daughter ramets. Another advantage of offspring developing close to the mother plant is that physiological integration allows support of clone parts growing in low-resource patches (Wang et al., 2009).
This may represent a trade-off with the spatial heterogeneity of resources.

| CON CLUS ION
We conclude that S. breviflora populations generate spatial aggregation with increasing grazing intensity and that offspring clusters are spread out to the surrounding area from the center of the mother plant. Our results suggest that spatial aggregation can enhance the ability of S. breviflora to tolerate grazing and that smaller isolated clusters are beneficial to the survival of this dominant species under heavy grazing. This may further affect ecosystem stability and sustainability in these grasslands.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

AUTH O R CO NTR I B UTI O N S
Guodong Han conceived and designed the experiments. Baolong Yan wrote the main manuscript text and performed the experiments.
Shijie Lv and Zhongwu Wang processed and analyzed the data.
Sarula Kang contributed in manuscript editing. All authors reviewed the manuscript.

DATA ACCE SS I B I LIT Y
All data (cover, density, height, standing crop and relative coordinates' data) from this manuscript are publically available in the FigShare database (https://doi.org/10.6084/m9.figshare.7326536).