Field study was done at the observation field in Naiman Desertification Research Station of the Chinese Academy of Sciences in the eastern part of Inner Mongolia, China (42°55′-42°56′N, 120°41′-120°42′E) (Figure 1). Naiman located in the Northern China’s agro-pastoral ecotone in the southwestern part of the Horqin Sandy Land, which is roughly 400 × 400 km in size and represents the most desertification threatened area in north China (Li et al., 2005). This area has a continental semi-arid monsoon climate in the temperate zone. The mean annual precipitation is 377 mm, and nearly 75% is concentrated in the growing season (May to August) (Zuo et al., 2012). The average annual temperature is 6.8°C. Vegetation cover in the study area ranged from 50% to 70%.

The experimental field consisted of four plots along a restoration gradient, i.e., Plot A, natural mobile dune grassland plot (NM); Plot B, post-non-grazing restoration plot (NGR); Plots C and D, post-light-grazing restoration plot (LGR) and post-heavy-grazing restoration plot (HGR). NM, NGR, LGR and HGR were long-term observation fields, and many studies have been carried out (Su et al., 2006; Zhao et al., 2009; Zhang et al., 2016). The NM Plot has been monitored since 2005 and has an area of 0.60 ha (200 × 30 m). The vegetation cover was less than 5%. The NM Plot was dominated by the annual herb Agriophyllum squarrosum. The NGR Plot was established and fenced in 1985 to assist with the natural recovery of the sandy grassland and covers 0.66 ha (200 × 33 m). The vegetation cover in the NGR Plot ranged from 30% to 50%. The NGR Plot was dominated by the perennial grasses Artemisia halodendron, Melissitus ruthenicus, and Euphorbia humifusa. The LGR and HGR Plot were established in 1992 to allow natural restoration of the 5-years lightly and heavily grazed lands, has 1.5 ha (200 × 75 m) each. The grazing intensities were 2 sheep units/hm² and 6 sheep units/hm², respectively. Vegetation cover in LGR Plot was more than 70%. The dominant plant species were annual herbs, including Artemisia scoparia, Setaria viridis, Salsola collina, and Corispermum hyssopifolium (Niu et al., 2018). Vegetation cover in the HGR Plot was about 60%. The HGR Plot is dominated by the perennial grass, Pennisetum centrasiaticum, S. viridis, Digitaria cilliaris, and Phragmites communis (Zhao et al., 2007).

The number of plants was determined in mid-August 2019, at which time coverage, height and frequency were measured and used to calculate their importance value. In each plot, 9 random quadrats (1 × 1 m) were used to measure vegetation characteristics. Young healthy leaves from 3 dominant plant species in each quadrat were arbitrarily sampled from plants and immediately stored with silica gel in zip lock plastic bags for later DNA extraction. Aboveground biomass (AGB) in each quadrat was harvested in August, standing biomass and litter separated, oven-dried at 65°C for 48 h and weighed (Zhao et al., 2016).

About 36 soil samples (0–10 cm; 9 samples per plot) were collected from the 4 recovering grassland plots. Soil samples were air-dried and hand-sieved through a 2-mm screen to remove roots and other debris. A portion of each air-dried soil sample was ground to pass a 0.1 mm mesh for soil nutrient analyses. Soil organic carbon (SOC) concentration was determined using the Walkley-Black dichromate oxidation procedure (Nelson and Sommers, 1982). Soil available nitrogen (SAN) was determined using the alkaline diffusion method, and soil phosphorus (SOP) was determined using the Olsen method (Olsen et al., 1954). Soil samples for molecular analysis were packed in polyethylene bags, immediately stored in portable refrigerator powered with a car battery and then stored in the lab at −80°C until processing.

Plant and soil total DNA was extracted using an AxyPrep genomic DNA mini kit (Axygen, Beijing, China) and Qiagen Dneasy Powersoil kit (Qiagen, Beijing, China), respectively. DNA was quantified spectrophotometrically; samples that yielded high quantities of good-quality DNA were used in consecutive experiments.

Soil bacterial and fungal communities were sequenced separately in 16S rDNA (amplified by primers 515F/806R) (Walters et al., 2016) and ITS rDNA (amplified by primers ITS5-1737F/ITS2- 2043R) (Zhou et al., 2020) (Appendix S1). Polymerase chain reactions (PCR) was performed in a 20 µL reaction volume, containing 10 ng of template DNA, 0.4 µL Pfu DNA polymerase (TransGen Biotech, Beijing, China), dNTPs (4 µL of 2.5 mM), 5 × reaction buffer (4 µL), each primer (0.8 µL of 5 µM). The amplification condition for bacteria are as follow: 95 °C for 3 min, by 27 cycles of 30 s at 95°C, 30 s at 55°C, 45 s at 72°C, and a final 10 min extension step at 72°C. The amplification condition for fungal are as: 94 °C for 1 min, followed by 30 cycles of 20 s at 94°C, 30 s at 56°C, 30 s at 72°C, and a final 5 min extension step at 72°C. High-throughput sequencing was performed on an Illumina NovaSeq platform (Illumina, San Diego, CA, United States) in Novogene Co., Ltd. (Beijing, China).

In this study, DNA barcoding markers were used to assess genetic diversity and variation in dominant plants. Several pairs of plant cpDNA and nrDNA primers designed by Hamilton (1999), Taberlet et al. (1991) and Sang et al. (1997), were used in the initial screening (Appendix S1). Marks trnL-F, matK and psbA-trnH were chose as primers for plants from Chenopodiaceae and Leguminosae amplification; ITS, rps16 and matK as primers for plants from Gramineae amplification; trnL-F, ITS and psbA-trnH as primers for plants from other families amplification. PCR was performed in a 25 µL reaction volume, containing 40 ng of genomic DNA, 1.0 U of Taq polymerase (Axygen Inc., Beijing, China), 3 mM MgCl₂, 500 µM each dNTP, 20 mM Tris-HCl (pH 8.3), 100 mM KCl, 0.3 µM each primer. The amplification condition was an initial denaturation step of 94°C for 3 min, followed by 30 cycles of 30 s at 94°C, 30 s at an appropriate annealing temperature, 1 min at 72°C, and a final 5 min extension step at 72°C. The PCR products were determined by 1.0% agarose gel electrophoresis. The amplification products were purified using AxyPrep PCR purification kit following the manufacturer’s protocol (Axygen Inc., Beijing, China). Purified plant DNA was sequenced by the MEIJI sequencing company in Shanghai, China, applying the PCR-primers as sequencing primers.

Two measures of plant species diversity from 4 recovering grassland plots were used: species richness (S_R) and Simpson’s diversity index (D). D was calculated from the formula, D = ∑ i S N i 2 where Ni = important value, which is the average of relative density, relative height, and relative biomass.

The number of species was taken as a measure of species richness. The measures of plant genetic diversity from recovering gradients plots were calculated with Arlequin 3.0 (Excoffier et al., 2005), including allelic richness (A_R) and expected heterozygosity (H_E), as they are analogous to S_R and D, respectively (Etienne, 2005; Evanno et al., 2009).

The raw reads of soil DNA sequencing were demultiplexed, quality-filtered by FASTQ, and merged by FLASH (Magoč and Salzberg, 2011). Quality-filtered sequences were aligned with SILVA alignment (Quast et al., 2013), and clustered into the operational taxonomic units (OTUs) with 97% sequence similarity using UPARSE (Edgar, 2013). Two measures of soil microbial diversity from recovering gradients plots were used: soil bacterial operational taxonomic unit (OTUs) and soil fungal richness.

The relationships between plant species diversity (S_R) and plant genetic diversity (A_R and H_E) and between plant diversity (P_D) and soil microbial diversity (S_D) were examined. Before analysis of variance (ANOVA) and correlation analysis, data were tested for normality to ensure normal distribution. Statistical significance was determined by one-way (plant species and soil microbial diversity index) or two-way (plant genetic diversity index) ANOVA followed by Tukey’s post hoc test. All statistical analyses were performed using SPSS 19.0 for windows (SPSS, Inc. Chicago, IL). Confirmatory analysis of structural equation model (SEM) was applied to the data to determine direct and indirect impacts of long-term ecological restoration on above-and belowground biodiversities. SEM analysis used lavaan packages in R version 3.5.2 (https://www.r-project.org/).

Thirty-three species were recorded (Table 1). Specifically, there are 9 species in natural mobile dune grassland plot (NM), 19 species in post-non-grazing restoration plot (NGR), 14 species in post-light-grazing restoration plot (LGR) and 18 species in post-heavy-grazing restoration plot (HGR). The results showed that herbaceous are dominated in all grassland plots. Biennials were dominant in NM and NGR, accounting for 92.40% and 78.22%, respectively. Perennials were dominant in LGR and HGR, accounting for 57.32% and 44.35% of the total importance value. According to the species importance value, the results showed that S. viridis and C. incertus were the dominant plant species in NM, A. scoparia and E. gmelini were the dominant plant species in NGR, E. stephanianum and C. squarrosa were the dominant plant species in LGR, S. viridis and P. centrasiaticum were the dominant plant species in HGR. One-way ANOVAs explained that grassland restoration were notable impacted on plant species richness (S_R) and Simpson’s diversity index (D) (p < 0.001) (Figure 2). The results indicated that the plant species richness from the NGR communities was significantly higher than that from the other recovering grassland plots (10 ± 1.1, p = 0.004), and the plant species richness from the NM communities was significantly lower than that from the other recovering grassland plots (4 ± 1.2, p < 0.05). The plant Simpson’s diversity index in the NM communities was remarkably higher than that in other recovering grassland plots (0.417 ± 0.036, p < 0.001) and exhibited no significant change in the NGR, LGR, and HGR communities.

Genetic diversity of dominant plant species from different recovering grassland plots is given in Figure 3. Two-way ANOVAs revealed that genetic diversity indices were affected significantly by different grassland plots, and the allelic richness (A_R) and expected heterozygosity (H_E) of dominant plant species in cpDNA was outstanding greater than those in nrDNA in each recovering sandy grassland plots. Among the restoration plots, the A_R of dominant plant species in nrDNA and cpDNA was varied from 1.931 to 2.590 and from 3.844 to 3.930, respectively. The A_R values of nrDNA and cpDNA in NGR were significantly higher than those in NM. The H_E of the dominant plant species in nrDNA and cpDNA varied from 0.2989 to 0.4822, and from 0.6924 to 0.7007, respectively, and the values of each recovering grassland plot were not significant (amplitude of difference up to 44.8%–70.5% in A_R, 81.9%–128.1% in H_E).

The bacterial operational taxonomic units (OTUs) were assigned into 37 phyla, 216 families and 487 genera. According to the annotated OTUs, Actinobacteria, Proteobacteria and Acidobacteria represented the most abundant phyla in all samples, comprising 21.87%–53.15%, 27.62%–45.37%, and 3.31%–18.60%, respectively, of the total OTUs in each treatment. These three phyla occupied more than 80% of each plot. The fungal OTUs were assigned into 16 phyla, 220 families and 388 genera. According to the annotated OTUs, Mortierellomycota and Ascomycota represented the most abundant phyla in all samples, comprising 23.59%–44.90% and 21.47%–44.94%, respectively, of the total OTUs in each treatment. These two phyla occupied more than 50% of each plot (Figure 4).

In different restoration grassland plots, one-way ANOVAs identified that community types had sharp influence on soil microbial richness (p < 0.001) (Figure 5). The results found that soil bacterial and fungal richness from NM communities were obvious fewer than recovering grassland plots (1641.9 ± 100.4, p < 0.001; 533.2 ± 16.6, p < 0.001). The soil bacterial richness from the HGR communities was greater than that from the NM and NGR communities (2714.1 ± 38.8, p < 0.05). However, the soil fungal richness demonstrated no clear change in the NGR, LGR and HGR communities (Figure 5).

The heterogeneous levels of genetic variability among different restoration grassland plots were detected. Composite allelic richness of cpDNA and nrDNA in NM, NGR, LGR and HGR communities varied from 1.02 to 1.78, 2.25 to 2.37, 1.72 to 2.32, 1.74 to 2.36, respectively. Composite expected heterozygosity of cpDNA and nrDNA varied from 0.34 to 0.78, 0.72 to 0.77, 0.72 to 0.76, 0.71 to 0.77, respectively. The heterogeneous levels of soil microbial diversity among different recovering grassland plots were revealed. Composite soil microbial shannon diversity of bacterial and fungal in NM, NGR, LGR and HGR communities varied from 4.728 to 6.253, 6.210 to 7.173, 6.616 to 7.090, 6.138 to 7.349, respectively. Composite plant diversity of species diversity and genetic diversity varied from 0.164 to 0.368, 0.099 to 0.496, 0.217 to 0.579, 0.147 to 0.355, respectively (Figure 6).

Correlation analyses disclosed that there were positive correlations in species diversity and genetic diversity (Figure 6). Data illustrated that there was a prominent correlation between S_R and A_R (r = 0.56, R ² = 0.31, p < 0.001), and a striking correlation between S_R and H_E (r = 0.33, R ² = 0.11, p = 0.045). The results also showed that there was a striking negative correlation between plant diversity (P_D) of different grassland types and soil microbial diversity (S_D) (r = -0.44, R ² = 0.19, p = 0.005) (Figure 6).

The final structural equation model (SEM) for projecting the direct and indirect effects of long-term restoration on S_R, A_R and S_D emerged when aboveground biomass (AGB), soil organic carbon (SOC), soil alkali-hydrolyzable nitrogen (SAN) and soil phosphorus (SOP) were the variables (P_Chi-squre = 0.914, P_RMSEA = 0.979, df = 11). The final model explained 38% of the variance of S_R, 57% of A_R, 52% of S_D, 49% of AGB, 87% of SOC, 47% of SAN and 69% of SOP (Figure 7). Long-term ecological restoration had outstanding direct active effects on AGB, SOC, SAN, SOP, A_R, S_R and S_D. There was an inactive correlation between aboveground and belowground biodiversity and biological and abiotic factors; that is, S_R was significantly negatively correlated with AGB, A_R was significantly negatively correlated with AGB, SAN and SOP, and S_D was negatively correlated with SAN but not significantly. The effect of SOC on aboveground plant diversity was mediated by AGB and SAN.

Long-term ecological restoration greatly affects the structure and function of grassland ecosystem, especially the vegetation characteristics. Different forms of restoration will have different influences. The results showed that the number of species in NGR, LGR and HGR was more than NM (Table 1), which might be due to ecological restoration had improved the original habitat conditions and increased the number of species suitable for living. The data pointed out that the proportion of herbaceous was higher than that of shrubs in all treatments, and perennial plants were the dominant species in LGR and HGR, and biennial plants were the dominant species in NM and NGR (Table 1). These results indicated that long-term ecological restoration affected the composition of life forms of grazing and natural grassland vegetation and changed the community structure. It is consistent with the research results of Lv et al. (2016) and Zhao et al. (2008).

In this study, long-term ecological restoration altered soil bacterial and fungal communities and different microbial taxa responded unequally. The results stated that in all restoration gradients, the dominant phyla of bacteria and fungi were Actinobacteria, Proteobacteria and Acidobacteria, Mortierellomycota and Ascomycota, respectively, which was consistent with relevant research datas in Horqin Sandy Land (Wang et al., 2018; Zhang et al., 2020). This indicated that soil microbial community composition was relatively stable in sandy grassland with different restoration modes. The abundance of Proteobacteria and Acidobacteria in the restoration plot (NGR, LGR, HGR) was higher than that in the NM. These results indicated that ecological restoration promoted the bacterial transformation from slow growth type to fast growth type (Xun et al., 2018; Zheng et al., 2022). The reason was that vegetation regrew after enclosed-fencing, soil nutrient content increased, and nutrient-loving bacteria were promoted to grow stably and multiply rapidly. The relative abundance of Mortierellomycota in NGR and LGR were much greater than that in NM and HGR, while that value of Ascomycota in HGR were much greater than that in other plots. This is the same as the research results of Wang et al. (2019) on the short-term grazing ban (14 and 19 years) in semi-arid steppe, that is, the diversity of soil fungi becomes richer with the increase of the years of enclosure. The reason may be that the increase of vegetation biomass and coverage in the process of ecological restoration changes the rate of soil nutrient conversion and migration, and accelerates the demand for saprophytic nutrient fungi in grassland ecosystem (Zhang et al., 2013) (Figure 4).

Ecological restoration measures are of great significance to the sustainable development of sandy grassland ecosystem. In this study, the results discovered that enclosure had a certain effect on aboveground plant species richness, allele diversity and belowground soil fungi and bacteria richness in sandy grassland (Figures 2, 3, 5). The species richness (S_R) from NGR was significantly higher than that of other restoration plots (Figure 2), because the original habitat conditions of NGR are more superior than other plots, which is more conducive to plant growth and community construction. The overall trend was consistent with that the aboveground plant diversity (S_R and A_R) and belowground soil microbiology diversity (bacterial OTUS and fungal richness) with NM were significantly lower than those of other treatments (Figures 2, 3, 5). This is due to the poor original habitat conditions in NM lead to low species richness and soil microbiology diversity, and habitat fragmentation leads to reduced allelic richness. The soil bacterial OTUS in HGR was outstanding more than that of other treatments, possibly because severe grazing had positive and lasting effects on soil bacterial community. There were no significant differences in Simpson’s diversity index (D), expected heterozygosity (H_E) and soil microbiology diversity values from NGR, LGR and HGR indicating that after years of restoration and conservation, the grazing grassland of different intensities in the early stage has been repaired well, and the habitat conditions are similar to the natural grassland in the same period.

The relationship between above- and belowground biodiversity in Horqin sandy grassland community was significantly affected by long-term ecological restoration, but the trend was not consistent. The data detected a positive correlation between plant species diversity and genetic diversity parameters, which was supported by many similar studies (He et al., 2008; Crawford and Rudgers, 2013). In this study, S_R was extremely significant vigorous correlated with A_R, and was significantly vigorous correlated with H_E, indicating that plant species richness determines the genetic differences of dominant species and plays a key role in the level of genetic diversity of a sandy grassland community. The results showed that there was a significant negative correlation between above- and belowground biodiversity in sandy grassland community, which suggested that plant species richness and genetic diversity index inhibited soil bacteria and fungi diversity in enclosure conditions (Figure 6).

Previous studies have shown that different management measures responded differently to biodiversity conservation (Lindenmayer et al., 2006; Kernecker et al., 2022). According to the SEM in the present study, long-term ecological restoration did active affect above- and belowground biodiversity directly in sandy grassland (Figure 7). However, biological and abiotic factors did inactive affect them directly. AGB was the main direct driving factor for S_R. SAN and SOP were the main direct driving factor for A_R. SAN was the main direct driving factor for S_D. SOC was affected plant diversity through AGB and SAN, and was an indirect driver of plant diversity.