# Towards the development of multifunctional molecular indicators combining soil biogeochemical and microbiological variables to predict the ecological integrity of silvicultural practices.

Vincent Peck, Liliana Quiza, Jean-Philippe Buffet, Mondher Khdhiri, Audrey-Anne Durand, Alain Paquette, Nelson Thiffault, Christian Messier, Nadyre Beaulieu, Claude Guertin, et al.

► To cite this version:

Vincent Peck, Liliana Quiza, Jean-Philippe Buffet, Mondher Khdhiri, Audrey-Anne Durand, et al. Towards the development of multifunctional molecular indicators combining soil biogeochemical and microbiological variables to predict the ecological integrity of silvicultural practices. Microbial Biotechnology, 2016, 9 (3), pp.316-29. 10.1111/1751-7915.12348. pasteur-01351533

HAL Id: pasteur-01351533

https://riip.hal.science/pasteur-01351533

Submitted on 3 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Towards the development of multifunctional molecular indicators combining soil biogeochemical and microbiological variables to predict the ecological integrity of silvicultural practices

Vincent Peck,<fnref n="1" /> Liliana Quiza,<fnref n="1" /> Jean-Philippe Buffet,<fnref n="1" />

Mondher Khdhiri,<fnref n="1" /> Audrey-Anne Durand,<fnref n="1" />

Alain Paquette,<fnref n="2" /> Nelson Thiffault,<fnref n="2" />,<fnref n="3" />

Christian Messier,<fnref n="2" />,<fnref n="4" /> Nadyre Beaulieu,<fnref n="5" />

Claude Guertin<fnref n="1" />,* and Philippe Constant<fnref n="1" />,**

fication of the soil samples into four distinct groups categorized as a function of their deviation from baseline ecological conditions. According to this classification model, simple trenching was the approach that represented the lowest ecological risk potential at the microsite level. No relationship was observed between MSP method and soil bacterial community structure as assessed by high-throughput sequencing of bacterial 16S rRNA gene; however, indicator genotypes were identified for each multifunctional soil class. This is the first identification of multifunctional molecular indicators for baseline and disturbed ecological conditions in soil, demonstrating the potential of applied microbial ecology to guide silvicultural practices and ecological risk assessment.

## Summary

The impact of mechanical site preparation (MSP) on soil biogeochemical structure in young larch plantations was investigated. Soil samples were collected in replicated plots comprising simple trenching, double trenching, mounding and inverting site preparation. Unlogged natural mixed forest areas were used as a reference. Analysis of soil nutrients, abundance of bacteria and gas exchanges unveiled no significant difference among the plots. However, inverting site preparation resulted in higher variations of gas exchanges when compared with trenching, mounding and unlogged natural forest. A combination of the biological and physicochemical variables was used to define a multifunctional classification of the soil samples into four distinct groups categorized as a function of their deviation from baseline ecological conditions. According to this classification model, simple trenching was the approach that represented the lowest ecological risk potential at the microsite level. No relationship was observed between MSP method and soil bacterial community structure as assessed by high-throughput sequencing of bacterial 16S rRNA gene; however, indicator genotypes were identified for each multifunctional soil class. This is the first identification of multifunctional molecular indicators for baseline and disturbed ecological conditions in soil, demonstrating the potential of applied microbial ecology to guide silvicultural practices and ecological risk assessment.

## Introduction

Tree plantations are gaining increased interest to protect natural forests, restore ecosystem services and meet various social needs. For instance, more than 15% of current wood production is supported by tree plantations and this contribution is expected to rise in the future (Paquette and Messier, 2009). Mechanical site preparation (MSP) is a common practice to improve seedling performance in tree plantations dedicated to intensive wood production and forest regeneration. This approach, applied following clear-cut or variable retention logging, involves the utilization of heavy machineries to break soil structure in order to improve soil physical conditions limiting tree growth and control competing vegetation (Löf et al., 2012). Soil scarification and mounding are the most usual MSP methods. Scarification consists in mixing organic and mineral upper layers of soil by trenching. This technique increases soil aeration and temperature, while favouring nutrient availability by accelerating nitrogen mineralization and limiting the invasion of competing vegetation (Prévost, 1992; Thiffault and Jobidon, 2006). Mounding consists in creating elevated planting spots where soil is excavated and deposited on the ground, next to the ditch, while inverting the soil horizons to get

Received 1 September, 2015; accepted 15 December, 2015. For correspondence. *E-mail Claude.Guertin@iaf.inrs.ca; Tel. (450) 687-5010 ext. 4234; Fax (450) 686-5501. **E-mail PhilippeConstant@iaf.inrs.ca; Tel. (450) 687-5010 ext. 4117; Fax (450) 686-5501.

Microbial Biotechnology (2016) 9(3), 316–329
doi:10.1111/1751-7915.12348

**Funding Information** This research was supported by a grant from the Natural Sciences and Engineering Research Council of Canada (NSERC) Engage Program to PC (grant EGP 463612-14).

<footnote n="1">¹INRS-Institut Armand-Frappier, 531 boulevard des Prairies, Laval, Québec Canada H7V 1B7.</footnote>

<footnote n="2">²Centre d'étude de la forêt, Université du Québec à Montréal, Case postale 8888, succursale Centre-ville, Montréal, Québec Canada H3C 3P8.</footnote>

<footnote n="3">³Direction de la recherche forestière, Ministère des Forêts, de la Faune et des Parcs, 2700 Einstein, Québec, Québec Canada G1P 3W8.</footnote>

<footnote n="4">⁴Institut des Sciences de la Forêt Tempérée (ISFORT), Université du Québec en Outaouais (UQO), 58 rue Principale, Ripon, Québec Canada J0V 1V0.</footnote>

<footnote n="5">⁵Produits Forestiers Résolu, 2419 Route 155 sud, La Tuque, Québec Canada G9X 3N8.</footnote>
mineral soil on the top and an organic layer at the bot-
tom (Sutton, 1993). As an alternative to mounding to
increase worker safety for following operations such as
cleaning and thinning, excavated soil can also be placed
back into the original ditch, resulting in inverting site
preparation. Mounding treatments provide the beneficial
effects on soil aeration, temperature, competing vegeta-
tion control and nutrient availability observed with scarifi-
cation and are particularly suitable to get seedling spots
free from water logging conditions in wet areas (Kabrick
et al., 2005; Simon et al., 2013). Mounding and inversion
using an excavator are also expected to be less damag-
ing to the environment than regular trenching, disturbing
a lesser proportion of the treated area (especially inver-
sion) and avoiding the creation of linear trenches that
may cause soil erosion (Buitrago et al., 2015).

Increased soil temperature, followed by the stimulation of organic carbon mineralization and nutrient release are the prevailing soil disturbances caused by MSP (Jandl et al., 2007). In contrast to plant diversity, which is expected to be resilient to this management practice (Haeussler et al., 2004), the resilience, resistance or vulnerability of soil microbiome to these marked soil disturbances remains unknown. This question is of critical importance in ecosystem-based forest management where establishment of tree plantations must exert minimal alteration or even restore ecosystem services sustained in natural forests (Martin et al., 2014). Because microorganisms play a crucial role in global biogeochemical cycles and closely interact with vegetation through nutrient transfers as well as water retention, the composition of soil microbial communities can be seen as an indicator of soil ecosystem functioning. For instance, high microbial biodiversity was shown to promote soil multiple ecological functions including plant diversity and nutrient cycling (Wagg et al., 2014), resistance to the invasion and survival of allochthonous pathogen species (van Elsas et al., 2012) as well as resistance and resilience of microbial community structure to certain environmental stress (Tardy et al., 2014).

The overarching objective of this study was to com-
pare the metabolic activity and composition of soil bacte-
rial community between recently clear-cut and site
prepared plots within an hybrid larch (Larix × marschlinsii
Coaz) plantation and unlogged natural mixed boreal
forest conservation areas to assess the impact of differ-
ent MSP on soil ecosystem functioning a few years
following treatments. The rationale for this approach is
the identification of soil biogeochemical processes and
microbes affected by MSP to select silvicultural practices
offering the best early environmental performance at the
microsite level. We tested the hypothesis that both the
conversion of natural forest to a hybrid larch monoculture
and the intensification of MSP treatments reduce the

activity and the taxonomic diversity of soil bacterial
community in the early years after this conversion. We
finally combined soil biogeochemical and microbiological
data sets to explore the relevance of microbial molecular
indicators to predict the environmental impact and guide
silviculture practices. In contrast to previous investigations
comparing single parameter such as soil enzyme activi-
ties between reference sites and managed forests (e.g.
Staddon et al., 1999), we have combined multiple biogeo-
chemical variables as a metric to classify soil under an
original multifunctional system and searched for indicator
16S rRNA gene sequences restricted to specific soil mul-
tifunctional classes. Biogeochemical variables were
selected to get a broad classification system, including
gaseous exchange involving taxonomically diverse and
specialized microbial guilds. Under this framework, we
expected to identify indicators for baseline or disturbed
ecological conditions in soil and to explore the potential of
soil bacterial community monitoring as a promising
approach to predict the ecological integrity of soil in the
early stage of intensively managed tree plantations.

Results

Soil physicochemical properties, gaseous exchanges
and abundance of bacteria

Conversion of unlogged natural forest to a hybrid larch
plantation caused no significant change in measured soil
physicochemical properties at the microsite level
(Table 1). Soil carbon and nitrogen concentrations were
significantly and positively correlated (Pearson,
*P* < 0.0001) and C:N stoichiometry ranged from 18 to 32
among samples. All tested soil represented a net sink
for H₂ and CO and net source for CO₂ (Table 1). Even
though gaseous exchanges were not affected by MSP
procedures, some treatments induced high variations
between replicates. In general, the coefficient of variation
(CV) of gas exchanges increased as a function of the
intensification of the MSP techniques (Fig. 1). For
instance, H₂, CO and CO₂ exchanges measured in soil
samples collected in excavated mounds displayed CV
ranging from 58% to 62% while the same variables
measured in unlogged natural forest samples were char-
acterized by CV of 7–32%. Trace gas turnover showed
no significant relationship with soil carbon and nitrogen
concentrations (Pearson, *P* > 0.05) but a negative corre-
lation was observed between soil respiration and pH
(Pearson, *P* = 0.03). H₂ uptake rates increased as
function of CO uptake (Pearson, *P* < 0.0001) and CO₂
production (Pearson, *P* = 0.02) activities. CO₂ production
rates were not related to measured CO uptake activities
(Pearson, *P* > 0.05). The abundance of bacteria was
proportional to CO₂ production rate in soil (Pearson,
*P* = 0.01).
<table><caption>**Table 1.** Physicochemical properties, trace gas exchanges, bacterial 16S rRNA gene abundance and bacterial diversity in soil. Averages are represented with standard deviations in parentheses<sup>a</sup>.</caption><thead><tr><th>Sites</th><th>Carbon (%)</th><th>Sand (%)</th><th>Silt (%)</th><th>Clay (%)</th><th>Nitrogen (%)</th><th>pH</th><th>H<sub>2</sub> oxidation (nmol g<sup>-1</sup> h<sup>-1</sup>)</th><th>CO oxidation (nmol g<sup>-1</sup> h<sup>-1</sup>)</th><th>CO<sub>2</sub> production (µmol g<sup>-1</sup> h<sup>-1</sup>)</th><th>16S rRNA gene (copies g<sup>-1</sup>)</th><th>Alpha diversity (Shannon)</th><th>Alpha diversity (ACE)</th></tr></thead><tbody><tr><td>Unlogged natural</td><td>4.9 (1.1)</td><td>83 (3)</td><td>13 (3)</td><td>4 (1)</td><td>0.22 (0.04)</td><td>4.3 (0.1)</td><td>6.3 (0.4)</td><td>4.5 (0.6)</td><td>4.1 (1.3)</td><td>7.9 (1.3) × 10<sup>9</sup></td><td>6.2 (0.08)</td><td>2950 (242)</td></tr><tr><td>Simple</td><td>5.4 (1.3)</td><td>86 (3)</td><td>13 (3)</td><td>1 (1)</td><td>0.27 (0.06)</td><td>4.6 (0.2)</td><td>5.8 (1.6)</td><td>4.5 (2.0)</td><td>3.1 (1.4)</td><td>4.9 (3.6) × 10<sup>9</sup></td><td>6.3 (0.2)</td><td>2895 (565)</td></tr><tr><td>Double intensive</td><td>4.1 (1.2)</td><td>84 (2)</td><td>14 (2)</td><td>2 (0)</td><td>0.21 (0.06)</td><td>4.7 (0.2)</td><td>5.2 (1.4)</td><td>3.8 (1.7)</td><td>1.5 (0.5)</td><td>2.0 (2.6) × 10<sup>9</sup></td><td>6.2 (0.3)</td><td>2617 (732)</td></tr><tr><td>Inversion</td><td>3.0 (1.4)</td><td>76 (2)</td><td>22 (2)</td><td>2 (0)</td><td>0.13 (0.04)</td><td>4.7 (0.2)</td><td>3.0 (2.3)</td><td>2.1 (1.4)</td><td>1.4 (0.2)</td><td>1.5 (1.4) × 10<sup>9</sup></td><td>6.1 (0.5)</td><td>2268 (985)</td></tr><tr><td>Mound</td><td>4.3 (0.8)</td><td>85 (8)</td><td>14 (8)</td><td>1 (1)</td><td>0.17 (0.03)</td><td>4.6 (0.2)</td><td>6.6 (4.1)</td><td>4.8 (1.8)</td><td>1.1 (2.3)</td><td>1.8 (2.1) × 10<sup>9</sup></td><td>6.0 (0.2)</td><td>2247 (877)</td></tr></tbody></table>

a. No variable showed significant difference between the treatments (ANOVA, *P* > 0.05).

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Coefficient of variation (CV) for gaseous exchanges. The bars represent the average and standard deviation of CV measured for H₂, CO and CO₂.</figcaption></figure>

## Multifunctional soil classification

Soil physicochemical properties, gaseous exchanges and bacteria abundance variables were utilized to define a multifunctional soil classification. Under this classification approach, the distribution of each variable was considered to compute a distance matrix measuring the association between soil samples (Fig. 2A). Four different multifunctional classes were identified in the clustering analysis. The level of disturbance characterizing each class was defined on the basis of their Euclidean distance from the unlogged natural forest plots. Firstly, the soil sample constituting class I (M-B) represented the MSP treatment that resulted in the most intense disturbance of baseline ecological functions. Secondly, soil samples included in class II (S-A, S-B) belong to the category of MSP treatments that caused slight deviations from baseline soil ecological functions. Soil samples belonging to class III (N-B, N-A, N-C) correspond to baseline of soil ecological functions. Finally, class IV (S-C, I-A, I-C, I-B, M-A, M-C, D-A, D-B and D-C) encompass soils for which the MSP treatments caused important alteration of baseline ecological functions at the microsite level. The unlogged natural forest was the sole condition for which replicated composite samples exhibited treatment-specific, conserved multifunctional profile. Taken together, this classification model indicates that soil multifunctional profile observed in inversion and mound excavations plots were those showing the strongest deviation from unlogged natural forest at the microsite level. With the exception of one replicate (S-C), simple trenching (S-A and S-B) was the treatment exerting the lowest incidence on soil multifunctional profile (Fig. 2A).

A principal component analysis was computed to represent the position of sampling sites in a reduced space defined by the measured variables. The ordination space
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

defined by the first two components explained 56% of
the variation observed (Fig. 2B). Five variables defined
the reduced space represented by both axes (Fig. 2C).
The gradient associated to CO₂ production rate
explained the distribution of plots along the first compo-
nent, while abundance of bacteria and H₂ uptake rate
defined the distribution along the second axis. Soil CO₂
respiration was the most preponderant variable responsi-
ble for the clustering of S-A and S-B trenching plots with
unlogged natural forest samples and the clustering of S-
C with mounding, inverting and intensive trenching plots.

**Fig. 2.** Multifunctional comparison of soil samples according to their physicochemical (C, N, C:N ratio, pH, texture), trace gas exchanges (H₂, CO, CO₂) and abundance of bacteria (16S rRNA gene abundance; labelled as 'Biomass' in the plot) profiles. (A) UPGMA agglomerative clustering of soil samples according to a Euclidean distance matrix calculated with standardized variables. The grey circles denote the nodes delineating the four multifunctional classes significantly discriminated by SIMPROF permutation procedure (*P* < 0.05). Colour labels show the assignation of the soil samples to their multifunctional class (red; class I, green; class II, blue; class II and black; class IV). The scale bar represents the Euclidean distance in the dendrogram. (B) Principal component analysis showing the distribution of sampling sites in a reduced space defined by soil physicochemical properties, gaseous exchanges and abundance of bacteria. The colours used to present soil samples correspond to the clusters identified in the UPGMA (A). (C) Variables defining the distribution of soil samples along the first and the second axis are represented along with the equilibrium circle of descriptors showing the contribution of variables to the formation of the reduced space. The detection limit of the qPCR assay was utilized to estimate the abundance bacteria in sample M-A for which the low yield of the DNA extraction procedure precluded qPCR and bacterial 16S rRNA gene profiling (see the Material and Methods section for more details).

Indeed, CO₂ efflux was 3.9 µmol g(dw)⁻¹ h⁻¹ in S-A and S-B, while a value of 1.6 µmol g(dw)⁻¹ h⁻¹ was measured in S-C. The sample M-B displayed the highest H₂ and CO uptake rates, explaining its position in the ordination space (Fig. 2B).

Soil bacterial community taxonomic structure

Quality control, classification and equalization of the 16S
rRNA gene sequence libraries yielded 5451 bacterial
OTU (97% identity threshold). Overall, 50% of the
sequences belonged to *Proteobacteria* mostly repre-
sented by *alpha*- (74%) and *delta-Proteobacteria* (11%).
The *Acidobacteria* and *Actinobacteria* were the two other
phyla dominating the bacterial communities with 19%
and 9% relative abundance respectively (Fig. S1). Nei-
ther the conversion of unlogged natural mixed forest to a
hybrid larch monoculture nor MSP treatments caused
significant alteration at the microsite level in bacterial
OTU richness as evaluated by the Shannon diversity
index and ACE estimator (Table 1). In general, beta
diversity defined as the variability in OTU composition
among replicated plots measured by multivariate disper-
sion showed more variations in the larch plantation than
in unlogged natural mixed forest conservation areas
(Fig. S2). Agglomerative clustering of the samples
according to their microbial community profile showed
that soil samples collected in different treatments could
not be discriminated on the basis of OTU composition
(Fig. 3A). Unlogged natural forest clustered together with
M-B, while all other clusters were composed either of
unique plots or plots originating from different MSP pro-
cedures. A redundancy analysis (RDA) was performed
to infer the relationship between 16S rRNA gene profiles
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Comparison of soil samples according to their 16S rRNA gene profile. (A) UPGMA agglomerative clustering of soil samples derived from a matrix of Euclidean distance calculated after Hellinger transformation of OTU (97% identity threshold) absolute abundance. The grey circles denote the nodes delineating the four groups of samples significantly discriminated by SIMPROF permutation procedure (*P* < 0.05). The scale bar represents the Euclidean distance in the dendrogram. Colour labels show the assignation of the soil samples to their multifunctional class (red; class I, green; class II, blue; class II and black; class IV). (B) Parsimonious RDA triplot of Hellinger-transformed OTU absolute frequency matrix explained by soil pH and C:N ratio. Only the 14 OTUs displaying extreme distribution in the reduced space are depicted for clarity. These OTUs are identified in the legend with colour bars discriminating α-Proteobacteria (black), β-Proteobacteria (blue), δ-Proteobacteria (red) and other phyla (orange), as determined using the Greengene reference database V13_8_99 (McDonald et al., 2012). The colour labels used to present soil samples in the RDA triplot correspond to the clusters identified in the UPGMA (Fig. 4A). The sample M-A is absent due the low yield of the DNA extraction procedure for this soil (see the Material and Methods section for more details).</figcaption></figure>

and environmental variables (Fig. 3B). The most parsimonious model to explain variation in the distribution of 16S rRNA gene sequences included soil C:N stoichiometry and pH. The other variables being redundant to soil C:N and pH, their addition to the analysis increased the

variance inflation factor unduly, and they were therefore ignored in the analysis. The first two canonical axes explained 49% of the total variance of bacterial OTU distribution. Significance of the RDA was confirmed with 1000 permutations of OTU data matrix (*P* = 0.001). Soil pH played an important role in the dispersion of the samples along the first axis, while C:N discriminated the samples along the second.

### Microbial molecular indicator for soil multifunctional classes

The bacterial community profiles derived from the soil samples were classified within the four classes defined by the clustering of soil according to the multifunctional classification model (Fig. 2A). In total, 693 OTU were ubiquitously detected in all soil samples (Fig. S3). However, the search for indicator OTU unveiled coherence between the distribution of several members of the rare biosphere comprising less than 0.1% of the bacterial communities and soil multifunctional classes (Fig. 4). Indicator OTU encompassed a broad taxonomic diversity and *Proteobacteria* was the only phylum represented in the four multifunctional classes. In the case of multifunctional classes represented by more than one indicator, the OTU displaying the highest abundance was identified as representative indicator. Two indicators are of particular interest for this study. Firstly, the OTU 3283 (classified as a member of the order *Myxococcales*) was considered as a disturbance indicator because it was detected in soil samples for which the multifunctional classification (class IV) diverged from baseline ecological functions observed in unlogged natural mixed forest (Fig. 4). Furthermore, distribution of this candidate disturbance indicator was consistent with the environmental variables defining the multifunctional soil classification model. According to the unconstrained principal component analysis (PCA) ordination analysis, soil CO₂ respiration was a preponderant variable discriminating soil samples classified in class IV from those belonging to class II and class III (Fig. 2). As expected, the abundance of 16S rRNA gene sequences belonging to OTU 3283 showed a negative correlation with CO₂ soil respiration (Table 2). Secondly, the OTU 398 (classified as a member of the order *Rhodospirillales*) was considered as an indicator for baseline ecological conditions in soil because it was only detected in unlogged natural mixed forest conservation areas (Fig. 4). Distribution of this OTU was related to soil pH, bacteria abundance and soil clay content (Table 2). This result is in agreement with the contribution of these variables to the distribution of unlogged natural mixed forest soil samples in the PCA utilized to identify the factors defining the structure of the multifunctional classification model (Fig. 2B and C). In
<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Identification of soil multifunctional molecular indicators. The heatmap shows the absolute abundance of OTUs detected in the soil samples categorized into the four multifunctional classes previously defined (Fig. 2A). Colour bars show the assignation of the soil samples to their multifunctional class (red; class I, green; class II, blue; class II and black; class IV). Taxonomic assignment of the OTU was done using the Greengene reference database V13_8_99 (McDonald *et al.*, 2012). Representative indicators are highlighted with bold characters and are identified with an arrow and asterisk (← *).</figcaption></figure>
<table><caption>**Table 2.** Spearman correlation between the abundance of the representative indicator ribotype selected for the multifunctional classes identified in this study and soil biological, physical and chemical variables (see Table 1).</caption><thead><tr><th rowspan="2">Variables</th><th colspan="4">Selected multifunctional indicators</th></tr><tr><th>Class I<br>OTU 2465 (Legionellales)</th><th>Class II<br>OTU 1984 (Opitutales)</th><th>Class II<br>OTU 398 (Rhodospirillales)</th><th>Class IV<br>OTU 3283 (Myxococcales)</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>0.34 (0.24)</td><td>0.03 (0.92)</td><td>0.25 (0.38)</td><td><strong>-0.67 (0.009)</strong></td></tr><tr><td>CO</td><td>0.05 (0.88)</td><td>-0.41 (0.15)</td><td>0.01 (0.98)</td><td>-0.53 (0.053)</td></tr><tr><td>CO<sub>2</sub></td><td>0.14 (0.64)</td><td>0.16 (0.60)</td><td>0.49 (0.08)</td><td><strong>-0.77 (0.001)</strong></td></tr><tr><td>Bacteria abundance</td><td>-0.11 (0.70)</td><td>0.30 (0.29)</td><td><strong>0.68 (0.007)</strong></td><td>-0.46 (0.101)</td></tr><tr><td>C</td><td>-0.46 (0.10)</td><td>0.22 (0.44)</td><td>0.21 (0.47)</td><td><strong>-0.59 (0.026)</strong></td></tr><tr><td>N</td><td>-0.36 (0.21)</td><td>0.30 (0.29)</td><td>0.21 (0.47)</td><td><strong>-0.62 (0.017)</strong></td></tr><tr><td>C/N</td><td>-0.46 (0.10)</td><td>-0.45 (0.11)</td><td>0.40 (0.16)</td><td>-0.10 (0.74)</td></tr><tr><td>pH</td><td>0.13 (0.65)</td><td>-0.02 (0.95)</td><td><strong>-0.73 (0.003)</strong></td><td><strong>0.61 (0.020)</strong></td></tr><tr><td>Sand</td><td>0.43 (0.12)</td><td>0.44 (0.11)</td><td>0.06 (0.84)</td><td><strong>-0.72 (0.003)</strong></td></tr><tr><td>Clay</td><td>-0.07 (0.80)</td><td>-0.44 (0.12)</td><td><strong>0.81 (0.0005)</strong></td><td>-0.13 (0.65)</td></tr><tr><td>Silt</td><td>-0.39 (0.16)</td><td>-0.38 (0.19)</td><td>-0.21 (0.47)</td><td><strong>0.77 (0.001)</strong></td></tr></tbody></table>

Spearman rho correlation coefficients are presented with the significance levels (α-value) in brackets. Bold characters represent significant correlation.

contrast to these two indicators, OTU for which the distribution was specific to multifunctional class I (OTU 2465; *Legionellales*) or class II (OTU 1984; *Opitutales*) did not show any significant correlation with the variables measured in this study (Table 2).

Two PCR assays were designed to challenge the indicator for baseline multifunctional conditions (OTU 398; class III) and the disturbance indicator (OTU 3283; class IV) identified through the indicator species statistical analysis (Table S1). Baseline indicator was detected in the three samples comprising class III (N-A, N-B, N-C), but weak PCR amplification signal also was detected for sample M-B (class I) and sample I-A (class IV), where no reads belonging to OTU 398 where retrieved from the high-throughput 16S rRNA gene sequencing analysis (Fig. S4A). On the other hand, disturbance indicator was detected in three out of the five samples from which reads assigned to OTU 3283 were detected in the sequencing analysis (Fig. S4B).

## Discussion

Monitoring of early growth of hybrid larch (*Larix × marschlinsii* Coaz) seedlings over two growing seasons following planting demonstrated that trenching (simple and double), mounding and inverting site preparation resulted in undistinguishable growth performance (Buitrago et al., 2015). This observation led to the conclusion that simple trenching represents the most economically attractive silvicultural practice for early growth of hybrid larch. Because soil microorganisms are at the core of key ecological functions such as nutrients transfer and biogeochemical cycles, we investigated the impact of MSP on the metabolic activity and composition of soil microbial communities. Our hypothesis that both the conversion of natural forest to a hybrid larch

monoculture and the intensification of MSP treatments reduce the activity and the taxonomic diversity of soil microbiome in the early years of this conversion was not verified. Neither the microbial activities nor the diversity indices showed significant difference between the MSP techniques and natural unlogged forest. Beside the examination of environmental variables and the distribution of microorganisms in soil, the originality of our approach was the development of a multifunctional soil classification system along with molecular indicators to assess the potential ecological risk of MSP techniques on ecosystem functions. The integration of biogeochemical and microbiological variables – instead of considering variables individually – was shown to be essential to determine the response of soil bacterial community to silvicultural practices.

An important question in conservation biology is the number of species that must be protected to ensure ecosystem functioning. A milestone in the field of biodiversity and ecosystem functioning was the observation that the minimal number of species required to support multiple functions is much higher than estimates derived from the analysis of individual processes (Hector and Bagchi, 2007; Gamfeldt et al., 2008). As ecosystems are conserved for their multiple natural services, these studies highlighted the need for integrated approaches to set reasonable conservation objectives (Balvanera et al., 2014). In this study, we selected multiple gas turnover as a subset of ecosystem functions. H₂ and CO soil uptake activities are catalysed by different guilds of bacterial species that contribute to mitigate the global emissions of these climate-relevant trace gases in the atmosphere. High-affinity, H₂-oxidizing bacteria are mostly represented by specialized *Actinobacteria* and some representatives of the *Proteobacteria*, *Chloroflexi* and *Acidobacteria* (Constant et al., 2011; Meredith et al.,
2013). These bacteria are responsible for 80% of the global sink of atmospheric H₂ (Constant *et al.*, 2009; Ehhalt and Rohrer, 2009). Soil survey showed the importance of carbon and nitrogen contents to explain the activity of this functional group in soil (Gödde *et al.*, 2000; Khdhiri *et al.*, 2015). Carboxydovore bacteria, encompassing the *Actinobacteria*, *Proteobacteria* and *Chloroflexi*, are responsible for the soil uptake of atmospheric CO (King and Weber, 2007; King and King, 2014; Quiza *et al.*, 2014). Even though no systematic investigation about the impact of soil nutrients on CO uptake rate has been reported, organic matter-rich soils displayed higher potential CO uptake activity than desert soils along vegetation transects (Weber and King, 2010). The activity of carboxydovore bacteria contributes to 15% of the global sink of atmospheric CO, while scavenging CO produced by biological and abiotic reactions in soil. In contrast to H₂ and CO, soil CO₂ respiration is supported by a broad diversity of bacteria and fungi thriving in soil using a heterotrophic growth metabolism. In addition to the availability of soil organic matter, factors influencing soil respiration rate are complex, including soil nutrient stoichiometry (Drake *et al.*, 2013). Including H₂, CO and CO₂ exchanges in multifunctional classification model was, therefore, a relevant approach to consider important ecosystem services that are provided by taxonomically diverse microorganisms. Addition of further environmental variables such as pH, C, N, abundance of bacteria and soil texture defining soil biogeochemical structure provided a supplementary dimension to the multifunctional classification system. Consideration of other environmental variables such as organic matter and nutrient turnover in the future is not precluded as this could refine the soil multifunctional classification system to satisfy other ecosystem functions prioritized in sustainable forest management.

In contrast to the comparative analysis of individual environmental parameters that unveiled no significant impact of MSP due to the variance of the observations, the multifunctional classification approach led to the identification of MSP treatment plots for which ecological functions differed from baseline conditions as measured in unlogged natural forest nearby. According to this classification, double trenching, mounding and inverting site preparation were the less sustainable MSP techniques at the microsite level. The low environmental performance of mounding and inversion plots at this scale was also supported by the higher coefficients of variation in the gaseous exchanges measured in these soils than those observed in trenching and unlogged natural mixed forest (Fig. 1). This finding, indicating that three distinct ecological services were sensitive to mounding and inverting, is probably due to the heterogeneity caused by the heavy machinery used for site preparation. The

procedure involved the transfer of excavated soil and was prone to large variations due to the presence of rocks, resulting in variable volume and depth of the excavations in addition to potential variation in the horizontal and lateral distribution of displaced organic and mineral soil horizons and ground microtopography. Our results are relevant at the microsite level, i.e. local environmental conditions directly influencing planted seedling physiology and growth. Although mounding and inverting intensively disturb the soil profile, thus having significant local impacts on soil functions, trenching furrows affect a higher proportion of the planted site. Determining which treatment is best in terms of overall ecological impacts will require further work, including scaling up and interpreting these impacts at the stand level, taking into account the proportion of disturbed soil in each treatment.

Soil microbial diversity was neither related to MSP treatment nor multifunctional classes identified in this study. Soil pH, rather than MSP treatments, was the most important parameter to explain the composition of soil microbial communities (Fig. 3). Indeed, variation in pH was shown as a preponderant variable explaining the composition of soil microbiome and any silvicultural treatment resulting in marked pH alteration are expected to influence the distribution of dominant lineages (Lauber *et al.*, 2009). The most important aspect of soil microbial community analysis was the clustering of the 16S rRNA gene profiles according to the multifunctional soil classification model, resulting in the identification of indicator OTU. Those indicators represented members of the rare biosphere in soil, for which the incidence of environmental variables on their distribution could not be observed using ordination techniques parameterized with the whole microbial community profiles (Fig. 3). Two relevant indicators were identified. The OTU 3283 (*Myxococcales*) detected in soil samples that were divergent from the baseline conditions observed in unlogged natural mixed forests could be used as a diagnostic tool to assess the impact of MSP on ecosystem services. On the other hand, the OTU 398 (*Rhodospirillales*) was associated to baseline soil multifunctional attributes. Bioindicators of soil quality involving exoenzyme activity measurements (Staddon *et al.*, 1999) and other approaches involving the monitoring of arthropods (Pearce and Venier, 2006) and small mammals (Pearce and Venier, 2005) as bioindicators of sustainable forest management have already been developed. In contrast to these metrics requiring relatively laborious monitoring efforts, molecular diagnostic tools, specific to selected indicator OTU, are fast and easily integrated with abiotic and biological factors in the environment. In the context of the hybrid larch plantation investigated here, monitoring of the identified indicators, especially the indicator for
baseline multifunctional soil attributes, could be used for environmental certification and ecological risk assessment after adequate scaling up addressing key requirements of ecological indicators.

The molecular indicators identified in this study are relevant at the stand level, with the abiotic and biotic conditions that prevailed in soil at the time of sampling. The experimental plan was not designed to address the consistency of molecular indicators over time, space and sites. Further investigations are needed to challenge the indicators in other hybrid larch as well as other fast-growing tree plantations and assess their specificity (Bartell, 2006). Longer term studies are also needed to assess whether spatial and temporal variations of soil biogeochemical processes in intensive silviculture can be related to alterations of molecular indicator distribution profiles. Indeed, it is expected that indicator OTU identified in this study are not universal due to the impact of abiotic and biological gradients shaping soil microbial communities in soil. Nevertheless, our data show that the distribution of soil microbial communities is not random, with the distribution of some members restricted to specific multifunctional soil classes, at the stand level. Combination of soil bacterial, fungal and archaeal community taxonomic and functional profiles through metagenomic surveys should be considered in future attempts to develop and validate molecular indicators for sustainable forest management. Finally, PCR diagnostic assays targeting the baseline and stress indicators clearly demonstrated that a composite of indicator OTU rather than a single indicator would be needed for predicting soil multifunctional classes (Fig. S4). This could be achieved through the adaptation of the indicator species modelling approach used to select the minimal subset of species necessary to predict species richness of invertebrate and vertebrate assemblages to microbial ecology (Thomson *et al.*, 2007; Azeria *et al.*, 2009).

In conclusion, we showed that MSP treatments influence the overall signature of soil biogeochemical structure at the microsite level, suggesting that mounding and inverting site preparation could pose higher potential ecological risk for ecological functions in soil than trenching at the local scale. Further work is needed to scale up and interpret these results at the landscape level and over longer period of time, taking into account the proportion of disturbed soil in each treatment and the overall productivity of each treatment. In effect, the lower ecological risk potential of a less productive treatment could be cancelled by the need of converting a larger area into fast-growing plantations to make up for the lesser productivity. Chronosequences, spatial variability, and different tree plantation types will need further attention in future investigations to elucidate the short term impact of MSP on soil biogeochemical structure and

challenge the molecular indicators identified in this study. Because of functional redundancy in microorganisms, these studies would be essential to challenge and define limitations of the indicators. Nevertheless, this work demonstrates the relevance of applied ecology to evaluate the sustainability of silvicultural practices.

## Experimental Procedures

### Sampling site

The study site was located near La Tuque (Québec, Canada; 47° 37' 19" N, 72° 49' 55" W), about 250 km north of Montréal. The experimental area, dominated by balsam fir (*Abies balsamea* (L.) Mill.), paper birch (*Betula papyrifera* Marsh.), yellow birch (*B. alleghaniensis* Britt.), red maple (*Acer rubrum* L.) and black spruce (*Picea mariana* (Mill.) BSP) was harvested in 2009 (clear-cut with 5% retention) prior to the installation of an experiment to investigate the impact of different MSP techniques on the growth performance of hybrid larch planted in April 2010 (Buitrago *et al.*, 2015). The experimental design consisted of a complete block design comprising three replicated blocks (Fig. 5). Briefly, each block was separated into four plots randomly assigned to different MSP treatments encompassing trenching and mounding. Trenching was performed with rotary discs mixing mineral and organic soil horizons. This treatment was applied either as a simple (simple plots; composite samples S-A, S-B, S-C) or double trenching (double intensive plots; composite samples D-A, D-B, D-C) corresponding to an increasing gradient of soil mixing (Buitrago *et al.*, 2015). An excavator was used for the mounding treatments consisting in inverting the soil horizons to place mineral soil on the top and the organic layer at the bottom. Excavated soil was either placed back to the original ditch (inversion plots; composite samples I-A, I-B, I-C) or on soil surface next to the ditch generated by the excavation (mound plots; composite samples M-A, M-B, M-C). Seedlings were planted in the hinge position of the trenching furrows or near the centre of excavated mounds. Finally, plots consisting of non-harvested, non-planted natural mixed forest areas (retention areas located within the three experimental blocks, and consisting in approximately 500–600 m² plots) were used as reference in this study (unlogged natural plots; composite samples N-A, N-B, N-C). These unlogged forest soils were selected as reference plots based on third-party forest certification criteria (e.g. Forest Stewardship Council), where reference forests in the vicinity of forest management units are used to measure the impact of management plans on ecosystem integrity.

The experimental area was visited in July 2014 for soil sample collection. In total, 10 replicated soil samples were collected in each plot (3 blocks × 5 plots × 10
<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Schematic representation of the sampling design. Five treatment plots were replicated in three blocks. Soil samples were collected at 10 locations from each plots (illustrated with grey trees). Soil was collected in proximity of the stem, at one of the four cardinal points. Replicated soil samples were pooled to obtain one composite sample per replicated plot. The (native) treatment consisted of non-harvested, non-planted natural mixed forest areas located within the experimental blocks.</figcaption></figure>

replicates = 150 samples in total). Individual trees for which information regarding early growth parameters was available (Buitrago et al., 2015) were first randomly selected and the A-horizon (0–10 cm) was collected within a 15 cm radius of the stem at one of the four cardinal points also determined randomly (Fig. 5). Each replicate sample consisted of ~200 g soil stored on ice in Whirl-Pak® bags. Composite samples (3 blocks × 5 plots = 15 composites) were prepared on the site less than 6 h after collection. Approximately, 100 g soil from each replicate were thoroughly mixed in a plastic bucket and transferred to plastic bags. All the samples were then immediately stored at 4°C. The 15 composite samples were processed within 3 weeks for DNA extraction, nutrient analyses and gas exchanges measurements. Soil texture was determined using the hydrometer method, and particle size distribution assigned soil samples to the loamy sand textural class (Elghamry and Elashkar, 1962). Because soil samples were collected in a 15-cm radius of the stem, the results reported in this study are representative of a microscale of the whole experimental larch plantation landscape area.

## Physicochemical analyses

Composite soil samples were dried at 20°C for 48 h, homogenized (2 mm sieve) and pulverized with mortar and pestle. Soil pH was analysed in suspensions using a 1:2.5 (w/v) soil-to-water ratio (MP220 pH-meter; Mettler Toledo, Mississauga, ON, Canada). Total soil carbon and nitrogen contents were determined using an elemental combustion system (ECS 4010; Costech Analytical Technologies Inc., Valencia, CA, USA). Analyses were performed with 18 to 123 mg pulverized soil samples and certified atropine standard ($C_{17}H_{23}NO_3$; Costech Analytical Technologies Inc.) containing 70.86% total carbon and 4.84% total nitrogen was used for calibration. Relative errors were lower than 2.8% for total carbon and 1.1% for total nitrogen analyses as observed with repeated analyses of atropine standard.

## Gaseous exchanges

Composite samples were dried at 20°C for 48 h and homogenized (2 mm sieve) and soil water content was adjusted to 20% water holding capacity. A defined amount of composite soil sample (20 g) was transferred in a 500 ml (nominal volume) Gibco® glass bottle (Wheaton Industries Inc., Millville, NJ, USA) fitted with foam plugs to allow gaseous exchanges between soil and atmosphere, while avoiding microcosm contamination with airborne particles. Soil microcosms were then transferred to an environmental chamber (MLR-350®; Sanyo, Osaka, Japan) and incubated 3 days in the dark, at 25°C and 50% relative air moisture. $H_2$, CO and $CO_2$ soil-to-air exchanges were measured using gas chromatography assays. Briefly, soil microcosm foam plugs were replaced with gastight caps equipped with butyl septa. A defined volume of two air mixtures containing either 469 ± 9 ppm $H_2$ (GTS-Welco, Minersville, PA, USA) or 508 ± 10 ppm CO (GTS-Welco) was injected to the static headspace of the microcosms, resulting in $H_2$ and CO levels of 2.5–3 ppmv. Headspace samples (10 ml) were collected with a Pressure Lok® gastight glass syringe (VICI® Precision Sampling Inc., Baton Rouge, LA, USA) and injected through the injection port of a gas chromatograph equipped with a reduction gas detector (ta3000R; Ametek Process Instruments®, Newark, DE, USA). The first-order oxidation rates were calculated by integrating the $H_2$ and CO mole fraction time series measured over a 1-h period, using at least five concentration points for data integration. Soil $CO_2$ respiration was measured in the soil microcosms fitted with gastight caps and flushed with ambient laboratory air. Linear $CO_2$ mole fraction time series was measured over a 72-h period, with four concentration points for data integration using a gas chromatograph equipped with a
methanizer and a flame ionization detector (7890B GC System; Agilent Technologies, Mississauga, ON, Canada). Sealed microcosms were incubated in the environmental chamber during the measurements. The reproducibility of gas analyses was assessed before each set of experiments by repeated analysis of three certified standard gas mixtures: (i) 2.13 ± 0.11 ppmv H₂ balance air (GTS-Welco), (ii) 2.05 ± 0.10 ppmv CO balance air (GTS-Welco) and (iii) 610 ± 12 ppmv CO₂, 5 ± 0.1 ppmv CH₄, 1 ± 0.02 ppmv N₂O balance nitrogen (Agilent Technologies), and standard deviations were < 5%. No significant gaseous exchanges were observed for blank measurements of empty microcosms. Because of the occurrence of simultaneous production and consumption of trace gases in nature and their dependence on temperature and moisture, gaseous exchanges presented in this study must be considered as potential activities. Gaseous exchanges were expressed in mol of gas h⁻¹ g⁻¹ of soil on a dry-weight basis; with soil water content measured using standard gravimetric method.

## Abundance of bacteria

The abundance of bacteria was estimated by quantification of 16S rRNA gene of bacteria in soil by a universal bacterial qPCR assay (Fierer *et al.*, 2005). Genomic DNA was extracted from 0.5 to 3.0 g(dw) composite soil samples using a combination of chemical and mechanical cell lysis procedure (Constant *et al.*, 2008). DNA was precipitated with 2 volumes 96% ethanol and polyvinylpolypyrrolidone spin column was used for final purification (Berthelet *et al.*, 1996). Purified DNA extracts from two technical replicates were pooled (200 µl in total) and kept frozen at -20°C before qPCR. The reactions were performed using 5 µl of diluted genomic DNA (1:100 and 1:500) and no quantification bias due to Taq polymerase inhibitors was observed. The qPCR assay was based on a standard curve prepared by using trilicate 10-fold dilutions of PCR-amplified standard DNA. Genomic DNA of *Burkholderia xenovorans* LB400 served as template for 16S rRNA gene standard DNA. PCR products were purified (E.Z.N.A. Cycle Pure Kit, Omega Bio-Tek®, Norcross, GA, USA) and quantified using the Quantifluor™ dsDNA System (Promega, Madison, WI, USA) according to the instructions of the manufacturers. Standard curves encompassing 10² to 10⁸ copies µl⁻¹ of standard DNA were prepared and displayed linear relationship between the signal and the logarithm copy number with reaction efficiencies of 0.99 (I² = 0.99). The Perfecta SYBR Green Fast Mix (Quanta Biosciences®, Gaithersburg, MD, USA) was used for the qPCR performed in a Rotor-Gene 6000 qPCR cycler (Corbett Life Science®, Concorde, NSW, Australia). For

unknown reasons, the extraction procedure was not successful for the sample M-A. Neither the utilization of FastDNA SPIN kit for Soil® (MP Biomedicals, Solon, OH, USA), modifications of the extraction buffer, nor increased amount of soil in the lysis procedure improved the yield of the extraction for sample M-A. This sample thus was not considered for qPCR and bacterial 16S rRNA gene ribotyping analyses.

## Bacterial 16S rRNA gene sequencing

PCR amplification of the V6-V8 regions of 16S rRNA, libraries preparation and Illumina MiSeq 250 bp paired-ends sequencing reactions were performed by the technical staff of Centre d'Innovation Génome Québec et Université McGill, resulting in 6 274 978 raw sequences. Paired-end reads were merged using the software Flash (Magoč and Salzberg, 2011) with minimum and maximum overlap length between the two reads of 20 and 250 bases respectively. The maximum proportion of mismatched base pairs tolerated in the overlap region was 30%. The 6 073 727 merged reads were processed using the software UPARSE (Edgar, 2013). Briefly, sequences were truncated to a uniform length of 420, representing the length of more than 97.5% of the sequences in the database. Reads with a low-quality score were removed using 2.0 as the maximum expected error value. The remaining 4 884 667 high-quality reads were de-replicated, sorted by size and singletons were removed. The remaining unique reads were clustered into 8248 operational taxonomic units (OTU) with the UPARSE-OTU greedy clustering method using a 97% identity threshold. The UPARSE-REF algorithm detected and removed 98 006 chimeric sequences during clustering procedure. Furthermore, 191 chimeric OTU were removed with the software UCHIME run against ChimeraSlayer 'gold' reference database (Edgar *et al.*, 2011). The final database contained 4 228 736 sequences. Libraries were normalized to the sequencing effort of the smallest 16S rRNA gene library (166 040 sequences) to avoid biases in comparative analyses introduced by sampling depth. The software QIIME version 1.8.0 (Caporaso *et al.*, 2010) was used to perform equalization of the libraries and to eliminate OTU comprising less than eight representative sequences, corresponding to a threshold of 0.005% of the total number of reads per library. QIIME was also used to pick and align one representative sequences for each OTU to assign a taxonomic classification using the Greengene reference database V13_8_99 (McDonald *et al.*, 2012). The resulting OTU table comprising the abundance and taxonomic affiliation (phylum, class, order and family levels) of the OTUs in the samples was utilized to compute alpha diversity (i.e. species richness with Ace and
Shannon indices) and the multivariate dispersion of the OTU as a measure of beta diversity (Anderson *et al.*, 2006) with the packages 'fossil' and 'vegan' implemented in the software R (Oksanen *et al.*, 2012; Vavrek, 2012). Raw sequences were deposited to the Sequence Read Archive of the National Center for Biotechnology Information under the Bioproject PRJNA280109.

## Statistical analyses

Statistical analyses were performed using the software R (R Core Development Team, 2008). The impact of MSP treatments on soil biogeochemical properties was tested using one-way analysis of variance and *post hoc* Bonferroni-corrected *t*-test. Shapiro-Wilk normality test was applied to assess normal distribution of data before the ANOVA. Distribution of soil pH, carbon content and nitrogen content followed a normal distribution, while the distribution of other variables was normalized by logarithmic ($H_2$, CO, $CO_2$ exchanges) or square root (abundance of bacteria) transformations. Normalized variables also were used for Pearson correlation analysis. Cluster analysis was computed using the function 'hclust' in the package 'stats' (R Core Development Team, 2014) to explore sampling site similarities defined by variations in gas exchange ($H_2$, CO and $CO_2$), pH, abundance of bacteria as well as total carbon and total nitrogen contents measured in soil. The environmental variables were standardized before the analysis by subtracting individual values by the average and dividing them by the standard deviation. A Euclidean distance matrix was used to generate a UPGMA agglomerative clustering of the samples according to their biogeochemical profile. The identification of statistically different soil multifunctional classes was done by performing 999 permutations of the environmental variables data set separately across the samples and comparing the observed similarity score of each cluster against the expected values under the null hypothesis using the similarity profile tool (SIMPROF) implemented in the package 'clustsig' (Clarke *et al.*, 2008). PCA was used to explore sampling sites partitioning in a reduced space defined by environmental gradients. Meaningful ordination axes whose eigenvalues were larger than the average of all eigenvalues were selected for biological interpretation. Equilibrium circle of descriptor with the radius $\sqrt{d/p}$ (where *d* is the number of dimensions of the reduced space: 2 and *p* is the total space: 10) was plotted to identify variables significantly contributing to the axes defining the position of sampling sites. Discrimination of the samples according to their ribotyping profile was performed after Hellinger transformation of the OTU table to avoid unduly relationships between explanatory variables and 16S rRNA gene composition supported by the high weight of rare species (Legendre and Gallagher,

2001). UPGMA clustering analysis was conducted to compare the samples according to their ribotyping profile. RDA was computed using standardized environmental variables with the package 'vegan' (Oksanen *et al.*, 2012), according to the comprehensive procedure described by Borcard and colleagues (Borcard *et al.*, 2011). The most parsimonious constrained model to explain bacterial ribotyping profile was obtained by forward selection of the environmental variables (Blanchet *et al.*, 2008) and permutation tests (n = 1000) were performed to assess the significance of the RDA. Indicator OTU characterizing UPGMA agglomerative clustering of the samples according to their multifunctional classification were identified using the indicator species analysis procedure implemented in the package 'indicspecies' (Dufrêne and Legendre, 1997). Minimal significance level (alpha) of the indicator OTU was 0.05, tested against 999 random permutations of samples among the biogeochemical clusters.

## Acknowledgements

This research was supported by a grant from the Natural Sciences and Engineering Research Council of Canada (NSERC) Engage Program to PC (grant EGP 463612-14). The work of VP was supported by a NSERC - Undergraduate Student Research Awards. The authors are grateful to the personnel staff of *Centre d'Innovation Génome Québec et Université McGill* for preparation of 16S rRNA gene libraries and sequencing services.

## Conflict of Interest

All the authors have no conflict of interest to declare.

## References

* Anderson, M.J., Ellingsen, K.E., and McArdle, B.H. (2006) Multivariate dispersion as a measure of beta diversity. *Ecol Lett*, **9**: 683-693.

* Azeria, E.T., Fortin, D., Hébert, C., Peres-Neto, P., Pothier, D., and Ruel, J.-C. (2009) Using null model analysis of species co-occurrences to deconstruct biodiversity patterns and select indicator species. *Diversity Distrib*, **15**: 958-971.

* Balvanera, P., Siddique, I., Dee, L., Paquette, A., Isbell, F., Gonzalez, A., *et al.* (2014) Linking biodiversity and ecosystem services: current uncertainties and the necessary next steps. *Bioscience*, **64**: 49-57.

* Bartell, S.M. (2006) Biomarkers, bioindicators, and ecological risk assessment- A brief review and evaluation. *Environ Bioindicators*, **1**: 60-73.

* Berthelet, M., Whyte, L.G., and Greer, C.W. (1996) Rapid, direct extraction of DNA from soils for PCR analysis using polyvinylpolypyrrolidone spin columns. *FEMS Microbiol Lett*, **138**: 17-22.
Blanchet, F.G., Legendre, P., and Borcard, D. (2008) Forward selection of explanatory variables. *Ecology*, **89**: 2623-2632.

Borcard, D., Gillet, F. and Legendre, P. (2011) *Numerical Ecology with R*. New York: Springer.

Buitrago, M., Paquette, A., Thiffault, N., Bélanger, N., and Messier, C. (2015) Early performance of planted hybrid larch: effects of mechanical site preparation and planting depth. *New Forest*, **46**: 319-337.

Caporaso, J.G., Kuczynski, J., Stombaugh, J., Bittinger, K., Bushman, F.D., Costello, E.K., et al. (2010) QIIME allows analysis of high-throughput community sequencing data. *Nat Meth*, **7**: 335-336.

Clarke, K.R., Somerfield, P.J., and Gorley, R.N. (2008) Testing of null hypotheses in exploratory community analyses: similarity profiles and biota-environment linkage. *J Exp Mar Biol Ecol*, **366**: 56-69.

Constant, P., Poissant, L., and Villemur, R. (2008) Isolation of *Streptomyces* sp. PCB7, the first microorganism demonstrating high-affinity uptake of tropospheric H₂. *ISME J*, **2**: 1066-1076.

Constant, P., Poissant, L., and Villemur, R. (2009) Tropospheric H₂ budget and the response of its soil uptake under the changing environment. *Sci Total Environ*, **407**: 1809-1823.

Constant, P., Chowdhury, S.P., Hesse, L., Pratscher, J., and Conrad, R. (2011) Genome data mining and soil survey for the novel group 5 [NiFe]-hydrogenase to explore the diversity and ecological importance of presumptive high affinity H₂-oxidizing bacteria. *Appl Environ Microbiol*, **77**: 6027-6035.

Drake, J.E., Darby, B.A., Giasson, M.-A., Kramer, M.A., Phillips, R.P., and Finzi, A.C. (2013) Stoichiometry constrains microbial response to root exudation- insights from a model and a field experiment in a temperate forest. *Biogeosciences*, **10**: 821-838.

Dufrene, M., and Legendre, P. (1997) Species assemblages and indicator species: the need for a flexible asymmetrical approach. *Ecol Monogr*, **67**: 345-366.

Edgar, R.C. (2013) UPARSE: highly accurate OTU sequences from microbial amplicon reads. *Nat Meth*, **10**: 996-998.

Edgar, R.C., Haas, B.J., Clemente, J.C., Quince, C., and Knight, R. (2011) UCHIME improves sensitivity and speed of chimera detection. *Bioinformatics*, **27**: 2194-2200.

Ehhalt, D.H., and Rohrer, F. (2009) The tropospheric cycle of H₂: a critical review. *Tellus B*, **61**: 500-535.

Elghamry, W., and Elashkar, M. (1962) Simplified textural classification triangles. *Soil Sci Soc Am J*, **26**: 612-613.

van Elsas, J.D., Chiurazzi, M., Mallon, C.A., Elhottova, D., Kristúfek, V., and Salles, J.F. (2012) Microbial diversity determines the invasion of soil by a bacterial pathogen. *Proc Natl Acad Sci USA*, **109**: 1159-1164.

Fierer, N., Jackson, J.A., Vilgalys, R., and Jackson, R.B. (2005) Assessment of soil microbial community structure by use of taxon-specific quantitative PCR assays. *Appl Environ Microbiol*, **71**: 4117-4120.

Gamfeldt, L., Hillebrand, H., and Jonsson, P.R. (2008) Multiple functions increase the importance of biodiversity for overall ecosystem functioning. *Ecology*, **89**: 1223-1231.

Gödde, M., Meuser, K., and Conrad, R. (2000) Hydrogen consumption and carbon monoxide production in soils with different properties. *Biol Fertil Soils*, **32**: 129-134.

Haeussler, S., Bartemucci, P., and Bedford, L. (2004) Succession and resilience in boreal mixedwood plant communities 15-16 years after silvicultural site preparation. *For Ecol Manage*, **199**: 349-370.

Hector, A., and Bagchi, R. (2007) Biodiversity and ecosystem multifunctionality. *Nature*, **448**: 188-190.

Jandl, R., Lindner, M., Vesterdal, L., Bauwens, B., Baritz, R., Hagedorn, F., et al. (2007) How strongly can forest management influence soil carbon sequestration? *Geo derma*, **137**: 253-268.

Kabrick, J.M., Dey, D.C., Sambeek, J.W.V., Wallendorf, M., and Gold, M.A. (2005) Soil properties and growth of swamp white oak and pin oak on bedded soils in the lower Missouri River floodplain. *For Ecol Manage*, **204**: 315-327.

Khdhiri, M., Hesse, L., Popa, M.E., Quiza, L., Lalonde, I., Meredith, L.K., et al. (2015) Soil carbon content and relative abundance of high affinity H₂-oxidizing bacteria predict atmospheric H₂ soil uptake activity better than soil microbial community composition. *Soil Biol Biochem*, **85**: 1-9.

King, C.E., and King, G.M. (2014) Description of *Thermogemmatispora carboxidivorans* sp. nov., a novel carbon monoxide-oxidizing member of the *Ktedonobacteria* isolated from a geothermally-heated biofilm, and analysis of carbon monoxide oxidation by members of the *Ktedonobacteria*. *Int J Syst Evol Microbiol*, **64**: 1244-1251.

King, G.M., and Weber, C.F. (2007) Distribution, diversity and ecology of aerobic CO-oxidizing bacteria. *Nat Rev Micro*, **5**: 107-118.

Lauber, C.L., Hamady, M., Knight, R., and Fierer, N. (2009) Pyrosequencing-based assessment of soil pH as a predictor of soil bacterial community structure at the continental scale. *Appl Environ Microbiol*, **75**: 5111-5120.

Legendre, P., and Gallagher, E. (2001) Ecologically meaningful transformations for ordination of species data. *Ecologia*, **129**: 271-280.

Löf, M., Dey, D., Navarro, R., and Jacobs, D. (2012) Mechanical site preparation for forest restoration. *New Forest*, **43**: 825-848.

Magoč, T., and Salzberg, S.L. (2011) FLASH: fast length adjustment of short reads to improve genome assemblies. *Bioinformatics*, **27**: 2957-2963.

Martin, B., Marc, L., Nelson, T., Alain, P., Luc, L., Louis, B., et al. (2014) Issues and solutions for intensive plantation silviculture in a context of ecosystem management. *For Chronicle*, **90**: 748-762.

McDonald, D., Price, M.N., Goodrich, J., Nawrocki, E.P., DeSantis, T.Z., Probst, A., et al. (2012) An improved Greengenes taxonomy with explicit ranks for ecological and evolutionary analyses of bacteria and archaea. *ISME J*, **6**: 610-618.

Meredith, L.K., Rao, D., Bosak, T., Klepac-Ceraj, V., Tada, K.R., Hansel, C.M., et al. (2013) Consumption of atmospheric hydrogen during the life cycle of soil-dwelling actinobacteria. *Environ Microbiol Rep*, **6**: 226-238.

Oksanen, J., Blanchet, F., Kindt, R., Legendre, P., Minchin, P., O'Hara, R., et al. (2012) Vegan: community ecology
package. R package version 2.0-4. URL http://cran.r-project.org/package=vegan.

Paquette, A., and Messier, C. (2009) The role of plantations in managing the world's forests in the Anthropocene. *Front Ecol Environ, **8**: 27–34.

Pearce, J., and Venier, L. (2005) Small mammals as bioindicators of sustainable boreal forest management. *For Ecol Manage, **208**: 153–175.

Pearce, J.L., and Venier, L.A. (2006) The use of ground beetles (Coleoptera: Carabidae) and spiders (Araneae) as bioindicators of sustainable forest management: a review. *Ecol Indic, **6**: 780–793.

Prévost, M. (1992) Effets du scarifiage sur les propriétés du sol, la croissance des semis et la compétition: revue des connaissances actuelles et perspectives de recherches au Québec. *Ann For Sci, **49**: 277–296.

Quiza, L., Lalonde, I., Guertin, C., and Constant, P. (2014) Land-use influences the distribution and activity of high affinity CO-oxidizing bacteria associated to type I-coxL genotype in soil. *Front Microbiol, **5**: 1–15. doi:10.3389/fmicb.2014.00271.

R Core Development Team (2008) R: A language and environment for statistical computing. Vienna, Austria: Computing RfS.

R Core Development Team (2014) The R Stats Package - version 3.1.3. URL https://stat.ethz.ch/R-manual/R-devel/library/stats/html/hclust.html.

Simon, B.-G., David, P., Christian, M., and Nicolas, B. (2013) Root production of hybrid poplars and nitrogen mineralization improve following mounding of boreal Podzols. *Can J For Res, **43**: 1092–1103.

Staddon, W., Duchesne, L., and Trevors, J. (1999) The role of microbial indicators of soil quality in ecological forest management. *Forest Chron, **75**: 81–86.

Sutton, R.F. (1993) Mounding site preparation: a review of European and North American experience. *New Forest, **7**: 151–192.

Tardy, V., Mathieu, O., Lévêque, J., Terrat, S., Chabbi, A., Lemanceau, P., et al. (2014) Stability of soil microbial structure and activity depends on microbial diversity. *Environ Microbiol Rep, **6**: 173–183.

Thiffault, N., and Jobidon, R. (2006) How to shift unproductive *Kalmia angustifolia* – *Rhododendron groenlandicum*

heath to productive conifer plantation. *Can J For Res, **36**: 2364–2376.

Thomson, J.R., Fleishman, E., Nally, R.M., and Dobkin, D.S. (2007) Comparison of predictor sets for species richness and the number of rare species of butterflies and birds. *J Biogeogr, **34**: 90–101.

Vavrek, M.J. (2012) Fossil: Palaeoecological and palaeogeographical analysis tools - R Package version 0.3.7. URL http://cran.r-project.org/web/packages/fossil/index.html.

Wagg, C., Bender, S.F., Widmer, F., and van der Heijden, M.G.A. (2014) Soil biodiversity and soil community composition determine ecosystem multifunctionality. *Proc Ntl Acad Sci, **111**: 5266–5270.

Weber, C.F., and King, G.M. (2010) Distribution and diversity of carbon monoxide-oxidizing bacteria and bulk bacterial communities across a succession gradient on a Hawaiian volcanic deposit. *Environ Microbiol, **12**: 1855–1867.

## Supporting information

Additional Supporting Information may be found in the online version of this article at the publisher's web-site:

**Fig. S1.** Distribution of the ribotypes classified at the phylum taxonomic level.

**Fig. S2.** Beta diversity of the ribotypes (OTU defined at 97% identity threshold) as assessed using multivariate dispersion measure (Anderson *et al.*, 2006).

**Fig. S3.** Relative abundance of ubiquitous ribotypes detected in all soil samples (OTU were classified at the phylum level).

**Fig. S4.** PCR detection of (a) OTU 398 and (b) OTU 3283 identified as potential indicator for soil samples characterized by baseline multifunctional attributes and soil samples that are divergent from baseline attributes respectively.

**Table S1.** Oligonucleotides and PCR conditions utilized to detect bioindicator for baseline soil multifunctional attributes and soil samples that are divergent from the baseline multifunctional attributes.