Environmental monitoring and sample collection were carried out monthly aboard the MSV John Strickland at station SI03 (48° 35.500 N, 123° 30.300 W) in Saanich Inlet, B.C. From February 2008 to July 2010 a total of 288 water samples were collected at 16 high-resolution depths (10, 20, 40, 60, 75, 85, 90, 97, 100, 110, 120, 135, 150, 165, 185, and 200 meters) spanning oxic (>90 μmol O₂ kg⁻¹), dysoxic (90-20 μmol O₂ kg⁻¹), suboxic (20-1 μmol O₂ kg⁻¹) anoxic (<1 μmol O₂ kg¹) and sulfidic water column compartments (Wright et al., 2012). Samples were processed and analyzed as previously reported for the Saanich Inlet time-series (Zaikova et al., 2010; Capelle et al., 2015; Torres-Beltrán et al., Unpublished data). Briefly, water samples were collected from Niskin or Go-Flow bottles for dissolved O₂ and CH₄, nutrients [Nitrate (NO₃), Nitrite (NO₂), Ammonium (NH₄), Silicon dioxide (SiO₂), Phosphate (PO₄), Hydrogen sulfide (H₂S)], DNA, and RNA. In addition, conductivity, temperature, and depth were measured using a Seabird SBE19 CTD-device (Sea-Bird Electronics Inc., Bellevue USA), with a PAR and O₂ sensor attached. CTD was also used to measure dissolved O₂ throughout the water column. Dissolved gases and nutrient measurement protocols have been previously reported for the Saanich Inlet time-series (Zaikova et al., 2010; Capelle et al., 2015; Torres-Beltrán et al., Unpublished data), and data is available through Dryad Digital Repository (www.dryad.org).

Biomass to generate full-length small subunit (SSU) rRNA gene clone library sequences and pmoA gene library sequences (Stilwell, 2007), and metagenomic datasets was collected on February 2006 and February 2010 respectively, at six depths (10, 100, 120, 135, 150, and 200 meters) and filtered with an in-line 2.7 μm GDF glass fiber pre-filter onto a 0.22 μm Sterivex polycarbonate cartridge filter. Biomass to generate SSU rRNA pyrotag datsets for microbial community composition profiling (2L) and filtered directly onto a 0.22 μm Sterivex polycarbonate cartridge filter from high-resolution depths collected between May 2008 and July 2010. Biomass for RNA analysis (2L) was collected on February 2010 at six depths (10, 100, 120, 135, 150, and 200 m), and filtered with an in-line 2.7 μm GDF glass fiber pre-filter onto a 0.22 μm Sterivex filter within 20 min of shipboard sample collection.

Environmental DNA was extracted from Sterivex filters as previously described (Wright et al., 2009; Zaikova et al., 2010). Briefly, after defrosting Sterivex on ice, 100 μl lysozyme (0.125 mg ml⁻¹; Sigma) and 20 μl of RNAse (1 μl ml⁻¹; ThermoFisher) were added and incubated at 37°C for 1 h with rotation followed by addition of 50 μl Proteinase K (Sigma) and 100 μl 20% SDS and incubated at 55°C for 2 h with rotation. Lysate was removed by pushing through with a syringe into 15 mL falcon tube (Corning) and with an additional rinse of 1 mL of lysis buffer. Filtrate was subject to chloroform extraction (Sigma) and the aqueous layer was collected and loaded onto a 10K 15 ml Amicon filter cartridge (Millipore), washed three times with TE buffer (pH 8.0) and concentrated to a final volume of between 150–400 μl. Total DNA concentration was determined by PicoGreen assay (Life Technologies) and genomic DNA quality determined by visualization on 0.8% agarose gel (overnight at 16 V).

Total RNA was extracted from Sterivex filters using the mirVana Isolation kit (Ambion) (Shi et al., 2009; Stewart et al., 2010) protocol modified for sterivex filters (Hawley et al., Unpublished data). Briefly, after thawing the filter cartridge on ice RNA later was removed by pushing through with a 3 ml syringe followed by rinsing with an additional 1.8 mL of Ringer's solution and incubated at room temperature for 20 min with rotation. Ringer's solution was evacuated with a 3 ml syringe followed by addition of 100 μl of 0.125 mg ml⁻¹ lysozyme and incubated at 37°C for 30 min with rotation. Lysate was removed from the filter cartridge and subjected to organic extraction following the mirVana kit protocol. DNA removal and clean up and purification of total RNA were conducted following the TURBO DNA-free kit (ThermoFisher) and the RNeasy MinElute Cleanup kit (Qiagen) protocols respectively.

To initially assess microbial community diversity in the Saanich Inlet water column full-length bacterial small subunit ribosomal (SSU) rRNA gene sequences datasets were generated from the February 2006 LV samples as previously described (Zaikova et al., 2010). Briefly, DNA extracts from 10, 100, 120, and 200 m samples were amplified using SSU rRNA primers targeting the bacterial domain: B27F (5′- AGAGTTTGATCCTGGCTCAG) and U1492R (5′-GGTTAC CTTAGTTACGACTT) under the following PCR conditions: 3 min at 94°C followed by 35 cycles of 94°C for 40 s, 55°C for 1.5 min, 72°C for 2 min and a final extension of 10 min at 72°C. Each 50 ml reaction contained 1 μl of template DNA, 1 μl each 10 mM forward and reverse primer, 2.5 U Taq (Qiagen), 5 ml 10 mM deoxynucleotides, and 41.5 μl 1X Qiagen PCR Buffer. Clone library construction and screening. The SSU rRNA gene amplicons were visualized on 1% agarose gels in 1X TAE and purified using the MinElute PCR Purification Kit (Qiagen) according to the manufacturer's instructions. Approximately 4 μl of each purified SSU rRNA gene product was cloned into a pCR4-TOPO vector using a TOPO TA cloning kit for sequencing (Invitrogen) and transformed by chemical transformation into Mach- 1-T1R cells according to the manufacturer's instructions. Transformants were transferred to 96-well plates containing 180 mlLB_kan50 and 10% glycerol and grown overnight at 37°C prior to storage at −80°C. Cloned inserts were amplified directly from glycerol stocks with M13F (5′-GTAAAACGACGGCCAG) and M13R (5′- CAGGAAACAGCTATGAC) primers using the SSU rRNA gene PCR protocol. Bidirectional end sequencing was performed on a Sanger platform at the Department of Energy Joint Genome Institute (DOE-JGI; Walnut Creek, CA).

To survey microbial community structure and dynamics throughout the complete water column profile over time, DNA extracts from the samples from May 2008 to July 2010 were used to generate SSU rRNA gene pyrotag datasets. Pyrotag libraries were generated by PCR amplification using multi-domain primers targeting the V6-V8 region of the SSU rRNA gene (Allers et al., 2013): 926F (5′-cct atc ccc tgt gtg cct tgg cag tct cag AAA CTY AAA KGA ATT GRC GG-3′) and 1392R (5′-cca tct cat ccc tgc gtg tct ccg act cag- <XXXXX>-ACG GGC GGT GTG TRC-3′). Primer sequences were modified by the addition of 454 A or B adapter sequences (lower case). In addition, the reverse primer included a 5 bp barcode designated <XXXXX> for multiplexing of samples during sequencing. Twenty-five microliter PCR reactions were performed in triplicate and pooled to minimize PCR bias. Each reaction contained between 1 and 10 ng of target DNA, 0.5 μl Taq DNA polymerase (Bioshop inc. Canada), 2.5 μL Bioshop 10 x buffer, 1.5 uL 25 mM Bioshop MgCl₂, 2.5 μL 10 mM dNTPs (Agilent Technologies) and 0.5 μL 10 mM of each primer. The thermal cycler protocol started with an initial denaturation at 95°C for 3 min and then 25 cycles of 30 s at 95°C, 45 s at 55°C, 90 s at 72°C and 45 s at 55°C. Final extension at 72°C for 10 min. PCR products were purified using the QiaQuick PCR purification kit (Qiagen), eluted elution buffer (25 μL), and quantified using PicoGreen assay (Life Technologies). SSU rRNA amplicons were pooled at 100 ng DNA for each sample. Emulsion PCR and sequencing of the PCR amplicons were sequenced on Roche 454 GS FLX Titanium at the DOE-JGI, or the McGill University and Génome Québec Innovation Center.

Pyrotag sequences were processed using the Quantitative Insights Into Microbial Ecology (QIIME) software package (Caporaso et al., 2010). To minimize the removal of false positive reads all 3,985,489 pyrotag sequences generated from the 288 samples were clustered together. Reads with a length shorter than 200 bases, ambiguous bases, and homopolymer sequences were removed prior to chimera detection. Chimeras were detected and removed using chimera slayer provided in the QIIME software package. Sequences were then clustered into operational taxonomic units (OTUs) at 97% identity using uclust with average linkage algorithm. Prior to taxonomic assignment singleton OTUs (OTUs represented by one read) were omitted, leaving 69,051 OTUs. Representative sequences from each non-singleton OTU were queried against the SILVA database (Pruesse et al., 2007) and the Greengenes database (Desantis et al., 2006) using BLAST (Altschul et al., 1990).

Extracted total DNA from the February 2006 samples were used to identify bacterial particulate methane monooxygenase (pMMO) subunit β (pmoA) presence in the Saanich Inlet water column. Briefly, pmoA sequences were PCR amplified (Stilwell, 2007) using gene-specific forward and reverse primers A189F (5′-GGNGACTGGGACTTCTGG) (Holmes et al., 1995) and mb661R (5′- CCGGMGCAACGTCYTTACC) (Costello and Lidstrom, 1999) and the following PCR profile: 30 cycles of 96°C for 25 s, 54°C for 45 s, and 72°C for 50 s. Each 50 ul reaction contained l μl of template DNA, 2 μl each 10 uM forward and reverse primer, 2.5U Taq (Qiagen), 4 μl 10 mM deoxynucleotides, and 40.5 μl lx Buffer (QiagenTaqPolymeraseKit). pmoA amplicons were purified using the MinElute PCRPurificationKit (Qiagen, CA), cloned into a pCR4-TOPO vector using a TOPO T A cloning kit for sequencing (Invitrogen, Carlsbad C A), and transformed by chemical transformation into Mach-1-T1R cells according to the manufacturer's instructions. Cloned inserts were amplified directly from glycerol stocks for fingerprint screening using the common 4-base cutter Rsa I (Invitrogen, CA). Restriction patterns were visually inspected and unique patterns selected for Sanger sequencing through the McGill University and Genome Quebec Innovation Centre (Montreal, Quebec, Canada) (Stilwell, 2007). Four libraries were constructed from 10, 100, 120, and 200 m, from which a total of 33 representative sequences were identified.

Extracted total DNA (6 samples) and RNA (6 samples) from February 2010, corresponding to the 10, 100, 120, 135, 150, and 200 m depth intervals, were used to generate metagenomic and metatranscriptomic datasets at the DOE-JGI following the protocols for library production and sequencing and assembly previously described for the Saanich Inlet time-series (Hawley et al., Unpublished data).

A total of 6 assembled metagenomes and 6 assembled metatranscriptomes were analyzed using MetaPathways V2.5.1, an open source pipeline for predicting reactions and pathways using default settings (Konwar et al., 2013) (https://github.com/hallamlab/metapathways2/wiki). For each gene, reads per kilobase per million mapped (RPKM) was calculated as a proportion of the number of reads mapped to a sequence section, normalized for sequencing depth and ORF length (Konwar et al., 2015). RPKM values and relative abundance to the total number of reads were used to describe the abundance of pmoA genes and transcripts.

To generate a reference phylogenetic tree for methanotrophic bacteria, full-length SSU rRNA gene sequences from Saanich Inlet, and diverse environmental reference sequences affiliated with methanotrophic bacteria were first aligned and compared. Full-length SSU rRNA gene sequences from Saanich Inlet were edited manually using Sequencher software V4.1.2 (Gene Codes Corporation) and imported into the full-length SSU SILVA database (www.arb-silva.de) and aligned to the closest relative. Sequences affiliated with known methanotrophs were extracted from the dataset. Methanotroph reference sequences from Saanich Inlet in addition to 53 reference sequences for methanotrophic bacteria including cultured type I and II methanotrophs, OMZ representatives, mussel symbionts and envioronmental clones were clustered at 97% identity using mothur v.1.19.0 (Schloss et al., 2009). A total of 88 sequences affiliated with Methylococcales were recovered representing 1.3% from the total Saanich Inlet full-length SSU rRNA sequences generated. Sequences clustered at 97% identity resolved into 11 distinct clusters, 5 of which contained most of the identified sequences (93%). Representative sequences for the most abundant clusters revealed 4 subgroups with phylogenetic similarity to environmental representatives of type I methanotrophs: Methylococcaceae (Mou et al., 2008), putative methanotrophic group OPU3 (Hayashi et al., 2007; Tavormina et al., 2010, 2013), environmental seafloor clones (Santelli et al., 2008), and methanotrophic symbionts (Streams et al., 1997; Dubilier et al., 2008; Petersen and Dubilier, 2009). Representative sequences for the most abundant (represented by more than one sequence) clusters were identified using the get.oturep command in mothur and were included in the phylogenetic tree. Representative sequences were aligned using the SSU SILVA database and imported into the ARB software (Ludwig et al., 2004) for tree distance matrix and alignment generation using the ARB parsimony tool. ARB sequences were exported to Mesquite (V.2.0) and edited manually. A maximum likelihood phylogenetic tree was inferred by PHYML (Guindon et al., 2005) using an GTR model of nucleotide evolution where the parameter of the gamma distribution, the proportion of invariable sites and the transition/transversion ratio were estimated for each data set. The confidence of each node was determined by assembling a consensus tree of 1000 bootstrap replicates.

To further resolve diversity of methanotrophic bacterial OTUs, we recruited pyrotag OTU sequences to full-length SSU rRNA gene sequences described above. Pyrotag sequences affiliated with methanotrophs based on BLAST-comparison in QIIME were re-clustered with SSU rRNA gene tree reference sequences using a 97% identity cut-off in mothur (Schloss et al., 2009). In addition, blastn was used to query representative pyrotag sequences from clusters against full-length SSU rRNA reference tree sequences. Only hits with a perfect match across the full length of a query sequence were retrieved, and the number of pyrotags mapping to all sequences in each cluster was summed. Clusters represented by one pyrotag sequence were not used in downstream analyses. Representative pyrotag sequences were aligned in ARB to reference tree sequences, imported to Mesquite for manual edition, and finally, included in the phylogenetic tree inferred by PHYML using the parameters detailed above.

To generate a reference phylogenetic tree for PmoA, conceptually translated and annotated ORFs from the metagenomic and metatranscriptomic datasets were manually extracted from the functional annotation table <ORF_annotation_table.txt> in the <results/annotation table> output directory. Sequences were aligned and compared to diverse environmental and reference PmoA sequences. A total of 60 PmoA sequences, retrieved from metagenomic (34 sequences) and metatranscriptomic (26 sequences) datasets, were clustered over a range of identity thresholds using the UClust algorithm (USEARCH V6.0) with 52 reference sequences including Saanich Inlet pmoA gene libraries, cultured and environmental sequences affiliated with Type I and II methanotrophs, and novel pmoA phylotypes found within the SAR324 clade, Verrucomicrobia and Candidate Methylomirabilis oxyfera NC10. Reference sequences also included ammonium monooxygenase subunit α (AmoA). The 97% identity threshold was selected based on resolution of the OPUs and symbiont groups. Cluster representative sequences were aligned using the Multiple Sequence Comparison by Log- Expectation (MUSCLE) method (EMBL-EBI), and was manually curated in Mesquite. A maximum likelihood phylogenetic tree was inferred by PHYML (Guindon et al., 2005) using an WAG model of amino acid evolution where the parameter of the gamma distribution, the proportion of invariable sites and the transition/transversion ratio were estimated. The confidence of each node was determined by assembling a consensus tree of 1000 bootstrap replicates.

Pyrotag datasets were normalized to the total number of reads per sample, and environmental parameter data were transformed to the same order of magnitude so that each variable had equal weight. Hierarchical cluster analysis (HCA) was conducted to identify groups associated with discrete water column compartments. In addition, OTUs were correlated using nonmetric multidimensional scaling (NMDS) with environmental parameters. Hierarchical cluster and NMDS analyses of microbial community compositional profiles were done using the pvclust (Suzuki and Shimodaira, 2015) and MASS (Venables and Ripley, 2002) packages in the R software (Rcoreteam, 2013) with Manhattan Distance measures, and statistical significance to the resulting clusters was computed as bootstrap score distributions with 1000 iterations and NMDS stress value ≤0.05.

Multi-level indicator species analysis (ISA) using the indicspecies package (De Cáceres and Legendre, 2009) in the R software (Rcoreteam, 2013) was performed to identify OTUs specifically associated with different water column compartments defined by HCA. The ISA/multi-level pattern analysis calculates p values with Monte Carlo simulations and returns indicator values (IV) and p-values with α ≤ 0.05. The IVs range between 0 and 1, where indicator OTUs considered in the present study for further community analysis shown an IV ≥ 0.7 and p-value ≤ 0.001.

A multivariate regression analysis (Fox and Weisberg, 2011) was conducted on time-series pyrotag data to infer significant correlations between OTUs affiliated with methanotrophic bacteria while controlling for the effect of depth. Analysis was conducted in the R environment (Rcoreteam, 2013). Parameter estimates were calculated with least squares fit between relative abundances and depth measurements. Only representative OTUs affiliated with methanotrophic bacteria were regressed. Statistical significance of correlations was determined with bootstrapping (1000 iterations), and using a Type-I error rate of 5%.

Univariate regression analyses were conducted on time-series pyrotag data to infer correlations between OTUs affiliated with methanotrophic bacteria counts and environmental variables (O₂, CH₄, NO3−, NO2−, and H₂S). Choice of model per OTU was determined by AIC testing (Akaike, 1998), using a combination of negative binomial regression (Venables and Ripley, 2002), and zero-inflated negative binomial regression (Zeileis et al., 2008) employed due to many zeroed observations (average: 54%). Parameter estimates and statistical significance were calculated in the R environment (Rcoreteam, 2013), using a Type-I error rate of 5%.

To generate a robust network emphasizing co-occurrences between prevalent OTUs in water column compartments defined by HCA rather than individual depth intervals, the Bray-Curtis and Spearman's rank correlations were used. Correlation coefficients were calculated using CoNet (Faust et al., 2012) with OTU abundance as count data. First, an OTU matrix derived from the pyrotag taxonomic analysis was transformed into presence-absence data to remove OTUs with less than 1/3 zero counts, leaving a matrix of 780 OTUs for all samples. Next, to construct ensemble networks, measure-specific thresholds set to 0.6 were used as a pre-filter and edge scores were computed only between clade pairs. To assign statistical significance to the resulting scores, edge and measure-specific permutation and bootstrap score distributions with 1000 iterations each were computed. p-values were tail-adjusted so that low p-values correspond to co-presence and high p-values to exclusion. After merging, p-values on each final edge were corrected to q-values (cut-off of 0.05). The positivity or negativity of each relationship was determined by consensus voting over all integrated data sources. Finally, only edges with at least two supporting pieces of evidence were retained.

The final edges matrix was visualized as a force directed network using Cytoscape 2.8.3 (Shannon et al., 2003). Network properties were calculated with the “Network Analysis” plug-in. Nodes in the co-occurrence network corresponded to individual OTUs and edges were defined by computed correlations between corresponding OTU pairs. The layout revealed distinct modules, which persisted after lowering the correlation coefficient cut-off for edge creation to 0.90 reinforcing the robustness of the network. Edges from modules were selected and visualized as sub-networks using the tool Hive Panel Explorer (https://github.com/hallamlab/HivePanelExplorer/wiki) (Perez, 2015). HivePlotter allows for edge selection based on interactions within specific OTUs such as those affiliated with methanotrophic bacteria.

The SSU rRNA gene sequences reported in this study have been deposited in the NCBI under the accesion numbers: GQ346856, GQ349233, GQ347199, HQ163221, GQ350623, GQ349295. The SSU rRNA pyrotag sequences reported in this study have been deposited in the NCBI under the accession numbers: SAMN03387532–SAMN03387915. Metagenomes reported in this study have been deposited in the NCBI under the accession numbers: SAMN05224436, SAMN05224437, SAMN05224442, SAMN05224443, SAMN05224447, SAMN05224451. Metatranscriptomes reported in this study have been deposited in the NCBI under the accession numbers: SAMN05238739, SAMN05238741, SAMN05238743, SAMN05238745, SAMN05238748, SAMN05238751.

We monitored changes in water column conditions corresponding to the progression of stratification (from winter through mid-summer) and deep water renewal (late summer into fall) events over a 2-year period (May 2008–July 2010). As previously described for OMZs (Wright et al., 2012), in this study we define water column compartments on the basis of O₂ concentration ranges: oxic (>90 μM O₂), dysoxic (90-20 μM O₂), suboxic (20-1 μM O₂) and anoxic (<1 μM O₂). Between May and August 2008, as water column stratification peaked, suboxic conditions intensified. This intensification corresponded with increasing levels of CH₄ and hydrogen sulfide (H₂S) below 150 m consistent with the development of deep-water anoxia. Two CH₄ concentration peaks were observed, at the subsurface (20–100 m) ranging between the 20–120 nM, and at in the deep basin increasing steadily below 150 m to a maximum of 800 nM at 200 m in late July. Over the same time interval, the concentration of H₂S ranged between 2 to 8 μM in the anoxic bottom waters (150–200 m). The concentration of NO3− between the surface and 100 m ranged between 5 to 20 μM, decreasing rapidly between 100 and 135 m before reaching a minimum of ≤1 μM in anoxic bottom waters. The concentration of NO2− between surface and 135 m ranged between <0.1 and 0.25 μM reaching a maximum 0.3 μM in anoxic bottom waters (150–200 m) (Supplementary Figure 1).

The beginning of 2008 deep-water renewal occurred between the end of July and early September, continuing through November. During this time interval, oxygenated nutrient-rich waters flowed over the sill, displacing anoxic bottom waters upwards and disrupting the redox gradient established during spring and summer months. Between September and October, dissolved O₂ was observed throughout the water column, although upwards shoaling O₂ and NO3− depleted bottom waters produced an intermediate suboxic layer between 100 and 135 m. Concomitantly, CH₄ concentrations increased transiently in the suboxic transition zone (120–135 m) ranging between 180 and 700 nM while decreasing below 135 m to 0 nM in bottom waters. In November, O₂ and NO3− and NO2− concentrations continued to increase above 100 and below 135 m with intensification of water column O₂ deficiency within intervening depth intervals (Supplementary Figure 1).

Beginning in December 2008 water column O₂ deficiency intensified below 100 m consistent with stratification. However, in contrast to previous studies at Saanich Inlet, our observations indicate a significantly weak renewal occurred in fall 2009 as indicated by no measurable increase in O₂ or NO3− in deep basin waters. This phenomenon extended O₂ deficiency below 100 m resulting in CH₄ accumulation through the summer of 2010, and corresponded with anoxic conditions below 135 m (O₂ and NO3− concentrations equal to 0 μM combined with high levels of H₂S). Over this period, the concentration of CH₄ between surface and 100 m ranged between 20 and 500 nM increasing with depth to reach a maximum 1250 nM in anoxic bottom waters in July 2010. The concentration of H₂S ranged between 2–20 μM in the anoxic waters. Over this period, the highest NO2− concentration (1.8 μM) was observed at 120 m in July 2010 (Supplementary Figure 1).

To identify microbial agents driving methane oxidation in the Saanich Inlet water column we analyzed SSU rRNA gene pyrotag sequences from 288 samples collected at 16 depths (10–200 m) over a 2 year time period between May 2008 and July 2010. Results revealed consistent microbial community partitioning as commonly observed in stratified ecosystems where O₂ deficiency is associated with redox-driven niche partitioning (Alldredge and Cohen, 1987; Shanks and Reeder, 1993; Wright et al., 2012). For instance, hierarchical cluster analysis (HCA) resolved three major groups or clusters (AU ≥ 70, 1000 iterations) associated with oxic (group I), dysoxic-suboxic (group II), and anoxic (group III) water column conditions (Figure 1A). Consistent with oxycline formation, O₂ (R² = 0.80) and NO3−(R² = 0.63) were negatively correlated with groups II and III (Figure 1B), while positively correlated with NH₄ (R² = 0.26) and H₂S (R² = 0.52) (Figure 1B).

Overall, microbial community composition was dominated (relative abundance >1%) (Rappé and Giovannoni, 2003) by OTUs affiliated with ubiquitous and abundant taxonomic groups previously identified in marine O₂ deficient environments (Field et al., 1997; Fuhrman and Davis, 1997; Brown and Donachie, 2007; Tripp et al., 2008; Lavik et al., 2009; Walsh et al., 2009; Zaikova et al., 2010; Walsh and Hallam, 2011; Wright et al., 2012) including SAR11, SAR324, Nitrospina, SUP05, Marine Group A and Methylophilales within the bacterial domain, and Thaumarchaeota within the archaeal domain (Figure 2). Consistent with previous observations in OMZs (Stevens and Ulloa, 2008; Glaubitz et al., 2010), OTUs affiliated with Methylococcaceae and Methylomonas, both canonical type I methanotrophs, were identified among the rare biosphere (relative abundance <1%) (Sogin et al., 2006). Additionally, OTUs affiliated with methylotrophic bacteria such as Methylobacteriaceae and Methylophaga, were also identified among the rare biosphere (Figure 3A).

Methanotroph diversity in Saanich Inlet was determined based on recruitment of representative OTU sequences to full-length SSU rRNA gene reference sequences (full description of reference sequences used in Methods). Initially, we identified all OTUs with a taxonomic assignment affiliated with methanotrophic bacteria using BLAST-based comparisons conducted in QIIME queried against the Silva and Greengenes reference databases. Sequences affiliated with these OTUs shared 90% identity with Methylococcales reference sequences in Silva, and 90% identity with Methylococcaceae and Methylomonas reference sequences in Greengenes. Subsequently, further analysis on representative methanotrophic OTUs was conducted using full-length SSU rRNA gene sequences as fragment recruitment platforms. A total of 3804 sequences affiliated with type I methanotrophs were recovered, representing 0.095% of the total pyrotag sequences generated. Sequences clustered at 97% identity resolved into 66 distinct OTUs, 6 of which contained 75% of total methanotroph sequences. Similar to full-length SSU sequences, pyrotag OTUs revealed 4 subgroups with phylogenetic similarity to cultured and environmental representatives of type I methanotrophs: Methylococcaceae (4%) (Mou et al., 2008), putative methanotrophic groups OPU1 (22.55%) and OPU3 (26.23%) (Hayashi et al., 2007; Tavormina et al., 2010, 2013), and methanotrophic symbionts (46.89%) (Streams et al., 1997; Dubilier et al., 2008; Petersen and Dubilier, 2009; Figure 3). The most abundant OTUs were related to putative methanotrophic OPU1 (OPU_01 = OTU55333; 16.9%), OPU3 (OPU3_01 = OTU6504; 8.4%), and methanotrophic symbionts (Symbiont_01 = OTU39693; 45.5%) (Figure 3B), comprising ~70% of total methanotroph sequences found in the pyrotag datasets. Based on this information we focused our analysis on these three groups.

Population dynamics of OPU1_01, OPU3_01 and Symbiont_01 were determined throughout the water column between May 2008 and July 2010. OPU1_01 and OPU3_01 were more abundant under dysoxic (<90 μM O₂) and suboxic (≤20 μM O₂) conditions while Symbiont_01 was more abundant under suboxic and anoxic (≤3 μM O₂) conditions (Figure 4A). Specifically, OPU1_01 showed two abundance peaks: (1) after the 2008 renewal (September–November) reaching up to 0.24% relative abundance under suboxic conditions (150 and 165 m), and (2) during water column stratification (March and April 2009) reaching up to 0.15% relative abundance under oxic conditions (10–85 m). OPU3_01 peaked during stratification periods (July–August 2008 and March–April 2009) reaching up to 0.15% relative abundance under oxic and dysoxic conditions (10–97 m), and during the extended stratification in July 2010 reaching up to 0.2% relative abundance under suboxic conditions (120–135 m, Figure 4B). Symbiont_01 also showed two abundance peaks: (1) during the 2008 stratification (May–July) reaching up to 0.23% relative abundance, and (2) after the 2008 renewal (November–January 2009) reaching up to 0.6% relative abundance under oxic conditions (40–97 m), and 0.55% under suboxic—anoxic (165–200 m) conditions (Supplementary Figure 2). Cumulative abundance for methanotroph OTUs throughout the water column indicated co-occurrence and peak abundance under dysoxic-suboxic water column conditions over time, suggesting that the highest CH₄ oxidation activity is carried out between these depth intervals in Saanich Inlet (Supplementary Figure 2 and Figure 3).

Observed OTU distribution patterns for OPU1_01, OPU3_01 and Symbiont_01 became increasingly compartmentalized during the extended stratification period. For instance, OPU1_01 was observed under dysoxic conditions (90–110 m) and OPU3_01 under dysoxic-suboxic conditions (100–135 m), while Symbiont_01 under suboxic-anoxic conditions (135–200 m) (Figure 4). These observations point to redox-driven niche partitioning among methanotrophic bacteria with the potential to mediate CH₄ oxidation in different water column compartments. Given these distribution patterns in relation to measured geochemical profiles i.e., OPU3 and NO2− (Figure 4B), we hypothesized the use of alternative terminal electron acceptors in the CH₄ oxidation process.

Given the distribution of OPU1_01, OPU3_01 and Symbiont_01 OTUs in the Saanich Inlet water column we were interested in determining functional potential and activity of these groups. This was determined based on the number of conceptually translated pmoA genes and transcripts found throughout the water column in February 2010 and the recruitment of water column PmoA sequences to selected reference sequences (see Methods). Initially, we identified all sequences with a functional assignment affiliated with PmoA using the MetaPathways functional annotation table output. Similar to observations made for methanotroph OTU abundance, PmoA sequences recovered in the metagenomic and metatranscriptomic datasets were rare, representing ≤0.001% (34 sequences) and ≤0.018% (26 sequences) from the total number of predicted proteins in the metagenomic and metatranscriptomic datasets, respectively. Representative sequences in pmoA clone libraries were mostly related to OPU3 (Hayashi et al., 2007), and methanotrophic symbionts within the type I methanotroph clade. However, all February 2010 metagenomic, and metatrasncriptomic PmoA representative sequences clustered at 97% similarity were related to OPU3 (Hayashi et al., 2007) (Figure 5). Corresponding metatransciptome RPKM values for pmoA in February 2010 showed differential expression across the redox transition zone under dysoxic-suboxic conditions. Transcript expression, depicted as dots sized based on RPKM values (Figure 5), was higher at 100 m (mean RPKM = 124, SD = 12.10) than at 120 (mean RPKM = 15.25, SD = 7.33) and 135 m (mean RPKM = 17.05, SD = 0.58) potentially indicating the boundaries for OPU-mediated CH₄ oxidation in the Saanich Inlet under water column stratification conditions (Figure 5).

To better constrain niche partitioning among methanotrophic bacteria and the potential use of alternative terminal electron acceptors, i.e., NO3− and NO2−, we conducted multivariate linear and beta regression analyses on the time-series data. Multivariate regression allowed us to minimize the possible linear effect of depth on OTU distribution while beta regression allowed us to use environmental and relative abundance data to model variable correlations (Kieschnick and Mccullough, 2003; Ferrari and Cribari-Neto, 2004). Multivariate analysis showed significant positive correlation (p < 0.001) between OPU1_01 and OPU3_03 over stratification periods. However, negative correlation (p < 0.01) between OPU1_01 and symbiont_01 was also observed during these time periods. In addition, significant positive correlation (p < 0.0001) between OPU3_01 and Symbiont_01 was also observed during transitioning to stratification periods (Supplementary Figure 2) (Table 1). Negative binomial regression results on the time-series data indicated OPU1_01 distribution was significantly and negatively correlated with O₂ (p < 0.05) and CH₄ (p < 0.001) and weakly correlated with NO3−, NO2− and H₂S. OPU3_01 distribution was significantly and negatively correlated with CH₄ (p < 0.001) and NO2− (p < 0.05), and weakly correlated with NO3− and H₂S (Table 2, Figure 4). Symboint_01 distribution was weakly significantly and negatively correlated with CH₄, H₂S, NO3−, and NO2− (p < 0.001), and weakly correlated with O₂ (Table 2).

To identify characteristic OTUs occurring under specific water column oxygen conditions, multi-level indicator species analysis (ISA) was conducted based on groups resolved in HCA. A co-occurrence network was then constructed to identify potential interactions with methanotrophic OTUs. Multi-level ISA identified OTUs affiliated with one carbon (C1) utilizing microorganisms including Methylophilales, Methylophaga, SAR324, Verrucomicrobia and Planctomycetes. SAR324 and Verrucomicrobia were indicators for oxic and anoxic water column conditions while Methylophilales and Methylophaga were indicators for suboxic and anoxic water column conditions. The indicator OTUs affiliated with Planctomycetes were evenly distributed throughout the water column (Supplementary Tables 1–3). OPU1_01 was identified as an indicator species for dysoxic-suboxic water column conditions, while no OTUs affiliated with OPU3_01, or Symbiont_01 were identified as indicator species for a particular water column condition. However, OPU3_01 and Symbiont_01, as well as OTU20751, affiliated with the OPU3 group, were identified as indicators for combined oxic and dysoxic-suboxic conditions (Supplementary Table 2).

Co-occurrence analysis based on Bray-Curtis and Spearman correlation values among OTUs, resulted in a microbial network (Figure 6A) partitioned in modules similar to hierarchical clustered groups associated with oxic (I), dysoxic-suboxic (II and III), and anoxic (IV) water column compartments (Supplementary Figure 1). Interestingly, significant positive correlations (CV > 0.6, p < 0.001), shown as sub-networks (Supplementary Figure 1, Table 5), were observed among OPU1_01, OPU3_01 and Symbiont_01 with indicator OTUs affiliated with potential C1 utilizing microorganisms such as Methylophaga, Methylophilales, SAR324, Verrucomicrobia and Planctomycetes, and other ubiquitous OMZ microbes including representative taxa such as Marine Group A, Nitrospina, and SUP05 (Figure 6B, Supplementary Table 4). To provide further evidence for potential methanotroph interactions observed in sub-networks we compared microbial OTUs to SSU rRNA gene sequences recovered from previous CH₄ microcosm experiments using Saanich Inlet suboxic waters (Sauter et al., 2012) (Supplementary Figure 2). Most bacterial OTUs (83% of OTUs in sub-networks) were found affiliated (≥80% identity) to microcosm SSU rRNA sequences that were enriched and active after CH₄ addition, and related to Bacteroidetes, Marine Group A, Planctomycetes, Sphingomonadales, Nitrospina, Methylophilales, Methylophaga, SUP05 and Verrucomicrobia (Supplementary Table 5) reinforcing co-occurrence network observations.

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. The reviewer MD and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.