Soil Layers Matter: Vertical Stratification of Root-Associated Fungal Assemblages in Temperate Forests Reveals Differences in Habitat Colonization

Ectomycorrhizal and saprotrophic fungi play pivotal roles in ecosystem functioning. Here, we studied the vertical differentiation of root-associated fungi (RAF) in temperate forests. We analysed RAF assemblages in the organic and mineral soil from 150 experimental forest plots across three biogeographic regions spanning a distance of about 800 km. Saprotrophic RAF showed the highest richness in organic and symbiotrophic RAF in mineral soil. Symbiotrophic RAF exhibited higher relative abundances than saprotrophic fungi in both soil layers. Beta-diversity of RAF was mainly due to turnover between organic and mineral soil and showed regional differences for symbiotrophic and saprotrophic fungi. Regional differences were also found for different phylogenetic levels, i.e., fungal orders and indicator species in the organic and mineral soil, supporting that habitat conditions strongly influence differentiation of RAF assemblages. Important exceptions were fungal orders that occurred irrespective of the habitat conditions in distinct soil layers across the biogeographic gradient: Russulales and Cantharellales (ectomycorrhizal fungi) were enriched in RAF assemblages in mineral soil, whereas saprotrophic Polyporales and Sordariales and ectomycorrhizal Boletales were enriched in RAF assemblages in the organic layer. These results underpin a phylogenetic signature for niche partitioning at the rank of fungal orders and suggest that RAF assembly entails two strategies encompassing flexible and territorial habitat colonization by different fungal taxa.


Introduction
Fungi are a remarkably diverse group of organisms on Earth, playing pivotal roles in terrestrial biogeochemical cycles [1,2]. In temperate forests, saprotrophic, symbiotrophic, and pathotrophic fungi have distinct tasks in carbon and nutrient cycling [3]. Saprotrophic fungi are mainly responsible for decomposing plant litter in the forest floor and for maintaining carbon cycling, whereas symbiotrophic fungi such as mycorrhizal fungi colonize roots and mine the soil for mineral nutrients, which they deliver to the plant in exchange for photosynthetically derived carbohydrates [3]. Pathogenic fungi cause diseases and thereby may eventually structure the composition of the vegetation [4]. A small change in these microbial structures can significantly impact matter fluxes and ecosystem functioning [5].
Because of their eminent roles in nutrient cycling in forests, the assembly processes of belowground fungal communities have received increasing attention. In general, fungal assemblages are not occurring at random [6,7]. The composition of soil fungal communities is mainly driven by habitat conditions [8][9][10]. Environmental factors with strong effects on the fungal community composition in soil of temperate forests include temperature, precipitation, soil properties, vegetation, etc., [8,[11][12][13][14].
In addition to forest soil, roots are an important habitat for belowground fungi [15]. The fungi colonize the rhizoplane and grow as saprotrophs or interact with living tissue as endophytes, pathogens or symbionts [15]. Based on their habitat, these fungi have been defined as root-associated fungi (RAF) [16] and are considered as critical components of the plant microbiome [15]. Like soil fungal assemblages, the RAF composition in temperate forests is shaped by various biotic and abiotic environmental conditions including tree species, soil pH, soil moisture, availability of mineral nitrogen and phosphorus and the soil C/N ratio [17][18][19][20][21]. However, along large-scale environmental gradients these factors have lower impact on the composition of RAF than on the composition of fungi in soil [17]. Therefore, it was suggested that the RAF community composition is influenced by additional factors than the soil residing fungi [17], for instance by root properties, especially root N contents [18,20]. Soil depth is a further important factor structuring fungal communities [19,[22][23][24][25][26]. The forest floor consisting of decaying organic material is characterized by high organic carbon (C) and N contents and a wider C/N ratio, whereas the mineral top soil contains narrower C/N ratios and generally decreasing availabilities of nutrient resources [27,28]. The transition from the organic layer to mineral top soil often shows a sharp border with an associated drastic shift in soil chemistry [27,28]. The structuring gradient of soil layers has mainly been studied for soil fungi, showing that saprotrophic fungi were more abundant in the uppermost soil layers composed of fresh litter and decomposing organic material, whereas mycorrhizal fungi dominated deeper in the mineral soil [19,22,29,30]. The influence of soil horizons on fungal composition was even stronger than that of other environmental factors [19,26]. Fine scale analyses of vertical distribution of fungi have also been conducted in different soil depth and showed vertical niche partitioning of ectomycorrhizal fungi on roots [31], fungal hyphae [32], and both [33]. Despite the huge impact of soil horizons on structuring fungal assemblages, studies on the stratification of RAF in organic and mineral soil are scarce and entirely lacking for large-scale biogeographic regions.
The goal of this study was to enhance our understanding of fungal niche partitioning by studying commonalities and differences of the composition of RAF communities in the organic layer and mineral top soil across different biogeographic regions. For this purpose, we used 150 well-established plots of a large-scale infrastructure for biodiversity studies (Biodiversity Exploratories: https://www.biodiversity-exploratories.de/en/, accessed on 10 May 2017) in typical European forests mainly composed of beech (Fagus sylvatica) and conifers (Picea abies, Pinus sylvestris, [34]). The Biodiversity Exploratories are located in the northeast, the middle and the southwest of Germany, encompassing areas of about 450 to 1300 km 2 . The northeast is characterized by drier and warmer climate and sandy soils, whereas the middle and the southwest have lower temperatures, higher soil moisture and silty or loamy soils [34][35][36]. We used this set-up to investigate the fungal assemblages of roots in organic and mineral soil. We know from our previous studies that the tree roots in our forest plots are massively colonized by ectomycorrhizal fungi [37] but how different fungal guilds on roots are affected by soil strata is unknown. Here, we hypothesized that (i) the species richness and relative abundance of symbiotrophic fungi is higher than that of saprotrophic fungi in RAF assemblages, irrespective of the soil layer. We further hypothesized that (ii) the taxonomic community composition of symbiotrophic and saprotrophic fungi differs significantly between organic layer and mineral soil and shows lower turnover for symbiotrophic than for saprotrophic fungi between the soil layers. We hypothesized that (iii) RAF patterns indicate response traits either to soil strata or to regional habitat conditions. To test the third hypothesis, we compared the responses of phylogenetically related taxa (at the rank of fungal orders) to mineral and organic soil in different biogeographic regions. Moreover, we expected to discover root-associated indicator species for distinct soil layers and for different biogeographic regions representing territorial or flexible behaviour.

Study Sites Characteristics
This study was conducted in 150 experimental forest plots in the Biodiversity Exploratories (https://www.biodiversity-exploratories.de/en/, accessed on 10 May 2017, [34]). The Biodiversity Exploratories are located in three geographic regions: Schorfheide-Chorin (SCH) in the northeast, Hainich-Dün (HAI) in the center and Schwäbische Alb (ALB) in the southwest region of Germany. The plots locations are indicated in Supplementary Materials (Table S1). The soil types vary among the regions, with silty soils in ALB, loamy soils in HAI and sandy soils in SCH [36] (Table 1). The plots are located in 53-to 141-year-old forest stands mainly composed of Fagaceae (Fagus sylvatica or Quercus sp.) and Pinaceae (Picea abies and Pinus sylvestris) [34,38]. Additional soil and climatic characteristics have been compiled in Table 1.

Root and Soil Sampling from Organic and Mineral Soil Layers
In each region, roots and soil were sampled in 50 experimental plots [34] in May 2017. In each plot, two transects of 40 m length from north to south and from east to west were established and soil was collected at a distance of 4.5, 10.5, 16.5, 22.5, 28.5, 34.5 and 40.5 m along the transects as described previously [20]. At each sampling point, organic layer (Oe horizon = moderately decomposed organic material with a proportion of visible plants residues between 16 and 66% and Oa horizon = strongly decomposed organic materials with less than 16% visible plant residues) and mineral top soil (to a depth of 10 cm) were sampled separately, resulting in 14 samples per soil layer. In each plot, the samples per layer were combined to one sample of the organic and one sample of the mineral soil. Subsequently, the soil was sieved. Fine roots (<2 mm in diameter) were collected and washed in tap water. Approximately 1 g of fine roots were frozen in liquid nitrogen in the field and stored at −80 • C. The remaining roots were collected, dried (60 • C) and the biomass was recorded (Supplementary Materials Table S1). Aliquots of soil samples (approximately 300 g) were stored at 4 • C.

Determination of Soil Chemical Properties
Carbon and nitrogen: soil samples were dried at 40 • C for two weeks and ground to a homogenous fine powder with a ball mill (Retsch, Type MM400, Haan, Germany). Aliquots of dry soil powder (10 mg organic soil, 30 mg mineral soil) were weighed in in 4 mm × 6 mm tin capsules (IVA Analysentechnik, Meerbusch, Germany) using a super-micro balance (S4, Sartorius, Goettingen, Germany). Total soil carbon (C) and nitrogen (N) were measured by dry combustion in a CN analyzer "Vario Max" (Elementar Analysensysteme GmbH, Hanau, Germany). Acetanilide (71.09% C, 10.36% N) was used as the standard.
Phosphorus and basic cations: potentially plant available phosphorus (P sol ) was extracted according to the method of Bray and Kurtz [39]. Approximately 100 mg of dry soil powder was mixed with 150 mL of Bray extraction solution (0.03 N NH 4 F, 0.025 N HCl). The suspension was shaken slowly for 60 min on a rotary shaker and subsequently filtered through phosphate-free paper filters (MN 280 1/4 125 mm, Macherey-Nagel, Düren, Germany). P sol was measured in the filtrate by inductively coupled plasma-optical emission mass spectroscopy (ICP-OES) (iCAP 7000 Series ICP-OES, Thermo Fisher Scientific, Dreieich, Germany). To determine the potassium (K), calcium (Ca) and magnesium (Mg) contents, approximately 40 to 50 mg soil powder was extracted in 25 mL of 65% HNO 3 (Merck, Darmstadt, Germany) for 12 h at 160 • C [40]. The suspension was filtered (filter papers MN 640 w, width 90 mm, Macherey-Nagel, Düren, Germany) and the filtrate used for element determination by ICP-OES (iCAP 7000 Series, Thermo Fischer Scientific). The calibration was performed with standards of 1 g L −1 (Einzelstandards, Bernd Kraft, Duisburg, Germany). The cations per sample (mmol g −1 dry soil) were added and the sum of cation was used in further analyses.

DNA Extraction and Polymerase Chain Reaction
Frozen fine roots that had been stored at −80 • C were milled with a Retsch ball mill (Type MM400, Retsch GmbH, Haan, Germany) at a frequency of 30 s −1 for 3 min in liquid nitrogen. The genomic DNA was extracted from the frozen root powder with the innuPREP Plant DNA Kit (Analytik Jena AG, Jena, Germany), following the manufacturer's instructions. We used the internal transcribed spacer (ITS) region 2 for fungal identification as recommended by Horton and Bruns [41]. To amplify the fungal ribosomal ITS2 region in roots, we used ITS3_KYO2 [42] as the forward primer (GATGAAGAACGYAGYRAA) and ITS4 [43] as the reverse primer (TCCTCCGCTTATTGATATGC). The primers contained adapter sequences for MiSeq sequencing (Illumina Inc., San Diego, CA, USA). The polymerase chain reactions (PCR) were conducted as described elsewhere [20]. The size of the PCR products was determined after electrophoresis in 2% agarose gels (Biozym Scientific GmbH, Hessisch Oldendorf, Germany) with a DNA standard (Thermo Scientific TM GeneRuler TM 1kb DNA Ladder, Life Technologies GmbH, Darmstadt, Germany). The PCR products were purified with a magnetic bead-based Magsi-NGS PREP kit (Steinbrenner Laborsysteme GmbH, Wiesenbach, Germany). PCR products were stained using GelRed (0.01 µL mL −1 , GelRed TM Nucleic Acid, Biotium Inc., VWR International GmbH, Darmstadt, Germany). The PCR products were visualized with an FLA-5100 Fluorescence Laser Scanner (Raytest GmbH, Straubenhardt, Germany) and an Aida Image Analyser v. 4.27 (Raytest GmbH). Purified PCR products were quantified using a Qubit dsDNA HS assay Kit in a Qubit 3.0 Fluorometer (Thermo Fischer Scientific, Dreieich, Germany). Amplicon sequencing was conducted on the Illumina MiSeq platform using the MiSeq Reagent Kit v3 (Illumina Inc., San Diego, CA, USA) at the Göttingen Genomics Laboratory, Germany.
Applying this pipeline, we obtained a total number of 2.6 and 4.0 million sequences for the root samples from the organic and mineral soil layers, respectively. The OTU table was rarefied to 3890 reads per sample (minimum number of counts in a sample) by employing the function amp_subset_samples() from the 'ampvis2 package implemented in R [53]. Rarefaction resulted in a total number of 2046 and 2147 distinct OTUs for roots from the organic and mineral soil, respectively. Fungal OTUs were functionally annotated as symbiotroph (SYM), saprotroph (SAP) and pathotroph (PAT) using the program 'FunGuild' [54].

Data Processing and Statistical Analysis
The statistical analyses were performed in R version 4.0.3 [55]. Data distribution and homogeneity of the variances were inspected visually using histograms and residual plots. Data were logarithmically transformed to meet the criteria of normal distribution and homogeneity of variances, where necessary. Generalized linear models (Poisson regression, chi-square test) were used with the function glm() from 'lme4 package to investigate datasets with count data (fungal species richness, sequence reads) [56]. Linear models were used to investigate the datasets with continuous data (C, N, C/N, P sol , sum of cations) with the function lm() from 'lme4 package [56]. Differences among variables were calculated using the Anova() function from the 'car' package. Pairwise differences between different groups were compared with a post hoc test (HSD Tukey's honestly significant difference) using the function glht() from the 'multcomp' package [57]. Differences were considered to be significant when p ≤ 0.05.
We investigated the dissimilarities of RAF communities by non-metric multidimensional scaling (NMDS) ordination. We performed NMDS analyses with the function metaMDS() in 'vegan' package [58] with three dimensions, 100 iterations and Bray-Curtis as the dissimilarity measure. Effects of the soil layers on fungal community composition were tested using Analysis of Similarities (ANOSIM with 999 permutation steps) with the function anosim() as implemented in 'vegan' package.
We used "generalized additive models for location scale and shape" (GAMLSS) with a zero-inflated beta (BEZI) family (GAMLSS-BEZI) from the 'metamicrobiomeR' package [59] to compare the relative abundance of the fungal orders between organic and mineral soil layers. Fungal orders were filtered applying the threshold of >0.5% of the mean sequence abundances using the function taxa.filter() to exclude low-abundant fungal orders from further analysis. The relative abundances of the root fungal orders from the organic and mineral soil layers were compared using the function taxa.compare(). The p-values were adjusted using the default method 'False Discovery Rate' (FDR). The analyses were conducted for each of the three study regions separately and then a meta-analysis was performed using the function meta.taxa() to estimate the overall effects across the three regions with "region" as random factor. In the meta-analysis, the model pooled adjusted estimates and standard errors with inverse variance weighting and the DerSimonian-Laird estimator from each region to determine the overall effects across regions. Heatmap and forest plot were generated using the function meta.niceplot(). The heatmaps show log-odds ratio (log(OR)), i.e., effect sizes indicating enrichment or depletion of fungal orders between organic and mineral soil layer. The fungal orders were grouped according to their inferred ecological functions as SYM or SAP. A fungal order was regarded as SYM, when >95% of sequences per order were annotated as SYM or as SAP when >95% of the sequences of the given orders were assigned as SAP. If a fungal order could not be assigned to SYM or SAP, it was classified as multiple functional mode such as SAP+PAT, SAP+SYM, etc.
To determine the β-diversity (Sørensen dissimilarity), the abundance-based OTU data matrices were transformed into presence-absence (1 or 0) matrices with the function beta.temp() from the 'betapart' package [60]. Total β-diversity was further partitioned into species turnover, i.e., replacement of taxa between organic and mineral soil layer and nestedness, i.e., loss of taxa between organic and mineral soil layer using 'betapart'. Pairwise differences of each response variable (e.g., overall beta diversity, nestedness and turnover components) and each fungal group (e.g., all fungi, SYM and SAP) were compared separately among the regions, using a post hoc test (HSD Tukey's honestly significant difference). Further, paired rank sum test was used to compare the different fungal functional group in each region.
Bipartite networks were built to evaluate associations of fungal taxa (OTUs) with organic and mineral soil layers, using the 'bipartite' package [61]. We included OTUs with relative abundances >0.1% of the fungal sequences in this analysis. The plotweb() function was used to visualize bipartite network plots. To determine fungal indicator taxa for roots in the organic and mineral soil layer, we used the multipatt() function from the 'indicspecies' package [62]. Further, the significance difference of the enrichment of the SYM and SAP taxa between the soil layers were tested using the function fischer.test(). Data were visualized using 'ggplot2' package [63] in R.

Differences in Soil Chemistry among Different Biogeographic Regions Are Larger in the Mineral Topsoil Than in the Organic Layer
Across the three regions, which spanned a geographic distance of about >800 km, the organic soil was characterized by higher concentrations of carbon, nitrogen, basic cations and potentially soluble P than the mineral top soil ( Table 2). In the organic layer, C/N ratios exhibited relatively stable values of approximately 22 across the three regions (Table 2). In the mineral soil, C/N ratios varied among the regions and were approximately 30% lower in ALB and HAI than in SCH (Table 2). Mean P sol concentrations varied among the regions, about 1.5-fold in the organic and 2-fold in the mineral soil ( Table 2). The basic cation concentrations varied among the regions, about 3.6-fold in the organic and 7.7-fold in the mineral soil (Table 2). Overall, regional differences were stronger in the mineral top soil than in the organic layer (Table 2). Table 2. Soil chemical properties in different soil layers and regions. ALB = Schwäbische Alb, HAI = Hainich-Dün, SCH = Schorfheide-Chorin. Carbon (C) (g kg −1 DW); nitrogen (N) (g kg −1 DW); ratio of CN; phosphorus (P sol ) (mg kg −1 DW); Cations: sum of potassium (K), calcium (Ca) and magnesium (Mg) (mmol kg −1 DW). Data indicate means ± SE (n = 50). Linear models were used to compare the means of the element between the regions. Significant differences of the means are shown in bold. Pairwise differences of the nutrient elements between organic layer and mineral soil were compared using a post hoc test (HSD Tukey's honestly significant difference). Different letters denote significant differences between soil layers and within the regions.

Strong Taxonomic Differentiation of Root-Associated Fungal Assemblages between Organic and Mineral Topsoil
Our analyses of RAF in 300 root samples from two different soil layers in three biogeographic regions resulted in a rarefied data set containing 1.2 Mio fungal sequences, which clustered into 2537 different fungal OTUs. The mean OTU richness per plot ranged from 158 to 188 for the roots in the organic layer and from 106 to 157 in roots from the mineral top soil (Table S2). For both soil strata and all fungal guilds, lowest fungal richness was found in SCH compared to HAI or ALB (Table S2).
SAP and PAT fungi on the roots exhibited higher OTU richness in the organic than in the mineral layer (Figure 1a-c), while the richness of the SYM fungi was higher in the mineral than in the organic soil (Figure 1a-c). Similarly to OTU richness, the number of SYM sequences was higher on roots in the mineral than in the organic layer, while SAP and PAT sequences were enriched on roots from the organic compared to the mineral layer (Figure 1d-f). Overall, OTU richness of PAT and their relative abundances (3% of the total number of sequences) were low compared to SYM (54%) or SAP (31%) (Figure 1). ). Generalized linear model with Poisson regression and chi-square test was used to compare the fungal groups between soil layers. Pairwise differences of the functional groups of root-associated fungi between organic layer and mineral soil were compared using a post hoc test (HSD Tukey's honestly significant difference). Different letters indicate significant differences of the means at p ≤ 0.05.
In each region, the RAF assemblage in the organic layer was clearly separated from that in the mineral soil (Figure 2a-c), showing dissimilarities of fungal communities on roots in different soil layers. These separations were significant (ANOSIM, p and r values in Figure 2a-c). Significant separations between the organic and mineral soil were also found for SYM and SAP, the main functional fungal groups colonizing roots (SYM: Figure 2d-f, SAP: Figure 2g-i). The goodness-of-fit (r values in Figure 2d-i) was greater for SAP than for SYM assemblages, indicating stronger dissimilarities for SAP than for SYM assemblages on roots in different soil layer. Overall, our data show a strong vertical differentiation of the taxonomic and functional composition of the RAF between organic and mineral soil layers ( Figure 2).

β-Diversity of Root-Associated Fungal Taxa between Organic and Mineral Soil Shows Regional Differences
We found significant differences for the β-diversity of the fungal community composition among the three regions, which were caused by lower β-diversity of RAF in SCH than in the other two regions (Figure 3; Table S3). We partitioned total β-diversity into the turnover component, i.e., replacement of taxa between organic and mineral soil layers, and nestedness, i.e., loss of taxa between organic layer and mineral soil. We found significant differences in RAF turnover among the regions, irrespective of the fungal group analysed with higher turnover in ALB and HAI than in SCH (Figure 3; Table S3). Nestedness of all fungi was higher in SCH than ALB or HAI (Figure 3). The enhanced higher nestedness in SCH was caused by enhanced nestedness of SAP since the SYM nestedness was unaffected ( Figure 3). In general, the extent of the turnover accounted for approximately 80% of the β-diversity of RAF assemblages (Figure 3). We further compared the turnover and nestedness components of SYM and SAP fungi in each region separately. We found significant differences of SYM and SAP fungal turnover and nestedness in HAI (turnover: test statistics = 2.10, p = 0.04; nestedness: test statistics = 3.21, p = 0.001) and SCH (turnover: test statistics = 3.21, p = 0.001; nestedness: test statistics = 3.20, p = 0.001) but not in ALB (turnover: test statistics = 1.65, p = 0.10; nestedness: test statistics = 1.01, p = 0.31). , and black sections represent the mean fungal nestedness (β SNE ). Pairwise differences of each variable were compared separately among the regions using a post hoc test (HSD Tukey's honestly significant difference). Different letters represent significant differences of the means at p ≤ 0.05. Capital letters refer to significant differences for the overall beta diversity and small letters for the turnover and nestedness components of the beta diversity.

Phylogenetically Related Fungal Groups Show Divergent Responses to Organic Layer and Mineral Soil
We classified OTUs according to fungal orders and selected 19 orders, each representing at least 0.5 % or more of the total number of sequences (Table S4). The selected orders accounted together for 85 % of the total number of sequences ( Figure S1). Eight of the 19 fungal orders were assigned as SYM, six as SAP, two as SYM+SAP and three as SAP+PAT guilds (Table S5). In most regions, we found significant differences in the relative abundance of the selected RAF orders between organic and mineral soil layers (Figure 4a). However, the responses did not always show the same direction in different regions (Figure 4a), thus, resulting only in eight significantly affected orders across all study regions (Figure 4b). The orders of Polyporales and Sordariales, which contained only SAP fungi, were enriched on roots from the organic layer compared to those from the mineral soil (Figure 4a,b). Significant enrichment on roots from the organic layer was also observed for the Hypocreales and Pleosporales, orders, which were composed of saprotrophic and pathotrophic fungi (Figure 4). Russulales and Cantharellales, orders containing SYM fungi, were enriched on roots in the mineral soil (Figure 4a,b). In contrast to Russulales and Cantharellales, Boletales, also a SYM order, was enriched on the roots from the organic layer (Figure 4a,b).
Heterogeneous enrichment patterns between different soil layers were observed for SYM fungal orders (Thelephorales, Atheliales, Pezizales and Mytilinidales) as well as for SAP fungal orders (Chaetothyriales and Tubeufiales). In most cases, RAF orders in ALB and HAI showed opposing responses compared to SCH (Figure 4a). . Effect sizes of the changes of root-associated fungal orders between organic to mineral soil. The most abundant fungal orders with >0.5% of the total number of sequences were included. Data are displayed as a heatmap of the log-odds ratio (log(OR)), which indicates the abundance in organic layer relative to that in the mineral soil. Data show the results for three regions (a) and the forest plot of pooled estimates across all three regions with 95% confidence intervals (b). ALB = Schwäbische Alb, HAI = Hainich-Dün, SCH = Schorfheide-Chorin. Red and blue colours indicate the enrichment of the root-associated fungal order in the organic layer and mineral soil, respectively. The order of Cystofilobasidiales was not present in the SCH region. Significant differences at p < 0.05 are denoted with * and at p < 0.001 with **. Significant pooled log-odds ratio (OR) estimates with p < 0.05 are shown in red colour in the forest plot (b).

Fungal Indicator Taxa in Organic Layer and Mineral Soil
Bipartite network association analysis was performed to determine the degree of habitat preference of distinct taxa (OTU based) in the organic layer and mineral soil ( Figure 5). Upper nodes refer to fungal OTUs, and lower nodes refer to organic layer (red) and mineral soil (blue). The width of the upper nodes reflects the relative abundance of fungal OTUs from both soil layers. Line widths represent the relative abundances of the fungal OTUs from the respective soil layer. The statistically significant associations of root-associated fungal OTUs are shown in red for the organic layer and in blue for the mineral soil in the upper nodes. Black nodes represent non-significant fungal OTUs associations between the organic and mineral soil layers.
In ALB, 82 fungal taxa were significantly associated with roots with equal numbers of taxa (41) being significantly associated with the organic and mineral soil (Figure 5a). A similar pattern was observed for HAI (OTUs: organic = 44 and mineral = 43, Figure 5b) but not in the SCH region (OTUs: organic = 31 and mineral = 18, Figure 5c). Across the regions, we found that 30 of 39 SAP indicator species (77%) were present on roots from the organic layer and that 52 of 63 SYM indicator species (81%) were present on roots in the mineral soil layer (Table S6). The enrichment of SYM in the mineral and of SAP indicator taxa in the organic layer was significant (p < 0.001, Fisher's exact test).

Stratification of Root-Associated Fungi by Organic Layer and Mineral Topsoil
Here, we characterized the assembly patterns of fungi associated with roots in two major soil strata, the organic layer and the mineral top soil across a large biogeographic scale. In agreement with studies on soil-localized fungi [22,26,[64][65][66][67][68][69], we found a clear separation of the RAF communities according to soil layers in each of the three biogeographic areas. In line with targeted analyses of mycorrhizal fungi colonizing the roots [19,29,70,71], we found a clear separation of symbiotrophic as well as of saprotrophic RAF communities by soil strata. Pathogenic fungi were rare (relative abundance of 1.3% to 4.0%), most likely occurring by chance, as suggested previously [18]. Overall, the consistent separation of the fungal groups according to soil strata indicate that soil quality is a dominant factor across a large biogeographic scale.
Saprotrophic fungi usually prefer the organic layer, where they degrade organic material to obtain C and contribute to the decomposition of organic substances in forest soil [72,73]. Mycorrhizal fungi mine the mineral soil for inorganic compounds and rely on host-derived C, and therefore may reside in older litter layers and mineral topsoil [74]. Our results show that soil C, N, C/N ratio, P sol and cations differed strongly between the organic layer and mineral top soil, similar to other studies [19,27,28]. Thus, the decline in soil nutrients from the organic layer to the mineral topsoil would support the well-known shift from saprotrophic to symbiotrophic fungi [19,22,75], if the RAF assemblage was only driven by soil properties. However, in support of our initial hypothesis, the relative abundance of symbiotrophic fungi was higher on roots in both soil strata, most likely because mycorrhizal fungi have direct access to carbohydrates from their host, irrespective of the soil layer and therefore have an advantage over saprotrophic fungi colonizing the rhizoplane. Still, strong influence of soil layers on the composition of the fungal assemblages was evident, as the richness of saprotrophic fungi was higher under nutrientrich conditions (organic layer), while the richness of symbiotrophic fungi was higher in nutrient-poor conditions (mineral soil). These results are in line with the results on soil fungi in boreal forests [22,75] and arbuscular mycorrhizal fungi in temperate forests [19]. Whether the shifts in the richness are the result of different nutritional acquisition strategies of saprotrophic and symbiotrophic fungi or of antagonistic relationships requires further studies. Future studies should also address the impact of root necromass on RAF diversity.
Since we conducted our study in different biogeographic regions, we were able to detect regional effects on the vertical distribution of RAF. This is a novel aspect shedding light on fungal assemblies. For example, RAF richness and β-diversity were lower in the region with drier, more acidic and nutrient-poor forest soil than in regions with cooler, moist climate and nutrient-rich soil. This can be explained by the classical ecological theory that environmental stress conditions result in reduction of species diversity, selecting species that can tolerate harsh condition [76]. Mycorrhizal fungi are known to be sensitive to elevated temperature and drought [77][78][79][80][81]. Therefore, the climatic differences might have negatively affected the richness of the symbiotrophic fungi in the SCH region compared to ALB or HAI, where the soils have a higher water-holding capacity [36]. In addition, more fertile soil conditions, as present in ALB and HAI, may have favoured higher mycorrhizal diversity, in line with previous studies [12,18,20,82].
At the plot level, we observed a low nestedness and high turnover of RAF between the soil strata but the turnover of symbiotrophic fungi was lower than that of saprotrophic fungi. This result was expected (hypothesis ii) since mycorrhizal fungi colonizing roots are directly influenced by roots of their host. Host effects are known from studies showing strong impact of vegetation [12,14,37], tree identity [70,[83][84][85] and root chemistry [20,86]. Host effects may stabilize mycorrhizal patterns at larger geographic scales [17]. Here, we found lower turnover of symbiotrophic than of saprotrophic fungi within the RAF assemblies. Therefore, the present results indicate that root traits may drive symbiotrophic differently from saprotrophic fungi by stabilizing the communities. However, it should be noted that the differences between the turnover of saprotrophic and symbiotrophic fungi were relatively small compared to the overall turnover, underpinning strong abiotic habitat impact on the RAF.

Indicator Taxa and Phylogenetic Groups Uncover Different Strategies of Root-Associated Fungi
An important goal of this study was to test the hypothesis that RAF patterns indicate response traits either to soil strata or to regional habitat conditions. We chose a novel approach by grouping the fungal taxa according to their phylogenetic rank at the order level and measuring the effect size imposed by soil strata on in the fungal abundance in different regions and across all regions. Phylogenetic community structures are known to carry information on ecological assemblages because the relatedness of members in a community suggests similar ecological requirements and functions of phylogenetically related taxa [87,88]. For example, Kohler et al. [89] and Nagy et al. [90] showed phylogenetic relatedness of fungal orders based on gene counts for certain saprotrophic traits such as degradation of cellulose or lignin. Using gene counts, a large fraction of variance of potential fungal traits for N and P transformation was explained at the level of the subphylum or phylum [91]. Mycorrhizal species (identified on roots by morphotyping-Sanger sequencing) showed stronger phylogenetic clustering in drier and acidic habitats than under cooler and humid soil conditions [37]. Here, we discovered that distinct RAF orders showed a territorial behaviour as they were always predominant in a certain soil stratum, irrespective of differences in the environmental conditions among the region, whereas other fungal orders showed flexible behaviour with changing effect sizes depending on the regional habitat conditions.
Our results indicate that the behaviour whether a fungal group was territorial or flexible distinguished saprotrophic and symbiotrophic orders in RAF assemblages. Shared territorial mycorrhizal fungal orders occurred in the organic layer (Boletales) and in mineral soil (Russulales and Cantharellales), whereas we did not find any saprotrophic fungal order nor any saprotrophic indicator species that was shared in the mineral soil among the three regions. As expected the saprotrophic groups were either flexible or enriched in the organic layer, the latter group including Polyporales, Pleosporales, and Sordariales. Polyporales (some causing brown-rot of timber) are well known for their efficient lignolytic capabilities to degrade deadwood [92]. They were also enriched in the litter and organic soil surface in other temperate forests [21,69,93]. At the species level, Calycellina fagina was a shared indicator taxon in the organic layer across all study regions. Members of this genus are usually growing on dead-wood and plant matter in forest floors and play a role in decomposition processes [94,95].
Thelephorales, Sebacinales (both ectomycorrhizal fungi), Agaricales and Helotiales (mainly mixed symbiotrophic/saprotrophic fungi with small contributions (<4%) of pathogenic fungi) were major fungal orders in RAF assemblages with flexible effect sizes. This result indicates that dominant members of these orders are driven by the specific habitat conditions in each region. On the contrary, Russulales were consistently enriched in the mineral soil, suggesting territorial behaviour of the members of this order. Russulales are highly abundant ectomycorrhizal species in temperate beech forests [37,84,96]. All known members of the Russulales exhibit a hydrophilic contact exploration type [97], absorbing nutrients from the surrounding soil [98]. At the taxon level, Russula sp. was a shared indicator species on roots in mineral soil across all regions. Whether the availabilities of mineral nutrients foster root colonization with Russulales is unclear because fertilization experiments with P or N showed divergent effects on the abundance of Russulales (P fertilization: decrease [21,99]), and no effect of N fertilization in boreal forests [100,101]). Russulales have been classified as nitrophilic fungi but they are also capable of producing extracellular enzymes that could degrade organic matter in litter and soil [97]. Recent genome analyses indicated that they lost most enzymes required for cell wall degradation but maintained Mn peroxidases and chitinases [102]. Their involvement in decomposition is still questionable since they are relatively inactive in the period of organic matter decomposition [103]. Similar to Russulales, we found positive effect sizes for Cantharellales in the mineral layer. Members of this order also contain class II peroxidases capable of degrading complex organic compounds [89,104]. Perhaps these enzymes are beneficial for degradation of recalcitrant organic compounds still present in the mineral soil.
Boletales were the only mycorrhizal fungal order that showed significantly higher abundance in the RAF of the organic layer than in mineral soil across all regions. The relative abundance of Boletales increased after P fertilization [21]. In our studied forest soils, P sol and N were higher in the organic layer than in mineral soil. Many Boletales species are characterised by long-distance rhizomorphs and can explore soil far beyond the rhizosphere [97]. Therefore, Boletales are considered as beneficial in nutrient-limited forest ecosystems. They can access distant resources [105] and, thus, meet the P demand of trees [106]. Almeida et al. [106] detected an increase in hyphal biomass of Imleria badia (Boletales) that accessed apatite (a recalcitrant P source) in N-fertilized soil but not in N-fertilized soil amended with soluble P sources, underpinning a role of these taxa for P nutrition.
At the species level, we detected Laccaria amethystina (Agaricales) as an indicator taxon on roots in the organic layer across all three regions, emphasizing the wide spread of this ectomycorrhizal fungal species across temperate forest ecosystems. Laccaria amethystina is among the most common species that colonize roots of oak and beech tree species in temperate deciduous and deciduous-coniferous forest ecosystem in Europe [107]. Besides Laccaria amethystina, related species such as Laccaria maritima were detected on beech roots in our studied forests [37,84]. Species from the Laccaria genus (e.g., Laccaria bicolor) have little capacity to degrade organic matter [108], which renders the occurrence of Laccaria amethystina as an indicator species in the RAF community in organic soil surprising. However, Laccaria amethystina is known as an ammonia utilizing fungus [109]. Therefore, we speculate that Laccaria sp. in the organic layer may benefit from N released by the degradation of organic material.
The consistent enrichment patterns of distinct RAF orders and species in response to soil layers across the regions support that ecological niche partitioning strongly influenced the differentiation of root-associated fungal community structures. Our results imply that RAF assembly entails two strategies encompassing flexible and territorial habitat colonization by different fungal taxa. However, it is important to take into account that this aggregated behaviour does not reflect collective response patterns of the species within a given order. This caveat is reflected by the identification of 39 saprotrophic and 63 symbiotrophic indicator taxa among which only three were confined to a distinct soil stratum and were shared among the regions. Other indicator taxa occurred only in one distinct region or were shared among the regions with higher similarities in soil resources [36]. Taking together individual variation of taxa in the RAF assemblage and relatively stable regional responses at a higher phylogenetic level, we suggest that this behaviour could lead to emergent properties of ecosystems that are more than the sum of their individuals.

Conclusions
Our results support the ecological concept that resource partitioning and phylogenetically conserved properties determine the ecological communities. A distinct response of RAF orders and indicator taxa that were specific to soil layer and region indicated that habitat conditions strongly influence the differentiation of the RAF community structure. The results also support a phylogenetic signature for niche partitioning since several fungal orders were enriched in distinct soil layers across a large-scale biogeographic gradient, irrespective of the habitat conditions. Saprotrophic fungal orders and indicator taxa were preferentially enriched in the organic layer and mycorrhizal orders and species in the mineral soil. However, our knowledge is still limited on fungal traits and fungal nutritional preferences. To better understand the role of RAF in shaping ecosystem properties, further studies on the interaction of mycorrhizal and saprotrophic fungi are required.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/microorganisms9102131/s1, Figure S1: Relative abundance of root-associated fungi compiled at the order level in the organic layer and mineral soil and in different regions. Table S1: Location of the 150 study plots and fine root biomass. Table S2: Observed species richness (OTU-based), estimated richness (Chao1), Shannon diversity (H ) and evenness (E H ) of root-associated fungi in different soil layers across three biogeographic regions. Table S3: Beta-diversity of root-associated fungal assemblages of all (OTUs), symbiotrophic (SYM), and saprotrophic (SAP) fungi between organic layer and mineral soil. Table S4: Relative abundance, number of sequences and number of taxa (OTU-based) of the root-associated fungal orders present in the organic layer and mineral soil in the Biodiversity Exploratory forest plots. Table S5: Classification of the top nineteen rootassociated fungal orders according to functional groups. Table S6: Root-associated symbiotrophic and saprotrophic indicator taxa in the organic layer and mineral soil in three biogeographic regions. Funding: The research was funded by the DFG Priority Program 1374 "Infrastructure-Biodiversity-Exploratories" (Po362/18-4 and Po362/ . We acknowledge the support by the Open Access Publication Funds of the Göttingen University. Data Availability Statement: All data are available in the BExIS database (https://www.bexis.unijena.de) under the following accession numbers (data owner): Soil chemistry-26228 and 26229 (Polle), root-associated fungal taxonomic and functional assignment-30973 and 30974 (Polle) and total fine root dry mass in organic layer and mineral topsoil-31048 and 31049 (Polle).