H₂-saturation of high affinity H₂-oxidizing bacteria alters the ecological niche of soil microorganisms unevenly among taxonomic groups

Soil sample

Soil sample was collected in an agricultural land located in St. Claude (Québec, Canada) on the south shore of the St. Lawrence River (45.6809°N, −71.9969°W). The field is managed with fallow, potatoes and maize crop rotation. Potato seedlings (approximately 10 cm height) were present on the site during sampling in July 2013. The top layer of the A-horizon (0–10 cm depth) was collected, stored in plastic bags at 4 °C and processed within a week. Soil was air-dried for 48 h in the laboratory and sieved (2 mm mesh size) through a vibratory sieve shaker AS 200^® (Retsch GmbH, Haan, Germany) before preparation of soil microcosms. Soil was classified as sandy clay loam according to soil textural class parameters identified with the hydrometer method (Bouyoucos, 1936). Soil pH was determined with 1:2 soil-water suspensions with an Accumet pH-meter (Fisher Scientific, Hampton, NH). Total carbon (3.1 ± 0.3%) and total nitrogen (0.3 ± 0.0%) content were determined using an elemental combustion system using the protocol described in Khdhiri et al. (2015).

Controlled H₂ exposure in dynamic soil microcosm chambers

Dynamic microcosm chambers were designed to expose soil to controlled levels of H₂ (Fig. S1). Microcosm chambers consisted of 0.9 L cell culture flasks (Corning, Tewksbury, MA, USA) equipped with a rubber stopper fitted with two 1/8” outside diameter PTFE (Teflon) tubes: the first supplied gas mixture to the microcosm chamber and the second was vented to the atmosphere. Synthetic gas mixtures supplied the microcosm chambers at a flowrate of 40 ml min⁻¹, resulting in a dynamic headspace with a residence time of approximately 22 min in the enclosures. Gas mixtures were bubbled in water before entering in microcosms to prevent soil dryness. Two different H₂ treatments were applied in parallel incubations. The first treatment (designated microcosms eH₂(a) and eH₂(b)), named elevated H₂ (eH₂), consisted of exposing soil to a dynamic headspace comprising 525 ppmv H₂ in synthetic air, simulating H₂ concentrations detected around N₂-fixing nodules in soil (Hunt & Layzell, 1993). The second was a control treatment (designated microcosms aH₂(a) and aH₂(b)), named atmospheric H₂ (aH₂), where soil was exposed to a dynamic headspace comprising 0.54 ppmv H₂ in synthetic air, representing H₂ concentrations found in the atmosphere. Both treatments were replicated, resulting in four soil microcosms in total. Each microcosm chamber contained 200 g_(dw) soil at the beginning of the incubation period. This amount of soil ensured sufficient material for the monitoring of microbiological and physicochemical variables throughout the incubation period, while avoiding diffusion limitation of the H₂ soil uptake rate measurements. Diffusion limitation was precluded since preliminary experiments using the same microcosm setup showed proportional H₂ uptake activity as a function of the amount of soil in the chamber using 150, 200 and 250 g_(dw) soil samples (raw data file). Soil water content was adjusted to 20% water holding capacity before incubation for 10 days at 28 °C in the dark. Synthetic gas mixture supply was continuously maintained, with the exception of routine soil subsamples collection and high affinity H₂ oxidation rate measurements (see below). A decrease in soil pH was observed over the course of the incubation (from pH 5.9 ± 0.2 to 5.1 ± 0.1) in all microcosms, without distinction between aH₂ and eH₂ treatments. Soil moisture was monitored using standard gravimetric method and maintained at 26 ± 5% throughout the incubation period. Blank (empty) microcosms incubated prior to the experiment did not show any H₂ oxidation or production activity.

H₂ uptake activity

High affinity H₂ oxidation activity was routinely monitored throughout the incubation of the soil microcosms. Briefly, microcosms were disconnected from their respective gas supply, flushed 5 min with a synthetic gas mixture (0.54 ppmv H₂) and tightly closed with rubber septum caps. The 5-min flush was shown sufficient to avoid residual H₂ degassing that would otherwise have led to an underestimation of H₂ oxidation rate. A defined volume of H₂ gas mixture (525 ± 10 ppmv H₂ GST-Welco, PA, USA.) was injected to obtain an initial concentration of approximately 3 ppmv in the static headspace. Decrease of the H₂ mixing ratio was monitored as a function of time by analyzing aliquots (10 ml) of the headspace air in a ta3000R gas chromatograph equipped with a reduction gas detector (Ametek Process Instruments^®, DE, USA.) as previously described (Khdhiri et al., 2015). Considering the low level of H₂ added in the headspace, H₂ uptake reflected the activity of high affinity HOB, knallgas bacteria showing low affinity for H₂ being unable to use these trace amounts of H₂ (Conrad, Aragno & Seiler, 1983). The biphasic kinetic parameters governing H₂ oxidation activity in soil (i.e., _(app)K_m and _(app)V_max) were measured at the end of the 10-day incubation period after the addition of specified amounts of pure H₂ into the headspace of the microcosms, as previously described (Schuler & Conrad, 1990).

DNA extraction, qPCR and high-throughput sequencing of PCR-amplified bacterial 16S rRNA gene

Soil subsamples (approximately 10 g per subsample, without replacement) were collected in the four microcosms after 0, 1, 3, 5, 7 and 10 incubation days to investigate the taxonomic structure of microbial communities. The extraction of total genomic DNA was performed using the FastDNA SPIN kit for Soil^® (MP Biomedicals, Solon, OH, USA.). DNA was eluted in 100 µL nuclease-free water. Quality of the DNA was then examined on agarose gels and samples were quantified using Quantifluor dsDNA System^® (Promega, Fitchburg, WI, USA). The V4 region of bacterial 16S rRNA gene was PCR amplified (Table S1) and sequenced on an Illumina MiSeq 2000 instrument as multiplexed libraries to generate paired-end reads (2 × 250 bp). Quantification of bacterial 16S rRNA and hhyL genes was done by qPCR following Khdhiri et al. (2015).

Quality control of sequencing reads

Our sequencing analysis pipeline using QIIME 1.8.0 (Kuczynski et al., 2011) software was based on the Itagger pipeline (Tremblay et al., 2015). V4 region Illumina sequencing reads were first processed using the split_library function implemented in QIIME for samples demultiplexing and forward primer removal. The assembly of paired-end reads was performed using FLASH software version 1.2.10 (Magoc & Salzberg, 2011). The software DUK (Li, Copeland & Han, 2011), a Kmer-based sequence matching tool, was then used to remove common sequencing contaminants as well as PhiX sequences added as control. Afterwards, sequences were trimmed to remove low quality ends and a general quality control was also applied. All reads containing at least a single N (ambiguous base), an average quality score of less than 30 or more than 15 bases with a quality score under 20 were discarded. The minimum and maximum read length were respectively set to 200 and 300 bases for the assembled paired-reads. Reverse primers were also removed using the truncate_reverse_primer function implemented in QIIME. Raw sequences were deposited in the Sequence Read Archive of the National Center for Biotechnology Information under the Bioproject PRJNA295403.

Clustering and taxonomic identification of sequencing reads

Processed 2,554,543 high quality reads were clustered into OTUs (Operational Taxonomic Units) using the function pick_otus implemented in QIIME along with the USEARCH (Edgar, 2010) clustering algorithm. Sequences dereplication was performed at 100% sequence identity and clusters were formed by denoising unique sequences at 99% identity. Singletons were removed to avoid diversity bias. Remaining clusters were filtered using UCHIME (Edgar et al., 2011) chimera filter in de novo mode, followed by a reference-based filter step against the Gold reference database. These 2 chimera scanning steps allow a better removal of chimeras by removing sequences flagged as chimeras by any of the two methods. Chimera-checked clusters were then clustered into OTUs at 97% identity threshold also using the function pick_otus with the USEARCH software. Taxonomic identification of OTUs was performed using the assign_taxonomy function along with the naïve Bayesian Ribosomal Database Project (RDP) classifier 2.7 (Wang et al., 2007). Taxonomy was assigned based on the Greengenes taxonomy with the 16S rRNA Greengenes reference database (DeSantis et al., 2006). A rarefaction step was applied to the OTU libraries to standardize all libraries, by random subsampling, to the lowest amount of sequences (74,316 reads) to avoid bias introduced due to unequal sequencing efforts of the samples. This rarefied OTU table comprising the frequency distribution of the OTUs in each sample was used in downstream statistical and co-occurrence network analyses.

Statistical analyses

Statistical analyses were performed using the software R version 3.0.2 (R Development Core Team, 2008). Pairwise comparison of H₂ uptake activities in soil microcosms incubated under aH₂ and eH₂ exposure treatments was tested using one-way analysis of variance (ANOVA) and post hoc Tukey test. Shapiro–Wilk normality test was applied to confirm normal distribution of data before the ANOVA. The package “nlme” was used to compute nonlinear regression for modeling the first-order decay of hhyL gene abundance determined by qPCR. Discrimination of the samples according to their ribotyping profile was performed with the package “vegan.” Hellinger transformation of the OTU table (97% identity threshold) was computed due to the presence of zeros in the OTU table (Legendre & Gallagher, 2001). Euclidean distance matrix was used to generate an UPGMA agglomerative clustering of the samples. Identification of statistically different clusters was done by performing 999 permutations of the OTU table dataset 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” 1.1 (Clarke, Somerfield & Gorley, 2008). Pairwise comparison of relative abundance of OTUs having a higher contribution than average to explain the two dimensions of the PCA space in soil microcosms incubated under aH₂ and eH₂ exposure treatments was tested using Kruskal–Wallis and post hoc Tukey test. Pairwise comparison of relative abundance of all OTUs in soil microcosms incubated under aH₂ and eH₂ exposure treatments was tested using the likelihood ratio test using the package “edgeR.” Covariation among OTUs during the incubation under controlled H₂ levels was analyzed by correlation networks using the package “WGCNA” 1.41 (Langfelder & Horvath, 2008). A detailed methodology for the computation of correlation networks is provided in Text S1.

Impact of H₂ exposure on the distribution and activity of HOB

Agricultural soil microcosms were incubated under two different H₂ level treatments. The first treatment simulated soil exposure to atmospheric level of H₂ (aH₂; 0.54 ppmv), while the second simulated high affinity hydrogenases substrate saturation (eH₂; 525 ppmv). High affinity H₂ oxidation rates were not significantly different (ANOVA, P > 0.05) between the two pairs of microcosms at the beginning of the incubation (Fig. 1). The influence of H₂ exposure could already be observed after 24 h. From hour 24 to day 10, a decline of high affinity H₂ uptake rate was observed in microcosms exposed to eH₂, while an increase of this uptake rate was measured in microcosms exposed to aH₂ levels. This trend was maintained over the course of the incubation, with an oxidation rate of 6.2 ± 1.9 and $0.3\pm 0.1\mathrm{nmol}{\mathrm{g}}_{\left(\mathrm{dw}\right)}^{-1}{\mathrm{h}}^{-1}$ in aH₂ and eH₂ treatments at the end of the experiment, respectively. This response of H₂ oxidation rate was accompanied by an alteration in the abundance of presumptive high affinity HOB in soil. The abundance of hhyL gene determined by qPCR did not change over the course of the incubation under aH₂ (exponential regression, P = 0.06), while a slight but significant exponential decay (exponential regression, P = 0.02) was observed under eH₂ (Fig. 2A). Indeed, the physiological role of high affinity hydrogenases differs within various taxonomic groups of bacteria. In streptomycetes, the enzyme is primarily expressed in spores to support a seed bank under a mixotrophic survival energy mode (Constant et al., 2010; Constant, Poissant & Villemur, 2008; Liot & Constant, 2015), while Mycobacterium express the enzyme in the exponential and stationary phases for mixotrophic growth as well as survival (Berney & Cook, 2010; Greening et al., 2015b; Greening et al., 2014b). Under axenic cultivation conditions, eH₂ level would be expected to favour an increase of Mycobacterium biomass and persistence of streptomycetes spores. It should be noted that our qPCR data cannot differentiate between active and inactive cells but the decreasing trend of hhyL copy number in soil exposed to eH₂ suggests that H₂ was not sufficient to promote growth and persistence of high affinity HOB.

H₂ soil exposure exerted a significant impact on the kinetic parameters governing H₂ oxidation activity in soil. At the end of the incubation, microcosms exposed to aH₂ level displayed _(app)K_m of 40 ± 5 ppmv (Fig. 3A), while microcosms exposed to eH₂ showed lower affinity towards H₂ with an _(app)K_m of 838 ± 163 ppmv (Fig. 3B). Similar observations were obtained by Dong & Layzell (2001), where soil exposed to elevated levels of H₂ (600 ppmv) displayed an _(app)K_m of 1,028 ppmv H₂, while soil exposed to low H₂ level (0.55 ppmv) was characterized by an _(app)K_m of 40 ppmv. The biphasic kinetic parameters characteristics of low and high affinity HOB were also computed using Eadie-Hofstee plots to highlight high- and low-affinity H₂ uptake activities in both H₂ treatments. The H₂ oxidation activity was mainly catalyzed by bacteria demonstrating intermediate and high-affinity in aH₂ treatment (Fig. 3A, insert), while H₂-oxidation activity by bacteria displaying low- and high-affinity for H₂ were detectable in the eH₂ treatment (Fig. 3B, insert). Low H₂ exposure supported the metabolism of intermediate- to high-affinity HOB, such as some Actinobacteria species of Streptomyces, Rhodococcus and Mycobacterium (Berney et al., 2014; Constant et al., 2011b; Golding et al., 2012; Liot & Constant, 2015; Meredith et al., 2014; Schäfer, Friedrich & Lenza, 2013). Triggering of low-affinity H₂ oxidation activity under eH₂ treatment might be explained by knallgas bacteria such as Proteobacteria species encompassing Ralstonia, Variovorax and Bradyrhizobium using H₂ for mixotrophic growth (Rittenberg & Goodman, 1969), as they are known to exhibit a H₂ uptake threshold ranging from 1 to 200 ppmv H₂ (Conrad, 1996).

Impact of H₂ exposure on soil bacterial community structure

Soil subsamples were collected in the microcosms during the incubation period to evaluate the temporal variation of soil bacterial community structure. H₂ treatment did not influence the abundance of 16S rRNA gene in soil (Fig. 2B). This is in contrast with another study that reported more than twofold increase of the amount of bacterial cells upon elevated H₂ exposure, but a mistake in the reported H₂ concentration units by the authors (i.e., the use of a gas mixture of 2,000 ml H₂ per L of synthetic air) impairs a sound comparison with the present investigation (Stein et al., 2005). High-throughput sequencing of PCR-amplified 16S rRNA gene unveiled that bacterial communities were dominated by Proteobacteria (34%), Acidobacteria (20%), Actinobacteria (10%), Verrucomicrobia (8%) and Bacteroidetes (5%) (Fig. S2). Species richness was not affected by H₂ exposure (ANOVA, P > 0.05), with Shannon indices of 9.38 ± 0.05 and 9.34 ± 0.17 for aH₂ and eH₂ treatment respectively. An UPGMA hierarchical clustering analysis was computed to compare bacterial community profiles from soil subsamples collected in microcosms exposed to different H₂ concentrations. The four microcosms encompassed the same cluster before the incubation, indicating high similarity of their initial ribotyping profile (Fig. 4A). No coherent impact of H₂ exposure on microbial community structure was observed after 1 and 3 incubation days. Ribotyping profiles obtained for sub-samples collected after 5 and 7 incubation days were separated with confidence in different clusters, suggesting a response of microbial communities to H₂ exposure. This discrimination was however transient as ribotyping profiles corresponding to the last incubation day could not be discriminated according to their respective H₂ treatments (Fig. 4A). The transient response was likely due to the fact that H₂ alone was not sufficient to maintain a sustained change in microbial community structure. Decrease of complex carbon sources and nutrients used along with H₂ for mixotrophic growth in soil is a potential explanation for the convergence of the ribotyping profiles after 7 incubations days.

A PCA was computed to identify OTUs contributing to the clusterization of ribotyping profiles. The ordination space defined by the first two components explained 38.6% of the variation observed (Fig. 4B). Three OTUs defined an important proportion of the reduced space represented by both axes (Fig. 5). The relative abundance of OTU 52 (classified as a member of the class Betaproteobacteria) contributed to distinguish the four soil subsamples collected before the incubation from the 20 other subsamples collected after 1, 3, 5, 7 and 10 incubation days along the first axis as it decreased after the beginning of the incubation period (Kruskal–Wallis, P < 0.05). Two OTUs classified as members of the order Bacillales (OTUs 65 and 299) contributed to discriminate both H₂ treatments in the ordination space after 5 and 7 incubation days (Kruskal–Wallis, P < 0.05). These OTUs were more abundant in microcosms exposed to eH₂ until day 10, where they reached a similar relative abundance in aH₂ and eH₂ treatments. The transient response of these OTUs in eH₂ treatment could be a direct consequence of H₂ exposure as representatives of the genus Bacillus possess putative membrane-bound type 1a and 1d [NiFe]-hydrogenases (Greening et al., 2016).

Conflicting results were obtained in two different studies reporting the impact of H₂ soil exposure on bacterial community structure based on 16S rRNA terminal restriction fragment profiles analysis. H₂ exerted no incidence on ribotyping profile in soil exposed to 250 nmol H₂cm⁻³h⁻¹ (500 ppmv H₂ in artificial air, added at 45 ml min⁻¹), while a single ribotype related to Mycobacterium increased in soils upon elevated H₂ exposure (Osborne, Peoples & Janssen, 2010). On the other hand, an exposure rate of 33 nmol cm⁻³ h⁻¹ (79 ppmv H₂ at 100 cm³ min⁻¹) exerted a significant influence on bacterial community profile, with an increase of T-RFLP peaks belonging to γ-Proteobacteria and a decrease of peaks belonging to Actinobacteria and α-Proteobacteria upon elevated H₂ exposure (Zhang, He & Dong, 2009). A likelihood ratio analysis, fitting a negative binomial generalized log-linear model to the sequencing data (McCarthy, Chen & Smyth, 2012; Robinson, McCarthy & Smyth, 2010), unveiled that distribution of 958 OTUs was influenced by H₂ exposure. The relative abundance of OTUs responding to H₂ treatment ranged between 0.001 and 1.8%, suggesting an incidence of H₂ on members of the rare biosphere and abundant taxa. The influence of H₂ was uneven among different taxonomic groups as different representatives of the same taxa (i.e., OTUs classified at the phylum and the order taxonomic levels) were found to be favored or disfavored in response to eH₂ exposure (Table 1). Even though the uneven response among taxonomic groups impairs prediction of H₂ exposure on metabolic functions in soil (Langille et al., 2013), this observation is sufficient to demonstrate that previous investigations relying on low-resolution community fingerprinting techniques have considerably underestimated the response of soil microbial communities to H₂ exposure.

Impact of H₂ exposure on the co-occurrence of OTUs

Correlation networks were computed to investigate the impact of H₂ exposure on the covariation of OTUs throughout the incubation period and identify OTUs for which the distribution is influenced by H₂. Two separate networks were computed: the first for elevated H₂ exposure (eH₂ network) and a second for low H₂ exposure (aH₂ network). Both networks contained 43 modules and module preservation statistics demonstrated that module composition significantly changed across H₂ treatments (Text S1). The significances of measured H₂ oxidation rate in explaining network structure were 0.33 and 0.98 in aH₂ and eH₂ networks, again pointing out a significant impact of H₂ on soil microbial community structure. Indeed, module eigengenes of five modules were significantly correlated with high-affinity H₂ oxidation rate time series in eH₂ network, while no module was related to the activity in aH₂ network (Fig. S3). Together, these 5 modules represent 1,140 OTUs, which is more than a third of the whole eH₂ network (3,154 OTUs). OTUs belonging to one module were hindered by eH₂ exposure, while members of the four other modules were favored by eH₂ treatment. The OTUs belonging to these five modules were members of the rare biosphere as well as more abundant ribotypes, with relative abundance ranging between 0.001% and 4.3% encompassing Proteobacteria, Chloroflexi, Acidobacteria as well as other phyla (Fig. S4). Clustering of these OTUs at the class and order levels unveiled that none of the taxonomic groups were restricted to the five modules correlated with H₂ oxidation rate, supporting the previous observation that response to H₂ exposure and the distribution of hydrogenase genes is uneven within each taxonomic group. The 38 modules in eH₂ network whose eigengene showed no significant correlation with H₂ oxidation rate are ecologically relevant observations since they represent indirect impacts of H₂ exposure on soil microbial communities and highlight a previously overlooked role of H₂ in shaping potential microbe-microbe cooperation and competition interactions. Deciphering the impact of these complex interactions on soil biogeochemical processes was beyond the scope of this study but will deserve attention in future investigations.