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 [17]. 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 [18, 19]. 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 [18, 19]. In autumn, usually by October, this hypoxic water mass disappears as the vertical convection becomes activated by cooling of the sea surface [20].
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. [15]). 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 of 2017 and 2020. Each sediment core was sliced at a depth of 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 sets 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 [7]. 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 20h 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. [7]. 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), dissolved in isopropanol for the filter and methanol for the sediment. We used the following equation to calculate the whole INT reduction rate (WIR):

where A and B are the amounts of INTF extracted from the filter and the sediment, respectively; Dw 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 WIR measurement. Oxygen concentrations in the 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. 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 [21] in R (version 3.5.0), assuming a Gaussian distribution, an identity link function, and a thin plate spline.
For the sediment fixed with formaldehyde, we also measured the rate of chemical oxygen consumption (COC) and the non-biological (chemical) reduction of INT (CIR). 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 stored 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 CIR was determined in the same way as described for WIR. 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. 2).

The model coefficients are b0 and b1, μ 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 [7]. The WIR rate was estimated using the method of Wada et al. [7]. Briefly, triplicate cores were used for the measurement of WIR. Triplicate or duplicate cores, collected in 2012 and 2013, were used for the measurement of CIR 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 CIR 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. [7]. 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. INTF (µmol g-1 day-1) was calculated by combining the INTF values in the overlying water and those of the sediment as described in the above section (See “Comparison between the extent of a tetrazolium dye reduction and oxygen consumption of the sediment samples (Slurry experiments)”). In order to simplify the calculation of INT reduction rate per unit mass of sediment (µmol g-1 day-1), 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 [8] 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. [7], 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) [8, 15] 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. 3).

The model coefficients are b0, b1, b2, and b3, μ 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” [22], “broom” [23], and “ggplot2” [24] 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 [25], 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 [26]. Following quality control processes, inclusive of chimera removal, the removeBimeraDenovo, filterAndTrim and dada function in the “DADA2” package version 1.6 [27] were used, with default settings. Finally, the amplicon sequence variants (ASVs) [28] were classified from kingdom to genus level using the SILVA 132 SSU Ref NR 99 database [29], resulting in the construction of an ASV table with read counts for the ASVs in all samples. In the visualization step, the “phyloseq” [30] 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 [31]. Thereafter, the distance matrix was calculated based on the Bray–Curtis dissimilarity, and clusters were calculated using Ward’s method. The “vegan” package [32] 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).