Quantifying population-specific growth in benthic bacterial communities under low oxygen using H218O.

The benthos in estuarine environments often experiences periods of regularly occurring hypoxic and anoxic conditions, dramatically impacting biogeochemical cycles. How oxygen depletion affects the growth of specific uncultivated microbial populations within these diverse benthic communities, however, remains poorly understood. Here, we applied H218O quantitative stable isotope probing (qSIP) in order to quantify the growth of diverse, uncultured bacterial populations in response to low oxygen concentrations in estuarine sediments. Over the course of 7- and 28-day incubations with redox conditions spanning from hypoxia to euxinia (sulfidic), 18O labeling of bacterial populations exhibited different patterns consistent with micro-aerophilic, anaerobic, facultative anaerobic, and aerotolerant anaerobic growth. 18O-labeled populations displaying anaerobic growth had a significantly non-random phylogenetic distribution, exhibited by numerous clades currently lacking cultured representatives within the Planctomycetes, Actinobacteria, Latescibacteria, Verrucomicrobia, and Acidobacteria. Genes encoding the beta-subunit of the dissimilatory sulfate reductase (dsrB) became 18O labeled only during euxinic conditions. Sequencing of these 18O-labeled dsrB genes showed that Acidobacteria were the dominant group of growing sulfate-reducing bacteria, highlighting their importance for sulfur cycling in estuarine sediments. Our findings provide the first experimental constraints on the redox conditions underlying increased growth in several groups of "microbial dark matter", validating hypotheses put forth by earlier metagenomic studies.

Introduction

Benthic microbial communities living in estuarine ecosystems play an important role in global biogeochemical cycles, because they drive organic matter decomposition, nutrient regeneration, and influence water column dissolved O2 concentrations [1, 2]. Benthic oxygen depletion is typical in estuarine habitats, where degradation of organic matter is 100–1000 times higher than corresponding values in the water column [3]. The resulting hypoxia impacts both fisheries by increasing fish mortality [4], as well as climate by facilitating increased fluxes of the greenhouse gas nitrous oxide [5].

Estuarine sediments have complex microbial communities composed primarily of uncultured lineages catalyzing aerobic, micro-aerophilic, and anaerobic metabolic pathways that impact carbon, nitrogen, and sulfur cycling [6]. Quantifying growth of specific bacterial populations is challenging, yet critical for understanding of ecosystem resilience and response to change [7]. The structure of microbial communities in estuarine habitats can exhibit resistance to environmental perturbation [8], but the levels of activity within populations can change dramatically in response to changing nutrients and oxygen levels, with clear impacts on biogeochemical cycles [9].

The activity of growing microbial populations in environmental samples can be quantified using quantitative DNA-stable isotope probing (qSIP) with H218O as a passive tracer [10]. Oxygen atoms from H218O are incorporated into DNA during genome replication, which when combined with quantitative PCR (qPCR) and high-throughput sequencing of 16S rRNA genes can be used to quantify activity of growing populations within complex microbial communities from environmental samples [11, 12]. Relative to energy-rich 13C-labeled carbon substrates, labeled water is a passive tracer of cell growth, whereby DNA replication generates a new DNA strand that will contain 18O atoms in the presence of labeled water [12]. The amount of 18O incorporated into the total DNA pool is correlated with growth rates [11-13], showing that 18O labeling occurs primarily during growth via DNA replication [12]. The degree of atomic incorporation can then be used as a quantitative proxy for growth [13]. qSIP with H218O has been applied to terrestrial habitats [14-16], including freshwater sediments [17]. But, to our knowledge, H218O qSIP has not yet been applied to quantify growth in estuarine sediments or under anoxic conditions.

We used H218O qSIP [12] to quantify for the first time population-specific growth dynamics in benthic bacterial communities in response to changing redox conditions. This allowed us to test hypotheses regarding the potential physiology of several groups of uncultivated microbial groups put forth by earlier metagenomics studies, which suggested an adaptation to low oxygen and anoxic aquatic environments [6, 18]. Our results showed that hypoxia and euxinia selected for specific phylogenetic groups of uncultivated bacteria whose metabolic activity was increased, providing evidence of their optimal redox conditions for growth. Notably, establishment of benthic anoxia coincided with increased growth from numerous uncultivated groups of sulfate-reducing bacteria (SRB) that were dominated by the Acidobacteria, which should thus be more closely considered as an important SRB group impacting sulfur cycling in estuarine sediments. Our experimental findings validate prior hypotheses put forth by metagenomics studies indicating micro-aerophilic and anaerobic lifestyles for many groups of “microbial dark matter (MDM)”.

Materials and methods

Sampling

Surface sediment samples were collected in July 2016 from 1 m water depth in Sage Lot pond, a coastal lagoon connected as a sub-estuary to Waquoit Bay (Cape Cod, Massachusetts). Sage Lot pond is a small (surface area 0.17 km2) shallow (ca. 2 m maximum depth) lagoon surrounded by dense vegetation including salt marshes and seagrasses [19, 20]. Sage Lot pond exhibits phytoplankton chlorophyll concentration up to 90 mg L-1 when nitrogen inputs increase [19]. These eutrophic conditions lead to frequent benthic anoxic events [19].

Experimental setup

We added sea salts (30 mM MgCl2 ٠6H2O, 16 mM MgSO4 ٠7H2O, 2 mM NaCO3, 10 mM KCl, 9 mM CaCl2, 450 mM NaCl) to 99% H218O (Sigma-Aldrich, St. Louis, MO, USA) in order to create 18O-labeled artificial seawater (ASW). As a control, ASW was also created using diethyl pyrocarbonate (DEPC)-treated (sterile, nuclease free) water. Both waters were filter sterilized (0.2 μm). One milliliter of either 18O-labeled or -unlabeled (control) ASW was added to 2 g of wet surface sediment from Sage Lot Pond in 20 mL sterile glass vials containing sterile oxygen sensor spots (PreSens Precision Sensing). The oxygen sensor spot was positioned at the sediment-seawater interface to measure benthic O2 concentrations, and additional sensor spots were placed in the headspace of two flasks to measure gaseous O2 levels throughout the incubation. Incubations were set up in biological triplicate for each timepoint (7 day and 28 day). The water content of the sediments was 15% (±1%), and thus the final concentration of H218O in the H218O incubations was roughly 66%. After addition of labeled and unlabeled ASW, flasks were crimp sealed with gas tight gray butyl rubber stoppers. All flasks contained ca. 15 cm of oxygenated headspace and were incubated in the dark for 7 and 28 days at 8 °C. Dissolved oxygen was measured noninvasively using a fiber optic oxygen sensor (PreSens, Regensburg Germany) ca. 0.5 cm above the sediment-water interface as described previously [21]. Oxygen measurements were also performed on autoclaved sediments as a killed control. DNA from the samples was extracted and quantified from the replicate incubations at the beginning (T0), 7 days, and 28-day timepoints as described previously [22].

Density gradient centrifugation and gradient fraction

DNA samples were prepared for density gradient centrifugation according to previously defined protocol for qSIP [23]. In brief, density gradient centrifugations were carried out in a TLN-100 Optima MAX-TL ultracentrifuge (Beckman Coulter, Brea, CA, USA) near-vertical rotor at 18 °C for 72 h at 165,000 × g. In all, 50 µl of DNA spanning from 0.5 to 1.5 µg [24] was added to a solution of cesium chloride (CsCl) and gradient buffer (0.1 M Tris, 0.1 M KCl and 1 mM EDTA) in order to achieve a starting density of 1.70 g mL-1 in a 3.3-mL polyallomer OptiSeal tubes (Beckman Coulter, Brea, CA, USA). After ultracentrifugation, the density gradients were fractionated into 15 equal fractions of 200 µl from the bottom of polyallomer OptiSeal tubes by using a syringe pump and fraction recovery system (Beckman Coulter, Brea, CA, USA). The density of these fractions was measured with an AR200 digital refractometer (Reichert Analytical Instruments, Depew, NY, USA). DNA was precipitated from the fractions using two volumes of polyethylene glycol with 2 µl (10 mg mL-1) glycogen and precipitated overnight at room temperature. DNA was pelleted by centrifugation (13,000 × g; 40 min), washed with 70% ethanol, and resuspended with 30 µl molecular-grade (DEPC-treated) water. DNA was quantified fluorometrically using a Qubit 4 fluorometer (Thermo Scientific).

qPCR, 16S rRNA gene, and dsrB gene sequencing

Universal primers targeting the V4 hypervariable region of 16S ribosomal RNA (rRNA) genes were used in qPCR to determine density shifts of key genes (16S and dsrB) for each incubation. We used a version of the 16S rRNA gene 515F primer with a single-base change (in bold) to increase the coverage of archaea (515F-Y, 5′-GTGYCAGCMGCCGCGGTAA-3′; [25]). All qPCR reactions were carried out as described previously using the Eppendorf EpMotion 5070 pipetting robot that has <5% technical variation [23]. Each density fraction was also screened using qPCR for SRB with primer pairs targeting the dissimilatory sulfite reductase β-subunit genes (dsrB) according to a previously published assay [26-28]. We chose to focus on the dsrB gene because a large database exists for dsrB sequences recovered from environmental samples that we could compare our data against [29] (www.microbial-ecology.net/download). qPCR standards consisted of 10-fold dilution series of the genes of interest that were PCR amplified from the sample at 40 cycles using the same primers. Prior to creating the dilution series, the correct size of amplified standard was confirmed via gel electrophoresis, gel extracted, and quantified with a Qubit. Reaction efficiencies in all qPCR assays were between 90 and 110% with r2 > 0.98 for the standards. dsrB amplicons were cloned and sequenced via sanger sequencing from density fractions at the 28-day incubation timepoint that exhibited 18O labeling, in the density range 1.70–1.71 g mL-1.

Two 16S PCR amplicons from each density fraction (technical replicates to reduce PCR bias) were pooled and sequenced on the Illumina MiniSeq as described previously [30]. To account for the influence of contamination, we included barcoded aerosol (laboratory dust) and kit reagents (DNA extraction blanks) samples.

dsrB amplicons were gel extracted and cloned using the TOPO TA cloning kit (Invitrogen, Life Sciences) according to the manufacturer instructions. A total of 132 clones were picked, the insert size confirmed via PCR, and those clones having the correct dsrB size (89 clones) were Sanger sequenced bidirectionally. The forward and reverse Sanger reads were used to create dsrB contig sequences in CodonCode Aligner version 8.0.2 (CodonCode Corporation, MA, USA).

Bioinformatic analysis

The Illumina reads were quality trimmed and assembled using USEARCH version 10.0.240 with the default parameters [31] resulting in 6.8 million quality checked V4 reads. Reads were then de novo clustered at 97% identity using UPARSE; OTUs represented by a single sequence were discarded [32]. Taxonomic assignments were generated by QIIME 1.9.1 [33] using the implemented BLAST method against the SILVA rRNA gene database release 132 [34]. After that, only operational taxonomic units (OTUs) >12 sequences in total in each replicate for the control and SIP-labeled fractions were selected for further study [23, 35]. OTUs detected in the contaminant datasets were removed from all downstream analysis if the total number of sequences in the contaminant sample was greater than the experimental sample. Working with this “cleaned” dataset, 598 OTUs and 523 OTUs from 7 days and 28 days incubations were used for downstream analyses.

Observed excess atom 18O fractions (EAFs) were calculated for each taxon as described previously [13] using a qSIP workflow embedded in the HTS-SIP R package [36]. To calculate the bootstrap confidence intervals (CI) for significant isotopic incorporation, bootstrap replicates (n = 1000) were run with the HTS-SIP R package [36]; an OTU was considered as having isotopic incorporation (true positive) if the lower CI was >0 [13].

Phylogenetic analyses were performed in SeaView [37] following alignment with MUSCLE [38]. Maximum likelihood (ML) with selected substitution model as general time reversable (GTR) was performed with PhyML version 3.0 [39]. dsrB gene translation was performed using EMBOSS Transeq [40]. W-IQ-TREE (http://iqtree.cibiv.univie.ac.at) was used to find the best model using Model Finder [41, 42], which resulted in LG + G4 model. Trees were visualized and edited using iTOL [43]. Statistical analyses and plots were performed using R.Studio Version 3.3.0 [44]. Blomberg’s K [45] and Pagel’s λ [46] tests for significantly non-random phylogenetic distributions of growing patterns from qSIP were calculated on all OTUs (labeled and unlabeled) using the phylosignal R package [47]. Both indices test species’ traits under a Brownian motion model (BM) of trait evolution; that is whether or not the distribution of traits across different phylogenetic groups is random or non-random. The BM assigns a 0 value to indicate phylogenetic independence (random phylogenetic distribution of traits) and values close to 1 for a strong phylogenetic signal (non-random phylogenetic distribution of traits) [48-50].

The sizes of the growing and dying fractions of each population, and their rates of change, was calculated using a model developed to determine population growth and mortality rates from 18O-qSIP data [51]. For all calculations, bootstrap resampling of replicates within each treatment was used to reproduce the uncertainty and 90% CIs were estimated. These calculations were performed in R using the code at https://bitbucket.org/QuantitativeSIP/qsip_repo [51]. Sequence data were entered in the NCBI Short Read Archive under BioProject ID PRJNA498588.

Results

Dissolved oxygen measurements

The sediment-water interface was well-oxygenated at the beginning of the incubation (70% atm. saturation), which declined exponentially during the first week until reaching 0% atm. saturation after day 5 (Fig. 1a). This was not observed in the killed control, showing that the rapid drawdown of benthic O2 was due to respiration. Small fluctuations in the oxygen measurements in the killed control were likely due to temperature fluctuations of the incubator itself (±1 °C), since the non-invasive fiber optic oxygen sensor spots are temperature sensitive [21]. By the end of the experiment (28 days), the sediments had turned from a gray-brown color to black, indicating the presence of iron–sulfur minerals (e.g., FeS and FeS2), and upon opening the vials sulfide could be smelled. Thus, while the headspace contained oxygen at the beginning of the experiment, the sediments had become anoxic and sulfidic (euxinic) by the end of the incubation period.

Fig. 1

a Benthic O2 concentrations during the incubation and killed (autoclaved) control. b, c Quantification of 16S rRNA gene copies across CsCl density gradient fractions after 7 (b) and 28 days (c). 18O water replicates are represented by blue solid lines with triangles and unlabeled replicates (control) are represented by red dashed lines with circles. The y axis represents the relative abundance of 16S rRNA genes quantified with qPCR, normalized to maximal abundance across all density fractions

qSIP of 16S rRNA genes

18O labeling of 16S rRNA genes was observed at 7 and 28 days, with 16S rRNA genes at 28 days exhibiting a higher degree of labeling compared with 7 days (Figs. 1b, c). The decrease in oxygen (Fig. 1a) was mirrored by a decrease in 16S rRNA gene copies: the total number of 16S rRNA gene copies per gram wet sediment decreased during the first week, from 2.5 ( ± 0.12) × 108 at T0 to 1.4 ( ± 0.09) × 108 at 7 days, and then decreased further after 28 days to 1.1 ( ± 0.04) × 108. This indicated net microbial death with a fraction of the community maintaining growth during the incubation.

The composition of microbial populations at 7 and 28 days were markedly different (analysis of variance; F = 3991, p < 0.001), but dominated by the same phyla in nearly equal proportion (Figure S1). In total, 443 OTUs were detected at both timepoints, whereas 235 OTUs were detected at only a single timepoint (Figure S1). In all, 128 OTUs were 18O labeled after 7 days, which increased to 395 OTUs after 28 days (Figure S1). Of the 18O-labeled OTUs at day 7, Bacteroidetes were the most abundant taxa with 67 OTUs (49.4% of the 18O-labeled OTUs), followed by 39 OTUs affiliated with Proteobacteria (42.1% of the 18O-labeled OTUs) (Fig. 2 and S1). On the other hand, the 395 18O-labeled OTUs at day 28 were comprised mostly of Proteobacteria (157 OTUs, 61.4% of 18O-labeled OTUs), followed by Bacteroidetes (94 OTUs; 18.1% of 18O-labeled OTUs) and Planctomycetes (53 OTUs; 8.71% of 18O-labeled OTUs) (Fig. 2 and S1).

Fig. 2

OTU-specific shifts in the median atom fraction excess (18O) of OTUs with 90% confidence interval (CI). OTUs were colored by phylum. OTUs that do not have a 90% CI overlapping with 0 are considered to be 18O labeled

Although the scope of this study is to determine the growing microorganisms, non-growing cells due to dormancy and or slow growing cells are represented in our results as those OTUs that did not become labeled. Accordingly, a total of 248 OTUs constituted non-growing microbes (unlabeled) in both incubations (Fig. 3). The non-growing or slow growing bacterial groups were dominated by OTUs affiliated with sulfate-reducing lineages in the Deltaproteobacteria (42 OTUs), Planctomycetes (38 OTUs), Epsilonbacteraeota (7 OTUs), and Spirochaetes (9 OTUs) (Fig. 3).

Fig. 3

Phylogeny of bacterial taxa detected and their extent of 18O labeling at 7 and 28 days. The inner circles correspond to excess atom fraction (EAF) 18O values of labeled taxa after 7 and 28 days of incubation. The outer ring of the heatmap represents the EAF change between the timepoints. The numbered and colored circles represent the growth mode of the OTUs and the legend defines the growth mode of categories explained in result section. Bottom panel shows phylogenetic signal tests (Blomberg’s K and Pagel’s λ) and corresponding p-values for labeled taxa and growth mode categories

Growing bacterial OTUs affiliated with MDM [18] candidate phyla Aegiribacteria (1 OTU, 0.08 EAF) and Patescibacteria (1 OTU, 0.04 EAF) were detected only in the 7-day incubations (Fig. 2 and Table 1). Labeled OTUs affiliated with MDM groups after 28 days included Latescibacteria (10 OTUs, 0.24 ± 0.04 EAF) and Calditrichaeota (1 OTU, 0.24 EAF). The only MDM group that was 18O labeled at both 7- and 28-day timepoints was the candidate phylum WPS-2 (Fig. 2 and Table 1).

Table 1

A summary of growth by uncultivated taxa, previously designated as “microbial dark matter”, under various redox conditions

Group	Growth mode	EAF values (number of OTUs)	r-Value (net population growth rate)	b-Value (rates of reproduction)	d-Valuea (mortality rate)
Aegiribacteria	Micro-aerophilic	0.08 (n = 1)	-0.202	0.081	-0.288
Latescibacteria	Aerotolerant anaerobic	0.013 ± 0.008 (n = 10)	-0.227	-0.038	-0.191
Gracilibacteria	Micro-aerophilic	0.04 (n = 1)b	-0.172	-0.041	-0.136
Calditrichaceae	Aerotolerant anaerobic	-0.007 (n = 1)	-0.214	-0.018	-0.197
WPS-2	Facultative anaerobic	0.1 (n = 1)b	-0.339	0.079	-0.421
28 days of incubation
Gracilibacteria	Aerotolerant anaerobic	0.17 ± 0.073 (n = 2)b	-0.136	0.023	-0.171
Latescibacteria	Aerotolerant anaerobic	0.19 ± 0.023 (n = 14)c	-0.373	0.019	-0.392
WPS-2	Facultative anaerobic	0.2 (n = 1)b	-0.422	0.014	-0.436

ar, b, and d values correspond to the rates per day (d-1)

bSignificantly growing microorganisms based on qSIP (i.e., lower boundary of bootstrap is >0)

cMost of the Latescibacteria grew (10 OTUs out of 14)

The change in oxygen concentrations over the course of the experiment allowed us to group OTU growth into five categories based on the pattern of 18O labeling at 7 days (micro-oxic conditions) and 28 days (anoxic conditions) (Fig. 3) [1]. Micro-aerophilic growth was defined as 18O-labeled OTUs detected only at day 7, and not at 28 days [2]. Anaerobic growth was defined as 18O-labeled OTUs detected only after establishment of euxinic conditions at day 28, and not at 7 days [3]. Facultative anaerobic growth was defined as 18O-labeled OTUs overlapping between both timepoints [4]. Aerotolerant anaerobic growth was defined as OTUs detected at both timepoints, but only 18O labeled during anoxic and sulfidic conditions [5]. Dormant or slow growing microorganisms were defined as OTUs that were not 18O labeled at either timepoint.

18O labeling of dsrB genes

Bacterial dsrB genes exhibited 18O labeling only after development of euxinic conditions sampled at 28 days of incubation, with peak DNA buoyant density (BD) of 1.71 (±0.008) g mL-1, which was greater than the control where the peak was 1.686 (±0.003) (Fig. 4a). This corresponds to an increase in the atomic enrichment percentage of >20%, which is typically regarded as the threshold for significant isotopic labeling [52]. In total, 89 dsrB sequences were obtained by molecular cloning from the density fractions that exhibited peak 18O labeling (Fig. 4a). Most of the 18O-labeled dsrB sequences were affiliated with novel groups of uncultivated Acidobacteria SRB (61 sequences, 82% of total), whereas the remainder of 18O-labeled dsrB sequences (13 sequences, 18% of total) were related to Deltaproteobacteria SRB (Fig. 4b). The 18O-labeled Acidobacteria dsrB genes include a clade of four dsrB sequences with close relation to a novel clade of SRB originally described as “novel dsrB Group IV” from the Guaymas hydrothermal vent [53], raising the possibility that this group of Acidobacteria contains both thermophilic and mesophilic SRB.

Fig. 4

a Quantification of dsrB gene copies across CsCl density gradient fractions after 28 days. 18O water replicates are represented by blue solid lines with triangles and unlabeled replicates (control) are represented by red dashed lines with circles. The y axis represents the relative abundance of dsrB genes quantified with qPCR normalized to maximal abundance across all density fractions. b Phylogenetic tree of 18O-labeled dsrB genes including their most similar sequences from the NCBI nr database, bold sequences indicate those from this study. Collapsed clades (triangles) show the number of 18O-labeled dsrB gene sequences contained within the clade. Black circles at nodes represent bootstrap support of 90%, gray circles represent bootstrap support from 70 to 90%, and white circles represent bootstrap support from 70 to 50%

Growth and death dynamics of 18O-labeled populations

After 7 days, most genera exhibited gross reproduction and three OTUs exhibited significant net production (defined as 90% CI in the growth/death model not overlapping 0), affiliated with Sva0081 sediment group (Desulfobacteraceae), Desulfobacterium catecholicum, and SB-5 family of Bacteroidetes (Fig. 5). The rate of mortality per genus was higher at day 28 than day 7, indicating that establishment of euxinic conditions caused the majority of cells per genus to die faster than they grew (Fig. 5). However, despite the higher net mortality rates, many exhibited relatively low positive gross reproduction rates (Fig. 5) indicating that a smaller proportion of individuals per population were actively growing. This is consistent with the result that most of the OTUs were 18O labeled at this timepoint. We also calculated the whole-assemblage turnover estimated via qSIP using the developed model of Koch et al. [51]. The seven-day incubation had an average community turnover value of 0.28 d-1 (90% CI: 0.219–0.36 d-1), whereas 28 days of incubation had an average community turnover value of 0.371 d-1 (90% CI: 0.315–0.476 d-1).

Fig. 5

Population growth rates (r) at 7 and 28 day timepoints. After 28 days, all taxa exhibited mortality rates greater than reproduction rates. Points indicate bootstrapped medians and bars show 90% confidence intervals for each OTU. The positive net production rates here are considered as statistically significant increase in the populations if the bootstraps are not crossing the zero. Open circles: rates of reproduction, black filled circles: mortality rates, blue filled circles: net reproduction rates. The only groups that exhibiting significant net production at day 7 are indicated with asterisks (see legend)

Discussion

Dissolved oxygen has declined in ocean water in the past five decades due to the increase in global temperature [54], resulting in the expansion of oxygen minimum zones (OMZs) in the oceans [55]. In coastal settings, increased human activity such as high fertilizer use has caused widespread eutrophication and recurrent bottom water anoxia that influences the benthos [19, 56, 57]. The effects of such benthic anoxia on the growth of diverse bacterial populations driving elemental cycles is poorly understood. We used 18O-qSIP to quantify for the first time to quantify the growing bacterial populations in response to benthic oxygen depletion, including many groups for which there currently exists no cultured representatives.

Assessing effects of incubation conditions

Although the sediment surface became oxygen depleted after 5 days of incubation, there was an oxygenated headspace, and O2 measurements of the headspace confirmed oxygen throughout the 28-day incubations at concentrations of 90–80% atmospheric saturation (data not shown). The gradual depletion of oxygen at the benthic interface during the first week thus indicates a shoaling of the oxic–anoxic transition zone within the sediments into the overlying ASW. Benthic anoxia is a common feature of the sampled environment in Waquoit Bay (Cape Cod, USA), where increased nitrogen input to the watershed through atmospheric deposition, fertilizer, and wastewater has led to an increase in primary productivity and recurrent bottom water anoxia [19, 56]. This phenomenon is also known to occur in the marine environment, for example, in the Benguela upwelling system where summer time water column stratification results in a shoaling of the sediment oxic–anoxic transition zone several meters upwards into the bottom waters and where sulfide accumulates to high levels on the continental shelf [58]. Thus, our experimental conditions are relevant to naturally occurring benthic habitats in estuaries and continental shelf settings that are especially prone to the development of anoxia. Dormancy by anaerobic bacteria during periods of oxygenation at the sediment surface may explain how they survive until favorable anoxic conditions are re-established, for example, after burial deep below the subseafloor in anoxic sediments [59].

Given the detection limit of the fiber optic oxygen measurements (ca. 0.5 % atmospheric saturation), we cannot rule out the presence of trace amounts of dissolved oxygen within the overlying seawater or sediment. Nevertheless, the appearance of black color throughout the sediments and the strong smell of sulfide at the end of the 28-day incubation strongly indicated anoxic conditions in the sediments by the end of the experiment. The strong smell of sulfide suggests that the black color reflected formation of amorphous iron–sulfur compounds, caused by a reaction between oxidized Fe(III) with HS- [60]. As we sampled the entire sediment slurry for our qSIP analysis, we likely sampled a mixture of both anoxic and micro-aerophilic habitats that were present as a steep redox gradient within the flasks. Therefore, the active microbes detected at 7 days probably represent micro-aerophilic bacteria, whereas active microbes sampled at 28 days were living under predominantly anoxic and sulfidic conditions.

It is likely that anoxic microsites likely developed within the sediment in the first 7 days, potentially complicating the categorization of aerobic, facultatively anaerobic, and aerotolerant growth in OTUs detected at this timepoint. Anaerobic SRB can remain active in oxidized marine sediments within anoxic microniches [61], but some SRB are aerotolerant since they can grow (albeit slightly) in the presence of oxygen [reviewed in [62]]. The total area occupied by anoxic microsites should be much smaller compared with the rest of the bulk sediment sampled that experienced oxygen during the first 7 days. Thus, the dominant OTUs detected at day 7 are unlikely to be those living in anoxic microsites. But, it is possible that some strictly anaerobic populations survived at low abundance under oxygen exposure during the first 7 days within anoxic microsites.

It is also likely that after the 28 days of incubation the label turned over, and a second generation of microbial cells became labeled not only from 18O water but also from 18O-labeled organic compounds that formed earlier. In this case, some of the labeling detected at the 28-day timepoint could have been due to assimilation of 18O-labeled organic compounds. However, the concentration of any 18O transferred to organic matter would be diluted within the much larger (unlabeled) dissolved organic matter pool. This would then be further diluted several orders of magnitude by the 18O label in the water that was present at a molar concentration. Thus, the potential assimilation of 18O-labeled organic compounds is unlikely to affect the conclusion that the degree of 18O labeling is a measure of assimilation of 18O from water, and thus activity, in growing populations.

Phylogenetic grouping of redox-specific activities

The oxygenated concentrations during the first week and development of anoxic and sulfidic conditions at 28 days allowed us classify 18O-labeled OTUs detected between these two timepoints into four categories of growth (Fig. 3) [1]; micro-aerophilic growth [2], facultative anaerobic growth [3], anaerobic growth, and [4] aerotolerant anaerobic growth (see Results for category definitions). The non-growing microorganisms were also considered as an additional category [5]: dormant or slow-growers. We recognize that because oxygen was present in the headspace throughout the incubation, strict anaerobic growth cannot be unequivocally assigned to the OTUs. But, given the increase in their growth later in the incubation after the onset of euxinic conditions at the sediment-water interface were established, we interpret this as an indicator of anaerobic growth. We also recognize that our designation of facultative versus aerotolerant anaerobic growth is arbitrary, but is used here to differentiate between those OTUs labeled at day 28 present also at day 7 that were, or were not, 18O labeled at day 7.

OTUs exhibiting micro-aerophilic growth had a relatively weak phylogenetic distribution across the phylogeny (λ = 0.38) (Fig. 3). In contrast, 18O-labeled OTUs exhibiting anaerobic growth corresponded to non-random phylogenetic groupings (λ = 0.9 and 0.87, respectively), implying that traits conferring anaerobic growth were conserved in the sampled communities. 18O-labeled organisms detected at both timepoints (facultative anaerobes) did not display a significant phylogenetic pattern (λ = 0.16), suggesting that facultative anaerobic growth was not a conserved trait in our sampled communities.

Populations exhibiting micro-aerophilic growth

A total of 47 OTUs affiliated with orders Flavobacteriales and Chitinophagales within the phylum Bacteroidetes were the most active growing bacteria after 7 days, maintaining metabolic activity under suboxic conditions. Similarly, seven OTUs affiliated with order Anaerolineae (Chloroflexi) were 18O labeled at day 7 (0.08 ± 0.008 EAF). Furthermore, 26 OTUs affiliated with Gammaproteobacteria (mainly orders Cellvibrionales and Thiotrichales) were labeled at day 7, and demonstrated relatively high EAF values (0.102 ± 0.031) indicating micro-aerophilic growth (Fig. 3).

18O-labeled OTUs affiliated with candidate class Gracilibacteria were observed solely at day 7, consistent with a micro-aerophilic growth (Table 1). Gracilibacteria have a cytochrome/quinol oxidase [18], most specifically cytochrome bd, which is implicated in ameliorating oxidative stress effects [reviewed in [63]]. Cytochrome bd could thus help to explain the higher growth of Gracilibacteria under low oxygen conditions.

Populations exhibiting anaerobic growth

There were 256 OTUs that exhibited a pattern of 18O labeling consistent with aerotolerant anaerobic growth (Fig. 3). They are apparently capable of surviving in the presence of oxygen during the beginning of the incubation, but their growth was maximized under anoxic conditions. The majority of these OTUs were affiliated with the Gammaproteobacteria (69 OTUs; 0.18 ± 0.06 EAF) and Deltaproteobacteria (19 OTUs; 0.29 ± 0.02 EAF). The highest 18O-labeled OTU was affiliated with genus Zhongshania that contains the facultative anaerobic heterotrophic species Zhongshania aliphaticivorans SM-2T [64]. In addition, one of the 18O-labeled Gammaproteobacteria OTUs was affiliated with the JTB255/Woesiaceae clade, which have been identified as the most dominant dark carbon-fixing microbes with a capacity to oxidize reduced sulfur compounds in anoxic and suboxic coastal sediments [65]. The genus Sandaracinus in Deltaproteobacteria known to degrade complex polysaccharides [66], also exhibited aerotolerant anaerobic growth. At day 28 after the onset of euxinic conditions, 10 OTUs affiliated with the Latescibacteria were detected that grew anaerobically (Table 1). This is consistent with their proposed fermentative mode of metabolism in anoxic sediment and water columns [18, 67].

OTUs affiliated with the known SRB genera [68] Desulforhopalus, Desulfosarcina, Desulfobulbus, Desulfopila, Desulfobacter, Desulfotignum, and Desulfatitalea were well-represented comprising 9% (42 OTUs) and 8% (26 OTUs) of the total sequences at days 7 and 28, respectively. However, only one SRB OTU was 18O labeled, which occurred at day 7 and was affiliated with the SRB genus Desulfobulbus (0.10 EAF). The relatively small number of growing Deltaproteobacteria SRB is low compared with the more numerous populations of sulfate-reducing Acidobacteria that have higher 18O EAF values (Figs. 3, 4). This indicates that these anaerobic Acidobacteria SRB were growing faster compared with the Deltaproteobacteria SRB.

Acidobacteria belong to several newly discovered groups of SRB, showing that dissimilatory sulfur metabolism is more widespread than previously thought [29, 69]. For example, Acidobacteria with a dissimilatory sulfur metabolism have been recently identified in acidic peatland [69, 70] and a DNA-SIP study [71] showed activity of dsrAB-containing organisms derived from Acidobacteria [70]. Our study shows that in addition to acidic peatland, dsrB carrying Acidobacteria grow in anoxic estuarine sediments with a relatively fast rate. Their increased growth rate and activity compared to Deltaproteobacteria SRB shown here implies that they should have a large impact on dissimilatory sulfur cycling under anoxic conditions.

Our finding that the majority of deltaproteobacterial sulfate reducers had minimal growth is in line with the previous reports that their mean in situ doubling times are on the order of months to years [72]. Metabolic activity (e.g., rRNA synthesis) in natural microbial communities is typically followed by cell division [16], but metabolic activity of non-growing organisms can also influence biogeochemical cycles [72]. For example, the increase in activity of SRB Candidatus Desulfosporosinus infrequens can occur independent of cell growth-associated processes [73]. Thus, the low degree of 18O labeling in abundant deltaproteobacteria SRB seen in our study (Fig. 3) may relate to a different ecophysiological strategy (e.g., slow growth) compared with the faster growing Acidobacteria SRB.

Within the Planctomycetes, members of the uncultured OM190 clade were abundant (44% total Planctomycetes) and only became 18O labeled at day 28 (Fig. 3), indicating anaerobic growth. Representatives of the OM190 clade closely related to anaerobic ammonia oxidizing bacteria (anammox) were also detected in hypoxic estuarine surface sediments in the East China Sea [74]. The anaerobic growth of OM190 clade organisms seen here is consistent with a potential anammox metabolism.

The 18O-labeled Actinobacteria were dominated by OTUs most closely related to the enigmatic Sva0996 actinobacterial clade first described from marine sediments [75], the “Candidatus Actinomarinales” [76], and Rhodococcus. While aerobic growth of Rhodococcus is well known as it relates to hydrocarbon degradation [77], our results showing facultative anaerobic growth indicate that Rhodococcus has potential to degrade hydrocarbons also under anoxic conditions in sediments. A facultative anaerobic lifestyle also explains why Rhodococcus are often found in deep subseafloor anoxic marine sediments [78].

The 18O-labeled Verrucomicrobia were dominated by OTUs affiliated with the uncultured DEV007 clade, first reported from the Elbe River in Germany (unpublished data, accession number: AJ401107). 18O-labeled OTUs affiliated with DEV007 were labeled at 7 days and became increasingly labeled after 28 days and the establishment of anoxic conditions (Fig. 2). This is consistent with biogeographic surveys that have detected this group in anoxic estuarine sediments (unpublished data, accession number: JN672646), OMZs (unpublished data, accession number: MG875625), and marine sediments [79]. However, the DEV007 clade is also found in oxic seawater attached to particles [35]. Our 18O labeling results showing a facultative anaerobic growth of the DEV007 clade explains this wide biogeographic range.

An OTU affiliated with the WPS-2 (Writtenberg Polluted Soil) clade was the only MDM group exhibiting facultative anaerobic growth (Fig. 2). The WPS clade was first described in a study of polychlorinated biphenyl-polluted soil in Germany [80] and was since detected in a wide range of oxic and anoxic environments [81]. The facultatively anaerobic growth shown here could potentially explain the ability of the WPS-2 group to survive in a large number of habitats with widely varying redox states.

Although total microbial abundance decreased over the incubation by nearly an order of magnitude, 16S rRNA genes became increasingly enriched in 18O (Fig. 1c). This raised the possibility that within populations, a high number of cells were dying while a smaller number were growing. In order to investigate this possibility further, we applied a model [51] that uses the 18O-qSIP data to calculate the number of 18O-labeled and -unlabeled 16S rRNA genes per OTU and the changes in their ratio over time to estimate rates of gross reproduction, mortality, and net production of individual OTUs (Fig. 5).

After seven days, only three OTUs exhibited net production (defined as net production 90% CI not overlapping zero), which were affiliated with the Sva0081 sediment group (Desulfobacteraceae), Desulfobacterium catecholicum, and SB-5 family of Bacteroidetes. Our analysis suggests that the Sva0081 marine benthic group (MBG), a putative group of SRB that are an important sink of acetate [82] and H2 [83] in coastal marine sediments, was one of the fastest growing populations under micro-aerophilic conditions at day seven (Fig. 5). The increased micro-aerophilic growth of Sva0081-MBG individuals is consistent with their higher acetate uptake rates in oxic–suboxic surface layer sediment slurries [82]. As Sva0081-MBG are H2 oxidizers [83], it is possible that the fast growing Sva0081-MBG individuals at the oxic–anoxic interface were also utilizing fermentation-derived H2 diffusing upwards from the deeper anoxic sediments. Desulfobacterium catecholicum, which exhibited net production at day seven (Fig. 5), is a mesophilic SRB that can also perform dissimilatory reduction of nitrate to ammonium [84]. This would explain its higher growth rate under the micro-aerophilic conditions at day seven during which nitrate was possibly still present. The SB-5 group of Bacteroidetes belongs to the Cytophagales and was previously detected in a benzene degrading, sulfate-reducing consortium [85]. The high net production under micro-aerophilic conditions indicates this group was growing faster compared with most other OTUs (Fig. 5).

After 28 days, most of the OTUs were 18O labeled (Fig 2). However, the rate of mortality per genus was higher at day 28 than day 7, indicating that establishment of euxinic conditions caused the majority of cells per genus to die faster than they grew (Fig. 5). However, despite the higher the net mortality rates, many exhibited relatively low positive gross reproduction rates. This indicates that a smaller proportion of individuals within each OTU were actively growing, while the majority of individuals were dying. This small fraction of active cells per population explains how nearly all OTUs became 18O labeled by day 28, despite of the increase in total microbial mortality.

Larger mortality compared with production at the end of the 28-day incubation could be partly explained by a decrease in electron donors (e.g., organic matter) and acceptors (e.g., nitrate, Fe(III), and sulfate). Indeed, under normal conditions organic matter and dissolved ions from seawater (e.g., sulfate) would be slowly but continuously entering the sediment to fuel new microbial growth. In our incubations, however, nutrients were constantly being depleted without replacement. Some of the inactive populations could also be explained by an increase in dormant cells due to unfavorable environmental conditions [9].

Conclusions

Our findings provide the first experimental evidence demonstrating the redox conditions promoting growth in several groups of uncultured “microbial dark matter”, validating hypotheses put forth by earlier metagenomics studies. The data help to explain previously observed biogeographic patterns for many uncultivated groups of bacteria that tend to correlate with anoxic or low oxygen conditions in aquatic habitats. This information could be helpful to guide future cultivation efforts for groups of ubiquitous, yet uncultured, bacterial taxa.

Fig S1