Grazing Altered the Pattern of Woody Plants and Shrub Encroachment in a Temperate Savanna Ecosystem

Ulmus pumila-dominated temperate savanna is an important tree-grass complex ecosystem in the Otindag sand land, northern China. To date, few investigations have been undertaken on the spatial patterns and structure of this ecosystem and its driving factors under different grazing pressures. The objective of our study therefore is to explore whether grazing has affected the population structure/pattern of woody plants and shrub encroachment in a temperate savanna ecosystem. Results indicate that species richness and seedlings decreased with increasing grazing pressure. An increase in grazing pressure did not significantly affect adult-tree density, but it hindered the normal regeneration of U. pumila seedlings, further inducing population decline. U. pumila seedlings had a more significant aggregated distribution than juvenile or adult trees. The adult and juvenile trees had an aggregated distribution at the small scale and a random distribution at the large scale. Shrubs also showed a significant aggregated distribution. No clear effect on the spatial patterns of adult trees was observed; however, there was a noticeable effect for juveniles and seedlings under different grazing pressures. U. pumila seedlings had a positive association with their juveniles and Spiraea aquilegifolia, but a negative association with Caragana microphylla. Shrub encroachment occurred with decreasing grazing pressures. In conclusion, overgrazing led to the decline of U. pumila population, but the decrease in grazing pressure increased shrub encroachment in the temperate savanna ecosystem. Moderate grazing management may be a better way to enhance the stability of U. pumila population and reduce shrub encroachment.


Introduction
Savanna patchy ecosystems are predominantly structured by complex interactions between the dominant trees (or shrubs) and their herbaceous plants found beneath or around trees [1][2][3]. Savanna ecosystems, which are found in temperate and tropical regions in Africa, Australia, Asia, and America, account for about 12% of the world's total area [4]. In China, temperate savanna ecosystems are mainly distributed in the ecotones of arid and semi-arid areas [5]. As a typical temperate savanna ecosystem, Ulmus pumila-dominated savanna is a climate climax community with a characteristic tree-grass two-phase structure that is composed of U. pumila-dominated woody plants and herbaceous understory vegetation. These ecosystems are widely distributed in the Otindag and Horqin sand lands The water table has been deepened to 1.0-1.5 m below ground level [34]. The two sites in this study have similar climate, geography, and soil conditions.
(a) (b) Figure 1. The location (a) and climate (b) of the study area.

Data Collection
To investigate the effects of different grazing methods on the spatial distribution of U. pumila, sample plots were selected in a flat area with a relatively low topographic relief (<1 m), thus eliminating the effect of habitat heterogeneity that is caused by the undulating dunes. To examine the effect of different grazing pressures on the U. pumila-savanna ecosystem, three rangelands were selected: (i) lightly-grazed rangeland: for 18 years, this area was grazed in the winter and not disturbed in other seasons; (ii) moderately-grazed rangeland: for 20 years, this area was fenced in the summer to exclude livestock and seasonal grazing occurred in the winter and spring, while it was mowed in the autumn; and, (iii) heavily-grazed rangeland: for 25 years, this area was public rangeland that was accessed all year for grazing.
Six sample plots with a size of 100 × 100 m (numbered as Plots 1-6) were established in the different grazing pressure areas in September, 2016 ( Figure 1). Plots 1, 3, and 5 were located in the Sanggendalai Township and Plots 2, 4, and 6 were located in the Haolaihure Township. Plots 1 and 2 were located in the lightly-grazed rangeland; Plots 3 and 4 in the moderately-grazed rangeland; and, Plots 5 and 6 were in the heavily-grazed rangeland ( Figure 2). To ensure consistency in the study, each site had similar conditions (such as soil type, geomorphology, and precipitation). Each sample plot was further divided into 25 subplots (20 × 20 m), and all woody plants in each subplot were measured and identified; trees with a diameter at breast height (DBH) ≥1 cm were marked. The species names, DBH, total height, clear-bole height, and crown diameter of all woody plants (trees and shrubs) were recorded. All elm seedlings (DBH < 1 cm) in each plot were recorded, including the heights and coordinates. The coordinates of woody plants in each plot were recorded using an Electronic Total Station with the southwestern corner of the plot as the origin. Additionally, in the center of each subplot, a small plot with a size of 1 m × 1 m was selected where all herbaceous vegetation was identified; the area of cover and number of species were also recorded. After recording this information, all grasses in aboveground were harvested to determine their biomass. The spatial distribution of woody plants is shown in Figure 2 and Table 1 details the plot characteristics.

Data Collection
To investigate the effects of different grazing methods on the spatial distribution of U. pumila, sample plots were selected in a flat area with a relatively low topographic relief (<1 m), thus eliminating the effect of habitat heterogeneity that is caused by the undulating dunes. To examine the effect of different grazing pressures on the U. pumila-savanna ecosystem, three rangelands were selected: (i) lightly-grazed rangeland: for 18 years, this area was grazed in the winter and not disturbed in other seasons; (ii) moderately-grazed rangeland: for 20 years, this area was fenced in the summer to exclude livestock and seasonal grazing occurred in the winter and spring, while it was mowed in the autumn; and, (iii) heavily-grazed rangeland: for 25 years, this area was public rangeland that was accessed all year for grazing.
Six sample plots with a size of 100 × 100 m (numbered as Plots 1-6) were established in the different grazing pressure areas in September, 2016 ( Figure 1). Plots 1, 3, and 5 were located in the Sanggendalai Township and Plots 2, 4, and 6 were located in the Haolaihure Township. Plots 1 and 2 were located in the lightly-grazed rangeland; Plots 3 and 4 in the moderately-grazed rangeland; and, Plots 5 and 6 were in the heavily-grazed rangeland ( Figure 2). To ensure consistency in the study, each site had similar conditions (such as soil type, geomorphology, and precipitation). Each sample plot was further divided into 25 subplots (20 × 20 m), and all woody plants in each subplot were measured and identified; trees with a diameter at breast height (DBH) ≥1 cm were marked. The species names, DBH, total height, clear-bole height, and crown diameter of all woody plants (trees and shrubs) were recorded. All elm seedlings (DBH < 1 cm) in each plot were recorded, including the heights and coordinates. The coordinates of woody plants in each plot were recorded using an Electronic Total Station with the southwestern corner of the plot as the origin. Additionally, in the center of each subplot, a small plot with a size of 1 m × 1 m was selected where all herbaceous vegetation was identified; the area of cover and number of species were also recorded. After recording this information, all grasses in aboveground were harvested to determine their biomass. The spatial distribution of woody plants is shown in Figure 2 and Table 1 details the plot characteristics.

Stand Structure
When grazing pressures increased from lightly-grazed to heavily-grazed conditions, herbaceous species respectively decreased from 20 and 21 species (Plots 1 and 2) to 14 and 15 species (Plot 5 and 6), and aboveground herbaceous biomass declined from 127.63 and 225.85 g·m -2 (Plots 1 and 2) to 80.84 and 103.02 g·m −2 (Plots 5 and 6) ( Table 1). Three shrub species with a large number of individuals (C. microphylla, S. aquilegifolia, and S. linearitipularis) were observed in the lightly-grazed plots (Plots 1 and 2), while only C. microphylla and S. aquilegifolia were observed in the moderately-grazed plots (Plots 3 and 4). These shrub species were not observed in the heavily-grazed plots (Plots 5 and 6).
For U. pumila, the number of seedlings significantly decreased from 106 and 101 seedlings/ha (Plots 1 and 2, respectively) to six and nine seedlings/ha (Plots 5 and 6, respectively) from the lightlygrazed plots to the heavily-grazed plots. In the moderately-grazed plots, U. pumila seedlings accounted for 84.0% and 66.4% (Plots 3 and 4, respectively) as compared to the lightly-grazed plots. Tree densities also declined from the lightly-grazed plots (214 and 187 trees/ha, Plots 1 and 2, respectively) to the heavily-grazed plots (100 and 95 trees/ha, Plots 5 and 6, respectively). However, the average DBH was highest in the heavily-grazed plots and lowest in the lightly-grazed plots, which may be due to fewer small trees being present in the heavily-grazed plots. For adult U. pumila trees, there was not a significant change in tree density between the different plots.
According to the polynomial fitting result (Figure 3a-d, polynomial order = 4), U. pumila trees had a bimodal DBH distribution under lightly-and moderately-grazed pressures (Plots 1-4). Seedlings (without DBH) and adult trees (DBH ≥ 20 cm) were more abundant in the moderatelygrazed plots (Plots 3 and 4) than juvenile trees (DBH < 20 cm), with a mixed unimodal distribution and reversed J shape. In the lightly-grazed plots, U. pumila were mainly distributed at smaller DBH

Population Structure
In a specific environment, the response of the age class of woody plants to its environment is consistent with the response of the DBH to it environment [36], meaning that the DBH of woody plants can reflect the age class of the population. Therefore, the diameter classes in our study were used to explain the age structures of the main tree species on the basis of the age dependence of tree diameter. According to the correlation between the diameter and age of U. pumila in the Otindag sand land [37], all of the trees were divided into 10 categories according to their DBH: Class I (without DBH), Class II (DBH < 5cm), Class III (5 cm ≤ DBH < 10 cm), Class IV (10 cm ≤ DBH < 15 cm), Class V (15 cm ≤ DBH < 20 cm), Class VI (20 cm ≤ DBH < 25 cm), Class VII (25 cm ≤ DBH < 30 cm), Class VIII (30 cm ≤ DBH < 35 cm), Class IX (35 cm ≤ DBH < 40 cm), and Class X (DBH ≥ 40 cm). The number and percentage of U. pumila individuals in each DBH class were calculated and mapped in each plot. Gaussian and polynomial models were used for data testing and curve fitting to infer population structure in the future.
To analyze the effect of grazing on herbaceous plants, the number of herbaceous species was counted from 25 samples in each plot. The average biomass and standard error in each plot was calculated from 25 observed herbaceous samples.

Point Pattern Analysis
Spatial pattern analysis categorized U. pumila into three classes by size: seedlings (without DBH), juvenile trees (DBH < 20 cm), and adult trees (DBH ≥ 20cm). If the number of individuals in a category was very small, this category was excluded from spatial pattern analysis.
Point pattern analysis, as proposed by Ripley [38] and improved by Diggle [39], is an effective tool in ecological studies. Ripley's K-function, expressed as the expected number of points in a circle of radius r centered at an arbitrary point divided by the intensity λ of the pattern, has been widely used in plant analysis to characterize the spatial patterns of species [40,41]. However, as the function is a cumulative measure incorporating information from both large and small scales, the function can confound large-scale effects with small-scale effects [40,42,43].
The pair correlation function g(r), which arises if the circles of Ripley's K-function are replaced by rings, provides the expected number of points at distance r from an arbitrary point, divided by the intensity λ of the pattern. The function g(r) is a probability density function for the interpretation of neighborhood density and it is sensitive to small-scale effects that can eliminate the cumulative effect that is caused by the K(r) function. Similar to Ripley's K function, the function g(r) includes both univariate and bivariate statistics. The univariate statistic is used to analyze the spatial pattern of one object, while the bivariate statistic is used to analyze the spatial association of two objects (pattern 1 and pattern 2) [40]. In this study, the univariate g(r) was used to analyze the spatial pattern of woody plants and the bivariate g 12 (r) was used to analyze the spatial associations between the U. pumila trees and shrubs, as well as between the different growth stages (i.e., seedling, juvenile, and adult stage) of U. pumila trees.
In our study, the impacts of habitat heterogeneity on point pattern analysis, including soil and topography, were eliminated by locating all plots on flat areas. Therefore, we used the null model of complete spatial randomness (CSR) as a null hypothesis to check the level of significance. In the univariate analysis, if g(r) higher than the upper confidence limit indicates aggregation, g(r) lower than the lower confidence limit indicates regularity [40]. In the bivariate analyses, when considering the fact that large trees may impact the distribution pattern of small trees within their area of influence (competition), we used a bivariate g function for large trees and small trees or trees and shrubs using both toroidal shift and the antecedent condition null model options [40]. This tests whether one pattern (small trees) is influenced by a second pattern (large trees) and it assesses whether there are more (or fewer) small trees in the neighborhood of large trees than expected under a random distribution of small trees [44]. In analyzing the interaction between trees and shrubs, we examined the spatial association between the two species with an independent null model [40]. If values of g 12 (r) are above the upper (or below the lower) limit of the confidence envelope, this indicates that pattern 2 is positively (or negatively) associated with pattern 1 at a given distance r. If values of g 12 (r) are within the estimated upper and lower envelopes, this indicates that there is no interaction between patterns 1 and 2.
Point pattern analysis was performed using Programita 2014 software (Leipzig, Germany) [40]. To eliminate the influence of edge effects on the point pattern analysis, spatial scale was assigned as the half length of plots (0-50 m) with a step of 1 m. The 99% confidence envelopes for each null model were calculated by running 199 Monte-Carlo simulations. Figures were plotted using ArcGIS 10.2 (Redlands, California, USA) and Origin 8.0 (Northampton, MA, USA).

Stand Structure
When grazing pressures increased from lightly-grazed to heavily-grazed conditions, herbaceous species respectively decreased from 20 and 21 species (Plots 1 and 2) to 14 and 15 species (Plot 5 and 6), and aboveground herbaceous biomass declined from 127.63 and 225.85 g·m −2 (Plots 1 and 2) to 80.84 and 103.02 g·m −2 (Plots 5 and 6) ( Table 1). Three shrub species with a large number of individuals (C. microphylla, S. aquilegifolia, and S. linearitipularis) were observed in the lightly-grazed plots (Plots 1 and 2), while only C. microphylla and S. aquilegifolia were observed in the moderately-grazed plots (Plots 3 and 4). These shrub species were not observed in the heavily-grazed plots (Plots 5 and 6).
For U. pumila, the number of seedlings significantly decreased from 106 and 101 seedlings/ha (Plots 1 and 2, respectively) to six and nine seedlings/ha (Plots 5 and 6, respectively) from the lightly-grazed plots to the heavily-grazed plots. In the moderately-grazed plots, U. pumila seedlings accounted for 84.0% and 66.4% (Plots 3 and 4, respectively) as compared to the lightly-grazed plots. Tree densities also declined from the lightly-grazed plots (214 and 187 trees/ha, Plots 1 and 2, respectively) to the heavily-grazed plots (100 and 95 trees/ha, Plots 5 and 6, respectively). However, the average DBH was highest in the heavily-grazed plots and lowest in the lightly-grazed plots, which may be due to fewer small trees being present in the heavily-grazed plots. For adult U. pumila trees, there was not a significant change in tree density between the different plots.
According to the polynomial fitting result (Figure 3a-d, polynomial order = 4), U. pumila trees had a bimodal DBH distribution under lightly-and moderately-grazed pressures (Plots 1-4). Seedlings (without DBH) and adult trees (DBH ≥ 20 cm) were more abundant in the moderately-grazed plots (Plots 3 and 4) than juvenile trees (DBH < 20 cm), with a mixed unimodal distribution and reversed J shape. In the lightly-grazed plots, U. pumila were mainly distributed at smaller DBH classes. DBH had a reverse J distribution, suggesting compatibility with a trend of expanding population size. According to the Gaussian fitting result (Figure 3e,f, Gaussian number peak = 1), the DBH distribution of U. pumila in heavily-grazed plots had a unimodal distribution; U. pumila trees thrived in the large diameter Class VI (20 cm ≤ DBH < 25 cm) at Plot 5 and in Class VII (25 cm ≤ DBH < 30 cm) at Plot 6. However, Classes II, III, and IV in Plots 5 and 6 were not present. This result suggests a declining U. pumila population under heavy-grazing pressures. Therefore, the effects of grazing on the population structures of U. pumila were mainly reflected in the damage of juvenile trees and seedlings, with grazing limiting the normal regeneration of the U. pumila population and accelerating the process of their decline.

Spatial Pattern and Association of U. Pumila
Analysis of the univariate spatial pattern of U. pumila using the CSR null model while accounting for antecedent conditions showed that all life stages exhibited small scale aggregations. As spatial scale increased the level of aggregation decreased, and it tended to have a random distribution (

Spatial Pattern and Association of U. Pumila
Analysis of the univariate spatial pattern of U. pumila using the CSR null model while accounting for antecedent conditions showed that all life stages exhibited small scale aggregations. As spatial scale increased the level of aggregation decreased, and it tended to have a random distribution (Figure 4).

Spatial Pattern of Shrubs
Univariate spatial analysis of shrubs found that C. microphylla had significant aggregation at the 0-14 m scale in Plot 1 and at the 0-12 m scale in Plot 2. S. aquilegifolia had significant aggregation at the 0-18 m scale in Plot 1 and at the 0-22 m scale in Plot 2. These two shrubs tended to have a random spatial distribution at all other scales in the two plots ( Figure 6).

Spatial Association of Ulmus Pumila and Shrubs
Bivariate spatial analysis of U. pumila trees and shrubs using the CSR null model while accounting for antecedent conditions found that adult trees had no significant associations with C. microphylla and S. aquilegifolia at all scales in Plots 1 or 2. In Plot 1, juvenile trees had significant negative associations with C. microphylla at the 0-10 m scale and no significant associations with S. aquilegifolia at all scales. Seedlings had significant negative associations with C. microphylla at the 0-22 m scale and a positive association with S. aquilegifolia at the 0-17 m scale. In Plot 2, juvenile trees had no significant associations with C. microphylla and S. aquilegifolia at all scales, seedlings had no significant association with C. microphylla at all scales, and S. aquilegifolia had a significant positive association with seedlings at the 0-9 m scale (Figure 7).

Spatial Association of Ulmus Pumila and Shrubs
Bivariate spatial analysis of U. pumila trees and shrubs using the CSR null model while accounting for antecedent conditions found that adult trees had no significant associations with C. microphylla and S. aquilegifolia at all scales in Plots 1 or 2. In Plot 1, juvenile trees had significant negative associations with C. microphylla at the 0-10 m scale and no significant associations with S. aquilegifolia at all scales. Seedlings had significant negative associations with C. microphylla at the 0-22 m scale and a positive association with S. aquilegifolia at the 0-17 m scale. In Plot 2, juvenile trees had no significant associations with C. microphylla and S. aquilegifolia at all scales, seedlings had no significant association with C. microphylla at all scales, and S. aquilegifolia had a significant positive association with seedlings at the 0-9 m scale (Figure 7).

Effects of Grazing on Stand Structure
Takahashi et al. [45] noted that the population structure of a species reflects its regeneration processes. When compared to the spatial structure of a stand, insights into the population dynamics and future trends of a stand can be obtained [46,47]. In an arid or semi-arid area, the most important limiting factor is precipitation and species regeneration usually occurs during periods with abundant precipitation or sufficient water habitats [48]. These ecosystems are more vulnerable to human disturbances (grazing, land use/cover change, etc.) relative to other ecosystems [31,49]. Understanding what drives grazing impacts on population structure is critical for sustainability [50].

Effects of Grazing on Stand Structure
Takahashi et al. [45] noted that the population structure of a species reflects its regeneration processes. When compared to the spatial structure of a stand, insights into the population dynamics and future trends of a stand can be obtained [46,47]. In an arid or semi-arid area, the most important limiting factor is precipitation and species regeneration usually occurs during periods with abundant precipitation or sufficient water habitats [48]. These ecosystems are more vulnerable to human disturbances (grazing, land use/cover change, etc.) relative to other ecosystems [31,49]. Understanding what drives grazing impacts on population structure is critical for sustainability [50]. In this study, our results indicated that an increase in grazing pressure resulted in a decrease in the number of herbaceous and shrub species, shrub density, the number of U. pumila juvenile trees and seedlings, and grass biomass, which was not consistent with the study by Gaitán et al. [50], who found that grazing pressure negatively impacted the cover of palatable grasses and species richness but did not affect shrub cover. This indicates that grazing pressure significantly affected the stand structure of the U. pumila-dominated temperate savanna. Our results were similar to those of Li et al. [51], who also found that grass biomass and density of U. pumila seedlings significantly decreased with increased human disturbance.
A comparison of U. pumila density under different grazing pressures showed that many juveniles and seedlings were present under lightly-grazed conditions (Plots 1 and 2), and the size distributions of U. pumila under this condition had a reversed J shape (Figure 3). This result suggests the presence of favorable conditions for the establishment and survival of seedlings and for the continuous regeneration of the population. However, U. pumila populations under heavily grazed conditions (Plots 5 and 6) had a unimodal DBH distribution ( Figure 3); very few juveniles were identified, which might not be offset by the longer lifespan of the adult trees. At each grazing pressure (Plots 1-6), adult U. pumila trees showed no obvious difference in density, suggesting that these trees had very strong anti-disturbance capability to grazing. Under moderately-grazed conditions (Plots 3 and 4), U. pumila could regenerate to a certain extent but the population density may decrease in the future. The continuation of heavy-grazing may lead to the disappearance of the U. pumila-dominated temperate savanna ecosystem Thus, over-utilization (for example, over-grazing and mowing) of the U. pumila-dominated temperate savanna may result in a decline of the U. pumila population in the Otindag sand land.

Ulmus Pumila Spatial Pattern and Influence of Grazing
Information about the regeneration of new individuals, the mortality of adults, and the overall demographics of a population following natural or human disturbance can be gained from investigating the spatial patterns of a species. This information can aid in understanding the ecological phenomena hidden in the spatial patterns [7,52]. Aggregation at a certain scale is a common pattern of species distribution in nature, which can result in the species resisting severe environments due to cluster effects; this can enable a stable population to continue to reproduce [53]. In our study, seedlings had a more aggregated distribution than juveniles or adult trees, occurring at scales <8-17 m for seedlings (Plots 1-4), <2-3 m for juvenile trees (Plots 1-4), and <1-2 m for the adult trees (Plot 1-6). This result suggests a random distribution for the adult U. pumila trees. Juvenile U. pumila trees had aggregation at the smaller scales and random distributions at the larger scales.
Interestingly, canopy diameter (Table 1) was about 7-8 m for the adult trees and 3-4 m for the juvenile trees, respectively, indicating that U. pumila seedlings may be successfully germinating under the canopies of adult trees. In arid and semi-arid areas, competition among seedlings for moisture and nutrients is relatively weak, a finding that has been frequently observed and is a response to patchy rainfall events in semi-arid savannas [25]. The aggregation pattern of U. pumila seedlings reflects their resistance to wind erosion and adaptation to the wind-sand environment, which may facilitate their growth and reduce the effects of sandstorms, thus improving their survival [19,54,55]. With the growth of U. pumila trees, increasing competition for water and soil nutrients led to self-thinning of the population. Adult and juvenile trees therefore had a random pattern at large scales. However, the aggregations for the adult and juvenile trees were evident at the small scales, which may be related to the initial distribution of the seedlings and sprouting reproduction due to human destruction.
In the lightly-grazed plots where seedlings were mostly unaffected by grazing, the aggregation distribution of seedlings may have been the result of seed dispersal or environmental heterogeneity [26,56,57]. The general characteristics of U. pumila seed shape and size (small and samara-shaped) reduce falling speed and enhance their dispersal capacity by wind. However, dispersal is often spatially limited, which leads to intra-specific aggregation [58]. In our study area, the prevailing wind-direction was from the northwest in the flowering season, which caused a spatial aggregation of seedlings to the southeast of the adult or juvenile trees (Figure 2). Directed dispersal can therefore generate apparent facilitative patterns of seedlings [54]. The aggregation may also have been the result of environmental heterogeneity, which causes uneven distribution of environmental variables both spatially and/or temporally [54]. In our study, a "fertility island" could have been created around the U. pumila trees due to plant litter and animal manure [35,59], which favored the growth of seedlings [60,61]. In addition, other investigations have found that micro-topography also affects the spatial patterns of species [62,63]. However, our plots were established in flat areas with a relatively low fluctuation (<1 m) to eliminate the effects of micro-topography. Therefore, the effects of micro-topography on the patterns of seedlings can be ignored.
Seasonal grazing and mowing in the moderately-grazed plots may lead to seedlings and juveniles being destroyed by the nonselective mowing activity in the autumn or a decline in density due to livestock feeding, thus further reducing spatial aggregation. This hypothesis may explain why numerous seedlings (predominantly germinated in the current year) were present in aggregations in August. Human disturbance therefore interrupts the chain of normal regeneration of U. pumila trees, which may eventually result in the breakdown of the U. pumila population in a semi-arid savanna ecosystem.
Overgrazing of livestock in heavily-grazed plots resulted in an extremely low density of seedlings and over-trampling from livestock affected tree regeneration. Due to this, seedling patterns in the heavily-grazed plots were not analyzed in this study. A comparison of patterns of U. pumila trees under different grazing pressures showed that there were no clear effects on the patterns for the adult trees. However, the effects were obvious for juveniles and seedlings. This implies that most juveniles are produced through coppice (sprouting from old stock) and can withstand trampling.

Spatial Associations among Woody Species
The spatial association of seedlings with juveniles was mainly positive at distances up to 4 m and the neutral at other distances. This finding suggests that seedlings were more likely to survive if their distance from other juveniles was less than 4 m. However, no significant spatial association between adults and seedlings at any scale were identified, presumably including water competition or grazing pressure. The areas under tree crowns are favored areas for animals to rest as they provide shelter from the summer heat and intense ultraviolet radiation. The soil under the canopy may also be destroyed due to trampling and recumbency, and at times a "bare soil circle" may be formed under trees [64]. All of this can lead to the hindrance of seedling survival. At the same time, adults had a significant positive association with juveniles on smaller scales, but not at other scales. We assume therefore that the associations between adults and juveniles may be caused by sprouting reproduction due to human destruction or originating from the associations between juveniles and seedlings.
The results for tree and shrub associations showed that there was almost a neutral association between the U. pumila adults and juveniles, and shrub species (including C. microphylla and S. aquilegifolia) (Figure 7). This finding suggests that the distance separating shrubs and trees had no significant effect on where they grew. Seedlings of U. pumila were negative or neutrally associated with C. microphylla, and positive associated with S. aquilegifolia (Figure 7). This finding may relate to water competition between the shrubs and seedlings. The association between U. pumila seedlings and C. microphylla shrubs suggests that C. microphylla may have adverse effects on the survival of U. pumila seedlings in areas with a high density of C. microphylla, i.e., competition correlation exists. However, in areas where the density of C. microphylla was lower, a significant correlation between C. microphylla and U. pumila seedlings did not exist, suggesting that there was no obvious competition for soil resources. The positive association between U. pumila seedlings and S. aquilegifolia shows that S. aquilegifolia may facilitate the regeneration of the U. pumila population.

Influence of Grazing on Shrub Encroachment
Although shrubs are often the dominant species in arid and semi-arid areas, their encroachment into rangelands has been considered an important indicator of degradation [65][66][67], further resulting in a decrease in grass productivity and biodiversity [68]. Although it is generally accepted that overgrazing is one of the main causes of shrub encroachment [69,70], recent findings have questioned the importance of overgrazing on shrub encroachment. Many studies have found that the reduction of grazing pressure after long-term grazing may be an important cause of shrub encroachment [71][72][73][74][75], which were consistent with our study findings. In the Otindag sand land, C. microphylla and S. aquilegifolia were the main species of shrub and additional similar results were also found: shrub density increased with decreasing grazing pressure (Table 1), with shrubs predominantly occurring in lightly-grazed plots (Plots 1 and 2). However, no shrubs were found in the heavily-grazed plots (Plots 5 and 6), which may be due to the abnormal germination and the regeneration of shrubs caused by over-trampling from livestock. Overgrazing in the heavily-grazed plots may reduce the productivity of the rangeland and result in edible shrubs being consumed, thus inhibiting their subsequent growth. After removing or reducing grazing pressures, through seasonal enclosure or other protective measures, bare land formed by grazing may be initially colonized by shrubs. This may be one reason why shrub density was higher in the lightly-grazed plots than in the other plots. Additionally, mowing in the autumn can remove shrub seedlings in the moderately-grazed plots, thus preventing shrub encroachment. Results from our investigation suggest that shrubs occurring in the lightly-grazed plots may be a result of encroachment due to decreased grazing pressure, similar to findings from previous studies. This finding also supported the suggestion that a shift in grazing pressure after long-term overgrazing is a direct factor that triggers shrub encroachment [69,76].
The mechanism of shrub encroachment to a new environment is very complex, which is closely related to the properties and functions of plants and their environment (such as precipitation pattern change, climate warming, grazing, etc.) [77]. The domestic herbivores by grazing can promote the dispersal of woody plant seeds, and the selective feeding by livestock reduced the grass biomass/cover, which led to reduction in competitiveness of herbaceous plants, and increase in competitiveness of woody plants [78]. In our study, the seeds of C. microphylla and S. aquilegifolia were dispersed by the feeding-excretion way or sticking to fur. Under the overgrazing, the livestock have to eat invasive shrub seedlings with better palatability because of fewer herbaceous plants, which makes it impossible to colonize and establish their communities. However, due to stronger competitiveness in shrubs than herbaceous plants, the reduced grazing pressure makes the shrub seedlings easier to survive and establish the stable communities, which enhance the shrub encroachments. Additionally, drought and precipitation variability are also the important factors that lead to shrub encroachments by enhancing the ability of competition and settlement of shrubs [79]. In this study area, drought is increasing continuously since 1960, especially during 1999-2011 when drought was the most serious period historically. Due to its deeper roots and stronger drought-tolerant than herbaceous plants, shrubs can utilize the deeper soil water during the drought, which promotes the growth and distribution of invaded shrubs. Under the combined effect of grazing and climate, shrubs expand continuously their distributions through self-reproduction and competition with the herbaceous species [79].

Conclusions
Because vegetation in fragile habitats is sensitive to environmental change and prone to degradation, U. pumila-dominated savanna in the Otindag sand land has been seriously degraded due to climate change. In addition, an increase in rangeland grazing pressures in recent decades has led to serious land degradation, and even to desertification of areas in this region [51]. In the heavily-grazed rangeland, the U. pumila population has declined and grasslands have become degraded. In the moderately-grazed rangeland, the U. pumila population was relatively stable but it suffered from poor seedling regeneration. In the lightly-grazed rangeland, the U. pumila population was increasing and it had good seedling regeneration, however serious shrub encroachment was occurring. The combination of our results indicates that controlled grazing is urgently needed in the Otindag sand land area to promote ecological restoration of the degraded U. pumila-dominated savanna ecosystem. Our findings show that light grazing is the optimum method for rangeland use and protecting the savanna ecosystem from further degradation. However, a shift in grazing pressures to a light-grazing situation after long-term overgrazing may reduce the availability of rangeland; therefore, we propose that a moderate-grazing regime is required. Reformation of property rights of the commons (heavily-grazed rangeland) is needed to reduce rangeland pressures and to avoid further degradation [80,81]. This investigation provides important information for management practices and planning, and specific recommendations for the U. pumila-dominated savanna. It is important to avoid planting woody shrub species that have negative interspecific associations (C. microphylla vs. U. pumila seedling) in homogeneous habitats. The presence of S. aquilegifolia and U. pumila seedlings together shows that S. aquilegifolia may be a more suitable shrub vegetation type to promote savanna restoration in this area.

Conflicts of Interest:
The authors declare no conflict of interest.