This study was carried out at the Natural Reserve of Ganshenling (109°37'−109°41'E, 18°19'−18°24'N), which has a typical tropical lowland rainforest with a tropical marine monsoon climate. The annual average precipitation and temperature are 1,800 mm and 25.4°C, respectively. The dominant plant species in this area include Hopea reticulata, Vatica mangachapoi, Koilodepas hainanense, Rhodomyrtus tomentosa, Melastoma sanguineum, Scleria elata, and Blechnum orientale (Qi et al., 2014; Mao et al., 2021).

Four successional stages in the Natural Reserve of Ganshenling were selected: (1) G (109°39'−109°41'E, 18°21'−18°23'N); (2) S (109°39'−109°41'E, 18°22'−18°23'N); (3) SF40 (109°39'−109°40'E, 18°21'−18°23'N); and (4) SF60 (109°39'−109°40'E, 18°22'−18°23'N).

From each successional stage, five soil samples of 0–10 cm depth were collected by following the 5-point sampling method from five plots of 900 m² (30 m × 30 m) in August 2019 (Wang et al., 2020). These samples were divided into two groups: one was stored at −80°C until the measurement of soil AMF communities; the other group was used to measure soil physicochemical properties. These soil samples were subsequently air-dried. A 2-mm screen sieve was used to sieve air-dried soil samples for analyzing soil pH, SOC, SN, and soil phosphorus (SP). The SOC was measured with the help of the H₂SO₄-K₂Cr₂O₇ oxidation method (Bisutti et al., 2004). SN was measured using the Micro-Kjeldahl method (Bremner, 1960), and the available nitrogen was estimated with the help of a micro-diffusion method (Mulvaney and Khan, 2001). SP was measured using an inductively coupled plasma atomic-emission spectrometer (Thermo Jarrell Ash Corporation, Franklin, USA) and HNO₃-HClO₄ soil solution (McDowell and Sharpley, 2001).

Soil DNA was extracted from 250 mg soil samples with the help of a MO-BIO PowerSoil DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA, USA). DNA was analyzed and quantified with the help of a NanoDrop ND-2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). A standard concentration of AMF was used for PCR amplification.

For PCR amplification, AMV4.5NF (5′-AAGCTCGTAGTTGAATTTCG-3′) and AMDGR (5′-CCCAACTATCCCTATTAATCAT-3′) were used as primers. An AMF library was constructed using two-step PCR amplification. Target fragments were amplified using these primers. Then, an Illumina adapter sequence, a primer pad, and a barcode were added to both the forward and reverse target fragments. A 50-μL reaction volume was used for PCR of each sample, comprising 10 μL of 5x buffer, 1 μL of dNTP (10 m/M), 1 U of Phusion High-Fidelity DNA Polymerase (New England Biolabs, Ipswich, MA, USA), 1 μL of forward and reverse phasing primers, and 3–10 ng of target sequencing. The reaction volume was supplemented with 50 μL of RNase-free ultrapure water. An ABI9700 PCR amplifier (Applied Biosystems, Waltham, MA, USA) was used to amplify DNA. The holding temperature was set at 10°C. PCR products were run on an agarose gel and PCR products were purified from the gel with the help of an AxyPrep DNA Gel Recovery Kit (Axygen Biosciences, Corning, NY, USA) and pooled. FTC-3000TM real-time PCR instrument (Richmond Biotechnology Inc., Richmond Hill, Ontario, Canada) was used to quantify purified DNA and mixed it to an equal molar ratio to perform the second PCR amplification, for which the reaction volume changed to 40 μl: 5x buffer (8 μl), dNTP (1μL, 10 m/M), Phusion High-Fidelity DNA Polymerase (0,8 U, New England Biolabs, Ipswich, MA, USA), both forward and reverse primers (1 μl, 5 μmol/L), respectively, and DNA template (5 μl, 100 μg/ml) and RNase-free ultrapure water (5 μl).

Barcodes were used to identify sequences of different samples in a parallel sequencing test. The sequence with a complete barcode was regarded as valid and extracted to a dataset stored in the FASTQ format. Quality assurance and splicing were performed on such valid sequences and included the following steps: (1) Poorly overlapped and low-quality sequences (quality score <20 and window size <5) were removed with the help of the Trimmomatic program (Bolger et al., 2014); (2) paired reads were merged into a sequence, where the overlap of PEreads (>10 bp) was evaluated with the help of the FLASH program (Magoc and Salzberg, 2011); and (3) merged sequences were screened according to an error ratio (<0.2). Mothur V.1.39.5 was used for quality control of merged reads by removing sequences that were ambiguous, homologous, or chimeras generated during PCR amplification (Schloss et al., 2009). Further analysis of information on species and clustering of operational taxonomic units (OTUs) was based on the clean tags. OTUs were selected with the help of UPARSE with a cutoff of 97% similarity (Edgar, 2013). The dataset was screened for chimeric reads using UCHIME in a reference database mode (Edgar et al., 2011). Through USEARCH global, all sequences were compared to OTU representative sequences to obtain the OTU abundance table. The identity of OTU taxonomy was determined with BLAST (Ye et al., 2006) by comparing them to representative sequences in the MaarjAM database (Zhang et al., 2000; Opik et al., 2010). The taxonomy of an OTU was obtained for the species that had the highest similarity or identity with OTU representative sequences in the MaarjAM database. Some OTUs didnot contain any species information and were therefore marked as “Nomatch”. The library preparation and sequencing were performed by TinyGene Bio-Tech (Shanghai, China). Raw sequence data were submitted to NCBI's Sequence Read Archives (SRP399604).

We analyzed the overall AMF species richness and diversity (Shannon H' and Simpson indices) and equitability (Pielou J'). Following the classification proposed by Logares et al. (2014), we classified species' relative abundances within each successional stage as “dominant” (>1%), “common” (1–0.01%), and “rare” (<0.01%) (Galand et al., 2009).

OTUs that included an unknown genus and contained <5 sequences were not used for bioinformatics analysis. The species diversity was calculated in R package “vegan.” The principal coordinate analysis (PCoA) was performed using Bray–Curtis distances with 999 permutations with the help of R 4.2.2. The linear discriminant analysis effect size (LEfSe) and linear discriminant analysis (LDA) were performed with R 4.2.2 (Sirova et al., 2018; Nauta et al., 2020). The phylogenetic tree and heatmap of AMF were analyzed with the help of R 4.2.2.

Bray–Curtis-based PCoA showed that successional stages had different effects on AMF communities (Figure 1). The AMF communities clustered according to different successional stages, demonstrating that there were large differences between G and S, as well as between SF40 and SF60 (Figure 1).

The PERMANOVAs of Bray–Curtis distances also revealed significant differences among the AMF communities at different successional stages (P = 0.001). Additionally, our beta-diversity results were not significant (P = 0.149), indicating that our groups had the same dispersions. This fact also confirmed that our PERMANOVA results were reliable. Successional stages had a significant effect on the overall AMF community, which contributed to 33.7% of the total variation (Figure 1).

To further elucidate the trend in recovery of the AMF communities during successional stages, the similarity of community composition at the OTU level between different successional stages and SF60 was calculated. There was an average Bray–Curtis similarity of 35.6% between all AMF communities and SF60. The average Bray–Curtis similarity of AMF community between the different successional stages and SF60 was 21.0% for G, 16.0% for S, and 25.3% for SF40, respectively (Figure 2). This suggested that there was a trend of increasing similarity of soil AMF communities from S and G to SF40 compared to SF60.

The alpha-diversity indices that were based on phylogenetic species richness and phylogenetic diversity of soil AMF communities exhibited different trends (Supplementary Figure S1). AMF phylogenetic species richness and phylogenetic diversity of G were significantly lower than those of SF60 (P < 0.01). However, these parameters showed no significant differences between SF40 and SF60 (Supplementary Figure S1).

The dominant orders of AMF communities among all samples were Glomerales (42.8 ± 3.8%; means ± SD), followed by Paraglomerales (29.1 ± 5.1%), Archaeosporales (0.8 ± 0.3%), and Diversisporales (0.8 ± 0.2%). Together, these accounted for more than 73.5% of the total AMF sequences in all the samples (Supplementary Table S1). The relative abundance of Paraglomerales was significantly increased along the successional stages (Figure 3). Interestingly, the lowest relative abundances of Glomerales and Archaeosporales were found in SF40. However, these relative abundances had no significant differences among the successional stages (Figure 3; Supplementary Table S1). These findings indicated a partial recovery of AMF communities at the order levels (Paraglomerales increasing, while Glomerales and Archaeosporales were decreasing) through the 40a recovery after deforestation disturbance.

The AMF were mainly distributed in the orders Glomerales (42.8% of the total AMF reads) and Paraglomerales (29.1% of the total AMF reads). The relative abundances of Glomerales were no significant differences among the successional stages (Figure 3). In addition, no significant changes were observed at the family or lower levels, except for OTU13, OTU15, OTU23, and OTU28 (Supplementary Table S1).

SOC, SN, carbon/nitrogen ratio (C/N), and pH had significant differences among the successional stages (Supplementary Table S2). The species richness of common species, rare species, and all species of AMF exhibited a significantly positive relationship with SN (Figure 4). However, species richness of all species of AMF, species richness, Shannon H', Simpson, and Pielou J' indices of rare species of AMF showed a significantly negative relationship with C/N (Figure 4).

The LEfSe of AMF at successional stages were OTU10 for G; OTU6, OTU9, and OTU141 for S; OTU3 and OTU38 for SF40; and OTU2, OTU15, OTU23, and OTU197 for SF60, respectively (Figure 5).

The most dominant orders of AMF communities among successional stages were Glomerales and Paraglomerales. The distribution of OTUs of G was similar to that of S, whereas this was in contrast to those of SF40 and SF60 (Figure 6).

Significant differences were observed in AMF communities between G and S, as well as between SF40 and SF60. Previous studies demonstrated that changes in fungal composition occurred during field successions, where late-successional plant species were highly dependent on AMF as compared to the early-successional plant species, which were not observed in the study of Koziol and Bever (2015). As many the mid- and late-successional plant species were shown to have a high overall beneficial response to AMF, these species were also likely to be more sensitive to the presence of particular beneficial AMF in their environment (Maitra et al., 2021). The changes in AMF composition could mediate plant species turnover during successional periods, which could be particularly important for the establishment and subsequent dynamics of late-successional plants (Koziol and Bever, 2016). AMF have played an important role in secondary forest successions (Heilmann-Clausen et al., 2015). Differences in AMF species composition along the successional stages were also found in the eastern periphery of Amazonia (Reyes et al., 2019), which confirmed the results of Pereira et al. (2014) in Atlantic rainforest areas of Northeast Brazil and those of Leal et al. (2013) in western Amazonia.

In this study, the SOC, SN, C/N, and pH showed significant differences with succession, and the species richness of common species, rare species, and all species of AMF indicated a significantly positive relationship with SN. Whereas, the species richness of all species of AMF showed a significantly negative relationship with C/N. Some previous studies reported that SOC and SN were increased along with proceeding successions (Allison et al., 2007; Eisenhauer et al., 2010). In addition, C/N in S was significantly lower than those of other successional stages (Xue et al., 2017). Both SOC and SN were increased along with the succession from G to other stages, which stimulated microbial metabolisms and changed bacterial and fungal community structure by increasing soil amino sugars (Shao et al., 2017). The C/N could significantly affect the structure of the soil microbial community and the nitrogen supply of plants (Waring et al., 2013).

The soil nutrient availability had been proven to be the key determinant that regulated the assembly of AMF communities (Johnson et al., 2015; Jiang et al., 2018) and vice versa (Rillig et al., 2015). The SN could directly drive the change in the AMF community in mixed roots (Shi et al., 2022). The change in SN indirectly regulated the species composition of the AMF community by driving plant species richness at both the whole-community scale and the individual plant species scale (Shi et al., 2022).