Regional-Scale In-Depth Analysis of Soil Fungal Diversity Reveals Strong pH and Plant Species Effects in Northern Europe

Soil microbiome has a pivotal role in ecosystem functioning, yet little is known about its build-up from local to regional scales. In a multi-year regional-scale survey involving 1251 plots and long-read third-generation sequencing, we found that soil pH has the strongest effect on the diversity of fungi and its multiple taxonomic and functional groups. The pH effects were typically unimodal, usually both direct and indirect through tree species, soil nutrients or mold abundance. Individual tree species, particularly Pinus sylvestris, Picea abies, and Populus x wettsteinii, and overall ectomycorrhizal plant proportion had relatively stronger effects on the diversity of biotrophic fungi than saprotrophic fungi. We found strong temporal sampling and investigator biases for the abundance of molds, but generally all spatial, temporal and microclimatic effects were weak. Richness of fungi and several functional groups was highest in woodlands and around ruins of buildings but lowest in bogs, with marked group-specific trends. In contrast to our expectations, diversity of soil fungi tended to be higher in forest island habitats potentially due to the edge effect, but fungal richness declined with island distance and in response to forest fragmentation. Virgin forests supported somewhat higher fungal diversity than old non-pristine forests, but there were no differences in richness between natural and anthropogenic habitats such as parks and coppiced gardens. Diversity of most fungal groups suffered from management of seminatural woodlands and parks and thinning of forests, but especially for forests the results depended on fungal group and time since partial harvesting. We conclude that the positive effects of tree diversity on overall fungal richness represent a combined niche effect of soil properties and intimate associations.

Soil microbiome has a pivotal role in ecosystem functioning, yet little is known about its build-up from local to regional scales. In a multi-year regional-scale survey involving 1251 plots and long-read third-generation sequencing, we found that soil pH has the strongest effect on the diversity of fungi and its multiple taxonomic and functional groups. The pH effects were typically unimodal, usually both direct and indirect through tree species, soil nutrients or mold abundance. Individual tree species, particularly Pinus sylvestris, Picea abies, and Populus x wettsteinii, and overall ectomycorrhizal plant proportion had relatively stronger effects on the diversity of biotrophic fungi than saprotrophic fungi. We found strong temporal sampling and investigator biases for the abundance of molds, but generally all spatial, temporal and microclimatic effects were weak. Richness of fungi and several functional groups was highest in woodlands and around ruins of buildings but lowest in bogs, with marked group-specific trends. In contrast to our expectations, diversity of soil fungi tended to be higher in forest island habitats potentially due to the edge effect, but fungal richness declined with island distance and in response to forest fragmentation. Virgin forests supported somewhat higher fungal diversity than old non-pristine forests, but there were no differences in richness between natural and anthropogenic habitats such as parks and coppiced gardens. Diversity of most fungal groups suffered from management of seminatural

INTRODUCTION
Soil microorganisms such as bacteria, archaea and fungi play integral roles in soil ecosystem functioning. In particular, fungi act as principal decomposers of organic material, key root mycorrhizal symbionts of plants as well as harmful pathogens in forest, grassland and agricultural habitats. It has been argued that the diversity of soil organisms is of particular importance for ecosystem services and resilience to disturbances, with a potential to ameliorate stress caused by global change (Langenheder et al., 2010;Bardgett and van der Putten, 2014). Many of these soil organisms are vulnerable to shifts in land use, changing climate and ecosystem management (Van der Putten, 2013;George et al., 2019;Makiola et al., 2019;Sterkenburg et al., 2019); therefore, some of the most conspicuous rare species are protected (Dahlberg et al., 2010). However, fungal species that have large fruitbodies represent only a tip of the iceberg in the enormous diversity of soil mycobiome, much of which is cryptic and previously undocumented. Advances in the socalled environmental DNA (eDNA) analyses have enormously improved our understanding about the diversity and distribution of small soil organisms including fungi (Taberlet et al., 2018;Nilsson et al., 2019).
Virtually all members of the soil biota exhibit distinct distribution patterns from microscale to global scale in relation to soil parameters, vegetation and climate. These environmental effects are also linked to spatiotemporal patterns caused by seasonal climate variability as well as dispersal limitation at the landscape and continental scales Cavicchioli et al., 2019). While several global efforts have provided rough information about the distribution of soil microbes on the Earth, where climate typically drives large-scale patterns (Tedersoo et al., 2014;Bahram et al., 2018;Delgado-Baquerizo et al., 2018), we still lack comprehensive studies from landscape to regional scales that are more relevant to understanding the immediate effects of vegetation, land use and edaphic properties on soil organisms (but see Karimi et al., 2018;George et al., 2019). The regional geographic scale is important to consider, because this is typically little affected by historical dispersal limitation and climate (except montane systems), especially the covariation of climate with vegetation and soil properties -all inherent to larger geographic scales (Lomolino et al., 2006). Multiple highly focused studies have demonstrated the effects of land use, dominant plant species, soil properties, pollution and management by humans on soil microorganisms (Rousk et al., 2010a;Hartmann et al., 2012;Urbanova et al., 2015;Creamer et al., 2016;Sterkenburg et al., 2019). However, these studies do not enable us to understand the relative importance of these factors because of their limited ecological scope.
Using third-generation PacBio sequencing and 1251 standard composite samples, we provide a detailed view on the biotic and abiotic factors underlying the distribution of soil fungi, across different habitat and vegetation types, ecosystem management regimes and edaphic gradients in the northern Baltic region (Figure 1). Our main purpose was to evaluate differences in the response of taxonomic and functional groups of fungi to various abiotic and biotic variables by explicitly accounting for spatial and temporal autocorrelation, microclimatic differences and technical biases in a multi-year sampling effort. We tested the following hypotheses: (1) diversity of biotrophic groups responds strongest to dominant vegetation, whereas free-living groups are more affected by abiotic variables; (2) diversity of soil fungi is relatively lower in island and fragmented forest ecosystems, especially in small and distant habitats following the island biogeography theory; (3) anthropogenic ecosystems harbor lower mycobiome diversity than natural and particularly pristine habitats; (4) management of forests, seminatural woodlands and parks reduces soil fungal biodiversity.

Sampling and Pre-treatment
Within the period of July, 2011 to May, 2018, we collected 1251 composite soil samples comprising altogether > 50,000 individual samples in the projected area of 80,000 km 2 in the northern Baltic region (Estonia and North Latvia; Supplementary Table S1 and Figure 1). The sites were selected based on aerial photographs, forestry and protected areas databases and haphazard visitation to cover a great variety of habitats. We also focused on old and protected habitats and monocultures of tree stands to improve our capacity to detect the effects of age and dominant plant species. Selection of the sites followed the basic criteria: (1) homogeneous area with no visible ecotones (except islands) and vegetation gradients; (2) >30 m distance from other vegetation type or trees belonging to species not present in the plot (to avoid influence by roots or litter), cutting of no more than 40% of trees basal area; and (3) no legal constraints.
As described in Tedersoo et al. (2014), we established circular 2500-m 2 plots in homogeneous vegetation (avoiding ecotones) in forests, bogs, grasslands, croplands, parks and around abandoned buildings (ruins) of contrasting age and plant species composition (for ecosystem types and measured variables, see Supplementary Table S2). From each four quadrats of the plots, five trees located at least 8 m apart were randomly selected. From two opposing sides of each 20 trees per plot, 1-1.5 m from the tree trunk, soil cores (5 cm diam. to 5 cm depth) were collected using a sterilized PVC tube or sharp knife. In non-wooded habitats, randomly located spots were selected instead of trees for collecting samples. About 5-10 g of soil from the margins of each core was placed into a zip-lock plastic bag, with the interior of the soil core placed back into the hole. The material from all 40 cores per plot was pooled into the same bag, without removing fine roots or small stones (<1 cm diam.). The composite soil sample was laid on clean newspaper and air-dried as soon as possible, within at least 24 h since collection, at <40 • C in a dry room, with heat from the sun, oven, light bulbs, floor heating system or other type of heating body. The newspaper was replaced when becoming moist to speed up the drying process. In a few exceptional cases (<2%), the sampling area was smaller (small island habitats) or larger (sparse woodlands) than 2500 m 2 (for exceptions, see Supplementary Table S1).

Molecular Analysis and Bioinformatics
The dried composite soil samples were kept in zip-lock plastic bags and homogenized manually by intensively rubbing the bag and its contents by hands for 3 min. DNA was extracted from 2.0 g of soil dust using the PowerMax Soil DNA Isolation kit (Qiagen, Carlsbad, CA, United States) following manufacturer's instructions. The DNA extracts were further purified using FavorPrep TM Genomic DNA Clean-Up Kit (Favorgen, Vienna, Austria).
The purified DNA extracts were subjected to amplification with the universal eukaryotic primers ITS9mun and ITS4ngsUni that have been nominated for global microbiome analyses of eukaryotes at species-level resolution (Tedersoo and Lindahl, 2016;Tedersoo and Anslan, 2019). These primers amplify 170 bp of the 18S rRNA gene and full-length Internal Transcribed Spacer (ITS), which serves as an official barcode for fungi (Schoch et al., 2012;Nilsson et al., 2018) and performs well on many groups of soil animals and protists (Coleman, 2009;Pawlowski et al., 2012). We selected this partial 18S + ITS marker for approximating operational taxonomic units (OTUs) at roughly species level, increased taxonomic resolution and higher capacity to detect and remove artifacts compared with regular short ITS1 or ITS2 and partial 18S markers (Tedersoo et al., 2018b). In spite of a high singlepass error rate, advances in PacBio third-generation sequencing technology enable producing circular consensus sequences (CCS) of >3000 bp with high quality.
For amplification, the PCR mixture comprised 5 µl of 5x HOT FIREPol Blend Master Mix (Solis Biodyne, Tartu, Estonia), 0.5 µl of each forward and reverse primer (20 mM), 1 µl of DNA extract and 18 µl ddH2O. Thermal cycling included an initial denaturation at 95 • C for 15 min; 25-30 cycles of denaturation for 30 s at 95 • C, annealing for 30 s at 57 • C, elongation for 1 min at 72 • C; final elongation at 72 • C for 10 min and storage at 4 • C. PCR products were normalized for library preparation and sequenced on PacBio Sequel instrument using SMRT cell 1M, v2 LR; Sequel Polymerase v2.1 and Sequencing chemistry v2.1. Loading was performed by diffusion, one SMRT cell was used for sequencing with a movie time of 600 min and pre-extension time of 45 min.
PacBio CCS reads (minPasses = 3, MinAccuracy = 0.9) were generated using SMRT Link v 6.0.0.47841. Subsequent quality filtering was performed using PipeCraft 1.0 (Anslan et al., 2017) as described in Tedersoo and Anslan (2019). UNITE 7.2 (Kõljalg et al., 2013) served as a reference database for chimera filtering and identification. Using UPARSE (Edgar, 2013), sequences were clustered to OTUs at 98% sequence similarity to enable separation of closely related species while keeping lower-quality sequences and rare variants adhered to OTUs' centroids. Global singletons were removed. The entire workflow is given in Supplementary Table S3. For taxonomic assignments, we evaluated the 10 best hits and conservatively kept OTUs with conflicting best matches unidentified at the level of that taxonomic rank. The criteria for kingdom-level treatment included e-value < e −20 and sequence similarity > 70% (BLASTn values); that for phyla and class was e-value < e −50 and sequence similarity > 75%, whereas orders, families and genera were assigned at >80%, >85%, and >90% sequence similarity, respectively. Taxa falling below these thresholds were considered unknowns at particular taxonomic levels. Fungal taxa [taxonomy follows Tedersoo et al. (2018a) as updated in Wijayawardene et al., 2020] were used at the level of phyla (early diverging lineages), classes (Chytridiomycota, most Mucoromyceta) and orders (Dikarya, molds) to balance between taxonomic resolution and coverage.
To assign functions to fungal OTUs, we used two parallel approaches. First, we used the newly built FungalTraits database (Põlme et al., unpublished; beta version available at doi: 10.15156/BIO/807446) to assign OTUs to guilds (Supplementary Table S4) and EcM fungi further to lineages and exploration types. For genera with multiple lifestyles, we used the assignments based on annotations at the level of sequences and species hypotheses (SHs) as given in UNITE 8.1.

Environmental Variables
During sampling, we recorded all woody plant and dominant herb species. Except for ericoid subshrubs, we estimated the relative basal area using direct measurements  or visual estimates, for which participants were trained. For tree species, we considered presence when saplings were >1 m high. For non-wooded habitats, we estimated the proportion of plant species based on their coverage. For small-sized and rare species with no or minimal basal area, relative abundance of 0.1% was assigned. The age of vegetation was determined based on forestry databases, annual growth of trees based on branching, or circumference of largest trees at breast height (for deciduous trees > 200 cm circumference, 2 cm corresponds to 1 year based on increment cores). According to the evidence for mowing or grazing and growth of coppice, we estimated the management status for woodlands and parks (categories: managed, unmanaged and coppiced). Based on forestry databases or visual estimates, we determined the proportion of trees harvested and time since harvesting for selectively harvested or thinned stands. We considered old forests as virgin habitats, when these were at least 120 years old according to the oldest trees and were databased as old-growth forests or harbored large amounts of woody debris in all stages of decomposition and showed no visible evidence for historical human interference (e.g., no visible stumps).
Chemical properties of each composite sample were measured from 20 g of dried, homogenized material. Soil pH KCl , P, K, Mg, and Ca concentration as well as 15 N and 13 C abundances were determined following Tedersoo et al. (2012a). Concentrations and ratios of C/N, C/P and N/P were log-transformed for statistical analyses.
The climatic conditions in the region displayed little variation in mean annual temperature (5-7 • C) and precipitation (500-750 mm). Although we expected that this has little effect on soil biota, we included a set of microclimatic variables (Karger et al., 2017) in our analysis (Supplementary Table S2).
Based on geographical coordinates of sample centroids, we calculated Moran's Eigenvector Maps (MEMs) to represent geographic distance vectors (gMEMs; n = 319) as implemented in adespatial package (Dray et al., 2018) of R (R Core Team, 2019). We also calculated 74 temporal distance vectors (tMEMs) based on sampling dates. We treated year and month of sampling (February to April pooled; November to January pooled) as dummy variables and considered these and square-root-transformed linear time as covariates.
Phylogenies of woody plant species were adapted from Zanne et al. (2014). Based on an ultrametric tree, we calculated nearest neighbour (comdistnt), average neighbor (comdist) and phylosor distances for all samples as implemented in picante package (Kembel et al., 2010). Using adespatial package, we generated the respective 32 ntMEMs, 44 cdMEMs and 119 psMEMs based on these distances for use as covariates.
Based on estimates of basal area (or herb/shrub coverage), we calculated the relative abundance of EcM plants (Soudzilovskaia et al., 2020). Proportion of EcM plants and individual species were log-ratio-transformed (+0.01%) to vary from −4 (0.0%) to 0 (50.0%) to 4 (100.0%). The age of vegetation, number of woody plant species and EcM plant species sampled were used in both untransformed and square-root-transformed formats.
We also included dummy variables for collector, habitat type, selective harvesting or thinning, island type, management type (woodlands and parks only), presence of water bodies, and anthropogenic habitat (wild, village, urban). In more specific analyses and for illustration, these variables were used as multilevel categorical predictors (Supplementary Table S2).
As a proxy for fungal diversity, we used both OTU richness and Shannon diversity index. We also calculated OTU richness for dominant functional guilds, taxonomic groups and EcM lineages of fungi (Supplementary Table S4). All OTU richness measures were converted to residuals based on the average values of raw residuals taken from regression analyses of OTU richness vs. square-root-transformed sequencing depth and log-transformed sequencing depth. Averaging was used, because these functions tended to overestimate residual richness at either low or high sequencing depth. The relative sequence abundance-based proportions of fungal guilds, mold groups, EcM exploration types and individual OTUs were log-ratio-transformed.

Pre-statistics Data Quality Evaluation
Before statistical analyses, we sorted the OTU matrices by negative and positive controls to seek for potential laboratory and technical contamination. The initial dataset#1 included 3,060,546 quality-filtered sequences from 1310 samples that were clustered into 23,787 OTUs (Table 1). We removed 59 control samples and 40 samples with low sequencing depth (<500 reads). Based on analysis of sequence distribution in control samples, tag switching and non-targeted markers, 28 OTUs were specifically excluded.
The resulting dataset#2 (1211 samples, 2,932,102 sequences and 23,464 OTUs) comprised on average 2421.2 sequences (SD = 1242.6; min = 504; max = 12763). We observed that a considerable proportion of samples revealed a high proportion of reads from molds (rapidly growing opportunistic saprotrophs), either a single high-abundance OTU or multiple OTUs from one or several phylogenetic groups. Previous studies suggested that high mold abundance is at least partly related to poor condition or poor preservation of samples (Tedersoo et al., 2014;. Due to high climatic variability or small sample size, we did not address this issue in-depth previously. Molds assigned to Mortierellales, Umbelopsidales, Mucorales and Pezizomycotina (mostly Eurotiales) contributed to 6.9 ± 4.8% (mean ± SD), 3.2 ± 7.5%, 0.2 ± 0.6%, and 2.4 ± 9.0% of sequences per sample, respectively, with maximum values between 12 and 90%. The negative effect of molds was tested on residuals of total fungal OTU richness, EcM fungal OTU richness and proportion of EcM fungal sequences, because we expected the latter obligately mutualistic group to be most responsive to molds (Tedersoo et al., 2014). The relative abundance of Umbelopsidales molds and proportion of most abundant OTU of Umbelopsidales had strongest negative correlations with overall fungal richness (R = −0.437; P < 0.001 and R = −0.325; P < 0.001, respectively), followed by negative effects of relative abundance of the most common Pezizomycotina mold OTU (R = −0.213; P < 0.001) and Mortierellales OTU (R = −0.140; P < 0.001). Residuals of EcM fungal richness were negatively affected by relative abundance of the most common Umbelopsidales OTU (R = −0.300; P < 0.001), Mortierellales OTU (R = −0.222; P < 0.001) and Mucorales OTU (R = −0.164; P < 0.001) as well as relative proportion of Mortierellales sequences (R = −0.189; P < 0.001). Random forest modeling (see below) confirmed that Mortierallales and Eurotiales mold relative abundances are among the strongest predictors of fungal richness, EcM fungal richness and EcM fungal proportion (Supplementary Figure S1). Accordingly, we removed 36 samples that were heavily occupied by molds (maximum OTU relative abundance > 30% in Pezizomycotina, >15% in Umbelopsidales, or >15% in Mortierellales molds).
Further analyses of molds revealed that relative abundance of Mortierallales and Umbelopsidales peaked at low-pH and C-rich forest soils, whereas Mucorales showed occasional high values in croplands. Pezizomycotina mold abundance was unrelated to any biological variables, suggesting that these molds are the most opportunistic but avoid samples where Mucoromyceta molds, especially Mortierellales are abundant. As relative abundances of all mold groups had strong temporal and/or collector biases, we conclude that their high abundance is related to either unfavorable meteorological conditions during the time of sampling (e.g., freeze-thaw and rewetting of dry soils) and/or differential post-sampling sample spoilage (too slow or incomplete drying) in hands of collectors. More detailed information about the environmental and technical variables driving mold abundance is given in Supplementary Figure S2 and Supplementary Table S5.

Statistical Analyses
We used the mold-corrected dataset#3 for the main analyses (Table 1), including the proportion of molds as covariates because of their potential residual effects. First, we used the random forest machine learning algorithm to generate non-linear models for all richness variables using combined features of randomforest (Liaw and Wiener, 2002) and VSURF (Genuer et al., 2019) packages in R. Model evaluation was performed using 999 trees. Because these analyses provide no information about the type of fit or determination coefficient, we tested these pre-selected best predictors by fitting quadratic functions and general linear modeling for each dependent variable as implemented in Statistica 13 (TIBCO Software Inc., NY, United States). For EcM fungi, we removed 64 non-EcM plots and samples containing < 100 EcM fungal sequences (dataset#4). The above-described model selection procedure was used for richness residuals of EcM fungi and EcM fungal lineages and relative abundance of exploration types. We also calculated niche models for 50 most frequent EcM and non-EcM fungal OTUs using log-ratio-transformed relative abundances. We further used piecewise structural equation modeling to delve into direct and indirect relationships between OTU richness and the key potential predictors, as implemented in piecewiseSEM package (Lefcheck, 2016) of R. Richness residuals of fungi ( Figure 2) and major fungal guilds were extrapolated to a land use and forest map of Estonia at 250 × 250 m resolution based on available map layers. The estimates of stand productivity were obtained from State Registry of forest data and soil map (Environment Agency, Tartu 1 ). Proportions of tree species (Lang et al., 2018) and stand height, volume and stand density (Lang et al., 2012) were calculated by equations from the Regulation of Ministry of Environment (Regulations of the Forest Survey, 2018; Regulation of Minister of the Environment No 2. Enforced 16.01.2009.). The estimated soil acidity and humidity scalars were adopted from the Estonian Soil Classification for Postlithogenic normally developed mineral soils (Astover et al., 2012) by converting from specific average values of soil types 2 . The extrapolated fungal richness values were transferred to the Estonian base map 3 .
To infer the most influential variables affecting the composition of fungi, we constructed multivariate PERMANOVA models (Anderson et al., 2008), followed by supporting analyses using Non-metric Multidimensional Scaling (NMDS) in vegan package of R (Oksanen et al., 2019) and estimates of potential non-linear effects as implemented in General Dissimilarity Modeling (GDM; Manion et al., 2018). For these multivariate analyses, environmental variables were normalized, abundances of OTUs were Hellinger-transformed, and Bray-Curtis dissimilarity was used as a distance metric. We assessed three broad groups of fungi: (1) EcM fungi including log-ratio-transformed proportions of lineages in a separate analysis; (2) putative plant pathogen guilds combined; (3) saprotrophs combined (excluding taxa also classified as plant pathogens). In all cases, samples with <20 OTUs in the particular categories were excluded from analyses. For the NMDS graphs, we fitted the environmental predictors and calculated 95% confidence ellipses using the vegan functions envfit and ordiellipse.
The observed taxonomic and functional distribution of fungi corroborates global-scale patterns (Tedersoo et al., 2014;Egidi et al., 2019;Vetrovsky et al., 2019), with somewhat higher proportions of EcM fungi and early diverging lineages. EcM fungal dominance is attributable to prevalence of EcM treedominated forest habitats in our sampling and in North European temperate and boreal forest biomes in general, compared with most other global biomes, especially drylands (Soudzilovskaia et al., 2019). The relatively greater diversity of early-diverging lineages is at least partly related to the use of pan-eukaryotic primers and a marker comprising fulllength ITS and flanking 18S V9 variable region (see section "Methodological Implications"), but it may also be related to their greater beta diversity.

Diversity of Fungi and Major Fungal Guilds
By combining the random forest machine learning and general linear modeling (dataset#3), we found that the overall proportion of fungi among eukaryote sequences responded positively to EcM plant proportion (R 2 = 0.098; Figure 5) and negatively to soil P concentration (R 2 = 0.096). Soil pH was the main predictor for richness of most fungal functional and taxonomic groups (Table 2, Figure 6, and Supplementary Figure S3). Soil pH (unimodal relationship) and plant richness (positive effect) had strongest effects on fungal OTU richness and Shannon index (Figures 5, 6), with an additional effect of community age (sigmoid relationship). Structural equation models (SEMs) confirmed these patterns and suggested that vegetation age has positive effects on fungal richness directly and through increasing plant diversity (Figure 7).
Ectomycorrhizal fungal richness was most strongly related to soil Ca concentration (positive effect; R 2 = 0.260), followed by relative abundance of Betula spp. (R 2 = 0.213), Corylus avellana (R 2 = 0.136) and EcM plants taken together (R 2 = 0.135), and soil pH (R 2 = 0.084; all unimodal responses; Figure 7). A more detailed analysis showed that the relative abundance of EcM fungi reached a plateau at around 50% EcM plant proportion ( Figure 5). Conversely, richness of the AM Glomeromycota was negatively related to EcM plant basal area (R 2 = 0.237) and soil C/N ratio (R 2 = 0.024). Relative abundance of Glomeromycota declined linearly from 0 to 90% EcM plant basal area, with a sharper drop beyond 90% (Figure 5), whereas Glomeromycota richness decreased linearly throughout the range of EcM plant proportions. By contrast, relative abundance and richness of putatively AM Endogonomycetes (including a single rare EcM fungal OTU) increased at around 98% of EcM plant proportion ( Figure 5) that coincided with acidic pine forest habitats. There was a negative correlation between the abundance of Endogonomycetes and Glomeromycota (R = −0.205), which we attribute to their contrasting ecological preferences rather than competition. Although several Endogonomycetes OTUs matched closely to the AM fine root endophytes recovered in previous studies (Orchard et al., 2017), our results provide circumstantial evidence that most of these taxa are not obligate AM symbionts. Our results furthermore suggest that AM Glomeromycota are more sensitive to relative abundance of their host plants compared with EcM fungi, which maintain high biomass and richness at lower proportions of their host plants. These differences indicate that soil and ecophysiological processes such as slow nutrient cycling and neutral or positive plant-soil feedback and low bacterial diversity and activity characteristic of EcM habitats (Phillips et al., 2013;Tedersoo and Bahram, 2019;Bahram et al., 2020;Tedersoo et al., 2020), may also prevail in mixed forests.  The best GLM predictors are given in bold. All effects are significant (P < 0.001). 1 Symbols in continuous predictors: , near-linear increase; , near-linear decline; , sigmoid increase; , cumulative increase; ∩, unimodal. Subgroups of response variables are indicated in bold.
Richness of all saprotrophs taken together was best predicted by proportion of EcM plants (negative effect, R 2 = 0.261) and soil pH (peak at 5.5-7; R 2 = 0.083), but best predictors differed among guilds ( Table 2). Litter decomposers prevailed at near-neutral pH (pH 6-7; R 2 = 0.354) and responded negatively to proportion of EcM plants (R 2 = 0.088), particularly Picea abies (R 2 = 0.031), with synergistic effects between pH and tree composition (Figure 7). Soil Mg concentration (R 2 = 0.195) and proportion of EcM plants (R 2 = 0.071) best explained wood saprotroph richness, both with unimodal relationships. Diversity of white-rot decomposers responded negatively to soil δ 15 N (R 2 = 0.084) but positively to proportion of Alnus spp. (R 2 = 0.073) and vegetation age (lag at 100 years; R 2 = 0.047). Soil saprotrophs responded negatively to soil pH (R 2 = 0.132) and positively to age (lag at ca 140 years; R 2 = 0.074). Dung saprotrophs had a pH optimum at pH > 6 (R 2 = 0.419) and positive response to soil δ 15 N (R 2 = 0.072). The strong negative effect of EcM plant relative abundance on litter saprotrophs is attributable to reduced soil pH and competition between EcM and saprotrophic fungi for nutrients bound in soil organic material (Sterkenburg et al., 2018).
Leaf pathogens prevailed at weakly acidic to neutral soils (R 2 = 0.241) and low relative abundance of EcM plants (R 2 = 0.054), in particular Picea abies (R 2 = 0.044). Wood pathogens were positively related to woody plant richness (R 2 = 0.141) and soil pH (peak at 4.5-6.5; R 2 = 0.057). Animal parasites were most diverse in neutral soils (pH > 6; R 2 = 0.219) with high P concentration (R 2 = 0.063). Richness of OTUs classified as opportunistic human pathogens similarly prevailed at high pH (pH > 6; R 2 = 0.132) but exhibited strong temporal shifts (R 2 total = 0.128). Similarly, mycoparasite richness was mostly explained by various temporal vectors (R 2 total = 0.122), with weak unimodal responses to soil pH (R 2 = 0.032) and Ca concentration (R 2 = 0.026). Similarly to molds, mycoparasites and fungal decomposers take advantage  of nutrients released from dying mycelium, which is related to annual and intra-annual cycles, confirming the temporal patterns in various biomes Castano et al., 2018b;Burke et al., 2019). The positive effects of soil pH and negative effects of EcM vegetation on plant pathogens support the observations that pathogens drive plant community dynamics in AM-dominated forests (Bennett et al., 2017;Chen J. et al., 2019).
Yeasts were dominated by members of Basidiomycota (Tremellomycetes and Pucciniomycotina). The richness of this entire life form increased with vegetation age (R 2 = 0.081). Dimorphic yeasts, mostly represented by members of Herpotrichiellaceae, were negatively affected by relative abundance of Picea abies (R 2 = 0.061).
In Ascomycota, the orders Glomerellales, Hypocreales, Sordariales, Xylariales (all Sordariomycetes), Pezizales and Pleosporales were more diverse in soils with neutral pH, whereas Helotiales, Thelebolales (both Leotiomycetes) and Venturiales had pH optima at moderately or weakly acidic soils ( Table 2 and Supplementary Figure S3). Archaeorhizomycetales richness was associated with strongly acidic soils in pine-dominated forests. Interestingly, Saccharomycetales yeasts were most diverse at both ends of the pH gradient, confirming their realized niche in extreme habitats (Kurtzman et al., 2011). Richness of Chaetosphaeriales and Eurotiales was positively related to the relative abundance of Quercus robur (Supplementary Figure S4). In Basidiomycota, the Tremellomycetes orders Cystofilobasidiales and Filobasidiales were relatively more diverse at neutral pH, whereas Trichosporonales prevailed at moderately acidic habitats and Tremellales had no particular pH preference ( Table 2). Of Agaricomycetes, Trechisporales, Hymenochaetales, Polyporales and Tremellodendropsidales were associated with highly or moderately acidic soils, but non-EcM OTUs of Agaricales prevailed in weakly acidic soils (Supplementary Figure S3). The saprotroph/necrotrophic pathogen family Ceratobasidiaceae was most diverse at high δ 15 N in non-forest habitats.
Although plant diversity was one of the strongest predictors of overall fungal diversity, it contributed directly to OTU richness of non-EcM Agaricales (R 2 = 0.203) and Helotiales (R 2 = 0.062) only. SEM revealed that tree richness had a direct positive effect on fungal richness, which was further promoted by vegetation age, soil pH and Ca that all exhibited direct residual effects as well (Figure 7). These results suggest that the positive effect of plant richness previously observed on local  and regional (Yang et al., 2017) scales represents a synergistic interaction among multiple coevolutionary as well as biotic and abiotic niche differentiation processes that favor richness of both groups. On a global scale, historical processes, dispersal limitation and differential tolerance to climatic conditions have effectively blurred the local-scale causal links and correlation between plant and fungal richness (Tedersoo et al., 2014).
Somewhat unexpectedly, relative abundances of EcM exploration types responded poorly to soil nutrients, with strongest effects of pH instead. Soil pH had a positive effect on short-distance-delicate type (R 2 = 0.447) but a negative effect on contact type (R 2 = 0.299). Soil δ 15 N had a negative effect on medium-distance-fringe type (R 2 = 0.275), whereas mediumdistance-smooth type was most abundant at medium values of soil C/N ratio (R 2 = 0.140). Long-distance type increased linearly with vegetation age (R 2 = 0.123), but mat type displayed negligible response to any environmental variables. Distribution of exploration types in NMDS ordination space relative to other variables is given in Supplementary Figure S6.
The lineage-level taxonomic composition of EcM fungi is consistent with previous comprehensive root-tip-based surveys in temperate Eurasia and North America, where forests of Pinaceae or Fagales dominate (Tedersoo et al., 2012a;Miyamoto et al., 2018;van der Linde et al., 2018). The overall commonness of/inocybe lineage and low abundance of/cortinarius, the rhizomorph-forming group and most species-rich EcM fungal genus in the Nordic countries (Knudsen and Vesterhold, 2012), is somewhat unexpected. Sequence-based dominance of/inocybe belowground and its large abundance of OTUs certainly indicate that this genus is more species-rich than expected, perhaps on the account of multiple cryptic species that we observed particularly for Inocybe maculata and I. rimosa species complexes (Larsson et al., 2009). The richness of/cortinarius but also/hebelomaalnicola lineages may be severely underestimated by the chosen 98% sequence similarity threshold, as multiple closely related species in these groups differ by only a few bases across the entire ITS region (Schoch et al., 2012;Garnica et al., 2016;Grilli et al., 2016).
Exploration types of EcM have been viewed as ways of obtaining nutrients and plant carbon investments into mineral nutrition by differences in soil exploration, potential enzymatic activities and access to nutrients in different forms (Agerer, 2001;Jones et al., 2010;Tedersoo et al., 2012b). Across various ecosystems, our analyses revealed that relative abundance of only medium-distance-fringe and medium-distance-smooth types are related to 15 N abundance or C/N ratio, whereas other types were unresponsive or responded to soil pH. Although pH affects nutrient availability, the distribution of exploration types may simply reflect that of dominant EcM fungal genera (e.g., Russula for contact type) and entire lineages (e.g., /tomentellathelephora for medium-distance-smooth type), because the major morphological features of mycelium and rhizomorphs vary little at these taxonomic levels, being phylogenetically conserved. Hence, the confounding effects of phylogeny and exploration type need to be carefully disentangled.

Structure of Fungal Communities
Multivariate PERMANOVA modeling revealed that community structure of saprotrophs, pathogens (both dataset#3) and EcM fungi (dataset#4a) was most strongly driven by soil pH, which explained 11.4, 6.4, and 5.3% of variation, respectively (Supplementary Table S8). The proportion of Picea abies, soil P and Mg concentrations had consistently significant effects on composition of fungal guilds, but these variables explained only 1.2-2.5% and 0.8-1.6% and 0.4-0.8% of variation, respectively. General dissimilarity modeling (GDM), which accounts for non-linearity in continuous variables, supported these results, but suggested relatively stronger complementary effects of soil C/N ratio, δ 15 N and Mg concentration to the pH effect ( Figure 8A and Supplementary Figure S7). While the relative effect strength of soil pH was nearly linear for composition of saprotrophs and EcM fungi, it was cumulative for plant pathogens, indicating lower importance at near-neutral soils. Soil δ 15 N and C/N ratio effects were generally cumulative, suggesting that low N availability and N occurrence in mostly organic form are relatively more influential. Conversely, the effect strength of soil P and cation concentrations increased exponentially, indicating relatively stronger effects of surplus nutrients (cropland and urban habitats). Alternatively, these patterns may reflect an interaction of nutrients with soil pH, which determines their mobility and availability. The effect strength of tree species ranged from cumulative to linear to exponential, with a pronounced exponential effect of Alnus spp. on EcM fungal composition ( Figure 8A). This is in agreement with specialized communities of Alnus spp. relative to other host plants, which become more prominent in monospecific stands (Kennedy et al., 2015). Conversely, cumulative effects can be interpreted as initial positive effects of critical abundance of a tree species or nutrient level until a certain threshold, beyond which further increase has a negligible effect.
As we obtained samples from multiple monospecific and mixed stands of various Salix species, we tested the effect of Circles, 95% ordiellipses for tree taxa; arrows, fit with environmental variables and placement of EcM fungal lineages and exploration types in top-right pane (a-t, /amphinema-tylospora; hygr, g-h, /genea-humaria; /hygrophorus; m-pg, /marcelleina-peziza gerardii; p-a, /pachyphloeus-amylascus; p-g, /paxillus-gyrodon; p-s, /pisolithus-scleroderma; pst, /pseudotomentella; colors represent lineages belonging to different orders or groups of orders) and mid-right pane (LDD, long-distance; MDF, medium-distance fringe; MDS, medium-distance smooth; SDC, short-distance coarse; SDD, short-distance delicate; colors depict different categories of variables), respectively. Bottom-left pane of (A) indicates the cumulative effect of environmental variables in the relative scale of minimum (left) to maximum (right) values as revealed from GDM. In (B), symbols and their colors indicate different subgenera and species of Salix, respectively; cdMEMs and psMEMs represent phylogenetic eigenvectors indicating phylogeny effects on composition.
host, their phylogenetic associations and other environmental predictors on EcM fungal composition for Salix-dominated (relative abundance > 66%) data subset (n = 60). The best model suggested that Salix species, their phylogenetic relationships, soil pH, other soil variables and other vegetation variables explain 9.5% (F 2,48 = 8.2), 2.4% (F 3,48 = 3.2), 6.2% (F 1,48 = 6.9), 2.5% (F 3,48 = 3.1), and 2.1% (F 2,48 = 3.8) of variation, respectively. The NMDS ordination indicates that the main differences occur between Salix subgenus Vetrix (S. caprea) and subgenus Salix, with no clear differences among S. alba, S. pentandra and S. fragilis belonging to the latter group ( Figure 8B). Our analyses suggest that although the overall EcM fungal composition is driven by soil pH, individual tree genera and subgenera (in Salix) have a strong residual effect. We confirmed the effect of phylogenetic relatedness within the genus Salix but not at higher taxonomic levels . Ordination graphs show that genera of Pinaceae and Fagales occur in separate positions, with non-overlapping confidence intervals ( Figure 8A).
Our analyses indicate that soil pH is the main driver of composition of different fungal guilds, followed by the proportion of Picea abies. This uniform residual effect of one tree species across fungal guilds is unexpected, suggesting its effect on soil properties other than pH or a unique biotic environment through direct interactions with EcM fungi, saprotrophs and pathogens. These non-exclusive, perhaps synergistic mechanisms are plausible based on mycological literature (Põlme et al., 2018;Wang et al., 2019).

Environmental Niche for OTUs
We modeled the realized environmental niche for 50 most frequent OTUs of EcM fungi and 50 non-EcM fungi. Relative abundance of EcM fungal OTUs was mainly driven by soil pH (28 cases), including 15 unimodal and 10 positive relationships (Supplementary Table S9 and Supplementary Figure S8). Proportion of a tree host was the best predictor for 16 OTUs, all with positive effects. Of these, nine OTUs (mostly species of Amphinema, Tylospora, and Tricharina) were strongly associated with Picea abies, whereas proportion of Salix spp. predicted relative abundance of OTUs corresponding to Hyaloscypha finlandica, Tuber sp. and Cortinarius decipiens. Vegetation age was the strongest predictor for three OTUs, with negative effects on relative abundance of Cortinarius diasemospermus and two OTUs corresponding to Hyaloscypha finlandica species complex. Relative abundance of non-EcM OTUs was mainly determined by soil pH (36 cases), with positive (17) and unimodal relationships (12) prevailing (Supplementary Table S10 and Supplementary Figure S9). Nine OTUs responded strongest to relative abundances of tree genera, mostly Picea (3) and Pinus (2). Furthermore, two non-EcM OTUs displayed strong preferences for non-forest or park habitats. Among the dominant OTUs, there were almost no taxa that occurred in all habitats or in all types of forests indicating a certain level of niche differentiation. Congeneric OTUs responded similarly to soil pH but similarly (Amphinema spp.) or differently (Inocybe spp., Hyaloscypha spp.) to host plants, suggesting that abiotic niche is more conserved than biotic niche within genera.
To evaluate the validity of clustering, we assessed whether OTUs of Hyaloscypha finlandica-H. bicolor species complex have differentiated niches or represent potentially the same biological/ecological species with high marker sequence variation, by combining phylogenetic analyses, niche modeling and co-occurrence analysis (legend to Supplementary Figure S10). Phylogenetic analyses resolved OTUs belonging to the H. bicolor-H. finlandica complex poorly, except that OTU00823 branched at the base of the complex and OTUs corresponding to H. finlandica formed a monophyletic group within H. bicolor that remained paraphyletic. OTUs of H. bicolor featured similar preference to low soil pH, whereas OTUs of H. finlandica responded negatively to vegetation age but variably to soil pH and host plants such as Salix, Populus and Picea. OTUs of both groups displayed negative (and neutral in H. finlandica) within-group cooccurrence patterns. Surprisingly, several interspecific pairs of OTUs displayed strongly positive co-occurrence values (Supplementary Figure S10).
The niche analysis shows that OTUs of H. bicolor have similar ecological preferences to low pH, whereas OTUs of H. finlandica are typically pioneer species with broad pH tolerance and unexpected host preference. In spite of some ecological niche overlap, the C-score analysis revealed no positive co-occurrence (Supplementary Figure S10), suggesting that the closely related OTUs avoid each other, possibly by competition (H. bicolor) or both competition and niche differentiation (H. finlandica). The presence of some strongly positive co-occurrences among more distantly related OTUs suggests that cryptic species of H. bicolor and H. finlandica are not strongly competing with each other when they co-occur.
Taken together, the OTU-specific niche analysis provides evidence that several closely related OTUs matching the same species can respond differently to environmental predictors and hence be ecologically relevant. However, this may differ by taxonomic groups and our ability to define the niche in a multivariate space. We used 98% ITS sequence similarity threshold for OTU clustering, but the optimal, biologically most relevant thresholds depend on genera and particular species complexes. It is likely that these ecologically distinct, cooccurrence-avoiding sympatric OTUs correspond to biological species or groups of species, but biological species need to be determined based on mating tests or phylogenetic concordance of multiple genes in fungi (Taylor et al., 2006).

Habitat Type and Land Use
Our analyses of habitat type and land use effects (dataset#3; for definitions, see Supplementary Table S2) revealed that fungal richness and Shannon index were highest around ruins and in woodlands but distinctly lowest in bogs (Supplementary Table S11 and Figure 9). EcM fungi were most diverse in woodlands, followed by forests, parks and ruins, with lowest richness in croplands, grasslands, energy plantations and bogs. Conversely, AM fungi prevailed in grasslands, followed by croplands and wastelands, being almost absent from bogs. Root endophytes were, however, most diverse in bogs, followed by forests. Dung saprotrophs prevailed in grasslands, with distinctly lowest values in bogs and forests. Litter, wood and white-rot decomposers were most diverse around ruins. Unicellular fungi such as Aphelidiomycota, Monoblepharomycota, Neocallimastigomycota, Rozellomycota but also Zoopagomycota and unclassified fungi displayed greatest richness in energy plantations, whereas Spizellomycetes prevailed in croplands (Supplementary Table S11). Energy plantations were diversity hotspots also for Pleosporales and non-EcM Sebacinales. Croplands and grasslands were among the most preferred habitats for Zoopagomycota, Diversisporales, Glomerales, Pleosporales, Onygenales, Glomerellales, Hypocreales, Microascales, Sordariales, Auriculariales, and Ceratobasidiaceae. Ruins constituted richness hotspots for Helotiales, Thelebolales, Pezizales, Glomerellales, Xylariales and non-EcM Agaricales. Notably, diversity of only Hymenochaetales and Trechisporales peaked in forest habitats. Bogs represented the habitats of lowest richness for most fungal taxonomic groups, with croplands (Helotiales), energy plantations (Auriculariales, Trechisporales), wastelands (Monoblepharomycota, Hymenochaetales) and alleys (Zoopagomycota, non-EcM Sebacinales) also recorded as such for certain groups.
We attribute the differences in habitat types mostly to substantial differences in soil pH, δ 15 N and C/N ratio but also to tree community, particularly EcM vegetation (Supplementary Table S11). Woodlands and surroundings of ruins typically harbor a variety of microhabitats, which differ in pH and plant growth form, and high plant diversity that all favor high fungal richness. Energy plantations (P. x wettsteinii and Salix viminalis) were prominent for supporting high diversity of early diverging fungal lineages and pioneer groups of EcM fungi, which can be linked to high productivity and age effects, respectively. Similarly, the low richness of most fungal groups in bogs may be related to both low productivity and anaerobic conditions .
Wild, village and urban tree-dominated habitats differed in soil properties but also in tree composition, tree richness and overall EcM plant proportion (Supplementary Table S12). In particular, the village biome was associated with higher pH and δ 15 N as well as greater tree diversity than other biomes, which all have a positive effect on diversity of fungi and most groups therein. The relatively lower fungal diversity in forest habitats compared to village habitats is attributable to the low pH and low tree diversity of coniferous forests. Because the dummy variables of anthropogenic influence did not contribute to the best richness models, we conclude that the direct human impact on fungal diversity is generally low in natural and seminatural habitats in the anthropogenic gradient. Other authors have reported shifts in fungal community composition when including heavily polluted urban sites in this gradient (Timonen and Kauppinen, 2008;Jumpponen et al., 2010).

Island Ecosystems and Habitat Fragmentation
To test the effect of forest fragmentation including natural and anthropogenic islands on soil biodiversity, we excluded the following plots: non-forest habitat, vegetation age < 30 years, EcM plant basal area < 30%, selective harvesting < 5 years ago. This dataset#6 comprised 784 samples including 69 forest fragments and 107 islands (<20 km 2 ), which were categorized into aquatic islands (n = 37), bog islands (27) and fieldsurrounded islands (43). Island type and island size usually had a weak but significant effect on biodiversity of one third of the tested fungal groups. Contrary to our hypothesis, Shannon index and richness of fungi and most fungal groups were greater in island ecosystems than non-island ecosystems; in a majority  (V-X) animal parasites. In box plots, central lines, boxes, whiskers and circles represent medians, quartiles, 90% quantiles and outliers, respectively; letters above whiskers depict significantly different groups. In scatter plots, black and red lines respectively indicate linear regression lines when including and excluding non-fragmented habitats. Note that fragments larger than 1024 ha were not precisely measured and truncated to 1024 ha. Note the base-2 logarithmic scale in x-axis of scatter plots. of these cases, increasing island or forest fragment size had a negative effect on diversity (Supplementary Table S13 and Figure 11). For many groups, there was also a meaningful island type × patch size effect, indicating that size matters depending on the island type. Richness and Shannon index of fungi and most fungal groups were highest in field islands (64% of cases) but lowest in non-islands or bog islands (for example, AM fungi, various guilds of pathogens and saprotrophs, taxonomic groups of non-ectomycorrhizal fungi). Conversely, root endophytes, Hymenochaetales, Umbelopsidomycetes and Endogonomycetes displayed greater diversity in bog islands and non-islands compared with field islands.
Unlike most other fungal groups, EcM fungal richness did not differ among island types and non-islands, and there was no effect of island size or distance. In spite of this, EcM fungal lineages differed strongly in their preferences for specific island types and island size (Supplementary Table S13). While their preference for island type reflected their pH optima, their response to island size largely coincided with their inferred pioneer vs. competitor strategy. Groups with multiple pioneer species such as pezizalean lineages, /pisolithus-scleroderma and /tomentella-thelephora were relatively more abundant and diverse on smaller islands. Conversely, the /russula-lactarius, /piloderma, /clavulina, /cenococcum, /cantharellus, /amanita, and /amphinema-tylospora lineages were more common in large islands and non-islands (Supplementary Table S13).
Non-mycorrhizal fungi responded relatively more strongly to island size. Dung saprotrophs (R 2 = 0.083), Pleosporales (R 2 = 0.099), Glomerellales (R 2 = 0.081) and Sordariales (R 2 = 0.093) were most diverse on small islands. The only groups that preferred larger islands and non-islands were root endophytes, as well as Hymenochaetales, Umbelopsidomycetes and Endogonomycetes. In both EcM and non-EcM fungi, preference for field islands generally coincided with higher diversity in smaller islands, whereas preference for bog islands tended to reflect greater richness in larger islands and non-islands.
Excluding plots from non-island and non-fragmented habitats, island distance and distance × island type affected only a few fungal groups (Supplementary Table S13 and Figure 11). With increasing isolation, richness of total fungi, litter saprotrophs, wood saprotrophs, Capnodiales and Xylariales declined. The EcM fungal lineages /genea-humaria, /sebacina and /tomentella-thelephora were also less diverse in more isolated islands. In other guilds and taxonomic groups, reduction in diversity was statistically marginally significant or the effect was absent; no significant increases occurred.
The tendency for greater fungal richness in islands, particularly small islands, is in disagreement with previous work that demonstrated positive effect (Peay et al., 2010) or no effect (Li et al., 2020) of island size. The high diversity in small islands can be attributable to edge effect -a non-linear gradient of environmental conditions that extends up to 50 m in the habitat interior -which promotes overall habitat heterogeneity and niches (Crockatt, 2012). However, previous studies revealed decline in mycelium biomass, activity, reproductive effort and aboveground diversity toward the ecotone (Crockatt, 2012;Grilli et al., 2017). Our large plots cover various habitat conditions in the ecotone and interior of islands, increasing the number of sampled niches. We attribute the relatively low diversity of fungi in bog islands to low pH in these habitats, because the fungal groups that responded to bog island habitat coincided with groups most strongly driven by soil pH. Besides pH, island types and islands of varying size differed most strongly in soil C/P ratio, δ 15 N and relative basal area of Picea abies and Salix spp. (Supplementary Table S13) that are important predictors for many groups of fungi.
Fungal richness declined with increasing island distance, corroborating the results of Peay et al. (2010) in small tree islands separated by shrubs but not in a system of larger aquatic islands in China (Li et al., 2020). Of individual groups, representatives of both asexual microfungi and trufflelike and corticioid fruit-body formers had significant dispersal limitation, indicating no striking overall differences among dispersal strategies. Forest fragments displayed somewhat weaker island biogeography patterns compared with islands, but this is probably related to more recent history of isolation and short evolution as an island ecosystem (decades vs. centuries or millennia). Our results confirm that forest fragmentation reduces richness of certain poorly dispersed and old forest specialist groups (Grilli et al., 2017).

Old and Virgin Forests
Although vegetation age was an important predictor for several fungal groups, the effect usually leveled off at 80-140 years. Therefore, we tested if tree age per se and virgin conditions in particular affect diversity patterns in mature forest ecosystems. We separated forest plots aged > 100 years into virgin (n = 34) and non-virgin (n = 421) habitats (dataset#8) based on literature and expert evaluation. Virgin conditions and increasing age in >100-years old stands had weak but positive effects on richness of most fungal taxonomic and functional groups (including total fungal richness and Shannon index and EcM fungal richness), with or without accounting for other predictors (Supplementary Table S14). Both old-growth conditions (R 2 = 0.029) and age (R 2 = 0.016) promoted richness of white-rot decomposers. Similarly, richness of Capnodiales was greater in virgin forests (R 2 = 0.056), increasing with age (R 2 = 0.018). Virgin forests harbored greater richness of three major agaricomycete decomposers -Hymenochaetales (R 2 = 0.027), Polyporales (R 2 = 0.029) and Trechisporales (R 2 = 0.032), without a residual age effect. In EcM fungal lineages, we observed greater richness for /amphinema-tylospora (R 2 = 0.065), /cortinarius (R 2 = 0.026) and /piloderma (R 2 = 0.043) lineages in virgin forests. Considering habitat conditions, virgin forests exhibited greater amounts of topsoil C (R 2 = 0.047) and K (R 2 = 0.059) and proportion of Picea abies (R 2 = 0.043) but lower δ 15 N values (R 2 = 0.049). With age increase from 100 to 300 years, there was a slight increase in topsoil C (R 2 = 0.023) and K (R 2 = 0.023) concentrations and proportion of Quercus robur (R 2 = 0.063; Supplementary Table S14).
Virgin forests harbor large amounts of decaying wood in various stages of decomposition and develop gaps that are suitable for establishment of pioneer species (Harmon et al., 1986). Richness increase in white-rot decomposers and corresponding Hymenochaetales, Trechisporales and Polyporales (partly brown-rotters) orders probably benefit from greater amounts of woody residues that are buried in soil. Many wood-inhabiting species of Hymenochaetales and Polyporales are red-listed and considered indicators of primeval habitats (Lõhmus et al., 2018). As described above, diversity of Hymenochaetales and Trechisporales declined in island habitats, suggesting that a combination of limited dispersal capacity and preference for virgin forests renders these taxa particularly vulnerable to habitat loss and fragmentation.
In selectively harvested plots, the decline in overall fungal richness was related to decrease in EcM fungi and yeasts. Our results are thus partly consistent with a study in Northern Sweden that found 30-40% decline in EcM fungal diversity in a few years after partial harvesting (Sterkenburg et al., 2019) but in a disagreement with studies in Spanish Mediterranean forests revealing no effect (Castano et al., 2018a;Parlade et al., 2019). However, our study corroborates the long-term selective harvesting effect on composition of fungi and specifically EcM fungi in Mediterranean forests, which was not evident in Sweden. Compared to 30-40% basal area harvesting in these studies 1-5 years ago, our plots were harvested on average at 23.8% intensity 10.2 years previously. These harvesting effects can be attributed to three main factors -(i) mechanical disturbance and disruption of EcM fungal networks; (ii) increase in fresh stumps and dying root systems; and (iii) short-term rise in soil pH -that all favor proliferation of opportunistic groups at the expense of EcM fungi (Lindahl et al., 2010). The observed shifts in EcM fungal composition may be related to altered competitive balance due to losses of the carbohydrate source and massive disturbance that favor establishment of species from genera with pioneer strategy such as Hebeloma, Paxillus, and Tuber. Relative abundance of contact exploration type, mostly representing Russula and Tuber spp., was positively related to harvesting, which points to an advantage of low C demanding groups after declining carbohydrate sources (Leake et al., 2004).
With increasing time since selective harvesting, shade-tolerant pioneer species such as Sorbus aucuparia and Corylus avellana and many other small bushes (Supplementary Table S15) establish from seed or stump sprouts, potentially contributing to recovery of fungal richness besides the effect of aging per se. However, full recovery of selectively harvested plots may take several decades or centuries if long-lived EcM tree species are replaced by long-lived AM trees (Chen J. et al., 2019). Even without major vegetation shifts, recovery of EcM fungi from a clear-cut state may take 90 years based on model averages (Spake et al., 2015). Correspondingly, our results from this and the virgin habitat analysis extend these findings, indicating that selective harvesting and thinning may affect fungal composition for multiple decades.

Management of Parks and Woodlands
To maintain the semiopen structure, European woodlands and parks require consistent management by grazing or mowing and removal of coppicing sprouts. We tested whether management of seminatural woodlands and parks affects diversity of soil organisms (dataset#10: 40 wooded meadows, 17 wooded pastures, 100 parks). Wooded meadows showed higher fungal richness (R 2 = 0.113) and Shannon index (R 2 = 0.061) than pastures ( Table 3). This pattern also occurred in soil saprotrophs (R 2 = 0.101), especially Mortierellomycetes (R 2 = 0.166) and Umbelopsidomycetes (R 2 = 0.080) and Rhizophlyctidomycetes (R 2 = 0.141), and non-significantly so for most fungal taxonomic and functional groups. Of EcM fungal lineages, /cortinarius (R 2 = 0.213) and /russula-lactarius (R 2 = 0.050) had substantially greater richness in wooded meadows than pastures, but this pattern was non-significant for the overall EcM fungal and lineage richness.
Greater woody plant richness, establishment of conifers along with several other tree species and lower mechanical disturbance probably favor richness of most fungal groups in coppiced habitats. Because of contrasting management effects on different taxonomic and functional groups, patch-wise management and restoration of parks and woodlands probably promotes greatest fungal richness in these ecosystems. Retention of as many tree species as possible will probably benefit fungal biodiversity most.

The Prevailing pH Effect
One of the most remarkable findings of this study is the overwhelming effect of soil pH or its closely related variables (C/N ratio, δ 15 N, cation concentration) on the distribution of individual fungal OTUs, which collectively drive the observed patterns in diversity and composition of most guilds and higher taxa of fungi. While most groups had a distinct unimodal relationship with soil pH, many order-level taxa displayed preference for extreme pH values, suggesting phylogenetically conserved adaptation to such conditions (Rousk et al., 2010a) or the realized bimodal niche in yeasts, for example. The effect of pH is related to physiological capacity of organisms to tolerate high H + or OH − ions and these differ strongly across broad taxonomic groups (Rousk et al., 2010b). Overall, the strong effect of pH is characteristic to much of the soil microbiome, especially bacteria (Rousk et al., 2010a;Bahram et al., 2018). The average pH effect on non-EcM fungal OTUs (R 2 average = 0.195), EcM fungal OTUs (R 2 average = 0.081) and orders (R 2 average = 0.120) somewhat exceeds the regional-scale pH effect on bacterial genera (R 2 average = 0.144) and phyla (R 2 average = 0.072) across France (pH range 3.7-9.2; Karimi et al., 2018). Our study therefore extends the high impact of soil pH on diversity of fungi (Rousk et al., 2010a;Glassman et al., 2017) and its functional and taxonomic groups by removing the confounding effects of latitude, climate and historical effects that prevail in global-scale studies (Tedersoo et al., 2014;Egidi et al., 2019;Vetrovsky et al., 2019;Oliverio et al., 2020). Based on extrapolated soil chemical data, Vetrovsky et al. (2019) found a relatively low soil pH effect on fungal diversity, indicating the importance of using empirically determined values. Furthermore, our correlation analyses and SEMs suggest that soil pH affects fungal composition also indirectly via plant diversity and composition as well as proliferation of molds (Figure 7 and Supplementary Figure S11). EcM fungi and evergreen plants acidify soil by exuding organic acids and shedding recalcitrant litter (Cornelissen et al., 2001;Tedersoo and Bahram, 2019). Therefore, combination of certain tree species and bedrock type may result in extreme pH values that magnify the pH effect on soil microbiome. The extrapolated maps of soil mycobiome mostly follow habitat type and soil pH (Figure 2 and Supplementary Figures S12-S18). Soil pH also determines the solubility of minerals containing P and cations, and hence the availability of these nutrients to both microorganisms and plants. Apart from Ca and Mg concentration, we found only negligible effects of soil nutrient concentration on fungal richness, but the effects of interaction between nutrient concentration, δ 15 N and soil pH are difficult to test non-experimentally. Although soil Mg availability is important for plant chlorophyll production, stress resistance and intercellular transportation (Senbayram et al., 2016), little is known about its role in soil microbiome, providing no clear explanation for the prominent effect of soil Mg on richness of certain fungal taxa and ciliates (L. Tedersoo et al., unpublished).

Tree Species Effects
Besides the indirect effects of soil pH, tree species and genera had strong additional effects on fungal diversity and composition, which sometimes prevailed in niche analysis of individual OTUs, particularly EcM taxa. Consistent with our hypothesis that biotrophic groups respond more strongly to host plants than saprotrophs, tree species had 27 and 52% stronger effects on composition of plant pathogens and EcM fungi compared with saprotrophs. In terms of richness, plant species did not affect the diversity of soil and dung saprotrophs or mycoparasites and animal parasites, but had a weak effect on litter (3.1% variation) and wood (4.4% variation) saprotrophs and leaf pathogens (4.4% variation) -guilds that are associated with decomposing or living plant material. Yet, individual host trees collectively explained 35.0% variation in EcM fungal diversity. Any specific host effects on AM fungi were probably masked by the dominant effect of EcM plant proportion and our focus on woody plants exclusively. Although AM fungi lack host specificity, many AM plants associate with greater diversity of fungi than others and display differences in community composition (Põlme et al., 2018;Sepp et al., 2019). Values in bold are significant (P < 0.001). 1 ns, not suggestive at P > 0.05; 2 WM, wooded meadow; WP, wooded pasture; > / = / < indicate significant differences among groups; 3 C, coppiced; M, managed; UM, unmanaged; > / = / < indicate significant differences among groups.
The importance of individual tree genera on fungal composition was unrelated to the frequency of these taxa and differed remarkably across both plant and fungal groups. Pinus sylvestris was always associated with the most divergent composition stretched out along the primary axis in NMDS plots, with Picea abies also in a separate position (Figure 8 and Supplementary Figure S7). The high diversity of several unicellular fungal groups and protists (L. Tedersoo et al., unpublished) in hybrid poplar energy plantations is unexpected, but it may be related to high primary productivity and perhaps well-developed herbaceous understorey. The effects of Betula spp. and Populus tremula on any fungal functional group were highly similar, although both host trees associate with a number of specific EcM, saprotrophic and pathogenic fungi. Perhaps due to the mull soils, Corylus avellana and Tilia spp. also shared much of the fungal composition and plots dominated by AM angiosperms (Acer, Ulmus and Fraxinus species) overlapped in fungal composition. However, the AM Juniperus was associated with highly distinct saprotroph communities, which was not observed for plant pathogens. Based on PERMANOVA analysis, other tree taxa harbored distinct EcM fungal communities, but similar communities of saprotrophs and pathogens. Although several studies have indicated strong grouping of closely related plant species in fungal composition (Põlme et al., 2013;Tedersoo et al., 2013;Wu et al., 2018;Wang et al., 2019;Yang et al., 2019), the lack of grouping of closely related tree genera reflects the overall low importance of tree phylogeny above genus level.
Taken together, most overstory tree species affect community composition of EcM fungi, but only a few have significant influence on other fungal guilds. At the plot scale, greater overall tree diversity promotes fungal richness either directly or indirectly through positive and unimodal effects of individual tree taxa. In particular, the EcM conifers Picea abies and Pinus sylvestris are associated with strongest effects on diversity and composition, with part of these effects attributable to soil acidification. Therefore, both of these conifers can be regarded as non-redundant keystone species in the temperate and boreal forest belt of North and Central Europe by generating unique habitats suitable for the specialist host-associated and acidophilic fungal communities. Although pure stands of conifers were associated with reduced diversity of several taxonomic and functional groups, their presence adds much to the landscapescale soil microbial diversity.

Methodological Implications
The universal eukaryotic primers enabled us to capture the full diversity of kingdom Fungi and other eukaryotes as well, but the relative sequence abundance and identification power for protists and animals was much lower because of high fungal biomass in soil and poorly managed reference data of non-fungal ITS sequences. For these groups, several-fold deeper sequencing or use of group-specific primers and rigorous development of reference databases would certainly benefit species-level biodiversity analyses using environmental DNA.
The 18S-ITS metabarcoding by third-generation sequencing allowed us to use the entire variable ITS region for identification, which represents a significant improvement compared to short ITS1 or ITS2 barcodes in terms of taxonomic resolution (Garnica et al., 2016) and precision of identification (Tedersoo et al., 2018b), reducing the proportion of unidentified taxa by an order of magnitude compared to studies using short amplicons Chen L. et al., 2019;Sterkenburg et al., 2019). The vast majority of the well-delimited, abundant fungal species were represented by a single OTU, which contrasts to the occurrence of multiple co-occurring sister OTUs as revealed by short-read metabarcoding or metagenomics using 454 or Illumina platforms . The long barcode enables more accurate clustering, because random PCR and sequencing errors are less likely to accumulate over longer marker regions and sequences with a few errors are more easily incorporated into the main OTUs (Tedersoo et al., 2018b). Furthermore, the PacBio platform provides better quantitative output compared with the Illumina platform (Castano et al., 2020).
Our study represents perhaps the only oversampled microbiome investigation on a regional geographic scale, which is reflected by virtually stabilizing OTU accumulation curves (Figure 4). The high sample coverage enabled us to test for meager effects of variables (R 2 = 0.006 and R 2 = 0.002 correspond to P < 0.001 in GLM and PERMANOVA analyses, respectively), far beyond their potential biological importance. However, the leveling off rarefaction curves need to be interpreted in the context of sequence data filtering and similarity thresholds. By incorporating sequences into UNITE Species Hypotheses (Kõljalg et al., 2013), we observed that many global singletons, otherwise excluded from analyses, fell into existing UNITE SHs, representing truly rare species in the context of this study. For example, singletons included the only representative of the/acephala macrosclerotiorum EcM lineage and several species of Tuber and Helvella. However, most singletons were chimeric, of low-quality, or possessed > 5mer indels rarely observed in Sanger sequences. Based on the data set of Schoch et al. (2012), it is also clear that multiple biological species of fungi cannot be distinguished based on the 98% or even 100% ITS sequence similarity threshold, for example many groups in Penicillium and Fusarium (Seifert, 2009). Based on ITS metabarcoding alone, we cannot predict the diversity of these species-rich genera and the proportion of taxa that are missed or lumped with this approach. For the morphological species of Hyaloscypha with around 3% ITS sequence variation, we showed that OTUs at 98% similarity level are ecologically differentiated, but whether this is an optimal threshold or lumps together multiple other ecologically distinct OTUs, remains to be evaluated more broadly.
With respect to statistical analyses, PERMANOVA and GDM algorithms strongly complemented each other in multivariate analyses by recovering the major predictors and determining their value-dependent importance (Landesman et al., 2014;Glassman et al., 2017), respectively. Combining a machine learning algorithm for initial variable selection with quadratic GLMs proved useful for determining best models amongst hundreds of variables with linear or non-linear effects within minutes. SEMs, however, did not add as much information as in global-scale studies because of the exclusion of latitude and climatic predictors (as non-significant), which have strong direct effects and indirect effects via floristic and edaphic variables on soil microbiome over larger geographic scales (Tedersoo et al., 2014;Bahram et al., 2018). As a downside, our results can be directly extrapolated to a relatively narrow climate belt in North Europe, with more caution required when interpreting to Central European, North American or Asian conditions. Although comparable climatic regions in these continents harbor similar soil conditions and share many tree genera, many other genera differ and may have strong influence on soil biota similarly to P. abies and P. sylvestris in our study (Wu et al., 2018). However, we believe that many of the general biodiversity patterns recovered in our study apply to other temperate habitats with comparable gradients in soil properties.
The main technical issues emerging from this study are the temporal sampling effect, investigator bias and sample spoilage by molds. While these could be modeled and accounted for in our dataset, this is far more complicated for smaller datasets and especially global datasets, where samples from different regions are typically collected by distinct persons on a specific time point, generating a set of confounding geographic, temporal and investigator effects that cannot be disentangled. In our samples, both the temporal effect and investigator bias were mostly related to the relative abundance of molds, suggesting that molds may become problematic after certain weather events such as dryingrewetting or freeze-thaw (Blazewicz et al., 2014;Castano et al., 2018a;Burke et al., 2019) and depend on individual-specific care during sample pre-treatment. Therefore, we recommend that the abundance of molds, spatial and temporal sampling effects and investigator biases be at least considered in diversity models whenever such information is available.

CONCLUSION
Our multi-year, regional-scale investigation demonstrates that taxonomic and functional distribution of much of the fungal kingdom follows soil pH more than previously expected. The pH effect also prevailed in comparisons of habitat types and anthropogenic impact. SEMs and correlation analyses indicated that soil pH is related to soil Ca concentration and the abundance of most tree species, which also affected both richness and composition of fungal groups. Accounting for pH, woody plant species richness has an independent positive effect on overall fungal diversity, which integrates differential biotic and abiotic tree species effects on different fungal guilds and taxonomic groups. Partly through altering soil pH, the conifers Picea abies and Pinus sylvestris serve as important tree species with strongest influence on diversity and composition of most fungal groups in North European forest ecosystems, which renders these tree species as keystone components in the soil microbiome perspective.
The analyses revealed mixed support to our hypotheses. In line with the first hypothesis, biotrophic groups responded more strongly to dominant vegetation than free-living groups, but both were generally more affected by soil pH. In contrast to second hypothesis, diversity of soil fungi tended to be higher in island habitats due to the edge effect. However, as predicted, fungal richness declined with island distance and in response to forest fragmentation. In agreement with the third hypothesis, pristine habitats supported somewhat higher fungal diversity, but there were no differences between natural and anthropogenic habitats focusing mostly on parks and coppiced gardens. In support to the fourth hypothesis, diversity of most fungal groups suffered from management of seminatural woodlands and parks (coppice cutting and mowing) and forests (selective harvesting or thinning of trees), but especially for forests the results depended on fungal group and time since partial harvesting. With no effect of microclimatic variables, the North European soil mycobiome is not vulnerable to direct effects of climate change. However, increasing drought periods and pest or pathogen outbreaks may alter tree composition with further effects on soil microbes. Substantial differences in richness of fungi and several functional groups across habitat types and in response to management suggest that shifts in land use may greatly affect fungal diversity and functional guild composition.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
LT designed the study. LT, SA, RD, IH, RL, RR, ER, KA, TD, HT, KJ, IS, EO, SP, MM, KL, and AA collected the samples. SA, MB, KP, FB, AP, NH-D, VM, DG, RA, RP, IV, UK, and KA performed the analyses. LT wrote the manuscript with input from other authors. All authors contributed to the article and approved the submitted version.