The collection of samples used in this study was described in Padilla et al. (2016). Briefly, discrete water samples from depths of 30, 60, 90, 100, 120, and 165 m were collected using a hand-deployed Niskin bottle in late January from a site at the northern head of the coastal basin of GD, Costa Rica (Supplementary Figure 1). Microbial biomass was collected for molecular analysis by filtration of seawater (~1–3 L) through a glass fiber disc pre-filter (GF/A, 47 mm, 1.6 μm pore-size, Whatman) and a primary collection filter (Sterivex™, 0.22 μm pore-size, Millipore) using a peristaltic pump immediately after recovery. Filters were preserved in RNA stabilizing buffer (25 mM Sodium Citrate, 10 mM EDTA, 5.3M Ammonium sulfate, pH 5.2), flash-frozen in liquid nitrogen, and stored at −80°C. Replicate filters for DNA analysis were collected from similar water volumes following RNA collection, preserved with lysis buffer (50mM Tris-HCl, 40 mM EDTA, 0.73 M Sucrose), and stored at −80°C. Dissolved oxygen concentrations were measured with a Clark-type O₂ electrode (Revsbech, 1989) mounted on a hand-deployed CTD (Sea & Sun Technology).

Nitrite and dissolved methane concentrations at the time of biomass collection were first reported in Padilla et al. (2016) and are presented again here for context (Figure 2). Samples for nitrate were filtered (0.45 μm cellulose acetate) and frozen until analysis. Concentrations of nitrate + nitrite were determined using chemiluminescence after reduction to nitric oxide with acidic vanadium (III; Braman and Hendrix, 1989).

DNA was extracted from Sterivex filters (>0.2 μm biomass size fraction) using a phenol:chloroform protocol. Cells were lysed by adding lysozyme (2 mg in 40 μl of lysis buffer per filter) directly to the Sterivex cartridge, sealing the ends, and incubating for 45 min at 37°C. Proteinase K (1 mg in 100 μl lysis buffer, with 100 μl 20% SDS) was added, and cartridges were resealed and incubated for 2 h at 55°C. The lysate was removed, and the DNA was extracted once with phenol:chloroform:isoamyl alcohol (25:24:1) and once with chloroform:isoamyl alcohol (24:1) and then concentrated by spin dialysis using Ultra-4 (100 kDA, Amicon) centrifugal filters.

RNA was extracted from Sterivex filters using a modification of the mirVana™ miRNA Isolation kit (Ambion). Filter cartridges were thawed on ice, RNA stabilizing buffer was then expelled and discarded, and cells were lysed by adding Lysis buffer and miRNA Homogenate Additive (Ambion) directly to the cartridges. Following vortexing and incubation on ice, lysates were transferred to RNAase-free tubes and processed via acid-phenol:chloroform extraction according to the kit protocol. The TURBO DNA-free™ kit (Ambion) was used to remove DNA, and the extract was purified using the RNeasy MinElute Cleanup Kit (Qiagen).

Fragments of the 16S rRNA molecule were analyzed to assess the taxonomic composition in both the community DNA (gDNA) and RNA (cDNA) pools. gDNA amplicon data were generated previously (Padilla et al., 2016) and re-analyzed here to focus on the methanotrophic community. cDNA amplicon data were generated in this study, as follows. Approximately 50 ng of RNA from each sample was reverse transcribed into cDNA using the Superscript First-Strand synthesis system for RT-PCR (Invitrogen), following protocols of Campbell et al. (2009). Amplicon generation and sequencing was done using an established pipeline in our lab (e.g., Padilla et al., 2015, 2016). Briefly, amplicons were synthesized using Platinum® PCR SuperMix (Life Technologies) with primers F515 and R806, encompassing the V4 region of the 16S rRNA gene (Caporaso et al., 2011). Both forward and reverse primers were barcoded and appended with Illumina-specific adapters according to Kozich et al. (2013). Equal amounts of cDNA were used for each PCR reaction to avoid biases due to variable template concentrations (Kennedy et al., 2014). Thermal cycling involved: denaturation at 94°C (3 min), followed by 30 cycles of denaturation at 94°C (45 s), primer annealing at 55°C (45 s), and primer extension at 72°C (90 s), followed by extension at 72°C for 10 min. Amplicons were analyzed by gel electrophoresis to verify size (~400 bp, including barcodes and adaptor sequences) and purified using the Diffinity RapidTip2 PCR purification tips (Diffinity Genomics, NY). Amplicons from different samples were pooled at equimolar concentrations and sequenced on an Illumina MiSeq using a 500-cycle kit.

For both gDNA datasets (from Padilla et al., 2016) and cDNA datasets (this study), barcoded sequences were de-multiplexed, trimmed (length cutoff 100 bp), and filtered to remove low quality reads (average Phred score <25) using Trim Galore! (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Paired-end reads were then merged using FLASH (Magoč and Salzberg, 2011), with program parameters of a minimum average length of 250 bp for each read, minimum average length of 300 bp for paired read fragments, and a maximum allowable fragment standard deviation of 30 bp. The number of trimmed and merged reads per sample ranged from 3844 to 95,351 for both gDNA and cDNA amplicon pools. Chimeric sequences were detected by reference-based searches using USEARCH (Edgar, 2010) against the SILVA rRNA database 119. Amplicon sequences were then processed using “UPARSE” as implemented in USEARCH. Briefly, sequences were dereplicated, sorted by size, and clustered into operational taxonomic units (OTUs) at 97% sequence similarity. Sequences from all samples were then mapped back to the OTU file using “usearch_global,” and sample/OTU matrices were produced from the resulting output with the python script “uc2otutable.py.” Methanotrophic sequences were identified in the amplicon OTUs by BLASTN against SILVA 119, followed by keyword parsing to identify all sequences matching a known methanotroph-containing genus. Top database matches (hits) to these sequences were also identified via BLASTN against the nt database. The SINA aligner with the ARB output designation was used to align sequences representing all identified putative methanotroph OTUs in our data (n = 11; both gDNA and cDNA amplicons pooled), the sequences most closely related to these OTUs via BLASTN, and the methanotrophs isolated in culture in Knief (2015). The resulting “.arb” file was imported into ARB (Ludwig et al., 2004). A reference phylogenetic tree was compiled using near full-length sequences (1300 bp or greater) with the Feldstein correction, a custom filter based on the associated nucleotide-alignment, and bootstrap support of 1000 re-samplings. Short sequences (1300 or less), including the amplicon sequences, were then inserted into the near full-length reference tree using the parsimony tool in the ARB environment. Proportional abundances of identified methanotroph OTUs (Figure 2) were calculated after rarefaction based on the sample with the lowest number of reads (3844 reads) for both gDNA and cDNA. In rarefying the data some OTUs that were either in low abundance or unique to specific samples were removed; their contribution to the total (all samples) community is expected to be negligible.

A single depth, 90 m at the peak of measured OMZ methane oxidation rates (Padilla et al., 2016) and methanotroph abundance based on amplicon sequence data, was selected for metagenome sequencing. Community DNA was processed using the Nextera XT DNA Sample Prep kit and sequenced using a paired-end Illumina MiSeq 600 kit.

Metatranscriptomes from three depths within the anoxic zone (90, 100, 120 m) were analyzed to assess methanotroph transcriptional activity. Three datasets representing these depths were generated previously (Padilla et al., 2016) and were reanalyzed in this study. To increase the number of sequences representing the methanotroph-enriched 90 m depth, a fourth dataset was generated using an aliquot of the total RNA used in Padilla et al. (2016). For this aliquot, the Ribo-Zero™ rRNA Removal Kit for bacteria (Epicenter) was used to deplete ribosomal RNA (rRNA) sequences prior to sequencing. rRNA was not depleted from the RNA aliquots used to generate the Padilla et al. (2016) datasets. For all samples, cDNA was prepared for sequencing using the ScriptSeq™ v2 RNA-Seq Library preparation kit (Epicenter) and sequenced on an Illumina MiSeq using a 500 cycle kit. Sequencing statistics are listed in Supplementary Table 1.

Sequences were trimmed using the same methods as described above for the amplicon analysis. Quality trimmed forward and reverse sequences were merged and assembled into contigs using SPAdes 3.7.0 (Nurk et al., 2013) with the “-meta” option. The number of contigs, contig length, GC content, N50, and L50 assembly statistics were calculated with metaQUAST (Mikheenko et al., 2016; Supplementary Table 2). Contigs ≥ 500 bp were organized into genome bins based on tetranucleotide frequency and sequence coverage using MaxBin 2.0 (Wu et al., 2016). Bin completeness and contamination levels were estimated based on the representation of lineage-specific marker gene sets using CheckM (Parks et al., 2015). Assembly and genome bin statistics are in Supplementary Table 3. For each bin with ≥50% completeness and ≤10% contamination, the program Prodigal (Hyatt et al., 2010) was used to predict open reading frames (ORFs). Predicted ORFs were queried via BLASTP against the NCBI-nr database (April 2016; criteria for a significant match: amino acid ID ≥40% over ≥70% of the alignment, and bit score ≥50) and assigned to a prokaryotic KEGG Orthology (KO) identifier (Kanehisa and Goto, 2000) using the MetaGenome Analyzer 5 (MEGAN5; Hudson et al., 2011), with taxonomic classification in MEGAN5 via the Lowest Common Ancestor (LCA) algorithm based on the NCBI taxonomy using the default settings. Annotations are reported in Supplementary Tables 4, 5.

Based on high completeness and low contamination (see below), a single bin (bin 010) taxonomically affiliated with OPU3-like Methylococcales bacteria was selected for comparative analysis. The identity of this bin was further confirmed based on phylogenetic analysis of 107 concatenated single copy marker genes present in this bin and in the genomes (n = 36) of known methanotrophic taxa from the Gammaproteobacteria, Alphaproteobacteria, Verrucomicrobia, and the NC10 candidate phylum. Marker genes were identified using Hidden Markov Models (HMMs) via HMMER3 (http://hmmer.org/; Finn et al., 2011) with default settings. Identified genes were aligned using clustalW (Sievers et al., 2011), then concatenated using the alignment tool “Aln.cat.rb” in the enve-omics package (Rodriguez-R and Konstantinidis, 2016). The resulting alignment was inspected manually for errors and then used to generate a maximum likelihood phylogeny inferred with the Dayhoff substitution model with 100 bootstrap iterations using MEGA6 (Tamura et al., 2013). For each genome in the phylogeny, the presence/absence of key genes of dissimilatory and assimilatory nitrogen metabolism was evaluated via BLASTX using the KEGG Automated Annotation Server (KAAS). Each genome was individually assessed with the bi-directional best-hit option on default settings and the results were manually searched for subunits of nitrogen metabolism genes, with the results (presence/absence) mapped onto the concatenated marker gene phylogeny. Gene content in bin 010 relative to the 36 other methanotroph genomes included in this analysis was evaluated via reciprocal BLASTP using the two-way amino-acid identity script in the enve-omics packages (Konstantinidis and Tiedje, 2005) with a bit score threshold of 50, and match criteria of >35% amino acid identity over >65% of each gene (Supplementary Table 4).

Phylogenetic analysis was used to further examine the taxonomic affiliation of key genes of the denitrification process recovered in bin 010: narG and nirK. The amino acid sequence from each gene was aligned against a representative sequence set, identified based on top BLASTP results, using MUSCLE with default settings (Edgar, 2004). The alignments were manually inspected and used to generate maximum likelihood phylogenies using MEGA6.0 with the Dayhoff substitution model and 1000 bootstrap iterations. To further characterize taxonomic affiliations, the taxonomic annotation of all genes occurring on the same contig as key marker genes of denitrification and methanotrophy (nar, nir, pmo) was determined via BLASTX against NCBI-nr (Supplementary Figure 2).

Reads were filtered by quality and merged as described above. The program ribopicker (Schmieder et al., 2012) was used to remove rRNA sequences from metatranscriptome datasets. To characterize the transcript pool affiliated with the OPU3 genomic bin described above, transcripts were mapped to bin 010 ORFs using BLASTX with match criteria of >90% amino acid identity over >60% of the transcript length and bit score >50. To compare transcript levels among bin 010 ORFS, the number of transcripts with significant BLASTX matches per ORF was normalized to variation in ORF length, and then expressed as a proportion of total transcripts mapping to bin 010. To also characterize the transcript pool at the community level, merged non-rRNA datasets were queried via DIAMOND (Buchfink et al., 2015) against the NCBI-nr database (April 2016) using the sensitive search setting. DIAMOND-identified protein-coding transcripts were assigned to functional categories as above based on Kegg Ortholog (KO) identifiers and taxonomically classified according to the NCBI taxonomy using the LCA algorithm in MEGAN5. Transcript counts per KO were normalized to the total number of mRNA transcripts assigned to a prokaryotic KO. To show variations in contributions of bin 010 to community mRNA transcription with depth, the number of transcripts with significant BLASTX matches to bin 010 ORFs was expressed as a proportion of total prokaryotic mRNA transcripts assignable to a KO.

All sequence data, including the bin 010 assembly, have been submitted to the Sequence Read Archive at NCBI under the following BioProject ID's: PRJNA328797 and PRJNA277357.

Deep-coverage sequencing of 16S rRNA gene (gDNA) and transcript (cDNA) amplicons revealed an abundant, diverse, and active methanotroph community in the Golfo Dulce OMZ. All sequences that matched putative aerobic methanotrophs were classified as type Ia Gammaproteobacteria methanotrophs. Across all depths, 11 OTUs from the gDNA dataset were identified as Methylococcales by BLASTN. Eight of these OTUs remained after rarefaction (3844 reads) of the total gDNA and cDNA amplicon pools. Phylogenetic analysis classified these eight OTUs as most closely related to methanotrophs of the OPU3 and OPU1 clusters and the genus Methyloprofundus (Figure 1), with six of the eight belonging to the OPU3 clade. A single OPU3 OTU (GD_7) dominated the methanotroph sequence dataset, constituting ~35–75% of the Methylococcales pool in both the gDNA and cDNA datasets at all depths (Figure 2). Only one OTU from the OPU1 clade (GD_791) was detected and was confined primarily to the 60 m sample where O₂ concentration was ~25 μM. Together, these eight OTUs, which collectively represented 1.8% of the total gDNA combined over all depths, were undetected or at low abundance (0.4%) above or within the oxycline (30 and 60 m, respectively) but increased to 4–5% of the gDNA pool in the upper anoxic zone (90 and 100 m) where O₂ concentrations fell below detection, <50 nM, and nitrite accumulated (Figure 2). Notably, OTU GD_7, the dominant OPU3 bacterium, constituted 4% of the total microbial community (gDNA pool) at 90 m. These OTUs then decreased to <1% of the community deeper in the anoxic zone at 120 m, and at 165 m where nitrate and nitrite were depleted and methane accumulated. The proportional representation of these eight OTUs was ~2-fold higher in the cDNA compared to the gDNA datasets at all depths, indicating a transcriptionally active Methylococcales community.

Metagenomic analysis confirmed the potential for coupled methanotrophy and denitrification in the Golfo Dulce OMZ. Genomic DNA from within the secondary nitrite maximum at 90 m was deeply sequenced, assembled into contigs, and binned into OTUs resulting in 52 draft genomes (bins) ranging from 12.4 to 100% completeness and 0–86% contamination (Supplementary Table 3). Of these, bin 010 contained a single near-full length (1405 bp) 16S rRNA gene fragment identical to that of OTU GD_7 identified in the amplicon analysis (above). The contig containing this fragment consisted solely of the 16S rRNA gene and exhibited an average per-base coverage of 6.6, nearly identical to the bin-wide average (all contigs) of 6.7. The bin contamination level was 6% based on duplicated marker genes, a level consistent with that reported in other studies drawing conclusions of genome content based on binning [see Sekiguchi et al., 2015 (6–7%); Campanaro et al., 2016 (3–5%); Güllert et al., 2016 (5–10%)]. On average, duplicate marker genes in bin 010 shared 81% amino acid identity, within the genome-wide range observed for bacterial strains categorized in the same species (Konstantinidis and Tiedje, 2005), suggesting that contaminant sequences in this bin are from closely related taxa.

Given the high abundance of OTU GD_7 at >4% in both the gDNA and cDNA amplicon pools (Figure 2), suggesting a substantial contribution to OMZ community processes, bin 010 was selected for in-depth characterization. The genome represented by this bin, hereafter designated GD_7, was estimated to be 95.3% complete, with 2011 coding sequences distributed over 305 contigs and a total of 2,106,486 bp (Supplementary Table 2). The estimated GD_7 genome size, at 2.2 Mbp, is consistent with that of pelagic free-living bacteria (Raes et al., 2007; Shi et al., 2011), but substantially smaller than that of most other gammaproteobacterial methanotrophs (~4–5 Mbp; Figure 3). LCA-based taxonomic classification of GD_7 coding sequences indicated an affiliation with the order Methylococcales. Phylogenetic analysis of 107 housekeeping genes (concatenated) from available methanotroph genomes confirmed GD_7 as basal but most closely related to the type-Ia methanotrophs of the Methylococcales, a diverse group of bacteria with some members shown experimentally to couple methane oxidation to denitrification (Figure 3; Kits et al., 2015a,b).

Below, we discuss a subset of key functional genes recovered in the GD_7 bin. As the inferred genome is fragmentary and exhibits 6% contamination (above), confirmation that the recovered genes are co-localized on the same chromosome in OPU3 requires further analysis of a pure culture or a closed genome. Similarly, as the recovered GD_7 genome is not 100% complete, we cannot definitively conclude that specific genes are absent from this genome. However, the gene content described below, when contextualized relative to that of other methanotrophs, prior results, and the environmental conditions of Golfo Dulce, is strongly suggestive of key metabolic functions in GD_7.

The genome content of GD_7 indicates the potential for methanotrophic denitrification. GD_7 contains genes required for energy generation by methane oxidation, including those encoding the three subunits of particulate methane monooxygenase (pmoCAB) for methane oxidation to methanol, which together co-occur on a contig with diverse other genes having close homologs in other gammproteobacterial methanotrophs (inferred from top BLASTX matches; Supplementary Figure 2, Supplementary Table 4). Also present are genes for the lanthanide-dependent (xoxF) methanol dehydrogenase involved in methanol oxidation to formaldehyde, NAD(P)-dependent methylene-H₄MPT dehydrogenase (mtd) involved in formaldehyde oxidation to formate, and formate dehydrogenase (fdh) for formate oxidation to CO₂. The presence of XoxF type methanol dehydrogenases, but absence of the calcium-dependent MxaF type, has also been reported for pelagic non-methanotrophic methylotrophs with relatively small genomes (Giovannoni et al., 2008) and is consistent with recent studies suggesting the prevalence of XoxF compared to MxaF in natural systems (Keltjens et al., 2014; Ramachandran and Walsh, 2015).

GD_7 also contains two operons encoding a dissimilatory nitrate reductase (nar), as well as genes for synthesis of the essential nar molybdenum cofactor. Of the downstream genes of the denitrification pathway, only the nirK gene encoding nitric oxide (NO)-producing nitrite reductase was detected in GD_7. The two nar operons, as well as nir, were recovered on contigs alongside other functional genes with homology to those of gammaproteobacterial methanotrophs (Supplementary Figure 2). Phylogenetic analysis of marker genes nirK and narG, with the latter encoding the alpha subunit of Nar, also supports a methanotroph affiliation. One of the narG copies (ORF 1321) clusters with high support in a clade containing denitrifying gammaproteobacterial methanotrophs, including Methylomonas denitrificans (Figure 4A). Likewise, GD_7 nirK clusters with a sequence from the gammaproteobacterial Methylothermaceae and that of the nitrite-oxidizing bacterium Nitrospina gracilis, although the function of NirK in Nitrospina remains unclear (Lüker et al., 2013; Figure 4B). In contrast, the second nar operon may have an origin outside of the gammaproteobacteria, as analysis of the corresponding narG copy (ORF 1413) reveals a closer affiliation to narG from a genus of Chlorobi bacteria (Ignavibacterium) containing facultatively anaerobic chemoheterotrophic members (Figure 4A).

These results suggest that certain genes for partial denitrification in GD_7 share ancestry with those of gammaproteobacterial methanotrophs, while others, notably ORF 1413, may have been acquired horizontally (Figure 4A). Horizontal acquisition of denitrification genes in methanotrophs may be common, as the presence of nar and nir genes is independent of methanotrophic phylogenetic histories (Figure 3). Indeed, the co-occurrence of divergent denitrification genes in the same genome is not unusual. Many bacteria, for example, encode divergent nar copies (Philippot, 2002; Tsementzi et al., 2016), potentially as an adaptation to variable oxygen conditions (Iobbi-Nivol et al., 1990). Although multiple nar operons have yet to be reported in other methanotrophs, a marine methylotroph was shown to encode two nar enzymes, with one copy playing a role in regulating the other depending on oxygen availability (Mauffrey et al., 2015). Some denitrifying methanotrophs also encode the perisplasmic dissimilatory nitrate reductase (Nap), hypothesized to function primarily under nitrate-limiting conditions (Figure 3; Ferguson and Richardson, 2004); however, nap genes were not detected in GD_7. In contrast to Methylomonas denitrificans FJG1, Methylomicrobium album BG8, and Candidatus “Methylomirablis oxyfera,” GD-7 does not appear to harbor an NO-reductase (nor), and is therefore likely incapable of producing N₂O as an intermediate or end-product (Figure 3; Kits et al., 2015a,b). In M. denitrificans FJG1, M. album BG8, and potentially GD-7, nitrate is likely utilized for respiratory purposes. This is in contrast to NC10 members, which likely utilize nitric oxide (derived from nitrate and nitrite reduction) for dismutation with concomitant production of intracellular oxygen. Rather, our data identify a partial methanotrophic denitrification pathway in GD_7, with NO as the likely respiratory end product.

GD_7 also contains genes encoding key enzymes of C₁ carbon assimilation from methane, including those of the ribulose monophosphate (RuMP) pathway for incorporation of methanol-derived formaldehyde and regeneration of ribulose-5-phosphate (Supplementary Figure 3). Key enzymes of the serine pathway for formaldehyde assimilation were also detected (Table 1, Supplementary Tables 4, 5), with 10 of the serine cycle genes recovered on the same contig. Diagnostic enzymes of the ethylmalonyl-CoA (EMC) pathway (ethylmalonyl-CoA mutase, crotonyl-CoA reductase/carboxylase), used by some serine cycle-utilizing methanotrophs for glyoxylate regeneration (Chistoserdova, 2011), were not detected in GD_7. Neither were isocitrate lyase and malate synthase, suggesting glyoxylate generation does not occur via the glyoxylate cycle in GD_7. This pattern, the presence of serine cycle genes and absence of an EMC pathway and glyoxylate shunt, has been reported in other type 1 methanotrophs (Poehlein et al., 2013; de la Torre et al., 2015). In GD_7, all enzymes of the Citric Acid (TCA) cycle are present, as are genes of glycolysis and for pyruvate dehydrogenase for acetyl-CoA generation (Supplementary Table 4). However, GD_7 appears to lack the genes encoding the Entner-Doudoroff Pathway, which would be atypical of Type I methanotrophs that have been hypothesized to depend primarily on Entner-Doudoroff cleavage reactions over glycolysis (Kalyuzhnaya et al., 2013). However, we cannot rule out that these or other genes were not detected due to incomplete coverage (~95%) of this genome bin. Together, these patterns suggest C₁ incorporation from methane via the RuMP pathway, with acetyl-CoA synthesis via pyruvate for entry into the TCA cycle for biosynthesis.

Reciprocal BLASTP against all methanotrophic bacteria in Figure 3 identified 102 protein-coding genes in GD_7 that lack close homologs in other methanotroph genomes (based on the imposed BLAST criteria: bit score >50, with >35% amino acid identity across >65% of the gene). Of this divergent gene set, 47 encode hypothetical proteins (Supplementary Table 4). The remaining annotated genes are associated with diverse functions, including flagellum synthesis or regulation, phage tail synthesis, and diverse functions of amino acid or peptide metabolism, as well as an arginase potentially involved in urea generation. Genes in this divergent set may be candidates for follow-up studies to better understand niche-specific adaptations in GD_7.

To characterize the Methylococcales-like transcript pool, we examined transcripts from 90, 100, and 120 m that mapped with high identity to GD_7. To increase coverage across the bin, we analyzed an additional 90 m dataset from which rRNA was removed prior to sequencing. This step increased data yield but altered the abundance of mapped transcripts (expressed as a % of total prokaryotic mRNA; Figure 5), suggesting that rRNA subtraction may also have removed certain mRNA transcripts. This procedure does not affect our conclusions, however, as our goal was to characterize the range of genes actively transcribed by GD_7; variation in transcription among genes is presented in only a general sense (e.g., Figure 5). Nonetheless, the rRNA-subtracted dataset is not considered when evaluating changes across depths (below).

Metatranscriptome analysis confirmed that Methylococcales bacteria were transcriptionally active in the Golfo Dulce OMZ. Across all datasets, a total of 24,989 reads representing 764 genes were recruited via BLASTX to GD_7 contigs with high identity (>90% AAI; bit scores >50) across 60% of the transcript fragment (Supplementary Table 5). The proportional representation of these mapped transcripts peaked at 90 and 100 m at ~0.5% of total mRNA reads (non-rRNA-subtracted dataset) before decreasing to <0.1% deeper in the anoxic zone (120 m sample; Figure 5B). Indeed, at 90 and 100 m, the vast majority (~90%) of recruited transcripts share >98% AAI with genes of GD_7 (Figure 5A; Supplementary Figure 4). In contrast, at 120 m, the majority of mapped reads share <94% AAI. These patterns suggest discrete OPU3-like populations based on depth, with GD_7 confined primarily to the upper OMZ where both nitrite and nitrate are enriched compared to greater depths in the anoxic zone (Figure 2).

Transcripts involved in methanotrophy and denitrification were among the most abundant of those recruiting to GD_7, based on the mRNA-enriched 90 m dataset (Supplementary Figures 4, 5; Supplementary Table 5). Mapped transcripts included all genes of the aerobic methane oxidation pathway (Table 1), with pmoA and pmoC being particularly abundant, as well as genes of the RuMP pathway. Denitrification transcripts included narG and nirK, both of which were detected at 90 and 100 m, but absent from the 120 m dataset. Interestingly, only transcripts mapping to narG ORF 1413 were detected (Supplementary Table 5), suggesting differential regulation of Nar expression and therefore the potential for functional variation between the two enzyme variants (ORF 1321, 1413). Terminal oxidase genes, though present in GD_7, were not detected in the mapped transcript pool, suggesting that any oxygen being utilized is for methane oxidation rather than respiration. Genes for fatty acid metabolism were also among the most abundant transcripts recruiting to GD_7, as were those mediating flagellar biosynthesis and regulation. The high transcription of flagellar genes is consistent with prior reports of the importance of motility in type 1 methanotrophs from low oxygen environments (Danilova et al., 2016), suggesting a need to traverse redox gradients to access zones of optimal substrate (e.g., methane) or oxygen conditions, potentially including those at the micron scale on suspended or sinking particles. Enzymes involved in amino acid metabolism were also abundant in the GD_7 transcript pool (Supplementary Table 5, Supplementary Figure 5), with proline dehydrogenase being among the most highly transcribed GD_7-specific protein at 90 m. This enzyme is involved in proline catabolism, producing reductant in the form of reduced co-enzyme Q10 (i.e., ubiquinol). Although use of proline for reductant has been reported for diverse bacteria under both oxic and anoxic conditions (Deutch and Soffer, 1975; Barker, 1981), the potential contribution of amino acid catabolism to reductant pools in methanotrophs, including OPU3 members, remains uncharacterized. Genes encoding enzymes of the methylcitrate pathway (MCA) involved in propionate catabolism were also highly transcribed (Supplementary Table 5, Supplementary Figure 5). These included methylcitrate synthase, methylcitrate dehydrogenase and aconitate hydratase, although the MCA gene encoding methylisocitrate lyase was not detected. MCA genes and transcripts have been reported in aquatic betaproteobacterial methylotrophs, raising the hypothesis that propionate, which can be produced via demethylation of compounds such as dimethylsulfoniopropionate, may be an important carbon substrate in these bacteria (Kalyuzhnaya et al., 2008, 2010). The potential for propionate catabolism in OPU3 should be further evaluated.

Taxonomic classification of the bulk metatranscriptome based on lowest common ancestry (LCA) corroborated the transcriptional activity of a Methylococcales community in the OMZ. The majority (~80%) of LCA-identified Methylococcales transcripts could not be assigned to a genus (Supplementary Figure 6). Of those transcripts that could be classified, most are affiliated with Methylomonas, Methylobacter, or Methylomicrobium. The proportional abundance of total Methylococcales transcripts across depths was nearly identical to that of the transcript pool mapping to GD_7 (Figure 5B, Supplementary Figure 6), indicating that the vast majority of Methylococcales transcripts were affiliated with this OPU3 taxon and that its contribution to community transcription peaked at 90 and 100 m just below the oxycline. As in the GD_7 transcript pool, biochemical functions highly represented in the bulk Methylococcales transcripts were predominantly associated with methane oxidation and C₁ assimilation, with transcripts encoding pmoCAB among the most abundant (data not shown).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.