Heterogeneous microcommunities and ecosystem multifunctionality in seminatural grasslands under three management modes

Abstract Increasing attention has been paid to the relationship between biodiversity and ecosystem functioning (BEF) because of the rapid increase in species loss. However, over the past 20 years, most BEF studies only focused on the effect of species diversity on one or a few ecosystem functions, and only a few studies focused on ecosystem multifunctionality (i.e., the simultaneous provision of several ecosystem functions). Grassland ecosystems have important economic, environmental, and esthetic value; thus, this study focused on the heterogeneous microcommunities in grasslands under three management modes. The multifunctionality index (M‐index) was assessed at community and microcommunity scales, and the relationship between species diversity and multifunctionality was investigated. The communities were found to be respectively composed of one, three, and six microcommunities in grazing, clipping, and enclosure management, based on a two‐way indicator species analysis (TWINSPAN) and detrended correspondence analysis (DCA) for community structure. Biodiversity and soil indicators showed an apparent degradation of the grazing community, which had the worst M‐index. Clipping and enclosure communities showed no significant difference in biodiversity indices, soil variables, and M‐index; however, these indices were clearly different among microcommunities. Therefore, the microcommunity scale may be suitable to investigate the relationship between vegetation and multifunctionality in seminatural grassland ecosystems. Dominant species richness had more explanatory power for ecosystem multifunctionality than subdominant species richness, rare species richness, and the number of all species. Therefore, it is important to distinguish the role and rank of different species in the species richness–multifunctionality model; otherwise, the model might include redundant and unclear information. Communities with more codominant species whose distribution is also even might have better multifunctionality.


| INTRODUCTION
Global change is altering biodiversity worldwide at an unprecedented speed, resulting in unpredictable consequences for ecosystem functions (Cardinale et al., 2012;Grime, 1997;MacDougall, McCann, Gellner, & Turkington, 2013;Valencia et al., 2015). In particular, with the dramatic decline in biodiversity, researchers have increasingly focused on the relationships among the changed biodiversity, community structure, and ecosystem functioning (BEF) (Diaz & Cabido, 2001;Diaz et al., 2007;Gamfeldt, Hillebrand, & Jonsson, 2008;Grime, 1997;Hector & Bagchi, 2007;Tilman et al., 1997). Over the past 20 years, most studies only focused the effect of species diversity on one or a few ecosystem functions (Hector et al., 1999;Simova, Li, & Storch, 2013;Valencia et al., 2015;Waide et al., 1999), and few evaluated multiple ecosystem functions. However, ecosystems are primarily valued because they provide multiple functions and services simultaneously (i.e., multifunctionality) (Gamfeldt et al., 2008;Hector & Bagchi, 2007;Sanderson et al., 2004;Soliveres et al., 2014). Therefore, assessing how human management activity may impact multifunctionality is crucial to understand the ecological consequences of human management on seminatural grasslands (Byrnes et al., 2014;Maestre, Castillo-Monroy, Bowker, & Ochoa-Hueso, 2012;Soliveres et al., 2014;Valencia et al., 2015;Wagg, Bender, Widmer, & Heijden, 2014). Hector and Bagchi (2007) defined multiple ecosystem services or processes as ecosystem multifunctionality and quantified the effects of species diversity on multiple combined ecological processes for the first time. Gamfeldt et al. (2008) and Zavaleta, Pasari, Hulvey, and Tilman (2010) also defined and quantified ecosystem multifunctionality separately. They considered that multifunctionality was the ability to maintain multiple ecological functions or services simultaneously and that sustaining multiple ecosystem functions in grassland communities requires higher biodiversity. Maestre, Quero, et al. (2012) used 14 soil variables (reflecting ecological processes of the C, N, and P cycles) to synthetically evaluate ecosystem multifunctionality in global dryland, and the evaluation method and indicators in this study have been widely adopted in recent years by studies that assess ecosystem multifunctionality (Byrnes et al., 2014;Soliveres et al., 2014;Valencia et al., 2015;Wagg et al., 2014).
In China, the grassland area is 3.55 × 10 8 hm 2 , accounting for approximately 41.7% of the territorial area and having important economic, environmental, and esthetic value (Yan et al., 2014); however, human activities have seriously affected the plant distribution, community structure, and ecosystem functions of the natural grassland communities (Xu, Gao, & Wang, 2014). Therefore, investigating the change in the structure and ecosystem function of grassland communities is important for the sustainable management and development of grassland. Grazing, clipping, and enclosure are three main grassland management modes, and they result in different succession states for grassland communities and generate different interspecific and intraspecific relationships, causing heterogeneous community structure. Heterogeneous microcommunities (or dependent communities) are notably different in species composition and appearance from the background community they belong to and are composed of different types of plants. The formation and existence of microcommunities are largely dependent on the gramineous background community, and the microcommunities are a part of the horizontal structural differentiation of the overall community (Chen & Li, 2011;Song, 2001). Because of the differences in species composition and microhabitat conditions among different microcommunities, the change rules in multifunctionality would be more distinct at the microcommunity scale. Therefore, this study aimed to (1) analyze the heterogeneous microcommunity composition in a grassland ecosystem under three management modes according to community classification and ordination, (2) compare the change rules of multifunctionality at the microcommunity and community levels to determine the appropriate scale for researching multifunctionality, and (3) investigate the relationship between species diversity and multifunctionality to determine the relative importance of different species (dominant species, subdominant species, and rare species) for changes in multifunctionality.

| Study site
The experiment was conducted in Huihe National Nature Reserve in Hulunbeier Grassland, Inner Mongolia (118°48-119°45′E, 48°10-48°57′N), which has a total area of 3468 km 2 . The regional climate is a temperate continental monsoon climate. The annual average temperature is −2.4 to 2.2°C. The frost-free period is 100-120 days, and the annual average precipitation in 2008-2014 was 375.03 mm, of which 70% was concentrated between June and August. Zonal soil types are chernozem and chestnut (Li, Zheng, Ye, Xia, & Feng, 2014).

| Experiment design and sampling
In June 2008, we set four study sites of 1 ha per site in a local herder's pasture with uniform and smooth terrain that was used as a freegrazing pasture before 2008 ( Figure 1a). According to the most common grassland management modes in the local region, each site was randomly divided into two parts (each part had a size of 0.5 ha) to conduct clipping and enclosure management, respectively ( Figure 1b).
In the enclosure parts, no grazing or human disturbance occurred in the entire experiment period, whereas in the clipping parts, grass was cut once every year around August 20 (local grass-mowing time).
In August 2014 at the peak of biomass, we established one plot of 25 m × 37 m in the center of each clipping and enclosure part, respectively. In each plot, we set 15 quadrats of 1 m × 1 m per quadrat at equal intervals ( Figure 1c); thus, 60 clipping quadrats and 60 enclosure quadrats were recorded in total in four study sites.
Because our study sites were located in a herder's pasture (rented from a local herder) where grazing was forbidden according to local pasture utilization status, we were not able to conduct grazing treatment in the above-mentioned four study sites. However, as a reference, two grazing sites (50 m × 100 m per site) were set within the public grazing path (the road that livestock crossed to pasture) outside the pasture where the above-mentioned four study sites were set, and the grazing intensity was found to be approximately 62.3 standard sheep unit/km 2 (equivalent to the local actual free-grazing intensity).
As the grazing community showed highly a homogenized structure, we set 25 quadrats of 1 m × 1 m along the grazing path direction. In addition, the species-area curves of communities under three management modes (Appendix S5) showed that no new species emerged when the number of sampling quadrats reached 19; therefore, we only selected 20 grazing quadrats in the subsequent experiments. Then, we recorded the presence of all plant species and determined their density, coverage, mean height, and relative frequency (RFi) in each quadrat under the three management modes.

| Soil sampling and analyses
Soil cores (0-15 cm depth) were sampled using the quincunx sampling method to obtain a mixed sample for each quadrat. We thus collected 60 soil samples for both clipping and enclosure communities, and 20 samples for the grazing community. Maestre, Castillo-Monroy, et al. (2012), Maestre, Quero, et al. (2012) and Valencia et al. (2015) respectively selected nine to 14 soil variables to evaluate ecosystem multifunctionality. Here, we selected the following 12 variables to assess multifunctionality: pH, total nitrogen (TN), available nitrogen (AN), soil organic carbon (SOC), total phosphorus (TP), available phosphorus (AvP), soil moisture content (SMC), bulk density (BD), capillary moisture capacity (CMC), cation exchange capacity (CEC), capillary porosity (CP), and noncapillary porosity (NCP). These variables constitute a good proxy for processes such as nutrient cycling, biological productivity, and buildup of nutrient pools, which are important determinants of ecosystem functioning (e.g., water and soil conservation, soil respiration, and carrying of flora and fauna) in dry lands (Valencia et al., 2015). Most of these processes are also considered to support ecosystem services, as other types of ecosystem services depend on them (Isbell et al., 2011). Redundance analysis (RDA) (Figure 2) also confirmed that the selected 12 soil variables had a significant influence on the plant distribution (p = .0020, Monte Carlo test). The first two environmental axes explained 70.1% of the total information, and the first axis was mainly determined by factors such as CEC, SMC, CMC, BD, and AvP, which explained as much as 53.4% of the information content.

| Species importance values and biodiversity
The species importance values (IV) under different management modes were calculated as follows: (1) Study area and sites setting in Hui He National Nature Reserve, Inner Mongolia. Red circles represent 4 study sites in fenced-off area (a); each site was divided two parts, conducting clipping treatment and enclosure treatment, respectively. (b). In 2014, one plot of 25 m × 37 m was set in the center of each treatment, and 15 quadrats of 1 m × 1 m were set at equal interval in each plot. (c) A total of 60 clipping quadrats and 60 enclosure quadrats were recorded in four study sites. Blue circles represent two grazing sites within the public stock road outside fenced-off area where RD i , RC i , and RF i represent the relative species density, relative coverage, and relative frequency, respectively.
The biodiversity index includes the Simpson index, Margalef index, and Evenness index as described below (Song, 2001): where i is the ith species, S is the number of species, N is the total number of individuals of a species occurring within the quadrats, and P i is the relative density of the ith species.

| TWINSPAN clustering and DCA ordination
According to the quadrat data from the three management modes, a quadrat-species matrix based on the coverage of each species was first established and then classified into four levels using twoway indicator species analysis (TWINSPAN) in PC-ORD 6.0, and finally, the classification results for the quadrats and species were obtained. Furthermore, a quadrat-species matrix based on coverage was sorted and mapped using detrended correspondence analysis (DCA) in CANOCO 4.5. Combining the results of TWINSPAN clustering (Appendix S3) and DCA ordination (Appendix S2), we obtained the main microcommunities under the three management modes (Appendix S1).

| Multifunctionality index
Multifunctionality was estimated from all of the soil variables measured using the M-index of Maestre, Castillo-Monroy, et al. (2012). To obtain the M-index for each community, Z-scores were first calculated for each of the 12 soil variables estimated at each quadrat surveyed.
Raw data were normalized before calculations; a square root transformation normalized most of the variables evaluated. The Z-scores of the 12 soil variables were averaged to obtain the M-index in each quadrat. This index provides a straightforward and easily interpretable measure of the ability of different communities to sustain multiple ecosystem functions simultaneously (Byrnes et al., 2014). It is also statistically robust (Maestre, Castillo-Monroy, et al., 2012) and is being increasingly used when assessing multifunctionality (Quero et al., 2013;Bradford et al., 2014;Pendleton et al., 2014;Wagg et al., 2014). Moreover, the relatively large number of variables employed to calculate the M-index makes it relatively robust to outliers or atypical values.

| Multifunctionality index among microcommunities
Based on TWINSPAN and DCA results, we divided communities under three management modes into 10 microcommunities (Appendix S1) and selected the most typical quadrats from the 10 microcommunities mentioned above to compare M-index. The whole grazing community was treated as one microcommunity (marked as G1) containing 20 quadrats; the enclosure community included three microcommunities (marked as E1, E2, and E3), each type containing 10 quadrats, and the clipping community included six microcommunities (marked as C1, C2, C3, C4, C5, and C6), each type containing five quadrats. First, we square root-transformed the 12 soil variables and biodiversity indices to meet the required assumption of normality of the dependent variables in further statistical treatments. We used a mixed-model approach to compare the differences for each index between the clipping community and enclosure community, among which soil variables, biodiversity indices, and M-index were regarded as independent variables, management modes were the fixed factor, and the four study sites were random variables. Then, to compare

| Species abundance class
To facilitate the analysis and interpretation, each species was classified as a dominant species (relative abundance ≥10%), subdominant species (1% ≤relative abundance <10%), or rare species (relative abundance <1%) based on its relative abundance (density or coverage) in the whole community under the three management modes (Clark & Tilman, 2008). were not significantly different between the clipping and enclosure communities, and the data from the different sites also did not show significant aggregation (Table 1).

| Heterogeneous community composition and microcommunity types
No significant differences were observed in biodiversity indices between clipping and enclosure communities at the community level, but the species composition and community structure were not identical. The species rank abundance curve showed that the species number (5) of the first superiority rank in the enclosure community was obviously lower than that in the clipping community (10) (Figure 3). Biodiversity indices and soil variables were treated as dependent variables. Management mode (clipping and enclosure) was treated as a fixed-effect factor. "Site" (n = 4) was treated as random factor to address the nonindependence of quadrats in the same sites. *p < .05; **p < .001.
In the enclosure community, the dominant species, Leymus chinensis and Artemisia capillaries, were in the absolute dominant position and had an IV significantly higher than that of the other dominant species; additionally, species that had only one individual also occupied a certain proportion (Figure 3, Table 2). In the clipping community, the species IV showed little difference among different dominant species ( capillaries (Appendix S2, S3). Based on the DCA and TWINSPAN results, we divided the three communities into 10 microcommunities (Appendix S1).

| Soil variables, biodiversity, and multifunctionality based on microcommunities
There were no significant differences in soil variables, biodiversity indices, and M-index between clipping and enclosure management at the community scale, but significant differences were observed at the microcommunity scale (Figure 4, Tables 3 and 4). Among the three microcommunities in the enclosure community, E2 was dominated by codominant species, whereas E1 and E3 were dominated by monodominant species, that is, Leymus chinensis and Artemisia capillaries, respectively. Dominant species and subdominant species richness and evenness were higher in E2 than those in E1 and E3 (Tables 3 and 4), whereas the values of major soil variables including BD, AN, TN, TP, SOC, and CEC in E2 were also obviously better than those in E1 and E3; therefore, E2 had a better M-index than E1 and E3 (Tables 3 and   4). In the clipping community, the values of the soil variables CEC, TN, AN, TP, and SOC were all worst in C4; therefore, its M-index was the lowest. Because C1 had the highest content of AvP, AN, and SOC, it also had the highest M-index. Additionally, the number of dominant species in the clipping community in C1 was higher than that in C4, which had higher rare species richness (Table 3).

| Relationship between species diversity and ecosystem multifunctionality
M-index was positively correlated with all diversity indices, suggesting that increasing diversity can enhance ecosystem multifunctionality.

| Heterogeneous microcommunities under different management modes
Grazing could surpass other environmental factors and become the leading factor controlling plant community (Wang & Li, 1999 Therefore, those three plants were the dominant species in the grazing community, and the entire grazing community was almost entirely composed of these three plants (Wang & Li, 1999;Wang, Liang, Liu, & Hao, 2000), leading to a homogeneous community structure and the least heterogeneous microcommunities.
Intermediate disturbance with annual clipping increased biodiversity (Table 1), but the clipping interval, that is, a single growing season, was not long enough for sufficient interspecific competition.
Nonselective clipping also resulted in an equal probability and frequency of each species being disturbed, hindering the generation of a dominant species that could win the competition rapidly. After the grazing stopped, some increased plant populations aggregated to form patches during the restoration process in the clipping community.
Population patches were an effective organization form for adapting to interspecific competition, and each population partitioned its resources by occupying space (Wang, Liu, He, & Liang, 1996a). The population patches were unstable because one or a few species could not take full advantage of the resources in the space that they occupied; however, the process of species interinfiltration in units of population patches would last longer than species interactions in units of individual plants. Therefore, there were more microcommunity patches in the horizontal structure of the clipping community, and this "metastable T A B L E 4 The difference of biodiversity indices and M-indexes among different microcommunities in enclosure and clipping communities state" or "disturbance climax" was expected to last for a long time (Wang, Liu, He, & Liang, 1996b).
In the enclosure community, stratification was apparent; the microcommunities decreased, and many population aggregates began to fuse. Compared to the clipping community, in the enclosure community without interference, more species began to appear in population patches. Because their niche complementarity and interspecific positive interactions then enabled them to better utilize resources, the biomass increased (Cardinale et al., 2012;Hector et al., 1999). High productivity reduced the spatial heterogeneity of the limited resources, making the habitat more homogeneous, and then, interspecific competition became important. The species that had higher resource utilization efficiency (especially for light) thus competitively excluded the subdominant species to become the dominant species (Hautier, Niklaus, & Andy, 2009;Partel, Laanisto, & Zobel, 2007;Rajaniemi, 2002). This competitive exclusion allowed the population patches to gradually fuse, and the physiognomy showed a uniform horizontal structure. The first-level classification results of TWINSPAN indicated that the main two microcommunities were dominated by Leymus chinensis and Artemisia capillaries in the enclosure community. Leymus chinensis and Artemisia capillaries outcompeted the other species and had absolute dominant status, which also confirmed the preliminary fusion of population patches.
Thus, the enclosure community had fewer microcommunities than the clipping community.

| Soil variables and community structure under different management modes
Grazing, clipping, and enclosure are the three main grassland management modes in Inner Mongolia, and investigation of the change rule of species diversity and ecosystem multifunctionality is beneficial for sustainable grassland management. All species populations were suppressed by grazing pressure in overgrazed grassland, resulting in miniaturization of individual species and narrowing of niche breadth (Wang et al., 2000), and the photosynthetic area decreased (Carrera, Bertiller, & Larreguy, 2008). As a result, the community coverage, height, and productivity all decreased (Chen et al., 2013;Wang et al., 2015), and when the surface vegetation was damaged, the surface barren area and water evapotranspiration both increased (Wang, He, & Zhou, 2002). Additionally, livestock trampling on soil reduced the soil porosity, osmotic force, and water-holding capacity (Carrera et al., 2008;Wang et al., 2002;Zhang, Han, & Li, 2002). Decreased productivity also reduced litter return, leading to the decline in soil organic matter and nitrogen content (Carrera et al., 2008;Pei, Fu, & Wan, 2008;Wang et al., 2002Wang et al., , 2007 M-indexes were treated as dependent variables. The number of the species with different roles was treated as a fixed-effect factor. "Site" (n = 4) was treated as random factor to address the nonindependence of quadrats in the same sites. *p < .05; **p < .001. matter, soil aggregates and surface crust were damaged (Wang et al., 2010), which resulted in a decrease in clay particles and an increase in sand grains, consequently causing the plant-soil interface to lose balance with concomitant desertification (Wang et al., 2007 Maestre, Quero, et al. (2012) investigated the relationship between species richness and M-index in 224 ecosystems in arid regions worldwide and found that eight optimal models selected from 255 models all included species richness, indicating that species richness is indispensable for explaining variation in multifunctionality. In existing studies of the relationship between species richness and multifunctionality, because a species played a role not only in one ecosystem function but also in other ecosystem functions, the overlap ratio of species supporting different functions simultaneously was 0.2-0.5 (Hector & Bagchi, 2007). Therefore, with increased ecosystem functions evaluated, species richness was positively saturated with ecosystem functions (Gamfeldt et al., 2008;Hector & Bagchi, 2007). In our study, the number of all species and M-index also showed a significant positive correlation, but subdominant and rare species richness did not have high predictive power for the M-index; only dominant species richness had a better predictive power. Therefore, it is important to distinguish the role and rank of different species in the species richness-multifunctionality model; otherwise, the model might include redundant and unclear information. Thus, species richness would be a poor predictor for multifunctionality if the roles of different species were not distinguished to some degree.

| Relationship between biodiversity and multifunctionality
In addition, our identification of a positive correlation between evenness index and M-index was consistent with the results of Maestre, Quero, et al. (2012), who found that species with even distribution could make complementary use of resources more fully, thus increasing ecosystem multifunctionality. Maestre, Castillo-Monroy, et al. (2012) also found that the interaction of evenness with species richness had a significant effect on M-index, but the relative importance of evenness for multifunctionality was lowest when compared separately with the factor species composition or species richness. These results are not surprising, as the effects of the richness × evenness interaction on ecosystem functioning will largely be driven by individual species. When a dominant species has a strong influence on a given function, a negative relationship between evenness and this particular function would be expected because evenness is inversely proportional to dominance. Therefore, in our study, the M-index of monodominant microcommunities was significantly smaller than that of codominant microcommunities; for example, E1 and E3 were significantly smaller than E2 in the enclosure community. This demonstrates that communities with more dominant species whose distribution is also even might have better multifunctionality.

| CONCLUSION
1. The different disturbance intensities and intervals of the three management modes led to distinct differences in interspecific and intraspecific relationships and progress of succession; consequently, the overgrazing community had the fewest heterogeneous microcommunities, followed by the enclosure community, whereas the clipping community had the most heterogeneous microcommunities.

2.
The soil and vegetation showed apparent degradation in the grazing community, and its M-index was the lowest. There were no significant differences in diversity indices and M-index between the clipping and enclosure communities, but these indicators showed obvious differences among the different microcommunities. Therefore, the microcommunity level would be a more suitable scale to investigate the change rule and relationship between plant species and multifunctionality in seminatural ecosystems.