The study site was located in the Hongshixia Sandy Botanical Garden of Yulin City, Shaanxi Province (38°26′N, 109°12′E, 1080 m a.s.l.), which is on the southern edge of the Mu Us Desert (Supplementary Figure S1). It is characterized by the temperate semi-arid continental monsoon climate with an annual temperature, precipitation, and evaporation of 9.1°C, 385 mm, and 2914 mm, respectively. The frost-free period ranges from 160 to 170 days. The soil is mainly composed of aeolian sandy soil, which is loose and very susceptible to flow and erosion. The dominated arbor plants are P. sylvestris and P. tabuliformis. The shrub and herb include Caragana korshinskii, Lespedeza bicolor, Cares duriuscula, Rhamnus parvifolia, Salsola collina, and Potentilla anserine.

Field investigation and sample collection were conducted in August 2017, which is the peak of plant growth season and microbial activity period. We randomly set a 50 m × 50 m plot in P. sylvestris plantations across 26, 32, and 43 a, which represent half-mature (MUh), near-mature (MUn), and mature (MUm) forests, respectively (Table 1). There were no tending measures and visible disturbance by anthropogenic activities since the establishment of plantations.

At each plot, five healthy and widely spaced standard trees were selected as sampling objects. Loose litters, herbs, and the undergrowth humus layer were gently removed prior to digging out the fine root samples carefully along the base of the trunk. After that, three repeated fine samples were taken from each standard tree; then, the three samples were fully mixed. A total of 15 mixed root samples (3 plots × 5 mixed samples = 15 composite samples) were used for DNA extraction and molecular identification. The ectomycorrhizal infected roots were accompanied with soil placed into plastic bags to keep the root tips fresh. Root samples were kept cool after sampling and stored at 4°C. Soil samples were collected around the root, with three repeated soil samples collected in each point, and fully mixed into one (3 plots × 5 mixed samples = 15 composite samples). Soil samples from 0 to 20 cm depth were obtained to determine soil physico-chemical properties and enzyme activities. After passing through a 2 mm sieve, the soil samples were divided into two parts; one was stored at 4°C in a refrigerator to determine enzyme activities and the content of available nutrients in soils, and the other was taken straight back to determine soil physical properties, soil organic matter (SOM), and soil total nutrient.

Soil water content (SWC) was determined by the gravimetric method after drying the soil sample in an oven at 105°C for 12 h. Soil pH was measured using a pH meter (PHS-3E, INESA, Shanghai, China); the mixture of soil and water is 1:2.5. Soil organic matter (SOM) content was measured by the K₂Cr₂O₇–H₂SO₄ oxidation technique (Walkley and Black, 1934). Soil total nitrogen (TN), total phosphorus (TP), nitrate nitrogen (NO₃ ⁻-N), and ammonium nitrogen (NH₄ ⁺-N) were determined using an automated element analyzer (Smartchem 200, Italy). TN used the indophenol-blue spectrophotometric method, and TP used the Mo-Sh anti-colorimetric analysis method (John and Matt, 1970; Mason et al., 1999). Soil invertase (INV) activity, urease (URE) activity, and phosphatase (PHO) activity were assayed as described by Tabatabai (1994).

DNA was extracted and purified using the Powersoil® DNA Isolation Kit (Mo Bio Laboratories, Carlsbad, CA, United States), according to the manufacturer's protocol. The genomic DNA concentrations were checked on 1% agarose gels electrophoresis and quantified on a NanoDrop NC2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States). The fungal internal transcribed spacer (ITS) region was amplified on Eppendorf Mastercycler Gradient Thermocycler (Germany); for the ITS1 regions of fungi, the primers were ITS1F (5′-GGA​AGT​AAA​AGT​CGT​AAC​AAG​G-3′) and ITS2 (5′-TCC​TCC​GCT​TAT​TGA​TAT​GC-3′) (Amend et al., 2010). For each sample, an eight-digit barcode sequence was added to the 5′ end of the forward and reverse primers (provided by Allwegene Company, Beijing). The polymerase chain reaction (PCR) mixtures were analyzed on a Mastercycler Gradient (Eppendorf, Germany) using 25 μL reaction volumes, containing 2.5 μL 2 × Taq PCR MasterMix, 3 μL BSA (2 ng/μL), 1 μL Forward Primer (5 μM), 1 μL Reverse Primer (5 μM), 2 μL template DNA, and 5.5 μL ddH₂O. The cycling parameters were 95°C for 5 min, followed by 28 cycles of 95°C for 45 s, 55°C for 50 s, and 72°C for 45 s with a final extension at 72°C for 10 min. The PCR products were purified using an Agencourt AMPure XP Kit. Deep sequencing was performed on an Illumina Miseq platform at Allwegene Company (China, Beijing). After the run, image analysis, base calling, and error estimation were performed using Illumina Analysis Pipeline Version 2.6.

The raw data were first screened, and then the sequences shorter than 230 bp and with low-quality scores (≤20) were removed. After removing non-fungal sequence reads, the fungal sequences were clustered into operational taxonomic units (OTUs) (Blaxter et al., 2005) at a 97% similarity level using the Uparse algorithm of Vsearch (Version 2.7.1) (Edgar, 2010). The remaining high-quality sequences were queried against the GenBank (https://academic.oup.com/nar/article/50/D1/D161/6447240) non-redundant nucleotide database in NCBI using local BLASTn (https://blast.ncbi.nlm.nih.gov/). The MEGAN (Version 6.0) program was used to assign BLAST hits to taxa of the NCBI taxonomy. Non-EcM fungal OTUs were filtered out by FUNGuild (Version 1.0, http://www.stbates.org/guilds/app.php) (Zhang et al., 2015). The key EcM fungal species were analyzed by using the analysis of modular topological roles (within-module connectivity (Zi) ≥ and among-module connectivity (Pi) ≥ 0.62). All the identified sequences were submitted to the GenBank with the accession numbers of SUB7166485 (from MT229482 to MT229553) (Table 3).

Define “dominant species” as the EcM fungi detected in 80% of all samples, “frequent species” as the EcM fungi detected in 30% of all samples, and “rare species” as the EcM fungi detected in less than 30% of all samples. The important value of OTUs is represented by the sum of its relative abundance and relative abundance. α-diversity indices were calculated based on the OTU dates, including the number of identified genera (S), Shannon index (H′), Simpson’s index (d), and evenness index of species richness [Pielou (J)]. The differences in community structure of the EcM fungi were analyzed by non-metric multidimensional scaling (NMDS) analysis. Bray–Curtis dissimilarity index was calculated to examine the dissimilarity of the OTU community in different sites. Dissimilarity matrices of environmental data sets were calculated based on Euclidean distance. A redundancy analysis (RDA) was used to determine the correlation between the EcM fungal community (genus level) and soil properties. Subsequently, Spearman rank correlation was used to explain the relationship between the EcM fungal community (genus level) and soil properties, also based on Bray–Curtis distance. The correction heatmap was a visual correlation between all soil properties and the five most frequently recovered genera. After converting soil properties to a matrix, we performed the Mantel test to determine the significance of each environmental variable. The EcM fungal network was constructed based on random matrix theory using the OTU data, which present in at least eight (>50%) samples (Deng et al., 2012). All network analyses were performed using the MENA platform (http://ieg4.rccc.ou.edu/mena/main.cgi).

One-way analysis of variance (ANOVA) with the least-significant difference (LSD) and the Tukey test were performed to determine the significant difference among soil variables and EcM fungal diversity. All statistical analyses were performed in SPSS (Version 20.0), and p < 0.05 was considered statistically significant. The α-diversity indices, Bray–Curtis dissimilarity index, network analysis, correction heatmap, and Mantel test were implemented using the R (Version 3.5.1) program. RDA was performed using CANOCO (Version 4.5). The network diagram was visualized by Cytoscape (Version 3.9.0).

In P. sylvestris plantation, all the soil properties varied by stand age (p < 0.05, Table 3). Soil TN, SOM, NO³⁻, INV, and URE were significantly different with stand aging (p < 0.05). In general, soil in the Mu Us Desert was alkali to neutral. The maximum values of soil nutrient contents (SOM, TN, and TP) were found in half-mature forests. Soil URE and PHO enzymes were the most active in mature forests.

A total of 361,548 high-quality sequences were obtained from roots associated with P. sylvestris, assigned to 807 fungal OTUs. A total of 72 of these had an EcM lifestyle and were identified to the genus level and below (Table 2). With stand aging, the numbers of OTUs and rhiness of genus increased first and then decreased; the values of the remaining α-diversity indices gradually decreased. The impact of stand age on the EcM fungal community has not reached a significant level (p > 0.05, Figure 1).

Our studies found that the EcM fungal OTUs of P. sylvestris plantations in the Mu Us Desert were assigned to two phyla, four classes, 12 orders, 18 families, and 23 genera (Figure 2; Supplementary Table S1). Of all 72 EcM OTUs, Basidiomycota was the most abundant phylum (83%) with 63 OTUs, following by Ascomycota (17%) with nine OTUs. Among them, 13 and 14 species were detected for the dominant and common species, respectively. The rare species was the largest pool (45 species) accounting for 62.50% of the total species. At the genus level, Rhizopogon was the most dominant EcM fungus; it has the highest relative abundance in all age groups, especially in MUn and MUm (more than 50%, Figure 2). Besides, the genera of Tuber, Tomentella, and Inocybe were also widely distributed in all plantations. With the aging of stand, Rhizopogon increased first and then decreased, while Tomentella and Inocybe were just the opposite. Tuber was accumulated with stand age.

The Bray–Curtis distances of the EcM fungal community composition in three stand ages are shown and ordered by NMDS. The stress value is 0.1679 (<0.2), so the Bray–Curtis distance can better reflect the difference of community composition (Table 3). The coordinates of scatter points indicated that the EcM fungal composition varied in different stand ages (Figure 3). ANOSIM results confirmed that the difference of fungal community structure was obvious with stand age (R = 0.88, p > 0.05). Based on the area of the confidence ellipse and Bray–Curtis distances (Figure 3; Table 3), the EcM fungal community structures of the near-mature and mature forests were the most similar.

The molecular ecological network of the EcM fungi reflects the coexistence relationship of species in the environment (Figure 4). After removing the OTUs with average relative abundances less than 0.01%, a network that included 57 nodes and 435 edges was constructed. Of all these edges, 206 (47.36%) had shown positive corrections and 229 (52.64%) had shown negative correlations. In the EcM fungal network, the average degree, average clustering coefficient, and average path distance were 15.263, 0.571, and 1.875, respectively. The modularity calculation of the EcM network divides the whole network into four individual models (Supplementary Figure S2). OTU306 (classified as Rhizopogon mohelnensis), OTU43 (uncultured Peziza, Peziza sp.), OTU40 (classified as Clitopilus sp.), and OTU14 (classified as Clitopilus sp.) of Basidiomycota and OTU114 (classified as Rhizopogon mohelnensis), OTU56 (classified as Wilcoxina mikolae), and OTU36 (classified as Geopora arenicola) of Ascomycota had weak relationship with other species. We identified key EcM fungal species by using the analysis of modular topological roles. According to within-module connectivity (Zi) and among-module connectivity (Pi), OTU60 (Zi = 2.60 > 2.5, Pi = 0.625 > 0.6, classified as Hebeloma collariatum) was the core species of the model (Supplementary Figure S3).

The EcM fungal community composition through stand ages was driven by soil properties (Figures 5, 6; Supplementary Table S3). The redundancy (RDA) analysis revealed the impact of soil characteristics on the EcM fungal community, and the contribution rates of principal components 1 and 2 were 25.78 and 18.49%, respectively. All four axes in total contributed 67.3% of the explanation (Figure 5). Among all soil properties, NO₄ ⁺ content had the highest contribution rate (15.32%). The contributions of SOC, TN, INV, and stand age all exceeded 11%.

The correction heatmap and Mantel test between soil factors and major EcM genera were created by Spearman’s rank correlation (Figure 6; Supplementary Tables S2, S3). The INV had a significant positive effect on the EcM fungal community, while SOC, TP, and stand age have a significant negative effect on the major EcM fungi (p < 0.05 or p < 0.01). URE and NO₃ ⁻ had an excellent negative correlation to the EcM fungal community (p < 0.05), although their contribution rates were low.

EcM fungi are abundant and diverse in temperate and boreal forests (Steidinger et al., 2020), and the succession and variation of the underground EcM fungal community was a complex ecological process (Bardgett and van der Putten, 2014; Franklin et al., 2014). The diversity indices and community structure of EcM fungi changed with stand aging in the Mu Us Desert, but the effect of forest age was not significant. This was consistent with previous studies that the diversity of the fungal community was not correlated with the aging stand (Wang et al., 2020). As expected, the EcM fungal community was dominated by Basidiomycota rather than Ascomycetes, which is consistent with previous results (Geml et al., 2017). It is known that most of the known EcM fungi are Basidiomycota, which are diverse and adaptive to various environments (Boeraeve et al., 2018). Therefore, Basidiomycota considerably contributes to the ecological restoration in deserts.

The endemic genera were found across different stands, and the EcM fungal community structure was significantly different between three age stands. The distribution pattern of the EcM fungi indicates the succession of EcM fungi during forest development and the influence of stand age on it. The variations in plant community and soil characteristics caused by stand age altered the establishment and distribution of the EcM fungal community (Hayward et al., 2015). Rhizopogon was the dominant EcM fungal genera of P. sylvestris across all stand ages in the Mu Us Desert. Besides, Tomentella, Inocybe, and Tuber were also widely distributed in P. sylvestris plantations. These OTU-richest genera are also commonly detected in boreal coniferous forests (van Dorp et al., 2016). They are easily dispersed and widely adapted to the harsh environment (Ryberg et al., 2008; Bahram et al., 2011). Obviously, they are “multistage” fungi with an excellent competitive ability (Taylor and Bruns, 1999).

As the largest genus of hypogenous Basidiomycota, Rhizopogon were widely distributed and associated with pines; it is also the most common EcM fungi of P. sylvestris plantations in the Mu Us Desert. Especially, Rhizopogon proliferated in the high-carbon near-mature forest due to the high demand for carbohydrates (Mujic et al., 2014). Further, R. rubescens and R. mohelnensis contribute to the stress resistance; therefore, they conduce to the improvement in survival of conifer seedlings in a desert environment (Bahram et al., 2011). As saprotrophic EcM fungi, Tomentella prospers the nutrient cycling and is widely distributed throughout the Northern hemisphere forests (Han et al., 2017). Besides, Tomentella is also of great resistance to environmental stresses, supported by the presence in the degraded oak forest (Mosca et al., 2007). Inocybe is ecologically important in the forest succession by providing water and minerals for the host plants (Mason et al., 1983). In the Mu Us Desert, Inocybe were clearly decreasing from the half-mature forest to a mature forest. Soil moisture consumption is greatly responsible for this instance due to their preference of an arid environment, which has obvious significance in the desertification region. Furthermore, Cortinarius, Amphinema, and Russula were widely distributed in the mature forest. They are typical “late-stage” fungi, which are usually found associated with older forests and have strong competitive abilities (Hayward et al., 2015; Guo M. et al., 2020). The core fungi of the EcM fungal network are classified as Hebeloma, which is a frequent EcM genus, and its species is usually associated with a wide range of Pinaceae trees (Gryta et al., 2010).

Suillus is regarded as the most generic EcM fungi associated with pine, but it was diminished in the Mu Us Desert; this result is inconsistent with previous studies (Policelli et al., 2018). Geopora is considered to be an important mutual partner for host plants resisting the stress conditions. Notably, Russula and Peziza are also detected in the P. sylvestris plantations, which usually rely on fertile soil and prefer high pH, respectively. Lactarius, Hebeloma, and Lachnum accounted for a very small percentage; this was supported by the previous studies that they were not easily established in desertification regions (Desai et al., 2016). Curiously, some common drought-tolerant and barren-resistant EM fungal species were missed, such as Lycoperdon and Cenococcum (Jany et al., 2003).

The primary drivers of the EcM fungi community at different scales are not necessarily the same (Peay and Matheny, 2016). EcM fungi were strongly determined by the soil environment, and the edaphic factors are approved to be the major drivers of the EcM fungal community at a local scale (Miyamoto et al., 2015). Forest development is usually accompanied by variations in the soil system (Zhao et al., 2020); EcM fungi are highly sensitive to the environmental changes. Therefore, host-related variables strongly influence the EcM fungal community. The variations of soil properties with stand aging are of great importance to shape fungal communities (Truong et al., 2019). We found that most of the factors we measured in this study were proved to sufficiently affect the community structure of the EcM fungi. Among them, the soil enzyme activities of invertase and urease were the most major abiotic factors causing EcM fungal communities in cation regions. They were the most important enzymes for soil fertility and nitrogen metabolism. ECM fungi have different preferences and demands for nutrient strategies, and the nutrient competition in fungal communities alters the assembly of the EcM fungal communities (Yang et al., 2021). EcM fungi-induced soil nutrient cycling substantially contributes to the growth and health of host plants. Nitrogen is the foremost limiting nutrient in conifer forest systems (Sebastiana et al., 2018). In our study, N shift was clearly observed, that is, NH₄ ⁺ and NO₃ ⁻ concentrations significantly increased with aging stand. Meanwhile, EcM fungi (e.g., Peziza and Russlua) were negatively correlated with soil N availability. This finding confirms that EM fungal absorption capacity varies in response to the various forms of N, so N availability determines the community composition of EcM fungi (Kjøller et al., 2012). Furthermore, EcM fungi are good decomposers of soil organic matter and play a major role in the turnover and stabilization process of soil nutrients. The sufficient soil nutrients provide a suitable microenvironment for fungi, which is the basis of expanding and prospering the EcM fungal community.

Water shortage was considered to be an important abiotic factor; it was also the strongest abiotic candidate filter (Liu et al., 2021b). This result is consistent with previous studies that soil moisture is a key factor to alter the EcM fungi community composition in arid environments (Dean et al., 2015). Meanwhile, with the specificity of EcM fungi to drought stress, in return, variations in the EcM fungal community have direct impacts on the growth and development of host plants (Kennedy and Peay, 2007). Besides, considerable research studies have been devoted to come to the conclusion that soil pH is another key driver of the EcM fungal community (Kutszegi et al., 2015). The EcM fungal community changes significantly even with a very small variation of pH values. On one hand, soil pH affects the growth of the EcM fungi and mycelial metabolism. On the other hand, soil pH indirectly alters the EcM fungal community by changing nutrient availability (Bennett et al., 2017). Generally, EcM fungi prefer slightly acidic to neutral soils, and the growth of the EcM fungi is inhibited at a higher-pH environment. In the Mu Us Desert, the soil pH ranged from 7.27 to 7.58. The weakly alkaline soil leads to potential variations in the EcM fungal community.

EcM fungi can secrete nutrient-related enzymes to regulate the nutrient cycling. The active soil enzyme contributes to the N and P absorption and slows the decomposition of soil organic matter (Lindahl and Tunlid, 2014). Therefore, the interaction between EcM fungi and the soil system is affected by soil enzyme activity. Generally speaking, urease and phosphatase benefit the mineralization of nitrogen and phosphorus, thus providing available nitrogen and phosphorus for the hosts (Bidartondo et al., 2001). In this study, the urease activity is increasing with stand aging; it is beneficial to available nitrogen supply for host plants. This finding was also confirmed by the soil nitrogen in suit investigation. Invertase is related to the transformation of soil organic carbon, but with the increase in soil organic matter, invertase activities were significantly reduced. This reduction was also reported that only a few EcM fungi have saprophytic functions and are active in the decomposition of soil organic matter (Jonard et al., 2014).

EcM fungi and host plants are a community of shared common interests, and they get feedback from each other. The growth of P. sylvestris plantation is inseparable from the assistance of EcM fungi, and EcM fungi also benefit from the host plant (Zhu et al., 2006). Therefore, the interaction between EcM fungi and host plants is a potential but crucial indicator for the degradation of P. sylvestris plantation. Natural forests with diverse ages and hierarchical structures usually have abundant, diverse, and complex EcM fungal community compositions and structures (Grebenc et al., 2009). Compared to the original forest, artificial plantation is simple in structure and has poor stability, resulting in EcM fungi with low heterogeneity and diversity (Guo M. et al., 2020). Therefore, this simple EcM fungal community composition may be one of the reasons for the degradation of plantation. There is overwhelming evidence corroborating the notion that structure and diversity determine ecosystem function and stability (Huang et al., 2021).

As the most typical EcM fungi, Suillus has extreme host specificity of Pinaceae, especially P. sylvestris. Suillus forms symbionts with conifer species and substantially contributes to P. sylvestris invasion and seedling establishment by improving soil nutrient absorption (Policelli et al., 2018). Previous studies have demonstrated that Suillus is the dominant EcM fungal genus of P. sylvestris in nature forests (Truong et al., 2019). However, only one species of Suillus associated with P. sylvestris was found in the Mu Us Desert. Hence, we deduced that functional EcM fungi Suillus (e.g., Lycoperdon, Cenococcum) commonly found in nature forests contribute to the P. sylvestris plantation degradation. This assignment has explained the central importance of EcM fungi in P. sylvestris plantation degradation. In a desert environment, there are intensely interspecific and intraspecific competitive interactions for the limiting resources. Therefore, the simplification community structure and the loss of EcM fungal diversity and function eventually lead to the degradation of P. sylvestris plantation.