Vegetation drives the structure of active microbial communities on an acidogenic mine tailings deposit

Study area

The studied mine site is located in Abitibi-Témiscamingue, Québec (Canada), at the Southern limit of the boreal forest belonging to the bioclimatic domain Abies balsamea-Betula papyrifera. The forest is mainly composed of balsam fir (Abies balsamea), black spruce (Picea marinara), white spruce (Picea glauca), paper birch (Betula papyrifera) and Jack pine (Pinus banksiana) (Bergeron et al., 2004; Kneeshaw & Bergeron, 1998). The mine site is home to several plant species, which grow naturally on tailings and mainly include alder (Alnus incana ssp. rugosa), paper birch and spruce (Picea sp.) (Gagnon et al., 2020; Tardif et al., 2019). For this study, the mine site was organized into four classes constituting the vegetation density classes (VDC): VDC-1, VDC-2, VDC-3, and VDC-4 (Fig. 1). For comparison, two adjacent boreal forest areas were studied: VDC-5, which was near the mine site and VDC-6 at 20 km from the mine site.

Bulk soil and rhizosphere soil sampling

Sampling was conducted at the beginning of September 2015 and 2016. The number of sample replicates for plants (for rhizosphere and bulk soil) is detailed in Table S1. In VDC-2, 3 and 4, trees of each species (alder, paper birch and spruce) were chosen according to their availability and size, selecting primarily only small plants (approximately 0.5 m in height) (Table S1). These provided plant samples, and soil/tailings samples. To collect fifty gram soils/tailings samples in areas without plants (VDC-1), in the forest adjacent to the mine site (VDC-5) and in the distant far boreal forest (VDC-6), samples were randomly collected under the litter (if present) in the first 30 cm of soil. More specifically, bulk soil (and not rhizosphere soil) was collected in VDC-1, 5 and 6, due to the lack of plants (VDC-1) or the difficulty in VDC-5 and VDC-6 of harvesting the rhizosphere soil of mature trees in the forest. All plant and soil/tailings samples were kept at 4 °C until they arrived at the laboratory. Then, the roots were shaken in a beaker to collect the rhizosphere soil (approximately 20 g). The rhizosphere and bulk soil were frozen at −20 °C until RNA extraction.

Soil physico-chemical analysis

Soil pH, water content and elemental analyses by inductively coupled plasma mass spectrometry (ICP-MS) (except for total nitrogen and organic content) were performed previously (Gagnon et al., 2020). For elemental analysis, samples were digested using nitric acid (HNO₃). For plant tissues, this digestion resulted in complete mineralization of the samples. Digestion of soil samples with this strong acid (HNO₃) is an established method known to extract all metals that can become environmentally available, but does not extract metal included in the mineral fraction of recalcitrant minerals (e.g., silicates) (Osher et al., 2006; Tuovinen et al., 2015, U.S. EPA, 1996; Zhi Dang, Liu & Haigh, 2002). Digestions were carried out on a heating plate (DigiPREP Mini, SCP SCIENCE, Québec, Canada). Briefly, samples were crushed using a mortar and pestle in liquid nitrogen. Twenty-five mg of powder were transferred to 15 ml digiTUBES^® (SCP SCIENCE) to which 1 ml of concentrated nitric acid (HNO₃ Optima™, trace metal grade, Fisherbrand, Fair Lawn, NJ, USA) was added. Tubes were heated at 45 °C for 30 min. Then, another 1 ml of nitric acid was added, and the tubes were heated at 65 °C for 2 h. Digestates were brought to 10 ml with milliQ^® water (Millipore, Québec, Canada). Elemental quantitation was performed using a Thermo Fisher XSERIES 2 ICP-MS (Thermo Fisher Scientific, Bremen, Germany) and PlasmaLab software (Thermo Electron Corp. 2004. version 2.6.1.335. [PlasmaLab]. Germany: PC). Rhodium was used as an internal standard. The SLR5 certified water standard was used as an internal control to calibrate the ICP-MS. The studied metals were: sodium (²³Na), magnesium (²⁴Mg), aluminum (²⁷Al), phosphorus (³¹P), potassium (³⁹K), calcium (⁴⁴Ca), titanium (⁴⁷Ti), vanadium (⁵¹V), chromium (⁵²Cr), manganese (⁵⁵Mn), iron (⁵⁷Fe), cobalt (⁵⁹Co), nickel (⁶⁰Ni), copper (⁶⁵Cu), zinc (⁶⁶Zn), arsenic (⁷⁵As), selenium (⁸²Se), molybdenum (⁹⁵Mo), silver (¹⁰⁷Ag), cadmium (¹¹¹Cd), antimony (¹²¹Sb), barium (¹³⁷Ba), tungsten (¹⁸²W), thallium (²⁰⁵Tl) and lead (²⁰⁶Pb). Total nitrogen was determined by the combustion method (AOAC 990.03) and organic content was quantified by the calcination method (AGDEX 522, method PR-2) in an external laboratory (Groupe EnvironeX, https://www.labenvironex.com/en/). The results are presented in Table S2.

RNA extraction, library preparation and sequencing

Total RNA was extracted from environmental samples and then separated from DNA. RNA-based analysis was preferred over DNA-based analysis, as it provides data on microbial populations that are present and active, and not simply present (DNA).

Two grams of frozen bulk soil or rhizosphere was weighed into tubes containing 1.5 g of silica beads. Total RNA was extracted, purified and eluted as described in the RNA PowerSoil^® Total RNA Isolation kit (MO BIO Laboratories Inc., #12866-25, CA, USA). For removal of residual DNA, the RNA extracts were subjected to DNase treatment using the DNase I kit (Thermo Fisher Scientific, #18068015, Canada). Efficiency of DNA removal was verified through universal PCR amplification of the bacterial 16S rRNA genes and visualization of PCR product on 2% agarose gels to confirm the absence of DNA.

The 16S rRNA gene sequences (eubacteria and archaea) and the fungal nuclear ribosomal internal transcribed spacer (ITS) region were first amplified by RT-PCR using the OneStep RT-PCR kit (Qiagen, #220211, Canada). The specific primer sequences (Integrated DNA Technologies, Canada) for these regions (16S rRNA gene and ITS) are described in Table S3. These were chosen to identify bacterial and fungal classes and genera, as reported in the current literature (Kaiser et al., 2016; Sun et al., 2017; Urbanovà, Šnajdr & Baldrian, 2015; Žifčáková et al., 2016). For each region (16S rRNA gene and ITS), four pairs of primers with additional bases inserted between the sequencing primer overhang and the gene specific sequence (staggered pad) were used proportionally among the samples to create diversity in the sequence reads, thereby improving the quality of sequencing data. For amplification of the 16S rRNA region, a peptide nucleic acid (PNA) clamp was used to prevent amplification of eukaryotic RNA. PNA chloroplast and PNA mitochondrial blockers were used as proposed in Hong, Yang & Zhuang (2016).

Amplicons were visualized on a 2% agarose gel. The samples were purified using AMPure XP magnetic beads (Beckman Coulter Genomics, #A63881, CA, USA), and a second short amplification of 8 cycles was carried out with the Nextera XT index kit (Illumina^® #FC-131-1001, set A, B, C and D, CA, USA) to insert different index primers on the amplicons, according to the protocol in the Illumina “16S Metagenomic Sequencing Library Preparation” guide (Part #15044223 Rev. B). Unique codes were added to each sample by amplifying 5 µl of the purified PCR product with 25 µl of KAPA HIFI HotStart Ready Mix, 300 nM each Nextera XT Index Primer (Illumina^® Inc., CA, USA) and 10 µl UltraPure™ DNase/RNase-Free Distilled Water for a total volume of 50 µl. Thermal cycling conditions were as follows: 3 min at 98 °C, 8 cycles of 30 s at 98  °C, 30 s at 55 °C, 30 s at 72  °C, and a final elongation step of 5 min at 72  °C. Indexed amplicons were purified with magnetic beads as described above. Amplicon products were quantified using the Quant-it PicoGreen™ dsDNA Assay Kit according to the manufacturer’s instruction (Thermo Fisher Scientific, MA, USA) and combined in an equimolar ratio. The amplicon pool was verified with a Bioanalyzer 2100 (Agilent Technologies, CA, USA) to confirm the size of the amplicons and to visualize the potential presence of primer dimers or adapters, in which case an additional purification step of the amplicon pools with Select SPRI (Beckman Coulter, #B23319, CA, USA) beads was performed. Finally, the pool was denatured with 0.2 N NaOH, and the internal control PhiX was added at 5% (Illumina^® Inc., #FC-110-3001, CA, USA). The samples were sequenced on the MiSeq instrument (Illumina^® Inc., #SY-410-1003, CA, USA) using the MiSeq v2 500 cycle kit (Illumina^® Inc., #MS-102-2003, CA, USA) at the National Research Council Canada in Montréal, Québec (Canada).

Bioinformatics analysis

Sequencing results were analyzed using a previously described methodology (Tremblay et al., 2015; Yergeau et al., 2015). Briefly, both the 16S rRNA gene and ITS amplicons were separately filtered, assembled, trimmed and controlled for quality. Sequencing data was analyzed using AmpliconTagger (Tremblay & Yergeau, 2019). Raw reads were scanned for sequencing adapters and PhiX spike-in sequences and remaining reads were merged using their common overlapping part with FLASH (Magoč & Salzberg, 2011). Primer sequences were removed from merged sequences and remaining sequences were filtered for quality such that sequences having an average quality (Phred) score lower than 30 or one or more undefined bases (N) or more than 5 bases lower than a quality score 15 were discarded.

The generation of OTU tables was done with a three round clustering strategy. Quality controlled sequences were dereplicated at 100% identity. These 100% identity clusters were denoised at 99% identity using dnaclust v3 (Ghodsi, Liu & Pop, 2011). Clusters having an abundance of 3 or more reads were scanned for chimeras with UCHIME denovo and UCHIME reference-based algorithms, using the Broad Institute 16S rRNA gene Gold reference (http://microbiomeutil.sourceforge.net). The remaining clusters were clustered at 97% identity (dnaclust) to generate OTUs.

The taxonomic assignment of bacterial and fungal OTU results was performed with the RDP classifier (Bayesian classifier) with a training model constructed from the Greengenes database (version 13.5) (DeSantis et al., 2006; Wang et al., 2007). Fungal OTUs were classified using a training model constructed with the United Database (Kõljalg et al., 2013).

The RDP classifier gives a score (0 to 1) to each taxonomic depth of each OTU. Each taxonomic depth having a score ≥ 0.5 was kept to reconstruct the final lineage. Taxonomic lineages were combined with the cluster abundance matrix obtained above to generate a raw OTU table. From that raw OTU table, an OTU table containing bacteria only was generated. This latter OTU table was normalized using the edgeR R packages (Robinson, McCarthy & Smyth, 2010) to obtain a Counts Per Million (CPM) normalized OTU table. A multiple sequence alignment was then obtained by aligning OTU sequences on a Greengenes core reference alignment (DeSantis et al., 2006) using the PyNAST v1.2.2 aligner (Caporaso et al., 2010). Alignments were filtered to keep only the hypervariable region of the alignment. Alpha (observed species (Chao1 index)) and beta (Bray Curtis distances) diversity metrics and taxonomic summaries were then computed using the QIIME v1.9.1 software suite (Caporaso et al., 2010; Kuczynski et al., 2011) using the normalized OTU table.

Sequencing data are available on the NCBI Sequence Read Archive (SRA) portal under accession number PRJNA634113.

Statistical analysis

A matrix of environmental properties needed for the redundancy analysis (RDA) and the MANOVA was performed as described in Gagnon et al. (2020). Briefly, the species data matrix (y) was transformed by a centered log ratio since OTU data are considered compositional data (Macklaim, 2014). Since there were 174 bacterial classes and 1906 bacterial genera and 31 fungal classes and 889 fungal genera, for uniformity we chose to keep the 30 most abundant bacterial and fungal classes and genera for RDA. The RDA was performed using the “rda” function with the significant explanatory variables selected by the “ordiR2step” function. A second RDA was performed on the significant explanatory variables. The significance of the explanatory variable (red arrows) was verified using the “anova.cca” function (Table S4). The RDA are presented in type II scaling to allow a better interpretation of the correlation between matrices (x) and (y) (Legendre & Legendre, 2009).

Principal coordinate analysis (PCoA) was performed on normalized OTU tables. A distance matrix was created via the “vegdist” function, and PCoA was performed using the “pcoa” function. MANOVA was conducted to determine if significant differences in active bacterial populations and active fungal populations were present across VDC, or between plant species (terHorst, Lennon & Lau, 2014; Liu et al., 2015b; Wu et al., 2017; Warton & Hudson, 2004). Table 1 presents the p-values indicating the significance of differences observed between the active microbial populations between the various VDCs, and between plant species.

The microbial community profiles of the most abundant bacterial and fungal classes and genera were determined using the Shiny software (In-house R script) that allows the illustration of 20 classes.

To characterize the degree of change, and possible trends in microbial diversity across VDC, Chao1 boxplots were performed (GraphPad Prism 8 software for OS X, v8.1.2). We used the “Identify outlier” (method: ROUT, Q = 1%) function in GraphPad Prism to remove outliers from Chao1 datasets.

Statistical differences between vegetation density classes and plant species were analyzed in R using the “Kruskal.test”, “Dunn.test”, “lm” and “Anova” functions. Significant differences in bacterial richness between VDC and plant species were investigated using non-parametric Kruskal–Wallis followed by Dunn’s multicomparison test as the dataset did not meet the normality requirement after transformation. Significant differences in fungal richness between VDC were investigated using unbalanced ANOVA (sum of squares of type III) analysis followed by Tukey’s post-hoc test.

To determine what microbial populations were shared (vs unique) amongst VDC or plant species, the following online software was used: http://bioinformatics.psb.ugent.be/webtools/Venn/. The software sorted our metadata and OTU table (at the genus level), producing a list of most to least abundant populations. The 100 most abundant bacterial OTUs, and 100 most abundant fungal OTUs were then categorized into Venn tables that list which OTUs were shared or unique amongst VDC or plant species (Tables S5 and S6).

Vegetation density classes and plant species as drivers of microbial communities

Our results showed that the vegetation density classes (VDC) and plant species (p_s) significantly influenced the structure of the microbial (bacterial and fungal) communities present on the study site (Figs. 2 and 3). The vegetation density classes presented in this paper ranged from a total absence of plants (VDC-1) to a boreal forest, unaffected by human activities (VDC-6). In a mining context, acidic and poor-nutrient tailings make plant establishment difficult, leading to a low carbon concentration in soil, which was observed on tailings (VDC-1, 2, 3 and 4). These VDC had an average organic matter content of 2.1% compared to 22.39% for VDC-5 and 6 (Table S2). The presence of plants is known to change soil properties through the secretion of root exudates composed of organic acids, H ⁺, chelators, and phytosiderophores, which influence pH, the stability of aggregates, and nutrient availability (Anger & Caron, 1998; Gould et al., 2016; Philippot et al., 2009; Van der Putten et al., 2013). The absence of vegetation in VDC-1 leads to a different soil structure and lack of root exudates, which are truly important in the structure of the bacterial communities (Chaparro et al., 2013; Marschner et al., 2001). All these factors could explain why VDC and plant species were the main factors structuring microbial communities in RDAs, as well as the results observed in PCoAs where VDC-1 showed the greatest difference in the bacterial community structure in comparison to all other VDC (Figs. 2 and 3A and Table 1A). Our study also showed that fungal communities were different in the boreal forest soils (VDC-5 and 6) compared to those found on the tailings deposit (VDC-1 to 4) (Fig. 3B and Table 1B). This might be due to the differences in soil carbon levels. A study conducted by Liu et al. (2015a) as well as Sun et al. (2016) reported that carbon concentration in soil was indeed one of the major determinants of fungal community structure.

Several studies have shown that plant species is one of the major determinants of the microbial composition in the rhizosphere soil (Cavaglieri, Orlando & Etcheverry, 2009; Chaparro et al., 2013). Root exudates vary according to plant species, age of the plant, its state of stress, as well as its stage of development (Cavaglieri, Orlando & Etcheverry, 2009; Chaparro et al., 2013; Turner, James & Poole, 2013). It was also showed by Prescott & Grayston (2013), as well as Urbanovà, Šnajdr & Baldrian (2015) that trees mostly influenced composition of fungal and bacterial communities. Our results also showed differences in the rhizospheres of different plant species. PCoA results showed significant differences of bacterial communities in the rhizosphere soil of alder compared to birch and different fungal communities in the rhizosphere soil of alder compared to spruce (Fig. 3). We measured the highest relative abundance of Actinobacteria in alder rhizosphere (4.68%) in comparison to birch (3.30%) and spruce (3.74%) rhizospheres (Fig. 4 and Table S7A). As mentioned previously, plant species can harbour a rhizomicrobiome that is specifically beneficial (Berendsen, Pieterse & Bakker, 2012). The fact that alders establish a root nodule-forming symbiosis with Actinobacteria of the genus Frankia might explain the higher relative abundance of Actinobacteria in their rhizosphere. Additionally, this symbiosis requires specific nutrients (i.e., Mo, Mg, Fe and P) to support its nitrogen fixation activity, which might in turn influence the plant-microbe associations that tend to develop (Bélanger, Bellenger & Roy, 2013; Kabata-Pendias, 2001; Ngom et al., 2016).

Bacteria associated with vegetation density classes and plant species

This study revealed that the six most abundant classes of bacteria for each VDC and p_s were Gammaproteobacteria, Deltaproteobacteria, Alphaproteobacteria, Planctomycetacia, Acidobacteria, and Actinobacteria and represented approximately 50% of the microbial community classes (Fig. 4 and Table S7A). A study conducted by Mendez, Neilson & Maier (2008) reported that the main bacterial classes associated with iron-rich tailings with low pH (2.7) to moderate pH (5.7) were Gammaproteobacteria, Alphaproteobacteria, Acidobacteria, Actinobacteria, Nitrospira and Firmicutes, which corroborates some of our results. Hur et al. (2011) showed that zinc mine tailings with low pH, low organic matter as well as high Zn, Cd and Pb concentrations were dominated by Proteobacteria, Actinobacteria and Acidobacteria. The studies done by Hur et al. (2011) as well as Mendez, Neilson & Maier (2008) corroborate what was observed on the studied mine site, that Proteobacteria, Acidobacteria and Actinobacteria were the dominant taxa colonizing mine tailings.

As mentioned above, the rhizosphere soil is the nutrient rich, high biodiversity zone surrounding plant roots. The rhizosphere soil is a beneficial niche for many key microorganisms (Prashar, Kapoor & Sachdeva, 2014). In VDC-1 where there was little to no plant colonization and where pH was the lowest (3.9) compared to other VDC, we found the highest relative abundance of the chemoautotrophs Acidiferrobacter sp. and Leptospirillum sp. This could be expected since Acidiferrobacter sp. derive their energy through the oxidation of various inorganic compounds and Leptospirillum sp. does so through the oxidation of iron, which is present in high concentrations (31,000 ppm) in VDC-1 (Table S2) (Tyson et al., 2000). The low organic matter content of the tailings (average of 2.1% in (VDC1 to 4) (Gagnon et al., 2020) likely limited the development of heterotrophic microbial communities in tailings. This could have also contributed to the relative abundance of chemoautotrophs in VDC-1. Such chemoautotrophic organisms tend to be less active in the presence of plants because the latter secrete many carbon-rich molecules, which allows heterotrophs to dominate in their presence (Bais et al., 2006; Turner, James & Poole, 2013).

It is known that PGPR organisms facilitate plant growth in abiotic-stressed environments (Grandlic et al., 2008; Sasse, Martinoia & Northen, 2018; Thijs et al., 2016; Weyens et al., 2009; Zhuang et al., 2007). Berendsen, Pieterse & Bakker (2012), as well as Santoyo et al. (2016), showed that plants could select and shape their rhizomicrobiome as a function of their homeostatic states and for meeting their needs. Plants growing in metal-contaminated soils have different rhizomicrobiomes compared to plants growing in non-contaminated soils; they recruit beneficial microorganisms to support their growth and limit their metal uptake (Sasse, Martinoia & Northen, 2018; Thijs et al., 2016). Microorganisms in the rhizosphere soil can stimulate plant growth by providing nitrogen, solubilizing phosphorus, producing phytohormones (i.e., auxins, cytokinins), synthesizing siderophores, and protecting plants from pathogens (Glick, 2010; Prashar, Kapoor & Sachdeva, 2014).

We observed the presence of nitrogen-fixing, symbiotic organisms such as Bradyrhizobium sp., Frankia sp., and Rhizobium sp. in VDC-2, 3, and 4. These nitrogen fixers may have improved the competitiveness of plants harbored in these nitrogen-poor substrates (Bais et al., 2006; Ma et al., 2011; Mus et al., 2016; Wani, Khan & Zaidi, 2007) as VDC-2 to 4 contained less than 0.25 ppm N (Table S2). Bradyrhizobium sp. have been reported to benefit plants growing in metal-contaminated soil. For example, Bradyrhizodium sp. (vigna) was shown to improve growth, nodulation, nitrogen content, and reduce the uptake of Ni and Zn in green gram (Vigna radiata). This occurred through the production of indole-3-acetic acid, siderophores and ammonia (Ma et al., 2011; Wani, Khan & Zaidi, 2007).

Despite many PGPR traits reported in Pseudomonas sp. in contaminated soils (i.e., mine tailings), our study showed (Fig. 6E) that the relative abundance of Pseudomonas sp. was higher in boreal forest soils (VDC-5 and 6) compared to that of the tailings storage area (VDC-1 to 4) (Glick, 2010; Ma et al., 2011). Pseudomonas sp. are known for their production of organic acids that solubilize mineral phosphorus, and the production of phosphatase/phytase that hydrolyzes phosphate-organic compounds, increasing the availability of this limiting nutrient in soils (Hunter, Teakle & Bending, 2014; Richardson et al., 2009; Vacheron et al., 2013).

In our study, average phosphate concentrations were not systematically higher in boreal forest soils (VDC-5 and VDC-6; 90.8 and 57.9 ppm, respectively) compared to the tailings (VDC-1 to VDC-4; 18.7, 25.8, 62.2, and 35.1 ppm, respectively) however, higher phosphate concentrations did tend to occur in denser vegetation (VDC-3 to VDC-6) (Table S2). Soil phosphate concentrations in VDC-6 might be influenced by the rapid uptake of phosphate by mature trees, as the nutrient becomes available. Our observations are therefore not contradictory; the likely higher demand for phosphorus in the boreal forest vegetation could explain the higher abundance of Pseudomonas sp. in these soil (VDC-5 and 6), compared to mine tailings (VDC-1 to 4).

The nitrogen demand of alders is higher than that of other plant species. Despite other tree species, the absorption of nutrient by alder (N and P) is approximately in the same proportion, which could explain the high nitrogen demands by this species to support growth. To meet this high nitrogen demand, plant recruitment of beneficial microorganisms in the rhizosphere soil is primary (Lõhmus et al., 2006). Bradyrhizobium sp. were known for their symbiotic association with legume plant, and to promote nodulation in these plants. In addition, like other rhizobia, Bradyrhizobium sp. also have the capacity to fix atmospheric nitrogen and to make this nutrient available for other organisms (Saharan & Nehra, 2011). There is little literature on recruitment of Bradyrhizobium sp. by alder, but it is possible that this microorganism could be favorably recruited by alder (in comparison to birch and spruce) to improve its nitrogen acquisition. Furthermore, few authors have reported the possible role of Haliangium sp. as PGPR organisms. However, Kundim et al. (2003) as well as Ma et al. (2018) do mention that Haliangium sp. secrete antifungal molecules that limit the development of phytopathogens. The highest relative abundance of Bradyrhizobium sp. and Haliangium sp. in alder rhizosphere soil could also be explained by the different bacterial selection between different plant species (Lei et al., 2018; Urbanovà, Šnajdr & Baldrian, 2015).

In essence, our results indicate that bacterial richness was significantly lower in VDC-1 compared to all other VDC, including those of the natural environment (VDC-5 and VDC-6). Considering the increase of vegetation density as an indicator of environmental recovery, these results strongly suggest that bacterial diversity promptly recovers to natural levels following the establishment of vegetation.

Fungi associated with vegetation density classes and plant species

The ectomycorrhizal (ECM) fungi, arbuscular fungi (AM) and dark septate endophytes (DSE) are reported to enhance the phytoremediation potential of plants through elemental cycling and fungal-metal interactions (Deng & Cao, 2017; Pirttilä & Frank, 2011). ECM fungi include an estimated 6000 species (mostly Basidiomycetes), and establish symbiosis with only 5% of terrestrial plants, in few woody plant families and genera (i.e., Pinaceae, Fagaceae, Betula sp., Populus sp., and Alnus sp.), AM fungi comprise approximately 150 species of Zygomycetes (Deng & Cao, 2017). They are associated with herbaceous plants, and various woody plant families (Deng & Cao, 2017; Fortin, Plenchette & Piché, 2015; Landeweert et al., 2001; Mandyam & Jumpponen, 2005). The DSE fungi are reported to establish symbiosis with over 600 plant species, including plants that were not reported to be mycorrhizal. They are present in the rhizosphere soil and roots of plants colonizing metal-contaminated sites, and they are characterised as conidial and sterile fungal endophytes, which form melanised inter- and intra-hyphal structures (Mandyam & Jumpponen, 2005; Upson et al., 2009). Many fungal taxa are reported to establish symbiosis with plants colonizing metal-laden sites. For example, strains of genera Trichoderma sp., Fusarium sp., Aspergillus sp., and Cladosporium sp., colonized Portulaca plant, a heavy metal hyperaccumulator (Deng et al., 2014). Penicillum spp. and Trichoderma spp. are also the most frequently isolated fungi that attenuate heavy metal stress in plants (Babu et al., 2014a; Babu et al., 2014b; Deng & Cao, 2017; Khan et al., 2014).

ECM fungi can play a role in metal housekeeping through mechanisms such as precipitation, chelation, cell-wall binding, and the binding of metals by organic acids, polyphosphates, peptides and their transport through intracellular compartments (Colpaert et al., 2011). Moreover, glomalin (a glycoprotein) produced by some AM fungi increase heavy metal binding, reducing the uptake of heavy metals by the host plants (Bano & Ashfaq, 2013). In addition, ECM and AM fungi are involved in contaminant detoxification and mediate the nutritional status of heavy metals in plants (Barea, Azcon & Azcon-Aguilar, 2002; Harms, Schlosser & Wick, 2011; Luo et al., 2014; Thijs et al., 2016).

DSE are known to enhance mineral uptake of host plants, increase the utilization of various organic pools and modification of host water uptake. Indeed, the high melanisation of DSE allows them to resist severe drought and heat and increase plant growth under such conditions (Kilvin, Emery & Rudgers, 2013; Mandyam & Jumpponen, 2005; Pirttilä & Frank, 2011; Zhang et al., 2008).

Liu et al. (2015a) reported a shift in fungal communities as a function of organic carbon content in soils: the abundance of Agaricomycetes decreased when more organic carbon was present, and when there was an increase of Incertae sedis. Our results corroborate this; the relative abundance of Agaricomycetes in mine tailings (VDC-1 to 4) (average ∼20%) is higher than its relative abundance in boreal forest soil (VDC-5 and 6) (average ∼15%). We also observed, as Liu et al. (2015a), a higher relative abundance of Incertae sedis in boreal forest soils (VDC-5 and 6) (∼16%) compared to that in mine tailings (VDC-1 to 4) (∼5%) (Fig. 4 and Table S7B) . The higher organic matter content in boreal forest could also explain our PCoA results that showed a difference in fungal communities in VDC-5 and 6 in comparison to those in mine tailings soils (VDC-1 to 4) (Fig. 3B and Table 1B).

In our study, ECM Hypocrea sp. were associated more with plant rhizospheres from VDC-2 to 4, and could potentially be beneficial to plants on such sites. As an example, Morales-Barrera & Cristiani-Urbina (2015) found that Hypocrea tawa could be useful to detoxify Cr(VI)-contaminated wastewaters because of its capacity to reduce hexavalent chromium (Cr(VI)) to a much less toxic trivalent chromium form (Cr(III)). Hypocrea can solubilize metals (i.e., Cr, Ni, Cu and Pb) and, with other microorganisms, could immobilize these same heavy metals following their solubilisation (Congeevaram et al., 2007; Kumari et al., 2015). Based on the previous studies, Hypocrea sp. could have alleviated heavy metal stress, helping host plants grow on the tailings storage site.

Our study revealed a higher relative abundance of Mollisia sp. in VDC-1 and 2. This could be explained by the fact that Mollisia sp. is known for its high metal tolerance (Pirttilä & Frank, 2011), and that the low organic matter content in these VDC likely provided conditions allowing high bioavailability of metals (Table S2).

Our study also found a greater relative abundance of Zygomycota fungi (Dissophora sp. and Mortierella sp.) in boreal forest soil (VDC-5 and 6) compared to tailings (VDC-1 to 4). These observations corroborate those of Liu et al. (2015a) and Sun et al. (2016) according to which the relative abundance of Zygomycota (and specifically Mortierella sp.) increased with soil carbon content. The organic matter content of VDC-5 and 6 were on average 22.4%, while they were 2.1% on the tailings storage area VDC (Fig. 7 and Table S2).

The capability of plants species to shape their rhizosphere soil microbiome to adapt to various environmental conditions can lead to significant differences in the relative abundance of diverse fungi and bacteria (Berendsen, Pieterse & Bakker, 2012; Sasse, Martinoia & Northen, 2018). Prescott & Grayston (2013) and Urbanovà, Šnajdr & Baldrian (2015) reported that fungal communities in soil were largely explained by dominant trees. Urbanovà, Šnajdr & Baldrian (2015) reported that different proportions of arbuscular and ectomycorrhizal fungi were found under different tree species. In our study, we showed that birch rhizosphere present a higher relative abundance of the DSE fungus Mollisia sp. compared to the rhizosphere soil of alder and spruce (Fig. 9A). This finding is similar to that reported by Fernández-Miranda Cagigal (2017), where DSE were the dominant fungi of healthy fine root of, amongst other species, Betula pubescens.

Globally, we observed that fungal richness was not as dependent upon the presence or density of vegetation as was bacterial diversity. Fungal diversity did not differ significantly between VDC-1 to VDC-4, compared to that of boreal forest soils (VDC-5 and VDC-6). While bacterial richness in diversity increased following vegetation establishment, that of fungal communities remained relatively stable. These observations corroborate those of Sun et al. (2017) which found that bacterial diversity was 7.0- to 7.5-fold greater compared to fungal diversity in a recovering forest ecosystem.