Potential oxygen consumption and community composition of sediment bacteria in a seasonally hypoxic enclosed bay

Study site description and sampling procedure

Omura Bay is in Nagasaki Prefecture, western Kyushu, Japan. It comprises steep cliffs at the margins and a large, flat seafloor at a water depth of 15–20 m (Kato et al., 2003). The bay is extremely enclosed, is connected to the open sea (East China Sea) by just two narrow channels. Cold and dense oxygenated water flows into the bottom of the bay in winter through these channels (Nogami et al., 2000; Takahashi et al., 2009). In early summer, warmer and less dense water intrudes into the middle layer of the bay, leaving the bottom water below the depth of intrusion stagnant, consequently leading to hypoxia in the central region of the bay (Nogami et al., 2000; Takahashi et al., 2009). In autumn, usually by October, this hypoxic water mass disappears as the vertical convection becomes activated by cooling of the sea surface (Suzaki, Miyake & Nakata, 2013).

Oceanographic monitoring and sediment core sampling were conducted at a central site in Omura Bay (Station 21; 20 m water depth; 32°55.390′N, 129°51.350′E; see Mori et al., 2018b). During the summer months of 2012 and 2013 (Table 1), sediment samples were hand-collected with an acrylic pipe (length 31 cm, inner diameter 26-mm) whilst scuba diving. Surface areas with visible bioturbation were not included in the samples. Temperature, salinity, DO, and chlorophyll-a fluorescence in the water column at the sampling site were measured using a multi-parameter monitoring device (AAQ, JFE-Advantec Co, Kobe, Japan). All of the sediment core samples were stored at an in situ temperature and carefully transported to the laboratory within 3 h of collection, avoiding direct exposure to sunlight and other physical disturbances. For the DNA analysis, the top layers (0–5 mm) of three sediment cores were pooled and stored at −20 °C until DNA extraction was performed.

For a series of experiment to compare INT reduction and oxygen consumption (slurry experiments, see below), additional cores were collected using a gravity core sampler (Ashura, Rigosha Co. Ltd, Japan), equipped with three polycarbonate tubes (length 57 cm, inner diameter 82 mm), in August 2017 and 2020. Each sediment core was sliced at a depth 0–5 cm, the sliced sediments were pooled on board the training ship, Kakuyo-Maru, and then stored at −80 °C until further processing.

For the time-course experiments, another set of additional sediment cores were collected from the same sampling site of Omura Bay (St. 21) using an acrylic tube (length 31 cm, inner diameter 26 mm), whilst scuba diving on August 7 and September 8, 2014.

Comparison between the extent of a tetrazolium dye reduction and oxygen consumption of the sediment samples (slurry experiments)

The sediment samples collected in 2017 and 2020 were used to determine the relationship between INT reduction and oxygen consumption in sediment. The tetrazolium dye reduction assay was based on the reduction of INT by the dehydrogenase activity of sediment microorganisms and reduced compounds such as sulfide in the sediment (Wada et al., 2012). To reveal the quantitative relationship between INT reduction and oxygen consumption by the sediment samples, we performed a series of incubation experiments using sediment slurries. A freshly thawed slurry sample was prepared by suspending sediment (approximately 3 g) into 1 L of aerated and filter-sterilized (0.22 µm) artificial seawater (Tetra Marine Salt Pro, Tetra Japan, Tokyo, Japan). The sediment slurry was further divided into two, one was incubated at 20 °C for 20 h without aeration prior to the assay, and the other was directly used for the assay without such a pre-incubation. Afterward, the samples were stored at room temperature and aeration was provided for different periods of time, for up to 48 h, to adjust the amount of reduced compounds in the samples. The sediment slurries were aliquoted into three glass bottles (100 ml each) to which INT (final concentration 0.01%) was added. Each of these samples was gently mixed with a stirrer for 24 h at 26 °C under dark conditions. Following this incubation, the slurries were centrifuged at 1,610× g for 10 min, and the supernatants were filtered through a cellulose acetate membrane filter (pore size 0.22 µm, Advantec, A020A025A). Both the filter and the pellet (sediment) were stored at below −20 °C until analysis. The reduced form of INT (INT-formazan [INTF]) was extracted from the filters and sediments according to the method of Wada et al. (2012). The INTF concentration in the extract was calculated by applying a standard curve with five different known concentrations (ranging from 0.5 to 50 µM) of pure INTF (Sigma-Aldrich, USA), dissolved in isopropanol for the filter and methanol for the sediment. We used the following equation to calculate the whole INT reduction rate: (1)${\displaystyle \text{INTF}\mu \mathrm{mol}{\mathrm{g}}^{-1}{\mathrm{day}}^{-1}=\frac{A+B}{D_w}\times \frac{1}{T}}$

where A and B are the amounts of INTF extracted from the filter and the sediment, respectively; D_w is the dry weight of the sediment used for the extraction of INTF; and T is the time of incubation.

Actual oxygen consumption by the sediment slurries was measured directly using an optode sensor (FireSting O2 fiber-optic oxygen meter; Pyro Sciences, Aachen, Germany) in parallel with the whole INT reduction measurement. Oxygen concentrations in the glass bottles were recorded every 30 s. The sediment slurry was incubated in duplicate under the same conditions as described for the assay of the sediment cores (26 °C, under dark conditions, 24 h). Ultra-pure MQ water (10 mL) was added to each of the slurry samples instead of the aqueous INT solution. Oxygen concentrations in the bottles were recorded for 24 h after the initiation of measurement. To estimate the oxygen consumption rate, a generalized additive model (GAM) was fitted to the changes in DO during the incubation period (Eq. 2). (2)${\displaystyle \mathrm{y}=f\left(x\right)+b}$

Here, y is the DO, x is time, f(x) is the smoothing function, and b is the intercept.

Instantaneous consumption rates were determined by taking the first-order derivative of the fitted GAM. These rates were integrated to obtain the total consumption rate during the experiment. The GAM was fitted using the mgcv package (Wood, 2003) in R (version 3.5.0), assuming a Gaussian distribution, an identity link function, and a thin plate spline as the smoothing function. The upper limit for the degrees of freedom of f(x) was 10.

For the sediment fixed with formaldehyde, we also measured the rate of chemical oxygen consumption (COC) and the non-biological (chemical) reduction of INT. Five samples of non-pre-incubated or pre-incubated sediment slurries (20 °C for 20 h without aeration) were prepared for the assay. Afterward, the samples were fixed by adding 10 mL formalin (formaldehyde: 3.7% w/v final concentration). After 15 min, INT (final concentration 0.01% w/v) and ultra-pure MQ water were added to the fixed slurries. They were incubated at 26 °C under dark conditions for 24 h. The COC of the slurry samples was determined in duplicate using the optode, as described above. INTF was extracted from the slurries in triplicate, and chemical INT reduction rate was determined in the same way as described for whole INT reduction. In an additional experiment, sediment slurries (non-pre-incubation sample) were supplied with sterile air for up to 24 h before formalin fixation to adjust the concentration of reduced compounds in the sediment. A generalized linear model (GLMs) assuming a gamma distribution with a log-link function was applied to analyze the entire set of data (Eq. 3). (3)${\displaystyle \begin{array}{c}\hfill y\sim Gamma\left(\mu ,\sigma \right)\hfill \\ \hfill log\left(\mu \right)=b_0+b_{1}x\hfill \end{array}}$

The model coefficients are b ₀ and b ₁, µis the expected value on the log-link scale, x is the INTF, σ is the shape parameter, and y are the observations (dissolved oxygen consumption rate).

Potential sediment oxygen consumption rate in omura bay estimated by an INT reduction assay

Potential oxygen consumption by the sediment samples was estimated by applying a tetrazolium dye reduction assay method to the intact sediment cores, as the overlying water was artificially replenished with oxygen (Wada et al., 2012). The whole INT reduction rate was estimated using the method of Wada et al. (2012). Briefly, triplicate cores were used for the measurement of whole INT reduction rate. Triplicate or duplicate cores, collected in 2012 and 2013, were used for the measurement of chemical INT reduction at each sampling time. The overlying water was replaced with 40 mL of aerated and filter-sterilized (0.22 µm) artificial seawater and 5 mL of MQ water and then gently mixed by pipetting. For chemical INT reduction measurement, 5 mL formalin was added instead of MQ water. After 10 min, 5 mL of 0.1% INT solution (w/v) was added to the overlying water and again gently mixed by pipetting. Upon replacing the water, special care was taken to avoid disturbing the sediment surface. Core samples in the acrylic pipes from each sampling time were incubated in the laboratory at 26 °C under dark conditions for 24 h. Following incubation, the overlying water was siphoned into a sterile plastic tube and filtered through a cellulose acetate membrane filter (diameter 25 mm, pore size 0.22 µm, Advantec, A020A025A). The filters were stored below −20 °C until analysis. The sediment cores were vertically extruded and sliced into horizontal sections of 0–5 mm depth. The sediment slices were stored at less than −20 °C. INTF extraction was performed as described by Wada et al. (2012). The INTF concentration in the extract was calculated by applying a standard curve derived from five different concentrations (ranging from 0.5 to 50 µM) of pure INTF (Sigma-Aldrich), dissolved in isopropanol and methanol for the water and sediment samples, respectively. The concentrations of INTF (µmol g⁻¹ day⁻¹) was calculated by combining the INTF values in the overlying water and those of the sediment as described in the above section. In order to simplify the calculation of INT reduction rate per unit mass of sediment (µmol g⁻¹ day⁻¹), we assumed the upper most sediment of 0–5 mm depth was responsible for the whole INT reduction in the overlying water of the core samples. Besides, the published data of INT reduction in 2011 (Mori et al., 2015) was integrated into the following statistical analysis. For the 2012 samples, the solvent extracted sediment residue was pooled and stored in each sampling month (May through September) at −20 °C until DNA extraction.

Time-course experiments

In order to validate the incubation time for the INT-reduction assay of Wada et al. (2012), triplicate core samples were collected in 2014 in the same way as in 2012 and 2013, incubated in the laboratory at 26 °C under dark conditions and stopped after different incubation times (0, 2, 4, 6, and 24 h). Following incubation, INTF extraction and the calculation were made as described above.

Principal Component Analysis (PCA) of the environmental parameters and correlation analysis of the pWOC, pBOC, and pCOC activities

Previously published data for the in situ environmental parameters (temperature, salinity, DO, and chlorophyll-a fluorescence in the bottom water and the quantities of sediment organic carbon, and nitrogen) (Mori et al., 2015; Mori et al., 2018b) were compiled for the PCA. The PCA allowed us to reduce the multicollinearity inherent in the environmental variables, and the first three principal components (PCs) were used as covariates for the ensuing GLMs to examine correlations between with pWOC, pCOC, and pBOC. (Eq. 4). (4)${\displaystyle \begin{array}{c}\hfill y\sim Gamma\left(\mu ,\sigma \right)\hfill \\ \hfill log\left(\mu \right)=b_0+b_1\mathrm{PC}1+b_2\mathrm{PC}2+b_3\mathrm{PC}3\hfill \end{array}}$

The model coefficients are b ₀, b ₁, b ₂, and b_3,µis the expected value on the log-link scale, σ is the shape parameter, and y are the observations (pWOC, pBOC, and pCOC).

PCA and the subsequent GLM analyses were performed using the prcomp and glm functions in R (version 3.5.0). In the visualization steps, the “factoextra” (Kassambara & Mundt, 2020), “broom” (Robinson, Hayes & Couch, 2020), and “ggplot2” (Wickham, 2016) packages were used.

Illumina MiSeq 16S rRNA gene amplicon sequencing

Microbial genomic DNA in the sediment samples was extracted using an ISOIL DNA extraction kit (Nippon Gene, Tokyo, Japan), according to the manufacturer’s instructions. A total of ten genomic samples obtained from the surface sediment layer (0–5 mm depth) at the five sampling times in 2012, collected before and after incubation for the INT reduction measurement, were amplified for Illumina MiSeq 16S rRNA gene amplicon sequencing. Primers specific for the V3–V4 regions of the 16S rRNA gene, 341F and 805R (Herlemann et al., 2011), were used. The PCR was conducted in a 25-µL reaction mixture using a PCR buffer, 5 µl forward primer (1 µM), 5 µl reverse primer (1 µM), 12.5 µl 2X Kapa HiFi HotStart ReadyMix (Kapa Biosystems, Boston, MA, USA), and 10 ng of the extracted DNA. Then, 27 cycles of PCR with an initial step of 95 °C for 3 min were run, followed by 95 °C for 30 s, 55 °C for 30 s, 72 °C for 30 s, and finally, 72 °C for 5 min. The size of the PCR products was verified by agarose gel electrophoresis, and the PCR products were further purified with Agencourt AMPure XP (Beckman Coulter, Indianapolis, USA), according to the manufacturer’s instructions. A second PCR reaction was performed on the purified PCR products (2.5 µl) to index each of the samples. Two indexing primers (Illumina Nextera XT 72 indexing primers, Illumina, California, USA) were used per sample. The second PCR was carried out in a 25-µL reaction mixture using a PCR buffer, 5 µl forward index primer (1 µM), 5 µl reverse index primer (1 µM), 12.5 µl 2X Kapa HiFi HotStart ReadyMix (Kapa Biosystems), and 2.5 µl of the first PCR products. The second PCR products were then also purified (Agencourt AMPure XP; Beckman Coulter). Following purification, all samples were sent to the Bioengineering Lab. Co., Ltd. for massively parallel 300-bp paired-end sequencing, which was performed on an Illumina MiSeq analyzer. Processing and quality control of the reads was performed using R. For raw sequence data, the adapter sequences were trimmed using the “quasR” package (Gaidatzis et al., 2015). Following quality control processes, inclusive of chimera removal, the removeBimeraDenovo, filterAndTrim and dada function in the “DADA2” package version 1.6 (Callahan et al., 2016) were used, with default settings. Finally, the amplicon sequence variants (ASVs) (Callahan, McMurdie & Holmes, 2017) were classified from kingdom to genus level using the SILVA 132 SSU Ref NR 99 database (Quast et al., 2013), resulting in the construction of an ASV table with read counts for the ASVs in all samples. In the visualization step, the “phyloseq” (McMurdie & Holmes, 2013) and “ggplot2” packages were used to create stacked bar and plots. To delineate a hierarchical clustering plot, the ASV table was standardized based on coverage-based rarefaction (Chao & Jost, 2012). Thereafter, the distance matrix was calculated based on the Bray–Curtis dissimilarity, and clusters were calculated using Ward’s method. The “vegan” package (Oksanen et al., 2019) and hclust function were used to calculate the dissimilarity and clusters.

All associated raw data can be found at the DDBJ Sequence Read Archive under accession numbers DRA008748 (DRX177081–DRX177090).

Comparison between the extent of tetrazolium dye reduction and oxygen consumption of the sediment samples (slurry experiments)

The whole INT reduction rates in the sediment slurries ranged from 7.4 to 28.5 µmol INTF g⁻¹ day⁻¹, while actual oxygen consumption rates ranged from 44.3 to 390.5 µmol O₂ g⁻¹ day⁻¹ (n = 16). The chemical INT reduction ranged from 0.6 to 25.6 µmol INTF g⁻¹ day⁻¹, while oxygen consumption from the same samples ranged from 38.8 to 147.1 µmol O₂ g⁻¹ day⁻¹ (n = 25) (Fig. 1, Table S1). We found a significant relationship between INT reduction rate and oxygen consumption rate in Omura Bay sediment (Eq. 3 P <  0.0001, Fig. 1). If we estimate the values of the pWOC based on the GLM model, it ranged from 45.9 to 91.1 µmol O g⁻¹ day⁻¹, which can be converted to areal unit using density of sediment and the area of sediment core (45.3 to 56.1 mmol O₂ m⁻² day⁻¹). These pWOC values were comparable with the WOCs reported for the basins of other coastal areas, such as the East China Sea (3.6–17.6 mmol O₂ m⁻² day⁻¹, Song et al. (2016)), the northern Gulf of Mexico (0.82–56.4 mmol O₂ m⁻² day⁻¹, Rowe, Kaegi & Morse (2002)), and the Seto Inland Sea (3.6–17.6 mmol O₂ m⁻² day⁻¹, Nakamura (2003)). Therefore, the INT reduction method should work reasonably well to estimate SOC in coastal environments.

The extent of INT reduction has been proven to be an reliable proxy to estimate the oxygen consumption rate of planktonic microorganisms (Martínez-García et al., 2009; Villegas-Mendoza, Cajal-Medrano & Maske, 2019; García-Martín et al., 2019), despite the inherent toxicity of INT (Martínez-García et al., 2009; Villegas-Mendoza, Cajal-Medrano & Maske, 2015). Previous studies also showed a positive correlation between INT reduction and oxygen consumption for estuary sediments (Trevors, 1984) and terrestrial soils (Trevors, Mayfield & Inniss, 1982). However, none of them have addressed whole INT reduction and chemical INT reduction separately or determined their correlation coefficients with SOC. Thus, to the best of our knowledge, the present results are the first to demonstrate the relationship between INT reduction and oxygen consumption in marine sediment with a statistically significant GLM model that would allow us to estimate potential SOC based on the INT reduction rates.

Time-course experiments

We found that the quantity of INTF increased consistently for up to 24 h without appreciable decline in selected sediment core samples (Fig. S1). This was in stark contrast to what had been reported to the planktonic samples. Martínez-García et al. (2009) and others (Villegas-Mendoza, Cajal-Medrano & Maske, 2015; García-Martín et al., 2019) demonstrated that prolonged incubation times can cause an underestimation of the INT reduction rate due to breakage of microbial cells that accumulated INTF crystals inside their cellular membrane. The seemingly discrepancy between sediment and planktonic samples may be attributable to better extraction efficiency of INTF from sediment samples. Regardless of cell breakage, majority of INTF crystals formed in sediment microbial cells were likely to be retained on the sediment surface and were directly extracted from the sediment slices with methanol. In contrast, INTF crystals formed in planktonic samples could be lost if the size of INTF grain was smaller than that of pore size of the filters. Furthermore, the time for INT in overlying water to reach and exhibit toxicity on microbial cells in the sediment may be delayed due to a slow diffusion of INT across the boundary layer that would develop over the sediment surface under the static condition of the assay.

Temporal variability of potential SOC

Based on the GLM model obtained with the slurry experiments, the potential SOC in the intact core samples ranged from 45.3 to 56.1 mmol O₂ m⁻² day⁻¹. However, as most of the whole INT reduction in the core samples turned out to be far below the minimum value in the slurry experiments (Fig. 1, Table S1), we decided not to describe and discuss the results any further based on the conversion the INT reduction rates into SOC, but to use whole, chemical and biological INT reduction as proxies for potential WOC (pWOC_INT), COC (pCOC_INT), and BOC (pBOC_INT), respectively.

Figure 2 showed the three components of SOC at the central part of Omura Bay during early summer to late autumn in 2012 and 2013 (Figs. 2A–2C, Table 1) together with those in 2011 that had been already published (Mori et al., 2015) in order to make the downstream statistical analysis more robust. Both pWOC_INT and pCOC_INT tended to be greater in August and September than in other sampling months. According to the criteria used by Wright, Konwar & Hallam (2012), we defined hypoxic conditions as <90 µMO₂, and divided the DO level into four categories: Oxic (>90 µMO₂), dysoxic (90–20 µMO₂), suboxic (20–1 µMO₂), and anoxic (<1 µMO₂). The potential activity of SOC was highest under suboxic or anoxic conditions in the bottom water. As the highest pBOC_INT in 2012 and 2013 (2.82 µmol INTF g⁻¹ day⁻¹ in 2012 and 4.16 µmol INTF g⁻¹ day⁻¹ in 2013) coincided with the timing of the highest bacterial count in each year (1.7  × 10¹¹ cells/g sediment in 2012 and 3.2  × 10¹⁰ cells/g sediment in 2013, Table 1), it was likely that the increases in facultative anaerobes or aerobes contributed to the surge of the pBOC_INT. In five out of 20 cases, pCOC_INT exceeded pWOC_INT, and pBOC_INT showed negative values (Table 1, Figs. 2A–2C). This may be partly explained by the effect of formalin addition on SOC as described above, which lead to an enhancement of COC in the sediment.

In order to identify the environmental driver(s) for the variation of potential SOC, we first conducted PCA for all the environmental variables examined (DO concentration, temperature, and salinity in the bottom water within 0.5 m above the seafloor; and TOC and TON contents in the sediment), during 2011 through 2013 (Table 2, Fig. 3A) and tested the correlation between the obtained principal components (PCs) and the potential SOC (Table 3). It was indicated that the first three PCs explained 93% of the variability of the variables (Table 2). The results of PC1 suggest that TOC as well as TON in the sediment increased in the opposite directions from DO and temperature, reflecting a tendency of organic matter preservation in the sediment under hypoxia. On the other hand, the results of PC2 suggest that temperature increased in the opposite directions from salinity and DO, reflecting a temporal variation of those variables in hypoxic water mass overlying the sampling location. GLM analysis showed that only PC2 had a significant effect on pWOC_INT and pCOC_INT (Table 3, Figs. 3B, 3D). These results strongly suggest that temperature as well as DO in the bottom water acts as environmental drivers that control the potential SOC.

Bacterial community composition and its response to re-oxygenation

A total of 2,546 unique ASVs were identified in the surface sediment from the center of Omura Bay in 2012. The dominant taxa in the bacterial communities were affiliated with the Deltaproteobacteria, Gammaproteobacteria, and Bacteroidia (Fig. 4). Within the Deltaproteobacteria class, the family groups of Desulfobacteraceae (major genera Sva0081 sediment group and Desulfobacter) and Desulfobulbaceae (major genera Desulfopila and MSBL7), which belong to the dominant order Desulfobacterales, were predominant (Fig. S2). In contrast, family groups in Woeseiaceae (major genus Woeseia), Halieaceae (major genus Halioglobus), and Thiotrichaceae (which were not assigned to genus level) were predominant in the Gammaproteobacterial class (Fig. S3). Among the bacterial groups, Woeseiaceae was recently reported to be a cosmopolitan and abundant core member of the bacterial communities in marine sediments globally (Dyksma et al., 2016; Mußmann et al., 2017). These bacterial groups are frequently reported to be present in marine sediment surfaces (Kawahara et al., 2009; Zinger et al., 2011; Dyksma et al., 2016; Jochum et al., 2017). Therefore, the bacterial community composition in the surface sediment of Omura Bay demonstrated the ubiquity of the organisms (e.g., Desulfobacteraceae and Woeseiaceae) found in marine sediments.

According to our previous study (Mori et al., 2018a; Mori et al., 2018b), temporal shifts in the bacterial community in the surface sediment of Omura Bay were characterized by decreases in the relative abundance of the ASVs affiliated with Gammaproteobacteria toward less oxic conditions; conversely, those affiliated with Deltaproteobacteria exhibited the opposite results. The present results were essentially consistent with this trend, although there was a peak in the relative abundance of Deltaproteobacteria under the oxic condition in September 2012 (Fig.4). A strong typhoon, Sanba, had passed near Omura Bay on the day before we collected the sediment cores in September 2012. Due to the strong wind (>12 m/s on average) of the typhoon, water column of entire bay must have been mixed rigorously and oxygenated the bottom water. This mixing event is likely to be largely responsible for the contradiction between bottom water oxygen levels and the peak in Deltaproteobacteria. In contrast to a large fluctuation of DO in the water column under the typhoon, relative abundance of Deltaproteobacteria in the sediment seemed less affected.

Following the incubation (24 h) of sediment core samples, more than 90% of the ASVs (comprising 90.2%–96.6% of the total reads) did not substantially change their relative abundance (<0.1% of the abundance). Consistent with this notion, the ASVs were tightly clustered by sampling months (Fig. 5). However, it was found that several ASVs exhibited large increases or decreases during the incubation period. The majority of the increased ASVs belonged to the Arcobacteraceae, Fusobacteriaceae and Halieaceae family groups, while those of the decreased ASVs were affiliated with the Flavobacteriaceae, Saprospiraceae and Spirochaetaceae family groups (Fig. 6). The relative abundance of the Arcobacteraceae family group major genus Arcobacter increased from <1% to 12.9% within 24 h in the sediment cores that were retrieved under anoxic bottom water conditions (August 24, 2012) (Table S2). It was interesting to note that the increase in the relative abundance of Arcobacter ASVs coincided with the highest pBOC (Fig. 2). The Arcobacter comprise a group of physiologically diverse bacteria that exhibit anaerobic dissimilatory manganese reduction (Vandieken et al., 2012) or aerobic sulfide oxidation (Wirsen et al., 2002; Sievert et al., 2007). Given that these bacterial groups have been found in marine sediments under suboxic conditions and were reportedly engaged in the sulfur cycles (Wirsen et al., 2002; Sievert et al., 2007), the surge in pBOC_INT in August 2012 could be ascribed to the aerobic oxidation of reduced sulfur compounds by Arcobacter. Consistent with this interpretation, a recent report on the metatranscriptomics of anoxic marine sediment revealed Arcobacter was highly responsive to oxygenation event and upregulated expression level of the genes encoding for respiratory proteins (Broman et al., 2017). Therefore, it is very likely that the anoxic sediment of Omura Bay harbored an Arcobacter population capable of rapidly growing in response to the replenishment of DO. Besides Arcobacter, it is reasonable to assume the other bacterial groups whose relative abundance increased notably during the incubation period, such as Fusobacteriaceae and Halieaceae (Fig. 6), also contributed to the pBOC_INT.

On the other hand, it was difficult to associate the dominant bacterial groups with potential respiratory activities. Even though a majority of the predominant components of the bacterial community, such as the Woeseiaceae group, is likely to play a role in sediment metabolism, the present results provided no further insights to support this inference. One possibility for the contribution of dominant groups to the SOC is that highly abundant anaerobic sulfate reducing bacteria might have affected the seasonal dynamics of pBOC via their aerobic respiration. The sulfate reducing bacteria inhabiting the sediment surface generally possess abilities to cope with high oxygen stress, and some of those are capable of aerobic respiration (Cypionka, 2000; Brioukhanov, Pieulle & Dolla, 2010; Zhou et al., 2011), which possibly contributes to pBOC_INT during the oxygenation experiment_. Indeed, the maximum pBOC_INT observed in August 2012 was the timing when cells abundance was peaked and sulfate reducing Deltaproteobacteria, including the groups capable of aerobic respiration such as Desulfobacter (Dannenberg et al., 1992), was highly abundant in the surface sediment (Fig. 2B, Table 1, Fig. S2). It is not clear whether electrogenic sulfur oxidation by cable bacteria plays roles in Omura Bay, as we did not detect the cable bacteria genus (i.e., Candidatus Electrothrix). It will be worth investigating its presence using other methods, such as fluorescent in situ hybridization, in any future studies.