Our study area extends for 988 km² in the Wudalianchi Volcano (125°45′-126°30′ E; 48°30′-48°50′ N), which is a well-preserved volcano in Heilongjiang Province, China (Figure 1). The climate is characterized by a continental monsoon pattern, with a mean annual temperature of -5 °C and annual precipitation averaging 476.33 mm.

Soil samples were collected in the selected area in August 2023. To ensure sampling uniformity, samples were collected from locations with similar elevation and soil types: moss, herb, shrub, broadleaf forest, and mixed coniferous-broadleaf forest. The basic information of sample plots across different vegetation types is presented in Table 1. Within the plots of each vegetation type, three 10 m × 10 m quadrats were established, with approximately 20 m spacing between adjacent quadrats. When sampling each plot, tree trunks and plot edges were avoided as much as possible. We removed surface litter and discarded stones and roots from the soil sample. Soil samples were collected from the 0–5 cm layer using the “five-point sampling method”. Soils from the same quadrat were thoroughly mixed to form a composite sample weighing approximately 1000 g. All final composite samples were transported to the laboratory for analysis. The collected soil samples were sieved through a 2 mm mesh to remove roots and leaves. One portion was stored at -80°C for DNA extraction, and the other portion was analyzed for soil physicochemical properties.

The soil water content (WC) was determined using the drying method, soil pH (pH) was measured using a soil-water ratio of 1:5, and total nitrogen (TN) and total carbon (TC) were measured by the fully automated carbon and nitrogen analyzer Elementarvario ELIII (Elementar Analysensysteme GmbH, Langenselbold, Germany). Total phosphorus (TP) was determined using melted molybdenum, antimony, and scandium colorimetry.

Total microbial genomic DNA was extracted from soil samples using the E.Z.N.A.^® soil DNA kit (Omega Bio-tek, Norcross, GA, U.S.) according to manufacturer’s instructions. The quality and concentration of DNA were determined by 1.0% agarose gel electrophoresis and a NanoDrop 2000 spectrophotometer (Thermo Scientific, United States). The ITS region was amplified using primers ITS1F (5’-CTTGGTCATTTAGAGGAAGTAA-3’) and ITS2 (5’-GCTGCGTTCTTCATCGATGC-3’) by T100 Thermal Cycler PCR thermocycler (BIO-RAD, USA). The PCR reaction mixture including 4 μL 5 × Fast Pfu buffer, 2 μL 2.5 mM dNTPs, 0.8 μL each primer (5 μM), 0.4 μL Fast Pfu polymerase, 10 ng of template DNA, and ddH₂O to a final volume of 20 µL. PCR amplification cycling conditions were as follows: initial denaturation at 95°C for 3 min, followed by 27 cycles of denaturing at 95°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 45 s, and single extension at 72°C for 10 min, and end at 4°C. The PCR product was extracted from 2% agarose gel and purified using the PCR Clean-Up Kit (YuHua, Shanghai, China) according to manufacturer’s instructions and quantified using Qubit 4.0 (Thermo Fisher Scientific, USA). Purified amplicons were pooled in equimolar amounts and paired-end sequenced on an Illumina Nextseq2000 platform (Illumina, San Diego,USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw sequencing reads were deposited into the NCBI Sequence Read Archive (SRA) database (Accession Number: PRJNA1402311).

Raw FASTQ files were demultiplexed using an in-house Perl script, and then quality-filtered by fastp (https://github.com/OpenGene/fastp, version 0.19.6) and merged by FLASH (http://www.cbcb.umd.edu/software/flash, version 1.2.11) with the following criteria: (i) The reads were truncated at any site receiving an average quality score of <20 over a 50 bp sliding window, and the truncated reads shorter than 50 bp were discarded; reads containing ambiguous characters were also discarded; (ii) only overlapping sequences longer than 10 bp were assembled according to their overlapped sequence. The maximum mismatch ratio of overlap region is 0.2. Reads that could not be assembled were discarded; (iii) samples were distinguished according to the barcode and primers, and the sequence direction was adjusted, exact barcode matching, 2 nucleotide mismatch in primer matching.

The resulting high-quality sequences were then clustered into operational taxonomic units (OTUs) using UPARSE 7.1 with a 97% sequence similarity level. The most abundant sequence for each OTU was selected as a representative sequence. The OTUs assigned to spike-in sequences were filtered out and reads were counted. A Standard curves (based on read counts versus spike-in DNA copy number) for each sample were generated, and the quantitative abundance of each OTU in a sample was determined. The taxonomy of each OTU representative sequence was analyzed by RDP Classifier version 2.2 against the fungal ITS reference database (UNITE version 9.0) using a confidence threshold of 0.7. To account for differences in sequencing depth across samples, the OTU abundance table was normalized using a rarefaction approach prior to downstream statistical analyses.

Sequencing reads were normalized by rarefaction to an equal depth across all samples prior to downstream analysis to account for uneven sequencing effort and to ensure unbiased community comparisons. Following sequencing, Origin 2021, Excel 2021, and SPSS 20.0 were used for data organization and statistical analysis. Soil physicochemical characteristics, fungal community diversity indices, and community composition were analyzed using one-way analysis of variance (ANOVA), and group differences were compared using Tukey’s test. The alpha diversity (Chao index, Shannon index, and Simpson index) of the fungal community was calculated using the aforementioned OTU table were characterized. Beta diversity, calculated based on the Bray-Curtis distance to assess differences in microbial fungal structure among samples, was statistically tested for variability within treatments using PERMANOVA with 999 permutations (Adonis function) and visualized via Principal Coordinates Analysis (PCoA) plots. Redundancy analysis (RDA) was employed to reveal the relationships between the fungal community and key environmental factors. Statistical significance of the overall model and individual environmental factors was validated via 999 permutations (conducted in R software using the vegan package). Using FUNGuild v1.0 software, the potential ecological roles of soil fungi were examined. Spearman correlation analysis was used to examine the relationship between fungal functions and soil physicochemical properties.

Co-occurrence network analysis was constructed by Spearman correlations using the “corr.test” function in the “psych” R package. To enhance the ecological relevance and stability of the network, only OTUs with relative abundances > 0.01% of the total fungal OTUs were retained for construction, thereby excluding rare taxa that may introduce stochastic variations. The correlation matrix was filtered using a threshold of |r| > 0.6, which was determined based on random−matrix−theory methods. This threshold not only retains statistically significant correlations but also focuses on strong ecological associations, reducing the inclusion of weak, potentially spurious links. In the resulting network, each node in the network represents a fungal taxon, and the edges signify correlations among different nodes. Using the igraph package in R-4.3.1, a set of topological coefficients of each network, including network diameter, modularity, graph density, average degree, and average path length, were calculated to assess the complexity of the network. The network was visualized using Gephi 0.9.2.

During the transition from the M community to the C community, soil TC, total nitrogen TN, pH, N/P ratio, C/P ratio, and moisture content all exhibited a trend of initially increasing followed by a decrease. Moreover, soil TC (94.75 g/kg), TN (5.68 g/kg), and WC (43.59%) in the S community were significantly higher than in other vegetation types (p < 0.05, Figure 2). The C/N ratio showed an opposite trend with vegetation type, with the H community exhibiting a significantly lower C/N ratio (11.98) than other vegetation types (p < 0.05). While TP exhibited a trend of decreasing, it was highest in the M community (3.30 g/kg) compared to other vegetation types.

Across all vegetation types, the Chao and Shannon indices of soil fungal alpha diversity showed an N-shaped trend, peaking in the S community and falling in the B community (p < 0.05). The Simpson index, on the other hand, exhibited the reverse pattern, peaking in the B community and falling in the S community (Figure 3). The results indicate that there were significant differences in the diversity and richness of soil fungal communities in different vegetation types (p < 0.05).

PCoA based on Bray-Curtis dissimilarities revealed a significant separation of fungal community structures among the different vegetation types (PERMANOVA test with 999 permutations: R² = 0.92, p = 0.001; Figure 4). The first two principal coordinates (PC1 and PC2) explained 29.42% and 26.45% of the total variance, respectively, cumulatively accounting for 55.87%. In the ordination space, samples exhibited clear separation by vegetation type. The M community formed a distinct cluster in the first quadrant, whereas the C community clustered in the second quadrant; both groups were located on the positive side of PC2. Conversely, the B community was situated in the third quadrant, and the S and H communities co-clustered in the fourth quadrant, all on the negative side of PC2. This spatial separation indicates distinct soil fungal community structures across vegetation types. Specifically, the proximity of the S and H communities suggests a higher degree of compositional similarity between them compared to other vegetation types.

High-throughput sequencing yielded 1413798 valid fungal sequences, which were aggregated into 2821 operational taxonomic units (OTUs).

At the phylum level (Figure 5A), the dominant soil fungal communities were Ascomycota (43.39%-71.54%), Basidiomycota (5.36%-53.21%), and Mortierellomycota (2.38%-11.95%). Ascomycota predominated in both the H (69.25%) and C (71.54%) communities, with the highest abundance observed in C, which was significantly greater than in the S (58.55%), B (43.48%), and M (43.39%). Basidiomycota was the most dominant in the B (53.21%), with its abundance following the order: B > M > S > C > H. Mortierellomycota was the most dominant in the H community (11.95%), with a relative abundance significantly higher than in other types (p < 0.05).

At the genus level (Figure 5B), considerable compositional differences were observed across vegetation types. The dominant genera of fungal communities were unclassified_p_Ascomycota, Cortinarius, Mortierella, and Scleroderma. Unclassified_p_Ascomycota was the most dominant in the C community (19.49%), with its abundance following the order: C > S > B > H > M. Cortinarius showed the highest relative abundance in the B community (44.29%), but in M, H, and C, its relative abundance was less than 1%. Mortierella exhibited its highest relative abundance in the H community (11.57%), which was significantly greater than that in the other types (p < 0.05). In contrast, Scleroderma reached its maximum relative abundance in the M community.

Redundancy analysis (RDA) was constructed with log-transformed fungal phylum relative abundance data as the response variable and selected soil physicochemical properties (TP, pH, C/N ratio, N/P ratio, etc.) as explanatory variables. At the phylum level (Figure 6A), the first two axes collectively explained 81.99% of the variance in fungal community composition (RDA1: 62.15%; RDA2: 19.84%). TP (R² = 0.749, p = 0.001), pH (R² = 0.750, p = 0.001), C/N ratio (R² = 0.530, p = 0.001), and N/P ratio (R² = 0.415, p = 0.043) were identified as significant drivers of community structure (p < 0.05). Notably, Ascomycota relative abundance correlated positively with N/P ratio and pH, but negatively with TP. Conversely, Basidiomycota was positively associated with TP and C/N ratio, and it was significantly enriched in the M community.

At the genus level (Figure 6B), the first two RDA axes explained 43.11% and 32.25% of the variance, respectively, yielding a cumulative explanation of 75.36%. The C/P ratio (R² = 0.501, p = 0.014) was identified as the most influential driver of fungal community structure among the measured environmental factors. The relative abundances of Cladophialophora and Penicillium were significantly positively correlated with the C/P ratio, N/P ratio, and pH, and both genera were predominantly enriched in the S and C communities. In contrast, Mortierella and Talaromyces were primarily associated with the H community. Scleroderma and Laccaria were mainly found in the M community and showed a significant positive correlation with TP. Finally, Archaeorhizomyces and Cortinarius were predominantly restricted to the B community.

As shown in Figure 7, co-occurrence network analysis at the OTU level revealed distinct topological structures across vegetation types and identified key fungal taxa. Ascomycota and Basidiomycota were identified as the main keystone taxa, despite dynamic fluctuations in their relative abundances across vegetation types. In particular, the C community had the lowest relative abundance of Ascomycota (57.69%), while the M community had the greatest (64.07%). Basidiomycota, on the other hand, was comparatively lower in the H community (13.03%), but it was greater in the M (20.35%) and C (20.69%) communities.

Network topology parameters further uncovered pronounced structural disparities in fungal communities across vegetation types (Table 2). Vegetation types S, H, and C formed highly complex, tightly interconnected networks, each characterized by an average path length and network diameter of 1. This topological configuration implies that all nodes were either directly connected or linked via at most one intermediary, suggesting a high potential for efficient propagation of ecological information or resource exchange. Among these, network S exhibited the largest scale (487 nodes, 30,790 edges) and strongest connectivity (average degree: 126.45). Network C displayed the highest density (0.293), and its greater proportion of positive correlations (53.55%) indicated a community structured primarily by synergistic biotic interactions. In contrast, network M featured a distinct “high-modularity, low-connectivity” architecture. It yielded the highest modularity index (0.74) across all vegetation types, accompanied by a notably high proportion of positive correlations (60.63%). This pattern points to a community organized into multiple ecologically discrete modules, predominantly governed by cooperative interactions. Network B, by comparison, exhibited a sparsely connected and structurally simplistic topology. It contained the fewest nodes and edges, alongside the longest average path length (1.72) and largest network diameter (2) across all groups. These metrics align with weak direct interactions among fungal taxa and reflect the lowest overall network complexity relative to other vegetation-associated networks.

The ecological functions of fungal communities were predicted using the FUNGuild database. After filtering unassignable fungi, nine trophic modes were identified based on the top 10 relative abundances (Table 3). Symbiotroph was the most prevalent trophic mode in the B community (49.46%), showing the highest relative abundance among all vegetation types. Saprotroph was the predominant trophic mode in the H (32.97%), S (29.49%), and C (32.06%), although its relative abundance varied across these vegetation types. Furthermore, the relative abundance of Saprotroph-Symbiotroph was significantly higher in the H (11.72%) community than in other vegetation types, while Pathotroph was significantly less abundant in the B (0.37%) and C (0.99%) communities compared to the rest. In contrast, Pathotroph-Symbiotroph (8.59%) reached its highest relative abundance in the C community.

Correlation analysis (Figure 8) between fungal trophic types and soil physicochemical properties revealed that the relative abundance of Symbiotroph was significantly positively correlated with TP (p < 0.05) and significantly negatively correlated with pH (p < 0.01). The soil C/N ratio showed a strongly significant negative correlation with the Saprotroph-Symbiotroph (p < 0.001) and a significant positive correlation with the Pathotroph-Symbiotroph (p < 0.01). In contrast, Pathotroph exhibited a significant positive correlation with pH (p < 0.05) and a significant negative correlation with the C/N ratio (p < 0.05).

Based on the data obtained in this study, the results suggest that soil fungal alpha diversity (as indicated by the Chao and Shannon indices) is significantly higher in shrub communities (S) than in moss (M) or broadleaf forest (B) communities within this volcanic habitat. These diversity patterns appear to be closely linked to variations in key soil properties, including soil total carbon, and total nitrogen content. In accordance with the general understanding that environmental conditions shape fungal community composition (Shen et al., 2015) and that alpha diversity serves as an indicator of species abundance and richness (Li et al., 2022), our observations in this volcanic ecosystem align with broader ecological principles. Similar patterns have been reported by Miao et al. (2025), who observed that shrub communities promote litter input and root exudation, which in turn are associated with enhanced soil structure and nutrient retention. These improved soil conditions are consistent with a favorable environment for microbial growth and reproduction (Dai et al., 2022) and may contribute to creating diverse niches for fungi, which may reduce interspecific competition and consequently support higher alpha diversity of the fungal community (Yang et al., 2022). Furthermore, such conditions can sustain a greater variety of functional fungi (e.g., decomposers and symbionts), thereby being associated with enhanced overall community complexity (Xu et al., 2021). This is analogous to findings in other volcanic contexts, where ash deposition has been noted to influence fungal diversity through modifications of soil physicochemical properties (Xie et al., 2024). These findings reinforce the ecological importance of soil properties in shaping fungal diversity, and suggest that comparable vegetation–soil–microbe dynamics may also occur in other volcanic systems. Thus, the elevated alpha diversity in shrublands can be interpreted as a direct biological response to the improved soil habitat they engineer.

Vegetation type is a key determinant of soil fungal community characteristics (Cheng et al., 2024). In this study, Ascomycota (43.39%-71.54%) and Basidiomycota (5.36%-53.21%) are the dominant fungi communities of soil fungi in Wudalianchi. The volcanic lava platform is characterized by less water and good permeability, which are consistent with the characteristics of Basidiomycetes and Ascomycetes that prefer soil permeability (Ni et al., 2021). Ascomycota, in particular, are known to exhibit a competitive advantage in nutrient-limited environments (Zou et al., 2024), a pattern further supported by prior research (Siles and Margesin, 2017; Adamczyk et al., 2019). Cortinarius, Mortierella, and Scleroderma emerged as dominant fungal genera, with broad environmental adaptability and the ability to form mutualistic symbioses with plants. These traits are crucial for sustaining soil nutrient cycling and ecosystem stability (Truong et al., 2024). The genus Mortierella not only thrives in oligotrophic environments, a trait that facilitates its survival and proliferation (Hansen et al., 2019), but also exhibits notable stress tolerance (Ochoa-Henriquez et al., 2024). In broadleaf forests, the relative abundance of Cortinarius was significantly higher than in other vegetation types, a result consistent with the findings reported by Tuo et al. (2022). This distribution aligns with the specific habitat requirements of ectomycorrhizal fungi, as broadleaf forests generally feature deeper root systems, higher soil permeability, and more stable soil structure than moss, herb, or shrub communities, and typically experience less soil disturbance (Lann et al., 2024). These conditions provide suitable habitats and abundant organic substrates for ectomycorrhizal fungi such as Cortinarius (Lindahl et al., 2021). The highest relative abundance of Scleroderma was observed in the moss (M) community, a pattern consistent with the findings of Chen et al. (2023). This is primarily attributed to its role in enhancing phosphorus availability through the secretion of enzymes and organic acids, which improves the soil microenvironment of the moss community, promotes the accumulation of organic matter and phosphorus, and thereby sustains nutrient cycling within this community (Zhang et al., 2024).

Microbial network analysis serves as an effective tool for investigating interactions within microbial communities under changing environmental conditions (Tao et al., 2018). The complexity and topological characteristics of soil fungal co-occurrence networks reflect the spatial distribution patterns of microorganisms across different vegetation types while elucidating species coexistence relationships (Freitas et al., 2024). Li et al. (2024) discovered that Ascomycota and Basidiomycota often serve as keystone species in these networks, playing a central role in maintaining network stability and functionality—a finding that aligns with the results of the present study. Their remarkable capacity for organic matter decomposition, broad ecological adaptability, and capacity for symbiosis significantly underpin community stability and ecosystem functioning (Qu et al., 2024).

Abrego et al. (2020) found that in forest ecosystems, microbial communities exhibit an interaction pattern dominated by symbiotic associations, complemented by competitive interactions, which is consistent with the results of this study. Symbiosis reduces direct resource competition and mitigates interspecific exclusion by promoting resource complementarity, functional division of labor, and positive feedback mechanisms (Tedersoo et al., 2020), thereby contributing to ecosystem stability. Furthermore, fungal communities in these systems form complex network architectures (Tian et al., 2025a). Within nutrient-limited volcanic soils, the complexity of fungal co-occurrence networks varies significantly across vegetation types (Yang H.L et al., 2024). Fungal community diversity drives the enhancement of topological structure and the stabilization of key community functions (Tian et al., 2025b), ultimately shaping soil microbial communities with increased resilience to disturbance. This underscores that soil fertility and heterogeneity are primary drivers of microbial interaction complexity. This study found that the fungal co-occurrence network in shrub communities exhibits more nodes, more connections, and a higher average degree compared to other communities, indicating tighter connectivity and more robust interaction patterns. These structural characteristics of the network are closely related to their unique habitat conditions. Shrub communities predominantly occur in complex microtopographic areas such as lava fissures and folds (Liu et al., 2025), where microtopography promotes the accumulation of soil, moisture, and litter—a process that not only results in complex and uncertain distribution patterns of soil properties but is also critical during initial pedogenesis (Bonnet et al., 2023). This mechanism may contribute to the more robust nodes and connectivity pathways observed in fungal−bacterial interaction networks within shrub communities (Mao et al., 2014).

Soil fungal trophic types are closely associated with vegetation types (Hou et al., 2024). While Whitaker et al. (2015) reported that nutrient enrichment generally promotes saprotrophs while inhibiting symbiotrophs, our study revealed a divergent pattern under high TP conditions. Elevated soil P content increases the availability of inorganic nutrients, which reduces plants’ reliance on nutrients derived from organic matter decomposition. Consequently, this diminishes the resource base and competitive advantage of saprotrophs, compressing their ecological niche and leading to reduced diversity–a result consistent with our findings (Yang B. et al., 2024). In the M community, higher TP content slows the decomposition of soil organic matter (Clausing et al., 2021), thereby reducing the substrate availability for saprotrophs. Concurrently, the microbial community structure shifts toward symbiotroph dominance. These combined factors further reduce the competitiveness of saprotrophs. Consequently, the relative abundance of saprotrophs in the M community is significantly lower than in other vegetation types (Ma et al., 2021). This illustrates how specific soil chemical factors (e.g., P availability) can override general patterns to shape functional guild distributions.