Rhizosphere fungal community assembly varied across functional guilds in a temperate forest

Rhizosphere fungi play an important role in plant community dynamics and biogeochemical cycling. While the drivers of fungal community assembly have been studied in varied ecosystems, it is still unclear how these processes function for rhizosphere soil fungi in temperate forests. Furthermore, it is unknown whether the relative contributions of important determinants remain consistent or vary across fungal ecological guilds. This study used high-throughput next-generation sequencing to characterize the fungal communities of 247 rhizosphere soil samples from 19 tree species in a temperate forest within Northeast China. We aimed to investigate how three important determinants in temperate forests (host tree species, neighbouring plant communities, and edaphic properties) influence the community assembly of fungal functional guilds in the rhizosphere soil of trees. We found that host tree species contributed more to plant pathogens’ community composition than ectomycorrhizal fungi, and plant pathogens consistently showed higher host specialization than ectomycorrhizal fungi. Saprotrophs also showed high host specialization, which was mediated by the tree species’ effect on rhizosphere soil pH. Although neighboring plant communities contributed remarkably to richness of all fungal guilds, this effect on fungal composition varied across functional guilds, with stronger effect for biotrophic guilds (plant pathogens and ectomycorrhizal fungi) than for non-biotrophic guild (saprotrophs). Neighboring plant communities shaped the ectomycorrhizal community composition strongly in all samples regardless of host trees’ mycorrhizal type, whereas edaphic properties were the most important drivers for this guild in samples from only ectomycorrhizal-associated trees. Edaphic properties played an important role in shaping ectomycorrhizal and saprotrophic fungal compositions, indicating the importance of edaphic properties on the fungal functional guilds associated with the absorption and decomposition of nutrients. These results demonstrated that rhizosphere soil fungal community assembly determinants varied across fungal guilds, reflecting their different ecological functions in temperate forest ecosystems.


Introduction
Forest ecosystems harbour a large diversity of soil fungi that regulate plant community dynamics (Averill et al. 2014;Molina and Horton 2015) and biogeochemical processes (van der Heijden 2008; Clemmensen et al. 2013). Fungi drive many classical ecological phenomena in forests, e.g., conspecific density dependence (Janzen 1970;Connell 1971;Chen et al. 2019) as well as the observed relationships between soil fertility and plant community structure (Mao et al. 2019). Plants recruit and cooperate with soil fungi, especially mutualistic and pathogenic fungal species, which influence the fitness and competitiveness of host plants (Peay 2016;Kandlikar et al. 2019) and ultimately determine host plant survival and growth. Plants may also alter nutrient availability and soil fungal communities in their rhizosphere. These plant-soil feedback loops are important for plant populations and forest community dynamics (Bennett et al. 2017;Fujii et al. 2018). Therefore, exploring the local drivers of rhizosphere fungal community assembly in natural forests is required to better understand how these microorganisms influence plant-soil feedbacks, and ultimately predict how natural forest communities develop and respond to environmental change.
Plant pathogens, mycorrhizal fungi and saprotrophs are the most abundant fungal guilds that occur in soils and in particular in plants rhizosphere (i.e., the soil compartment that is directly influenced by the root system, henceforth, we call the rhizospheric fungi). Plant pathogens cause diseases and typically accumulate around adult trees by colonising conspecific or phylogenetically related heterospecific species, thereby facilitating the establishment of other species (Connell 1971). Therefore, plant pathogens play an important role in promoting the coexistence of diverse tree species (Connell 1971;Jia et al. 2020). Mycorrhizal fungi form mutualistic relationships with plants associating directly with their roots enhancing water and nutrient absorption and increasing host resistance to pathogens while in exchange fungi receive carbon photosynthesised by the plants (Smith and Read 2008). As the primary decomposers of organic matter, saprotrophs can regulate physicochemical cycles in the soil and influence plant growth by increasing belowground soil mineral nutrient content (Dighton et al. 1987). Generally speaking, plant pathogenic fungi and mycorrhizal fungi are part of the biotrophic guilds that depend on living plants to survive. In contrast, saprotrophic fungi are non-biotrophic, meaning they acquire all their nutrients from non-living plant material or surrounding soils. Due to their ecological and functional dissimilarities, these fungal guild communities may be shaped by different determinants (Nguyen et al. 2016b;Yang et al. 2019).
Several biotic and abiotic mechanisms can shape fungal community assemblies, such as host phylogeny , neighbouring plant communities (Chagnon et al. 2020;Cheng and Yu 2020), and edaphic properties (Glassman et al. 2017;Schappe et al. 2017). Edaphic properties, such as soil pH and nutrition, are considered primary factors that impact fungal community assembly (Dumbrell et al. 2010;Glassman et al. 2017;Davison et al. 2021). Host plants have also been shown to affect rhizosphere soil fungal communities (Becklin et al. 2012;Huang et al. 2014;Sweeney et al. 2021), since they directly interact with a myriad of fungi in the root rhizosphere. Likewise, rhizosphere fungal community might be indirectly influenced by plant-mediated changes in edaphic properties, such as altered soil pH in proximity to plant roots (York et al. 2016) and altered soil temperatures by litter cover (Weltzin et al. 2005). Moreover, neighboring plant communities can influence the rhizosphere fungal communities by constructing the above-and belowground environments around the hosts, where local light conditions, litter input, and secretion of chemicals are regulated (Badri et al. 2009;Neuenkamp et al. 2020;Kong et al. 2021). In addition, neighboring plants may act as fungal nurseries (Facelli et al. 2018) or alternative carbon sources (Moeller et al. 2015), which helps maintain the local fungal species pool. While the drivers of fungal community assemblage have been extensively studied, the relative importance of these determinants has been scarcely explored locally in natural forest ecosystems. Exploring contributions of these determinants to fungal community assembly in a natural forest ecosystem could further elucidate the mechanisms behind soil fungal community assembly. These studies would provide a better understanding of how trees interact with their plant neighbors and the surrounding soil by influencing rhizosphere fungi with different ecological functions. Previous studies have observed host preferences in plant pathogens, mycorrhizal fungi, and saprotrophs (Zhou and Hyde 2001;Wang et al. 2019). These non-random associations between fungi and their plant hosts are important mechanisms implicit in fungal niche partitioning (Dickie 2007), ultimately influencing fungal species richness and community composition. However, the degree of fungal specialization for tree species likely varies between fungal functional guilds and is influenced by the evolutionary strength of the host-microbe interaction. Long-term coevolution of direct interactions between plants and their associated pathogenic and mycorrhizal fungi has created strong host specializations of biotrophic fungi (Adamson and Caira 1994;Field and Pressel 2018). Saprotrophic fungal host specialization might be shaped by different rhizosphere habitat environments, since studies have shown the effect of plant species on soil properties (Calvaruso et al. 2011) and litter characteristics (Makita and Fujii 2015). Therefore, the assessment of fungal host specialization can help shed light on the host-dependent mechanisms behind fungal community assembly.
Although the factors shaping soil fungal communities have been widely investigated across ecosystems (Peay et al. 2013;Sweeney et al. 2021) and from local to global scales (Tedersoo et al. 2014;Glassman et al. 2017), relatively few studies have focused on the drivers of rhizospheric fungal communities associated with various co-occurring tree species across a steep gradient in edaphic properties at a local scale in natural temperate forests. Thus, this study aims to assess the importance of the host tree species, neighbouring plant communities, and edaphic properties in shaping tree rhizosphere fungal communities in a natural temperate forest. We collected rhizosphere soil samples from 19 major tree species in a temperate forest in northeast China, characterized the fungal communities, and analyzed the edaphic properties. We also retrieved information on the taxonomy and phylogeny of the host tree species and the characteristics of their neighbouring plant communities. With this information, we aimed to address the following questions: 1. What is the relative contribution of host tree species, neighbouring plant communities, and edaphic properties to the richness and composition of plant pathogens, ectomycorrhizal fungi and saprotrophs? We hypothesised that biotrophic guilds (i.e., plant pathogens and ectomycorrhizal fungi) would be mainly affected by host tree species and neighbouring plant communities, whereas non-biotrophic guilds (i.e., saprotrophs) are affected mainly by edaphic properties. 2. How does fungal specificity to tree species vary among the three functional guilds? We expected biotrophic guilds to exhibit higher specialization than non-biotrophic guilds. Moreover, we expected that phylogenetic relationships among host trees would shape the two biotrophic guilds more strongly, which would explain their higher host specialization.

Site description and soil sampling
The experiment was conducted in the Changbaishan Nature Reserve (42.38°N, 128.08°E) in northeast China. This region has a typical monsoon climate, with a mean annual precipitation of 700 mm and a mean annual temperature of 2.8 °C. Our research took place in a 25-ha broad-leaved Pinus koraiensis mixed forest dynamics plot (500 × 500 m and approximately 300 years). The plot elevation ranged from 791.8 m to 809.5 m (Wang et al. 2015). There were obvious dominant species in the plot, including Pinus koraiensis, Tilia amurensis, Quercus mongolica, and Fraxinus mandschurica (Additional file 1: Table S1). The understory also had a dense and diverse assemblage, including herbs and mosses. The herb assemblage was dominated by Meehania urticifolia, Anemone amurensis, and Cardamine leucantha (Jia et al. 2022). Most of these herbaceous species are known to form arbuscular mycorrhizas (AM), and other species do not associate with mycorrhizal fungi. This plot was established in 2004, where all free-standing trees with ≥ 1 cm in diameter at breast height (dbh) of 1.3 m have been mapped, identified to species, and measured every 5 years. Since 2009, more than 1300 individual trees with dbh ≥ 5 cm belonging to 27 tree species (representing all tree species in the plot) have been under growth monitoring with dendrometer bands. We collected rhizosphere soil from July 10 to August 10, 2019, from the 19 most abundant tree species under growth monitoring with dendrometer bands (Additional file 1: Table S1). We randomly selected 8-15 individuals per species and tried to maximize the range of their dbh within species, ensuring that the sampling sites were in different soil nutrient gradients (Mao et al. 2019; Additional file 1: Fig. S1). For each tree selected, loose debris was removed, and fine roots were excavated with a shovel by tracing from the thicker main roots carefully. The loose soil around the fine roots was gently removed, and only approx. 10-20 g soil firmly attached to the fine roots was collected using a brush. Soil samples from three different sub-sampling locations around each tree were collected and pooled together. The sub-sampling locations were always within 10 cm below the soil surface and within 10 m around the targeted tree. The shovel used for sampling was sterilized with 75% alcohol before each sample collection, and new disposable gloves and brushes were used between samples. We performed the aforementioned experimental operations in the field and all the samples were cooled with ice, brought back to the laboratory and stored at − 80 °C until further processing. To determine the soil physical and chemical properties, we collected approximately 200 g of bulk soil closely adjacent to the location, where the rhizosphere soil was collected. The soil samples used to determine soil physical and chemical properties were air-dried for approximately 1 month before subsequent analysis.

Measurements of soil physicochemical properties
Nine soil properties were measured according to Lu (1999): soil pH, soil organic matter, available nitrogen (N), phosphorus (P) and potassium (K), total N, P and K, dissolved organic carbon (Additional file 1: Table S2). The soil C/N ratio was calculated using the soil organic matter and total N values generated in this study.

DNA extraction and molecular analysis
Total DNA was extracted from 0.25 g of frozen soil using a FastDNA ® Spin Kit for Soil (MP Biomedicals, Inc. USA), according to the manufacturer's instructions. Primers ITS1F (5′-CTT GGT CAT TTA GAG GAA GTAA-3′; Gardes and Bruns 1993) and the ITS2 reverse primer (5′-GCT GCG TTC TTC ATC GAT GC-3′; Op De Beeck et al. 2014) were used to amplify the fungal ITS1 region of rRNA gene. Both forward and reverse primers were tagged with adapter, pad, linker, and 6 base pair (bp) barcode sequences. Polymerase chain reaction (PCR) amplification was carried out using a 20 μL reaction solution containing 4 μL of 5 × TransStart FastPfu buffer, 2 μL of dNTPs (2.5 mM), 0.8 μL of the forward primer (5 μM), 0.8 μL of reverse primer (5 μM), 0.4 μL of TransStart FastPfu DNA Polymerase and 10 ng of template DNA. Thermal cycling included an initial denaturation step at 95 °C for 1 min, then 35 cycles of 95 °C for 30 s, 55 °C for 30 s, and 72 °C for 45 s, followed by a final extension at 72 °C for 10 min. PCR products were sequenced on the Illumina MiSeq platform (Illumina Inc., San Diego, CA, USA).

Bioinformatic analysis
The raw gene sequencing reads were demultiplexed and quality-filtered using the software fastp (version 0.20.0). First, raw reads were truncated at any site with an average quality score < 20 over a 50 bp sliding window. The high-quality reads longer than 50 bp were resorted into sample groups based on the library prep barcoding, allowing up to 2 nucleotides to mismatch. The demultiplexed reads were merged by the software FLASH (version 1.2.7). Only overlapping sequences longer than 10 bp were assembled according to their overlapped sequence with the maximum mismatch ratio of 0.2. The ITS1 regions of the merged sequences were extracted using the ITSX software package (ITSx_1.1.2; Bengtsson-Palme et al. 2013). Subsequently singleton reads or reads with a length of < 100 bp were removed. The remaining reads were checked for chimeras using vsearch uchimeref in Quantitative Insights into Microbial Ecology software package (QIIME2; Caporaso et al. 2010), referencing the unified system for DNA-based fungal species linked to the classification (UNITE) database. This procedure led to 8,798,773 high-quality ITS1 sequences clustered into 12,224 distinct operational taxonomic units (OTUs) at a 97% sequence similarity threshold in de novo mode using function vsearch cluster-featuresde-novo in QIIME2. The most abundant sequence of each ITS1-based OTU was chosen as the representative FASTA sequence for BLAST against the international nucleotide sequence databases collaboration and UNITE database Version 8.2 (Abarenkov et al. 2020). Among the 247 total samples, we excluded one sample with < 10,000 ITS1 sequences and rarefied the remaining 246 samples to the minimum depth of 11,106 sequences to reduce potential amplicon sequencing biases using the function diversity core-metrics in QIIME2. According to Tedersoo et al. (2014) and the FUNGuild database (Nguyen et al. 2016a), fungal OTUs were assigned to functional guilds based on their genus-level taxonomic annotation. In the FUNGuild database, we only used the annotations with confidence ranking "Highly Probable" and "Probable". If one OTU was assigned to multiple guilds by FUNGuild, we classified it into "Others" in our study. We assigned all Glomeromycota as arbuscular mycorrhizal. Among these functional guilds, we listed plant pathogens, arbuscular mycorrhizal fungi (AMF), ectomycorrhizal fungi (EMF) and saprotrophs for subsequent analyses.
Other guilds, such as endophytes, animal pathogens, or lichenized fungi accounted for a very small proportion of our data set and were categorized as "Others". Sequencing data were deposited in the NCBI under the accession number PRJNA834583.

Statistical analysis
All statistical analyses were conducted in R v.3.6.1 (R Core Team 2020). Data sets were logarithmically or square-root transformed before analysis to improve our model residuals' normality and variance homogeneity. To assess the effect of tree species on edaphic variables, we used analysis of variance (ANOVA) to calculate the variation in edaphic variables captured by tree species. The ANOVA was executed using the lm function from the stats package (version 3.6.1). The rarefaction curves of each fungal guild were calculated using the specaccum function in the vegan package (version 2.5-6; Oksanen et al. 2013) to check whether our samples gave a typical representation of the entire soil fungal community at the sampling site. Fungal species richness was defined as the number of OTUs per sample. To address the variation in fungal composition, we calculated pairwise Bray-Curtis dissimilarity distances among samples (Faith et al. 1987). The sample-OTU matrix underwent Hellinger transformation beforehand (Legendre and Gallagher 2001) to downweight rare OTUs to analyze fungal composition. We visualized the dissimilarities in each fungal composition using nonmetric multidimensional scaling (NMDS) ordination performed by the metaMDS function in the vegan package.
The relative contributions of three sets of determinants (tree species, neighboring plant communities, and edaphic properties) to rhizosphere fungal community assembly were explored using variance partitioning. The most important explanatory variables were chosen using forward selection with the forward.sel function from the adespatial package (version 0.3-8; Dray et al. 2018). Then, the variation in fungal richness was partitioned by the three determinants using function varpart from the vegan package. Based on distance-based redundancy analysis (db-RDA; Legendre and Anderson 1999), this particular variance partitioning approach was applied to the fungal community composition. Principal coordinate analysis (PCoA) was conducted on the pairwise Bray-Curtis dissimilarity distances matrix using the cmdscale function from the stats package (version 3.6.1) to represent dissimilarities in the the fungal community composition. Next, as described earlier for species richness, the variance partitioning method was implemented. A variance partitioning analysis was also used to quantify the richness and composition of EMF when they occurred in the samples only from ectomycorrhizal (EM) trees.
The explanatory variables describing the three sets of determinants (i.e., tree species, neighbouring plant communities and edaphic properties) of the fungal community (i.e., species richness and community composition) were assembled as follows: (1) to assess how tree species could influence a fungal community, we used tree species as a categorical variable; (2) Given the heterogeneity and old age of the sampled forests, we defined an area within a 40 m radius of each target tree as the "neighborhood area" to assess the effect of neighboring plant communities. Neighbouring plant information for each sampled tree was extracted from the census data of our 25-ha forest dynamics plot in 2019. We used the crowding index to represent the conspecific, heterospecific AM, and heterospecific EM adult tree densities. These indexes were calculated based on basal area divided by the distance between conspecific, heterospecific AM, or heterospecific EM adult neighbours and the focal trees: where i is an individual tree; A can be conspecific, heterospecific AM, or heterospecific EM neighbouring effect. BA is the basal area of conspecific, heterospecific AM, or heterospecific EM individual tree. DISTANCE represents the distance between the conspecific, heterospecific AM, or heterospecific EM neighbour and the focal tree. The neighboring variables also included tree species richness, the PCo1 and 2 axes of the Hellinger transformed neighboring tree abundance matrix representing neighboring species composition, the PC1 and PC2 axes of community-level functional trait composition based on 13 species-level functional traits representing neighboring A = N 1 BA i /DISTANCE i functional traits, and two community-level phylogenyrelated indices: the Nearest Relative Index (NRI) and the Nearest Taxon Index (NTI) in neighboring communities (Additional file 1: Table S2); (3) Nine variables: soil pH, soil organic matter, available nitrogen (N), phosphorus (P), potassium (K), total N, P, K, dissolved organic carbon and the soil C/N ratio were used to represent edaphic properties. Each variable was scaled between 0 and 1 using the formula (x-x min )/(x max -x min ) to standardize their effects and make them comparable. Pairwise collinearity between variables was assessed using Pearson's correlation coefficient (r), with one of the pairs being removed from downstream analyses when r > 0.8. Collinearity was only tested between variables within hostrelated, neighboring, and soil-related groups but not among these groups to account for the possibility that different factors could have a partial collinear effect.
Interaction networks between target host trees and rhizosphere soil fungi were constructed to test whether fungal specialization to tree species varied among the three functional guilds. A plant-fungal interaction matrix was built for each fungal guild, respectively, with rows representing tree species and columns representing fungal OTUs, where each cell represented the frequency of observed associations between tree species and fungal OTUs. Given the low possibility of an association between EMF and non-ectomycorrhizal trees, only EM tree species were included in the plant-EMF interaction networks. The standardized Kullback-Leibler distance (d') index was used to quantify the degree of interaction specialization at the species level between fungal OTUs and tree species (Blüthgen et al. 2006) using the dfun function from the bipartite package (version 2.15; Dormann et al. 2008). d' represented the degree of specialization to tree species (i.e., fungal host specialization), and from the trees' point of view, d' represented the degree of their specialization to fungal OTUs. The index fell between 0 (lowest specialization) and 1 (highest specialization). We obtained 1,000 randomized interaction matrices by shuffling the species labels of the observed tree individuals and proceeded to calculate the standardized d' value of each tree species and fungal OTU as follows: d' standardized = [d' observed − Mean (d' randomized )]/SD (d' randomized ), where d' randomized are the d' values of the randomized interaction matrices. The Mann-Whitney test (Bauer 1972) was used to detect whether the standardized d' index of fungal OTUs significantly differed among the three guilds.
Finally, to test whether fungal species richness and composition were structured by tree species phylogenetic relationship, we calculated Blomberg's K value (Blomberg et al. 2003) and Pagel's λ value (Pagel 1999) using the phylosignal function in the phylosignal package (version 1.3; Keck et al. 2016) with the phylogenetic tree generated in a previous study conducted in the same field .

Edaphic properties
We observed large ranges in most edaphic variables in the tree rhizosphere (Table 1). For example, the value of pH ranged from 4.00 to 7.72, and organic matter ranged from 90.37 g/kg to 508.37 g/kg. For some edaphic variables, relatively large variation could be captured by tree species. The variables with highest variation explained by tree species were soil TP (ANOVA, r 2 = 55%, F = 17.72, P < 0.01), C:N ratio (ANOVA, r 2 = 46%, F = 12.63, P < 0.01) and soil pH (r 2 = 30%, F = 6.9, P < 0.01).

Fungal community taxonomical and functional composition
After the rarefication of the 246 samples, we obtained 5848 fungal OTUs, including 2,051 OTUs (35.1%) from the phylum Ascomycota and 2716 OTUs (46.4%) from Basidiomycota. The majority of the fungal OTUs (72.3%) were categorized into specific guilds successfully: 230 (3.9%) were plant pathogens, 1691 (28.9%) were ectomycorrhizal fungi, 1878 (32.1%) were saprotrophs, and 37 (0.63%) were arbuscular mycorrhizal fungi (Additional file 1: Fig. S3). Due to the low number of sequences assigned to AM fungi, we did not include this guild in subsequent downstream analysis. Our sample-based species accumulation curves for each guild did not reach asymptotes (Additional file 1: Fig. S2), suggesting that the entire community for each guild was not fully covered by our sampling. However, the assessment of differences in fungal species richness and composition in relation to host tree species, neighbouring plant communities, and edaphic properties was representative, because the rarefaction curves of different fungal guilds reached similar degrees of saturation (Additional file 1: Fig. S2).

Relative contributions of determinants to fungal richness and composition
The variables we used explained a relatively small proportion of the variance in fungal richness (< 40% Fig. 1A-C). For plant pathogens, edaphic properties alone explained 3.4% of the variations in richness, more than tree species or the neighbouring plant communities (Fig. 1A). For ectomycorrhizal fungi and saprotrophs, host tree species and neighbouring plant communities contributed more to the variation in fungal richness than edaphic properties (Fig. 1B, C). Especially for ectomycorrhizal fungi, host tree species and neighbouring plant communities contributed to most of the variation in fungal richness compared to the other fungal guilds, with 11% and 5.4% of the variation separately and 12.8% jointly. Regarding the ectomycorrhizal fungi sampled from EM trees, only the neighbouring plant communities explained a significant variation in fungal richness (10%; Fig. 3A). The relative importance of the three determinants of variation in the fungal community composition differed greatly among the three functional guilds (Figs. 1D-F, 2B). The plant pathogenic fungal community composition variation was primarily explained by host tree species (6.8%). The neighbouring plant communities explained most of the variation in ectomycorrhizal fungal community composition (9.9%), whereas saprotrophic fungal community composition variation was primarily explained by edaphic properties (24.4%). Edaphic properties alone explain 13.6% of the variation in EMF community composition associated with EM trees (Fig. 2B). However, the largest amount of the variation in the EMF community composition associated with EM trees (15.2%) was  (Fig. 2B).

Specialization in tree-fungi associations
Host specialization varied significantly between fungal guilds, with higher or similar specialization for plant pathogens (d' standardized -median = 1.83) and saprotrophs (d' standardized -median = 1.76), respectively, and significantly lower specialization for ectomycorrhizal fungi (d' standardized -median = 1.32, Fig. 3). Pinus koraiensis (d' standardize = 6.00) and Fraxinus mandshurica (d' standardize = 5.16) showed higher specialization in their  purple, edaphic properties. Overlapping circles represent the variation explained by those factors jointly. Portion with R 2 < 0.02 was not denoted. Digits in Venn diagrams represent variation explained by each portion; Digits in brackets represent their proportion in total explained variation correspondingly plant pathogens than other tree species. Betula platyphylla (d' standardize = 8.99), P. koraiensis (d' standardize = 7.28), and Tilia mandshurica (d' standardize = 6.69) exhibited greater saprotroph-related specialization. All five EM tree species showed significant higher levels of specialization in their associated ectomycorrhizal fungi, with the most specialized species being B. platyphylla (d' standardized = 7.17).

Discussion
In this study, we explored the rhizospheric fungal community assemblages in a temperate forest. The majority of the fungal functional guilds were identified as ectomycorrhizal and saprotroph. In our analysis, few sequences were identified as arbuscular mycorrhizal fungi, though there are a significant proportion of AM tree species in this temperate forest ecosystem (Mao et al. 2019). The absence of AM fungi could mainly be attributed to the sample-collection methodology we imposed. We collected the soil samples from root surface, but not the plant tissue, where AM fungi reside. The identification limitations of the ITS region on the AMF fungi might also cause this bias. Though the ITS region has the highest probability of successful identification for the broadest range of fungi, ITS may not be appropriate for AMF, and some of the saprotrophic groups (Nilsson et al. 2019). The 18S region would be targeted in future work to get more comprehensive knowledge of the full spectrum of the studied fungal community (Gorzelak et al. 2012).
We aimed to assess the importance of the host tree species, neighbouring plant communities, and edaphic properties shaping tree rhizosphere fungal communities in a natural temperate forest located northeast China. Our results demonstrated that, at the local scale, the relative importance of the ecological factors shaping these fungal communities in these ecosystems can differ across fungal functional guilds. We found evidence for differences in fungi and tree species associations between the two biotrophic guilds (i.e., plant pathogenic and ectomycorrhizal fungi). We observed that the effect of the host tree species on fungal composition and fungal host specialization was consistently higher for plant pathogens than ectomycorrhizal fungi. Although neighbouring plant communities had a common effect on fungal richness for each guild, the neighbouring plant effect was stronger in shaping the two biotrophic fungal compositions than saprotrophs, indicating that biotrophic fungal communities might be more sensitive to changes in the plant communities surrounding their host plant. In addition, neighbouring plant communities shaped the ectomycorrhizal community strongly in all samples, whereas this effect gives way to soil effect when ectomycorrhizal fungi associated with their EM hosts. Edaphic properties Fig. 3 Boxplots showing the distribution of standardized d' index for each fungal guild. The thick bar represents the median, the box represents the 25th and 75th quartiles, and the whiskers represent first and fourth quartiles. Single points are outlying values. PF: plant pathogens; EMF: ectomycorrhizal fungi; SAF: saprotrophs. We only included OTUs that presented at least three samples in any one of 19 tree species (for ectomycorrhizal fungi, 5 EM tree species). d' significantly differed between plant pathogens and ectomycorrhizal fungi, and between ectomycorrhizal fungi and saprotrophs, but not between plant pathogens and saprotrophs (Mann-Whitney test) Fig. 4 Fungal richness (A) and composition (B) are structured by tree phylogeny. Species level data with Pagel's λ and Blomberg's K measures of phylogenetic signal were showed above the respective chart. Fungal richness was represented by the mean richness per species; fungal composition was represented by the mean of first axis scores in db-RDA per species. Colours of tree species names: blue, EM tree species; Red, AM tree specie were the most influential factors shaping the composition of saprotrophs and ectomycorrhizal fungi associated with EM trees, indicating the importance of edaphic properties influencing fungal functional guilds associated with the absorption and decomposition of nutrients.

Tree species effects and fungal specialization
Host tree species contributed more to the composition of plant pathogens than to that of ectomycorrhizal fungi (Fig. 1D, E), corroborating empirical evidence that plant pathogens show stronger host preference than mutualists (Schroeder et al. 2019;Wang et al. 2019). Stronger host specialization was observed in plant pathogens than in ectomycorrhizal fungi, likely due to an underlying mechanism leading to pronounced tree species effects on plant pathogens, as reported in a previous study . The high host specialization of plant pathogens is believed to result from their co-evolutionary pathway which force fungi and plant hosts to adapt to each other constantly in an evolutionary "arms race" relationship to avoid each other's defences (Antonovics et al. 2013). Cooperation between plants and mycorrhizal fungi has depended more on mutual rewards and resource availability (Werner et al. 2015). Hence, the higher specialization of plant pathogens compared to ectomycorrhizal fungi may underline the relevance of co-evolutionary processes in shaping communities of plant pathogens (Gilbert and Webb 2007), whereas the coupling between mycorrhizal fungi and plant partners might be more dependent on the environmental context (Glassman et al. 2017;Ning et al. 2019). While our fungal specialization analysis supported these assumptions, the analysis of phylogenetic signals did not provide evidence for this. We found that ectomycorrhizal fungi rather than plant pathogens showed the strongest phylogenetic structure in the three guilds (Fig. 4), probably because ectomycorrhizal fungi were constrained in their co-occurrence with EM tree species, clustering in plant phylogenetic trees in our study (Fig. 4), In contrast, the associations with host trees were not restricted to a particular group of trees for plant pathogens.
Somewhat to our surprise, saprotrophs also showed a similarly strong specialization to tree species, like plant pathogens (Fig. 3), which is consistent with a recent study on tree rhizosphere fungi in subtropical forests (Chen et al. 2019). Unlike plant pathogens, the strong association between saprotrophs and their preferential tree species may be primarily shaped by the edaphic properties of the rhizosphere. We found B. platyphylla, with a significantly lower rhizosphere soil pH than other tree species, harboured the most specific saprotrophs but not significantly specific plant pathogens (Additional file 1: Fig. S6; Table 2). Since soil pH was the most important factor shaping saprotrophic community composition (Additional file 1: Table S5), B. platyphylla likely recruited saprotrophic fungal species adapted to acidic habitats.

Neighbouring plant communities
Neighboring plant communities had a common effect on fungal richness in all guilds (> 10% of the explained variations, Figs Tables S5, S6). This effect of neighbouring plant communities on fungal richness could be due to neighbouring plant communities acting as nurseries that stabilize local fungal species pools, a phenomenon that has been observed in arbuscular (Mony et al. 2021;Neuenkamp et al. 2021), ectomycorrhizal fungal communities (Peay et al. 2010) and plant pathogens (Hantsch et al. 2014). In this study, we provide further evidence for the contribution of neighbouring plant communities to the saprotrophic fungal richness. Guild-dependent effects of the neighbouring community on fungal composition were consistent with our hypothesis, where we expected that neighbouring plant communities could explain a certain percentage of the variation in two biotrophic fungal guilds but contributed very little to the saprotrophic fungal composition ( Fig. 1D-F). Our results align with previous studies showing that neighbouring plants played an important role in root-associated fungi, such as plant pathogens and ectomycorrhizal fungi (Hubert and Gehring 2008;Facelli et al. 2018;Cheng and Yu 2020). At the local scale, where multiple interactions occurred between neighbouring plant individuals, neighbouring trees were expected to have a marked effect on the assembly of rhizosphere fungal communities (Fichtner et al. 2017). For instance, soil-borne plant pathogens predominantly colonised conspecific or phylogenetic-related heterospecific plant species. Thus, neighbouring trees and their phylogenetic relatedness may have significantly shaped the colonization processes of pathogenic fungal communities from their soil surroundings and their local community composition (Rottstock et al. 2014;Cheng and Yu 2020). This was supported by our results revealing that the effects of conspecific individuals were the primary drivers shaping the communities of plant pathogens (Additional file 1: Table S3).
Neighbouring plant communities contributed substantially to the composition of ectomycorrhizal fungi in samples from all tree species (Fig. 1E), aligning with the idea that the surrounding vegetation shaped the overall pool of mycorrhizal fungi in the soil, from which the roots selected only a subset (Neuenkamp et al. 2021). Conversely, neighbouring community effects disappeared and were replaced by strong effects of edaphic properties when focusing only on the 'suitable' symbionts (i.e., EM fungi in the rhizosphere soil of EM trees) (Fig. 2B). Strong edaphic effects on mycorrhizal fungal communities have been reported in previous studies (Glassman et al. 2017;Van Geel et al. 2018), yet the observed differences in drivers of overall and 'suitable' mycorrhizal fungal communities pointed toward two underlying scale-dependent mechanisms. First, the general availability of mycorrhizal fungi in the rhizosphere of a tree was mediated by the neighbouring plant communities and their effects on the mycorrhizal fungal community (Neuenkamp et al. 2020;Mony et al. 2021). Second, the availability of the subset of 'suitable' symbionts representing the functional fungi was mediated by local conditions, such as the availability of nutrients, where trading of nutrients in exchange for carbon between fungi and their host plants can be expected. This assumption aligned with earlier findings of contextdependent interactions between plants and mycorrhizal fungi (Hoeksema et al. 2010;Johnson et al. 2010) and highlighted multiple scale-dependent, possible nested filtering mechanisms that shaped the biotrophic fungal communities in the tree rhizosphere.

Edaphic properties played an important role in fungal guilds associated with nutrition
Consistent with our hypothesis and reflecting previous studies (Ballauff et al. 2021;Huang et al. 2021), edaphic properties primarily influenced saprotrophic fungal composition compared to the other two guilds (Fig. 1F). Together with the fact that ectomycorrhizal fungal community composition was primarily shaped by edaphic properties (Fig. 2B), our results suggested that when fungi perform functions related to the absorption and decomposition of nutrients, edaphic properties became important for shaping their communities for both biotrophic and saprotrophic guilds.
Saprotrophs play a crucial role in biogeochemical cycles, since they excrete extracellular enzymes to depolymerize recalcitrant lignin and cellulose molecules (Baldrian 2008). The efficiency of their life activities has been shown to be strongly affected by abiotic factors, such as soil moisture (A' Bear et al. 2014) and soil pH (Hobbie and Gough 2004). In addition, soil nutrient content like nitrogen and phosphorus, which are required macronutrients for mycorrhizal fungi, could structure the community of the mycorrhizal guild. Furthermore, soil nutrient levels influence plant investment strategies, subsequently affecting the mycorrhizal fungi who partner with plants by trading nutrients and carbon (Treseder Tables S5; S6). This finding supported the notion that soil pH is an important determinant of microbial community structure (Hobbie and Gough 2004;Suz et al. 2014;Glassman et al. 2017;Ge et al. 2017;Arraiano-Castilho et al. 2021). This remarkable effect of soil pH on fungal communities also reflected the effects of soil nutrient availability, which can be particularly regulated by soil pH (Kluber et al. 2012). B. platyphylla had a significantly lower rhizosphere soil pH than other species, meanwhile, and was the most specialized species for both ectomycorrhizal and saprotrophic fungal communities (Additional file 1: Fig. S6; Table 2). This result suggested that tree species' effects on rhizosphere soil properties may shape the specialization of these two guilds.
Notably, the host tree, neighboring plant communities, and edaphic properties could not explain a high proportion of residuals. The driving factors for how the soil fungal community assembles in a natural forest ecosystem would be very complicated, including soil nutrients, aboveground plant status, and the microclimate caused by topology. The unexplained residuals might be caused by some variables that are not measured in this study (Chomicki et al. 2020). In addition, all the soil samples were collected from trees' rhizosphere but not roots, and the fungal communities without plant tissue may have more randomness than that inside (Yao et al. 2019).
Our results also demonstrated that a considerable component of the variability in fungal communities was represented by the joint effect of multiple determinants. This may be due to the following reasons: (1) edaphic properties affect both tree and fungal community composition in a similar way. (2) The edaphic properties, tree, and fungal community composition are reacting in a similar way to some other variables not taken into account (e.g., temperature, topography, etc.).

Conclusions
This study revealed that different primary determinants were driving the assemblages of plant pathogens, ectomycorrhizal fungi, and saprotrophs, consistent with the ecological characteristics and functions of the fungal guilds. The host tree species effect and host specialization of plant pathogens were higher than those of ectomycorrhizal fungi. Saprotrophs were also strongly shaped by host tree species and showed host specialization, yet this was mediated by the tree species' effect on rhizosphere soil pH. Neighbouring plant communities shaped the pool of available fungi for the host tree and significantly affected biotrophic fungal composition rather than saprotrophic fungal composition. Edaphic properties were the most critical determinants of ectomycorrhizal fungi associated with EM trees and saprotrophs, suggesting that when fungi performed functions related to the absorption and decomposition of nutrients, edaphic properties influenced the assembly of the fungal community.
Additional file 1. Table S1. Information of sampled tree species in our study. Table S2. Details of the variables used in our study. Table S3. Neighbouring and edaphic variables used in variance partitioning analyses for richness of each fungal guild. Table S4. Neighbouring and edaphic variables used in variance partitioning analyses for ectomycorrhizal fungal richness based on samples from 5 EM tree species (69 samples). Table S5.
Neighbouring and edaphic variables used in variance partitioning analyses for composition of each fungal guild. Table S6. Neighbouring and edaphic variables used in variance partitioning analyses for ectomycorrhizal fungal composition based on samples from 5 EM tree species (69 samples). Fig. S1. The individuals (dbh ≥ 5cm) distribution in our plot. Fig. S2 Species accumulation curves of fungi in rhizosphere soil. Fig. S3. Proportion of number of fungal OTUs and sequences assigned to major functional guilds and major fungal phyla. Fig. S4. Nonmetric multidimensional scaling (NMDS) ordination of the Bray-Curtis dissimilarity matrix of fungal communities. Fig. S5. Nonmetric multidimensional scaling (NMDS) ordination of the Bray-Curtis dissimilarity matrix of ectomycorrhizal fungal communities in the samples from EM trees. Fig. S6. Boxplot of variation in edaphic properties among tree species.