The field experiment was set up in two borrow pits ( Figure S1A, B ) in the North Shore region of Québec, Canada (lat. 49°51′40′′ N, long. 68°51′40′′ W) that were excavated between 1964 and 1968, leaving large stripped surfaces surrounded by boreal forest stands ( Figure S1C, D ). The harsh micro-environmental conditions, combined with the absence of topsoil, have prevented the natural recolonization of these sites over the years. The climate is humid and cold with a mean annual temperature of 1 °C and total precipitation of 950 mm. In 1998, a tractor was used to level the ground before planting. Jack pine (Pinus banksiana) seedlings obtained from the Centre de production de plants forestiers de Québec (provenance x05-33-92) were planted in June 1999 at a planting density of about 7,000 trees per ha. Each seedling received 25 g of Nutricote ® 18-6-8 fertilizer (Plant Products Co. Ltd., Brampton, ON, Canada, Table S1 ), with 80% of N being released within a 70-day period at 25 °C. The experimental design included two treatments and one control per block, each replicated three times, for a total of nine experimental units (hereafter plots) each covering 50 m² (N = 9) and planted with 35 trees ( Figure S2A-D ). Two complete blocks were established in the first borrow pit (BP1-1 and BP1-2, 0.4 ha, Figure S1C ), and one complete block in the second borrow pit (BP2, 0.2 ha, Figure S1D ). The treatments ( Figure S2A-D ) were as follows: (1) plots blanketed with locally produced RCW (dry weight = 0.157 kg per L, 3 kg per m², thickness = ~2-3 cm, ~150 kg ~ plot), leading to high-productivity plots (HPP); (2) plots where ~3 – 4 L of RCW was applied only at the seedling base at the time of planting with the same wood mulch thickness as in HPP but over a smaller area (~0.2 m² per tree, ~19 kg per plot), resulting in low-productivity plots (LPP); and (3) unamended control plots with no surviving trees, resulting in null-productivity plots (NPP).

In September 2014, three soil stations were installed in each experimental unit to determine the nutrient stocks in the soil and tree tissues. The soil was sampled volumetrically at depths of 0–10 cm and 10–20 cm using a 4.5-cm diameter steel cylindrical corer. Samples were air-dried and sieved through a 2-mm mesh prior to analysis. On the HPP, diameter at breast height (DBH = 1.3 m) was recorded for all trees. Two trees were selected in each experimental unit to cover the full range of DBH values found on all planted plots, and they were then harvested and their stumps excavated. The branches and foliage were separated from the trunks in the field. The samples were brought to the laboratory, where they were oven-dried and weighed according to the type of material sampled (foliage, branches, trunk and stump). Second-order polynomial equations were fitted to the DBH data to estimate dry mass for each of the four material types. The DBH of all the trees was recorded to estimate the biomass and the amount of nutrients that accumulated in the trees. The nutrient content of the tree tissues was determined from samples oven-dried at 70°C. Nutrient accumulation in the soil (total N and available P and exchangeable cations) and in the vegetation during the 16 years of ecosystem development in the HPP was estimated by assuming that the NPP represented the initial conditions for all sites since time zero measurements were not available ( Table 1 ). Given that these sites were very poor in nutrients and organic matter with no vegetation colonization ( Figures S2B, E ), it is reasonable to assume that the soil conditions in NPP did not evolve significantly over the 16 years of the experiment. RCW was produced from local small deciduous trees (Alnus, Betula) and its nutrient content was assessed from ten subsamples using the methods described below for vegetation. Total C content and total N content in soil and vegetation samples were analyzed using a LECO TruMac CNS analyzer (LECO Corporation, St. Joseph, MI, USA). Following calcination of plant tissues at 500°C and dilution with hydrochloric acid, cations and (P) were determined by inductively coupled plasma optical emission spectroscopy (ICP-OES) (Miller, 1997). A Mehlich III solution was used to extract base cations and available P (Carter and Gregorich, 2007), which were determined by ICP-OES. Finally, soil samples were also collected simultaneously with samples used to characterize the soil microbiome (see below) during the sampling campaign conducted in October 2016. Soil nutrients were measured as previously described, and pH was analyzed in distilled water and in 0.01M CaCl₂.

To characterize the soil microbiome, bulk soil, rhizosphere and roots were sampled at the end of October 2016. In the LPP and HPP, soil cores were collected at each cardinal point around three randomly selected trees (sampling unit), at a distance of 30 cm from the base of the trunks. In the NPP, soil cores were collected at each cardinal point around three randomly selected points because no trees had survived ( Figures S2E-G ). The roots in the soil cores from the HPP and LPP were separated from the bulk soil. Sieved bulk soil and root-attached soil (i.e., rhizosphere) were used to characterize the soil and rhizosphere compartments, respectively ( Figure S3 ). Samples were kept in a cooler during fieldwork and were stored at –80 °C in the laboratory until DNA extraction.

In the step prior to DNA extraction, root fragments were cleaned in 150-ml glass test tubes containing 15 ml of sterile phosphate buffered saline (PBS) solution (1X). Samples were vortexed for 10–15 s and sonicated 5 times for 30 s with a 30 s interval between each sonication. This step was repeated using new tubes and fresh PBS 1X solution. Cleaned root fragments were removed from the tubes, blotted dry and cut up. Up to 100 mg of roots were then weighed and placed in a 2-ml screw cap tube with a 5-mm stainless steel bead. The roots were flash frozen in liquid nitrogen and ground into a fine powder using a TissueLyzer II (QIAGEN, Valencia, CA, USA) at 26 Hz for 2 × 45 s. DNA was then extracted from the powder using the MO BIO PowerSoil DNA Isolation Kit (Mo Bio Laboratories Inc., Solana Beach, CA, USA). A total of 250 mg of starting material was used to isolate DNA from bulk soil and rhizosphere soil using the same extraction kit. Crude genomic DNA (gDNA) extracts were quantified using a Qubit Fluorometer 2.0 (Life Technologies, Burlington, ON, Canada) and the Qubit dsDNA HS assay kit. The concentration of gDNA was standardized by diluting each sample to 5 ng/μl (samples with concentrations lower than 5 ng/μl were not diluted).

Primer pair 515F-Y (5’-GTGYCAGCMGCCGCGGTAA-3’) and 926R (5’-CCGYCAATTYMTTTRAGTTT-3’) (Parada et al., 2016) was used to amplify the bacterial V4–V5 regions of the nuclear 16S ribosomal gene. The primer pair ITS9 (5’-GAACGCAGCRAAIIGYGA-3’) (Menkis et al., 2012) and ITS4 (5’-TCCTCCGCTTATTGATATGC-3’) (White et al., 1990) was used to amplify the fungal ITS2 region of the nuclear ribosomal gene. The amplification was performed in 25 µl of reaction mix in triplicate as follows: 9 μl of UltraPure™ Dnase/Rnase-Free distilled water (GIBCO, life technologies), 200 µM of each dNTP, 1.5 mM of Mg²⁺, 200 NT of each primer, 1 U of HotStarTaq Plus DNA Polymerase (Qiagen, Valencia, CA, USA) and ~ 12.5 ng of gDNA. Thermocycling conditions were as follows: initial denaturation step at 95°C for 5 min, 34 cycles at 94°C for 30 s, 50°C for 30 s, and 72°C for 1 min, and a final elongation step at 72°C for 10 min.

For high-throughput multiplexing and sequencing of amplicons, triplicate PCR products from each sample were pooled, purified and indexed. Residual dNTPs, primers and buffers were removed in a purification step using 81 μl of Agencourt AMPure XP beads (Beckman Coulter Inc., Indianapolis, IN, USA) according to the protocol described in Illumina’s 16S Metagenomic Sequencing Library Preparation Guide (Part # 15044223 Rev. B). Purified amplicons were indexed in 50 µl of the following reaction mixture: 10 μl of UltraPure™ Dnase/Rnase-Free distilled water, 5 μl of each Nextera XT v2 Index Primer (Illumina Inc., San Diego, CA, USA), 25 μl of KAPA HiFi HotStart Ready Mix (KAPA Biosystems, Wilmington, MA, USA), and 5 μl of the purified PCR product. Thermocycling conditions were as follows: initial denaturation step at 95°C for 3 min, 8 cycles of 95°C for 30 s, 55°C for 30 s, 72°C for 30 s, and a final elongation step at 72°C for 5 min. Indexed amplicons were purified using 56 μl of Agencourt Ampure XP beads as described above, quantified using a Qubit Fluorometer 2.0 with the Qubit dsDNA BR assay kit, and multiplexed using an equimolar ratio. Paired-end sequencing (2 × 250 bp) of the 16S/ITS pools was performed on an Illumina MiSeq System sequencer and the MiSeq Reagent Kit v2 (500 cycles) on the Nucleic Acid Solutions Group’s Illumina sequencing platform at the Aquatic and Crop Resource Development Research Centre, National Research Council Canada, Saskatoon, SK, Canada. The Illumina data generated in this study have been deposited in the National Center for Biotechnology Information’s (NCBI’s) Sequence Read Archive under accession number PRJNA775865.

Bacterial and fungal biomass in the rhizosphere and bulk soil and the abundance of the nifH gene, which encodes the nitrogenase enzyme and is responsible for N₂ fixation in bacteria, were measured using a qPCR assay. Primer pair 27F (5’-AGAGTTTGATCMTGGCTCAG-3’) and 518R (5’-TTACCGCGGCTGCTGG-3’) (Okubo et al., 2015) was used to amplify the 16S rRNA gene. The primer pair FungiQuant-F (5’-GGRAAACTCACCAGGTCCAG-3’) and FungiQuant-R (5’-GSWCTATCCCCAKCACGA-3’) (Liu et al., 2012) was used to amplify the 18S rRNA gene. The primer pair Ueda19F (5’-GCIWTYTAYGGIAARGGIGG-3’) and UedaR6 (5’-AAICCRCCRCAIACIACRTC-3’) (Angel et al., 2018) was used to amplify the nifH gene. The 10 μl reaction mix consisted of 5 μl of 2× QuantiTect SYBR Green PCR Master Mix, 2 μl of a 5 μM primer mix, 1 μl of Rnase-free water, and 1 μl of DNA extract (~ 0.5–5 ng of DNA). Thermocycling conditions were as follows: initial denaturation step at 95°C for 15 min; 40 cycles at 95°C for 15 s; 30 s at 54°C for the 16S rRNA gene, at 59°C for the 18S rRNA gene and at 52°C for the nifH gene; and 65°C for 90 s. Fluorescence measurements were performed at the end of the elongation step of each cycle, and the specificity of each reaction was verified by melting curve analysis. Standard curves were used to determine the copy number of each gene in the DNA extracted from the rhizosphere and bulk soil. The curves were calculated using a 10-fold dilution series of quantified DNA from pure cultures of Sinorhizobium meliloti 1607 (for the 16S rRNA and nifH genes) and Ophiostoma novo-ulmi H327 (for the 18S rRNA gene). For all assays, qPCR was performed on an ABI 7500 apparatus (Applied Biosystems; Thermo Fisher Scientific Inc., Waltham, MA, USA) using a QuantiTect SYBR Green PCR Kit (QUIAGEN, Valencia, CA, USA) according to the manufacturer’s instructions.

To generate microbial species tables using the raw demultiplexed sequences, paired-end sequences were first trimmed, truncated, assembled, filtered for chimeras, and denoised using the DADA2 plugin (Callahan et al., 2016) as established in QIIME 2 v2019.1.0 (Bolyen et al., 2019). Sequences with median quality scores below 37 were trimmed. The quality scores of the ITS2 sequences from soil samples decreased from positions 140–160 bp onwards, preventing the assembly of the forward and reverse sequences. Therefore, only the forward sequences were included in the analysis of the soil fungal community. For the libraries from root and rhizosphere samples, the 3’ends of the forward and reverse sequences were truncated, and the average length of the assembled sequences was 410–415 bp, with a 21–26 bp and 50 bp overlap in the 16S and ITS2 sequences, respectively.

Amplicon sequence variants (ASVs) were generated by clustering 16S sequences (pairwise similarity threshold of 100%), and operational taxonomic units (OTUs) were generated by clustering ITS sequences (pairwise similarity threshold of 97%) as proxies for bacterial and fungal species, respectively. The number of sequences used to train the error model was set to 2,000,000. ASVs or OTUs with a frequency of less than 0.1% of the mean sample depth were filtered out, which excluded the results caused by the MiSeq bleed-through between runs as reported by Illumina (https://github.com/LangilleLab/microbiome_helper/wiki/Amplicon-SOP-v2-(qiime2-2019.7)#41-filter-out-rare-asvs). De novo clustering was performed using VSEARCH (Rognes et al., 2016) in QIIME 2 to cluster the fungal ITS2 sequences (97% similarity threshold) into OTUs. Taxonomic information was assigned using a out Bayes classifier (Wang et al., 2007) that was pre-trained using the Greengenes (version 13_8, files clustered using a 99% similarity threshold, McDonald et al., 2012) and UNITE (version 8.0, released on 04/02/2020, Nilsson et al., 2019) databases for bacterial and fungal datasets, respectively. Each dataset was then filtered based on taxonomic criteria. A total of 184 ASVs that were identified as ‘Archaea’, ‘Chloroplast’ and ‘Mitochondria’ were discarded from the 16S dataset. A total of 115 unassigned and non-fungal OTUs were filtered out from the ITS2 dataset. Figure S4 shows how the number of samples, sequences, bacterial ASVs, and fungal OTUs changed during the bioinformatics processes, from raw data to ASV/OTU tables used for downstream community ecology analyses.

Rarefied ASV/OTU tables ( Figures S5, S6 ) were used to calculate diversity indices and relative abundances of bacteria and fungi. Alpha diversity (within-sample diversity) was calculated using the sample-size- and coverage-based approaches for the rarefaction and extrapolation of the ASV/OTU-generated Hill numbers (i.e., richness, q = 0; Shannon diversity, q = 1, the exponential of Shannon entropy; and Simpson diversity, q = 2, the inverse of Simpson concentration) in the R package for iNEXT v2.0.19 (Chao et al., 2014; Hsieh et al., 2016). The Hill numbers, or effective species numbers, are related to Renyi’s definition of a generalized entropy formula, and they differ only by an exponent q, which determines the sensitivity to rare (low q values) and abundant species (high q values) (Hill, 1973). Faith’s phylogenetic diversity (Faith and Baker, 2006) was calculated using a subset of the fungal dataset that included only OTUs identified as agaricomycetes. For the bacterial and agaricomycetes datasets, the QIIME 2 pipeline align-to-tree-mafft-fasttree in the phylogeny plugin was run to align the sequences using MAFFT (Katoh et al., 2002) and to produce a phylogenetic tree (root: midpoint) using FastTree (Price et al., 2010). The phylogenetic trees were then used as input to calculate Faith’s phylogenetic diversity. To generate the bacterial and fungal taxonomic profiles for the NPP, LPP, and HPP, the relative abundances of ASVs and OTUs were calculated, and ASVs and OTUs with a relative abundance of less than 1% were removed. Bacterial and fungal functions were predicted using the Functional Annotation of Prokaryotic Taxa (FAPROTAX) database (www.zoology.ubc.ca/louca/FAPROTAX) (Louca et al., 2016) and FUNGuild (http://www.funguild.org) (Nguyen et al., 2016), respectively. Only guild assignments that were ‘highly probable’ and ‘probable’ were considered to avoid overinterpretation of the ecological data (Nguyen et al., 2016).

To identify differentially abundant bacteria and fungi in LPP and HPP, a differential ranking procedure was applied to the non rarefied ASV/OTU matrices using Songbird v1.0.3 (Morton et al., 2019) as established in QIIME 2. Because Songbird is sensitive to sparsity, only ASVs and OTUs that were observed with a total abundance of more than 5,000 sequences across all samples and that were found in at least 10 samples from LPP and HPP were considered. Songbird allows for the identification of differentially abundant features while accounting for the compositional nature of microbiome data and accounting for confounding variables in the datasets. Differentials are defined as the estimated log-fold change in the abundance of features across different sample types. Multinomial regressions were performed while incrementing the number of covariates of interest (treatment, soil compartment, site) in the model. The convergence statistics (graph of the prediction accuracy of the model and graph of the loss function) generated by running multinomial regressions were visualized and analyzed using QIIME 2 Visualization. Parameters (covariates, number of iterations, number of samples assessed per training iteration, width of normal prior for the differentials) were adjusted to obtain a reasonable fit. The retained model was then validated against the null and baseline models. Qurro v0.7.1 (Fedarko et al., 2020), as implemented in QIIME 2, was used to visualize rankings and sampled plots and to explore feature rankings and log-ratios.

The bulk soil, rhizosphere and root datasets were combined prior to performing the core microbiome analyses. Core bacterial and fungal communities ( Figure S8 ) were calculated using the rarefied datasets (raw abundance) and the core and plot_core functions in microbiome R package v1.8.0 (Lahti and Shetty, 2019). Detection and prevalence thresholds were set to 0 and 70, respectively. Interactions between bacterial ASVs and fungal OTUs in LPP and HPP were modeled as Poisson log-normal (PLN) networks (Chiquet et al., 2019a; Chiquet et al., 2021) using the PLNnetwork function in R package PLNmodels v0.10.6. Compared to other models, the Poisson Log-Normal models have the advantage of accounting for confounding covariates and differences in sampling effort, while correcting for the compositional nature of count data, and integrating different data sets into the same network (i.e., bacterial and fungal communities in this study). Model selection was performed by incrementally adding covariates, and the model with the lowest Bayesian Information Criterion (BIC) was selected for the PLN inference. The relative log expression (RLE) normalization scheme was applied to compute scaling factors used as offsets during PLN inference and to account for differences in sequencing depth between samples.

To compare the means of Hills numbers (q = 0, q = 1, q = 2) and Faith’s phylogenetic diversity from each level of productivity, Poisson (q = 0) and linear (q = 1, q = 2, and Faith’s phylogenetic diversity) mixed ANOVA models were fitted separately for each compartment (soil, rhizosphere, root). ‘Block’ and ‘restoration treatments nested within block’ were defined as nested random intercepts, and ‘restoration and ecosystem productivity level’ as fixed effects, using the lme4 package (Bates et al., 2015).

Statistical differences in the average composition of bacterial and fungal communities between (1) productivity levels (NPP, LPP, HPP), (2) the compartments (soil, rhizosphere, root), and (3) their interactions were tested using permutational multivariate analysis of variance (PerMANOVA) based on Bray-Curtis dissimilarities calculated with log(x+1)-transformed abundance data (Anderson et al., 2008). Note that because only mineral soil was available in the NPP, this interaction term was calculated only for the LPP and HPP. To test for significant differences in multivariate dispersion between groups (i.e., beta diversity), permutational analyses of multivariate dispersion (PERMDISP) were used with log(x+1)-transformed abundance data and with the same factors used above for PerMANOVA analyses (Anderson et al., 2008). Finally, we used distance-based redundancy analyses (dbRDA) to analyze the relationship between multivariate bacterial and fungal abundance data (similarity matrix) and continuous soil chemical variables measured in the 20 cm topsoil layer. PerMANOVA, PERMDISP and dbRDA analyses were all performed using the Primer 7 software (Anderson et al., 2008).

To test whether concentrations of each soil nutrient differed among productivity levels, we first fitted linear mixed models with ‘block’ and ‘restoration treatments nested within blocks’ as nested random intercepts, and ‘productivity level’ as a fixed effect. Second, when the treatment effect was significant (or close to being significant), we performed Tukey’s multiple comparison test. Statistical analyses were performed in R v3.6.3 (R Core Team, 2020) using the packages nlme v3.1-152 (Pinheiro et al., 2020), multcomp (Hothorn et al., 2008), lsmeans v2.30 (Lenth, 2016), and lme4 1.1-27 (Bates et al., 2015).

A Poisson generalized linear mixed-effect model was used to assess the effect of forest restoration success and productivity on each gene copy number separately for each compartment, with ‘block’ and ‘restoration treatments nested within blocks’ as nested random intercepts and ‘productivity levels’ as a fixed effect. Pairwise comparisons of means were made using Tukey’s test. Statistical analyses were performed in the R package using lme4.

Data visualization (bar plots, non-metric multidimensional scaling [NMDS] biplots, core microbiomes, and networks) was performed using R and the following packages: corrplot v0.84 (Wei and Simko, 2017), ggplot2 v3.3.2 (Wickham, 2016), ggthemes v4.2.0 (Arnold, 2021), graphics v3.6.3 (R Core Team, 2020), gridExtra v2.3 (Auguie, 2017), igraph v1.2.5 (Csardi and Nepusz, 2006), knitr v1.29 (Xie, 2014; Xie, 2015; Xie, 2022), microbiome v1.8.0 (Lahti and Shetty, 2019), phyloseq v1.30.0 (McMurdie and Holmes, 2013), PLNmodels v0.10.6 (Chiquet et al., 2018; Chiquet et al., 2019a; Chiquet et al., 2019b; Chiquet et al., 2021), plyr v1.8.6 (Wickham, 2011), reshape2 v1.4.4 (Wickham, 2007), rmarkdown v2.3 (Xie et al., 2018; Xie et al., 2020; Allaire et al., 2022), and tidyr v1.1.0 (Wickham and Girlich, 2022).

Each of the three levels of forest restoration success and productivity produced a unique belowground microbiome ( Figure 1 ) with contrasting levels of microbial biomass ( Figure 2 ). The average taxonomic and functional composition of the bacterial and fungal communities differed markedly between the LPP and HPP ( Table 2 ), and this effect depended on the compartments (p< 0.01 for most interaction terms ‘Treatment × compartment’). The level of plot productivity (LPP vs. HPP) had a significant effect on the average taxonomic and functional composition within each compartment at p< 0.005 for most pairwise comparisons. The largest differences were observed in the rhizosphere for bacteria (SS = 40659; p< 0.001) and fungi (SS = 36472; p< 0.001), and in the bulk soil for fungal functions (SS = 12209; p = 0.006, Table 2 ). No significant or biologically relevant interaction was detected for the variance in the taxonomic and functional composition (i.e., beta diversity) of the bacterial and fungal communities ( Table S2 ). Therefore, the observed differences in community composition between plots are due to differences in average composition rather than differences in variance.

The taxonomic profiles based on relative abundance differed less at higher taxonomic levels ( Figure S9A, B ). Regarding the ecological functions, the predicted functional profiles of fungi were more similar between productivity levels than those of bacteria ( Figure S9C, D ). For bacteria, a decreasing trend was found for chemoheterotrophic functions between the HPP, LPP and NPP, while an increasing trend was found for nitrogen cycling. On the other hand, nitrogen fixation alone showed an opposite trend: it was predicted with a relative abundance of 1.74%, 0.42% and 0.65% in the HPP, LPP and NPP, respectively, and it was systematically higher in the HPP than in the LPP for the rhizosphere and root compartments. This trend was supported by the real-time qPCR analysis, which showed that the number of copies of nifH per gram of soil or root sample was 6 to 7 times greater in the HPP than in the LPP and NPP ( Figure 2 ). Finally, the bacterial and fungal biomass values obtained for the soil and rhizosphere ( Figure 2 ) in the HPP were 3 to 6 times higher than in the LPP, and up to 19 times higher in the HPP than in the NPP for soil fungal biomass (p< 0.01).

The microbiome composition in the HPP was associated with higher concentrations of C, Ca, K, and Mg and lower pH values ( Figures 1 , 3 ). The median concentration of these elements was about three times higher in the HPP than the LPP and NPP (p ≤ 0.06), while the sum of the vegetation and soil nutrient stocks was higher in the HPP than in the NPP ( Figure 4 , Table 1 ). C as well as N, P and K stocks were significantly greater in the HPP than the NPP (p ≤ 0.01). In the HPP, the tree tissues stored more C, N, K, and Mg than the soil and were therefore the main reservoir of accumulated carbon and nutrients, except for P. Phosphorus was mainly stored in the soil in the HPP ( x ¯ = 4.4 g/m^-2), while a similar amount of P was recorded in tree tissues in the HPP ( x ¯ = 2.9 g/m^-2) and soil in the NPP ( x ¯ = 2.6 g/m^-2). The average rate of nutrient accretion ranged from 1.2 kg ha^-1 y^-1 for Mg to 12.5 kg ha^-1 y^-1 for N ( Table 1 ) in HPP. Assuming an unlikely 100% retention over the 16 years of the experiment, the RCW amendment would account for 24%, 14%, 42%, 90%, and 56% of the accumulated N, P, K, Ca and Mg, respectively ( Table 1 ).

The levels of forest reclamation and tree productivity had contrasting effects on the alpha diversity of bacteria and of fungi ( Figure 5 ). The diversity indices obtained for soil bacteria were greater in the NPP than in the HPP and LPP (p< 0.05). For example, the effective number of species derived from Shannon entropy (q = 2) recorded in the soil samples decreased by 20% to 26% from the NPP to the HPP (NPP: 185 ± 49, LPP: 147± 31, HPP: 108± 35). However, this trend was not consistent between soil, rhizosphere, and root samples: the bacterial diversity in root and rhizosphere tended to increase with tree productivity, but the differences were not significant at p = 0.05. In contrast to bacterial diversity, the highest values of fungal diversity were recorded in the HPP, and this was consistent between the soil, rhizosphere, and root samples. However, only the species’ richness of the soil fungi was significantly higher (p< 0.05) in the HPP (77 ± 15) than in the LPP (48 ± 12) and NPP (45 ± 14).

Samples from LPP and HPP could be distinguished on the basis of the differentially abundant bacteria and fungi ( Figure 6 ). The functional screening of 16S-derived bacterial data was performed using FAPROTAX, which showed that decomposers (Sinobacteraceae), heterotrophic ammonia oxidizing bacteria (Xanthomonadaceae and Chitinophagaceae), and denitrifying bacteria (Comamonadaceae and Actinomycetales) were significantly more abundant in the HPP ( Table S3 ). The LPP provided a favorable niche for taxa related to the Mycobacteriaceae, Koribacteraceae, and Sinobacteraceae and as well as for Chitinophagaceae and Burkholderiaceae. Ectomycorrhizal fungi related to the genus Cortinarius dominated the top contributing OTUs in the HPP ( Figure 6B , Table S4 ). The log-ratio of Cortinarius OTUs classified in the top 10% of OTUs showed a clear separation between the LPP and HPP (inset boxplot, Figure 6B , Table S5 ), indicating that Cortinarius taxa were key players in the HPP, both in terms of diversity and abundance.

The core and differentially abundant taxa had a low connectivity in the LPP ( Figure 7 ). In the HPP, the network based on the core microbiome had approximately 10 times more edges (110 versus 12) than the network in the LPP. For the differentially abundant taxa, the number of edges was twice as high as in the HPP (42 edges) than the LPP (21 edges). Nine out of ten Cortinarius species were interconnected and two bacterial taxa related to the Chitinophagaceae were positively associated with six out of ten Cortinarius species. In the LPP, Lactarius sp. (OTU059) and Rhizopogon pseudoproseolus (OTU027) were most associated with the other differentially abundant fungi and bacteria. Two bacterial taxa related to the Sinobacteraceae (ASV0070) and the genus Burkholderia (ASV0036) were positively associated with Lactarius sp. (OTU059), Rhizopogon pseudoproseolus (OTU027) and Tomentella sp. (OTU046).

The contrasting levels of productivity observed in each experimental unit demonstrated the importance of applying RCW as a continuous layer over plots to be reforested in borrow pits. The wood mulch likely provided protection from the extreme abiotic conditions and its decomposition provided the nutrients necessary for seedling establishment. This allowed the seedlings to promote enhanced mineral weathering of the coarse sand particles found in the borrow pits through the recruitment of a specialized soil microbiome. Ectomycorrhizal fungi (van Schöll et al., 2008), bacteria (Uroz et al., 2009) and their interactions (Koele et al., 2009; Fontaine et al., 2016) contribute to accelerated mineral weathering rates. Their activity supported the development of trees measuring up to 6 m tall, 17 years after reclamation of the sand borrow pits began. The Hubbard Brook sandbox experiment showed that unvegetated soils had negligible weathering, while mineral weathering for K, Ca and Mg was high on vegetated plots (Bormann et al., 1998). The rate of mineral weathering required for ecosystem nutrient accretion appears to be high in the HPP compared to the long-term mineral weathering rates in granitic soils observed by Vadeboncoeur et al. (2014). For example, these authors reported Ca rates ranging from 0.1 to 4 kg ha^-1 y^-1, whereas Ca accretion rates in the HPP were approximately 9 kg ha^-1 y^-1 ( Table 1 ). Our results show significant nutrient inputs from mineral weathering in the upper portion of the soil, even after accounting for inputs from wood mulch and Ca in precipitation (which are likely less than 1 kg ha^-1 y^-1, Vadeboncoeur et al., 2014). In addition, soil acidification in the HPP, which was likely due to EM fungi, which accounted for the highest fungal biomass (Bahram et al., 2020), and root exudates and decaying soil organic matter, which released organic acids (Jones, 1998), also contributed to mineral weathering.

With respect to nitrogen, the ecosystem rate of N accretion in the HPP (12.5 kg ha^-1 y^-1) was six to twelve times higher than the estimated atmospheric N deposition rate for this region (1 and 2 kg ha^-1 y^-1 according to the Oak Ridge National Laboratory [ORNL] Distributed Active Archive Center [DAAC], data from July 15, 2016). The N inputs through the fertilizer applied at planting ( Table S1 ) in addition to that contained in the wood mulch ( Table 1 ) amount to at most to 50% of the N accretion, assuming an unrealistic level of 100% fertilizer retention in the HPP and noting that fertilizer was also applied to all planted trees in all treatments. Therefore, some of the N accretion in the HPP was likely the product of non-symbiotic N fixation (Johnson and Turner, 2014). For example, N-fixing activity has been associated with the ectomycorrhizal fungi Suillus variegatus, S. tomentosus and Tomentellopsis submollis on pine roots (Izumi et al., 2006; Paul et al., 2007). Our results showed that Suillus tomentosus (OTU007) and S. grevillei (OTU029) were part of the core fungal community in the HPP ( Figure S8 , Table S8 ). Similarly, families of free-living and plant symbiotic nitrogen fixers (Bradyrhizobiaceae, Burkholderiaceae, Methylocystaceae, Rhizobiaceae) were dominant in the HPP. The genera Bradyrhizobium and Burkholderia are generally the most abundant in soil microbiota (Shah and Subramaniam, 2018), and ASVs assigned to these genera dominated the core bacterial community in each plot. However, the number of copies of the nifH gene indicated that taxa from these groups thrived only in the HPP. EM fungi can utilize both inorganic and organic nitrogen, and efficiently immobilize nitrogen and redistribute it to their hosts through common mycorrhizal networks (He et al., 2003; Sterkenburg et al., 2018). In addition, EM fungi likely received more carbohydrates from the vigorous jack pines than from the stunted pines, which would increase mycelium biomass and nitrogen retention in mycelium (Näsholm et al., 2013).

Although the post hoc nature of the current study prevents a precise assessment of nutrient inputs (lack of information on soil and wood mulch conditions at time zero), our estimates of external nutrient sources including those from wood mulch and atmospheric deposition suggest that PSF contributes to the nutrient accumulation observed in HPP. The environmental conditions at our sites, with an extremely poor soils and a boreal climate are also conducive to PSF. In fact, PSF is more commonly reported for extreme environments (Ehrenfeld et al., 2005).

More fungi-fungi, bacteria-bacteria and fungi-bacteria interactions were observed in the HPP than the LPP. The soil microbiome found in the HPP was dominated by Cortinarius taxa and Chitinophagaceae. Cortinarius is the largest genus of ectomycorrhizal fungi (> 2,000 species) (Høiland and Holst-Jensen, 2000; Soop and Gasparini, 2011), and represents approximately 10% of the species forming ectomycorrhizae (Rinaldi et al., 2008). Cortinarius species often dominate the ectomycorrhizosphere in reforested plots (Stefani et al., 2009; Pickles et al., 2010; Villarreal-Ruiz and Neri-Luna, 2017; Haas et al., 2018). Their extraradical mycelium is classified as a medium-distance exploratory type (Agerer, 2001; Agerer and Raidl, 2004), in the fringe subtype. This particular type of mycelium forms fans and rhizomorphs that ramify and interconnect repeatedly (Agerer, 2001). Furthermore, species within Cortinarius are considered to be nitrophobic taxa (Lilleskov et al., 2001; Lilleskov et al., 2002; Haas et al., 2018). Thus, the oligotrophic soil environment in borrow pits provides suitable conditions for the proliferation of Cortinarius. Interestingly, the roots of stunted pines in the LPP were colonized by several EM fungi and were showed a positive correlation with bacterial taxa different from those observed in the HPP. No fungal or bacterial pathogens were detected; therefore, they are unlikely to be responsible for the stunted trees that survived in the LPP. EM fungi were detected even in the absence of host trees in the NPP, indicating the constant presence of substantial inoculum in the surrounding area. Sampling was performed in late autumn, and the phenological timing of the Basidiomycota fruiting bodies in late summer may have influenced the spore load in bare soil samples, through wind transport and deposition (Després et al., 2012; Manninen et al., 2014).

The core taxa and the differentially abundant taxa suggest the presence of abundant resources in the HPP. For example, two bacterial taxa from the family Chitinophagaceae were significantly more abundant in the HPP and positively correlated with the abundance of Cortinarius. Some species in the family Chitinophagaceae are predicted to have β-glucosidase activity and are responsible for degrading polymeric organic matter including chitin and cellulose (Bailey et al., 2013). The inferred network for the core taxa showed that ASV0101 (Acidisoma sp.), ASV0087 (candidate alphaproteobacterial order Ellin329), OTU069 (Saccharomycetales), and OTU047 (Mortierella) were putative keystone taxa in the HPP. Sugars, polyalcohols and some organic acids are the preferred growth substrates for taxa within the genus Acidisoma (Belova et al., 2009). The presence of fungal yeasts (OTU069) reflects the availability of simple compounds such as sugars (monomers and oligomers), released by the hydrolytic enzyme activities of fungal decomposers or by root exudation. The decomposition of RCW and leaf litter releases polysaccharides and lignin which are likely to stimulate the development of these keystone species. The soil environment in the HPP was also more favorable for chemoheterotrophic bacterial taxa, that derive energy from decomposition of soil organic matter. This is consistent with the highest bacterial biomass values found in the HPP. The bacterial and fungal biomass recorded in these plots reflects nutrient assimilation and active biogeochemical transformation (Mason-Jones et al., 2022).