Estimation of Fungal Diversity and Identification of Major Abiotic Drivers Influencing Fungal Richness and Communities in Northern Temperate and Boreal Quebec Forests

Fungi play important roles in forest ecosystems and understanding fungal diversity is crucial to address essential questions about species conservation and ecosystems management. Changes in fungal diversity can have severe impacts on ecosystem functionality. Unfortunately, little is known about fungal diversity in northern temperate and boreal forests, and we have yet to understand how abiotic variables shape fungal richness and composition. Our objectives were to make an overview of the fungal richness and the community composition in the region and identify their major abiotic drivers. We sampled 262 stands across the northern temperate and boreal Quebec forest located in the region of Abitibi-Témiscamingue, Mauricie, and Haute-Mauricie. At each site, we characterized fungal composition using Illumina sequencing, as well as several potential abiotic drivers (e.g., humus thickness, soil pH, vegetation cover, etc.). We tested effects of abiotic drivers on species richness using generalized linear models, while difference in fungal composition between stands was analyzed with permutational multivariate analysis of variance and beta-diversity partitioning analyses. Fungi from the order Agaricales, Helotiales, and Russulales were the most frequent and sites from the north of Abitibi-Témiscamingue showed the highest OTUs (Operational Taxonomic Unit) richness. Stand age and moss cover were the best predictors of fungal richness. On the other hand, the strongest drivers of fungal community structure were soil pH, average cumulative precipitation, and stand age, although much of community variance was left unexplained in our models. Overall, our regional metacommunity was characterized by high turnover rate, even when rare OTUs were removed. This may indicate strong environmental filtering by several unmeasured abiotic filters, or stronger than expected dispersal limitations in soil fungal communities. Our results show how difficult it can be to predict fungal community assembly even with high replication and efforts to include several biologically relevant explanatory variables.


Introduction
Fungi represents the most diverse groups on Earth after insects [1]. They are diverse in terms of shape, color, and lifestyle and they are distributed throughout the globe in all types of terrestrial and aquatic ecosystems, and even in extreme life conditions such as those prevailing in the Antarctic [2,3]. They play pivotal roles in ecosystem functioning by modulating nutrient cycling, organic matter decomposition, carbon storage, and plant nutrition through the formation of mutualistic symbiosis [4][5][6][7]. Several studies have shown that ecosystems functionality is altered when changes occur in soil fungal diversity [5,6] which is why they are used as a forest and soil health bioindicator [8,9]. In addition some species, especially edibles, are an important source of income and a significant component of the diet, particularly in developing countries and remote regions [10]. Unfortunately, fungal diversity can change because of time, climate, biota, topography, natural disturbance, or human caused perturbation and contamination [11][12][13][14]. For these reasons, there is interest in developing approaches to predict various facets of fungal diversity, and how it is likely to change over space and time in natural and managed ecosystems.
Fungal species richness, an intuitive component of fungal diversity, is often used to compare habitats, as species diversity is often assumed to reflect, at least partly, niche diversity when limiting similarity drive species coexistence [15]. Fungal richness has also been linked to key ecosystem parameters and functions such as litter decomposition rates [16] or soil and plant nutrient content [7]. Many abiotic variables can influence species richness. In the case of fungi, climatic variables such as mean annual precipitation and temperature are important drivers of fungal richness at global scale [17], but other variables such as vegetation cover [18], tree species [19], and soil texture [20] can have positive or negative impacts on richness.
Beyond the simple counting of species number (i.e., α-diversity), another way to characterize biodiversity is to study the composition of the community and how it changes across sites (i.e., β-diversity). Understanding compositional pattern of species helps researchers to understand different aspects of species interaction and ecosystem function [21,22]. β-diversity patterns can also guide conservation practices, for example by providing knowledge about the uniqueness of community composition in the landcapse [23]. Community dissimilarity can stem from two mutually exclusive, and thus additive, components: richness difference (i.e., nestedness) and compositional turnover/replacement [21]. Because these components are additive fractions of community β-diversity, they can easily be mapped in a ternary plot linking nestedness, turnover, and community similarity (i.e., 1-β) [22]. Such β-diversity partitioning analyses can help formulate more specific hypotheses regarding the potential drivers of community dissimilarity in a metacommunity [24].
Nevertheless, in Quebec, Canada, fungal diversity in northern temperate and boreal forests is generally poorly described. A reliable and exhaustive global description of the fungal diversity is needed for the protection and conservation of fungal species, and to further predict changes in fungal richness and community composition in our forests. This lack of data represents an obstacle for mycologist and for companies working with natural resources, and data is needed to improve forest fungal biodiversity conservation and management. The main aim of this study was to sample the northern temperate and boreal Quebec forests; to (1) make an overview of the fungal richness and the community composition; and (2) identify the main environmental drivers of these facets of fungal biodiversity.

Study Sites and Soil Sampling
The sampling design consisted of 262 sites that were located in three different regions in the province of Quebec (Canada). The sites were divided in 130 stands located in Abitibi-Témiscamingue (AT), 50 stands in Haute-Mauricie (HM), and 82 stands in Mauricie (MAU) (Figure 1). We used a forest inventory database composed of 75 abiotic variables and 33,000 stands, released in 2009, through the Ministry of Forest and Parks of the province of Quebec. This database served to identify and locate sites according to four ecological parameters used by stakeholders: the stand age, the dominant tree species, the level of drainage, and the surficial deposit. We selected sites to obtain the largest diversity of ecological parameters combinations that characterized each region. Sites also had to be accessible by road. Sites were sampled between June and September 2017. For each site, we used ropes to delimitate a representative area of 400 m 2 (20 × 20m). In each site, two 20 m transects were established 10 m apart from one another. The forest litter was removed from the surface and 6 soil cores of 20 cm depth were collected every 5m along the transects. Soil samples were pooled, mixed and a subsample was taken for DNA analyses and stored immediately on ice and at −20 °C at the end of the day. The rest of the soil samples were sieved and stored at different temperature depending of the future analyses (see below).

Clustering of Sites into Climatic Domains
Using georeferenced information of the 249 stands (positive at sequencing step), we extracted historical climatic data ranging from 1974 to 2008, from the site of the government of Canada (https://ouvert.canada.ca) (Table S1). Total degree-days (based on 0, 5, and 10 °C), mean length of the season without frost (T min ≤ 0, −1, −2, −3, and −4 °C), mean length of the growing season, average cumulative precipitation during the growing season, and the mean temperature during the growing season (calculated with the total degree-days based on 0 divided by the length of the growing season) were used for K-means clustering with kmeans function and cluster package in R software v.3.5.3 [25,26]. This resulted in three climatic domains distributed as follows: 70 sites, mostly localized at the north of Abitibi-Témiscamingue (exception of four sites in Haute-Mauricie and two in south of Abitibi), were grouped in domain 1 (NAT), domain 2 (MAU) contained 79 sites situated in Mauricie and domain 3 (SAT-HM) included 100 sites from the south of Abitibi-Témiscamingue (exception of five sites at the north) and Haute-Mauricie ( Figure S1).

Soil Analyses
Part of the soil sample from each site was dried at room temperature to measure pH (hydrogen potential) and Eh (redox potential) with the method proposed by Husson et al. [31]. Measurements were made using platinum M241Pt electrodes and a 10Mhoms resistance multimeter inside a Faraday cage. The measure was considered stable after 30 s. Because Eh and pH are strongly correlated values, Eh was standardized for a pH of 7 based on Equation (1) where the Eh is measure in volts, R the perfect gas constant (8.314471 JK −1 mol −1 ), F the Faraday constant (96485.3383 C mol −1 ), and T the temperature (Kelvin): Temperature and pH were taken with a pH electrode (HI 1292D HANNA) immediately after the Eh. Organic matter (OM) was measured by mass loss ignition (550 °C during 5 h). After heating, the soil was analyzed for texture using a Fritsch Analysette 22 microtech plus to get the content (%) of silt, sand, and clay.

Soil DNA Extraction, Amplification, and Sequencing
Soil for DNA was kept frozen in −20 °C until extraction and then DNA from 0.30 g of soil was extracted with the DNeasy PowerSoil kit (Qiagen, Montreal, QC, Canada) according to the manufacturer's instructions. 200 µL of 100 µM ammonium aluminium sulfate (AlNH4(SO4)2) was added, as an additional step, in the column to remove inhibitors such as humic acid and increase PCR efficiency on several types of soil (clayey, sandy) [32]. DNA was eluted in 100 µl. PCR reactions were performed in a volume of 25 µL containing 2.5 µL of 10 × standard Taq buffer (100 mM Tris-HCl, 500 mM KCl, 15 mM MgCl2), 0.5 µL of dNTPs (10 mM), 0.5 µL of primers (20 mM), 0.25 µL of standard Taq DNA polymerase, 2 µL of total soil DNA and 18.75 µL of nuclease free water. We amplified the ITS2 region with the universal fungal primers fITS7 and ITS4 with a CS1 (ACACTGACGACATGGTTCTACA) and CS2 (TACGGTAGCAGAGACTTGGTCT) tag added at the 5′ extremities [33,34]. Three different melting temperatures were used for each PCR reaction with the following thermocycling conditions: 34 cycles of 30 s at 94 °C, 30 s at three different melting temperatures of 48, 51, and 58°C, and 1 min 30 s at 72 °C. Reactions were preceded by a 3 min denaturation step at 94 °C and terminated with a 5 min elongation step at 72 °C. If no PCR band was obtained, a 1/10 dilution of the total soil DNA was prepared. Positive PCR from the same sites (3 PCR × 256 sites) were pooled and sent for sequencing. A second PCR was done to add adapters and barcodes (one barcode per site), libraries were quantified and the sequencing was performed with Illumina MiSeq PE 300 following the procedures of Quebec Genome Innovation Center of McGill University (Montreal, QC, Canada).

Bioinformatic Analysis
Raw sequences were demultiplexed and adapters and barcodes were removed for 249 sites (7 samples didn't work at sequencing). The sequences were uploaded to the Galaxy web platform v.1.1.2-806.1 and analyzed for quality using the FastQC Read Quality reports v.0.72 [35,36]. Then sequences were filtered and trimmed (trimleft = c (19,20), truncLen = c (230, 200), maxN = 0, maxEE = c (2, 2), minLen = 150) using dada2 package v. 1.7.9 in R [26,37]. Only sequences longer than 150 bp with a mean number of expected errors below 2 were kept, trimming these to keep only the first 230 bp for the forward reads and the 200 bp for the reverse reads. Paired-ends sequences were merged using Fastq-join (pN = 0, mN = 20) in the Galaxy web platform [38]. Then Fastq files were converted in Fasta files with FASTX-toolkit v.1.0.2 and ID sequences were changed for the unique number ID of the sites using Seqkit v 0.8.0 [39,40]. The open source software ITSx v. 1.1.1 was used to remove conserved regions (5.8S and 28S) may be present in our ITS2 amplicon [41]. Finally, dereplication procedure for each site was independently (sample level) applied using Vsearch v. 2.8.0 and then all files (249 files corresponding to 249 sites) were combined in one single Fasta file in order to obtain one OTU data file and compare fungal composition between all sites [42]. Then, the dereplication procedure was repeated, but across all sequences. We removed singletons (minuniquesize = 2) chimera sequences using the reference database (UNITE/UCHIME reference datasets v.7.2) and de novo chimeras, and we clustered sequences (id = 0.97) and performed taxonomic assignation (id = 0.98) with our local barcode database containing 169 barcodes from fungi collected in Canada (mostly in Quebec province) (Document file S1) and public reference databases (USEARCH/UTAX reference datasets v.7.2 from UNITE and 600,000 ITS2 sequences download from the National Center for Biotechnology Information Search database (https://www.ncbi.nlm.nih.gov/)) using Vsearch v. 2.8.0. Non fungal OTUs were removed.

Statistical Analysis
We obtained two databases, one with the abiotic variables for each site and the one with presence/absence OTUs by site. All statistical analyses were coded in R [26]. Taxa frequency bar chart was created using ggplot2 [43]. To identify major abiotic drivers, we correlated fungal occurrences in 240 sites (we excluded 9 sites with missing environmental variables) with 18 abiotic variables from these sites (Table S1). Percentage of sand was not used because of high correlation with clay and silt content. Mean temperature during the growing season and average cumulative precipitation during the growing season were used as climatic variable for each sites. Other climatic data were excluded from the statistical analyses because of high colinearity. Continuous numerical variables were standardized prior to analysis. Generalize linear model with negative binomial family error distribution was used to identify the drivers of fungal OTU richness, using the MASS and MuMIn R packages for model selection based on dredge and importance functions [44,45]. A pairwise site dissimilarity matrix based on Jaccard's index was computed for 240 sites keeping only 14 000 OTUs appearing in these 240 sites. A PERMANOVA analysis using the adonis2 function (vegan R package) was used to explain our fungal distance/dissimilarity matrix using our set of 18 abiotic variables [46]. Prior to running the PERMANOVA analysis, we verified the assumption of multivariate dispersion using the function betadisper (vegan R package) [46]. To verify the homogeneity of variance for the continuous predictors, we quantified the distance of each sample in the ordination to the overall metacommunity-wide centroid, and correlated this value to the continuous variable. Principal coordinate analysis (PCoA) was used to visualize fungal community composition among sites. Beta-diversity partitioning analyses were used to evaluate whether shifts in fungal community structure were mostly due to species replacement or differences in fungal richness among sites using the beta.div.comp function (coef = "J") from Legendre (2014) [21]. This beta-diversity partitioning could be visualized on a ternary plot using the R package Ternary [47]. Species indicator analysis was used to determine whether some species were preferentially associated with specific site characteristics [48].

Taxonomic Diversity
Illumina sequencing of our fungal amplicon resulted in 12,455,757 raw sequences which was pruned to 264,000 unique reads after bioinformatics treatment. Clustering at 97% yielded 14,755 OTUs (6687 singletons), out of which 4796 could be assigned to a known taxon. We merged OTUs (and their frequencies) if they received the same taxonomic assignation and removed OTUs that were assigned to non-fungal organism for a total of 252 removed OTUs. For the 249 sites, 14,503 OTUs, including 4544 OTUs with taxonomic assignation, were obtained.

Fungal Richness
OTUs richness varied between 26 to 825 OTUs across the sites with a mean of 193 OTUs per site. Lower OTUs richness was observed in domain SAT-HM compared with domain NAT and MAU ( Figure 3). We could also observe a tendency of OTUs richness to decrease in presence of higher moss cover and in older stand age but, to increase with high herb cover ( Figure 4A-C). OTUs richness also varied depending on the dominant tree species. Intolerant hardwood (IH) stands showed high fungal richness as opposed to jack pine stands (JP) ( Figure 4D). A generalized linear model was used first, with all the habitat variables, and used in dredge to performed model selection and model averaging. The dredge function generated a set of 262,144 possible combinations (Table S2). Model at the first rank (lowest AICc) included the variables stand age, average precipitation, herb cover, moss cover, soil pH and silt content, but all first six models were equivalent in weight and for Akaike's information criteria (AICc). Variables that were not included in all the first six models were deposit, drainage, EhpH7, litter, organic matter, and shrub cover. These variables also have low sum of Akaike weights, meaning they appeared in few models ( Figure 5). The moss cover and the stand age stood out as the most important predictors of fungal OTUs richness with a sum of Akaike weights close to one (appeared in almost all the models) ( Figure 5). Herb cover, average precipitation, silt content and lichen cover could also important predictors of fungal richness with a sum of Akaike weights over 0.5, indicating that they appeared in half of the possible models ( Figure 5).

Fungal Community Composition
Soil pH has a significant effect on fungal composition, followed by the average cumulative precipitation and stand age (Table 1 and Figure 6). For stand age, older stand (blue) seemed slightly nested into those from younger stand (red) ( Figure 6C). The multivariate dispersion analysis indicated that the assumption of variance homogeneity was violated for stand age (F = 3.68, p = 0.06) (Document file S2). However, given the robustness of the PERMANOVA test to this assumption, this should not translate as a major bias in our downstream analyses [49]. For the continuous variables, no statistically significant correlations with their distance to the centroid were found for soil pH and average precipitation, the two best predictors (Document file S2 and Figure S2).
We analyzed the fungal community structure with a simplex analysis of beta-diversity partitioning from pairwise comparisons of stands by using different indices that measure similarity (number of share species), relativized richness difference and relativized species replacement according to the R-script methods of Legendre et al. (2014) [21]. Values obtained were positioned into an equilateral triangle, the so-called ternary plot. High species turnover, perfect nestedness, and richness identity could be easily recognized because the points are positioned on a given side of the triangle. This analysis shows an alignment of the majority of the pairwise site comparison at the bottom of the ternary plot, which evidences species replacement as the major driver of community dissimilarity in this regional metacommunity ( Figure 7A). This could be due to a very high frequency of rare taxa, which could have inflated the turnover component of beta-diversity. Yet, a similar pattern of strong turnover was observed even when only the most frequent fungi (appearing in more than 30 sites) were retained for the analysis ( Figure 7B). Likewise, the pattern held constant when considering only the rarer fungal taxa ( Figure 7C).

Fungal Diversity in Temperate and Boreal Quebec Forest
Diversity in the regional species pool could be explained by the large environmental heterogeneity across our sampled sites and the sequencing depth used. Indeed, our study gathered data on highly contrasted sites differing in dominant tree species, geological substrate origin, drainage, climate, and abiotic soil parameters (Table S1). Many studies and reviews have supported the role of environmental heterogeneity in supporting greater species diversity through niche partitioning [50][51][52]. A study published by Tedersoo et al. (2014) also reported very large fungal richness; they sampled 365 sites across the world divided in eleven different biomes (not limited to boreal and temperate forest) and detected 44,563 nonsingleton fungal OTUs that corresponded to approximately half of the described fungal species on the globe [17]. They also observed that boreal and temperate forests were also more diverse in ectomycorrhizal and several other fungal guilds than other biomes. The three most frequent orders that characterized our regions were the Agaricales, Helotiales, and Russulales ( Figure 2). These three groups are generally dominant in boreal forest [11,[53][54][55]. The genus Lactarius and Russula, from the order Russulales, are abundant and diverse ectomycorrhizal fungi in northern forests, especially because there are host generalists [14,56,57]. In our data, the genus Russula and Lactarius represented 84% of the OTUs from the Russulales (Table  S3). Helotiales are common and frequently associated with roots in boreal forest [11,54], whereas the Agaricales are known to associate with bryophytes a dominant component of the boreal forest understory [58]. Agaricales and Russulales are also known to be abundant and diverse in temperate forests [59][60][61][62][63], which represented 32% of sites included in our study. In temperate forests, Agaricales are mostly represented by the genus Cortinarius and Inocybe [60,63]. Cortinarius genus represented 24% and Inocybe 14% of our OTUs belonging to the Agaricales (Table S3).
Soil was collected from June to September, which included a sampling done over two seasons (summer and beginning of fall). Fungal groups do not seem to change between season, although their abundance may change [64,65]. Active mycelium can decrease at the soil surface in the middle of summer, because of drought, but remains detectable [66]. Season should not have an impact on the global biodiversity retrieved. Unknown OTUs also represented a large proportion of the fungal diversity in our study, which was comparable to other studies using similar high-throughput sequencing methodologies [11,17,67,68]. Despite all the progress in fungal DNA barcoding, most of the fungal species in the world remains unknown for many reasons. Some species are not easily identifiable, difficult or impossible to growth in vitro, they are only represented by only one discovered specimen and in some case it is truly the first encounter [69,70].

Abiotic Drivers of Fungal Richness
We found that predictors of the fungal richness were the stand age and the moss cover. The first ranked model included the variables stand age, average precipitation, herb cover, moss cover, soil pH and silt content, but all first six models were equivalent in weight and for Akaike's information criteria (AICc) ( Table S2). Among these variables stand age and moss cover come out as the most important predictors of fungal OTUs richness ( Figure 5).
For moss cover, a decrease of OTUs richness was observed when the moss cover was higher ( Figure 4A). Fungi living in moss covered stand, would have considerable potential to degrade carbon substrates contained in moss residues and in peatlands [71]. Mosses are well known to be recalcitrant substrates and they can affect conditions growth by maintaining high water level and slowing down the nutrient cycle [72]. They also produce secondary metabolites which would make them difficult to digest for most decomposers [73]. Only a few specialist fungi would be able to decompose moss tissues [74,75]. This could explain a lower fungal richness compared to other type of stands.
In addition to moss cover, stand age come out as an important predictor ( Figure 5). A decline of OTUs richness was observed in our older stand ( Figure 4B). A general trend, observed in animals and plants, is the colonization of young stand, normally newly exposed, by opportunistic species followed by an increase of species richness [76,77]. In the case of fungi, generally mature coniferous forests can support a larger number of fungal species and provide a higher fungal richness [78,79]. Blaalid et al. (2012) also observed for root-associated fungi a high fungal richness in the recently exposed area (young stand), followed by a decrease to finish with an increase of richness in more established ecosystems [80]. In these studies, only one fungal guild (mycorrhizal fungi) was taken into account. In our case the total fungal composition was measured and as such included several fungal groups. Another study compared the presence of several fungal groups over three successional forest stages (2 years, 65 years, and 110 years) and they detected fewer different groups of fungi (species richness) in the older forest than in the young and intermediate forest taken together [81]. In our data, old forests are more than 80 years old (Table S1). Such a high richness in early successional stands contrasts with other findings in western Canada where fungal richness increased over time after clearcuts or fires [82]. Early-successional stands may have increased primary productivity and also potentially higher plant diversity and theses two factors can determine the diversity of soil fungal guilds [83]. Mature forest may not be able to support a large diversity of fungi because of loss and changes in nutrient and plant species. Forest productivity will reach a plateau around 40 years old for most of the tree species [84,85]. After the maximum growth, we can observe a decrease of the primary production, photosynthesis, nutrients, and an increase of dead woods [86].
Other predictors appeared in more than half of all possible models as herbaceous cover, lichen cover, average precipitation, and silt content ( Figure 5). In the case of herbaceous cover, an increase in OTUs richness was observed with higher herb cover ( Figure 4C). Tedersoo et al. (2016), identified that fungal richness was positively affected by herb cover [18]. Several groups of root-inhabiting fungi are related to herbaceous plants and a majority of herbaceous plants are dependent on mycorrhizae to survive especially in phosphorus poor soil [87,88]. Higher density of herbaceous cover could provide more colonizable roots for fungi and results in a higher fungal richness [89]. For lichen cover, no information was found to explain how lichen cover can directly affect fungal richness. In our study, we evaluated cover from macrolichen living on the soil surface. One explanation would be that the type of stand where there is high abundance of lichen is not favorable for fungi. Lichen can grow on several extreme environments: on dry and hot environments, on bare rock, walls, gravestones, roof, exposed soil surfaces, etc. In our study, high lichen cover (covering more than 40% of the area) was found in stand plantation dominated by jack pine. Pine stand is known to have lower fungal richness because of forest litter that could be more acid [18,19,90]. Tree species can differentially influence fungal richness, either because of their partner selectivity or because of the local soil conditions they create [19,90]. We can see this variation in OTUs richness depending on the dominant tree species in our sites, with especially low richness under jack pine stands ( Figure 4D).
One edaphic factor that appeared to be important predictor is silt content, with an importance over 0.50 and was included in the lowest AICc model (Table S2 and Figure 5). Fungal response to soil texture has been found elsewhere [91]. Higher global fungal richness was previously observed in fine-textured soils [20]. A reason for that may be the higher propension of fine particles to be bound within aggregates. It is well known that aggregates formation can create microenvironments for soil microbes and create refuges from predators (e.g., collembolas) in inaccessible micropores [92], but also protection for desiccation [93], heat [94], and pH fluctuation [95]. Such creation of microenvironments may have favored fungal niche partitioning locally, thus potentially increasing richness [96]. Providing more isolate microhabitat demonstrated an increase in microbial species richness [97].
Average cumulative precipitation during the growing season was also an important predictor (Table S2 and Figure 5). Tedersoo et al. (2014), previously identified mean annual precipitation as a strongest predictor with a positive effect on global fungal diversity [17]. They compared fungal richness between tropical regions and Africa, the latter had a lower fungal richness because mean annual precipitation is relatively lower. Another study noticed that the number of soil fungal taxa increases significantly with increasing mean annual precipitation [98]. An explanation of the role of precipitation on richness would be related to the fact that most of the mesophilic fungi fail to grow at less than 96% humidity. A study from Talley et al. (2002), showed that relative fungal richness increase in presence of high relative humidity (>50%) [99].

Drivers of Fungal Community Composition
Soil pH was the main driver of fungal community composition (Table 1 and Figure 6A). This is consistent with previous studies showing that soil pH is the primary driver of fungal community composition and it is a main determinant on fungal community composition across the globe [11,17,67,[100][101][102]. Fungi generally grow, without significant inhibition, on a wide pH ranges [103][104][105]. However, soil pH could drive fungal composition by changing the abundance of some fungal groups [106], on microorganism diversity [107], and soil nutrients [108]. Acidic soils are characterized by a low cation exchange capacity, particularly in magnesium and calcium, and then cations are not available to organisms, which may act as a strong environmental filter to few adapted and specialized species [109,110]. Soil pH is here confirmed as must be part of the toolkit for modeling the mushroom niche, for conservation and management purposes.
Average precipitation was also identified as important drivers of fungal community composition (Table 1 and Figure 6B). In a previous study, fungi directly responded to variation in precipitation with more abundant, diverse, and consistent communities under drought conditions and more variable communities with an increase in rainfall, but difference in the community composition was more evident when rainfall difference was over 200 mm [111]. For our sites, average cumulative precipitation range was between 489.36 mm to 638.03 mm (Table S1). Fewer sites (blue) were at the opposite of the first PCoA axis corresponding to sites with an average cumulative precipitation over 581 mm (Tables S1 and S4). Over a specific stress threshold, biotic factors such as competition will be more present and determine coexistence [96]. Hawkes et al. (2011) observed increase diversity in dry periods and suggested that drought stress moderates competition among fungal taxa [111]. The species have therefore acclimated to specific precipitation conditions and to a certain level of competition and then specific groups of species should be indicative of some stands depending on precipitation level. A species indicator analysis was used to identify some indicator species between the different levels of precipitation. Fourteen genera were identified as significant indicators of low average cumulative precipitation and only three as mid average cumulative precipitation, compared to fifty-two for high average cumulative precipitation (Table S5). Species belonging to Inocybe, Trichoderma, and Clitocybe are examples of significant indicators of high average cumulative precipitation (Table S5).
For stand age, old stands (blue) seemed slightly nested compared to those from younger stands (red) ( Figure 6C). Stands could be different in community composition, because of some species unique to an age-class [47,112]. We identified two genus as an indicator of old stand (80 years and more), Gloioxanthomyces and Urnula (Table S5). Among all the OTUs, 590 OTUs were only detected in old stand (Tables S1 and S3). Difference in fungal species composition among forest age classes was previously found in several studies. They explained the difference by changes in plant species composition, because plant communities evolves with time [113,114]. Additionally, accumulation of some recalcitrant plant litter is not favorable for the survey of some fungi [115]. One study found that fungal diversity reached a plateau in 26 years-old stand and community composition was similar for stand between 65-100 years, but differ among stand over 100 years-old [82]. Another explanation could be after the maximum growth of a stand, a decrease of the primary production, photosynthesis, nutrients, and an increase of dead woods is observed and then competition to survive could be more present. Models of fungal primary succession showed that few selected species are able to colonize new area and with time they are joined or replace with more competitive species [116,117].
The fungal community structure is however characterized by a high turnover rate across the region (Figure 7). This specifically means that there were no cosmopolitan generalist fungal taxa that were present in the vast majority of our sites. Such a high replacement or turnover of fungal OTUs could reflect either (1) narrow fungal niches and strong environmental filtering [118], (2) strong interspecific competitive exclusion, or (3) historical contingencies with priority effects allowing early colonizers to prevent the establishment of latecomers [21,22,24]. Without additional data, it is impossible to determine in our dataset what the major cause of such a strong species turnover was.

Conclusions
A global overview of the fungal diversity and richness was assessed across the northern temperate and boreal Quebec forests. It allowed to localize area with high fungal richness and identified abiotic predictors across the sampled region. This could be used as a tool to follow fungal richness through time and to anticipate future changes that could be damaging to forests. A high moss cover, old stand (over 80 years) and some dominant tree species (pine stands) are not suitable for a high fungal richness contrary to herb cover and average precipitation. The structure of the fungal community in the region was unknown and now is it characterized by a high turnover rate. Fungal composition is shaped by several abiotic drivers, but we identified soil pH, average cumulative precipitation, and stand age as the most important drivers among the sampled variables. This effort to use meta-barcoding will allow a better characterization of fungal biodiversity and community structure from soil samples across various stands in the province of Québec. Collectively, our results show how difficult it can be to predict fungal community assembly even with high replication and efforts to include several biologically relevant explanatory variables: this poses a significant challenge for fungal conservation and biodiversity management in a changing climate, as it seems that suitable fungal niches will be hard to track. The inclusion of functional traits and phylogeny could be a way to better understanding fungal biogeography [119].
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1. K-means clustering of 249 sites based on the historical climatic data extracted from the site of the government of Canada using the GPS localization. Group 1 (black) includes 70 sites, mostly localized at the north of Abitibi-Témiscamingue (exception of four sites in Haute-Mauricie and two in south of Abitibi) and represents the domain NAT. Group 2 (red) has 79 sites situated in Mauricie and named domain MAU. Group 3 (green) includes 100 sites from the south of Abitibi-Témiscamingue (exception of five sites at the north) and Haute-Mauricie and represents the domain SAT-HM. Figure S2. Plots representing values from soil pH (a) and (b) average precipitation and in y axis, their distances to the centroid. Document file S1. Local fungal ITS DNA barcode reference; Document file S2. Multivariate dispersion results; Table S1. Environmental database; Table  S2. Dredge results; Table S3. OTU database; Table S4. Ordination scores; Table S5. Results of species indicators analysis.
Author Contributions: G.L. contributed to the overall process of the experimental design, field sampling, laboratory manipulations, data analysis and manuscript drafting; P.-L.C. contributed with his expertise, made comments for the community analysis and ecological parts, and helped for R coding and statistical analysis; R.G.-T. participated to field sampling, executed the soil analysis, helped for R coding, and made comments; A.M. participated to field sampling, executed the pH-Eh soil analysis, and made comments; D.B. contributed with his technical support and expertise in geomatics; V.M. contributed with his expertise, made comments for the statistical and ecological parts, participated in field sampling, and contributed to the experimental design; H.G. supervised the project, participated to field sampling, contributed to the experimental design, reviewed the manuscript, made comments, and obtained funding.

Funding:
The study was funded by the Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG) and the Consortium de recherche et innovations en bioprocédés industriels au Québec (CRIBIQ).