Study area
This study took place in western Quebec (Canada) on the traditional territories of the Abitibiwinni (Anishnaabe), Mistissini (Cree), and Nemaska (Cree) Indigenous communities (Fig. 1). The Abitibiwinni First Nation’s territory is located near the Ontario border, mainly in the black spruce – feather moss bioclimatic domain, except for its southern part located in the balsam fir – paper birch bioclimatic domain (Bélisle and Asselin 2021). The Mistissini Cree Nation’s territory is located further north and is very large, with its southern part in the black spruce – feather moss bioclimatic domain and its northern part in the spruce – lichen woodland bioclimatic domain. The territory of the Nemaska Cree Nation is located at the northern boundary of the black spruce – feather moss bioclimatic domain. This territory has experienced frequent fires in the past decades, which led to younger and more open stands (Eeyou Planning Commission 2017).
Fruit and soil sampling
Samples were collected on the territories of the Abitibiwinni, Mistissini and Nemaska communities on August 2nd, 16-17th, and 20th 2022, respectively. The difference in timing between territories was planned to collect fruits at the same phenological stage, accounting for latitude differences. All samples were collected in the black spruce – feather moss bioclimatic domain (Fig. 1), more specifically in black spruce stands, to limit the influence of potential confounding environmental variables. A total of 24 sites were sampled, that is 8 in each territory: 3 under hydroelectric power lines, 3 near a mine (200 m or less), and 2 control sites at least 1 km away from these disturbances. The selection of the distance around mining sites for analysis was informed by findings from prior studies, which determined that the effect of mining activities on vegetation extended to an average distance of 200 m (Boisvert et al. 2021; Yin et al. 2023a). All sites were located at least 100 m from roads to limit edge effects, trampling, and disturbances due to traffic.
At each site, four fruit samples were collected from different plant individuals for analysis of (poly)phenol compounds. Three organic soil samples were also collected with an auger in a triangular fashion around the fruit sampling locations, for analysis of soil microbiome and soil properties. The three soil samples were combined into a composite sample representative of the soil around the sampled plants. All samples (fruit and soil) were preserved in a cooler immediately after collection, then transferred to a -80°C freezer after each sampling day until extraction.
Extraction and analysis of (poly)phenol compounds
Fruit samples were sent to INAF’s chemical analysis laboratory (Institute of Nutrition and Functional Foods, Laval University, Quebec City, Canada) where they were extracted for analysis of flavonoids and total (poly)phenol compounds.
Flavonoids
The proanthocyanidins (PACs) content of the samples was measured by phloroglucinolysis followed by a UPLC-UV-MS/MS analysis. Phloroglucinolysis is a process in which PACs are cleaved into their base units, flavan-3-ols, using phloroglucinol (Kennedy and Jones 2001). This allows to quantify PACs and to determine their polymerization degree. A catechin standard was used to quantify PACs by UPLC-UV. The identification of detected compounds was then confirmed by triple quadrupole mass spectrometry.
The flavonol content of the samples was measured using the same procedure as for PACs, albeit with the omission of phloroglucinolysis and the use of a quercetin-3-glucoside standard for the purpose of quantification.
Total (poly)phenols
Total (poly)phenolic content was determined according to the Folin-Ciocalteu method as described in Dudonné et al. (2015) using gallic acid as a standard. Extract solutions (20 mL in 20% methanol 0.1% TFA) were mixed with 100 mL of 10-fold diluted Folin-Ciocalteu reagent and 80 mL of sodium carbonate solution (75 g/L). After 1 h of incubation at room temperature, the absorbance was measured at 765 nm using a BMG Labtech Fluostar Omega microplate reader (Offenburg, Germany).
Physicochemical analyses
Before proceeding to physicochemical analyses, soil samples were oven-dried at 35°C for 3 days, then passed through a 5.6 mm sieve to remove large debris, and stored in a -80°C freezer until analysis. Soil samples were analyzed for mineral contents (calcium, aluminum, potassium, phosphorus, magnesium, boron, copper, iron, manganese, zinc, sodium), pH, cation exchange capacity (CEC), organic matter content, as well as nitrogen and carbon content. Minerals were extracted with a Mehlich III solution and quantified by ICP-Optical Emission Spectrometry with the appropriate standard for each mineral (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques du Québec 2014). Soil pH was determined in water with a pH-meter (Ministère de l’Environnement, de la Lutte contre les Changements Climatiques, de la Faune et des Parcs 2023). The organic matter percentage in the soil was calculated by loss on ignition (Ministère du Développement durable, de l’Environnement et de la Lutte contre les changements climatiques du Québec 2017).
Microbiome analysis
DNA Extraction
The DNA of the soil samples was extracted from 250 mg of sampled soil using the Qiagen DNeasy Powersoil Pro kit (QIAGEN 2022) following the manufacturer’s protocol. Two negative extraction controls with no soil were also processed following the same protocol. The extracted DNA was stored immediately in a -80°C freezer until further processing. DNA samples were then sent to Genome Quebec Innovation Center (Montreal, Canada) for amplification and sequencing.
Amplification and sequencing
The metabarcoding method was used to detect and quantify bacteria and fungi in the organic soil. For bacteria, the extracted DNA samples were amplified using the primers 515b-FwR1 forward (GTGYCAGCMGCCGCGGTAA) and 926-RvR2 reverse (CCGYCAATTYMTTTRAGTTT) (Parada et al. 2016). For fungi, the primers used were ITS-9F forward (GAACGCAGCRAAIIGYGA) and ITS4R reverse (TCCTCCGCTTATTGATATGC) (White et al. 1990; Ihrmark et al. 2012). PCR amplifications were performed with 5 minutes of initial denaturation at 95°C, 34 cycles (bacteria) or 40 cycles (fungi) of 30 seconds at 94°C, 30 seconds at 50°C, and 1 minute at 72°C, then a final elongation step of 10 minutes at 72°C. Prior to amplification and sequencing, DNA quality was checked using 1% agarose gel electrophoresis.
Amplicons were sequenced on the Illumina MiSeq platform for paired-end reads. A negative control was included in the sequencing for both bacteria and fungi, to ensure that the extraction step did not result in contamination of the samples.
Bioinformatic workflow
The DADA2 R-package was used to build amplicon sequence variants (ASVs) from the raw sequences (Callahan et al. 2016). For bacteria, sequence primers were removed by trimming the first 19 nucleotides of forward reads and first 20 nucleotides of reverse reads in DADA2. After checking reads quality, 16S reads were also truncated at position 260 for forward reads and 190 for reverse reads, as there was a drop in read quality after these points. For fungi, ITS sequences primers were removed with cutadapt (Martin 2011) prior to assembly with DADA2. As the ITS sequence length is variable, ITS reads of lesser quality were not truncated to ensure that forward and reverse read could merge, but ITS read quality was decent overall. For 16S and ITS analyses, reads were pseudo-pooled during the ASVs assembly step in order to allow for the detection of rare ASVs. External contaminants were then removed using the decontam R-package using the prevalence method (Davis et al. 2018). Taxonomy was assigned using the Silva v138.1 database formatted for DADA2 for the bacterial sequences (McLaren and Callahan 2021) and the UNITE v.9.0 database for the fungal sequences (Abarenkov et al. 2022).
Once ASVs were built, data was transferred into a phyloseq object with the phyloseq R-package for handling (McMurdie and Holmes 2013). As a quality check, we removed ASVs that were found in only one sample, and those that were found less than 10 times across all samples. A phylogenetic heat tree of the taxa found in the samples was built using the metacoder R-package for visualization of taxonomic diversity across all samples (Foster et al. 2017). Before downstream analyses, the library of each sample was repeatedly rarefied to the library size of the smallest sample by drawing random ASVs without replacement with 1000 repetitions with the mirlyn R-package (Cameron et al. 2021), to control for biases in richness induced by differences in library sizes between samples. As the analyses available in the mirlyn package are limited, the multiple libraries created by the repeated rarefaction were then condensed into a single phyloseq object for compatibility with further packages. To do so, a table of ASV abundance was constructed by rounding the mean abundance of each ASV after 1000 rarefactions for each sample.
Statistical analyses
All analyses were performed using the R software version 4.3.1 (R Core Team 2023). Bacteria and fungi were treated separately for all analyses.
Effect of disturbances and territory on (poly)phenol concentrations
The effect of disturbance type, territory, and their interaction on (poly)phenol concentrations was evaluated with a PERMANOVA using the adonis2 function of the vegan R-package (Oksanen et al. 2022). The random effect of the sampling sites was accounted for by constraining permutations: sampling sites could be permuted, but not samples between sites. A total of 9999 permutations were performed using Euclidean distances.
Effect of disturbances and territory on soil properties and microbiome
The effect of disturbance type, territory, and their interaction on soil properties was also analyzed with a PERMANOVA using the adonis2 function from the vegan R-package (Oksanen et al. 2022). Variables were normalized by scaling prior to the PERMANOVA, and 9999 permutations were performed using Euclidean distances. To further study the response of soil properties to disturbances, an ANOVA followed by a Tukey test was performed.
Differences in alpha diversity between disturbance types and territories were plotted using the plot_richness function of the phyloseq R-package. Non-metric multidimensional scaling (NMDS) ordination plots were also produced, using the vegan package (Oksanen et al. 2022) on Hellinger transformed data to visualize the relation between phylum abundance, soil properties and sampling sites. Finally, a linear discriminant analysis effect size analysis (LEfSe) was performed to detect taxa that were differentially abundant between disturbance types (Segata et al. 2011). Since this analysis has a high false discovery rate (Nearing et al. 2022), the LDA threshold was conservatively set to 3.5. The LEfSe was computed in R using the microbiomeMarker package (Cao et al. 2022).
Effect of soil properties on microbiome
The effect of soil properties on soil microbiome was evaluated using a distance-based redundancy analysis (db-RDA) using the vegan R-package (Oksanen et al. 2022). The explanatory variables used in the model were the following: pH, contents of nitrogen, carbon, phosphorus, magnesium, potassium, iron, copper, manganese, zinc, aluminum, calcium, and sodium, and CEC. Explanatory variables were scaled prior to analysis, and the dissimilarity matrix of microbiome abundance was calculated using the Hellinger distance. Significance of the db-RDA model, axis, and terms was then evaluated by permutation, using the anova function of the vegan R-package with 9999 permutations (Borcard et al. 2018). P-values were adjusted with the Benjamini-Hochberg correction (Benjamini and Hochberg 1995).
Effect of microbiome on (poly)phenol concentrations
(Poly)phenol concentrations were not scaled prior to analyses, as concentration was the best available proxy for (poly)phenol bioavailability.
The effect of soil bacterial abundance on (poly)phenol concentrations was evaluated. For each site, the average concentration of each (poly)phenol compound was calculated. The rarefied ASV abundances were summed to the phylum level, and phyla with less than 500 observations across all samples were removed from the analysis, as they were unlikely to meaningfully affect plant (poly)phenols. Then, a redundancy analysis (RDA) was used to investigate the effect of the abundance of each phylum on (poly)phenol concentrations using the vegan R-package. Significance of the model, axis and explanatory variables was then evaluated by permutation, with the anova function of the vegan R-package using 9999 permutations (Borcard et al. 2018). P-values were adjusted with the Benjamini-Hochberg correction (Benjamini and Hochberg 1995). If a significative effect was found at the phylum level, another RDA was conducted with the abundance of the classes of this phylum, and so on with inferior taxonomic levels. In addition, the Proteobacteria and Firmicutes phyla were further explored as they contain plant growth promoting bacteria (Bulgarelli et al. 2013; Youseif 2018; Getahun et al. 2020). The same procedure was followed to evaluate the effect of fungi abundance on (poly)phenols, except with a threshold of 300 observations across all samples, as fungi were less abundant in general. In addition, the Ascomycota and Basidiomycota phyla were also explored as they contain ericoid mycorrhiza (Dong et al. 2022). FUNGuild v1.1 (accessed on august 3rd 2023) was used to get fungi putative functional assignations in order to help discuss the results (Nguyen et al. 2015).
When an RDA highlighted a taxon as having a significant effect on the (poly)phenol profile, the effect of this taxon on individual (poly)phenol was further evaluated by constructing linear models.