In this study, rhizosphere soil samples were collected from both normal roots and cluster roots of Macadamia integrifolia across four major production regions in Yunnan Province–Changning (CN), Yingjiang (YJ), Lancang (LC), and Yunxian (YX)–during both the dry and rainy seasons. A total of 80 samples were obtained (4 regions × 2 seasons × 2 root types × 5 replicates), thereby establishing a multidimensional sampling framework incorporating seasonal, geographical, and root–type variation. These orchards exhibited distinct geographic and climatic characteristics (Figure 1; Table 1). As demonstrated in Table 1, each site exhibited significant seasonal fluctuations. For instance, in Changning, the air temperature increased from 10 to 17.4 °C in the dry season to 17.5–20.8 °C in the rainy season, while precipitation rose from 248 mm to 995 mm. Greater shifts were also observed in Yingjiang (14.2–21.8 °C and 186 mm vs. 22.5–24.4 °C and 1,366 mm), Lancang (13.2–21.8 °C and 195 mm vs. 23.5 °C and 1,429 mm), and Yunxian (13.8–24.4 °C and 273 mm vs. 26.9 °C and 1,094 mm). These seasonal contrasts were layered over clear spatial gradients. During the dry season, Changning exhibited the lowest temperatures, while Yunxian demonstrated the most extensive temperature range and the highest levels of precipitation. During the rainy season, Yunxian experienced the highest temperatures, while Lancang received the most precipitation, and Changning remained predominantly dry. The range of elevation was from 992 m in Lancang to 1,278 m in Yunxian. Collectively, these climatic and geographic gradients furnish a natural experimental framework for assessing how environmental variation influences the composition and dynamics of the macadamia rhizosphere mycobiome.

Sampling was conducted during two distinct seasons: the dry season (DS, November to April) and the rainy season (RS, May to October). In each orchard, sampling plots were meticulously selected based on flat terrain, uniform soil properties, and consistent cultivation and management practices. Each plot encompassed an area of 500 m × 500 m. All plots were exclusively planted with Macadamia integrifolia trees (cultivar O. C.) of similar age (≥8 years). At each site, rhizosphere soil was collected from two root types–normal roots (NR) and cluster roots (CR)–with five biological replicates per group, totaling 80 samples. Rhizosphere soil was collected by gently shaking roots to remove loose soil, then brushing adhering soil (2–5 mm from root surface) into sterile tubes. This follows standard protocols to ensure collection of root-influenced soil (Gong et al., 2024). The identification of root types was conducted in situ through the evaluation of tree morphological characteristics. Trees that were not exhibiting visible indications of disease were meticulously selected for the sampling process. A five-point method was employed to randomly collect soil tightly adhering (2–5 mm) to the root surface. Samples were promptly transferred to 50 mL sterile tubes, flash-frozen in liquid nitrogen, and stored at −80 °C for subsequent analysis. The implementation of sampling procedures was standardized across all regions.

The procedures for DNA extraction and sequencing were carried out in accordance with previously established protocols (Gong et al., 2024). Briefly, genomic DNA was extracted from rhizosphere soil samples using the CTAB/SDS method, visualised on 1% agarose gels, and diluted to a final concentration of 1 ng/μL. The fungal internal transcribed spacer (ITS1–ITS2) region was amplified using the primers ITS1F (5′–CTTGGTCATTTAGAGGAAGTAA–3′) and ITS2R (5′–GCTGCGTTCTTCATCGATGC–3′) with Phusion® High–Fidelity PCR Master Mix under standard conditions. The presence of PCR products was confirmed on 2% agarose gels. Subsequently, the samples were combined in equimolar concentrations, and the purification process was executed using the Qiagen Gel Extraction Kit. The preparation of sequencing libraries was undertaken using the NEBNext® Ultra II DNA Library Prep Kit, and the quality of the libraries was subsequently assessed with a Qubit 2.0 Fluorometer and an Agilent Bioanalyzer 2,100. Paired–end sequencing (250 base pairs) was performed on the Illumina NovaSeq platform.

Paired–end reads were demultiplexed using unique sample barcodes and trimmed to remove barcode and primer sequences. Merged reads were generated using FLASH (v1.2.11), and quality filtering was performed with fastp (v0.20.0) to obtain high–quality tags. The detection and removal of chimeric sequences was conducted using Vsearch (v2.15.0) against the UNITE fungal database, yielding a set of effective tags. Amplicon sequence variant (ASV) inference and denoising were conducted using DADA2 within QIIME2 (v2020.6), followed by taxonomic and functional annotation. Subsequently, the ASV abundance tables were normalized to the sample with the lowest sequencing depth. Alpha diversity indices, including Observed_features, Chao1, Shannon, Simpson, Goods_coverage, and Pielou_e, were calculated in QIIME2 (Bolyen et al., 2019). Beta diversity was assessed using Bray–Curtis distances, and community composition was visualized via principal coordinate analysis (PCoA) and non–metric multidimensional scaling (NMDS). To assess the impact of season, region, and root type on fungal community structure, permutational multivariate analysis of variance (PERMANOVA) was conducted, leveraging the Adonis function with 999 permutations.

The top ten most abundant fungal phyla and genera across all samples were used to generate stacked bar plots to visualize community composition trends. The identification of differentially abundant taxa was achieved through the implementation of Linear Discriminant Analysis Effect Size (LEfSe) with a threshold of LDA > 3.5 and p < 0.05, as well as through Random Forest modeling. The key genera were determined on the basis of their classification contribution and mean importance scores. The similarity percentage (SIMPER) analysis was employed to identify genera that contributed most to community dissimilarities across seasons and regions. The FUNGuild database was employed to perform functional prediction, thereby enabling the categorization of fungal taxa into ecological roles such as saprotrophic, pathogenic, or symbiotic. A comparative analysis was conducted, with the functional profiles (Modes) and ecological guilds (Guilds) being visualised.

A co–occurrence network was constructed based on Spearman correlation coefficients (ρ > 0.6 and p < 0.01) among ASVs in order to identify key interactive taxa within the rhizosphere fungal community. The calculation of network topology metrics, incorporating node degree and betweenness centrality, was conducted using the igraph package in R (v1.2.6). Subsequently, the network was rendered visually using Gephi. Within the network, nodes are used to represent individual ASVs, with edge colors indicating the direction of correlation (red indicating positive and blue indicating negative).

The differences in alpha and beta diversity between groups were analyzed using the Kruskal–Wallis test in QIIME2. Furthermore, bioinformatic analysis was conducted using the OmicStudio tools, accessible at https://www.omicstudio.cn/tool (Lyu et al., 2023).

After stringent quality control (Supplementary Table 1), an average of 79,799 high-quality reads per sample were retained (86.19% of raw reads, range: 62.31–96.39%). Q20 and Q30 scores exceeded 98.78 and 95.56%, respectively, with GC content (39.06–57.55%) consistent with fungal ITS regions. These metrics confirm the reliability of sequencing data for subsequent analyses.

A total of 6,563 fungal ASVs were identified across all samples. Venn diagram analyses revealed substantial differences among environmental groups (Figure 2). A mere 11.21% of ASVs were shared between the DS and RS, while 35.27 and 53.51% were unique, respectively. In a similar manner, CR and NR samples exhibited a mere 17.83% of ASVs in common. Within each site, combinations of season and root type also exhibited high ASV specificity. The findings indicate that seasonality, root morphology, and site conditions interact to shape a highly specialized and context–dependent rhizosphere mycobiome.

Alpha diversity varied across regions, seasons, and root types (Figures 3A,B), with 6,563 fungal ASVs identified. Chao1 (83.75–478.96) and Shannon (1.24–6.15) indices showed near-significant differences (Kruskal–Wallis, Chao1: p = 0.097; Shannon: p = 0.051), suggesting interactive effects of season, region, and root type. Rainy-season CR samples had higher diversity (Chao1: 285.3 ± 124.7; Shannon: 3.8 ± 1.4) than dry-season CR samples (Chao1: 172.4 ± 63.5; Shannon: 2.9 ± 1.1). Regional specificity was evident: e.g., DSYX_CR samples had higher Chao1 than NR, while RSYJ_CR samples had lower Shannon than NR. RS samples from Yingjiang and Lancang showed higher diversity, likely due to increased nutrient availability with precipitation, whereas extreme dry-season samples had reduced diversity, suggesting drought-driven community simplification.

PERMANOVA indicated season (p = 0.001, 3% variance), location (p = 0.001, 12% variance), and their interaction (p = 0.001, 5% variance) significantly shaped fungal community structure, while root type and its interactions did not (p > 0.05; Table 2). PCoA and NMDS (Stress = 0.10) confirmed these patterns, with clear separation of DS vs. RS samples and distinct clustering by location, while root-type effects were minimal (Figure 4). These results highlight season and geography as primary drivers of β-diversity.

Taxonomic analysis (Figure 5) showed Ascomycota dominated all groups (55–65%) at the phylum level, with Basidiomycota (20–30%) and Mortierellomycota (5–10%) as secondary phyla. Seasonal and spatial variations were more pronounced at lower taxonomic levels: e.g., RS samples had higher abundances of Basidiomycota genera (Agaricus, Hebeloma) and Ascomycota genera (Fusarium, Acremonium), while DS samples were enriched in Mortierella and Penicillium (Ascomycota). Yunxian (YX) had higher Basidiomycota, and Lancang (LC) had more Mortierellomycota. At the genus level, cluster roots (CR) were enriched in symbiotic taxa (e.g., Glomus-related), while normal roots (NR) had more saprophytes (e.g., Trichoderma). These patterns indicate stability at the phylum level but environmental responsiveness at finer scales.

Differential analyses identified key taxa associated with season and location (Figure 6). LEfSe (LDA > 3.5, p < 0.05) showed dry-season enrichment of Mortierella and Solicoccozyma, and rainy-season enrichment of Vishniacozyma, Filobasidium, and Fusarium (pathogenic; Figure 6A). Random Forest confirmed Penicillium, Hypocrea, Trichoderma, and Fusarium as top seasonal predictors (Figure 6B). SIMPER analysis revealed region-specific drivers: e.g., Mortierella and Penicillium drove seasonal differences in Changning; Fusarium and Acremonium distinguished CR vs. NR; Metarhizium and Aspergillus increased in rainy seasons in Lancang and Yunxian (Figure 6C).

Co-occurrence network analysis revealed environment-specific interaction patterns (Figure 7; Table 3). Dry-season networks were dominated by Ascomycota (Talaromyces, Penicillium) with high degree and betweenness centrality, while rainy-season networks featured Cladosporium (Ascomycota) and Mortierellaceae (Mortierellomycota) as keystone taxa. Regionally, keystone species varied: e.g., Phialomyces in Changning, Cladosporium and Talaromyces in Lancang, and Penicillium with Acremonium in Yunxian. These patterns highlight environmental shaping of network structure.

FUNGuild predictions (Figure 8) showed saprotrophic fungi dominated (50–55%), with a 10% higher proportion in rainy vs. dry seasons. Guild analysis revealed Unassigned and Soil_Saprotroph as primary groups, with seasonal fluctuations in functional groups like Animal_Pathogen–Plant_Pathogen_Soil_Saprotroph. Functional patterns (e.g., Pathotroph–Saprotroph–Symbiotroph) and guilds (e.g., Endo_Mycorrhizal, Wood_Saprotroph) varied with region, season, and root type: e.g., higher Animal_Pathogen–Plant_Pathogen_Soil_Saprotroph in DSCN_CR, and prominent Soil_Saprotroph in RSYX_NR. These patterns indicate environmental factors drive functional adaptation of rhizosphere fungi.

Our results strongly support the hypothesis that season and geographic location are primary drivers of rhizosphere fungal community structure. PERMANOVA analysis revealed significant effects of season (p = 0.001, 3% variance explained), location (p = 0.001, 12% variance explained), and their interaction (p = 0.001, 5% variance explained) on community structure, whereas root type alone had no significant effect (p = 0.072; Table 2). These findings align with the conclusion that climatic variables dominate fungal community composition (Bui et al., 2020) but extend it by identifying specific mechanisms through our data:

Seasonal effects: Alpha diversity indices (Chao1: 285.3 ± 124.7; Shannon: 3.8 ± 1.4) were significantly higher in rainy-season CR samples compared to dry-season CR samples (Chao1: 172.4 ± 63.5; Shannon: 2.9 ± 1.1; Figure 3), consistent with increased soil moisture and root exudates (e.g., sugars, amino acids) in the rainy season stimulating fungal proliferation (Jiang et al., 2023). Taxonomic analysis further showed that rainy-season samples were enriched in saprotrophic taxa such as Fusarium and Acremonium (Figures 5, 6), validating our prediction that high humidity accelerates organic matter decomposition and supports saprotroph growth. In contrast, dry-season dominance of Mortierella (Figure 6A) likely reflects adaptation to drought via accumulation of osmolytes like proline (Kaushal, 2019).

Regional effects: Geographic variation in dominant phyla—e.g., higher Mortierellomycota in Lancang (LC; low elevation, 992 m; high rainfall, 1,429 mm) and higher Basidiomycota in Yunxian (YX; high elevation, 1,278 m; Figure 5)—directly demonstrates environmental filtering of fungal taxa. SIMPER analysis identified region-specific drivers, such as increased Metarhizium in rainy-season Lancang (Figure 6C), reinforcing the role of local microclimates in shaping communities (Tao et al., 2020).

Our results partially support the hypothesis that cluster roots selectively enrich symbiotic fungi, with enrichment dependent on environmental conditions. While root type alone had no significant effect (PERMANOVA, p = 0.072), cluster roots (CR) showed 1.8- and 2.3-fold higher relative abundances of Glomus and Trichoderma, respectively, compared to normal roots (NR; Figure 5), consistent with the prediction that cluster roots recruit beneficial taxa via exudates (Zhao et al., 2021).

Notably, this enrichment was context-dependent: Chao1 indices were significantly higher in dry-season CR samples from Yunxian (DSYX_CR) compared to DSYX_NR (Figure 3A), but no such difference was observed in rainy-season Lancang (RSLC), likely due to dilution of root exudates by heavy rainfall (Figure 2B). This suggests cluster root effects are amplified in seasons/regions where root exudates remain concentrated (e.g., dry-season Yunxian), supporting our hypothesis that environmental conditions modulate root-fungus interactions.

Key taxa (Fusarium, Mortierella, Trichoderma, Penicillium) identified by LEfSe and Random Forest act as environmental indicators. Rainy-season enrichment of Fusarium (pathogenic) likely reflects favorable growth conditions under high humidity (Popovski and Celar, 2013), while dry-season dominance of Mortierella (drought-tolerant saprotroph) indicates adaptation to resource-limited conditions (Caravaca et al., 2024). FUNGuild analyses supported the hypothesis that functional guild composition shifts with seasonal and regional changes in community structure. Saprotrophs dominated in both seasons but were 10% more abundant in the rainy season (50–55% vs. dry season; Figure 8), consistent with increased organic matter decomposition under high moisture (Tanunchai et al., 2023). Conversely, dry-season samples showed higher proportions of multi-trophic guilds (e.g., Animal_Pathogen–Plant_Pathogen–Soil_Saprotroph; Figure 8B), reflecting fungal adaptation to resource scarcity via facultative parasitism.

Key taxa identified by LEfSe–e.g., rainy-season enrichment of Fusarium (pathogen) and Trichoderma (biocontrol agent; Figure 6A)–further validated functional differentiation under high humidity. Dry season dominance of drought–tolerant Mortierella (saprotroph; Figure 6A) aligned with our prediction of stress–driven functional specialization. Additionally, 35% of edges in rainy season co–occurrence networks were negative (Figure 7), indicating intensified resource competition among saprotrophs, directly linking community structure to functional dynamics.

Co-occurrence networks revealed environmental shifts in fungal interactions. Dry-season networks were dominated by Ascomycota (Talaromyces, Penicillium) with high connectivity, forming mutualistic modules (e.g., with Trichoderma) likely involved in cellulose decomposition and nutrient cycling (He et al., 2017, 2025; Samson et al., 2011). Rainy-season networks shifted to Cladosporium and Mortierellaceae, with 35% negative interactions indicating heightened resource competition—likely due to rapid organic matter decomposition limiting carbon availability. Regional variations in keystone species (e.g., Phialomyces in Changning, Cladosporium in Lancang) highlight geographic shaping of network structure (Liu et al., 2024; Zhang et al., 2017, 2018). Penicillium, a prevalent hub across networks, may act as a ‘resident’ core taxon, warranting functional validation.