The present study was conducted in two typical areas of L. gmelinii forests in the Greater Khingan Mountains, China ( Figure 1 , Table 1 ). One study area was located in the Aershan National Forest Park, Aershan (47°15′17″N, 120°17′34″E, shortened to Aershan), the southernmost area of L. gmelinii growing in China. The other was located in the Jinlin Forest Farm, Genhe (51°21′51″N, 121°20′31″E, shortened to Genhe), which is the core area for the distribution of L. gmelinii. The study area is located in a high-latitude permafrost distribution area in China (Ding et al., 2019). The winters in the study areas were long and cold, the summers were short and rainy, and spring and autumn were cool with less precipitation (Liu et al., 2020). The dominant tree species in these areas were L. gmelinii and Betula platyphylla, with a forest coverage rate of 80%. The zonal soil at both sites is Retisol (IUSS Working Group, 2015). The soil in Aershan is a typical black soil that is soft, moist, and nutritious, whereas the soil in Genhe is rich in gravel. A shovel was used to turn over the stones and obtain the soil in the stone cracks.

The rhizosphere and bulk soil samples were collected in July and October 2021, representing summer and autumn, respectively. At each sampling site, we established three 20 m × 20 m plots at least 100 m apart and collected five soil samples from five points (four corners and the center) as independent samples. Roots were extracted by digging the trees near the points to collect rhizosphere soil. In total, 120 samples (5 independent replicates× 3 plots × 2 compartments × 2 sites × 2 seasons) were used for DNA extraction. Additionally, 60 bulk soil samples were used to determine the soil physicochemical properties. Specifically, bulk soil from 0 to –20 cm was obtained, sifted through a 2 mm mesh sieve, mixed, and then packed into centrifuge tubes and sealed bags. The roots were stored in sealed bags. The samples in centrifuge tubes were stored with dry ice during the sampling process and transported to the laboratory stored at −80°C. The roots were stored with ice in the field and then transported to the laboratory stored at 4°C for the rhizosphere soil. The bulk soil was stored in sealed bags and transported at room temperature.

Rhizosphere soil was then separated from roots by shaking root and soil in 25 mL 0.1 M sterile phosphate buffer (7.1 g Na₂HPO₄ and 4.4 g NaH₂PO₄·H₂O were added to 820 mL deionized water; pH 7.0) in a 50 mL plastic conical centrifuge tube at 67 g for 5 min. The soil suspension was centrifuged at 4000 g for 15 min, and the obtained soil pellet was regarded as the rhizosphere soil and stored at −80°C until DNA extraction (Liu et al., 2021).

To determine soil properties, we followed the method described by Bao (2000). The soil in the sealed bags was dried naturally to determine the following physicochemical properties: soil organic carbon (SOC) was determined by the potassium dichromate volumetric method (Li et al., 2020), available phosphorus (AP) was determined by molybdenum antimony resistance colorimetry (Cai et al., 2018), pH was determined by the potential method (Li et al., 2020), total nitrogen (TN) was determined by Kjeldahl method (Shi et al., 2014), and the cation exchange capacity (CEC) was determined by the ammonium acetate exchange method (Liang et al., 2015).

Historical monthly weather data were obtained from WorldClim (https://www.worldclim.org/data/monthlywth.html#; accessed July 20, 2023), including average minimum temperature (°C), average maximum temperature (°C), and total precipitation (mm) in July and October 2021. The mean annual temperature and precipitation for 1970–2000 were obtained from WorldClim.

DNA was extracted using the DNeasy PowerSoil Pro Kit (Qiagen, Frankfurt, Germany). A NanoDrop NC2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) was used to quantify the DNA, and the quality of the DNA was detected using 1.2% agarose gel electrophoresis. The ITS1 region of fungi was amplified, and the primers were ITS5F (5-GGAAGTAAAGTAAAAGTCGTAAAGG-3) and ITS2R (5-GCTGCGTTCTTCATCGATGC-3) (White et al., 1990). PCR system (20 μL): 2 μL (2.5 mM) of dNTP. Of note, 1 μL (10 µM) forward primers and 1 μL (10 µM) reverse primers, 2 μL template DNA, 10 μL ddH2O, and 4 μL Fast pfu DNA polymerases. Circulating system: 95°C for 2 min; 30 cycles at 95°C for 30 seconds, 55°C for 30 seconds, and 72°C for 30 seconds; 72°C for 5 min. PCR amplification was performed using an Applied Biosystems 2720 Thermal Cycler (Thermo Fisher Scientific, Waltham, MA, USA). The Quant-iT PicoGreen dsDNA assay was performed using a Microplate Reader FLx800 (BioTek, Burlington, Vermont, USA). PCR amplicons were purified and recovered by adding Vazyme VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China), quantified with the fluorescent reagent of Quant-iT PicoGreen dsDNA assay kit (Invitrogen, Carlsbad, CA, USA), and then mixed in proportion to the sequencing amount. Pair-end 2 × 250 bp sequencing was performed using the Illumina MiSeq platform with a NovaSeq 6000 SP reagent kit at Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China). All raw sequencing data from this study have been deposited in the NCBI database with the Short Read Archive (SRA) accession number PRJNA1037725.

The early sequence processing in this experiment was based on QIIME 2–2021.2 (Bolyen et al., 2019). A Demux plug-in was used to split the samples. The DADA2 plug-in was used to perform quality control, such as filtering and noise removal; the positions at which the forward and reverse read sequences should be truncated due to a decrease in quality were 222 and 232, respectively, and the minimum and maximum sequence lengths were 230 and 438, respectively. The Vsearch plug-in was used to cluster the sequences into operational taxonomic units (OTUs) according to the 97% similarity principle. The Phylogeny plug-in was used to generate the phylogenetic trees. According to the fungus UNITE v.8.2 database (Nilsson et al., 2019), we used the Feature Classifier plug-in to annotate the species. Functional guilds and trophic modes were obtained from the FUNGuild database (Nguyen et al., 2016).

ITS is not a conserved region, making the constructed phylogenetic trees unstable. In this case, a hybrid tree based on the Silva 138 18S database and Unite 8.2 fungal dynamic database was constructed using q2-ghost-tree (Fouquier et al., 2016). The results of species composition and diversity analysis based on the OTU selected according to the ghost tree were basically consistent with those of the ITS tree, and the community assembly processes of the two phylogenetic trees were consistent. However, compared to phylogenetic trees based on ITS primers, the minimum frequency of samples decreased from 39,602 to 9,499, resulting in the loss of 76% of OTUs. Therefore, we decided to construct the ITS phylogenetic tree in this study.

The data processing part was mainly completed with R 4.1.2 (R Core Team, 2022) and SPSS 26.0 (SPSS Inc., Chicago, IL, USA). The vegan package provides tools for describing the community ecology. First, we used the rarefy command of the vegan package to standardize the OTU table according to the minimum sequence number of samples and constructed the rarefaction curve using the rare-curve command. Counts of different taxonomic categories were obtained using the Phyloseq package. Then we grouped taxa and functions which was not identified or whose average relative abundance was less than 1% into “others,” and defined dominant taxa as which a relative abundance more than 1%. The MicrobiotaProcess package was used to obtain alpha diversity indices. For alpha diversity analysis, the Shannon diversity index, Chao1 richness index, and Pielou evenness index were selected for one-way analysis of variance (ANOVA) and Tukey’s honest significant difference (HSD) test at P < 0.05 to compare species diversity. Beta diversity was tested by Principal coordinate analysis (PCoA) based on the OTU level, and the differences between fungal communities were tested by Permutational multivariate analysis of variance (PERMANOVA) of 999 permutations. Differential expression analysis was performed using the DESeq2 package based on the family level except for unassigned families. Comparisons were performed in four groups for each part to avoid the influence of other parts, and differentially expressed families (P < 0.05) were categorized as up regulated (log2FoldChange > 2) or down regulated (log2FoldChange < −2). A family that was up regulated or down regulated in two or more groups in each part was defined as differentially expressed in this part. We selected the null model described by Stegen et al. (2013) to classify community assembly processes. We calculated the β-nearest taxon index (βNTI) and modified Raup-Crick index (RC_Bray) using iCAMP package. βNTI values < −2 (homogeneous selection) or > +2 (heterogeneous selection) indicated deterministic processes; −2 < βNTI values <+2 indicated stochastic processes including homogenizing dispersal (RC_Bray < −0.95), dispersal limitation (RC_Bray > 0.95) and undominated (−0.95 < RC_Bray < 0.95). The Mantel test was performed on community structure and environmental variables. A heatmap was used to show Spearman’s correlation between the relative abundance of dominant phyla and environmental variables. The data were visualized using the ggplot2 package. All analyses of variance (ANOVA) were calculated using SPSS, with P < 0.05 considered statistically significant.

We only analyzed the soil properties in bulk soil, as we can only extract a small amount of rhizosphere soil for DNA extraction. There were significant differences in the physicochemical properties of the soil among the different variables ( Table 2 ). Soil total nitrogen (TN), SOC, pH, and CEC in Aershan were significantly higher than those in Genhe. Available phosphorus (AP) in Genhe was significantly higher than that in Aershan. Between the two seasons, there were no significant differences in TN, SOC, and CEC in Aershan, and no significant difference in pH in Genhe. In summary, compared with the season, the differences in soil physicochemical properties between the two sites were greater.

In total, 7,593,884 quality-filtered sequences were obtained from the 120 samples. After rarefying 39,602 sequences per sample ( Figure 2A ), 4,752,240 unique sequences were clustered into 7,788 operational taxonomic units (OTUs) with 97% identity. These OTUs belonged to 17 phyla, 60 classes, 138 orders, 362 families, and 864 genera.

We found that the differences in fungal composition were significant ( Figures 2B, C ), and there were also significant differences in functional guilds between the groups ( Figure 2D ). Basidiomycota, Ascomycota, Mortierellomycota, and Mucoromycota were the dominant phyla (average relative abundance > 1%). There were 13 dominant families, among which the families with a relative abundance more than 5% included Thelephoraceae, Mortierellaceae, Archaeorhizomycoaceae, and Inocybaceae. In terms of functional guilds, symbiotrophic fungi had an absolute advantage in the identified functions and the relative abundances of pathotrophic and saprotrophic fungi in Aershan were significantly higher than those in Genhe.

More families were differentially expressed between sites ( Figure 3C ), while there were fewer families with differential expression in different compartments and seasons ( Figures 3A, B ). Clavulinaceae, Hoehnelomycetaceae, Hydnangiaceae, Hygrophoraceae, Pyronemataceae, Rhizopogonaceae, Sporidiobolaceae, and Tritirachiaceae were significantly express families in bulk soil. Mytilinidiaceae, Russulaceae, Trichosporonaceae, and Tulasnellaceae were significantly expressed in rhizosphere soil. In July, only Archaeorhizomycetaceae and Trimorphomycetaceae were significantly expressed. As for October, Clavulinaceae, Cordycipitaceae, Leucosporidiaceae, Mortierellaceae, Mycosphaerellaceae, Pezizaceae, Schizoporaceae, and Vibrisseaceae were significantly expressed. A total of 69 families showed differential expression between these two sites, with 40 families dominating at Aershan and 29 families dominating at Genhe ( Table 3 ).

These results indicated that there were significant differences in community composition and functional guilds between different compartments, seasons, and sites.

The differences between the Shannon diversity, Chao1 richness, and Pielou evenness indices were significant in different groups, which consist of compartments, seasons, and sites. Shannon diversity and Pielou evenness indices of RGO were highest among these groups ( Figures 4A, C ). The group with the highest Chao1 richness index is BAJ ( Figures 4B ). Three-way ANOVA was used to analyze the effects of compartment, season, and site on alpha diversity. Among the single factors by Three-way ANOVA ( Table 4 ), the compartment had a significant effect on the Shannon (F = 10.179, P = 0.002) and Pielou (F = 12.58, P = 0.001) indices, whereas both season (F = 24.934, P < 0.0001) and site (F = 85.252, P < 0.0001) had a significant effect on the Chao1 index.

Principal coordinate analysis (PCoA) based on the Bray–Curtis distance ( Figure 4D ) indicated significant differences among multiple groups, especially between different sites (R² = 0.202, P = 0.001). The difference in fungal community structure between the rhizosphere and bulk soil in Genhe was relatively small, whereas the difference between the two seasons was significant. According to the results of the permutational multivariate analysis of variance (PERMANOVA), the grouping was credible (F = 9.088, P = 0.001).

According to the result of the β-nearest taxon index (βNTI, Figure 4E ), the values were dispersed around −2. Homogeneous selection was more important than other processes in Genhe fungal community assembly, whereas dispersal limitation was more important than other processes in Aershan ( Figure 4F ). Soil fungal communities in Aershan were dominated by stochastic processes with average relative importance of 70.48−84.76%, whereas those in Genhe were dominated by deterministic processes with average relative importance of 59.05−68.57%. The proportions of deterministic and stochastic processes in the fungal community assembly process significantly differed between the two sites.

According to the results, there were significant differences in alpha diversity and beta diversity of fungal communities in the soil of L. gmelinii. Differently dominant processes in the community assembly may cause the above results.

We performed mantle tests on the relationship between individual environmental factors and the fungal community structure of bulk soil and found that all soil properties, temperature, and precipitation significantly affected the fungal community structure ( Table 5 ). Among the environmental factors, TN (r = 0.734, P = 0.001) and pH (r = 0.174, P = 0.001) had the greatest effect on fungal community structure. Altogether, the correlation between the total environmental factors and community structure was not high (r = 0.174, P = 0.001).

In the correlation relationships between alpha diversity indices and environmental factors ( Figure 5A ), climatic factors had a significant effect on all three alpha diversity indices, and alpha diversity increased as temperature and precipitation increased. Additionally, the Chao1 evenness index was significantly influenced by the CEC, TN, pH, SOC, and AP. The effects of the CEC and pH on the Shannon diversity index were also significant.

To determine which taxa dominate the variation in fungal community structure in L. gmelinii forests under environmental changes, we evaluated the relationship between fungal families and the environment. In total, 16 of 17 dominant families were significantly correlated with environmental factors ( Figure 5B ), whereas the correlation between Hygrophoraceae and environmental factors was low and not significant. Russulaceae was significantly and negatively correlated with climatic factors. Inocybaceae were only significantly correlated with soil physicochemical properties. Mortierellaceae, Atheliaceae, Pyronemataceae, Pseudeurotiaceae, Hyaloscyphaceae, Tricholomataceae, and Thelephoraceae significantly correlated with AP, CEC, TN, pH, Tmax, and Tmin. A few of them were significantly correlated with Prec or SOC. Umbelopsidaceae was significantly and positively correlated with CEC, TN, and climatic factors. Archaeorhizomycetaceae and Serendipi,taceae were significantly correlated with environmental factors except Prec. Myxotrichaceae, Cortinariaceae, and Clavulinaceae were correlated with some soil physicochemical properties and Prec. The 17 dominant families were divided into two categories based on their correlations as a result of clustering.

Our results demonstrated that soil physicochemical properties and climatic factors played key roles in the diversity, composition and community structure of soil fungi.

Hosts produce exudates to regulate the rhizosphere microenvironment and provide carbon to microorganisms (Zhuang et al., 2020). Fungi combined with host roots to assist plants in absorbing carbon, nitrogen, phosphorus, and nutrients and play an important role in plant growth, development (Lu et al., 2018), and health (Berendsen et al., 2012). Many factors affect rhizosphere fungal communities, including vegetation type (Philippot et al., 2013), soil properties (Li et al., 2016), and biotic interactions (Paungfoo-Lonhienne et al., 2015). In our study, the fluctuations in Shannon and Chao1 indices were smaller in rhizosphere soil than in bulk soil. We supposed that the composition of fungal community in the rhizosphere was more stable. The alpha diversity of the fungal community in rhizosphere of L. gmelinii forest in Genhe was higher than that in bulk soil. However, the result was completely opposite in the soybean soil (Zhang et al., 2018). The selective effects of the root on specific taxa caused all the differences (Philippot et al., 2013). There were four families enriched in the rhizosphere soil of L. gmelinii. Russulaceae and Tulasnellaceae formed obligate mutualistic root symbionts with L. gmelinii, and they were also widely distributed in the mycorrhizae of orchid species (Li et al., 2021). Mytilinidiaceae and Trichosporonaceae were mainly involved in soil saprophytic processes and animal pathogenic transmission (Lourenço et al., 2020; Eagar et al., 2023). The families enriched in rhizosphere suggesting that those fungi were active in assisting nutrient uptake and transforming organic matter for L. gmelinii. Our results further confirmed that the host root system played a decisive role in soil fungi, and selected the specific fungal species. Exploring the status of rhizosphere fungal communities helps us understand the root health status of plants.

Seasonal changes in temperature and precipitation have important effects on aboveground and underground organisms. A previous study has shown that fungal communities have better adaptability and resistance to cold environments compared with bacterial communities (Buckeridge et al., 2013). Regardless of seasonal changes, Basidiomycetes and Ascomycetes are the main components of fungal communities (Xie and Yin, 2022). However, Zeng et al. (2020) reported that the relative abundance of Ascomycota decreases as the annual average precipitation increases. Previous studies have shown that saprophytic fungi dominate the soil during the cold season, whereas ectomycorrhizal fungi have a greater advantage during the growing season (Voříšková et al., 2014; Santalahti et al., 2016; Zeng et al., 2022). At our sampling sites, temperature and precipitation were concentrated in July; the temperature in October was below zero and the soil was relatively dry. High soil moisture promotes the decomposition of dead leaves and increases soil nutrients, thereby affecting the functional groups of the soil fungal community (Brockett et al., 2012). During the drier season, soil fungal abundance decreases compared to the wet season, but fungal diversity increases (Kivlin and Hawkes, 2016). In our previous study, variations in temperature and precipitation were important factors influencing seasonal changes in the fungal communities (Zhao et al., 2023). In fact, seasonal changes were masked by sites in this study.

Deterministic and stochastic processes were used to determine soil fungal community structure. Selection, drift, speciation, and dispersal are usually responsible for the deterministic fitness differences, species abundance, creates new species, and movement by the four processes, among which selection is a deterministic process and diffusion is a stochastic process (Vellend, 2010). The fungal community assembly in Aershan was dominated by stochastic processes, whereas it was dominated by deterministic processes in Genhe. The body size is crucial for dispersal (Luan et al., 2020), it is also an important stochastic factor that leads to the low similarity of the two communities in our study. The assembly process of different communities can cause differences in the community structure. Nevertheless, the role of stochastic processes in microbial community assembly remains unclear (Nemergut et al., 2013; Gao et al., 2020; Jiao et al., 2020). Combining niche breadth, deterministic processes have a strong impact on specialist species with narrow niche breadth, while the generalist species with wide niche breadth has stronger competitiveness (Graham et al., 2016). In general, the community dominated by stochastic process has more resources around it and less competitive pressure (Jiao et al., 2020). Based on this understanding of community assembly, we consider that the soil condition of Genhe was more limited, resulting in contemporary selection dominating the fungal communities. The soil conditions in Aershan were better, providing more opportunities for fungal communities to drift, diversify, and disperse.

Notably, contemporary environmental heterogeneity and historical events (geographic distance) were the two major factors affecting the biographical distribution of soil fungal communities (Ramette and Tiedje, 2007). The two sampling sites for L. gmelinii were 460 km apart, and significant differences were found in the soil physicochemical properties between the two sites, so we considered the effect of geographical distance. In previous studies, the geographical scale played an important role in fungal community structure (Tedersoo et al., 2014; Zhang et al., 2018; Adamo et al., 2021). He et al. (2017) even suggested that spatial heterogeneity has a more significant effect than temporal consistency.

In the biogeographical distribution of soil fungi, the contemporary environment played a more important role (Zeng et al., 2019). As important soil nutrient factors, C, N, and P not only have a profound effect on aboveground vegetation but also provide an environment for the living dynamics of soil fungi. Nutritional factors such as TN, AP, and SOC are important in the regulation of fungal communities (Lauber et al., 2008; He et al., 2016; Dang et al., 2017; Li et al., 2020; Luo et al., 2023). Long-term nitrogen enrichment decreases the relative abundance and diversity of ectomycorrhizal taxa; however, the saprotrophs respond positively to nitrogen enrichment. (Zhou et al., 2016; Morrison et al., 2016). Nitrogen deposition affects the soil fungal community in L. gmelinii forests annually and is a major cause of global environmental change (Galloway et al., 2008). SOC significantly affects the structure of soil fungal communities (Sul et al., 2013; Li et al., 2017). A few families in our study were significantly related to SOC; however, all related families belonged to saprophytic and symbiotic fungi, such as ericoid mycorrhizal fungi, orchid mycorrhizal fungi, and ectomycorrhizal fungi. They degrade cellulose (Li et al., 2017), and assist hosts in obtaining carbon (Suetsugu et al., 2021). Nevertheless, when the other four physicochemical properties were significantly correlated with alpha diversity indices and some dominant families, AP showed a significant opposite correlation ( Figure 5 ). Similar to our results for paddy fields in subtropical China, fungal richness and diversity indices declined significantly with increasing phosphorus rates (Liu et al., 2018). Bao et al. (2013) also found that the relationship between available phosphorus and fungal biodiversity showed a hump trend, slowly increasing before peaking and rapidly decreasing thereafter. As a non-nutritional factor, pH has a significant effect on soil fungal communities (Tedersoo et al., 2020). Compared with bacteria, fungi have a wider pH adaptation range (Rousk et al., 2010), but the dominant fungal taxa in different soil pH environments differ (Zhang et al., 2016). In our study, the correlation between the most dominant families and pH was high and extremely significant. CEC also affects the structure of fungal communities (Docherty et al., 2015). Some symbiotic taxa and the family Archaeorhiziomyceae were negatively correlated with pH, CEC, TN, and SOC. The relationship between these fungi and their hosts allows them to gain advantages in poor and acidic soils (Rosling et al., 2011; Vohník et al., 2016; Zhang and Sun, 2021). In addition to the factors studied, soil moisture content (Docherty et al., 2015) and salinity (Cao et al., 2021) also play important roles.

There were close relationships between TN, SOC, AP, pH, CEC and community structure, with TN and pH having the greatest effect ( Table 5 ). The soil total nitrogen of our samples ranged from 0.81 g/kg to 6.88 g/kg, and the pH values ranged from 3.71 to 6.23, which was consistent with the values in the study of Miyamoto et al. (2021). We did not conduct further experiments on its optimal gradient, but pH around 4.89 (4.89 ± 0.08) may be its superior pH by ANOVA. Combined with the results that soil total nitrogen was positively correlated with alpha diversity indices (including Shannon diversity index, Chao1 richness index and Pielou evenness index), we supposed that high soil total nitrogen is more advantageous for soil fungal communities. So appropriate nitrogen supplementation and soil pH regulation can maintain and regulate the health of L. gmelinii forest.