Indirect E ﬀ ects of Grazing on Wind-Dispersed Elm Seeds in Sparse Woodlands of Northern China

: Grazing leads to the reduction of biomass and plays a critical role in land degradation in arid and semiarid lands. However, the indirect e ﬀ ects of grazing on the ecosystem, e.g., the e ﬀ ect on seed dispersal, have not been well understood. In this study, we built an agent-based model (ABM) to simulate how grazing intensity a ﬀ ects the seed dispersal of elm trees, one of the native vegetation species of temperate woodlands in semiarid lands. The simulated results from the ABM and observed data from the real world were compared to assess the accuracy and validity of the ABM. The results show that elm seed densities in non-grazing, light, moderate, and heavy grazing lands were 74.97 ± 1.44, 57.63 ± 0.89, 37.73 ± 0.95, and 0.97 ± 0.05 seeds m − 2 , respectively—an apparently decreasing trend. Moreover, as grazing intensity increased, the values of nugget, sill, and partial sill decreased and the value of the ratio of nugget to sill increased. This study indicates that the grazing indirectly leads to the reduction of elm seed density and the increase of spatial heterogeneity of elm seed on the ground in sparse elm woodlands. Moreover, values of geostatistical indices from the ABM were not signiﬁcantly di ﬀ erent from ﬁeld observation data except for the ratio of nugget to sill. It shows that ABMs can reasonably replicate the spatial pattern of elm seed densities in the ﬁeld and thus are useful for simulating long-distance seed dispersal in sandy lands. This ﬁnding suggests that the indirect e ﬀ ects of grazing should be considered to e ﬀ ectively protect sparse elm woodlands.


Introduction
Grazing influences land use in arid and semiarid lands in various ways, ranging from changing vegetation structure to inducing land degradation [1,2]. Grazing influences vegetation by modifying individual growth, varying species abundance, and changing population dynamics [3]. Moreover, grazing regulates the vegetation spatial pattern by changing the spatial heterogeneity of soil nutrients and influencing both intra-and interspecies interaction [4][5][6][7][8][9]. However, the effects of grazing on seed dispersal have rarely been considered. It is vital to obtain a thorough understanding of the effects of grazing on seed dispersal. This could help us protect the ecosystem and combat land degradation in arid and semiarid lands.
Seed dispersal is a vital process linking the reproductive cycle of plants with the establishment of their offspring. Thus, it has various implications for species population and ecosystem biology, especially in arid and semiarid lands where winds primarily drive seeds to shape the pattern of seed dispersal and mainly determine the spatial structure of plant populations [10,11]. At the same time, the seed dispersal process is highly stochastic, and the underlying mechanism is notoriously challenging to quantify. In the current literature, field experiments and modeling approaches have been conducted to explore the influencing factors of seed dispersal [12,13].
The seed dispersal process might be influenced by grazing in direct and indirect ways. Directly, grazing promotes more seed production compared with exclosure lands [14,15]. Indirectly, the grazing animals eat and trample vegetation. Furthermore, vegetation cover regulates seed survival on the ground [16]. Thus, we hypothesize that grazing intensities regulate seed density patterns by changing vegetation coverage.
The effects of grazing on seed dispersal have not been well understood. This might be due to two reasons. One is that the effects of grazing on seed dispersal are difficult to test with field experiments. Firstly, the trace of seed dispersal, especially long-distance, is not easy to follow [12]. Secondly, heterogeneous patches of vegetation cover could cause nonuniform samples in designed experiments [17,18]. Additionally, the effects of grazing on seed dispersal in arid and semiarid lands are commonly intertwined with the wind, the main factor regulating seed dispersal [19][20][21]. This makes disentangling various driving factors extremely difficult.
Given the challenges and obstacles in field experiments, spatial simulating models have gained increasing attention in ecological modeling for their capability to incorporate spatial and temporal aspects of environmental factors and the dynamics of plants and animals [22]; in this case, the effects of grazing and wind. Agent-based models (ABMs), also known as individual-based models in the ecological modeling field, are spatially explicit models that take into account aspects of agents and their interactions, as well as heterogeneous environmental factors [23][24][25][26]. While ABMs are mostly used in simulating animal behaviors and patterns [27][28][29], recently, ABMs have been increasingly used for plant ecology; for example, in simulating species distribution [30], tree regeneration patterns [31], plant extensions [32], interspecies competition [33], and seeds dispersal [24,34]. However, the combination of ABMs and grazing has rarely been reported, especially for arid and semiarid regions, where grazing is the main driving force regulating vegetation and land-use type.
In this research, we built an ABM simulating seed dispersal that is associated with vegetation cover, grazing, and wind. In this study, the main aim is to explore the indirect effects of grazing on seed dispersal. Knowing the indirect effects of grazing on seed dispersal is vital for predicting the spatial pattern of seed distribution and evaluating the probability of recruitment in lands. The specific objectives are (1) to test our hypothesis, provided here, and (2) to verify the feasibility of the ABM on sparse woodlands. For that, firstly, we built an ABM to simulate seed dispersal under various wind frequencies and directions. Secondly, we compared the simulated results with the experimental results observed in the field. Thirdly, we used the ABM to evaluate the effects of grazing under various scenarios that consider the vegetation cover as heterogeneous environments.

Study Area
Horqin Sandy Land is located in a typical agriculture and pasture interlaced zone on the eastern edge of the Mongolian Plateau [35] (Figure 1). Horqin Sandy Land, with an area of 5.23 × 10 4 km 2 , is one of the largest sandy lands in China [36]. Due to the sandy and loose soil, sparse vegetation cover, and its semiarid climate condition, this area is highly prone to wind erosion and desertification. Exacerbated by ever-intensifying human activities such as livestock grazing and agricultural expansion and climate changes (e.g., the decreasing precipitation) in the past years, the risk of desertification and land degradation remains very high [37,38]. The continuous land degradation of Horqin Sandy Land attracts much attention from researchers and policymakers. A series of national ecological restoration projects, such as the Grain for Green Project and the Beijing and Tianjin Sandstorm Source Treatment Project, was implemented to fight for land degradation and preserve the fragile ecosystems [35]. January being the coldest month, averaging −14.0 °C. The mean annual precipitation is 340 mm, of which 70% falls in June, July, and August. The mean annual wind velocity is 4.4 m s −1 , and the number of gale days (>16 m s −1 ) ranges between 21 and 80. The windy season is from March to May, and the growing season begins in late April and ends in late September. Typical plants are Caragana microphylla Lam., Setaria viridis (L.) Beauv., Bassia dasyphylla, Chenopodium glaucum Linn., Chenopodium aristatum Linn., Lespedeza daurica (Laxm.) Schindl., and Pennisetum centrasiaticum Tzvel [6]. Sparse woodland is the native vegetation type and climax community in Horqin Sandy Land [14]. It plays an essential role in reducing wind erosion, promoting ecological stability, and supplying ecological services [14,39]. Sparse woodland constitutes three layers, i.e., trees, shrubs, and herbs. Additionally, Ulmus pumila L. (elm trees) is the key and defining species in sparse woodland. Elm trees provide various ecosystem services and play critical roles in helping to maintain the stability of the ecosystem of sparse woodland. For example, elm trees provide a seed shadow for regeneration and shade for grasses and livestock [15,16].
Sparse woodlands and elm trees, in particular, are, however, threatened by ever-increasing grazing pressure. Grazing leads to a decrease in the area of sparse woodland [14]. In addition, elm trees are also adversely influenced by grazing in the aspects of seedling growth, reproduction allocation, and spatial distribution [6]. More interestingly, grazing also changes vegetation coverage and subsequently influences elm seed dispersal, which plays a vital role in the renewal of elm trees in the wind season [40].

The Agent-Based Model
The ABM here was developed using NetLogo, a free software platform for building ABMs [41,42]. The model description largely follows the ODD protocol, i.e., overview, design concepts, and details [43]. More details can be found in Supplementary Material. The purpose of the model is (1) to ). This area is characterized by a continental semiarid monsoon climate. The average annual temperature is 6.3 • C, with July being the warmest month, averaging 23.0 • C, and January being the coldest month, averaging −14.0 • C. The mean annual precipitation is 340 mm, of which 70% falls in June, July, and August. The mean annual wind velocity is 4.4 m s −1 , and the number of gale days (>16 m s −1 ) ranges between 21 and 80. The windy season is from March to May, and the growing season begins in late April and ends in late September. Typical plants are Caragana microphylla Lam., Setaria viridis (L.) Beauv., Bassia dasyphylla, Chenopodium glaucum Linn., Chenopodium aristatum Linn., Lespedeza daurica (Laxm.) Schindl., and Pennisetum centrasiaticum Tzvel [6].
Sparse woodland is the native vegetation type and climax community in Horqin Sandy Land [14]. It plays an essential role in reducing wind erosion, promoting ecological stability, and supplying ecological services [14,39]. Sparse woodland constitutes three layers, i.e., trees, shrubs, and herbs. Additionally, Ulmus pumila L. (elm trees) is the key and defining species in sparse woodland. Elm trees provide various ecosystem services and play critical roles in helping to maintain the stability of the ecosystem of sparse woodland. For example, elm trees provide a seed shadow for regeneration and shade for grasses and livestock [15,16].
Sparse woodlands and elm trees, in particular, are, however, threatened by ever-increasing grazing pressure. Grazing leads to a decrease in the area of sparse woodland [14]. In addition, elm trees are also adversely influenced by grazing in the aspects of seedling growth, reproduction allocation, and spatial distribution [6]. More interestingly, grazing also changes vegetation coverage and subsequently influences elm seed dispersal, which plays a vital role in the renewal of elm trees in the wind season [40].

The Agent-Based Model
The ABM here was developed using NetLogo, a free software platform for building ABMs [41,42]. The model description largely follows the ODD protocol, i.e., overview, design concepts, and details [43]. More details can be found in Supplementary Material. The purpose of the model is (1) to simulate the seed dispersal of elm trees with wind directions in Horqin Sandy Land and (2) to evaluate the effects of grazing on seed dispersal under scenarios of grazing intensities, i.e., light, moderate, and heavy. The model framework is illustrated in Figure 2. The variogram is defined by the following parameters: nugget (C0), partial sill (C), sill (C0 + C), and range (a). These parameters perform differently in the spherical model, the Gaussian model, and the exponential model. These models are determined by the equations above. We used the geoR package in R to perform the geostatistical analysis [53]. The difference of mean in the spatial indices in this model and the field investigation data was determined with a t-test. The difference between grazing intensities was investigated with ANOVA and Turkey's HSD test. P-values less than 0.05 were considered significant in this study.

Model Validation
According to the results of the field investigation, the values of nugget, partial sill, sill, and the ratio of nugget to sill in the plots were 127.93 ± 118. 18 The difference of mean in nugget, partial sill, and sill was not significant between the results of the field investigation and this model (p > 0.05, Table 1). The difference of mean in the ratio of nugget to sill was significant between the results of the field investigation and this model (p < 0.05; Table 1).  We built the model with two entities, seeds, and vegetation cover. The seeds were modeled as agents, and their movements were driven by winds. Wind directions controlled the direction of the movement of seeds. The seed locations, quantified in a 2D coordinate system, i.e., x-and ycoordinates, changed depending on the wind, with eight directions. We recorded the locations of seeds at each moment and calculated the distance from the original point (x = 0, y = 0). The vegetation cover was modeled as patches within a regular bidimensional lattice. Vegetation coverage, defined as the percentage of the ground surface covered by herbs and shrubs, was assigned as an attribute of each patch, a grid cell. The relationship between vegetation coverage and seed density followed an experimental equation in the same study area [16], i.e., b = −0.131a 2 + 9.527a − 30.53 (R 2 = 0.446). Here, b indicates seed densities, and a indicates vegetation coverage.
The model was initialized with several parameters, including seed number, frequency of the eight wind directions, and vegetation coverage. The value of these parameters, except vegetation coverage, was taken from a previous study [40]. The value of vegetation coverage was randomly designed and varied within a range of 0-30%. We conducted the sensitivity analysis on the total number of patches and running time steps of the model, which might influence the output of the model.
The effects of grazing on seed dispersal were focused on changes in vegetation cover. Vegetation cover has a significant relationship with seed densities [16]. Meanwhile, grazing influences vegetation cover by various degrees, according to grazing intensities [44]. Thus, we could evaluate the effects of grazing on seed dispersal through changes in vegetation cover. The changes in vegetation cover reflect the changing grazing intensities under various scenarios.

Sensitivity Analysis
A sensitivity analysis is used to determine which parameters have more significant impacts on a model and whether the model has stability [45,46]. The high and low values of parameter sensitivity indices indicate whether a parameter change has a great or little effect on the model, respectively. The changes in the sensitivity index could be the criteria to set the optimal values of particular parameters. The sensitivity index (Sy i ) is expressed as ( [47,48]) In Equation (1), X t and Y t are the input parameters and output of the model at time t; dX t and dY t are the changes in X t and Y t . The output of this model is the seed density. The average sensitivity index (S) is calculated as Equation (2), where n is the number of model repetitions.
A small value of the sensitivity index of a particular parameter indicates that a change in the values of this parameter has little effect on the model. Conversely, the relatively large value of the sensitivity index of a particular parameter indicates that a change in the values of this parameter has significant effects on the model. We select the relatively small values of the sensitivity index. Here, the running time steps are set to 40.

The Field Investigation
At the end of the seed rain period, we randomly selected four plots (10 × 10 m) in enclosed sandy grasslands in our study area. For these plots, the dominant species was Chloris virgata Swartz, Setaria viridis (L.) Beauv., and Corispermum hyssopifolium L. The vegetation coverage in the four plots was 10%, 30%, 10%, and less than 5%, respectively (the vegetation coverage was relatively low because of decreased precipitation). In each plot, we counted the number of elm seeds in 100 subplots (1 × 1 m) separately. According to seed densities in the subplots, the geostatistical indices of the spatial patterns of seed dispersal were estimated.

Model Validation
How successfully the model-simulated results replicate the system behaviors of the real world, i.e., observations on the ground, validates the ABM; this is critical to assessing the credibility of an ABM [49]. We used an ABM to model the seed dispersal process by simulating the movement of individual seeds. However, the ultimate goal of the model is to understand the macro patterns of the targeted system, in this case, the spatial distribution patterns of the elm seeds. Thus, we used the pattern-oriented modeling approach (POM) to validate the simulated outputs of our ABM [50,51]. Given the abstract nature of ecological modeling, with a loss of information during the modeling process due to simplification and abstraction, POM focuses on a few relevant patterns in the real system. Patterns are macro-level characteristics of a complex system and often emerged from the essential underlying processes. Thus, POM is widely used to validate an ABM by testing how well an ABM can replicate system-level patterns. In other words, if an ABM can reasonably reproduce structurally realistic patterns, it is more likely that the model catches the key dynamics in the modeled system. As there were many potential patterns for the seed dispersal process (e.g., the maximum or average distance a seed travels, the shape of seed distribution of a single elm tree, density distribution of seeds), we here focused on the spatial patterns of seed dispersal and chose geostatistical indices as our pattern.
The results of the spatial pattern of seed densities were calculated with geostatistical methods and presented with geostatistical indices, including nugget, partial sill, sill, and the ratio of nugget to sill Land 2020, 9, 490 6 of 11 (see the data analysis section). Meanwhile, we obtained these geostatistical indices from the results of the field investigation, and we compared the geostatistical indices from the field investigation to the simulated results. The values of the geostatistical indices from the simulations fell into the interval (mean ± 1 standard error) of corresponding values from the investigated results, indicating that the model's results are consistent with the data observed on the ground.

Scenario Analysis
To further explore the effects of grazing on elm seed dispersal under intensity scenarios, we simulated the effects of grazing on seed dispersal with different intensities, including heavy grazing (HG), moderate grazing (MG), and light grazing (LG). The loss of vegetation cover caused by HG, MG, and LG was 85.90%, 45.90%, and 23.70%, respectively, when HG, MG, and LG was designed as six sheep ha −1 , four sheep ha −1 , and two sheep ha −1 in Horqin Sandy Land (fine wool sheep, 50 kg on average, 2-year-old) [44]. The values of loss in vegetation caused by grazing intensities were set as parameters in this model under different scenarios. The different scenarios of grazing intensities on seed dispersal could be simulated with the ABM.

Data Analysis
Geostatistical indices, including nugget, partial sill, sill, and the ratio of nugget to sill, were calculated with geostatistical models (Figure 2). In geostatistical models, one core concept is the variogram. A variogram characterizes the spatial continuity or roughness of a data set [52]. For each variable (Z), the variogram (γ) can be calculated from all sample pairs with lag distance h between locations (s) as The variogram is defined by the following parameters: nugget (C0), partial sill (C), sill (C0 + C), and range (a). These parameters perform differently in the spherical model, the Gaussian model, and the exponential model. These models are determined by the equations above. We used the geoR package in R to perform the geostatistical analysis [53]. The difference of mean in the spatial indices in this model and the field investigation data was determined with a t-test. The difference between grazing intensities was investigated with ANOVA and Turkey's HSD test. p-values less than 0.05 were considered significant in this study.

Model Validation
According to the results of the field investigation, the values of nugget, partial sill, sill, and the ratio of nugget to sill in the plots were 127.93 ± 118. 18 The difference of mean in nugget, partial sill, and sill was not significant between the results of the field investigation and this model (p > 0.05, Table 1). The difference of mean in the ratio of nugget to sill was significant between the results of the field investigation and this model (p < 0.05; Table 1).

Effects of Grazing on Seeds Dispersal
The elm seed densities in light, moderate, and heavy grazing lands were 57.63 ± 0.89, 37.73 ± 0.95, and 0.97 ± 0.05 seeds m −2 , respectively. Meanwhile, the elm seed density in lands without grazing was 74.97 ± 1.44 seeds m −2 . The difference in seed densities between grazing intensities was significant (p < 0.05).
The values of nugget, partial sill, and sill in the heavy grazing treatment were far less than in the other three treatments (Figure 3a). The values of nugget, partial sill, and sill decreased along with the increase in grazing intensities (Figure 3a). Meanwhile, the values of the ratio of nugget to sill increased along with the increase in grazing intensities (Figure 3b). The values of nugget, partial sill, and sill were significantly different between heavy grazing and other intensities (p < 0.05).

Effects of Grazing on Seeds Dispersal
The elm seed densities in light, moderate, and heavy grazing lands were 57.63 ± 0.89, 37.73 ± 0.95, and 0.97 ± 0.05 seeds m −2 , respectively. Meanwhile, the elm seed density in lands without grazing was 74.97 ± 1.44 seeds m −2 . The difference in seed densities between grazing intensities was significant (p < 0.05).
The values of nugget, partial sill, and sill in the heavy grazing treatment were far less than in the other three treatments (Figure 3a). The values of nugget, partial sill, and sill decreased along with the increase in grazing intensities (Figure 3a). Meanwhile, the values of the ratio of nugget to sill increased along with the increase in grazing intensities (Figure 3b). The values of nugget, partial sill, and sill were significantly different between heavy grazing and other intensities (p < 0.05).

Discussion
With an increase in grazing intensity, seed densities decreased from 58 to 1 seed m −2 , indicating that grazing is likely to reduce the seed density of elm trees in sandy lands. This finding supports our hypothesis that grazing intensities regulate seed densities. However, this result is not consistent with a previous study, where the seed density in grazing lands was more abundant than in exclosure lands [14]. In the previous study, only the direct effects of grazing on elm seed production were considered, as seed density was investigated in the seed rain period; however, the indirect effects of grazing on seed production were not considered. Different from the previous study, we explored the indirect effects of grazing on seed production, i.e., the effects of grazing on seed dispersal through the change in vegetation coverage. This study found that the indirect effects of grazing play a vital role in seed dispersal. Vegetation coverage influences elm seed densities because it regulates seed dispersal by wind [16]. Grazing mainly reduces vegetation cover through excessive foraging and trampling. Thus, these studies suggest the effects of grazing on the elm seed dispersal are more complex than what we thought and have an impact through direct and indirect ways.
With an increase in grazing intensity, values of nugget, sill, and partial sill decreased, but the value of the ratio of nugget to sill increased, indicating that grazing causes an increase in the spatial heterogeneity of seed densities. This finding confirms our hypothesis that grazing indirectly influences seed density patterns by regulating grass vegetation coverage on the ground. Moreover,

Discussion
With an increase in grazing intensity, seed densities decreased from 58 to 1 seed m −2 , indicating that grazing is likely to reduce the seed density of elm trees in sandy lands. This finding supports our hypothesis that grazing intensities regulate seed densities. However, this result is not consistent with a previous study, where the seed density in grazing lands was more abundant than in exclosure lands [14]. In the previous study, only the direct effects of grazing on elm seed production were considered, as seed density was investigated in the seed rain period; however, the indirect effects of grazing on seed production were not considered. Different from the previous study, we explored the indirect effects of grazing on seed production, i.e., the effects of grazing on seed dispersal through the change in vegetation coverage. This study found that the indirect effects of grazing play a vital role in seed dispersal. Vegetation coverage influences elm seed densities because it regulates seed dispersal by wind [16]. Grazing mainly reduces vegetation cover through excessive foraging and trampling. Thus, these studies suggest the effects of grazing on the elm seed dispersal are more complex than what we thought and have an impact through direct and indirect ways. With an increase in grazing intensity, values of nugget, sill, and partial sill decreased, but the value of the ratio of nugget to sill increased, indicating that grazing causes an increase in the spatial heterogeneity of seed densities. This finding confirms our hypothesis that grazing indirectly influences seed density patterns by regulating grass vegetation coverage on the ground. Moreover, this is also in line with the findings from previous studies [8,16,54]. This study emphasizes the intermediary role of vegetation cover in linking grazing and seed dispersal. Moreover, the spatial heterogeneity of seed densities changed by the grazing may intensify the spatial heterogeneity of the vegetation as seed dispersal provides seed sources for vegetation restoration [55].
The values of geostatistical indices from the ABM were not significantly different from the investigated results, except for one index (the ratio of nugget to sill). It is suggested that the model represents some of the key processes of the seed dispersal, and, consequently, the results of the simulation are consistent with the field observations. It shows that ABMs are useful for simulating the long-distance seed dispersal in sandy lands-a daunting task faced by researchers due to its highly stochastic nature. Long-distance dispersal is relatively difficult to trace in experimental studies [12]. Meanwhile, existing simulated models depend on various environmental and vegetation factors (including wind velocity, seed release height, and vegetation height) and their prior distributions, which can be difficult to parameterize [13]. Therefore, ABMs have their advantages when simulating seed dispersal as seed movements are modeled with simple rules [34].
The study also shows that ABMs can reasonably replicate the observations on the ground, that is, the spatial pattern of elm seed densities, which has profound implications on elm seed germination, seedling growth, and sapling growth. In this ABM, seed dispersal is mainly driven by winds and is influenced by vegetation cover. To keep the parsimony of the model, only a few key factors such as the direction and frequency of wind, vegetation cover, and grazing intensities are included. Factors not considered in this research, such as wind power, slope, and animal predation, may be considered in future research [56]. However, we believe the agent-based models should not include unnecessary complications, which hinder model output analysis and the complexity represented in the model [57]. We believe that the model designed in this study strikes a balance with its parsimony and its capability to replicate empirical observations.

Conclusions
The results show that a decrease of vegetation cover (mainly grass cover), as an indirect effect of grazing, has a negative effect on the elm seed dispersal; in other words, grazing causes a decrease of seed density on the ground. Alternatively, grazing leads to an increase in spatial heterogeneity in sparse elm woodlands. This finding supports our hypotheses and suggests that the indirect effects of grazing should be fully considered to protect sparse elm woodlands, which is the original vegetation type of temperate woodlands in semiarid lands in this region.
This research demonstrates that ABMs are a promising approach in studying the seed dispersal process, e.g., with their capability to fully consider plant and environmental factors at the same time. Such ABM-based seed dispersal models can help us reveal the key dynamics in the complex process of vegetation cover and land degradation in sandy lands.