Ectomycorrhizal fungal community varies across broadleaf species and developmental stages
Ectomycorrhizal fungi (EMF) form mutualistic relationships with plant roots, where their extramatrical mycelia extend many centimeters into soil, facilitating root nutrient (e.g., nitrogen) uptake and, in turn, receiving photosynthesized carbon from host plant , . This symbiosis is crucial for the survival, health, and growth of temperate forests , . For example, EMF hyphae enable plants to access a larger soil bank , secrete hydrolytic exoenzymes to mobilize nutrients from soil mineral–organic matter , and trigger plant defense mechanisms against pathogens and abiotic stressors . Broadleaf trees in temperate forests, particularly those of the Betulaceae, Fagaceae and Sapindaceae families, rely heavily on EMF symbiosis for establishment and growth . Investigating the assembly of the root EMF communities during host development is essential for understanding the survival strategies of these obligately ectomycorrhizal trees in temperate ecosystems . Previous studies document that EMF diversity increases as hosts mature. For instance, EMF richness of Quercus liaotungensis notably increases across sapling (1–3 years), juvenile (20–30 years) and mature stage (50–70 years) . Similarly, the Shannon index of Fagus sylvatica EMF rises as hosts mature from early (~ 20 years) to later stages (~ 100 years) . These shifts are attributed to the expansion of well-established root systems in adult trees, which allows for greater soil exploration and facilitates the recruitment of new fungal species via ectomycorrhizal networks (EMNs) shared by neighboring trees . Albeit typically increasing with host maturity, EMF richness and evenness may decrease at very late host stages due to the competitive pressures of conspecific neighbors, which can elevate mortality rates via species-specific pathogens – . In contrast, a study on Betula papyrifera suggests minimal changes in EMF diversity between the host stages of 5 and 65 years, probably because paper birch regenerates through stump sprouts and thereby maintains the host-specific EMF . In addition, Corylus avellana exhibits little variation in EMF community composition between the ages of 50–100 years and 100–200 years, with younger (< 50 years) stages differing more significantly . Furthermore, the EMF community compositions of young (20–30 years) and mature (50–70 years) Quercus liaotungensis are more similar to each other than to saplings (1–3 years) . These suggest that adult host trees serve as major key repositories of forest biomass and pivotal hubs in EMNs, facilitating local coexistence of EMF , . Despite extensive research on the EMF of broadleaf trees, few studies have examined the effects of host development across multiple coexisting tree species – . The zonal top broadleaf Korean pine mixed forests in Changbai Mountain are renowned for their high biodiversity and complex stand structure . For instance, the second poplar–birch forests, which developed following disturbances in the original broadleaf Korean pine forests, mainly comprise of Betula platyphylla and Populus davidiana . These forests serve as vital repositories of various EMF symbionts and are relatively short-lived, typically declining after 50–60 years . This study aims to investigate the succession patterns of root EMF communities across the sapling, juvenile, and adult stages of three typical broadleaf trees— Acer mono (Sapindaceae), Betula platyphylla (Betulaceae) and Quercus mongolica (Fagaceae). By utilizing Illumina MiSeq sequencing, we examine root EMF communities and identify specific EMF taxa associated with each host species and developmental stage. The correlations among the root EMF ASVs across different developmental stages were statistically analyzed. 

Database and major EMF taxa After quality control and denoising, we obtained 1,153,982 high-quality concatinated nonchimeric reads from 3,875,730 raw paired-end reads across 77 root samples. These sequences were mapped to 5,726 amplicon sequence variants (ASVs), of which 1,108 ASVs (579,294 sequences) were identified as EMF species (Supplementary Table ). These EMF ASVs spanned major phyla of Basidiomycota and Ascomycota, families such as Thelephoraceae, Russulaceae, Sebacinaceae, and Inocybaceae (Supplementary Fig. ). The dominant EMF genera (> 1%) were Tomentella , Russula , Inocybe , Sebacina , Tuber , Cortinarius , Amphinema , Humaria , Cenococcum , Piloderma , Trichophaea , Lactarius , Leccinum , Phaeohelotium , and Thelephora (Supplementary Fig. ). Alpha diversity of the EMF We compared the alpha diversity (such as Observed number of ASVs, Shannon, Simpson, Chao1, ACE, and Pielou indices) differences of root EMF by Kruskal‒Wallis test to discern the patterns of fungal richness as influenced by the maturation of hosts (Supplementary Table S3, Fig. ). Furthermore, the EMF alpha diversity indices (Observed number of ASVs, Chao1, ACE) were significantly influenced by the tree species at particular developmental stages ( P < 0.001, Supplementary Table S3). For example, from sapling to juvenile stage, EMF diversity remained relatively stable in Acer mono ( P > 0.05), marginally decreased in Betula platyphylla (0.05 < P < 0.1), whereas marginally increased in Quercus mongolica (0.05 < P < 0.1). Moreover, at the sapling stage, EMF diversity was significantly higher in Betula platyphylla than in Acer mono and Quercus mongolica ( P < 0.05). At the juvenile stage, EMF diversity ranked as Acer mono < Betula platyphylla < Quercus mongolica ( P < 0.05). Unlike those at the early stages, the EMF diversity in the adult hosts did not significantly vary across tree species ( P > 0.05, Supplementary Fig. S4). The above results implicate that as hosts mature, the EMF diversity increases and transitions from early-stage heterogeneity to late-stage convergence among the different host species. Beta diversity of the EMF Nonmetric multidimensional scaling (NMDS) and Permutational multivariate analysis of variance (PERMANOVA) demonstrated significant compositional shifts across stages ( P < 0.001, Fig. a), with adult communities showing lower dissimilarity among tree individuals ( P < 0.05, Fig. b; Supplementary Fig. S5). Particularly, unlike Acer mono and Quercus mongolica , the EMF communities in Betula platyphylla were not significantly different between the sapling and adult stages ( P > 0.05, Fig. ). Our findings indicate that the variance in the root EMF communities among the different hosts decrease markedly with host maturity, revealing a trend toward age-related convergence in the root EMF community composition. Cooccurrence network analysis of the EMF To explore the interaction patterns within EMF communities, we constructed a meta-community cooccurrence network based on Spearman’s correlation ( P < 0.05) and examined how these networks varied with host maturity and tree species. We identified stage-specific EMF taxa and categorized them into four clusters: sapling-, juvenile-, and adult-enriched ASVs, along with a set of common ASVs (Fig. ). Cooccurrence network demonstrated high degrees of positive correlations among EMF taxa (positive correlations exceeding 99%, Supplementary Table S4) in all three tree species: 6,290 correlations (edges) among 501 ASVs (nodes) for Acer mono , 7,471 edges and 631 nodes for Betula platyphylla , and 11,983 edges and 540 nodes for Quercus mongolica . The modularity values were greater than 0.4 (Supplementary Table S4), indicating that the networks exhibited distinct modular structures. The node-level topological features (degree, closeness centrality, and eigenvector centrality) were significantly greater at the adult stage compared to earlier stages ( P < 0.05, Supplementary Fig. S6). Furthermore, the adult-stage subnetworks encompassed more nodes and edges, higher average degree and graph density, but lower average path length and diameter, relative to its earlier stages (Supplementary Table S4). These findings demonstrate that the EMF communities are more aggregated and cohesive with host maturity, occupying a central position in cooccurrence networks. Stage-associated indicators of specific host species To identify key ASVs driving EMF community assembly, we conducted linear discriminant analysis (LDA) effect size (LEfSe) to detect EMF indicators of specific host species and stage. As for Acer mono , the adult stage was predominated by 20 EMF species primarily belonging to Atheliales and Pilodermataceae; the earlier stages were enriched in Pyronemataceae, Boletales and Russula vinosobrunneola . As for Quercus mongolica , the adult stage were dominated by 19 EMF species predominantly belonging to Atheliales, Pilodermataceae, Thelephora sp., and Cenococcum sp.; the sapling and juvenile stages were closely associated with Thelophora sp. and Elaphomyces papillatus . Betula platyphylla exhibited 11 EMF species enriched at the sapling stage and 7 EMF species at the later stages. Moreover, certain EMF indicators were common across tree species at specific development stages: the sapling stage was specifically abundant with Agaricales, Inocybaceae, Tylosporaceae and Amphinema ; the juvenile stage with Russula and Pyronemataceae; and the adult stage with Russula font - queri , Cortinariaceae and Suillaceae (Fig. , Supplementary Fig. S7). 

After quality control and denoising, we obtained 1,153,982 high-quality concatinated nonchimeric reads from 3,875,730 raw paired-end reads across 77 root samples. These sequences were mapped to 5,726 amplicon sequence variants (ASVs), of which 1,108 ASVs (579,294 sequences) were identified as EMF species (Supplementary Table ). These EMF ASVs spanned major phyla of Basidiomycota and Ascomycota, families such as Thelephoraceae, Russulaceae, Sebacinaceae, and Inocybaceae (Supplementary Fig. ). The dominant EMF genera (> 1%) were Tomentella , Russula , Inocybe , Sebacina , Tuber , Cortinarius , Amphinema , Humaria , Cenococcum , Piloderma , Trichophaea , Lactarius , Leccinum , Phaeohelotium , and Thelephora (Supplementary Fig. ). 

We compared the alpha diversity (such as Observed number of ASVs, Shannon, Simpson, Chao1, ACE, and Pielou indices) differences of root EMF by Kruskal‒Wallis test to discern the patterns of fungal richness as influenced by the maturation of hosts (Supplementary Table S3, Fig. ). Furthermore, the EMF alpha diversity indices (Observed number of ASVs, Chao1, ACE) were significantly influenced by the tree species at particular developmental stages ( P < 0.001, Supplementary Table S3). For example, from sapling to juvenile stage, EMF diversity remained relatively stable in Acer mono ( P > 0.05), marginally decreased in Betula platyphylla (0.05 < P < 0.1), whereas marginally increased in Quercus mongolica (0.05 < P < 0.1). Moreover, at the sapling stage, EMF diversity was significantly higher in Betula platyphylla than in Acer mono and Quercus mongolica ( P < 0.05). At the juvenile stage, EMF diversity ranked as Acer mono < Betula platyphylla < Quercus mongolica ( P < 0.05). Unlike those at the early stages, the EMF diversity in the adult hosts did not significantly vary across tree species ( P > 0.05, Supplementary Fig. S4). The above results implicate that as hosts mature, the EMF diversity increases and transitions from early-stage heterogeneity to late-stage convergence among the different host species. 

Nonmetric multidimensional scaling (NMDS) and Permutational multivariate analysis of variance (PERMANOVA) demonstrated significant compositional shifts across stages ( P < 0.001, Fig. a), with adult communities showing lower dissimilarity among tree individuals ( P < 0.05, Fig. b; Supplementary Fig. S5). Particularly, unlike Acer mono and Quercus mongolica , the EMF communities in Betula platyphylla were not significantly different between the sapling and adult stages ( P > 0.05, Fig. ). Our findings indicate that the variance in the root EMF communities among the different hosts decrease markedly with host maturity, revealing a trend toward age-related convergence in the root EMF community composition. 

To explore the interaction patterns within EMF communities, we constructed a meta-community cooccurrence network based on Spearman’s correlation ( P < 0.05) and examined how these networks varied with host maturity and tree species. We identified stage-specific EMF taxa and categorized them into four clusters: sapling-, juvenile-, and adult-enriched ASVs, along with a set of common ASVs (Fig. ). Cooccurrence network demonstrated high degrees of positive correlations among EMF taxa (positive correlations exceeding 99%, Supplementary Table S4) in all three tree species: 6,290 correlations (edges) among 501 ASVs (nodes) for Acer mono , 7,471 edges and 631 nodes for Betula platyphylla , and 11,983 edges and 540 nodes for Quercus mongolica . The modularity values were greater than 0.4 (Supplementary Table S4), indicating that the networks exhibited distinct modular structures. The node-level topological features (degree, closeness centrality, and eigenvector centrality) were significantly greater at the adult stage compared to earlier stages ( P < 0.05, Supplementary Fig. S6). Furthermore, the adult-stage subnetworks encompassed more nodes and edges, higher average degree and graph density, but lower average path length and diameter, relative to its earlier stages (Supplementary Table S4). These findings demonstrate that the EMF communities are more aggregated and cohesive with host maturity, occupying a central position in cooccurrence networks. 

To identify key ASVs driving EMF community assembly, we conducted linear discriminant analysis (LDA) effect size (LEfSe) to detect EMF indicators of specific host species and stage. As for Acer mono , the adult stage was predominated by 20 EMF species primarily belonging to Atheliales and Pilodermataceae; the earlier stages were enriched in Pyronemataceae, Boletales and Russula vinosobrunneola . As for Quercus mongolica , the adult stage were dominated by 19 EMF species predominantly belonging to Atheliales, Pilodermataceae, Thelephora sp., and Cenococcum sp.; the sapling and juvenile stages were closely associated with Thelophora sp. and Elaphomyces papillatus . Betula platyphylla exhibited 11 EMF species enriched at the sapling stage and 7 EMF species at the later stages. Moreover, certain EMF indicators were common across tree species at specific development stages: the sapling stage was specifically abundant with Agaricales, Inocybaceae, Tylosporaceae and Amphinema ; the juvenile stage with Russula and Pyronemataceae; and the adult stage with Russula font - queri , Cortinariaceae and Suillaceae (Fig. , Supplementary Fig. S7). 

To explore the dynamics of EMF communities in coexistent broadleaf trees at a regional scale, we analyzed 77 root tips from a chronosequence encompassing sapling, juvenile, and adult stages in temperate forests of Northeast China. EMF alpha diversity was significantly higher in adult trees than their younger counterparts. Furthermore, EMF community composition diverged to a greater extent during early host development, with mature trees hosting more converged EMF community composition. Cooccurrence network analysis confirmed that, as trees mature, EMF taxa became increasing integrated, displaying higher centrality, connectivity, and modularity. This suggests that mature trees foster more robust mutualistic networks, enhancing symbiotic interactions. Mature broadleaf trees tended to host EMF communities with lower alpha diversity but higher structural dissimilarity compared to saplings and juveniles, consistent with previous studies , . This can be attributed to the limited contact with EMF spores in undeveloped root systems of early-stage hosts . These young pioneers also face higher mortality rates due to environmental stressors such as nutrient depletion, moisture stress, herbivore damage, and wind disturbance . However, mature trees benefit from fully extended root systems that connect them to a broader range of EMF species, leading to more stable nutrient and water acquisition and increased EMF diversity – . These long-term obligate symbioses may streamline EMF partners and optimize nutrient-use-efficiency, facilitating adaptations to similar biotic and abiotic stressors over the host’s life span . Across premature host stages, the stochastic EMF recruitment , compounded by environmental stressors, disrupt the continuity of the mycorrhizal networks, increasing the divergence of the EMF communities . While we did not directly examine whether the EMF community of younger trees is nested within that of mature individuals, this remains a plausible hypothesis for future investigation. Network analysis demonstrated increased EMF integration in mature hosts, with higher modularity and centrality suggesting stabilized and interconnected mutualistic networks. This pattern aligns with prior findings on the stability of adult-hosted EMF communities , , as long-lived trees selectively recruit EMF to synchronize growth and facilitate community convergence. Even though stage-specific EMF did not always directly interact, they were consistently linked to common EMF species, especially in adult hosts. This suggests that the evolution of stage-specific EMF communities heavily rely on common species that mediate the development of mutualistic networks during host maturation . We found that positive associations among EMF species predominated (over 99%), which likely stems from niche overlap, ecological compatibility, and resource exchange. This promotes plant-soil feedbacks, supports plant maturation, and contributes to forest succession . The minimal negative correlations between stage-specific and common EMF species suggest low competition, facilitating the equilibrium and coexistence of EMF species . Several EMF taxa demonstrated broad host spectra across multiple developmental stages. For instance, Suillus is known as an obligate symbiont of Pinaceae and can significantly alters the morphology and hormone content in host seedlings , and it also benefits the growth of Quercus . The mixed forest composition at the current sites supports the likelihood of multihost colonization, highlighting a trend towards symbiotic promiscuity rather than specificity. Furthermore, the host preference of Suillus underscores the evolutionary constraints of EMF species , as its persistence in mature hosts reflects legacy effects of past host community dynamics . Moreover, the adult-stage indicator Plioderma , found in Acer mono and Quercus mongolica , is adapted to low nitrogen availability and typically associated with late-stage hosts – . This indicates that trees at very late developmental stages acquire resources via specific EMF species, thereby supporting their persistence across different succession stages . EMF species adopt different life strategies at various stages of host development. For example, Pyronemataceae, adapted to nutrient-poor environments , are suitable partners for Acer mono saplings that thrive in highly heterogonous niches . Additionally, Inocybe and Amphinema , preferring nutrient-rich and disturbed environments, aid the establishment of Betula platyphylla in post-disturbance sites . In contrast, Russula prefers mature forests, probably due to challenges with spore dispersal and germination at early developmental stages of hosts , . In addition, we also observed that Acer mono adults were predominated by Phaehelotium and Hebeloma , and Quercus mongolica adults by Cenococcum and Thelephora . All these families are optimized for short-distance exploration , maximizing hydrophilic hyphal extension to efficiently absorb nutrients and enabling energy-conserving late-stage growth. 

In summary, we revealed root EMF dynamics and succession across the various development stages of broadleaf trees in Changbai Mountain. The transition from early-stage heterogeneity to late-stage homogeneity illustrates adaptive strategies of symbiontic relationships as host trees mature. This suggests a natural selection process driving the convergence of symbiotic networks, which is crucial for the survival and success of heterospecifics in temperate forests. This adaptive mechanism helps select specific EMF partners as keystone taxa to the hosts at each stage. The EMF–host relationship is not strictly species-specific but exhibits overlap as hosts mature. This highlights the succession of EMF and contributes to understanding how these interactions evolve within forest ecosystems. These insights lay a framework for further investigations on the mechanisms underlying EMF dynamics and implications for the management and protection of forests. 

Experimental site The study was conducted in the National Nature Reserve on the northern slope of Changbai Mountain (41°43′ N to 42°26′ N, 127°42′ E to 128°17′ E), northeastern China. This region expriences a temperate continental monsoon climate, with an annual average temperature of 2.8 °C (January average: -13.7 °C; July average: 19.6 °C) and annual precipitation of approximately 700 mm, mainly occurring in June and August . The reserve, established in 1960, spans 200,000 ha and ranges in elevation from 740 m to 2691 m. Five vegetation zones are present: aspen-white birch ( Populus davidiana and Betula platyphylla ) forests, broadleaf Korean pine ( Pinus koraiensis ) mixed forests, spruce-fir ( Picea jezoensis and Abies nephrolepis ) forests, subalpine birch ( Betula ermanii ) forests, and alpine tundra , and the broadleaf Korean pine mixed forests, predominantly composed of Pinus koraiensis , Tilia amurensis , Quercus mongolica , Fraxinus mandshurica , and Acer mono . The three tree species investigated were: Acer mono (Sapindaceae), at an elevation of 850 ± 50 m; Betula platyphylla (Betulaceae), at an elevation of 900 ± 50 m; and Quercus mongolica (Fagaceae), at an elevation of 950 ± 50 m. These species grow in broadleaf Korean pine mixed forests, secondary poplar–birch forests, and secondary coniferous-broadleaved mixed forests, respectively. Tree classification and sampling To examine how tree ontogeny influences EMF communities, trees were classified into three developmental stages based on diameter at breast height (DBH) as a proxy for tree age , : saplings (DBH < 5 cm), juveniles (DBH 5–20 cm), and adults (DBH > 20 cm) (Supplementary Table ). Sampling was conducted in July 2019 at three distinct sites, with each site representing each host identity with all three developmental stages (Supplementary Fig. ) , . At each site, 10 individual trees per developmental stage were randomly selected, with at least 20 m between individuals. Fine roots were excavated in three directions along the taproots, and the roots from all directions were combined into a single sample for each individual tree. All the saplings were removed from the soil because of the limited number of root systems. A total of 90 root samples (three tree species × three stages × ten replicates) were collected. All saplings were removed due to limited root systems. Samples were transported on ice to the laboratory and stored at -20 °C before processing. DNA extraction Before molecular analysis, we washed away any residual soil clinging to the fine roots under running tap water and cut them into small root segments. Approximately 400 root tips were randomly obtained from each sample under a stereomicroscope and ground into powder with the aid of liquid nitrogen. We selected all the root tips of the sampled saplings (< 200) as a result of their undeveloped root systems. Total root genomic DNA was extracted from 0.5 g of powder via an FH Plant DNA Kit following the manufacturer’s instructions (Beijing Demeter Biotech Co., Ltd., Beijing, China). The quality and concentration of DNA were measured via the optical density (OD) at 260/280 nm on a Thermo Scientific NanoDrop 2000 spectrophotometer, and values ranging from 1.8 to 2.0 were selected for downstream polymerase chain reaction (PCR) amplification. PCR amplification and illumina sequencing The internal transcribed spacer 2 (ITS2) was amplified from the extracted DNA via the primers ITS3tagmix (5′-CATCGATGAAGAACGCAG-3′) and ITS4NGS (5′-TCCTSCGCTTATTGATATGC-3′) . PCRs were performed in 25-µL solutions containing 2.5 µL of DNA template, 12.5 µL of KAPA HiFi HotStart ReadyMix (Thermo Scientific, Waltham, MA, USA) and 1 µM of each primer. The PCR program began with initial denaturation at 95 °C for 3 min, followed by 25 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min. Under the same reaction conditions, we replaced the extracted DNA with nuclease-free water as a negative control to ensure successful amplification. Nextera XT Index primers (Illumina, San Diego, CA, USA) with N7 (S5) were employed as the forward (reverse) reads for the second round of PCR amplification. The 50-µL reaction mixtures comprised 5 µL of DNA, 10 µL of PCR-grade water, 25 µL of KAPA HiFi HotStart ReadyMix, 5 µL of N7xx primers and 5 µL of S5xx primers. The reaction was performed at 95 °C for 3 min, followed by 8 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min. The PCR products obtained were cleaned with AMPure magnetic beads (Beckman Coulter Life Sciences, Indianapolis, IN, USA). The purified PCR products were pooled at equimolar concentrations for amplicon library preparation and sequencing, which were completed on an Illumina MiSeq platform (Illumina, San Diego, CA, USA) at the Testing and Analysis Center, Institute of Applied Ecology, Chinese Academy of Sciences. Sequence processing Raw paired-end reads were analyzed via QIIME2 2020.11 . The DADA2 was used for quality filtering, denoising, chimera removal, and ASVs generation. Low-quality sequences were trimmed fro the 3’ ends, with forward and reverse sequences cut at position 280 and 220, respectively. Taxonomic assignment of fungal sequences was performed against the UNITE + INSDC fungal ITS database (version 8.3) . EMF taxa were identified according to genuslevel functional annotations from the Fungaltraits database . Due to uneven sequence depths between samples, rarefaction was used to normalize each ample to a minimum depth of 24,576 sequences, a threshold chosen to ensure statistical power and representativity s , . The final dataset consisted of 77 samples, after excluding low-quality samples. Raw sequencing data were deposited into the NCBI Sequence Read Archive (SRA) database (BioProject ID: PRJNA1051329). Statistical analysis Taxonomic assignments were visualized via a cladogram through the head_tree function in the metacoder package . The effects of developmental stage and tree species on the EMF community were tested using generalized linear mixed models implemented via the glmer function in the lme4 package (v1.1-35) . Alpha diversity indices (Chao1, ACE, Shannon, Simpson, and Pielou) were estimated with estimateR and diversity functions from the vegan package (v2.6-4). NMDS with Bray‒Curtis dissimilarity distances (metaMDS function) visualized community composition similarity. PERMANOVA (adonis2 function with 999 permutations) tested stage-driven compositional differences. The EMF diversity indices and community compositional differences among various developmental stages were examined as follows. Normality was evaluated via Shapiro-Wilk tests (shapiro_test in rstatix v0.7.2); Homoscedasticity via Levene’s test (levene_test). Where assumptions were violated (e.g., Simpson index: W = 0.872, p = 0.012), non-parametric Kruskal-Wallis tests (kruskal_test) with Dunn’s post-hoc comparisons (dunn_test) replaced one-way ANOVA. All p-values were adjusted via Benjamini-Hochberg false discovery rate control . Cooccurrence networks of the EMF communities were constructed for each tree species with the igraph package . The nodes represented the EMF ASVs, the edges represented correlations between EMF ASVs. Only robust correlations with Spearman’s correlation coefficients (r) > 0.6 or < -0.6 and Benjamini − Hochberg false discovery rate (FDR)-adjusted P values < 0.01 were retained to establish networks . To explore stage-discriminatory patterns, ASVs unique to a given stage were defined as stage-enriched ASVs. Pairwise paired t-test were applied to compare the relative abundance of each ASV among the different developmental stages, with no-significant differences indicating common ASVs. To reflect network aggregation and cohesion, we calculated network-level metrics, including average degree, clustering coefficient, and graph density . The average path length was derived to measure network connectivity and dispersion . Furthermore, the diameter was used to assess the network size and scale . At the node level, we calculated degree, betweenness, closeness, and eigenvector centrality. Degree denotes the number of connections per node; betweenness centrality measures a node’s role as a bridge; closeness centrality represents the average shortest path between nodes; and the eigenvector centrality indicates the influence of a node based on its neighbors’ centrality . Cooccurrence networks were visualized with the interactive platform Gephi 9.0 . To describe and verify distinct species at specific developmental stages, LEfSe analysis was adopted to discover stage-associated EMF indicators via classification comparison, biological consistency testing, and effect size estimation . First, we used the Kruskal‒Wallis test to identify EMF taxa with significant abundance differences across different stages. The Wilcoxon rank sum test method was subsequently employed to determine whether these significantly different taxa converged to the same classification within each stage. Finally, the LDA was applied, with an LDA score threshold of 4 and a significant threshold of P < 0.05 to identify stage-specific indicators. LEfSe analysis was conducted using the Galaxy platform . 

The study was conducted in the National Nature Reserve on the northern slope of Changbai Mountain (41°43′ N to 42°26′ N, 127°42′ E to 128°17′ E), northeastern China. This region expriences a temperate continental monsoon climate, with an annual average temperature of 2.8 °C (January average: -13.7 °C; July average: 19.6 °C) and annual precipitation of approximately 700 mm, mainly occurring in June and August . The reserve, established in 1960, spans 200,000 ha and ranges in elevation from 740 m to 2691 m. Five vegetation zones are present: aspen-white birch ( Populus davidiana and Betula platyphylla ) forests, broadleaf Korean pine ( Pinus koraiensis ) mixed forests, spruce-fir ( Picea jezoensis and Abies nephrolepis ) forests, subalpine birch ( Betula ermanii ) forests, and alpine tundra , and the broadleaf Korean pine mixed forests, predominantly composed of Pinus koraiensis , Tilia amurensis , Quercus mongolica , Fraxinus mandshurica , and Acer mono . The three tree species investigated were: Acer mono (Sapindaceae), at an elevation of 850 ± 50 m; Betula platyphylla (Betulaceae), at an elevation of 900 ± 50 m; and Quercus mongolica (Fagaceae), at an elevation of 950 ± 50 m. These species grow in broadleaf Korean pine mixed forests, secondary poplar–birch forests, and secondary coniferous-broadleaved mixed forests, respectively. 

To examine how tree ontogeny influences EMF communities, trees were classified into three developmental stages based on diameter at breast height (DBH) as a proxy for tree age , : saplings (DBH < 5 cm), juveniles (DBH 5–20 cm), and adults (DBH > 20 cm) (Supplementary Table ). Sampling was conducted in July 2019 at three distinct sites, with each site representing each host identity with all three developmental stages (Supplementary Fig. ) , . At each site, 10 individual trees per developmental stage were randomly selected, with at least 20 m between individuals. Fine roots were excavated in three directions along the taproots, and the roots from all directions were combined into a single sample for each individual tree. All the saplings were removed from the soil because of the limited number of root systems. A total of 90 root samples (three tree species × three stages × ten replicates) were collected. All saplings were removed due to limited root systems. Samples were transported on ice to the laboratory and stored at -20 °C before processing. 

Before molecular analysis, we washed away any residual soil clinging to the fine roots under running tap water and cut them into small root segments. Approximately 400 root tips were randomly obtained from each sample under a stereomicroscope and ground into powder with the aid of liquid nitrogen. We selected all the root tips of the sampled saplings (< 200) as a result of their undeveloped root systems. Total root genomic DNA was extracted from 0.5 g of powder via an FH Plant DNA Kit following the manufacturer’s instructions (Beijing Demeter Biotech Co., Ltd., Beijing, China). The quality and concentration of DNA were measured via the optical density (OD) at 260/280 nm on a Thermo Scientific NanoDrop 2000 spectrophotometer, and values ranging from 1.8 to 2.0 were selected for downstream polymerase chain reaction (PCR) amplification. 

The internal transcribed spacer 2 (ITS2) was amplified from the extracted DNA via the primers ITS3tagmix (5′-CATCGATGAAGAACGCAG-3′) and ITS4NGS (5′-TCCTSCGCTTATTGATATGC-3′) . PCRs were performed in 25-µL solutions containing 2.5 µL of DNA template, 12.5 µL of KAPA HiFi HotStart ReadyMix (Thermo Scientific, Waltham, MA, USA) and 1 µM of each primer. The PCR program began with initial denaturation at 95 °C for 3 min, followed by 25 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min. Under the same reaction conditions, we replaced the extracted DNA with nuclease-free water as a negative control to ensure successful amplification. Nextera XT Index primers (Illumina, San Diego, CA, USA) with N7 (S5) were employed as the forward (reverse) reads for the second round of PCR amplification. The 50-µL reaction mixtures comprised 5 µL of DNA, 10 µL of PCR-grade water, 25 µL of KAPA HiFi HotStart ReadyMix, 5 µL of N7xx primers and 5 µL of S5xx primers. The reaction was performed at 95 °C for 3 min, followed by 8 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 30 s, and a final extension at 72 °C for 5 min. The PCR products obtained were cleaned with AMPure magnetic beads (Beckman Coulter Life Sciences, Indianapolis, IN, USA). The purified PCR products were pooled at equimolar concentrations for amplicon library preparation and sequencing, which were completed on an Illumina MiSeq platform (Illumina, San Diego, CA, USA) at the Testing and Analysis Center, Institute of Applied Ecology, Chinese Academy of Sciences. 

Raw paired-end reads were analyzed via QIIME2 2020.11 . The DADA2 was used for quality filtering, denoising, chimera removal, and ASVs generation. Low-quality sequences were trimmed fro the 3’ ends, with forward and reverse sequences cut at position 280 and 220, respectively. Taxonomic assignment of fungal sequences was performed against the UNITE + INSDC fungal ITS database (version 8.3) . EMF taxa were identified according to genuslevel functional annotations from the Fungaltraits database . Due to uneven sequence depths between samples, rarefaction was used to normalize each ample to a minimum depth of 24,576 sequences, a threshold chosen to ensure statistical power and representativity s , . The final dataset consisted of 77 samples, after excluding low-quality samples. Raw sequencing data were deposited into the NCBI Sequence Read Archive (SRA) database (BioProject ID: PRJNA1051329). 

Taxonomic assignments were visualized via a cladogram through the head_tree function in the metacoder package . The effects of developmental stage and tree species on the EMF community were tested using generalized linear mixed models implemented via the glmer function in the lme4 package (v1.1-35) . Alpha diversity indices (Chao1, ACE, Shannon, Simpson, and Pielou) were estimated with estimateR and diversity functions from the vegan package (v2.6-4). NMDS with Bray‒Curtis dissimilarity distances (metaMDS function) visualized community composition similarity. PERMANOVA (adonis2 function with 999 permutations) tested stage-driven compositional differences. The EMF diversity indices and community compositional differences among various developmental stages were examined as follows. Normality was evaluated via Shapiro-Wilk tests (shapiro_test in rstatix v0.7.2); Homoscedasticity via Levene’s test (levene_test). Where assumptions were violated (e.g., Simpson index: W = 0.872, p = 0.012), non-parametric Kruskal-Wallis tests (kruskal_test) with Dunn’s post-hoc comparisons (dunn_test) replaced one-way ANOVA. All p-values were adjusted via Benjamini-Hochberg false discovery rate control . Cooccurrence networks of the EMF communities were constructed for each tree species with the igraph package . The nodes represented the EMF ASVs, the edges represented correlations between EMF ASVs. Only robust correlations with Spearman’s correlation coefficients (r) > 0.6 or < -0.6 and Benjamini − Hochberg false discovery rate (FDR)-adjusted P values < 0.01 were retained to establish networks . To explore stage-discriminatory patterns, ASVs unique to a given stage were defined as stage-enriched ASVs. Pairwise paired t-test were applied to compare the relative abundance of each ASV among the different developmental stages, with no-significant differences indicating common ASVs. To reflect network aggregation and cohesion, we calculated network-level metrics, including average degree, clustering coefficient, and graph density . The average path length was derived to measure network connectivity and dispersion . Furthermore, the diameter was used to assess the network size and scale . At the node level, we calculated degree, betweenness, closeness, and eigenvector centrality. Degree denotes the number of connections per node; betweenness centrality measures a node’s role as a bridge; closeness centrality represents the average shortest path between nodes; and the eigenvector centrality indicates the influence of a node based on its neighbors’ centrality . Cooccurrence networks were visualized with the interactive platform Gephi 9.0 . To describe and verify distinct species at specific developmental stages, LEfSe analysis was adopted to discover stage-associated EMF indicators via classification comparison, biological consistency testing, and effect size estimation . First, we used the Kruskal‒Wallis test to identify EMF taxa with significant abundance differences across different stages. The Wilcoxon rank sum test method was subsequently employed to determine whether these significantly different taxa converged to the same classification within each stage. Finally, the LDA was applied, with an LDA score threshold of 4 and a significant threshold of P < 0.05 to identify stage-specific indicators. LEfSe analysis was conducted using the Galaxy platform . 

Below is the link to the electronic supplementary material. Supplementary Material 1 Supplementary Material 2 

