Ganshiling Natural Reserve served as the site of the study (109°37’∼109°41′E, 18°19’∼18°24′N), a typical tropical lowland rainforest from 50 to 681 m, and the slope is <50°. It is a typical lowland hilly landform with a tropical marine monsoon climate on Hainan Island, China. The annual average precipitation and temperature are 1800 mm and 25.4°C, respectively. The soil parent material was granite and the main 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).

Three succession stages in Ganshiling Natural Reserve were studied: 1) the shrub forest (SC; 109°39’∼109°41′E, 18°22’∼18°23′N), was subjected to logging and formed by artificial disturbance such as construction of nature reserve, and the main plant species include: Aporosa dioica and Olea hainanensis; 2) the 40 years recovery secondary forest after logging ceased in 1979 (SF40; 109°39’∼109°40′E, 18°21’∼18°23′N), and the main plant species include: Garcinia oblongifolia and Vatica mangachapoi; and 3) the 60 years recovery secondary forest after logging ceased in 1959 (SF60; 109°39’∼109°40′E, 18°22’∼18°23′N), and the main plant species include: Vatica mangachapoi and Hopea reticulata. Detailed information of each sample were shown in Supplementary Table S1.

In August 2019, five soil samples of 0–10 cm depth were obtained in each plot from three plots of 900 m² (30 m × 30 m) in each succession stage. These samples were then divided into portions intended for AMF community (stored at −80°C) and soil property assessment. The soil moisture was measured by gravimetric method (Blake, 1965). A 2 mm screen sieve was used to sieve air-dried soil samples for analyzing soil pH, soil organic carbon (SOC), total nitrogen (TN), and total phosphorus (TP), as well as soil stoichiometric characteristics (e.g., C/N ratio, C/P ratio, and N/P ratio). In order to measure soil pH, soil suspensions were mixed with distilled water. An analysis of the soil TN and TP was conducted using sulfuric acid and potassium sulfate after digestion with sulfuric acid, and then followed with the protocol described by Kjeldahl (Bremner, 1960; Taylor, 2000).

At the same time, vegetation surveys and root samples were conducted on sample plots to obtain aboveground plant community characteristics. Using gravimetric method to measure root biomass. And the vegetation survey data were used to calculate the Shannon-Wiener index, Margalef index and Pielous index.

Using the MIO-BIO PowerSoil DNA Isolation Kit (MO BIO Laboratories, Carlsbad, CA, United States), soil DNA was extracted from 250 mg soil samples. The NanoDrop ND-2000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) was used to analyze and quantify DNA, which was prepared at a standard concentration for PCR amplification.

For PCR amplification, AMV4.5NF (5′-AAG​CTC​GTA​GTT​GAA​TTT​CG-3′) and AMDGR (5′-CCC​AAC​TAT​CCC​TAT​TAA​TCA​T-3′) were used as primers (Lumini et al., 2010). An AMF library was constructed using two-step PCR amplification. Target fragments were amplified using specific primers. Later, an Illumina adapter sequence, a primer pad, and a barcode were added into both the forward and reverse target fragments. A 50 ml reaction volume was used for PCR of each sample, comprising 10 ml of 5x buffer, 1 ml of dNTP (10 m/M), 1 U of Phusion DNA High-Fidelity polymerase (New England Biolabs, Ipswich, MA, United States), 1 ml of forward and reverse phasing primers, and 5–50 ng of target sequencing. The reaction volume was supplemented to 50 μl with RNase free ultrapure water. An ABI9700 PCR amplifier (Applied Biosystems, Waltham, MA, United States) was used to amplify AMF. Amplification products were kept warm at 10°C. With the assistance of an AxyPrepDNA Gel Recovery Kit (Axygen Scientific, Corning, NY, United States), electrophoresis was used on an agarose gel to pool and purify PCR products. FTC 3000TM real-time PCR (Fungyln Biotech, Richmond Hill, Ontario, Canada) was used to quantify purified DNA and mixed to an equal molar ratio to perform the second PCR amplification. When the second PCR amplification, the reaction volume change to 40 μl: 5x buffer (8 μl), dNTP (1 μl, 10 m/M), Phusion High-Fidelity polymerase (0.8 U, New England Biolabs, Ipswich, MA, United States), both forward and reverse phasing primer (1 μl), respectively, and DNA template (5 μl) and RNase free ultrapure water (5 μl).

Using different barcodes to identity sequences of different samples in the parallel sequencing test, the sequence with a complete barcode was regarded as a valid sequence and extracted to a dataset stored in the FASTQ format. Quality assurance and splicing was performed on the valid sequences, the specific process for which follows: 1) using the Trimmomatic program, poorly overlapped and low-quality sequences (quality score <20, window size <5) were removed (Bolger et al., 2014); 2) a merged pair of reads as a sequence was used to refer to the overlap relationship of PEreads (>10 bp) by the FLASH program (Magoc and Salzberg, 2011); and 3) merge sequences were screened according to error ratio (<0.2). By using mothur V.1.39.5, quality control of merge reads was performed by removing sequences that were ambiguous, homologous, or chimeras generated during PCR amplification (Schloss et al., 2009). Analysis of species information and operational taxonomic units (OTUs) clustering was based on the clean tags. OTUs were selected using UPARSE with a cutoff of 97% similarity (Edgar, 2013). The dataset was screened for chimeric reads using UCHIME in reference database mode (Edgar et al., 2011). Through USEARCH_global, all sequences were compared to OTU representative sequences to obtain the OTU table. The identity of OTU taxonomy was determined by using Blast to compared representative sequences with the Maarjam database (Zhang et al., 2000; Opik et al., 2010). The taxonomy of OTU was obtained for species which had the highest similar identity with OTU representative sequences in the Maarjam database, except for uncultivated environmental information. Some OTUs without species information were marked ‘Unidentified’. The library preparation and sequencing were performed by TinyGene Bio-Tech (Shanghai, China). Raw sequence data has been submitted to NCBI’s Sequence Read Archives (accession no. PRJNA884298).

In order to explore the AMF community dynamics, the change in distribution and composition of AMF were analyzed. Furthermore, the impact of environmental factor on AMF community dynamics were explored.

For bioinformation analysis, some OTUs that the classification results could not identify specific genus, and <5 sequences were not used. The species diversity was calculated in R package “vegan”. The non-metric multidimensional (NMDS) and principal coordinate analysis (PCoA) were analyzed using Bray-Crutis dissimilarity with R. The Venn diagram was used to demonstrate the distribution of OTUs during succession stages by R version 4.1.2.

Based on this dataset, the genus level was used to explore AMF composition differences which had the clearest taxonomy annotations. To explored the AMF community difference across groups, the OTU level was used for analysis. The top 10 relative abundances under genera and OTUs were selected for multi-group differential analysis.

Using R package ‘microbiomeMarker’, the linear discriminant analysis effect size (LEfSe) and linear discriminant analysis (LDA) were performed (Segata et al., 2011). Screening the OTUs which it frequency of occurrence >70% in samples for the following calculation. According to Spearman’s correlation coefficient (p-value ≤0.05), the dominant OTUs was selected to construct microbial network. The microbial network was using selected dominant OTUs as nodes and its strong correlations as edges. The screened of dominant OTUs and visualize data preparation were completed by R, and the microbial network visualized with Gephi (Sirova et al., 2018; Nauta et al., 2020).

By using R version 4.1.2, the mantel test, two-factor correlation heatmap and distance-based redundancy analysis (db-RDA) was analyzed. The significances of environmental factors with AMF were tested using Monte Carlo permutation ANOVA (Tamminen et al., 2018).

Among alpha diversity indices, the Chao1, ACE, and phylogenetic diversity (PD) of the soil AMF community showed significant differences between shrubland and secondary forest (p < 0.05). These indices showed an increasing trend with succession. Compared with other succession stages, SF60 group had the highest Chao1, ACE, and PD (Table 1). Conversely, Simpson and Shannon-Wiener indices showed no significant differences across succession stages (Table 1).

NMDS, PCoA, Venn, differential, and beta diversity analyses revealed significant differences in soil AMF communities between shrubland and secondary forest, but not between SF40 and SF60 (Figure 1A; Supplementary Figures S2, S3).

The Venn diagram showed that there were 270 shared OTUs in different groups, and the number of shared OTUs between SC and SF40 and between SC and SF60 were 46 and 39, respectively. However, the number of shared OTUs between SF40 and SF60 was 156 (Figure 1B). Kuklospora and Redeckera were just present in SC, and Diversispora were only present in secondary forests, respectively. In summary, AMF community distribution was influenced by succession, and the differences between shrubland and secondary forest were especially obvious.

As succession stages progressed, the composition of soil AMF community changed. As a result, Paraglomus and Glomus emerged as dominant genera, with Paraglomus exhibiting higher relative abundances in SF40. In contrast, the higher relative abundance of Glomus and Claroideoglomus was present in SC and SF60 (Figure 2A).

At the genera level, the multi-group differential analysis showed the variation of Paraglomus and Glomus were non-significant during succession. In addition, Ambispora and Rhizophagus were differential genera in succession, and Ambispora was significantly less abundant in shrubland than in secondary forest. In addition, Rhizophagus had the highest relative abundance in shrubland, and showed a decreasing trend with succession (Figure 2B).

At the OTUs level, the multi-group differential analysis showed that most of the Top 10 OTUs belonged to Paraglomus and Glomus, and only OTU15 belonged to Claroideoglomus. Multi-group differential analysis showed that the relative abundance of OTU1 and OTU4 increased along succession stages. The relative abundances of OTU5, OTU6, and OTU72 were greater in SC than in SF40 and SF60. OTU1, OTU2 and OTU3 were higher in SF40 than in SC and SF60 (Figure 3).

LEfSe analysis showed that Rhizophagus and Geosiphon were specific genus in SC, but Ambispora and Diversispora were specific genus in SF60 (Figure 4).

The correlation network analysis showed AMF community co-occurrence patterns were significantly different between succession stages. Although the nodes were similar among the three groups, the edges were significantly higher in SF40 (nodes = 54, edges = 272) than in SC (nodes = 56, edges = 145) and SF60 (nodes = 53, edges = 170). Furthermore, the SF40 network had higher network topological properties (average degree, network diameter, and clustering coefficient) than SC and SF60. Although most of modes in these networks were similar, the interaction between nodes were significantly different among groups (Figure 5).

In the Mantel test, aboveground plants community characteristics were the major factor in varieties of soil AMF community composition (r ² = 0.608, p = 0.004; Table 2). Likewise, the db-RDA analyses also showed that the Margalef index, Shannon-Wiener index, Pielou index, and root biomass of aboveground vegetation were drivers in changes of AMF community composition (Figure 6). In addition, the soil moisture, SOC, and total N—which belong to abiotic factors—also affect AMF community composition significantly (Supplementary Table S2). The composition of AMF communities could be explained by 54.2% of these factors. Among them, root biomass was the key factors which influenced AMF community composition (r ² = 0.789, p = 0.012; Figure 6; Supplementary Table S4).

The two-factor correlation heatmap of network attributes and environmental factors showed that C/N ratio significantly correlated to Harmonic geodesic distance (HD), network transitivity (Trans), and average clustering coefficient (avgCC), while root biomass (RB) and Margalef significantly correlated with total nodes (Figure 7).

With succession, the dominant vegetation gradually evolves from shrubs to trees, and it showed an increase trend in vegetation community complexity. Consistent with the variation of aboveground vegetation communities, AMF species richness and spectral diversity were increased during succession stages in this study (Table 1), and it could be confirmed that the AMF is closely related to host plants, which have host preferences (Janos, 1980; Xu T. et al., 2016; Li et al., 2020). It could be further found that the composition of AMF communities is highly similar in higher taxonomy levels during restoration succession, and it contrasted to vegetation communities composition had significant differences with succession (Leal et al., 2013; Reyes et al., 2019). For example, AMF community composition was similar at the genera level during succession in this study (Figure 2B), and the differences in AMF community composition were significant at the OTU level (Figure 3). In this study, Paraglomus and Glomus dominated in AMF community during succession, which benefits from its biological characteristics that have wide biogeographic distributions and vegetation cover (Figure 2A) (Lovelock et al., 2003). This study demonstrates that the importance of OTUs varied significantly between groups. This dynamic is associated with changes in the AMF niche and functionality at the OTUs level and it confirms that the host preference of AMF community was reflected in changes to the species importance and many coexisting microorganisms have functional redundancy at the OTU level (Louca et al., 2018). The niche analysis of keystone microbial taxa at different successional stages proof that the niche of these keystone taxa at each stage was clearly differentiated (Supplementary Figure S5). In generally, AMF functional redundancy and niche changes are likely response mechanisms to help AMF adapt the habitat changes during succession which consistent with previous studies (Varela-Cervero et al., 2015; Alguacil et al., 2016).

This study showed solid correlation between plant communities and soil AMF communities (r ² = 0.608, p = 0.001; Table 2). As an important factor, the host plants could affect the AMF community dynamic by altering the microhabitat (McGuire et al., 2012; Chai et al., 2019). With the change in vegetation community, the host plant reshapes the AMF community composition by secreting apoplast chemistry and root exudates (Goldmann et al., 2015; Hanif et al., 2019; Li et al., 2020). This effect further impacts the uptake of nutrients and water by vegetation (McGuire et al., 2012; Nakayama et al., 2019). In this study, the root biomass also has the strongest correlation with soil AMF community (r ² = 0.789, p = 0.012; Figure 6; Supplementary Table S4) because microbial activity and growth could be stimulated by roots (Wang et al., 2017). By altering root turnover and/or root exudation, the change in root biomass may affect resource availability for the soil AMF community (Philippot et al., 2013). As well as influencing soil microbial community nutrient uptake, root exudates act as signaling molecules and attractants, as well as inhibitors or repellents (Baetz and Martinoia, 2014). Therefore, the soil AMF community might be regulated by changes in root biomass during succession (Eisenhauer et al., 2017).

In addition, soil properties and AMF community composition correlated strongly (Figure 6). Decomposition rates of organic matter and nutrient cycling were directly influenced by soil properties, which contributed to a diverse community of soil AMF (Xu X. et al., 2016). The alternations in AMF community also react upon nutrient cycling (Phillips et al., 2013). As an important environmental factor, the soil moisture drive for soil AMF community characteristics (r ² = 0.689, p = 0.036; Figure 6; Supplementary Table S4). There is a direct relationship between soil moisture and nutrient effectiveness, which could impact on the soil microbial diversity and abundance through affect microbial nutrient uptake, and it could further affect in soil AMF community composition (Deepika and Kothamasi, 2015; Peng et al., 2021). In addition, soil moisture had clear legacy effects on soil properties, the nutrient uptake of plants, the formation of AMF, and mycorrhizal responsiveness, which would alter the synergistic relationship between AMF and host plants (Cavagnaro, 2016).

This study’s results showed that the most complex network occurs in SF40 (Figure 5). In the Ganshiling Nature Reserve, the aboveground vegetation suffered extensive deforestation in the 1980s and gradually recovered to become a tropical lowland rainforest after successive anthropogenic disturbances. With increasing restoration years, the intensity of disturbance to the successional stages which were selected for this study varied. With further explored, it could be found that SF40 suffer intermediate disturbance frequency compared with other. According to intermediate-disturbance hypothesis (IDH), the SF40 with intermediate disturbance frequency could keep the most active AMF community. Moderate ecological scale disturbances keep off interspecific competition and provide supportive environment for microbial species which about both early and late succession to survival (Grime, 1973; George et al., 2019).

Moreover, consistent with existing research, the two-factor correlation heatmap showed that changes in AMF community structure were closely related to soil C/N ratios (Figure 7; Wan et al., 2014). It is well established that the Trans and avgCC positively correlated with AMF community activity and HD showed negative correlation with it. In general, a significant positive correlation between change in soil C/N ratio and AMF community activity as succession progressed. Soil C/N ratio responds to soil nutrient status by influencing nutrient mineralization and organic matter decomposition. It has been shown that the soil C/N ratio affects microbial biomass and activity by altering the nutrient effectiveness and quality of the soil substrate in which microorganisms survive (Demoling et al., 2007).