Large-diameter trees and deadwood correspond with belowground ectomycorrhizal fungal richness

Large-diameter trees have an outsized influence on aboveground forest dynamics, composition, and structure. Although their influence on aboveground processes is well studied, their role in shaping belowground fungal communities is largely unknown. We sought to test if (i) fungal community spatial structure matched aboveground forest structure; (ii) fungal functional guilds exhibited differential associations to aboveground trees, snags, and deadwood; and (iii) that large-diameter trees and snags have a larger influence on fungal community richness than smaller-diameter trees. We used MiSeq sequencing of fungal communities collected from soils in a spatially intensive survey in a portion of Cedar Breaks National Monument, Utah, USA. We used random forest models to explore the spatial structure of fungal communities as they relate to explicitly mapped trees and deadwood distributed across 1.15 ha of a 15.32-ha mapped subalpine forest. We found 6,177 fungal amplicon sequence variants across 117 sequenced samples. Tree diameter, deadwood presence, and tree species identity explained more than twice as much variation (38.7% vs. 10.4%) for ectomycorrhizal composition and diversity than for the total or saprotrophic fungal communities. Species identity and distance to the nearest large-diameter tree (≥ 40.2 cm) were better predictors of fungal richness than were the identity and distance to the nearest tree. Soil nutrients, topography, and tree species differentially influenced the composition and diversity of each fungal guild. Locally rare tree species had an outsized influence on fungal community richness. These results highlight that fungal guilds are differentially associated with the location, size, and species of aboveground trees. Large-diameter trees are implicated as drivers of belowground fungal diversity, particularly for ectomycorrhizal fungi.


Background
Trees, snags, and deadwood interact directly and indirectly with fungi by providing resources aboveground and belowground and modifying the microenvironment (Smith and Read 2010). These interactions vary with tree species, size, age, and distance, which have all been implicated as influencing fungal community composition (Teste and Simard 2008;Birch et al. 2021). Tree species vary in the amount and nutrient content of their litterfall (Vogt et al. 1995;Keane 2008), root exudates (Vives-Peris et al. 2020), longevity, and life-history strategies, all of which can alter the environment and influence the assembly of fungal communities (Birch et al. 2022b). In addition to species identity, individual trees' size, location, and age can also influence fungal communities (Nguyen et al. 2016b;Rudawska et al. 2018;Wasyliw and Karst 2020). Trees can also directly host symbiotic and saprotrophic fungi on and within their roots, and thus directly influence fungal communities (Smith and Read 2010). Because these tree-provided resources vary in quality and quantity, the aboveground structure and composition of forests likely influence the belowground community of fungi.
Large-diameter trees contain a disproportionate amount of aboveground biomass within forests ) and deadwood (Lutz et al. 2020, b). The largest trees may have a distinctive influence on belowground fungal community structure through greater root length and area (Bauhus 2009), higher respiration relative to photosynthesis (Tang et al. 2014), or through their longevity (Birch et al. 2021). The influence of a tree on the environment may decline with increasing distance, with the strength of the distance-decay relating to the size of the tree and the density of the surrounding forest. Thus, we expect that large-diameter trees might have a large influence on belowground fungal communities, even at great distances. Large-diameter trees could support unique fungal communities owing to specific fungi promoting greater growth and survival or because large-diameter trees select different fungi due to altered nutrient requirements and larger rooting extents that change with tree ontogenesis (Fig. 1). The greater resource requirements, rooting extents, and potential for differential bi-directional selection between fungi and large trees could conceivably select for different ectomycorrhizal fungal communities relative to small trees. Likewise, large-diameter trees are likely to release larger amounts of litterfall and necromass (Maguire 1994) relative to smaller trees which could support saprotrophic communities. In either case, we would expect that largediameter tree location, size, and species may have consequences for belowground fungal communities.
The spatial extent over which fungal community composition is similar, i.e., autocorrelated, is important when parsing the abiotic drivers of plant composition and structure (Kim and Shin 2016), as the inclusion of spatial autocorrelation can improve identification of the environmental variables influencing communities. Likewise, understanding the drivers of forest and fungal community composition necessitates identification of the extent over which they vary in space. Many soil fungal communities are autocorrelated between 2-3 m, though some species can have spatial autocorrelation upward to 17 m, indicating large genets (Lilleskov et al. 2004;Boraks et al. 2021). Trees, deadwood, and abiotic variables, such as pH, vary with space (Lutz et al. 2013(Lutz et al. , 2021a, and if fungal communities are influenced by these features, we would expect congruent spatial autocorrelation between aboveground forest structure and belowground fungi. For example, the diversity of tree litter corresponds with Fig. 1 Conceptual diagram of how fungi may promote (path A) or be selected by large-diameter trees (path B) with time. Under path A the same fungal community endures from tree recruitment to old-growth and promotes the tree survival and growth. Under path B the tree selects different fungal partners as it progresses from recruitment to old-growth saprotrophic diversity and may alter saprotrophic community composition due to saprotroph specialization in species specific carbon and litter (Hanson et al. 2008;Otsing et al. 2018).
The subalpine forests of the southwestern United States represent a diverse and spatially isolated collection of mixed-conifer forests. These forests frequently encompass multiple coexisting overstory Pinaceae species as well as aspen and tall shrubs. Decomposition rates are slow, so deadwood remains for centuries (Kueppers et al. 2004). Through a spatially intensive survey of soil fungal communities, we sought to identify fungal communities in a high-elevation, old-growth forest. Specifically, we sought to test if (H1) fungal communities had spatial autocorrelations that closely matched aboveground forest structure, (H2) fungal functional guilds would exhibit differential associations to aboveground tree, snag, and downed wood structure, and (H3) that large-diameter tree and snags would have a larger influence on fungal community diversity than would smaller-diameter trees.

Site description
The study was located within the subalpine, mixedconifer forest at the Utah Forest Dynamics Plot (UFDP) located within Cedar Breaks National Monument, Utah, USA (37.66°N, 112.85°W) (Fig. 2). The UFDP was established according to the protocols of the Smithsonian Forest Global Earth Observatory network (Davies et al. 2021) with all living and dead stems ≥ 1.0 cm diameter at breast height (DBH; 1.37 m) in 15.32 ha being identified, mapped and annually assessed for mortality (Lutz 2015). This survey included 25,478 living trees, 4400 snags, and 5,310 pieces of deadwood originating from the main stems of trees ≥ 10 cm diameter (Lutz et al. 2021a). The UFDP contains 11 overstory tree species with maximum species ages ranging from 220-1400 years old (Furniss et al. 2017;Birch et al. 2021). Deadwood was divided into five decay classes (Harmon et al. 2008). Species-specific tree-ring chronologies have been collected for all tree species within the UFDP with more than 15 individuals (Anderson-Teixeira et al. 2022).
We established a sampling grid within a 1.15-ha section of the UFDP characterized by a diverse, old-growth closed-canopy forest. We selected this area due to the diverse composition and structure of the aboveground forest and similarity in aspect across a contiguous area. This area spans an elevation of 3037-3089 m with a mean annual temperature of 3.9 ºC and mean annual precipitation of 768 ± 119 mm, of which 506 ± 137 mm occurs as snow. Soils are categorized within the Mollisol (Argicryolls) and Inceptisol (Haplocryepts) orders and are alkaline (pH 7.0-8.1), calcareous, and have high base cation content. Soil textures within the sampling area are generally categorized as clay to clay loams.
The sampling area is composed of Abies bifolia A. Murray bis (Rocky Mountain subalpine fir), Populus tremuloides Michaux (quaking aspen), Picea pungens Engelmann (Colorado blue spruce), Picea engelmannii var. engelmannii Parry ex Engelmann (Engelmann spruce), Pinus longaeva D. K. Bailey (Great Basin bristlecone pine), and Pinus flexilis E. James (limber pine), with low abundances of Pseudotsuga menziesii var. glauca (Mayr) Franco (interior Douglas-fir), Pinus ponderosa Douglas ex Lawson & C. Lawson (ponderosa pine), Abies concolor (Gordon & Glendinning) Hildebrand (white fir), Juniperus communis Linnaeus (common juniper), and Ribes cereum Douglas (wax currant) in the adjoining forest (Additional file 1: Table S1) (Furniss et al. 2017). The forest has had no recorded management actions and has multiple cohorts originating from natural regeneration following patchy and infrequent high-severity fire with Fig. 2 The location of soil sampling (a) within the Utah Forest Dynamics Plot (UFDP), Utah, USA. Samples were taken 3 m away from monumented pins (b) placed within a 20 × 20 m grid at the UFDP the most recent fire likely occurring in 1802 (Lutz et al. 2021a). Live tree ages in the dataset range from 2-year seedlings to 1400 years old (Birch et al. 2021). Characteristic disturbances include infrequent, high-severity fire, drought, bark-beetle outbreaks, high wind events, and growing season frost (Lutz et al. 2021a;Birch et al. 2022a). The 1.15-ha sampling area is structurally diverse (Fig. 3) with both small (< 5 m) and large (~ 10 m) canopy gaps distributed across the multiple cohorts of mixed conifers.

Belowground sampling
To sample fungi, we collected soil samples across a 115 × 100 m grid overlaid across the monumented 20 × 20 m grid present in the UFDP (Fig. 2a). We sampled at distances of 3.0 m from these grid-points at fixed azimuths of 45º, 135º, 225º, and 315º ( Fig. 2b). We used a soil knife (Zenport; Sherwood, Oregon) to sample soil in a 4 × 4 cm area to a depth of 25 cm. We sanitized the knife using a 10% bleach solution between each sample and froze samples at -20 ºC before processing. At each sampling point, we also measured, to the nearest centimeter, the depth of the litter, fermentation, and humus layers (LFH layer). The LFH layer comprises the organic (O) horizon that overlies the uppermost mineral soil (A horizon).
To assess the concentration of mineral soil nutrients (phosphorus, ammonium, nitrate, calcium, iron, potassium, magnesium, manganese, effective cation exchange capacity, base saturation, and sodium) and pH, we sampled soils to 10 cm in a semi-random grid across the entire UFDP (Additional file 2). Soil nutrient concentrations were measured following the analytical methods of Baldeck et al. (2013). Non-nitrogen nutrients were extracted using a Mehlich III solution and analyzed using an atomic emission-inductively coupled plasma (AE-ICP, Perkin Elmer Inc., MA, USA). Nitrogen was extracted as NH + 4 and NO − 3 with a 2 M KCI solution and analyzed using an auto-analyzer (OI FS 3000, OI Analytical, College Station, TX, USA).

Sample preparation
To prepare soil samples for DNA extraction, we thawed each sample and sieved them through progressively finer sieves up to 2 mm to remove fibrous organic matter. We bleached and dried sieves between each sample to avoid cross-contamination. We freeze-dried samples until they had reached a stable weight. We pulverized freeze-dried samples to a fine powder using a Tissuelyser 2 (Qiagen, Hilden, Germany) and used MoBio DNeasy PowerSoil Kits (Qiagen, Hilden, Germany) to extract DNA following the manufacturer's instructions. After extraction of DNA, we iteratively amplified and visualized the DNA using gel-electrophoresis (Additional file 2).

Sequencing of fungal rDNA
We used MiSeq sequencing targeting the ITS2 region of fungal rDNA using fITS7 (Ihrmark et al. 2012) and ITS4 (White et al. 1990) primers. We amplified extracted rDNA through a two-step polymerase chain reaction (PCR). We ran the first PCR using 25 μL sample volumes with 2.5 μL DNA template, 7.8 μL water, 12.5 μL Platinum SuperFi PCR MasterMix (Thermo Fisher Scientific, Waltham, USA), 0.125 μL Bovine Serum Albumin, and 1 μL of each of the forward and reverse primers (95 ºC for 3 min; 95 ºC for 30 s, 50 ºC for 1 min, 72 ºC for 30 s, 35 cycles; and 72 ºC for 10 min). The samples were pooled and cleaned using Sera-Mag Select beads (GE Healthcare, Chicago, USA) and washed in two 80% ethanol washes. We used Nextera Index primers with N7 (N5) for the forward (reverse) reads for the second PCR. For the second PCR, a 25 μL volume was used with 5 μL of water, 12.5 μL of Platinum SuperFi PCR MasterMix, 2.5 μL N7xx primers, 2.5 μL N5xx primers, and 2.5 μL DNA template (95 ºC for 3 min; 95 ºC for 30 s, 55 ºC 30 s, 72 ºC for 30 s, 72 ºC for 5 min, for 8 cycles). Sequencing was conducted at the University of Alberta Molecular Biological Sciences Unit, with a MiSeq sequencing platform.

Bioinformatics
We denoised, quality-checked, and quality filtered the demultiplexed sequences in the R package dada2 1.18.0 (Callahan et al. 2016) to create amplicon sequence variants (ASV). We filtered sequences to remove those with ambiguous bases and limited the maximum number of expected errors to two per sequence. Following identification and removal of chimeric sequences (Additional file 1: Table S2), we assigned taxonomy from the UNITE database (Nilsson et al. 2019) with the IDTAXA function in the DECIPHER 2.18.1 R package (Wright 2016).
After assigning taxonomy, we analyzed the results using the phyloseq 1.38.0, vegan 2.5-7, ggplot2 3.3.5, and ggpubr 0.4.0 R packages (McMurdie and Holmes 2013; Wickham 2016; Oksanen et al. 2020;Kassambara 2020). We had several duplicated samples that were amplified as safeguards against the loss of low-quality samples. In four cases the duplicates and originals were successfully amplified. We merged the duplicate samples by averaging sequence abundance. After merging, we visualized the distribution of sequence abundance and rarefied to 12,500 sequences to account for uneven sampling depth (Additional file 3: Fig. S1). We removed four samples that were below the 12,500-sequence cutoff.
To assign functional guilds to ASV, we used the FUN-Guild database (Nguyen et al. 2016a), which assigns functional guilds to fungi based on the taxonomic assignment of the sequence. Common guilds include "Ectomycorrhizal", "Saprotrophic", or some combination of guilds as in "Ectomycorrhizal-Saprotrophic" when fungi are capable of existing across multiple functional guilds. Because of unknown or unresolved taxonomic assignments, many ASVs are classed as having an 'Unknown' guild. Additionally, many sequences are assigned to multiple potential guilds and were classified as 'Mixed EM' if one of the potential guilds was ectomycorrhizal. We included mixed guilds within the Saprotrophic classification if they included at least one known saprotrophic function, i.e., "wood-saprotroph-ectomycorrhizal".

Analyses
We used kriging to interpolate measurements of pH, LFH depth, and soil nutrients, across our sampling area. We used the 'powerTransform' function from the car 3.0-7 R package (Fox and Weisberg 2018) to normalize the variables. After transformation, we excluded nitrate and base saturation, which were still strongly non-normal. We generated semivariograms for each normalized variable and assessed the spatial autocorrelation of each. We used the 'krige' function to interpolate the remaining measurements and then transformed each nutrient to their original units (mg kg −1 ) for analysis (Additional file 3: Fig. S2). For the LFH, we included an additional 160 depth measurements of LFH (n total = 295) that were taken across the sampling area (Lutz et al. 2021a).
To assess spatial autocorrelation of fungal communities and functional guilds we used the 'mantel.correlog' function to plot mantel correlograms with 999 permutations (McMurdie and Holmes 2013). Distances were constructed using Euclidean distances between sampling locations. After assessing the entire fungal community, we split apart the fungal dataset into ectomycorrhizal and saprotrophic guilds and ran separate correlograms of each guild's community.
To compare the spatial autocorrelation of belowground fungi with aboveground trees and deadwood, we assessed the spatial autocorrelation of all living trees and mapped deadwood. Deadwood biomass calculations were conducted using the R packages rgdal version 1.4-8, rgeos version 0.5-2, sp version 1.3-2, and sf version 0.8-1 (Pebesma and Bivand 2005;Pebesma 2018;Bivand et al. 2019a, b). To assess the spatial autocorrelation of living trees, we binned trees within 5 × 5 m cells and summed the basal area, by species, within each cell. After binning trees, we calculated the basal area of each tree species and standardized the community dataset using a Hellinger transformation. After standardization, we detrended the tree data using a linear regression against the centroid coordinates of each cell. Finally, we calculated the Euclidean distance of the regression residuals to calculate the mantel correlogram for the tree community.
To assess the spatial autocorrelation of deadwood, we calculated the biomass of all mapped deadwood (Additional file 2). We aggregated piece-level values of deadwood biomass within 5 × 5 m cells divided across the 1.15-ha sampling area. For pieces that lay across cell boundaries, we virtually separated them into new fragments, each wholly contained in a single cell. After summation, we standardized the dataset with a Hellinger transformation, detrended the data using a linear regression, and calculated the Euclidean distance of the regression residuals to calculate the mantel correlogram for deadwood.
We used the "specaccum" function in the vegan 2.5-7 R package to generate species accumulation curves by fungal guilds (Oksanen et al. 2020). We calculated Shannon's diversity using the 'estimate_richness' function in the R package phyloseq 1.38.0, on un-rarefied community data (McMurdie and Holmes 2013).
To assess the influence of soil nutrients on fungal community richness, we aggregated nutrients using a principal coordinate analysis with a correlation matrix to collapse nutrients into a primary and secondary component which collectively explained 82.7% of the variation within nutrients (Additional file 1: Table S3, Additional file 3: Fig. S3).
To assess the influence of aboveground forest structure, nutrients, and topography on the total and guildspecific communities, we used a permutational analysis of variance (PERMANOVA) with 999 permutations. We derived elevation from 0.5 × 0.5 m LiDAR (UGRC 2018), with landscape profile curvature calculated using ArcGIS 10.7 (ESRI, Redlands, CA, USA). To characterize the tree and snag community around each sample, we selected the maximum extent of their spatial autocorrelation of each fungal community (e.g., total, ectomycorrhizal, saprotroph) as a buffer in which to assess aboveground tree, snag, and downed wood characteristics. We calculated the basal area of all trees and snags and weighted each by the distance from the sampling location before summing the basal area, by tree and snag species. For deadwood, we summed the biomass of each species of deadwood within the specific extent of the sample. Finally, we retrieved the point-estimates of the nutrient principal components, LFH depth, landscape curvature, and elevation for each sample.
We used a random forest model to identify key aboveground factors influencing belowground fungal richness. We chose a random forest model as it accommodates non-normal data and can assess non-linear relationships between the dependent and independent variables. We used the randomForest 4.6-1.4 R package (Liaw and Wiener 2002) and ran separate random forest models for total, ectomycorrhizal, and saprotrophic community richness. With richness as our dependent variable, we tested both the identity and distance of the nearest tree and nearest large-diameter tree (> 40.2 cm DBH), by tree species, as our independent variables. The 40.2 cm DBH threshold (Additional file 3: Fig. S4) was selected as an a priori threshold as it is the diameter in the UFDP above which trees contain half of the aboveground biomass ) and this metric can thus be used to set similar thresholds for large-diameter trees in other forests. In all models, we included the first two nutrient components derived from a principal component analysis (PCA), as well as elevation, LFH depth, and landscape curvature as abiotic covariates. We used the 'randomForest' function with 999 trees (Liaw and Wiener 2002) and default settings for the number of splits. We plotted partial dependence plots for each variable and graphed the increase in mean square error (MSE) for each variable.
Statistical analyses were conducted within the R Studio 2021.09.1 interface for R 4.2.1 (R Core Team 2020; R Studio Team 2021).

Fungal composition, richness, and function
After dada2 quality control, we had 2,305,722 sequences from 126 samples with a mean sequence length of 309 ± 3.5 bases (SE). After rarefaction and merging of duplicate samples, we had 117 samples and found 6,177 fungal ASVs with a mean of 155 ± 5 ASVs per sample (Additional file 3: Fig. S5). Mean Shannon's diversity was 4.0 ± 0.0 with a range of 2.7-4.9 per sample. A total of 2207 ASVs (35.7%) representing 13.4% of sequences were not assigned a taxonomic level of Phylum or below. Species accumulation curves indicate that we sampled > 90% of the fungal community and that 99 samples were needed to sample 90% of all ASVs (Additional file 3: Fig. S6). Ascomycota represented 52.1% of all sequences and had 2003 ASVs. In contrast, Basidiomycota represented 31.6% of all sequences and 1,446 ASVs. The most read-abundant ASVs were within the ectomycorrhizal genera of Inocybe, Wilcoxina, Tricholoma, Geopora, and Sebacina (Fig. 4). The largest single functional guild was ectomycorrhizal (597 ASVs; 25.5% of sequences) while 4349 ASVs (70.4%) representing 54.8% of sequences were not assigned a functional guild (Additional file 1: Table S4). Fungi capable of acting as saprotrophs accounted for 956 ASVs and 12.9% of total sequences. The most read-abundant ASVs representing saprotrophs were unidentified fungi within the Thelephoraceae and Pseudeurotiaceae families.

Spatial autocorrelation of fungi, trees, and deadwood
Within the sampling grid, composition of the entire fungal community showed weak spatial autocorrelation at distances from 4.5 m to 5.8 m (Fig. 5a). In contrast, ectomycorrhizal fungal communities showed consistent autocorrelation from 0.6 to 16.1 m (Fig. 5b). Similarly, saprotrophic fungal communities had spatial autocorrelation from 0.6 to 19.9 m (Fig. 5c). Both forest tree composition and deadwood biomass were positively autocorrelated from 0.6 to 30.5 m (Additional file 1: Table S5, S6, S7). Negative spatial autocorrelation existed for both trees (> 30.0 m) and deadwood (> 30.5 m) (Additional file 1: Table S8,

Influence of aboveground forest structure on fungal composition and richness
After identifying the maximum distance of spatial autocorrelation of each fungal guild, we used this value to test the distance-weighted influence of all tree species within each guilds' neighborhood. Total fungal community dissimilarity was associated with the distance-weighted basal area of living P. tremuloides and R. cereum within 6.0 m (Table 1). Ectomycorrhizal fungal community dissimilarity was associated with soil nutrients and living distance-weighted basal area of A. bifolia and A. concolor within 16.0 m. Finally, saprotrophic fungal community dissimilarity was associated with soil nutrients, depth of the LFH layer, elevation, living distance-weighted basal area of P. longaeva, P. menziesii, and R. cereum along with dead distance-weighted basal area of A. bifolia and P. tremuloides within 20.0 m (Table 1). Random forest models testing the influence of the nearest tree and nearest large-diameter tree varied substantially, by fungal guild, in their variation explained. Total fungal richness was better explained by models with distance to the Fig. 4 The total proportional read abundance for the 30 most read-abundant fungal genera within the Utah Forest Dynamics Plot, Utah, USA nearest large-diameter tree (10.4% variation explained) than by models with distance to the nearest tree (9.5%). The variables with the greatest increase in mean square error, a measure of variable importance, were distance to the nearest large-diameter dead A. bifolia, large-diameter living P. flexilis, and large-diameter living P. engelmannii (Additional file 1: Table S10). Partial dependency plots indicate that fungal richness was higher by nearby (0-20 m) large-diameter dead A. bifolia, depressed near large-diameter live P. flexilis (0-20 m), and large-diameter living P. engelmannii (0-20 m) (Additional file 3: Fig.  S9).
Ectomycorrhizal fungal richness was better explained by distance to the nearest large-diameter tree (38.7% of variation explained; Fig. 6) than by distance to the nearest tree (34.7%). Fungal richness increased at close distances (0-25 m) to large-diameter living P. engelmannii, increased with increasing values of soil component 2, and was depressed at close distances (0-25 m) around largediameter dead P. engelmannii (Fig. 6b). In contrast, saprotrophic fungal richness was poorly explained by distance to the nearest large-diameter tree (− 11.0%) and distance to the nearest tree (− 15.8%). The most important variables for saprotrophic fungal richness were distance to the nearest large-diameter live A. bifolia, elevation, and large-diameter dead P. flexilis (Additional file 3: Fig.  S10).

Discussion
In support of our first and second hypotheses, we found that tree and belowground fungal communities were associated with one another and that specific aboveground variables were associated with ectomycorrhizal and saprotrophic fungal communities. Our third hypothesis was supported by the total and ectomycorrhizal fungal community richness, which had more variation explained by the nearest large-diameter tree than by the nearest tree. Cumulatively, our results highlighted a remarkably rich and unexplored fungal community that was associated with distinct abiotic and biotic variables operating at differing scales within the forest. Further research is needed to fully parse the mechanisms through which these variables influence fungal community composition and richness. However, our methods highlight the utility of pairing rigorous belowground sampling with spatially explicit tree and deadwood data.

Fungal guilds and forest neighbors: unique spatial patterns and associations
Composition of fungal communities responded to both forest composition, soil nutrients, and topography, though at different spatial scales. Both ectomycorrhizal and saprotrophic fungal communities were autocorrelated beyond 15 m, nearly threefold greater distance than the total community. The autocorrelation of the total community is similar, though greater, than the 2-3 m distances commonly associated with fungal species (Lilleskov et al. 2004;Boraks et al. 2021). In contrast, the spatial autocorrelation of ectomycorrhizal and saprotrophic fungal guilds was far greater, though many meters short of the spatial autocorrelation of the forest (30.0 m) and deadwood (30.5 m). One potential reason for the offset between spatial scales of the aboveground forest and belowground fungal guilds is that individual trees most strongly influence their immediate 2. The living and dead basal area of specific tree species was associated with both the composition and richness of each fungal guild. However, ectomycorrhizal fungal composition and richness were better explained by our data than were the entire or saprotrophic communities. Our results mirror other studies which have found that tree species more strongly influenced ectomycorrhizal fungi than saprotrophic fungi (Nguyen et al. 2016b;Chen et al. 2019). Interestingly, the current survey detected more than three times as many fungal ASVs compared to our previous, targeted root-and soil samples collected around P. flexilis, P. longaeva, and P. menziesii at the UFDP (Birch et al. 2021(Birch et al. , 2022b. Because ectomycorrhizal fungal richness increases with aboveground diversity, particularly at the genus-level (Gao et al. 2013), the diverse forest within the 1.15 ha likely drives a portion of the high fungal richness.
Composition of ectomycorrhizal fungal communities were better explained by all models than was the total community. This is likely due to divergent responses to environmental variables across the various functional guilds that are obscured by the poor taxonomic resolution of the entire community. Of note, P. longaeva was not strongly associated with ectomycorrhizal fungal composition or richness, despite its dominance across our sampling area. Only the basal area of living A. bifolia and A. concolor were significant in explaining ectomycorrhizal fungal community composition. Because P. longaeva is widely dispersed in the 1.15 ha, the locally rare tree species, such as A. bifolia, or A. concolor, may be hosting fungal communities that are not well-suited for P. longaeva. Supporting this idea, the most important variables explaining ectomycorrhizal fungal richness were the distance to the nearest P. engelmannii, and soil nutrients. While P. engelmannii was once abundant in the forest, after a Curculionidae Dendroctonus rufipennis Kirby (spruce beetle) outbreak in the 1990s, they are present mainly at small diameters across the landscape and uncommon within the UFDP (DeRose and Long 2007). The current scarcity of P. engelmannii and its association with fungal richness, suggests that locally rare tree species are important for promoting belowground fungal diversity.

Decomposing the relationships of saprotrophic fungi and the aboveground forest
In all cases, the saprotrophic fungal community had less variation explained by the model variables than did the ectomycorrhizal fungal community. Notably, the saprotrophic fungal composition was associated, in part, with living P. longaeva, the dominant tree species within the forest. Saprotrophic fungal composition was also influenced by the basal area of dead A. bifolia and P. tremuloides. Litter and wood from A. bifolia have a longer residence time than Picea wood (Bigler and Veblen 2011), while P. tremuloides likely decays faster than the other species (Angers et al. 2012). Thus, the differing residence times and parent material of the litter and deadwood may alter the surrounding saprotrophic communities. The turnover time for deadwood carbon has been estimated at being more than 580 years in similar forests (Kueppers et al. 2004), highlighting the enduring legacy of deadwood on saprotrophic communities and forest nutrients. Live and dead tree species were better predictors of saprotrophic richness than were soil nutrients or the depth of the LFH. This is surprising, as snags represent a relatively small (36 Mg ha −1 ) component of the dead biomass within the UFDP with most (62 Mg ha −1 ) provided by litter, duff, and deadwood (Lutz et al. 2021a). In addition to the differences in the size of biomass pools between snags and the LFH and deadwood pools, standing dead trees will have a lower root surface area in contact with the soil as a result of root decay and thus have a smaller portion of their biomass accessible to soil-associated saprotrophs. While our saprotrophic-richness models included the depth of the LFH, a surrogate for saprotroph-available biomass, we did not characterize the relative contribution of each tree species litterfall to the standing LFH layer, which influences saprotrophic communities (Otsing et al. 2018). In general, microbes break down organic matter more quickly when it has low molecular diversity, low spatial heterogeneity, and when organic matter deposition occurs regularly (Lehmann et al. 2020). In Fig. 6 The increase in mean square error (a) and partial dependence plots (b) from a random forest model assessing the influence of the nearest large-diameter (> 40.2 cm diameter at breast height) tree species on ectomycorrhizal fungal species richness at the Utah Forest Dynamics Plot, Utah, USA contrast, the UFDP has several tree species which differ in deadwood residence times by many orders of magnitude and has sporadic, high-volume deposition of woody material from bark-beetle outbreaks or fire every several centuries (Veblen et al. 1994). Constitutive defenses of wood vary widely among UFDP tree species (Runyon et al. 2020;Ott et al. 2021) and may drive known differences in dead wood residence time. The total amount of organic material present may have lesser importance for saprotrophic fungal diversity than its composition (i.e., variability in the types of plant litter deposited), which in turn influences the composition of the resulting organic matter and its physical and chemical stabilization in the soil matrix (Schmidt et al. 2011;Lehmann et al. 2020). Assessing these characteristics of the LFH layer and soil organic matter were beyond the scope of this study but represent directions for future research. The presence of living trees likely indicates the presence of local ectomycorrhizal fungi which may compete with saprotrophic fungi and thereby reduce the decomposition rates of soil organic matter (Gadgil and Gadgil 1971;Fernandez and Kennedy 2016). This apparent competition between ectomycorrhizal and saprotrophic fungi is variable between forests, but may further reduce the ability of saprotrophic fungi to access and decompose organic matter.
Soil nutrients were associated with compositional changes for ectomycorrhizal and saprotrophic fungal communities, but not for the entire fungal community, indicating different proclivities for soil nutrients across the full fungal community. The richness of all fungi was associated with soil component 2, which was correlated with higher manganese and potassium, and lower ammonium. Manganese peroxidase is used by white-rot fungi to decompose lignin (Berg et al. 2015) and is a limiting factor for saprotrophic fungi seeking to decompose organic material (Kranabetter 2019(Kranabetter , 2021. Potassium is a key plant nutrient and ectomycorrhizal fungi can participate in both weathering and transport of potassium to plants in exchange for carbon (Dominguez-Nuñez et al. 2016). Finally, ammonium is a source of nitrogen, a critical nutrient for plants and fungi (Smith and Read 2010). Each of these nutrients could play a role in shaping specific fungal guilds depending on local resource limitations, at least as observed at the time of sampling.

Legacy of forest disturbances on soils, forests, and fungi
The structure and composition of the aboveground forest and the soil environment are a partial product of past disturbances. Disturbance to the overstory composition such as beetle-kill or fire can differentially select for or against certain tree species, thereby influencing the amount, time-scale, and chemistry of organic matter inputs to soil, which in turn provide the substrate for decomposition (Štursová et al. 2014;Williams et al. 2016). The characteristics of these inputs also differ among disturbances; for example, beetle infestations may increase plant litter and deadwood. In contrast, fires can consume deadwood and LFH horizons, as well as organic matter in surface mineral soils (Miesel et al. 2018). Fires can also alter the soil physical environment (Quigley et al. 2019) and the composition of microbial communities (Adkins et al. 2020). Because tree size and age are closely related at the UFDP (Birch et al. 2021), the largest-diameter trees are likely to also be the oldest trees and have survived multiple disturbances. Thus, large-diameter trees at the UFDP may represent survival of tree roots through multiple disturbances over time, leading to greater diversity of soil organic matter and fungal communities.

Fungi and old trees: who selects whom?
Perhaps one of the most interesting questions raised by this and other studies of mycorrhizal fungal communities associated with old trees (Birch et al. 2021) is whether large, old trees create unique fungal communities, or whether specific fungal communities allow trees to persist to old ages and large sizes (Fig. 1). From one perspective, as trees become older, they may be able to be colonized by a more diverse fungal community as a function of random colonization over time. However, as trees age they may support co-evolved fungi that have adapted to late-successional forest conditions. Previous work at the UFDP found that deterministic selection overwhelmingly contributed to fungal community assembly (> 80%) on the roots and in soils around Pseudotsuga menziesii around our sampling area and that drift and dispersal limitations accounted for < 20% of community assembly (Birch et al. 2022b). This suggests that deterministic selection, either from bi-directional tree-fungi selection, or environmental filtering are likely the dominant factors influencing fungal community assembly.
Fungal communities become increasingly dissimilar to one another with tree age, even out to 1400 years old (Birch et al. 2021) and at least a portion of this is likely due to the effect of species richness (Jiang et al. 2021) or within-species genetic diversity (i.e., ploidy of Populus; Bishop et al. 2019). The age of the forest at the UFDP is a remnant signature of multiple fires and intermittent drought over the previous millennia that disturbed heterogeneous patches within the forest (Additional file 3: Fig. S8). A combination of size and age-related changes within trees are likely responsible for explaining the relationship between large trees and fungal richness. Likewise, specific fungal communities may also promote the health and persistence of large-diameter trees which in turn promote forest structural heterogeneity (Lutz et al. 2013). Old-growth trees are of special conservation concern globally (Lindenmayer and Laurance 2016), and the preservation of large, old trees may help preserve unique belowground fungal communities. We used a largediameter threshold above which trees contained half of the forest biomass which can be easily adapted to identify and test the influence of large-diameter trees in other forests . Our results suggest that management actions that conserve large-diameter trees may be an effective way of promoting belowground ectomycorrhizal fungal richness, particularly if trees are spaced within 0.0-16.0 m. The exact density, species, and size of large-diameter trees which are associated with fungal community composition will likely vary by forest and requires further research. Management aimed at promoting structural diversity and deadwood retention has been previously linked with improved fungal richness (Dove et al. 2021;Tomao et al. 2020) and we expect that structurally complex forests will promote a variety of microhabitats and assembly processes that may improve fungal richness.

Conclusion
Hidden below an old-growth, globally unique forest was a diverse, largely unknown fungal community. Fungal functional guilds differed markedly in their spatial structure and the environmental variables explaining their composition and richness. These results highlight the utility of pairing spatially explicit tree and wood mapping with a rigorous sampling of belowground fungi. Of key importance, we found that the location and species of large-diameter trees were excellent predictors of ectomycorrhizal fungal richness and fair predictors of overall fungal richness. This work adds further evidence to the body of literature supporting the outsized ecological influence that large-diameter trees have on their ecosystems and the critical importance of conserving these vulnerable giants.
Additional file 1: Table S1. The basal area, number of stems, and range of ages for each of the woody species within the 1.15 ha forest sampled for fungal communities in Utah, USA. Table S2. The number of fungal sequences at each step of the DADA2 quality filtering pipeline. Table  S3. The component loadings (Comp.1-Comp.5) of soil nutrients from a principal coordinate analysis at the Utah Forest Dynamics Plot, Utah, USA. Table S4. The proportional sequence abundance (%) and number of amplicon sequence variants (ASVs) for the fungal functional guilds within the Utah Forest Dynamics Plot, Utah, USA. Table S5. The Mantel correlation coefficients between distance and community composition of all fungi at the Utah Forest Dynamics Plot, Utah, USA. Bolded values are significant (p <0.05). Table S6. The Mantel correlation coefficients between distance and community composition of all ectomycorrhizal fungi at the Utah Forest Dynamics Plot, Utah, USA. Bolded values are significant (p <0.05). Table S7. The Mantel correlation coefficients between distance and community composition of all saprotrophic fungi at the Utah Forest Dynamics Plot, Utah, USA. Bolded values are significant (p <0.05). Table  S8. The Mantel correlation coefficients between distance and community composition of tree biomass at the Utah Forest Dynamics Plot, Utah, USA. Bolded values are significant (p <0.05). Table S9. The Mantel correlation coefficients between distance and community composition of all deadwood biomass at the Utah Forest Dynamics Plot, Utah, USA. Bolded values are significant (p <0.05). Table S10. The rank of variables in the increase in mean square error, in Random Forest models for fungal diversity, by fungal guild, at the Utah Forest Dynamics Plot, UT, USA Additional file 2. Methods for removal of PCR inhibitors, biomass calculation for deadwood, and sampling of soil nutrients.
Additional file 3: Figure S1. The distribution of fungal sequence count for soil samples taken at the Utah Forest Dynamics Plot, Utah, USA. The red horizontal line denotes the 12,500 sequence cutoff used to rarefy the sequence abundances. Figure S2. The interpolated mineral soil nutrients at the Utah Forest Dynamics Plot, Utah, USA. Figure S3. A principal coordinate analysis (PCA) of soil nutrients at the Utah Forest Dynamics Plot, Utah, USA. Soil nutrients were interpolated across the plot and the values at 117 sampling locations were used for this PCA. Figure S4. The diameter at breast height (1.37 m; cm) for all tree species with > 10 individuals at the Utah Forest Dynamics Plot, Utah, USA. Figure S5. The distribution of (A) observed fungal amplicon sequence variants (ASVs) and (B) Shannon's Alpha diversity for soil samples taken at the Utah Forest Dynamics Plot, Utah, USA. The red horizontal lines in each panel mark the mean value. Figure S6. The species accumulation curves for all fungal amplicon sequence variants (ASV) (a), ectomycorrhizal fungi (b), and saprotrophic fungi (c) at the Utah Forest Dynamics Plot, Utah, USA. Figure  S7. A fitted variogram of the semi-variance of the depth of the LFH layer (litter + duff ) from 295 measurements across the Utah Forest Dynamics Plot, UT, USA. Figure S8. The estimated age of the forest at the Utah Forest Dynamics Plot, Utah, USA. Tree age was estimated using kriging interpolation of 127 trees. Figure S9. The increase in mean square error (a) and partial dependence plots (b) from a random forest model assessing the influence of the nearest large (>40.2 cm diameter at breast height) tree species on total fungal species richness at the Utah Forest Dynamics Plot, Utah, USA. Figure S10. The increase in mean square error (a) and partial dependence plots (b) from a random forest model assessing the influence of the nearest large (>40.2 cm diameter at breast height) tree species on saprotrophic fungal species richness at the Utah Forest Dynamics Plot, Utah, USA.