Community assembly mechanisms and succession processes significantly differ among treatments during the restoration of Stipa grandis – Leymus chinensis communities

Understanding community assembly mechanisms is helpful to predict community dynamics. To explore which community assembly mechanism(s) drive(s) the grassland restoration in semi-arid region, we investigated the relationships between plant trait and species relative abundance (SRA), and estimated community functional diversity indices for each community under different treatments (enclosure, grazing and mowing treatment) in a restoration region of Stipa grandis – Leymus chinensis communities in the northern China from 2010 to 2012. There was a high fraction of significant relationships between trait value and SRA, suggesting that niche theory structured the grassland restoration in this region. The functional richness was higher and the functional divergence was lower in the enclosure community than that in the grazing or mowing community, and significantly positive plant height - SRA relationship was found in the enclosure community. These findings demonstrated that limiting similarity based on niche theory was more important in structuring the enclosure community and that environmental filtering based on niche theory played a more important role in driving the grazing or mowing community. Only the factor of year significantly affected the functional evenness (FEve), and the lowest FEve in 2011 implied that the relatively lower precipitation could enhance the effect of limiting similarity on community assembly in the semi-arid grassland.


En-Gr-Mo-En-Gr-Mo-En-Gr-Mo-En-Gr-Mo-En-Gr-Mo-En-Gr-Mo-En-Gr-Mo-En-Gr-Mo-En-Gr-Mo-
Both the Ln RR (the natural-log-transformed response ratio) and the plasticity of plant height across treatments were significantly (P < 0.05) affected by species and Y × S interaction. Both the Ln RR and the plasticity of plant height through time were significantly (P < 0.05) affected by year, species and T × S interaction, and the Ln RR of plant height through time was also significantly (P < 0.05) affected by Y × S interaction. Both the Ln RR and the plasticity of SLA across treatments were significantly (P < 0.05) affected by all these factors (year, treatment, species and their interactions) except T × S interaction. Both the Ln RR and the plasticity of SLA through time were significantly (P < 0.01) affected by all these factors (year, treatment, species and their interactions) except the factor of year and Y × T × S interaction (Tables 1, 2).
The CWM of plant height was significantly (P < 0.001) affected by treatment, year and their interaction. It was significantly (P < 0.05) lower within the grazing community than that within any of the other two communities in 2010 and 2012, and was significantly (P < 0.05) higher within the enclosure community than that within any of the other two communities in 2011 (Fig. 1a). The CWM of SLA was significantly (P < 0.001) affected by treatment, year and their interaction, and it was significantly (P < 0.05) lower within the enclosure community than that within any of the other two communities in any observation year. Compared with the CWM of SLA within the mowing community in the same year, the performance within the grazing community was significantly (P < 0.05) lower in 2011 and significantly higher in 2012 (Fig. 1b). The effects of treatment, year and their interaction on the CWM of one thousand grain weight was not significant (P > 0.05) (Fig. 1c).
Response of species relative abundance (RSRA). The SRA and the RSRA across treatments were significantly (P < 0.05) affected by species, Y × S, T × S and Y × T × S interactions (Table 2). There were four species, S. grandis, C. squarrosa, Allium tenuissimum and Salsola collina, whose RSRA to grazing (grazing vs enclosure) was significantly positive in at least one observation year or no difference with zero. There were three species, L. chinensis, Carex korshinskyis and Chenopodium glaucum, whose RSRA to grazing was significantly negative in at least one observation year or no difference with zero. The RSRA of Artemisia frigida to grazing was significantly positive in 2011, significantly negative in 2010, and no difference with zero in 2012. The significantly positive RSRAs to mowing (mowing vs enclosure) were shown by Agropyron cristatum, C. squarrosa, and A. frigida in 2011 and 2012, and Allium tenuissimum in 2011. The significantly negative RSRAs to mowing were shown by L. chinensis in any observation year, Carex korshinskyi, Salsola collina and Ch. glaucum in 2010 (Table 3).
The RSRA through time was significantly (P < 0.001) affected by treatment, species, T × S and Y × T × S interactions. The significantly (P < 0.05) positive RSRAs to low precipitation (2011 vs 2010) were shown by S. grandis under the grazing treatment and L. chinensis under the enclosure or grazing treatment. The significantly (P < 0.05) negative RSRAs to low precipitation were shown by S. grandis and C. squarrosa under the enclosure treatment, Ch. glaucum and Ch. strictum under any of these treatments ( Table 3). The significantly (P < 0.05) positive RSRAs to high precipitation (2012 vs 2010) was shown by L. chinensis under any of these treatments, Ch. glaucum under the enclosure or mowing treatment, and Ch. strictum under the enclosure treatment. The RSRAs of C. squarrosa to high precipitation was significantly (P < 0.05) negative under any of these treatments (Table 3).  Year (Y)  www.nature.com/scientificreports www.nature.com/scientificreports/ the mowing community in 2011, and significantly (P < 0.05) negative in the enclosure community in 2012. No significant relationship was found between the one thousand grain weight and the SRA in this study (Table 4).
Relationship between response of plant trait and RSRA. The Ln RR of plant height was significantly (P < 0.05) positively correlated with the RSRA to grazing in 2010 and 2012, with the RSRA to mowing in 2011 and 2012, with the RSRA to low precipitation under the enclosure treatment, and with the RSRA to high precipitation under any of these treatments. The plasticity of plant height was significantly (P < 0.05) positively correlated with the RSRA to grazing in 2012, with the RSRA to mowing in 2011, and with the RSRA to high precipitation under any of these treatments. Both the Ln RR and the plasticity of SLA were significantly (P < 0.05) negatively correlated with the RSRA to grazing in 2011, with the RSRA to mowing in 2011, and with the RSRA to high precipitation under the enclosure treatment (Table 5). community functional diversity. The functional richness and functional divergence were significantly (P < 0.05) affected by treatment, year and their interaction. The functional richness was significantly (P < 0.05)  www.nature.com/scientificreports www.nature.com/scientificreports/ higher in 2011 than that in 2012, and the highest functional richness and the lowest functional divergence were found within the enclosure community in 2011. Only the factor of year significantly (P < 0.05) affected the functional evenness, and the lowest functional evenness was found in 2011 (Fig. 2).

Discussion
The present study found that the percentage of significant plant trait -SRA relationships was higher than random expectation (5%), indicating that niche theory played more important role in affecting the community restoration in the study region. The environmental adaptability or tolerance of these co-occurring species could provide evidence that not all individuals were functionally equivalent and support the niche theory. The RSRAs across treatments were significantly different among the co-occurring species, implying the difference of adaptability and tolerance to disturbance among these species, which could cause different restoration outcomes due to the change of competition superiority of co-occurring species 17,23 .
In the enclosure community, the relative abundance of S. grandis or L. chinensis was significantly higher than that (<5%) before 2003, and the relative abundances of S. grandis and L. chinensis were the first two highest among all these species in 2012 (Table 1), implying the dominance of S. grandis and L. chinensis in this community. That is to say, the composition and structure in the enclosure community was close to the climax community in this region. The CMW of SLA was lower within the enclosure community than that within any of the other communities, and the SLA of the dominant native species, S. grandis and L. chinensis, was relatively lower than that of any of the other species, implying that the relatively lower SLA has the higher functional adaptability in this region. In semi-arid region, species with low SLA and therefore long-lived leaves can accumulate a greater mass of leaf, and long mean residence time of nutrients, made possible by leaf longevity, permits a progressively larger share of nitrogen pools to be sequestered 24 .
However, in the grazing/mowing community, the abundance of C. squarrosa was highest (>30%) among all these species (except in the mowing community in 2012), indicating that C. squarrosa but not S. grandis or L. chinensis was the dominant species. More time was needed for the grazing or mowing community to reach the status of S. grandis -L. chinensis community. The significant differences of RSRAs through time with zero suggested that these communities were more dynamic, or hysteresis may be occurring in response to the grazing/  Table 5. Pearson correlations between responses of plant height or specific leaf area and the corresponding RSRAs across treatments or through time. *, ** and *** indicate significant relationships at the 0.05, 0.01 and 0.001 level, respectively.
www.nature.com/scientificreports www.nature.com/scientificreports/ mowing disturbance 25,26 . The RSRAs of L. chinensis to grazing or mowing were significantly lower than zero while the RSRAs of C. squarrosa to grazing or mowing were significantly higher than zero except one that was no difference with zero, implying the relatively higher tolerance of C. squarrosa to grazing/mowing than L. chinensis. It has been reported that C. squarrosa showed higher tolerance to environmental disturbance or stress than S. grandis in the Inner Mongolia Steppe 27 . As a result, the relatively lower recovery (competitiveness) of L. chinensis population and the relatively higher competitiveness of C. squarrosa population in the grazing or mowing community than in the enclosure community. A disturbance could damage or destroy some vegetations, reducing the uptake of light, water and nutrient. Once resource uptake goes down for a time in an ecosystem, there will be more resource available and the competition among co-occurring species will decrease, thus the community is particularly vulnerable so that invasion or the dominant species replaced by high resource-capture species happens very easily 28 . Contrary to species with low SLA, species with high SLA, take C. squarrosa as an example, can capture resource rapidly, in sequence, become a superiority to occupy the gaps in a disturbance community such as a grazing or mowing community. Consequently, the disturbance could lead to increasing the competition superiority of C. squarrosa, and then increasing its SRA. In fact, when there is no disturbance, C. squarrosa is an inferior competitor; and even if disturbance exists, C. squarrosa is still an inferior competitor in high precipitation year 29 , as the results in Table 3 showed in this study. The significantly different of characteristic (adaptability, tolerance, SLA, plant height and so on) between C. squarrosa and S. grandis or L. chinensis could explain not only the degradation processes in this study region but also the restoration processes under different treatments.
Even though niche theory affected the community assembly in this region, the limiting similarity based on niche theory was more important in structuring the enclosure community and the environmental filtering based on niche theory played a more important role in driving the grazing or mowing community 3,10,30 . The FRic within the grazing/mowing community was lower than that within the enclosure community or that reported by other researchers 31,32 , therefore, it could be viewed as evidence of environmental filtering affecting the grazing/mowing community 3 . While the relatively higher FDiv within the enclosure community could be considered as evidence of limiting similarity among co-occurring species 10 . Moreover, the significantly positive plant height -SRA relationships in the enclosure community demonstrated that the higher species had more advantage, supporting that light competition was the driving force as community height and canopy density increased over years of enclosure and that limiting similarity and niche differentiation by light competition affected the community structure and processes in the enclosure community 33 . However, the dramatic inter-annual changes of the correlation coefficient or significance under the grazing/mowing treatment emphasized the importance of environmental filtering by inter-annual climatic changes in driving the community assembly 34 . For example, in the mowing community, the correlation coefficient between the SLA and the SRA was 0.602 and more higher than zero in 2011 but −0.126 in 2012 (Table 4), and similar results have been reported by many researchers 4,35 . From the restoration outcomes under different treatments, we conjectured that the limiting similarity played more important role in community closer to the climax community in the study region.
The present findings indicated that the factor of year was a significant factor interacting with restoration treatments in affecting the community dynamics in the semi-arid grassland. However, only the factor of year significantly affected the FEve, and the lowest FEve was found in 2011. According to viewpoint of Pakeman 11 , low level of FEve can be indicative of sites where the limiting similarity may be important in structuring the communities. Our finding implied that the relatively lower precipitation in 2011 enhanced the effect of limiting similarity on community assembly. Seed mass is quite a good indicator of a cotyledon-stage seedling's ability to survive various hazards, representing a species' chance of successful dispersing a seed into an establishment opportunity 36 . The relationship between one thousand grain weight and SRA was negative in 2011 but positive in 2012, suggesting the trend of community succession and the important effect of year-to-year (climatic) fluctuation on recruitment.
Our results supported that the restoration succession of the semi-arid grassland in northern China was sensitive to environmental conditions, including restoration treatments or year-to-year fluctuation, and demonstrated that long-term consecutive observations of plant traits and estimating FD indices are vital for understanding www.nature.com/scientificreports www.nature.com/scientificreports/ community assembly mechanisms. RSRAs across treatments or through time are helpful to calculate trajectories of community process and judge community state. Our findings provided a more comprehensive understanding of why some species are more abundant in one particular habitat but less abundant in another habitat, but also demonstrated the complexity of community assembly mechanisms and processes of grassland restoration in the northern China.

Methods
Study region. The study region (44°10′~44°20′N, 116°20′~116°30′E, 1110~1154 m a.s.l.) located at Xilingol League, the middle of Inner Mongolia Steppe of China where S. grandis -L. chinensis community is the typical and climax vegetation type. Before 2003, S. grandis -L. chinensis community in this region degraded heavily due to overgrazing by cattle or/and sheep, the relative abundance of C. squarrosa was more than 30% while that of S. grandis or L. chinensis was less than 5%. That is to say, C. squarrosa community replaced S. grandis -L. chinensis community. The soil type is typical chestnut soil (Chinese Soil Taxonomy System) and Calcic luvisol (FAO-UNESCO). The mean annual temperature is 0.8 °C, ranging from −39.9 °C in January to 37.4 °C in July, with approximately 106 frost-free days per year. Mean annual precipitation measured over the last 30 years is 263.5 mm 22 , falling mainly from July to September. experimental design. In 2003, over-grazing was forbidden and three different restoration treatments were carried out in the study region. One-third was grazing plot, where light grazing was performed with approximately 3 sheep/ha from May to October. Two-thirds area was enclosed, and separated into enclosure plot and mowing plot. In the enclosure plot, there was no human-being disturbance; and in the mowing plot, the aboveground plants was mowed once in every autumn, leaving stubble height about 3~5 cm. In at the end of July, avoiding any overlap with the sampling plots set up before. Four 1 m × 1 m quadrats per sampling plot were randomly arranged with the constraint that it was at least 0.5 m from the edge to avoid marginal effect, totalling 288 sampling quadrats during these three years.
The species were not present at all quadrats, species that appeared in higher than 10% sampling plots were thought of as common species and the functional traits of common species were collected for further analysis. These species presented more than 80% of the total above-ground biomass and vegetation cover, therefore, according to the viewpoint of Cornelissen et al. 37 , the result may provide useful information to scale-up the trait values from species to the community level.
A database of functional traits for each species observed during the plant composition surveys using standard methods 38 . We measured two traits from plants in the sampling plot (plant height and SLA), one trait in the study region (seed mass), and two categorical traits from literature (functional groups [perennial graminoids, perennial forbs and annual plants] and seed diaspore type [seed, caryopsis and achene]). We selected these traits because they are main predictors for RSRA across treatments or through time 2,36 . The plant species density and the plant height of 3~5 mature individuals per species were recorded within each quadrat. Mature, fully expanded leaves of each species were collected to measure SLA which was calculated as the ratio of leaf area to dry leaf mass 37 , three replicates per sampling plot. The leaves were scanned for measuring leaf area before being dried at 80 °C for 48 h and weighed for biomass to the nearest 10 −4 g. For species with needle-like leaves or small leaves, four to six full leaves were thought of as one sample. The SLA values of both A. frigida and Carex korshinskyi showed a high random error, therefore, they were not used for further analysis. Fresh diaspores (referring to the seed and dispersal structures) were collected from May to September in 2010-2012. Diaspores were collected from as many individuals as possible (n > 20) in the study region. Structures having the function of contributing to dispersal were not included as part of seed mass 39 . One thousand grain weight was measured for each species as the characteristic of seed mass, ten replicates per year. We did not measure the one thousand grain weight by treatment, because there were not enough diaspores within the grazing community. Data analysis. Based on data collected from the same sampling plot, the mean values of plant height and SLA were calculated for each species, and the SRA of each species was calculated as the ratio of a given species' density to the total species' density. The CWM trait value was calculated with equation, where SRA ij is the SRA of the species i in sampling plot j, T ij is the mean trait value of the species i in sampling plot j, and CWM j is the community-weighted trait of sampling plot j. The Multi-trait FD indices, including FRic, FEve and FDiv, were calculated based on the species abundance and the 'Gower' distance matrices of the five traits (species functional group, diaspore type, plant height, SLA and one thousand grain weight). FRic represents the amount of functional space occupied by a species assemblage. FEve corresponds to how regularly species abundances are distributed in the functional space. FDiv defines how far high species abundances are from the center of the functional space. These three facets of FD are complementary and, taken together, describe the distribution of species and their abundances within the functional space 30 . If a plant trait is absent from one sampling plot, it is treated as missing data. If a species is absent from one sampling plot, it had a value of zero. These indices were calculated by using the 'FD' package 40 41 . T G , T M and T E were the performance of plant height or SLA of species i measured in the same sampling belt in the same year in the grazing, mowing and enclosure plots, respectively; and T L , T H and T N were the performance of plant height or SLA of species i calculated in the same sampling belt under the same treatment in 2011, 2012 and 2010, respectively.
The RSRAs across treatments or through time obtained for each species are symmetric around zero. The difference significance with 0 for each species' RSRA was examined by t-test, where a significantly higher value than zero, a significantly lower value than zero and no difference with 0 indicate a positive response, a negative response and a non-significant response across treatments or through time, respectively. In addition, the effects of treatment, year, species and their interactions on plant traits, traits' responses, SRA and RSRA were analyzed by linear mixed model restricted maximum likelihood (REML) with species as a fixed factor, treatment and year as random factors. The REML analyses allow the specification of fixed and random model structures for the analysis of unbalanced data set. Four variables, including SRA, plant height, SLA and one thousand grain weight, were log-transformed to meet the normality. Moreover, the effects of year, treatment and their interaction on the three FD indices (FRic, FEve and FDiv) and the CWM of each trait were analyzed by general linear model -univariate analysis with treatment and year as fixed factors and the FRic being log-transformed. At last, the relationships between each plant trait (plant height, SLA and one thousand grain weight) and the SRA, and between responses of plant height or SLA and RSRA were analyzed by Spearman correlation analysis. All these analyses were performed by using SPSS 21.0 (IBM, Chicago, IL, USA).