This study was performed on the Changbai Mountain (126^°55′-129^°00′E, 41^°23′-42^°36’N) in northeastern China. The local climate is a typical continental temperate monsoon climate with a daily average temperature of 5.9°C during the growing season (June to September). The average annual precipitation during the growing season can reach 958 mm. The mean annual precipitation and temperature is approximately 600 mm and 4°C, respectively.

To study the elevational pattern of soil fungal composition and diversity, we selected pure D. angustifolia population at 1,690, 1,810, 1,910, and to 2020 m a.s.l. along an elevational gradient on the north-facing slope of the Changbai Mountain. During October 1 to 7, 2018, three independent plots (10 m × 10 m) were set up in D. angustifolia population at each elevation. Ten to fifteen soil samples (0 ~ 20 cm organic layer) were sampled and pooled for each plot, using a sterile soil drill (5 cm in diameter, 20 cm deep). After removing the surface litter and humus layer, approximately 1 kg soils for each plot were collected. The soil samples were sieved (2 mm mesh) to remove stones, visible roots and residues and other debris, and divided into two sub-samples: one stored at −80°C for sequencing, and the other one stored at 4°C for soil physicochemical properties.

A soil-water (deionized water) (1:2.5 w/v) suspension was shaken for 30 min prior to measuring the pH with a pH meter (Thermo Scientific Orion 3-Star Benchtop, Cambridge, United Kingdom). Soil moisture content (SMC) was measured by comparing the fresh wet weight with the dry weight after drying at 120°C for 24 h. Soil organic carbon (SOC) and the total nitrogen (TN) content were measured using an elemental analyzer (Elementar, Langenselbold, Germany). Ammonium (NH₄⁺-N) and nitrate (NO₃⁻-N) nitrogen content were measured using a continuous flow analysis system (SKALAR SAN++, Breda, the Netherlands). The total phosphorus (TP) content was measured using a spectrophotometer, and available phosphorus (AP) content was measured using the colorimetric method upon extraction with 0.5 M NaHCO₃. The total potassium (TK) content was measured by digesting the soil with concentrated hydrofluoric acid, and available potassium (AK) content was extracted by acetic acid and ammonium leaching method. The extracted TK and AK content were determined using inductively coupled plasma atomic emission spectrometry (ICP-AES-7500, Shimadzu, Japan). Soil microbial biomass C (MBC) and biomass N (MBN) content were measured with a TOC analyzer (TOC-LCPH, Shimadzu, Japan). Soil mechanical compositions (Sand, Silt, Clay) were determined according to the method of Zhang et al. (2021). Three independent replicates per sample were performed for all the soil physicochemical properties.

Using the MOBIO Power Soil Extraction Kit (Mo Bio Laboratories, Carlsbad, CA, United States), soil total DNA was extract from 1 g of fresh soil according to the manufacturer’s instructions. The DNA was diluted in TE buffer (DNA Elution Solution-Ultra Pure Water). The DNA quantity and quality were detected using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, United States).

Fungal ITS rRNA region was amplified using primers ITS1 (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and ITS2 (5′- GCTGCGTTCATCGATGC -3′) (Fouts et al., 2012). A 6-bp barcode sequence unique to each sample was added to the primers for distinguishing multiple samples. The PCR reaction was performed in triplicate in a 25 μL mixture containing 2.5 μL of TransStart Buffer, 2 μL of dNTPs, 1 μL of each primer (10 ng/μL), and 30 ng of template DNA. The PCR conditions were as follows: pre-denaturation at 94°C for 5 min, 30 cycles of dunaturation at 94°C for 30 s, 55°C for 30 s, 72°C for 45 s, and a final extension at 72°C for 10 min. The PCR products were inspected by 2% agarose electrophoresis, and were purified using the AxyPrep DNA purification kit. Three independent PCR replicates per sample and then three PCR samples were pool at equal amount and PE300 paired-end sequenced on the Illumina Miseq v3 platfrom (2 × 300 bp). The raw sequences were uploaded to the Sequence Read Archive (SRA) database and accession number was SUB10527794.

Sequences were analyzed using QIIME (version 1.8 ¹) software on the Allwegene bioinformation cloud platform. ² The original PE reads were quality filtering following criteria: if the mean score < 20 or the length < 200 bp, and the ambiguities sequence were removed. The forward and reverse reads merged using PEAR software (version 0.9.8). The chimeras removed using Usearch (version 7.1). ³Operational taxonomic units (OTU) were generated at a similarity level of 97% using the UPARSE algorithm (Edgar, 2013). Taxonomic analysis was performed on the representative sequences of OTU, with a confidence threshold of 0.7, and the UNITE 8.2 fungi database was used for comparison. ⁴

Before further analysis, the sequences were normalized according to the lowest number of sequences for a single sample. The sequences occurring <1% were classified into “others.” Alpha diversity indexes (Richness, Chao1, Shannon index, PD=Phylogenetic diversity) were calculated in QIIME1 based on OTU table. Principal coordinate analysis (PCoA) was finished using R software (v.3.2.5, R Development Core Team, 2016) “Vegan” package based on Bray–Curtis dissimilarity at OTU level. Rarefaction curve was finished using R software (v.3.2.5, R Development Core Team, 2016) “microeco” package. Redundancy analysis (RDA) was performed using R software (v.3.2.5, R Development Core Team, 2016) “microeco” package based on OTU table and soil physicochemical parameters. Fungal functional guilds were assigned by using FUNGuild v1.0 and the differences of guilds among elevations were performed by one-way analysis of variance (ANOVA), Duncan test (p < 0.05). One-way analysis of variance (ANOVA) was used to detect the difference of soil physicochemical parameters among elevations using SPSS software (version 22.0). Pearson correlation analysis among fungal community composition (both phyla and genera) and soil physicochemical properties was also performed using SPSS software (version 22.0). Permutational multivariate analysis of variance (PERMANOVA) was used to test the differences in soil fungal composition among the four elevations, using Bray-curtis distance matrices (i.e., the adonis2 function of the vegan package) (Zhong and Fu, 2022).

All soil physicochemical parameters, except SOC and sand content, were significantly different (p < 0.05) among the four elevations (Table 1). Soil pH, and the content of SMC, nitrate nitrogen and available potassium declined with increasing elevation, while other soil physicochemical properties did not show a clear tendency (Table 1). Silt ranged from 0.5% (1,910 m a.s.l.) to 1.6% (2,020 m a.s.l.), and the clay ranged from 4.3% (1,690 m a.s.l.) to 5.1% (1,910 m a.s.l.).

The rarefaction curve (Figure 1) tended to flatten, indicating that the sequencing number was sufficient and reasonable to cover the fungal communities. The Chao1, richness, phylogenetic diversity (PD), and Shannon index of the soil fungi were significant different (one way ANOVA, p < 0.01) among elevations along the elevational gradient (Table 2). The alpha diversity indices including Chao1, richness, PD_whole_tree, and Shannon-Wiener decreased from 1,690 m to 1800 m a.s.l., then increased to 1910 m a.s.l., followed by decreasing tendency up to 2,020 m a.s.l.

Pearson correlation analysis showed that the soil fungal alpha diversity (richness and PD index) were significantly positively correlated with the content of soil NO₃⁻-N, but significantly negatively correlated with the content of soil NH₄⁺-N and TN (Table 3). PD index of soil fungi was significantly positively correlated with soil pH, while the Shannon-Wiener index of soil fungi was significantly negatively correlated with the content of soil TN (Table 3).

The ordination showed that the fungal communities were clearly separated by elevation (PCo1, 44.63% and PCo2, 25.36%) (Figure 2). The soil fungal beta diversity differed significantly among elevations (Figure 2, Anosim R = 0.97, p < 0.01). Moreover, there were significant differences in species composition between any two out of the four elevations (Supplementary Table S2, p < 0.05).

The soil fungal community composition (both at the phyla and genera level) was significantly different among the four elevations (Figure 3 and Supplementary Tables S1, S3). All the obtained sequences belonged to 8 phyla. The prevailing phyla was Basidiomycota (53% relative abundance), Ascomycota (26%), and Mortierellomycota (19%) across all the soil samples (Figure 3A). The relative abundance of these dominant fungal phyla changed remarkably with elevation (Supplementary Table S1), showing that the relative abundance of Basidiomycota had an increased tendency but Ascomycota showed a decreased tendency with increasing elevation (Figure 3A and Supplementary Table S1).

The dominant genera with a relative abundance of >1% were Mortierella (18.7% relative abundance), Russula (13.3%), Inocybe (7.7%), Archaeorhizomyces (4.7%), Laccaria (4.2%), Tricholoma (3.2%), Leotia (2.5%), Cortinarius (2.5%), Clavulina (2.4%), Entoloma (2.1%), Elaphomyces (2.0%), Piloderma (1.6%), Tomentella (1.3%), Solicoccozyma (1.1%), Gymnomyces (1.1%), Sebacina (1.0%) (Figure 3B). Russula was most abundant at the highest site (2,020 m a.s.l.), while Inocybe was most abundant at the lowest elevational site (Figure 3B). Similar to the phyla, the relative abundance of these most dominant fungal genera also changed remarkably with increasing elevation (Supplementary Table S3).

RDA revealed that soil properties (i.e., NH₄⁺-N, NO₃⁻-N, TN, SOC, MC, pH, AK, TK, and TP) were the key environmental factors that shaped the soil fungal community (Figure 4). The first two axes of the RDA accounted for 37.57% of the total variance. As shown in Figure 4, MC (p < 0.05), pH (p < 0.05), NH₄⁺-N (p < 0.05), NO₃⁻-N (p < 0.05), TN (p < 0.05) significantly influenced the fungal community. The soil fungal community structure at 2,020 m a.s.l. was significantly positively correlated with NH₄-N and TN, while that at 1,690 m a.s.l. and 1,800 m a.s.l. was significantly positively correlated with soil moisture content, pH, AK and NO₃^—N (Figure 4).

At the phyla level, the abundance of phyla was closely correlated with some certain, phyla-specific soil physicochemical factors (Table 4). For the two dominant phyla, for example, the abundance of Ascomycota was significantly negatively correlated with the content of soil TN, NH₄⁺-N, and AP (Table 4), while the abundance of Basidiomycota was significantly negatively correlated with soil TK (Table 4).

At the genus level, the abundance of genus was also closely correlated with some certain, genus-specific soil physicochemical factors (Table 4). For instance, the abundance of Ilyonectria, Inocybe and Archaeorhizomyces was significantly positively correlated with SMC, while the abundance of Russula was significantly negatively correlated with the content of SMC (Table 4).

The dominated functional groups of fungi were ectomycorrhizal, Endophyte, Undefined Saprotroph, Plant Pathogen, Animal Pathogen, Wood Saprotrophc, Lichenized, Ericoid Mycorrhizal and Arbuscular Mycorrhizal (Supplementary Table S4). Except Ericoid Mycorrhizal, Lichenized and Plant Pathogen, other soil fungal functional groups were significantly different among elevations along the elevational gradient (Supplementary Table S4). The absolute abundance of Arbuscular Mycorrhizal and Lichenized was highest at 1,690 m a.s.l., and the absolute abundance of Animal Pathogen, Ericoid Mycorrhizal and Wood Saprotroph was highest at 1,810 m a.s.l. while the absolute abundance of Ectomycorrhizal was highest at 2020 m a.s.l. (Supplementary Table S4).

In this study, we found significant differences in soil physicochemical properties among elevations along the elevational gradient from 1,690 m to 2020 m a.s.l. on the Changbai Mountain (Table 1), this is consistent with the results of Zong et al. (2014) who investigated the soil physicochemical properties in D. angustifolia population along elevational gradients on the same mountain. The climate changes dramatically with increasing elevation on mountains, which causes markedly changes in soil environment such as biogeochemical cycling and soil nutrients (Zong et al., 2014). In our study, for example, the soil water content and pH value were significantly higher at 1,690 m than at 2,020 m a.s.l., which indicates that soil moisture holding capacity decreases with the invasion of D. angustifolia. We found a significant lower value of soil organic carbon at higher elevation (2,020 m) and at lower elevation (1,800 m), which is consistent with results of previous studies (Meng et al., 2018; Luo et al., 2020).

Soil fungi play an important role in biogeochemical cycling and ecological process (Buee et al., 2009; Russo et al., 2012), but previous studies have focused on bacterial diversity and composition, with only a few focusing on fungi in mountain ecosystems. Therefore, we have only limited information about changes in fungi with elevation (Djukic et al., 2010; Adamczyk et al., 2019; Zhou et al., 2021). Our results showed that the alpha diversity (Shannon index, Chao1 index, Richness and PD index) of soil fungi changed significantly with elevation, showing that the alpha diversity decreased with increasing elevation (Table 3). However, this pattern does not seem to be widespread. The distribution pattern of soil fungal diversity in the literature showed declining, humped, U-shaped, or no change with increasing elevation. Yang et al. (2017) reported that soil fungal diversity decreased monotonically from 700 m to 2,600 m a.s.l. across various ecosystems on the Changbai Mountain. Ping et al. (2017) reported that soil fungal diversity in Pinus koraiensis forest showed a hollow curve’s pattern along an elevational gradient from 699 m to 1,044 m a.s.l. on the Changbai Mountain.

Our study indicated that the alpha diversity of soil fungi showed a decreasing trend with increasing elevation. Shen et al. (2014) showed that the Chao1 index of soil fungi on the Changbai Mountain was not correlated with elevation, but it had a strong correlation with soil pH. Ni et al. (2018) found that soil fungal Chao1 index increased with increasing elevation, and soil C/N was the most important environmental factor determining the soil fungal Chao1 index along an elevational gradient from 2,000 m to 2,500 m a.s.l. on the Changbai Mountain. However, Zhang et al. (2022) found that both the Chao1 and Shannon index of soil fungi decreased significantly from 2,785 m a.s.l. to 4,578 m, and the climate and soil properties had opposite effects on them in Tibetan Plateau.

In our present study, soil TN, NH₄⁺-N, and NO₃⁻-N were significantly correlated with fungal richness (r = −0.53, p < 0.05; r = 0.57, p < 0.05; r = −0.59, p < 0.01) along the elevational gradient. As the soil N content decreased with increasing elevation (Zong et al., 2014), resulting thus in a reduction in soil fungal diversity found in the present study. Similar studies have also demonstrated that soil nutrients are an important environmental factor affecting the distribution pattern of soil fungal diversity (Newsham et al., 2016; Yang et al., 2017).

In our present study, the dominant phyla are Basidiomycota and Ascomycota, which is consistent with findings of previous studies on fungal composition along elevational gradients on Changbai Mountain. However, the elevational patterns of Basidiomycota and Ascomycota differed among studies. Our study found that Basidiomycota showed an increasing trend, while Ascomycota showed a decreasing trend with elevation from 1,690 m to 2,020 m a.s.l. However, Shen et al. (2014) and Ni et al. (2018) found that Basidiomycota showed an overall increasing trend, while Ascomycota’s did not change or tended to decrease with elevation on Changbai Montain. Ping et al. (2017) found that the abundance of Basidiomycota first decreased from 699 m to 937 m a.s.l., then increased from 937 m to 1,177 m a.s.l., while the abundance of Ascomycota increased from 699 m to 937 m a.s.l., and then decreased from 937 m to 1,177 m a.s.l. Yu et al. (2019) and Zhang et al. (2021) found that Ascomycota and Mortierellomycota are the predominant fungi in Tibetan grassland communities. This difference in dominant phyla may be closely related to the fact that the vegetation (or plant ecosystem) and soil physicochemical properties differ significantly between the two regions (Jin et al., 2018; Zong and Fu, 2021).

Soil fungal composition is affected by multiple factors such as vegetation composition, soil physicochemical properties, and microclimate along elevational gradients (Ni et al., 2018). A previous study carried ou on the Changbai Mountain found that the aboveground vegetation composition was closely related to the fungal composition, and dominant plant species significantly affected the fungal composition (Ni et al., 2018). In our present study, the dominant Basidiomycota and Ascomycota were mainly affected by soil physicochemical properties because the vegetation population (i.e., D. angustifolia population) did not change along with the elevational gradient. The fungal Basidiomycota and Ascomycota are involved in the soil organic metabolism (Luo et al., 2021) and thus their abundances were significantly determined by soil organic matter content as a result of decomposition of plant residues (Li et al., 2019).

Our results also showed that the functional guilds of soil fungi at high elevation were mainly ectomycorrhizal fungi and plant pathogens. Ni et al. (2018) also reported that the dominant functional fungi were ectomycorrhizal fungi in the alpine tundra on the Changbai Mountain. Similarly, Timling et al. (2014) reported that the dominant functional guilds in the arctic were the ectomycorrhizal fungi. However, in this study, the alpine tundra population of D. angustifolia is an invasive plant, and the effect of the original tundra species on the soil has not completely replaced by the invasive plant, so that the soil fungal functions were still ectomycorrhizal fungi and plant pathogens.