Anaerobic methane oxidation in a coastal oxygen minimum zone: spatial and temporal dynamics.

Coastal waters are a major source of marine methane to the atmosphere. Particularly high concentrations of this potent greenhouse gas are found in anoxic waters, but it remains unclear if and to what extent anaerobic methanotrophs mitigate the methane flux. Here we investigate the long-term dynamics in methanotrophic activity and the methanotroph community in the coastal oxygen minimum zone (OMZ) of Golfo Dulce, Costa Rica, combining biogeochemical analyses, experimental incubations and 16S rRNA gene sequencing over 3 consecutive years. Our results demonstrate a stable redox zonation across the years with high concentrations of methane (up to 1.7 μmol L-1 ) in anoxic bottom waters. However, we also measured high activities of anaerobic methane oxidation in the OMZ core (rate constant, k, averaging 30 yr-1 in 2018 and 8 yr-1 in 2019-2020). The OPU3 and Deep Sea-1 clades of the Methylococcales were implicated as conveyors of the activity, peaking in relative abundance 5-25 m below the oxic-anoxic interface and in the deep anoxic water respectively. Although their genetic capacity for anaerobic methane oxidation remains unexplored, their sustained high relative abundance indicates an adaptation of these clades to the anoxic, methane-rich OMZ environment, allowing them to play major roles in mitigating methane fluxes.

Introduction

Microbial processes are the main source of the potent greenhouse gas methane, CH4, to the atmosphere, where it accounts for about one‐quarter of the increase in radiative forcing since the preindustrial era (IPCC, 2013; Etminan et al., 2016). Marine methane emissions originate mainly from coastal and shelf environments, yet despite recent modelling and meta‐data analysis efforts, estimates of methane emissions from such systems remain uncertain (4–27 Tg CH4 yr−1; Rosentreter et al., 2021; Weber et al., 2019), due to sparse measurements, high spatial and temporal variability, and changing conditions as a result of human influence (Rosentreter et al. 2021; Weber et al., 2019). Notable areas of marine methane accumulation (tens of nanomolar to tens of micromolar) include anoxic basins and fjords as well as coastal upwelling systems, particularly those characterized by oxygen‐depleted subsurface waters known as oxygen minimum zones (OMZs; Capelle et al., 2019; Kessler et al., 2006; Naqvi et al., 2010; Sansone et al., 2001; Thamdrup et al., 2019). Here, benthic archaeal methanogenesis is suggested to be the major methane source (Sansone et al., 2001; Chronopoulou et al., 2017), but molecular data have additionally revealed the potential for pelagic methanogenesis (Padilla et al., 2016). It is, however, unclear to what extent methane fluxes from anoxic waters are mitigated by microbial oxidation, and the microorganisms and metabolic pathways potentially involved remain elusive. Nonetheless, to understand these processes is imperative given the predicted expansion of marine oxygen‐depleted waters (Breitburg et al., 2018).

Microbial methane oxidation can be performed both aerobically and anaerobically. The latter uses alternative oxidants such as sulfate, iron, manganese, nitrate and nitrite, with nitrate and nitrite being the most energetically favourable (see Guerrero‐Cruz et al., 2021 and references therein). Overlapping distributions of oxidized nitrogen compounds and methane in oxygen‐depleted waters, for example in the eastern tropical North Pacific (ETNP) OMZ and Golfo Dulce, Costa Rica, suggest a niche for microbes coupling anaerobic methane oxidation to reductive nitrogen transformations (Padilla et al., 2016; Padilla et al., 2017; Thamdrup et al., 2019). Measuring rates of anaerobic metabolisms in these waters requires extensive protocols to minimize oxygen contamination that could inhibit anaerobic pathways (De Brabandere et al., 2012). Thamdrup et al. (2019) utilized procedures developed to minimize oxygen contamination and were able to quantify rates of anaerobic methane oxidation in the ETNP OMZ between 0.014 and 0.12 nmol L−1 d−1. Highest rates were located between the nitrite and methane maxima and rates were consistently inhibited with 0.5 μmol L−1 of oxygen (Thamdrup et al., 2019). Whether anaerobic methane oxidation plays a comparable role in other OMZ systems is still unknown as are the temporal dynamics of the process.

Molecular analysis performed in anoxic marine waters indicates that diverse microbes may be responsible for linking the methane and nitrogen cycles. In the OMZs of the Golfo Dulce and ETNP, close relatives of ‘Candidatus Methylomirabilis oxyfera’ of the NC10 clade have been shown to be present and transcriptionally active (Padilla et al., 2016; Thamdrup et al., 2019). Members of the NC10 are hypothesized to oxidize methane aerobically by producing intracellular oxygen through the dismutation of nitric oxide (NO) from nitrite reduction, which allows them to thrive in anoxic waters despite the obligate need for oxygen (Ettwig et al., 2010). A coupling of anaerobic methane oxidation to NO in the ETNP OMZ was further indicated by rate measurements being inhibited by the presence of the NO scavenger PTIO (Thamdrup et al., 2019). Further links to the nitrogen cycle in the ETNP OMZ have been indicated by the recovery of transcripts affiliated with ‘Candidatus Methanoperedens nitroreducens’ of the ANME‐2d clade, known to couple anaerobic methane oxidation to nitrate reduction (Haroon et al., 2013; Thamdrup et al., 2019).

Another group of methanotrophs detected in anoxic marine waters are members of the Gammaproteobacterial order Methylococcales, which are traditionally considered as aerobic methanotrophs (Hayashi et al., 2007; Padilla et al., 2017; Thamdrup et al., 2019). Diverse marine clades of Methylococcales have been identified based on the phylogenetic marker gene particulate methane monooxygenase (pmo) and are referred to as operational pmo units (OPUs, Tavormina et al., 2008). Among these, the clade OPU3 has been shown to be abundant in the ETNP OMZ, peaking in abundance at <4 μM of oxygen (Tavormina et al., 2013), and in the anoxic core of the Golfo Dulce OMZ (Padilla et al., 2017), the latter hinting at a potential role in anoxic waters. A metagenomic study performed in Golfo Dulce showed that an abundant member of the OPU3 clade could perform partial denitrification to NO (Padilla et al., 2017) and the analysis of a related isolate strain likewise demonstrated the ability to link methane oxidation with denitrification to nitrous oxide (N2O) when oxygen concentrations were <50 nmol L−1 (Kits et al., 2015). However, neither of the studies demonstrated methane oxidation completely independent of oxygen, since oxygen still appears to be required for the initial methane oxidation step. It, therefore, remains to be determined how members of OPU3 either obtain oxygen for methane oxidation or otherwise metabolize within anoxic waters.

Here, we present data from a 3‐year (2018–2020) study of methane oxidation in the coastal OMZ of Golfo Dulce. Our aim was to explore the temporal and spatial dynamics of anaerobic methane oxidation in a biogeochemical and microbial context. To do so, each year we quantified rates of methane oxidation in the anoxic core of the OMZ, measured relevant biogeochemical parameters, and investigated the methanotroph community through 16S rRNA gene amplicon sequencing. We found that (i) the biogeochemical zonation in the OMZ, which had high concentrations (up to 1.7 μmol L−1) of methane, was stable between years. (ii) The OMZ supported high rates of anaerobic methane oxidation, (iii) most likely performed by members of Methylococcales, although the exact pathway remains unclear.

Results

The stations sampled in Golfo Dulce in 2018, 2019 and 2020 were adopted from Thamdrup et al. (1996) (Fig. 1A). In 2018, all sampling was performed at Sta. 1 at the head of the bay, whereas in 2019 and 2020, a six‐station transect of chemical profiles was sampled along the central axis of the bay (Sta. 1 through Sta. 3a) in addition to high‐resolution profiles of water column chemistry, biomolecular samples and methane oxidation rates at Sta. 1 and Sta. 2 (Fig. 1A; Table S1). We refer to the chemical data collected during sampling of the 2019 and 2020 transects as transect data, while the depth profiles of chemical‐, molecular‐ and rate data additionally obtained from Sta. 1 and Sta. 2 are referred to as high‐resolution profiles (HR1 and HR2 respectively). The depth profiles sampled at Sta. 1 in 2018 are referred to as HR1a and HR1b, collected on January 29th and February 2nd, respectively.

Fig. 1

Study area and distribution of chemical parameters along the central axis of Golfo Dulce, Costa Rica.

A. Map of stations sampled in the Golfo Dulce with stations sampled for high‐resolution rate profiles in red. Concentration of oxygen from the standard range oxygen sensor in 2019 (B) and trace range oxygen sensor (<2 µmol L−1) in 2020 (C) (no standard range oxygen data were available in 2020). Concentrations of nitrate (D and E), nitrite (F and G), and methane (H and I) measured in January 2019 (left panels) and January 2020 (right panels) respectively. Black line indicates the oxic–anoxic interface.

Biogeochemical profiles along Golfo Dulce

The transects along the bay revealed a consistent vertical redox zonation, which was largely conserved between the 2 years sampled (2019 and 2020; Fig. 1). Oxygen concentrations decreased from a shallow mixed surface layer to below the detection limit of the trace range oxygen sensor (<5 nmol L−1) at an oxic–anoxic interface at mid‐water depth. In 2019, the interface deepened from 87 m (specific density, σ , 25.52 kg m−3) at the head of the bay to 108 m (25.61 kg m−3) at Sta. 3, before dipping to 145 m (25.68 kg m−3) at the outermost station (3a). In 2020, the interface remained at 100–110 m (25.46–25.50 kg m−3) along the transect (Fig. 1B and C). Due to challenges in calibration of the standard range oxygen sensor only data from the trace range oxygen sensor (<2 μmol L−1) were available during the transect in 2020. Nitrate was depleted in the surface waters in both years and then increased to a maximum of 13–19 μmol L−1 at 50–70 m depth. Nitrate concentrations decreased gradually from the lower oxycline to the sediment–water interface, reaching 3–5 μmol L−1 at 180–190 m depth in 2019, and ≤1 μmol L−1 in 2020 (Fig. 1D and E). A primary nitrite maximum was apparent at 30 m in 2019 (Fig. 1F). This feature was not evident in 2020 but may have fallen between sampling depths. As is characteristic to anoxic OMZs, the concentration of nitrite increased below the oxic–anoxic interface (Fig. 1F and G). In both years, the highest nitrite concentrations (0.6–1 μmol L−1) were measured at the head of the bay (Sta. 1–2) between 130 and 150 m. In 2019, nitrite at these stations remained high (>0.5 μmol L−1) from its peak to 180–190 m depth; in 2020, nitrite decreased with depth in the OMZ (Sta. 1‐2a). At the mouth of the bay, maximum nitrite concentrations were lower (0.4–0.6 μmol L−1) and occurred at deeper depths (180 m). The anoxic waters were further characterized by high concentrations of methane, which in both years increased approximately linearly from the oxic–anoxic interface (~0.005 μmol L−1) to the sediment–water interface. Methane concentrations peaked at 0.7–1.1 μmol L−1 at 180–190 m across stations, with the lowest values at the mouth of the bay (Fig. 1H and I). No sulfide was detected in any of the years. Ammonium concentrations below the oxic–anoxic interface were consistently less than 100 nmol L−1, and the profiles showed no distinct features with depth (data not shown).

High‐resolution profiles of water column chemistry and methane oxidation

High‐resolution depth profiles of water column chemistry were sampled at Sta. 1 (HR1) and Sta. 2 (HR2) in 2018 (Sta. 1 only), 2019 and 2020 and reflected an apparent steady state with little variation in chemical zonation between years (Fig. 2). Similar to the profiles recorded during sampling of the transect in 2019 and 2020 (Fig. 1), nitrate decreased with depth in the OMZ, while nitrite and methane accumulated. Notably, the highest concentration of methane in the OMZ (1.7 μmol L−1) was measured in 2018 at Sta. 1 (Fig. 2D). Despite the stability of the OMZ chemical zonation, the position of the oxic–anoxic interface fluctuated over time. In 2018, the interface at Sta. 1 deepened from 114 m (25.61 kg m−3) to 127 m (25.64 kg m−3) over 4 days (HR1a on January 29th and HR1b on February 2nd) and similar shifts were also observed in 2019, as the interface at Sta. 1 deepened by 16 m between transect sampling (87 m, 25.52 kg m−3) and high‐resolution sampling (103 m, 25.55 kg m−3) 3 days later.

Fig. 2

Depth profiles of oxygen (A, E) (inset shows oxygen concentrations from the trace range oxygen sensor), nitrate (B, F), nitrite (C, G) and methane (D, H) concentrations at Sta. 1 (top, A‐D) and Sta. 2 (bottom, E‐H). HR1a and HR1b indicate repeated sampling of Sta. 1 in 2018, on January 29 and February 2 respectively (Sta. 2 was not sampled in 2018).

The 3 years of high‐resolution profiles of methane oxidation activity measured with minimal oxygen contamination at Sta. 1 and Sta. 2 showed that rate constants of methane oxidation, k, determined as the fraction of methane tracer turned over per unit time, were generally low at oxic depths but increased sharply below the oxic–anoxic interface (Fig. 3B, F, J). Above the oxic–anoxic interface, k was ≤1 yr−1 in 2019 and 2020, but ranged from 3 to 15 yr−1 in 2018, with the highest values measured near the interface (100–115 m). In general, k was high in 2018 compared to the two following years. The profiles from 2019 and 2020, when both high‐resolution stations were sampled, showed that k at Sta. 2 was generally higher than at Sta. 1. In 2018, k peaked at 50 yr−1 at 130 m depth during HR1a and at 57 yr−1 at 180 m 4 days later during HR1b. In the following 2 years, three out of the four profiles (HR1 and HR2 2019 and HR1 2020) showed high values of k in a primary peak corresponding with the increase in nitrite concentrations (k, 5 to 8 yr−1 at 110 to 145 m depth) and again in a secondary maximum occurring at the deepest depth (k, 3 to 8 yr−1 at 180–190 m depth), where methane concentration was highest (Fig. 3F and J). In contrast at Sta. 2 in 2020 (HR2), the highest activity occurred at 150 m (23 yr−1) (Fig. 3J).

Fig. 3

Depth profiles of methane concentrations (A, E, I), rate constants of methane oxidation (k; B, F, J – note variability in x‐axis), methane oxidation rates (C, G, K – note variability in x‐axis), and relative abundances of Methylococcales (D, H, L). Data from 2018 (top, A–D), 2019 (middle, E–H) and 2020 (bottom, I–L). Solid lines with circles show Sta. 1 (HR1) and dashed lines with diamonds show Sta. 2 (HR2). HR1a and HR1b (A–C) indicate repeated sampling of Sta. 1 in 2018 and only HR1b is shown in panel D (no Methylococcales abundance data are available for HR1a). Error bars represent the standard error and open symbols represent rate constants/rates that are not significantly different from zero. Horizontal lines mark the position of the oxic–anoxic interface, where solid lines on panels E–L show the position during HR1 at Sta. 1 and dashed lines show the position during HR2 at Sta. 2, while on panels A–D, solid lines show the position for HR1a and dashed lines show the position for HR1b. Rate constants, k (B, F, J) and rates of methane oxidation (C, G, K) were measured in incubations with 3H‐CH4 tracer and <100 nmol L−1 oxygen. Methyloccocales abundances were obtained with Illumina (D, L) and PacBio sequencing (H) respectively. Full‐depth profiles (0–200 m) can be seen in Fig. S2.

The kinetics experiment performed in 2019 showed that methane oxidation rates, calculated by multiplying k by the methane concentration, in the OMZ increased approximately linearly with added methane up to 1.5 μmol L−1 (Fig. 4). While no signs of saturation were present in incubations from the OMZ, incubations from the oxycline (95 m) fitted a Michaelis–Menten model (K  = 0.63 ± 0.1 μmol L−1 and V max = 16.2 ± 1.4 nmol L−1 d−1) only indicating saturation at the highest methane concentration applied (~1.6 μmol L−1), well above methane concentrations observed at that depth (~0.02 μmol L−1). Our observations thus confirmed the assumption of first‐order reaction kinetics for methane oxidation for the range of concentrations found at the different depths (Reeburgh et al., 1991; Valentine et al., 2001; Thamdrup et al., 2019) and justified using k and in situ methane concentrations to calculate methane oxidation rates.

Fig. 4

Kinetics of methane oxidation in incubations with varying methane additions from the oxycline, A (95 m), and the OMZ core, B (130 m) and C (180 m). Vertical error bars represent the standard error of the slope obtained from time series incubations while horizontal error bars represent the standard deviation of added methane concentrations. Michaelis–Menten kinetics were fitted to the data in panel A (95 m; K  = 0.63 ± 0.1 μmol L−1; V max = 16.2 ± 1.4 nmol L−1 d−1). Regression lines are shown for panels B and C (130 m: r 2 = 0.99, p << 0.001; 180 m: r 2 = 0.99, p << 0.001). See Fig. S1 for a magnified view of incubations with <0.6 μmol L−1 methane.

Rates of methane oxidation from depths above the oxic–anoxic interface were always ≤0.2 nmol L−1 d−1 (Fig. 3C, G, K). Rates in the OMZ core were highest in 2018, increasing from 0.2 nmol L−1 d−1 at the oxic–anoxic interface to 243 nmol L−1 d−1 at 180 m depth. In comparison, rates in 2019 and 2020 increased to a maximum of 31 nmol L−1 d−1 (170 m depth). In all 3 years, rates of methane oxidation correlated with the in situ concentration of methane (r 2: 0.6–0.9; p < 0.01).

Methanotroph and methanogen community composition

The microbial community composition was analysed by sequencing of 16S rRNA gene amplicons, using two different sequencing methods across years: Illumina with universal primers in 2018 and 2020, and PacBio with bacteria‐specific primers in 2019. The relative abundance of Archaea varied between 10% and 35% of the microbial community across all samples from the OMZ in 2018 and 2020. As detailed below, the two methods yielded very similar results with respect to methanotroph phylogeny and relative abundance.

Sequencing of amplicons revealed 31 methanotroph‐related ASVs across the 3 years in Golfo Dulce. These ASVs were classified as members of the Gammaproteobacterial order Methylococcales and clustered into three environmental clades: OPU1 (five ASVs), OPU3 (23 ASVs) (Hayashi et al., 2007; Tavormina et al., 2013) and ‘Deep Sea‐1’ (three ASVs), with the latter including the genus Methyloprofundus sedimenti (Tavormina et al., 2015) (Fig. 5). We did not identify any sequences affiliated with anaerobic methane oxidizers of the nitrite‐reducing NC10 clade (Ettwig et al., 2010), the nitrate‐reducing ANME‐2d clade (Raghoebarsing et al., 2006; Haroon et al., 2013), or any Alphaproteobacterial methanotrophs. Vertical profiles of Methylococcales relative abundance were similar across years. Methylococcales abundance at HR1 was negligible at oxic depths, increased to 0.5%–7.8% of amplicons in the oxycline (Fig. 3D, H, L), and peaked at 9%–11% in the anoxic OMZ core (~5–25 m below the oxic–anoxic interface). In 2019 and 2020, abundances decreased with depth below their peak, declining to ~1.5% at 180 m just above the sediment. In 2018, a second Methylococcales peak (9%) was present at 180 m depth.

Fig. 5

16S rRNA gene phylogeny and depth distribution (coloured circles) of methanotrophs in Golfo Dulce in 2018 (black), 2019 (blue) and 2020 (orange). Consensus phylogeny based on full (1461 bp, GD19) and partial (253 bp, GD18 and GD20) 16S rRNA sequences, included sequences have relative abundances >0.05% of total amplicons. Circles show relative ASV abundance (%) with depth. Methanotroph OTUs from Golfo Dulce in 2015 (labelled with star, Padilla et al., 2017) are included for comparison. Scale bar represents nucleotide substitutions. Dotted line indicates the position of the oxic–anoxic interface. Bootstrap values >50 (grey) and >70 (black) are shown.

The three identified clades of Methylococcales differed in their depth distribution (Fig. 5). The distribution of OPU3 mirrored that of the total Methylococcales community, with OPU3 sequences comprising 60%–100% of the Methylococcales community at all depths. The OPU3 community was dominated by a single OPU3 ASV in all datasets with a nucleotide identity of >99% between years (GD18‐ASV 4, GD19‐ASV 4 and GD20‐ASV 7 in 2018, 2019 and 2020 respectively) (Fig. 5). ASVs of the OPU1 clade (GD18‐ASV 651, GD19‐ASV 1694, GD20‐ASV 637 and GD20‐ASV 779) occurred in lower relative abundances (≤35% of Methylococcales sequences) and were only observed in samples from the oxycline. In contrast, ASVs of the Deep Sea‐1 clade (GD18‐ASV 65, GD19‐ASV 309 and GD20‐ASV 85) appeared almost exclusively at depths below the oxic–anoxic interface, reaching up to 40% of the Methylococcales community in samples from 125 to 180 m (Fig. 5).

In Illumina‐based datasets from 2018 and 2020 containing archaea, we additionally searched for methanogen‐related 16S rRNA gene sequences and detected sequences belonging to the methanogenic Euryarchaeal class ‘Candidatus Methanofastidiosa’, which based on metagenomic and transcriptomic analyses have been indicated to perform methanogenesis through reduction of methylated compounds (Nobu et al., 2016; Zhang et al., 2020; Wang et al., 2021) (Fig. S3). The vertical distribution of ‘Ca. Methanofastidiosa’ differed between 2018 and 2020. In 2018, ‘Ca. Methanofastidiosa’ constituted <0.07% of total sequences when detected, with no obvious relationship with depth. In 2020, the group increased with depth to reach 1.2% in the OMZ core (130–170 m depth; Fig. S4).

Discussion

The 3 years of data on biogeochemical zonation, methane oxidation rates and microbial community composition demonstrate that the anoxic waters of Golfo Dulce serve as a stable niche for Gammaproteobacterial methanotrophs and support high rates of anaerobic methane oxidation. Below we discuss methane dynamics within the bay and how the biogeochemical setting may shape the structure and function of the methanotrophic community, and thereby ultimately its role in mitigation of methane emissions in the coastal OMZ.

Long‐term stability of the redox zonation in Golfo Dulce

The three expeditions to the Golfo Dulce OMZ (2018–2020) revealed a consistent vertical redox zonation across years (Fig. 2), as well as along the bay (Fig. 1). The accumulation of nitrite at anoxic depths indicated active nitrate reduction. We detected nitrate at all anoxic depths in 2018 and 2019, while it was depleted (LOD 0.5 μM) at 180 m at the innermost four stations in 2020 (Fig. 1E, Sta. 1‐2a). The lack of detectable sulfide suggested that nitrate and nitrite were the preferred electron acceptors over sulfate at the depths examined. These conditions are similar to those found in open‐ocean OMZs (Ulloa et al., 2012).

The redox zonation observed in 2018–2020 was broadly similar to patterns shown in previous studies in Golfo Dulce since the first surveys along the bay in March 1969 (Richards et al. (1971). This prior work further includes a mapping of the transect in January 1994 (Thamdrup et al., 1996; Ferdelman et al., 2006) as well as profiles from selected stations in November 2001 and January 2015 (Dalsgaard et al., 2003; Padilla et al., 2016). One notable difference is that the previous studies have typically detected low levels of sulfide (≤7 μmol L−1) in deeper waters at the head of the bay with the most recent report from 2015 finding sulfide at and below 150 m at Sta. 1 (up to 6.6 μmol L−1; Padilla et al., 2016). However, in 1994, sulfide was only detected in one profile from Sta. 1 and then disappeared a few days later (Ferdelman et al., 2006). Thus, the differences between prior results and those from our study may potentially reflect minor short‐term variability rather than a long‐term trend. Overall, the similarity of data obtained over more than 50 years indicates a remarkable stability of the Golfo Dulce OMZ.

Distribution and potential sources of methane

Methane consistently accumulated with depth in the OMZ with highest concentrations near the sediment–water interface (Figs 1H, I and 2D, H). The concentrations of methane in Golfo Dulce (≤1.7 μmol L−1) are high compared to those of open ocean OMZs (typically ≤0.1 μmol L−1, Chronopoulou et al., 2017; Jayakumar et al., 2001; Sansone et al., 2001) but similar to concentrations observed in the seasonally anoxic fjord of Saanich Inlet (1.3 μmol L−1; Capelle et al., 2019) and anoxic basins such as the Black Sea and Cariaco (up to 17 μmol L−1; Kessler et al., 2006 and references therein). The single previous study of methane concentrations in Golfo Dulce from January 2015 also reported an increase with depth, but only to 80 nmol L−1 at 180 m (Padilla et al., 2016). At that time, bottom waters contained 6.6 μmol L−1 sulfide, in contrast with our study where sulfide was not detected. The higher concentrations of methane in 2018–2020 can therefore not be explained by more reduced conditions in the OMZ. We see no obvious explanation for this difference and further measurements are needed to determine whether such variability is a recurring phenomenon.

The accumulation of methane towards the bottom waters could suggest that the sediments are the main methane source to the OMZ. No studies to date have analysed the biogeochemistry of methane in the sediments of Golfo Dulce. In general, upper layers of marine sediments contain low concentrations of methane, as the supply of sulfate from seawater typically favours sulfate reducers over methanogens in shallower sediment layers, and strongly attenuates the upward flux of methane from deeper layers through anaerobic oxidation in the sulfate methane transition zone (Iversen and Jorgensen, 1985; Knittel and Boetius, 2009). During the survey of Golfo Dulce in 1994 (Thamdrup et al., 1996) high concentrations of sulfate (>24 mM) were measured in the upper 30 cm of the sediments, with no sign of decline with depth (B. Thamdrup unpubl. res.). A large benthic methane flux, therefore, seems unexpected considering that, in continental shelf sediments, methane concentrations in the sulfate zone are usually <0.01 mmol L−1 and sediments represent a negligible methane source to the water column (Iversen and Jorgensen, 1985; Niewöhner et al., 1998; Beulig et al., 2019).

Other methane sources that could contribute to the high concentrations observed in the Golfo Dulce OMZ are seeps/vents and methanogenesis in the water column. Seeps and vents are common on the Costa Rican Pacific continental margin (Mau et al., 2006; Sahling et al., 2008) and methane seeps have been observed in Golfo Dulce. The currently documented seeps are located at water depths ≤10 m (Berrangé, 1987; Wild et al., 2015), yet undiscovered methane seeps could be located deeper in the bay. Another possibility is that methane is produced in situ in the water column. This hypothesis is supported by the previous recovery of transcriptionally active methanogenic archaea in the Golfo Dulce OMZ (Padilla et al., 2016) and our current observation of the putative methyl‐reducing methanogens ‘Ca. Methanofastidiosa’ (Nobu et al., 2016; Zhang et al., 2020) in our 16S rRNA gene dataset (Figs S3 and S4). While the descriptions of ‘Ca. Methanofastidiosa’ are from an anaerobic digester and mangrove sediment, we are not aware of published reports of their occurrence in OMZ waters. However, screening of the NCBI database identified a ‘Ca. Methanofastidiosa’ sequence recovered from the Eastern Tropical South Pacific OMZ (Belmar et al., 2011). Phylogenetic analysis subsequently revealed a clade within the ‘Ca. Methanofastidiosa’ that exclusively contained sequences from pelagic habitats, including those from OMZs (Fig. S3). In Golfo Dulce, methylotrophic methanogens would likely have to compete for substrates against nitrate/nitrite respiring methylotrophs in the water column and might therefore be confined to reduced micro‐environments such as inside sinking particles (Beck et al., 2014), including faecal pellets from anoxic gut environments of zooplankton (Stief et al., 2017; Wäge et al., 2020). There was a notable difference in the relative abundance of ‘Ca. Methanofastidiosa’ between 2018 and 2020 (Fig. S4). Since the biogeochemical conditions remained similar and the sequencing procedure and analysis were the same between these 2 years, we see no obvious cause for the abundance difference. Nonetheless, the presence of ‘Ca. Methanofastidiosa’ in the water column, and especially its high abundance in 2020, suggests this taxon may play a role in methane production in the OMZ. We suggest that both benthic and pelagic sources should be evaluated further through biogeochemical experiments and analysis of molecular capabilities in order to understand the origin of methane in the Golfo Dulce OMZ.

A role for anaerobic methane oxidation

The consistent observation of a steep increase in rate constants and rates below the oxic–anoxic interface indicates that the anoxic, nitrate‐, nitrite‐ and methane‐rich environment in the OMZ of Golfo Dulce is a niche for anaerobic methane oxidation. In our incubations, which were carried out with minimal oxygen contamination (<100 nmol L−1 oxygen, see Experimental procedures), the activity from above the interface was far lower (k < 1 yr−1) than in the OMZ (1–57 yr−1), with the exception of Sta. 1 in 2018 (HR1a), where rate constants in the oxycline ranged from 9 to 15 yr−1 (Fig. 2A). These higher values might be due to greater mixing around the oxic–anoxic interface, as suggested by the 13‐m deepening in interface depth between sampling dates 4 days apart, which could have transported methanotrophs from deeper depths to the oxycline. Unfortunately, 16S rRNA gene sequence data were not available for HR1a in 2018, leaving us unable to explore features of the methanotroph community that might further explain the rate variation.

Rate constants are assumed to scale with the population size of microbes involved in methane oxidation. The observed depth distribution of rate constants, therefore, suggests a community of methanotrophs specialized for OMZ conditions. The activity observed near the oxic–anoxic interface might still be influenced by occasional oxygen supply from mixing or photosynthesis (e.g. Tiano et al., 2014; Garcia‐Robledo et al., 2017). Based on the 16S rRNA gene sequence data, Cyanobacteria were present below the oxic–anoxic interface in Golfo Dulce at relative abundances between 0.5% and 1% (Fig. S5d). However, light (>0.1 μmol photons m−2 s−1) was never observed to penetrate below the interface (Fig. S5C and G). The highest relative abundance of Cyanobacteria occurred in the shallowest samples, suggesting that the Cyanobacteria observed in the OMZ were sinking and likely inactive, as no chlorophyll was detected below the oxic–anoxic interface either (Fig. S5B andF). Therefore, photosynthesis is unlikely to serve as an internal oxygen source for methane oxidizers in this OMZ (incubations were also carried out in the dark). While mixing is likely to supply oxygen occasionally to microbes in the upper OMZ, the likelihood decreases with the distance to the oxic–anoxic interface, and we never detected oxygen deeper than 110 m at our experimental stations. Thus, our finding of the highest rate constants measured at depths between 130 and 190 m suggests an anaerobic metabolism for methane oxidizers.

Contamination with low levels of oxygen appears inevitable when performing laboratory incubations with OMZ waters (Ganesh et al., 2015; Thamdrup et al., 2019). Continuous monitoring of oxygen in our incubations constrained oxygen to low nanomolar amounts (<100 nmol L−1; see Experimental procedures), which could potentially sustain microaerophilic methane oxidation as observed in lakes (Blees et al., 2014; Oswald et al., 2015) and coastal waters (Steinle et al., 2017). Still, the methane oxidation rates of up to 243 nmol L−1 d−1 observed in our incubations could in several cases not be sustained by the maximum oxygen concentrations in incubations, given the expected oxygen to methane ratio of 2:1 for aerobic methane oxidation (Naguib, 1976), or even with a decreased ratio if parts of the metabolism are carried out anaerobically (Dam et al., 2013; Kalyuzhnaya et al., 2013). Moreover, the methanotrophs would face competition for any available oxygen from other members of the microbial community, such as oxygen‐respiring heterotrophs. Previously measured oxygen consumption rates from the Golfo Dulce OMZ reached as high as 1000 nmol L−1 d−1 (Garcia‐Robledo et al., 2016), which in our incubations would largely deplete an oxygen pool of <100 nmol L−1 within the incubation period. In the kinetic experiments rates were linear throughout the incubation period, thus we saw no indications of a shift from aerobic to anaerobic methane oxidation. Taken together these lines of evidence strongly indicate that methane oxidation activity observed in the OMZ core of the Golfo Dulce was not dependent on oxygen and was anaerobic, as also concluded for similar experiments in the ETNP OMZ (Thamdrup et al., 2019).

Rate constants of anaerobic methane oxidation from the OMZ core in 2018 (6–57 yr−1) were high compared to those measured the subsequent 2 years (1–23 yr−1, Fig. 3B, F, J). The difference might in part be driven by the steeper methane gradients observed in 2018 (Fig. 2D) that suggest a higher substrate flux (the density gradient in the OMZ core did not change between years, Fig. S6) that would be expected to sustain a more active community (Brune et al., 2000). Notably, incubations from the OMZ core in 2018 were carried out with lower methane concentrations than in situ, while in situ concentrations were re‐established in the incubations in 2019 and 2020 (see Experimental procedures). However, the methane concentration does not affect the calculation of k (see Experimental procedures) and based on the first‐order kinetics observed in the OMZ (Fig. 4; S1) and typically assumed for methane oxidation (e.g. Valentine et al., 2001), this difference should not affect the determination of rates. Alternatively, as the radiolabel constituted a larger fraction of the methane in the incubations from 2018 than in subsequent years, a faster turnover of the tracer relative to the unlabelled pool could potentially have led to higher rates in 2018. While we cannot fully exclude such a kinetic effect, we note that previous comparisons of methane oxidation rates obtained with 3H and 14C labelled methane respectively, report similar rates with the two tracers in systems with rates in the range observed in Golfo Dulce (Mau et al., 2013; Pack et al., 2015), arguing against such a bias. In addition to the inter‐annual variability, the activity exhibited some variability within years, despite relatively stable biogeochemical conditions (Fig. 3B, F, J). We see no clear explanation for this variability, which suggests that unresolved chemical, physical or biological dynamics play a role in modulating methane oxidation rates.

Overall, rate constants from the Golfo Dulce OMZ (1–57 yr−1) are similar to estimates from vent/seep influenced waters (≤73 yr−1; Chan et al., 2019, ≤25 yr−1; Heintz et al., 2012) and coastal systems (1–31 yr−1, Steinle et al., 2017), which are environments characterized by high methane turnover rates. Values were 1–2 orders of magnitude higher than those estimated for the open ocean ETNP OMZ (0.2–3.3 yr−1; Pack et al., 2015, ~1 yr−1; Thamdrup et al., 2019). This observation is similar to relative differences in rates observed between the coastal OMZ of Golfo Dulce and the oceanic OMZs for the anaerobic processes denitrification and anammox (Dalsgaard et al., 2003; Bulow et al., 2010; Dalsgaard et al., 2012) and likely reflects the overall higher areal productivity of coastal zones (Sigman and Hain, 2012) and the resulting higher input of substrates that sustain efficient elemental cycling relative to the open ocean (Kalvelage et al., 2013).

The high rate constants of anaerobic methane oxidation suggest this process could be a substantial methane sink in Golfo Dulce. An additional sink is mixing, which has been suggested to be driven primarily by the horizontal advective exchange of water that occurs due to the inflow of intermediate Pacific water (Richards et al., 1971; Ferdelman et al., 2006). Ferdelman et al. (2006) using a salt balance model estimated the replacement time of water in the deep basin of Golfo Dulce to be 1–2 months. Averaging measured rate constants over anoxic depths, we estimate a turnover time for methane in the OMZ core of ~12 days in 2018 and between 1 and 4 months in 2019–2020, indicating that anaerobic methane oxidation and mixing serve as near equal methane sinks in the OMZ.

Our kinetics experiment showed that anaerobic methane oxidation in Golfo Dulce followed first‐order reaction kinetics (Fig. 4), which is generally assumed to hold for methane oxidation in marine waters (Reeburgh et al., 1991; Valentine et al., 2001; Thamdrup et al., 2019). First‐order kinetics allowed us to use rate constants and in situ methane concentrations to calculate methane oxidation rates. The steep increase in methane concentration with depth in the OMZ resulted in a similar trend in anaerobic methane oxidation rates, which were typically highest near the sediment–water interface (Fig. 3C, G, K). The rates measured in the OMZ spanned a wide range (0.02–243 nmol L−1 d−1) but nonetheless fell in the upper half of an even broader set of aerobic methane oxidation rates estimated from various marine waters (10−5 to 103 nmol L−1 d−1; Mau et al., 2013), where the variability is largely attributed to variations in methane concentrations. Thus, anaerobic methane oxidation appears at least as efficient as the aerobic equivalent.

Methylococcales as the dominant methane oxidizers

Members of Methylococcales were the only known methane oxidizers recovered from the Golfo Dulce water column in all years sampled (Fig. 5). Relative abundances of this group peaked each year (6%–11% of amplicons) below the oxic–anoxic interface between 115 and 130 m (Fig. 3D, H, L). In 2018 we observed a second peak in relative abundance at 180 m (9%), which could potentially support the particularly high methane oxidation activity observed at that depth in 2018, although relative abundances cannot be translated directly to absolute abundances. In 2019 and 2020, the highest proportional representation of Methylococcales (115–130 m, 6%–9%) roughly corresponded to the primary peak in rate constants (110–130 m, 5–8 yr−1; Fig. 3). In contrast, the secondary peak in methane oxidation activity (3–8 yr−1) at 180–190 m depth was not reflected in Methylococcales abundance, suggesting that other taxa may have contributed to methane oxidation at that depth in 2019/2020. The high relative abundance of Methylococcales and their dominance of the aerobic methanotroph community in Golfo Dulce align with work from the same site in 2015 by Padilla et al. (2017). In contrast, in 2015 members of the nitrite‐dependent NC10 clade were also detected in Golfo Dulce through qPCR (up to 0.006% of total 16S rRNA sequences; Padilla et al., 2016) and may have escaped detection in this study due to low abundance. We did not identify any other groups of methanotrophs in our 16S rRNA gene datasets, but the possibility for overlooked players cannot be ruled out.

The Methylococcales community each year was dominated by a single ASV from the OPU3 clade (Fig. 5). The three sequences were highly similar between years and closely related to the OTU (GD_7) that dominated methanotroph sequences in Golfo Dulce in 2015 (Padilla et al., 2017). Although these members of OPU3 typically peaked in abundance between 5 and 25 m below the oxic–anoxic interface, they had a high proportional representation (50%–100% of methanotrophs) across all depths and on both sides of the interface. In contrast, the distribution of OPU1 and Deep Sea‐1 ASVs showed redox driven niche separation, where OPU1 was confined to samples from the oxycline, while members of Deep Sea‐1 only appeared in samples from anoxic depths (with one exception; Fig. 5). Similar redox driven niche separation of Gammaproteobacterial methanotrophs has been observed in stratified lakes (Mayr et al., 2020).

A niche separation between the two OPU groups is consistent with earlier observations of the two clades in the region. In the OMZ off the Costa Rican coast, OPU3 was observed to primarily inhabit the centre and edges of the OMZ, while OPU1 dominated in the oxycline below the OMZ (Tavormina et al., 2013). Previously in Golfo Dulce, OPU3 peaked in abundance inside the OMZ while OPU1 was only present in the oxycline (Padilla et al., 2017). Members of Deep Sea‐1 are commonly found in methane seep‐influenced surface sediments and as endosymbionts in seep‐associated fauna (Redmond et al., 2010; Raggi et al., 2013; Tavormina et al., 2015), but have also been observed in OMZ waters of the Eastern Tropical Pacific (Hayashi et al., 2007) and in the anoxic water column of the Black Sea (Glaubitz et al., 2010).

Our molecular and experimental evidence point to Methylococcales as the main, if not only, conveyors of anaerobic methane oxidation in the Golfo Dulce OMZ. This aligns with evidence from other anoxic aquatic settings, including transcriptional activity of OPU3 members in the anoxic core of the ETNP OMZ (Thamdrup et al., 2019) and repeated observations of Gammaproteobacterial methanotrophs in deep anoxic lakes (Blees et al., 2014; Oswald et al., 2016; Naqvi et al., 2018). The chemical zonation indicates nitrate and nitrite as the most favourable electron acceptors in the OMZ (Fig. 2). Indeed, Padilla et al. (2017) showed that the dominating OPU3 OTU in Golfo Dulce in 2015, GD_7, contained and transcribed genes for dissimilatory nitrate and nitrite reduction to NO and suggested that this group could link partial denitrification to methane oxidation (Padilla et al., 2017), an ability that the closely related OPU3 members identified in 2018–2020 likely share with GD_7.

A coupling between methane oxidation and denitrification (to N2O) has been demonstrated experimentally and genetically for the strain Methylomonas denitrificans of the Methylococcales (Kits et al., 2015). However, a means for members of Methylococcales, including the OPU3 clade and Methylomonas denitrificans, to circumvent the need for oxygen for the initial methane activation step by the oxygen‐requiring particulate methane monooxygenase (pmo) has not been identified. While Padilla et al. (2017) could not exclude the possibility that the OPU3 population was sustained to at least some degree via oxygen through a fluctuating oxic–anoxic interface or transient inputs, our finer spatial resolution, demonstration of long‐term stability, and rate measurements all point to an anaerobic lifestyle. The detection of the Deep Sea‐1 clade at anoxic depths in the Golfo Dulce (Fig. 5) suggests that this group may be even less likely to rely on oxygen, although the ability for anaerobic methane metabolism has not been investigated in this clade. It thus seems likely that these groups possess the ability for anaerobic methane oxidation, and further investigation of the genetic potential of these clades should have a high priority.

Taken together, our work indicates that the Golfo Dulce OMZ contains a stable redox zonation with high concentrations of methane compared to open‐ocean OMZs, and an active methane cycle that includes high rates of methane oxidation and, potentially, also pelagic methanogenesis. Based on the consistent steep increase in methane oxidation rate constants below the oxic–anoxic interface measured from incubations containing less than 100 nmol L−1 oxygen, our results strongly indicate that methane oxidation in the OMZ is anaerobic. Our data also suggest that methane oxidation is primarily conducted by specific members of the OPU3 and Deep Sea‐1 clades of the Methylococcales, with the former being present at particularly high relative abundances, though the exact pathway for anaerobic methane oxidation in the two clades remains to be described. Contributions from other unidentified taxa are also possible, including the NC10 clade reported in a previous investigation (Padilla et al., 2016), although this clade seems to be in too low abundance to explain the high rates seen in bottom waters, where the relative abundance of Methylococcales is comparatively low. Thus, from our 3 years of investigation, we hypothesize that members of the OPU3 and Deep Sea‐1 clades are adapted to the stable, anoxic, methane‐rich environment in the OMZ, and substantially reduce the methane flux from the OMZ.

Experimental procedures

Study site, sample collection and chemical analysis

Golfo Dulce is a coastal basin on the Pacific coast of Costa Rica with a water depth of ~200 m and anoxic bottom waters (Richards et al., 1971; Thamdrup et al., 1996). Sampling in the Golfo Dulce was performed in January–February of 3 consecutive years, 2018–2020. Stations were adopted from Thamdrup et al. (1996) and positions and sample types are listed in Table S1. All sampling in 2018 was performed at Sta. 1 at the head of the bay, whereas in 2019 and 2020, a six‐station transect of biogeochemical profiles was sampled along the central axis of the bay (Sta. 1 through Sta. 3a) in addition to high‐resolution profiles of water column chemistry, biomolecular samples and methane oxidation rates at Sta. 1 and Sta. 2. The chemical data collected during sampling of the 2019 and 2020 transects are referred to as transect data, while the high‐resolution chemical‐, molecular‐ and rate depth profiles obtained from Sta. 1 and Sta. 2 are referred to as HR1 and HR2 respectively. The depth profiles from Sta. 1 in 2018 are referred to as HR1a and HR1b, measured on January 29th and February 2nd respectively.

Oceanographic parameters were measured with a CTD (conductivity, temperature and depth‐measuring device), equipped with a PAR sensor, a fluorometer (Sea & Sun Technology), and two optical oxygen sensors, a standard range sensor (Borisov et al., 2011) and a trace range sensor, with a detection limit of ~5 nmol L−1 produced according to Larsen et al. (2016). Water from discrete depths was collected from hand‐deployed Niskin bottles to measure sulfide, methane, nitrate, nitrite, and ammonium concentrations, as well as for 3H‐CH4 incubations and molecular sampling. All water was sampled through gas‐tight viton tubing. Samples for sulfide were collected in 12 ml glass vials (Exetainers® Labco, UK) filled from the bottom with three volumes overflow and no headspace. Sulfide concentrations were determined by the methylene blue method (Cline, 1969), with samples either amended immediately with reagent and analysed within 4 h of collection (2020) or stored with 100 μl 50% wt./vol. zinc chloride until analysis (2018 and 2019). Water for methane quantification was sampled into 60 ml serum bottles with three volumes overflow. A 5 ml air headspace was inserted followed by 0.9 ml 6 N HCl for preservation and the bottle was crimp‐sealed with a butyl rubber septum. Methane was analysed via headspace analysis on a gas chromatograph (Thermo Trace 1300 equipped with an FID detector). In 2018, the samples were stored between 1 and 2 months before analysis, whereas in 2019 and 2020, methane was analysed on the day of sampling after vigorous shaking to equilibrate the water and gas. Methane concentrations were calculated using the temperature‐ and salinity‐dependent distribution coefficient (Wiesenburg and Guinasso Jr., 1979). Samples for nitrite and nitrate were filtered through 0.2 μm cellulose acetate filters and analysed spectrophotometrically (García‐Robledo et al., 2014). Nitrite was analysed within 5 h of collection and samples for nitrate were frozen until analysis by reduction to nitrite with acidic vanadium (III) (García‐Robledo et al., 2014). Samples for ammonium were collected in 12 ml glass exetainers, amended with orthophthaldehyde reagent and analysed fluorometrically (Holmes et al., 1999).

Methane oxidation rate incubations

Rates of methane oxidation were determined in incubations with 3H labelled methane (Valentine et al., 2001) following the procedure described in Steinsdóttir et al. (2022) to measure rates of methane oxidation in OMZ waters with minimal oxygen contamination. Briefly, water was sampled from the Niskin bottle through gas‐tight viton tubing into 250 ml glass bottles with ≥2 volumes overflow. Bottles were closed bubble‐free with deoxygenated butyl rubber stoppers (De Brabandere et al., 2012) and stored dark and cool until returned to the land‐based laboratory (< 5 h) where they were immediately placed in an incubator at in situ temperature (16°C) for re‐acclimation (~1 h). Bottles were removed from the incubator and a headspace was introduced before the water was purged with helium through a glass frit for 15 min (flow rate 0.4 L min−1). Water was dispensed through glass and viton tubing into 12 ml glass vials (Exetainers Labco, UK) that were filled with >2 volumes overflow and immediately closed with deoxygenated chlorobutyl rubber septa. A helium headspace (1 ml in 2018–2019 and 2 ml in 2020) was introduced into each vial, which subsequently underwent two cycles of 30 s of vigorous shaking followed by 30 s of flushing the headspace with helium (flow rate 0.5 L min−1) to drive out traces of oxygen. The final concentration of oxygen in incubations was always <100 nmol L−1 as measured from a subset of vials mounted with oxygen optodes (details in ‘Monitoring oxygen in incubations’ section). Samples were returned to the incubator for re‐acclimation to in situ temperature (~1 h) before being injected with 10 μl 3H‐CH4 tracer in N2 gas and, additionally in 2019 and 2020, with unlabelled methane as methane saturated water to restore in situ concentrations, via gas‐tight syringes. The unlabelled methane stock was prepared as methane saturated water in a gas‐tight glass vial. Additions of methane saturated water to incubations ranged between 5 and 50 μl depending on the target concentration. Further details on tracer preparation and methane additions are provided below (see ‘Verification of methane additions, scintillation counting and calculations’). Incubations were terminated as triplicates by injection of 100 μl of either 50% wt./vol. zinc chloride (2018 and 2019) or saturated mercuric chloride (2020), either immediately after addition of tracer or after 24 h of incubation. The switch to mercuric chloride was made as we observed lower counts in killed controls with mercuric chloride (1023 ± 214 dpm ml−1 compared to 3961 ± 379 dpm ml−1 in the killed controls with zinc chloride) and thereby obtained better sensitivity. Incubations were carried out in the dark at in situ temperature (16°C). The samples killed immediately were stored under the same conditions as live incubations, and hence served as killed controls to account for any potential abiotic transfer of 3H between the 3H‐CH4 tracer and water.

The kinetics of methane oxidation were investigated in a time course experiment (terminated at 0, 6, 12 and 24 h, n = 2) for three depths (95, 130 and 180 m) at Sta. 1 in 2019. The in situ concentrations of methane at these depths were 17, 216 and 855 nmol L−1 respectively. For each depth we tested six different concentrations of methane by injecting only tracer or tracer and additional methane‐saturated water. Methane concentrations in incubations covered a range of 18–1662 nmol L−1 (18 ± 3, 71 ± 7, 122 ± 10, 283 ± 10, 510 ± 65 and 1662 ± 87 nmol L−1) as determined in a subset of vials via gas chromatography.

Monitoring of oxygen in incubations

The concentration of dissolved oxygen in incubations was measured from a subset of the glass vials mounted with optical oxygen sensors for non‐invasive measurements. Vials with optodes were kept under a helium atmosphere for 24 h prior to incubation to remove any traces of oxygen dissolved in the silicon used to mount the optode. In 2018 and 2019, the measurement procedure followed that described in Thamdrup et al. (2019), where vials were mounted with trace range optode spots (TROXSP5, Pyroscience GmbH) and measurements were performed every 4–6 h using a Firesting fibre‐optic oxygen meter (Pyroscience GmbH; LOD ~ 50 nmol L−1). In 2020, vials were mounted with in‐house manufactured optodes with a LOD of ~0.5 nmol L−1 (Lehner et al., 2015) that were calibrated via multipoint calibration at in situ temperature (e.g. Larsen et al., 2016). Measurements were carried out continuously (every 2 min) during the 24‐h incubations.

Verification of methane additions, scintillation counting and calculations

The tracer mixture (3H‐CH4 in N2) was prepared within 3 weeks of use from 20 Ci mmol−1 single‐labelled 3H‐CH4 aliquots (American Radiolabeled Chemicals, St. Louis, MO, USA). Each aliquot was diluted with N2 to a final activity of ~0.05 mCi ml−1. The final tracer mixture was stored over alkaline ascorbate (0.1 M sodium ascorbate in 0.1 M NaOH) as an oxygen scrubber to avoid oxygen contamination. Based on a previous report, 10 μl tracer injections were expected to result in a final concentration of ~75 nmol L−1 methane in the incubations due to the presence of unlabelled methane in the tracer (Thamdrup et al., 2019). This would approximate the highest expected in situ methane concentrations previously reported from Golfo Dulce (up to 80 nmol L−1; Padilla et al., 2016). However, analysis of methane concentrations after the field campaign in 2018 revealed that relative to the previous reports 10 μl injections of tracer yielded lower final methane concentrations (10–30 nmol L−1) while in situ methane levels were higher (up to 1.7 μmol L−1). As a result, incubations from above the oxic–anoxic interface had 3–10 times higher methane concentrations than in situ while those from the OMZ core in 2018 had up to 50‐fold lower concentrations than in situ. In 2019 and 2020, in situ methane concentrations were analysed during each field campaign and used to calculate additions of unlabelled methane required to restore the in situ methane pool after purging. Thus, incubations from below the oxic–anoxic interface in 2019 and 2020 were carried out with in situ methane concentrations, while concentrations in most incubations in 2018 and those from above the oxic–anoxic interface in the subsequent years differed from in situ. However, as previously discussed (Thamdrup et al., 2019), the first‐order reaction kinetics of methane oxidation in seawater (Ward et al., 1987; Reeburgh et al., 1991; Thamdrup et al., 2019), which were also shown to apply in Golfo Dulce in the present study (see Results; Fig. 4), allow us to determine rate constants and to estimate in situ methane oxidation rates using the calculated rate constants and in situ concentrations of methane (see calculations below).

The addition of unlabelled methane and tracer to incubations was verified in all killed control incubations. After termination of the final timepoint (24 h), a fraction of headspace gas (0.5 ml in 2018–2019 and 1 ml in 2020) was transferred from unpurged killed controls to 3 ml glass vials (Exetainers Labco, UK) pre‐filled with N2 purged 0.5% wt./vol. zinc chloride solution, which served as storage vials. The amount of tracer (3H‐CH4) in each storage vial was measured from the headspace gas following Thamdrup et al. (2019) and was used to calculate the activity of methane dissolved in the water of the incubation vial () by applying the temperature and salinity dependent distribution coefficient (Wiesenburg and Guinasso Jr., 1979) and accounting for differences in temperature, salinity and pressure between vials. The concentration of unlabelled methane in the incubations was determined from the remaining gas in the storage vial by injecting 250 μl of the headspace gas into a gas chromatograph. For each replicate treatment, an average was calculated and used in the calculation of rate constants.

The activity of 3H in water () was measured by scintillation counting of 4 ml sample from the incubation vials after all 3H‐CH4 had been stripped from the sample by purging for 15 min with N2 (Valentine et al., 2001; Thamdrup et al., 2019). Purging was carried out within 18–26 h of terminating the final time point of the incubation experiments, with initial and final samples from the same depth being purged simultaneously. in the kinetics experiments increased approximately linearly with time (no time lags, or flattening off) and the rate constant, k, was calculated from the slope of the linear regression of plotted over time, by k =  × t −1 ×  . When first‐order reaction kinetics apply (Chan et al., 2019; Thamdrup et al., 2019), the fraction of methane oxidized per unit time equals the rate constant, thereby rates of methane oxidation were calculated as the product of k and the in situ concentration of methane [k × [CH4] , (Reeburgh et al., 1991; Valentine et al., 2001]. The LOD for rates was ~0.001 nmol L−1 d−1. A t‐test was applied to determine if rates were significantly different from zero. Data analysis was performed in R (Venables and Smith, 2003). Graphics were made using ODV (Schlitzer, 2015) and the ggplot2 package in R (Wickham, 2016).

DNA sampling, extraction and analysis

Microbial biomass for DNA analysis (16S rRNA gene amplicons) was collected and extracted following Padilla et al. (2016). Briefly, water (>1 L) was collected from Niskin bottles into pre‐rinsed Cubitainer® carboys and immediately filtered by peristaltic pump through 0.22 μm Sterivex™ filters. Filters were subsequently filled with 1.8 ml nucleic acid stabilizing buffer (RNAlater™ in 2018–2019 and PBS buffer in 2020), capped and flash‐frozen in liquid nitrogen. Filters were transferred to −18°C (up to 4 weeks) before storage at −80°C until extraction. DNA was extracted from the filters as described previously (Padilla et al., 2016; Padilla et al., 2017). In short, cells were lysed by adding 40 μl lysis buffer [40 mM EDTA, 50 mM Tris (pH 8.3), and 0.73 M Sucrose] containing 2 mg lysozyme directly to each Sterivex cartridge before sealing it and incubating for 45 min at 37°C on a rotating table. Proteinase K (1 mg in 100 μl lysis buffer, with 100 μl 20% SDS) was added and cartridges were re‐sealed and incubated rotating for 2 h at 55°C. Lysate was removed before DNA was extracted once with phenol:chloroform:isoamyl alcohol (25:24:1, pH 8) and once with chloroform:isoamyl alcohol (24:1). Finally, DNA was concentrated by spin dialysis using Ultra‐4 (100 kDa, Amicon) centrifugal filters.

Taxonomic composition was assessed by amplicon sequencing of 16S rRNA genes. In 2018 and 2020, we used the primer pairs F515 (GTGYCAGCMGCCGCGGTAA, Parada et al., 2016) and R806 (GGACTACNVGGGTWTCTAAT, Apprill et al., 2015) modified from Caporaso et al. (2011) to amplify the variable fragment (V4) of the 16S rRNA gene, whereas in 2019, we applied the bacteria‐specific PacBio primers 27F/1492R (AGRGTTYGATYMTGGCTCAG, RGYTACCTTGTTACGACTT) to synthesize near full‐length 16S rRNA amplicons (Callahan et al., 2019).

In 2018 and 2020, amplicons were generated by using GoTaq® Green PCR Master Mix (Promega) with the F151/R806 primers. Both forward and reversed primers were barcoded and appended with Illumina‐specific adapters (Kozich et al., 2013) and PCR cycle conditions were as described previously (Patin et al., 2018). Amplicons from different samples were pooled at equimolar concentrations and sequenced using paired‐end Illumina MiSeq 500 cycle kit at Georgia Tech. Sequence results were imported into DADA2 (Callahan et al., 2016; Callahan et al., 2019) and the default Illumina pipeline was used to infer amplicon sequence variants (ASV) and create ASV sequence files and ASV tables. In 2019, PCR setup and cycle conditions were as described in Earl et al. (2018), and near full‐length 16S rRNA gene amplicons were sequenced on PacBio Sequel at the Next‐Generation Sequencing facility of the Vienna Biocenter Core Facilities. PacBio tools (https://github.com/PacificBiosciences/pbbioconda) pbccs, bam2fastq and lima were used to find circular consensus sequences, convert bam to fastq files and demultiplex the near full‐length 16S rRNA sequences. ASV inference was conducted using DADA2 following the protocol for PacBio data (Callahan et al., 2019). Taxonomic classification was performed simultaneously on all three datasets (2018, 2019 and 2020) via the RDP classifier in DADA2 using the SILVA nr 138 reference taxonomy (accessed in March 2020).

Two phylogenetic trees were constructed: one of methanotrophs (Fig. 5) and one of methanogens (Fig. S3). The partial and near‐full length 16S rRNA sequences representing the identified methanotrophs (relative ASV abundance >0.05%) were aligned with a selection of closely related environmental sequences (Tavormina et al., 2013; Padilla et al., 2017) obtained from the SILVA SSU database and methanotrophs isolated in culture (Knief, 2015) using mafft‐linsi. The 16S rRNA sequences representing the identified methanogens were aligned with related environmental sequences (Nobu et al., 2016) using MUSCLE (Edgar, 2004). Based on each respective alignment, a phylogenetic tree was inferred based on maximum‐likelihood with IQTree v 1.6.12 (Nguyen et al., 2015). The model finder selected TIM3 + F + R8 and TVMe+R2 as substitution models for methanotroph and methanogen alignments respectively (Kalyaanamoorthy et al., 2017). Branch support was calculated with 500 and 1000 non‐parametric bootstrap replicates for the methanotroph and methanogen tree respectively. Amplicon data were visualized using the ggtree package in R (Yu et al. 2017). The sequence data are available at GenBank under the accession number PRJNA771562.

Author Contributions

B.T., H.G.R.S. and L.A.B. designed experiments with support from F.J.S. A.D.B, B.T., E.G.R., H.G.R.S, L.A.B, S.M. and Z.A.P. carried out sampling and lab work in Costa Rica. Sequencing was carried out by A.D.B. and Z.A.P. H.G.R.S. and S.M. analysed the data with support from B.T., C.S. and L.A.B. B.T., H.G.R.S. and L.A.B. wrote the manuscript with contributions from all authors.

Data Availability

Sequence data have been submitted to GenBank under accession number PRJNA771562. All other data are available at zenodo.org (DOI: 10.5281/zenodo.6377989).