Consortia of low-abundance bacteria drive sulfate reduction-dependent degradation of fermentation products in peat soil microcosms

Dissimilatory sulfate reduction in peatlands is sustained by a cryptic sulfur cycle and effectively competes with methanogenic degradation pathways. In a series of peat soil microcosms incubated over 50 days, we identified bacterial consortia that responded to small, periodic additions of individual fermentation products (formate, acetate, propionate, lactate or butyrate) in the presence or absence of sulfate. Under sulfate supplementation, net sulfate turnover (ST) steadily increased to 16–174 nmol cm–3 per day and almost completely blocked methanogenesis. 16S rRNA gene and cDNA amplicon sequencing identified microorganisms whose increases in ribosome numbers strongly correlated to ST. Natively abundant (⩾0.1% estimated genome abundance) species-level operational taxonomic units (OTUs) showed no significant response to sulfate. In contrast, low-abundance OTUs responded significantly to sulfate in incubations with propionate, lactate and butyrate. These OTUs included members of recognized sulfate-reducing taxa (Desulfosporosinus, Desulfopila, Desulfomonile, Desulfovibrio) and also members of taxa that are either yet unknown sulfate reducers or metabolic interaction partners thereof. Most responsive OTUs markedly increased their ribosome content but only weakly increased in abundance. Responsive Desulfosporosinus OTUs even maintained a constantly low population size throughout 50 days, which suggests a novel strategy of rare biosphere members to display activity. Interestingly, two OTUs of the non-sulfate-reducing genus Telmatospirillum (Alphaproteobacteria) showed strongly contrasting preferences towards sulfate in butyrate-amended microcosms, corroborating that closely related microorganisms are not necessarily ecologically coherent. We show that diverse consortia of low-abundance microorganisms can perform peat soil sulfate reduction, a process that exerts control on methane production in these climate-relevant ecosystems.


Nucleic acid extraction
Soil samples were ground in liquid N₂ and separated into approx. 300 mg aliquots, followed by total nucleic acids extraction as described previously (Leininger et al., 2006). Nucleic acids were purified with the OneStep PCR inhibitor removal kit (Zymo Research, Irvine, CA, USA) and split into two fractions. Each fraction was either treated with RNase (RNase ONE ribonuclease, Promega, Fitchburg, WI, USA) followed by sodium acetate precipitation (Sambrook and Russell, 2001) or DNase (TURBO DNA-free kit, Thermo Fisher Scientific, Waltham, MA, USA) to obtain RNA-free DNA or DNA-free RNA, respectively. PicoGreen dsDNA and RiboGreen RNA assays (Thermo Fisher Scientific) were used for nucleic acids quantification. cDNA was obtained using the SuperScript III first-strand synthesis kit (Thermo Fisher Scientific) with either random hexamer primers for subsequent amplicon sequencing or the primer DSP821R (Pester et al., 2010) for subsequent qPCR analysis of Desulfosporosinus 16S rRNA (Pester et al., 2010).
16S rRNA gene and cDNA amplicon sequence data processing 34.4 million sequences were obtained from the native soil and all microcosm incubations at days 5, 26, and 50. Quality control, trimming, paired-end assembly, and clustering of obtained MiSeq reads was performed with the iTagger pipeline version 1.1 using the default 16S rRNA V4 region configuration (http:/ /bitbucket.org/berkeleylab/jgi_itagger/). In short, DUK (http:/ /duk.sourceforge.net/) was used for removal of PhiX control sequences, sequencing library adapter dimers, and other known contaminants (for example, human sequences), followed by trimming of the PCR primer with cutadapt (Martin, 2011). Paired-end assembly was performed using FLASH (Magoč and Salzberg, 2011) and nonmerged/single-end reads were discarded. Final quality control consisted of filtering with a cutoff of 2 errors (below defined quality threshold) per 100 nt. All sequences represented by at least two identical reads were used for clustering into operational taxonomic units (OTUs) using UPARSE (Edgar, 2013). To avoid insufficient clustering, an iterative approach was applied by using a 99%, then 98%, and then 97% sequence identity threshold. Singletons that were removed in previous steps were assigned to existing OTUs by iTagger, if possible, or discarded. Reference-based chimera filtering was performed with UCHIME (Edgar et al., 2011). 5.4%, 27.9%, 4.5%, 0.1% of all reads were removed during quality control, UPARSE clustering (including de novo chimera filtering), singletons filtering, and referencebased chimera filtering using UCHIME, respectively. Alpha diversity metrics (observed OTUs, Chao1, ACE, and Good's coverage) were determined (Supplementary Figure S3a) using R (R Core Team, 2015) and phyloseq (McMurdie and Holmes, 2013).

Sulfate reduction-associated bacteria in peat soil
Supplementary Information

Phylogenetic tree reconstruction
The phylogeny of selected OTUs was analyzed by calculating trees in ARB (Ludwig et al., 2004). Reference sequence RAxML trees were inferred based on 817 (Telmatospirillum) or 843 (Desulfosporosinus) alignment positions without applying a conservation filter (Stamatakis, 2014). Tree topology was evaluated by bootstrap analysis (100 resamplings) with RAxML (Stamatakis, 2014). Shorter V4 amplicon sequences were added to the reference trees by using ARB's parsimony interactive tool.
(1) Sulfate effect: for each substrate and time point, sulfate-amended microcosms were compared to their corresponding non-stimulated controls (i.e., did the tested OTU respond to external sulfate stimulation).
(2) Substrate effect: substrate-amended microcosms were compared to the nosubstrate-controls. This was done separately for each time point and substrate under either sulfatestimulated or unstimulated conditions (i.e., under which substrate-amendments did the tested OTU respond).
(3) Temporal effect: within each treatment, the 26 d and 50 d time points were compared to the 5 d time point (i.e., did the relative 16S rRNA (gene) abundance of the tested OTU increase in time). Pairwise comparisons consisted of three replicates each, with the exception of the propionate-and sulfate-amended microcosms, where one replicate was excluded from all analyses because of its inconsistent sulfate turnover ( Figure 1). Before each differential abundance test, OTUs were filtered. To be included, OTUs had to be reliably detected, i.e., be represented by more than 10 reads in at least two out of six samples. This resulted in the statistical testing of 1491 out of 7435 OTUs.
OTUs were considered responsive if all three types of differential abundance tests were significant (FDR-corrected p-value <0.05). OTUs significantly more abundant under sulfate-stimulated than under unstimulated conditions were denoted as sulfate-stimulated responders. Inversely, OTUs significantly more abundant in the no-sulfate microcosms were denoted as sulfate-deterred responders. Each response was always assigned to a certain substrate (significant substrate effect) and could be an intermediate response (after 25 days), and late response (after 50 days), or both (significant temporal effect). We also tested for sulfate-stimulation and time-dependent response with endogenous substrates by omitting the substrate effect test in the no-substrate-control microcosms.
For snapshot analyses of environmental 16S rRNA using only one time point, it has been pointed out before that the amount of ribosomal RNA is not only influenced by current conditions but also by life history (past) or ecological strategy (future) (Blazewicz et al., 2013). Our experimental design largely excluded such effects. Since we tested for increases in ribosome abundance over time, a legacy from past activity expressed as a constantly high ribosome content (for example, Sobek et al., 1966) could be excluded. An increase of the ribosome content by inactive populations entering starvation and dormancy (for example, Sukenik et al., 2012) could be excluded as well by omitting OTUs that responded unspecifically, i.e., irrespective of sulfate and/or specific substrate additions. qPCR targeting the genus Desulfosporosinus and Bacteria/Archaea qPCR assays targeting the 16S rRNA genes of the genus Desulfosporosinus and most Bacteria and Archaea were performed as described previously (Pester et al., 2010). Primer 1389F was modified to cover more Archaea (5′-TGY ACA CAC CGC CCG T-3′). In addition we analyzed the 16S rRNA cDNA of the genus Desulfosporosinus by reverse transcription using the primer DSP821R, followed by qPCR (Pester et al., 2010). To estimate Desulfosporosinus ribosome and genome abundance per cm³ of fresh soil, 16S rRNA (gene) copies per ng of extracted DNA/RNA were normalized against the amount of soil used for nucleic acid extraction, based on following assumptions: 100% nucleic acids extraction efficiency, 100% reverse transcription efficiency, Desulfosporosinus rrn operon copy number of 8.75 (Pester et al., 2012a), a soil bulk density of 0.29 g cm⁻³ (Goldberg et al., 2008), and a water content of 78% (this study). Total Bacteria and Archaea gene copy numbers per cm³ of fresh soil were calculated analogous.

qPCR of species-level dsrA OTUs
Abundances of the two dsrAB OTUs 1 and 2 in the incubations were measured by using established dsrA-targeted qPCR assays (Steger et al., 2011). Gene copy numbers per cm³ of fresh soil were calculated analogous to the qPCR assay for Desulfosporosinus.

Global distribution and abundance of Desulfosporosinus in wetlands
We gained insights into the global distribution and abundance of Desulfosporosinus species by identifying sequences in public databases with ≥97% 16S rRNA identity to the Desulfosporosinus OTUs 0051 and 0228. We used blastn (Camacho et al., 2009) and the 253 nt long OTUs sequences to query GenBank (Benson et al., 2013). To identify closely related sequences in the short read archive (SRA) (Leinonen et al., 2011), we initially choose proxy sequences for each OTU that were >1400 nt long and 100% identical over the V4 region of the 16S rRNA (GenBank accession numbers EU981233 and KJ650684, respectively). This made it possible to survey sequence information of non-V4 region amplicon data sets. We applied these proxy sequences in the integrated microbial NGS platform (http:/ /www.imngs.org/) to retrieve all SRA samples that harbored ≥97% identical reads. The metadata of all GenBank and SRA matches was manually screened to identify all sequences that originated from wetlands (including rice paddy and permafrost soils). If sufficient metadata was available, we classified source environments into peatland, permafrost, rice paddy, or other soils. We distinguished between sequences obtained directly from soils (i.e., pristine soils, rice paddy fields, or restored wetlands) or from soils that were manipulated (for example, microcosms incubations, SIP enriched samples, or enrichments/isolates). Locations of samples without available geographical coordinates were estimated based on other information given (for example, country/region of origin). For all demultiplexed samples, we also retrieved the relative abundance of the 16S rRNA gene reads that were closely related (≥97% sequence identity) to the two Desulfosporosinus OTUs.

CO₂ and CH₄ production
We used a second set of microcosms that were identical in incubation conditions to the setup used for the molecular analyses to monitor total CO₂ and CH₄ production under the various substrate and sulfate scenarios. The amount of total organic carbon mineralized (sum of produced CO₂ and CH₄) in the individual microcosms as well as the fraction of organic carbon added by external substrates is presented in Supplementary Table S1. Surprisingly, added substrates did not result in a surplus amount of organic carbon mineralization. A possible explanation could be that microorganisms using either internal or external substrates competed for the same limiting factor, which was however not the electron donor. All provided substrates are most likely easier to degrade than internal substrates, which have to be extracted from humic matter. Short-chained fatty acids and lactate can be only degraded by anaerobic respiration (for example, using nitrate, Fe(III), or sulfate as electron acceptor) or by secondary fermenters, which have to cooperate typically with methanogens and have a lower energy yield as compared to anaerobically respiring microorganisms (Schink and Stams, 2006). Therefore, a plausible limiting factor could be the various electron acceptors, especially Fe(III). It has been shown that iron reduction is an important anaerobic carbon mineralization process in the studied peatland (Küsel et al., 2008). In addition, iron reduction is a primarily surface-bound process posing a potentially limiting factor to the accessibility of this electron acceptor.
In control incubations without sulfate amendment, methanogenesis accounted for 0.1-0.2% of the total carbon flow during organic matter mineralization over the total incubation time of 27 days (Supplementary Table S1), which is again most likely explained by the presence of active Fe(III)reducing microorganisms (Küsel et al., 2008). In line with these low activities, methanogenic archaea did not respond significantly by increasing their 16S rRNA or 16S rRNA genes over time in these methanogenic controls. A primer bias seems unlikely because the used primers cover 87% of all Euryarchaeota sequences in the ARB-SILVA database (Methanobacteriales 93%, Methanococcales 84%, Methanocellales 89%, Methanomicrobiales 91%, Methanosarcinales 89%, Methanopyrales 0% -contains only five sequences, Thermoplasmatales 85%; Quast et al., 2013) and we detected five OTUs within the Methanosarcinales and five OTUs within the Methanocellales.

Mock community
A mock community consisting of the activated sludge clones H42, H29, H28, H13, and H44 (GenBank accession numbers AF234715, AF234692, AF234749, AF234737, and AF234743, respectively) was mixed at a ratio of 76%, 18%, 5%, 0.7%, and 0.09% relative abundance, respectively (Herbold et al., 2015), and included in one of the Illumina MiSeq sequencing runs. Clones were mapped to the representative OTU sequences using blastn and applying a threshold of 97% sequence identity to accommodate potential sequencing errors, resulting in an unambiguous assignment. Log-transformed expected and measured abundances for each mock community clone were compared by fitting a linear regression (R Core Team, 2015), which resulted in a significant slope of 1.0 (p-value=0.01, R²=0.90, Supplementary Figure S4a) and suggests that relative abundance shifts can be reliably measured by our amplicon sequencing approach. The five most abundant OTUs were assigned to the original mock community clones (Supplementary Figure S4b). The sixth most abundant OTU with only 7 reads was assigned to Sulfate reduction-associated bacteria in peat soil Supplementary Information Escherichia/Shigella and was most likely a contamination (Salter et al., 2014). For comparison, the least abundant mock community clone had 119 reads.

Phylogenetic analysis of Desulfosporosinus and Telmatospirillum OTUs
Phylogenetic analysis of the two responsive Desulfosporosinus OTUs and the three responsive Telmatospirillum OTUs verified their placement within the genera Desulfosporosinus (Supplementary Figure S7a) and Telmatospirillum (Supplementary Figure S7b), respectively. The pairwise similarity of the two Desulfosporosinus OTUs was 96.8% and thus within the genus-level threshold of 94.5% (Yarza et al., 2014). Interspecies similarity of Desulfosporosinus pure cultures reaches a minimum of 92.5%. OTU0051 is identical (over the V4 region) with eight clones from the heavy DNA fraction of sulfateamended SIP incubations of the same peatland (Pester et al., 2010). In contrast, OTU0228 had a maximum sequence similarity of 98% to clones from the SIP study.
Pairwise sequence similarities between the three Telmatospirillum OTUs varied between 94.5% and 96.8%. OTU0062, which responded stronger in sulfate-stimulated incubations as compared to unstimulated control incubations, was phylogenetically clearly separated from Telmatospirillum OTU0029, which showed the opposite response. In addition, OTU0062 was identical (over the V4 region) to 7 SIP clone sequences identified previously in the heavy DNA fraction of sulfate-amended as well as unamended SIP incubations of the same peatland (Pester et al., 2010). OTUs 0029 and 0273 had a maximum sequence similarity of 95% and 96% to clones from the SIP study, respectively.

qPCR analysis of abundant species-level dsrAB OTUs
Previous studies of this peatland identified a large diversity of reductive-type dsrAB sequence variants (coding for subunit A and B of the dissimilatory sulfite reductase), which serve as functional marker genes for SRM but are also present in some sulfite reducers, degraders of organosulfonates, and syntrophic microorganisms that likely lost the capability to reduce sulfate (Müller et al., 2015). In the studied peatland, so far 53 different dsrAB OTUs at the approximate species level were identified with most of them affiliating to deep-branching lineages with no cultured representative (Pester et al., 2012b). Several of these unidentified dsrAB OTUs represent autochthonous and abundant members of the microbiota in this peatland (Steger et al., 2011). We monitored the abundance of two of these OTUs (designated dsrAB OTU 1 and 2; Steger et al., 2011) in our incubations by qPCR. However, there was no clear response to sulfate amendment under any of the tested substrate scenarios (Supplementary Figure S6c). These results are in line with findings from a previous DNA-SIP study, which could not detect selective ¹³C-labeling of these novel dsrAB variants in sulfate-treated microcosms (Pester et al., 2010). These novel deep-branching dsrAB lineages either derive from non-SRM or the tested substrates might not have been in the metabolic range of SRM harboring these novel dsrAB variants.
Sulfate reduction-associated bacteria in peat soil Supplementary Information

Supplementary Tables
Supplementary Table S1 Amount of mineralized (CO₂ + CH₄) and supplemented carbon after 27 days as observed in parallel incubations (complete time series data shown in Supplementary Figure S2c). Mean and 95% confidence intervals are given in µmol.
Percentages are relative to total mineralized carbon. Overview of the experimental setup. Sulfate and substrate amendment regime is indicated with colored bars. Supernatant sampling for quantification of sulfate and individual substrates is indicated with black bars. Arrows indicate microcosm sampling and molecular analysis workflows or pH measurements. Additional parallel incubations for measurements of CO₂ and CH₄ are not shown. These were set up with the same sulfate and substrate amendment regime but were stopped after 27 days and gas from the headspace was sampled at days 0, 4, 11, 18, 25, and 27 to follow CO₂ and CH₄ production.
Sulfate reduction-associated bacteria in peat soil Supplementary Information

Yes No
Temporal changes in relative ribosome and genome abundance of selected responsive OTUs (i.e. mentioned in the main text) in incubations with different substrates. OTUs are sorted by their number. Solid lines and symbols depict sulfatestimulated microcosms whereas dashed lines and open symbols depict controls without additional sulfate. Diagonal crosses depict the abundance in the native soil, which was sampled from parallel peat soil subsamples and plotted at day 0 in all panels. Halos around symbols indicate significantly higher abundance in the unstimulated or sulfatestimulated microcosms, respectively, as compared to their respective controls and day 5. Data points drawn directly on the x axis represent relative abundance values of zero. Maximum likelihood 16S rRNA trees of the genera Desulfosporosinus (a) and Telmatospirillum (b). Representative V4 amplicon sequences (bold) of identified Desulfosporosinus and Telmatospirillum OTUs were added to the tree by using ARB's Parsimony Interactive tool. Solid circles indicate ≥90% and open circles indicate ≥70% bootstrap support. The bars represent 0.05 inferred sequence divergence. GenBank accession numbers of reference sequences are given in parentheses. All SIP clones (Pester et al., 2010), clone SBIb49 (Lüdecke et al., 2010), and clone P6K11f (unpublished) are from the same peatland as the OTUs identified in this study.