Grazing enhances species diversity in grassland communities

In grassland studies, an intermediate level of grazing often results in the highest species diversity. Although a few hypotheses have been proposed to explain this unimodal response of species diversity to grazing intensity, no convincing explanation has been provided. Here, we build a lattice model of a grassland community comprising multiple species with various levels of grazing. We analyze the relationship between grazing and plant diversity in grasslands under variable intensities of grazing pressure. The highest species diversity is observed at an intermediate grazing intensity. Grazers suppress domination by the most superior species in birth rate, resulting in the coexistence of inferior species. This unimodal grazing effect disappears with the introduction of a small amount of nongrazing natural mortality. Unimodal patterns of species diversity may be limited to the case where grazers are the principal source of natural mortality.

species richness in Inner Mongolian grassland sites 18 . Thus, explanations of why and how an increase in grazing intensity promotes the coexistence of many more species are needed.
Previous studies have modeled the coexistence of vegetation. The lottery model revealed the importance of environmental variability in promoting coexistence [19][20][21] . However, lottery models eventually lead to competitive exclusion over a very long time scale 22 . Therefore, a lottery model is not appropriate for testing the effects of grazing because of the eventual extinction of all but one species. We need a spatial component in the model to achieve the coexistence of plant species as in lattice models 22 .
Using a lattice model, Kakishima et al. examined the contribution of animal seed dispersers to the coexistence of tree species 23 . Tubay et al. applied this lattice model to build a general model of plant communities and reported the importance of spatial heterogeneity in determining the coexistence of many species 22 . These two models are able to avoid competitive exclusion as a stable state. Therefore, we can test the effects of grazing by using these lattice models. Here, we introduce grazing effects into a lattice model of a grassland community to examine how species diversity is affected by grazing intensity. We modify the lattice community model of Tubay et al. to develop a two-layer lattice model of a grassland community under grazing (with an animal layer) 22 . Using our model, we examine the effects of animal grazing on grasslands and investigate the underlying mechanism of these effects. The highest species diversity in the grasslands is observed under an intermediate grazing intensity. In contrast, species diversity decreases with an increase in grazer density. Furthermore, this unimodal grazing effect disappears with the introduction of a small amount of nongrazing natural mortality. We discuss when, why and how the highest diversity is achieved at an intermediate grazing intensity. We also discuss the implications of our results for the management of grassland communities.

Results
First, we simulate a lattice grassland model with 20 plant species under various levels of grazing intensity (G) without nongrazing natural mortality, i.e., d i = 0 ( Fig. 1, Supplementary Figure 1). The temporal dynamics exhibit large fluctuations in the number of individuals of each species over a long period of time ( Fig. 1a1,b1,c1). In contrast, this number is fairly stable over a short timescale (Fig. 1a2,b2,c2). Based on the short-term temporal dynamics (Fig. 1a2,b2,c2) and snapshots of the plant distributions ( Fig. 1a3,b3,c3), the number of surviving plant species is highest at an intermediate grazing intensity (G).
We vary the basal grazing intensity G to examine the effects of grazing (Fig. 2, Supplementary Figures 3, 4, Supplementary Table 1). The number of species is highest at an intermediate grazing intensity (G) for a given fecundity B and interspecific grazing difference G′ ( Fig. 2; Supplementary Figure 2b,c,e,f). The interspecific grazing difference G′ is the difference in animal preference or selectiveness in grazing, e.g., palatable herbs or unpalatable grass. To examine the effects of a trade-off between birth rate and grazing intensity, we consider the following opposite relations between the species-specific expected birth rate B i and grazing rate g i : (1) the trade-off case: species with a high B i have a high g i (grazing is subject to a trade-off) and (2) the opposite case: species with a high B i have a low g i . In both cases, the number of surviving species exhibits unimodal behavior in response to the grazing intensity G. The number of coexisting species at peak diversity is larger under the first case (with the trade-off) than under the second case (the opposite) (Fig. 2).
If there is no interspecific grazing difference (G′ = 0), then the number of surviving species decreases monotonically with an increase in the grazing intensity G (Supplementary Figures 2a,d, 3a, 4a,c). However, if a slight difference (e.g., G′ = 0.0005) is introduced, then the number of species responds unimodally (Fig. 2, Supplementary Figures 2b,c,e,f, 3b, 4b,d). Interestingly, the peak height (the number of surviving species) decreases as the interspecific grazing difference G′ is further increased (differences in colored lines in Fig To understand the effect of grazing on population density, we examine the composition of surviving species (Fig. 4). In the first (trade-off) case, the species with the lowest grazing rate dominates at a low grazing intensity (Fig. 1a3), and that with the highest birth rate dominates at a high grazing intensity (Fig. 1c3), while various species survive at low densities at an intermediate grazing intensity (Fig. 1b3). In contrast to the second case (opposite to the trade-off case) ( Fig. 4i-p), however, the results of the first case (with the trade-off) reveal a large fluctuation when the grazing intensity G is high (Fig. 4a-h). Here, which species survive depends on the simulation run (Fig. 4c,d,g,h). In any case, the results converge, and monotonic behavior is recovered as G′ decreases to zero (Supplementary Figure 2a,d). Accordingly, the results indicate that the species dependence of grazing G′ is important for unimodality but the trade-off between the birth rate and grazing intensity is not essential.
We also examine the effects of dispersal distance and lattice size on the number of surviving species (Supplementary Figure 5). We find no effects of dispersal distance on surviving species (Supplementary Figure 5a). Furthermore, the number of surviving species increases with the lattice size (Supplementary Figure 5b).

Discussion
The current report demonstrates that the empirically observed promotion of diversity by grazing is possible. This paper presents the first possible mechanism to explain the enhancement of species diversity by grazing. Our result confirms that the coexistence of plant species is possible under grazing in grasslands (Fig. 1). The results reveal unimodality in which the peak species diversity occurs at an intermediate grazing intensity (Fig. 2). When there is a trade-off between the birth rate and grazing rate, the number of coexisting species is large (more than 10 species) (Fig. 2a-c). Therefore, this trade-off promotes coexistence by balancing competitiveness among all plant species. We also test the case in which the grazing rate is negatively correlated with the birth rate. Surprisingly, unimodality in response to grazing intensity is still observed (Fig. 2d-f). In this case, the intermediate grazing rate allows less superior species to persist by occupying the new vacant sites produced by grazers (Fig. 4j,n), while it is not strong enough to exclude them (as in Fig. 4l,p). Thus, the coexistence of species is mechanically promoted because (weak) grazing produces vacant sites available for other species. The number and composition of coexisting species vary depending on simulation runs ( Fig. 4a-d,i-l), while on average they are determined by the birth rate and grazing mortality rate ( Fig. 4e-h,m-p). The trade-off between the birth rate and grazing rate is not necessary for the unimodality in response to grazing intensity.
The promotion of coexistence by grazing disappears when nongrazing mortality is introduced (Fig. 3d-f, Supplementary Figures 2a, 2d, 3). Therefore, an increase in diversity due to grazing does not occur when grazing is not the principal source of mortality in grasslands, i.e., when the nongrazing mortality is not negligible. This prediction should be tested empirically in future field studies. The interspecific grazing difference G′ also plays an Which species survive under grazing depends on the interspecific relationship between the birth rate and grazing mortality rate (Fig. 4). In the trade-off case, inferior species (with low birth rates) exclude superior species (with high birth rates) (Fig. 2a-c). In the opposite case, superior species dominate the grassland (Fig. 2d-f). The composition of surviving species is highly variable in the trade-off case but stable in the opposite case (Fig. 4). These results should be verified empirically in managed grasslands. Furthermore, the current simplified assumptions about grazing rate g i may be modified explicitly to include the behavior of grazing animals, e.g., the frequency dependence in food-plant selection.
The effects of herbivores on plant diversity in grasslands have been analyzed by using empirical data [24][25][26] . The proposed grazing model integrated with a microhabitat locality model 22 should be sufficient as a canonical model with which to investigate the functional mechanisms of herbivore grazing effects mathematically. In terms of maintaining and increasing diversity, our results are in line with those of some empirical studies 2,15-18,27,28 . For example, Komac et al. stated that maintaining an adequate grazing intensity (avoiding both abandonment and overgrazing) is necessary to preserve diversity in grasslands 15 . Mu et al. also found that light and moderate grazing promote plant biodiversity 16 . In addition, Chen et al. showed that heavy grazing has an indirect influence on biodiversity 18 . Furthermore, Chen et al. based on empirical observation inferred that grazing type and vegetation structure that affect spatial variation are the reasons for the high species richness in the Qinghai alpine meadow 2 . The two reasons are in line with our models and results. Moreover, by modifying the lattice setting, the current model can be applied to other terrestrial communities while considering predators in the system. However, this mechanism of grazer or predator effects should be directly verified by using empirical data in the future.
We should also note that this grazing model presents similar results with the intermediate disturbance models 13,19,20,29,30 . In these models, there may be several specific underlying mechanisms promoting coexistence. www.nature.com/scientificreports www.nature.com/scientificreports/ Qualitatively the current grazing model appears to conform to these models. However, whether the current model is considered as one such specific mechanism or not is a future issue.

Methods
Lattice Model. Following the lattice models of Tubay et al. 22 , we analyzed a modified lattice Lotka-Volterra (LV) competition model of a grassland community with grazing: i i The parameters b i , d i and g i denote the birth rate, nongrazing death rate and grazing mortality rate of plant species i, respectively. The model of Tubay et al. represents only the situation where plant species compete for space (i.e., direct sunlight and soil) 22 . In the modified model, we add another two-dimensional lattice (100×100) of grazers that randomly feed on plants. The grazer site Y, i.e., the site occupied by a population of grazers, moves randomly to feed on plants. Death by grazing (Eq. (3)) occurs only in site Y. The density I y of cells occupied by grazers (Y) and the initial density I i of plant cells (X i ) are kept constant. The stronger the grazing intensity is, the higher the g i is. One site is either vacant or occupied by a single plant species. In contrast, one grazer site may represent multiple grazing animals, depending on the grazing intensity. Nongrazing mortality is assumed negligible, i.e., d i = 0, unless otherwise stated (Fig. 3).
Following Tubay et al. 22 , we assume that the birth rate is given by  www.nature.com/scientificreports www.nature.com/scientificreports/ We consider the following two cases of species dependence of the grazing mortality g i : Trade-off case: i . The parameters G and G′ are the (basal, species-independent) grazing intensity and the species-dependent factor, respectively. In the first (second) case, the grazing intensity is strongest for the weakest (strongest) of 20 species.
We analyze the effect of grazing mortality (g i ) on the number of surviving plant species according to the following simulation procedure: 1. Initial distributions a. Individuals of plant species i are randomly distributed over the square-lattice cells with the initial density I i ( = 0.03). b. Grazer cells are also randomly distributed over the square-lattice cells with a given density I y . 2. Reaction dynamics a. Birth process (Eq. (1)): Two cells are chosen randomly. If the cells are X i and O, then cell O is changed to X i with probability b i . Otherwise, the cells remain unchanged. b. Death process (Eq. (2)): One cell is chosen randomly. If the cell is X i , then it is changed to O with probability d i . Note that this process is omitted when d i = 0. c. Grazing process (Eq. (3)): One cell is chosen randomly. If the cell is occupied by a grazer (Y) and species i (X i ), then X i is changed to O with probability g i . Otherwise, the cell remains unchanged. 3. Repetition process a. Steps 2a-2b and 2c are repeated L×L and k×L×L times, respectively, where L×L is the total number of square-lattice sites and the number of grazing episodes is fixed at k = 3. As a default, we set L = 100. This whole process is called a Monte Carlo step. 4. Process 3 is repeated for a specific number (e.g., 10,000 times) of Monte Carlo steps.

Data Availability
Simulation program is provided in Supplementary Information. All data are derived from the simulation program.