Characterization of microbial associations with methanotrophic archaea and sulfate-reducing bacteria through statistical comparison of nested Magneto-FISH enrichments

Sample collection and Magneto-FISH

iTag Magneto-FISH enrichments were conducted using a large scale (1 L) incubation of methane seep sediment from Hydrate Ridge North (offshore Oregon, USA) collected in September 2011 at 44°40.0$2^{\prime }\mathrm{N}$ 125°6.0$0^{\prime }\mathrm{W}$, from a water depth of 775 m using the ROV JASON II and the R/V Atlantis. Marine sediment was collected using a push core to sample a sulfide-oxidizing microbial mat adjacent to an actively bubbling methane vent. A sediment slurry from the upper 0–15 cm depth horizon of the push core was prepared with 1 volume N₂ sparged artificial seawater to 1 volume sediment, overpressurized with methane (3 bar) and incubated at 8 °C in a 1 L Pyrex bottle capped with a butyl rubber stopper until subsampling for Magneto-FISH.

In February 2015, incubation samples were immediately fixed in 0.5 ml sediment aliquots in 2% paraformaldehyde (PFA) for 3 h at 4 °C. The samples were washed in 50% phosphate-buffered saline (PBS): 50% EtOH, then 75% EtOH: 25% DI water, and resuspended in 2 volumes (1 ml) 100% ethanol. Samples were centrifuged at 1,000 × g for 1 min between wash steps. After fixation, the Magneto-FISH method first described by Pernthaler et al. (2008) and further optimized by Schattenhofer & Wendeberg (2011) and Trembath-Reichert, Green-Saxena & Orphan (2013) was used. Briefly, a liquid CARD-FISH reaction was followed by immunomagnetic bead incubation coupled with anti-fluorescein attaching magnetic beads to CARD-FISH hybridized aggregates. Samples were then held against magnets and the sediment matrix was washed away, retaining target cells and physically associated microbes in the magnetic portion as described in Trembath-Reichert, Green-Saxena & Orphan (2013). Four previously published FISH probes were used targeting a range of Deltaproteobacteria and Methanomicrobia (Table 1). A subset of three 0.5 ml aliquots was also immediately frozen before fixation (unfixed bulk sediment), and another three aliquots were frozen after fixation (fixed bulk sediment) for bulk sediment comparison with Magneto-FISH enrichments. Sediment for MSMX-Eel_932 Magneto-FISH metabolic gene analysis was fixed and washed onboard in September 2011, as described above. See methods flow chart provided in Fig. S1.

iTag amplification

For iTag sequencing, ten Magneto-FISH enrichments were performed in parallel using the FISH probes DSS_658 (triplicate), MSMX-Eel_932 (triplicate), SEEP-1a_1441 (duplicate), Delta_495a + Delta_495a competitor (duplicate). Magneto-FISH enrichments and bulk sediment samples were resuspended in 650 µl solution PM1 and transferred to silica tubes from the PowerMicrobiome RNA Isolation Kit (MoBio). This kit was chosen based on manufacturer recommendation for formalin-fixed sediment samples, with the added capability to co-elute RNA if desired. 6.5 µl of beta-mercaptethanol was added, and samples were mechanically lysed in a bead beater (FastPrepFP120; Thermo Electron Corp., Waltham, MA, USA) for 45 s at setting 5.5 and incubated at 65 °C for 3.5 h. The remaining steps in the PowerMicrobiome RNA Isolation Kit were followed according to manufacturer instructions (starting at step 5) without any DNase procedures, and eluting in a final volume of 60 µl ultrapure water. DNA extracts were quantified using a Qubit Flurometer and HS dsDNA kit (Invitrogen; Table S1). All but one Magneto-FISH sample had DNA concentrations below detection (<0.5 ng/µl); however, all samples yielded PCR amplicons when viewed on a gel after initial pre-barcoding PCR (30 cycles).

iTag samples were prepared with Earth Microbiome Project (EMP) primers 515f and 806r (Caporaso et al., 2012). An initial amplification of 30 cycles with primers lacking the barcode, linker, pad, and adapter was performed for all samples, in duplicate. Duplicate PCR reactions were pooled and reconditioned for 5 cycles with barcoded primers, for a total of 35 cycles. A master mix of 2 X Q5 Hot Start High Fidelity Master Mix (NEB) and 10 µM forward and reverse primers was prepared for a final volume of 15 µl per sample, with 1 µl DNA template. PCRs had an initial 2 min heating step at 98 °C, followed by cycles of 10 s 98 °C, 20 s 54 °C, and 20 s 72 °C, and finished with a final extension of 2 min at 72 °C. PCR negative controls, substituting ultrapure water for DNA template, were amplified for 40 cycles total. We note that these are not the official recommended reagents or PCR conditions from the EMP, but internal lab tests showed that for 6 out of 9 mock community taxa, recovered sequence relative abundances were more accurate when using Q5 polymerase rather than the recommended Hot Start MasterMix (5-prime). EMP primers were chosen for iTag for cross-comparison between studies, though there is known primer bias within this universal primer set (Parada et al., 2015) and sequencing reactions will always have some inherent variability.

Mock communities

Four mock communities were prepared with a range of relative proportions of nine common methane seep taxa (Table S2). Full-length 16S rRNA gene plasmids from each taxa listed were quantified by Qubit. Taking into account the plasmid’s nucleotide composition and length in order to calculate its molecular weight, plasmids were quantitatively combined in known volumetric fractions to achieve a range of desired mock community compositions. These combined plasmid mixes were diluted to ∼1 ng/µL and then prepared according to the same iTag methods as all other samples.

iTag sequence processing

We followed the mothur Standard Operating Procedure (SOP) for Illumina MiSeq sequencing of the 16S rRNA gene V4 region, accessed May 2015 and using methods described in Kozich et al. (2013) with UCHIME chimera checking (Edgar et al., 2011). A concatenated file of the mothur version of separate archaeal and bacterial SILVA 119 databases (Quast et al., 2013) was used for alignment and classification. Unfixed Bulk Sediment 1 only returned 8% of the average DNA concentration of the other two samples. (Table S1). This sample was removed from statistical analyses because it fails to be a representative of the unfixed bulk sediment community baseline. The mock communities were processed following the “Assessing Error Rates” section of the mothur SOP to compute sequencing error rates and spurious OTU rates (Table S4). Additional analysis demonstrating sequence processing did not selectively remove ANME-2c sequences and relative sequence abundances recovered with iTag sequencing of mock communities are provided in Tables S3 and S2, respectively.

Using R version 3.1.3 (R Core Team, 2015), an average number of sequences per OTU was calculated from unfixed bulk sediment samples (2 and 3). All OTUs with an average relative sequence abundance below 0.1% in the unfixed bulk sediment were identified and removed from all samples using mothur. 135 unique OTUs remained out of 25,354. We also verified that after the 0.1% cutoff was applied, no negative control contaminant OTUs remained. The top 20 OTUs amplified from the no template negative control were classified as, in order of sequence abundance: Sphingomonas*; Planctomyces*; Escherichia-Shigella*; Staphylococcus; Roseomonas*; Pir4_ lineage; Delftia*; Macrococcus; Myxococcales;0319-6G20;unclassified; Planctomyces; Enhydrobacter; Sphingobium*; Caenispirillum; Bacillus*; Pseudoxanthomonas*; Peptoniphilus; Lysobacter; Salinicoccus; Propionibacterium.* Reagent contaminant genera discussed in Salter et al. (2014) are denoted by (*). All samples (including mock community and negative controls) were submitted to the SRA under the accession SAMN03879962, BioSample: SAMN03879962, Sample name: PC47 (5133-5137) mixed slurry.

Gene libraries of the Magneto-FISH samples were prepared as in Trembath-Reichert, Green-Saxena & Orphan (2013) using the primers and annealing temperatures listed in Table 1 and TOPO TA Cloning Kit for Sequencing with pCR4-TOPO Vector and One Shot Top 10 chemically competent Escherichia coli (Life Technologies). All full-length 16S rRNA gene sequences were aligned by the SINA online aligner (v1.2.11) (Pruesse, Peplies & Glöckner, 2012) and added using maximum parsimony to the SILVA 119 database (Quast et al., 2013) for classification. A taxonomy-based count table was prepared (sequences per taxa, per sample) and all taxa absent from the bulk sediment library were removed from Magneto-FISH enrichment libraries (for parity with iTag contaminant removal processing). Functional gene sequences were translated using the EMBOSS online translation tool (Li et al., 2015), then added to ARB (Ludwig et al., 2004) databases for phylogenetic placement and classification. Sequences were submitted to NCBI under the following accession numbers: AprA ( KT280505– KT280517), DsrA ( KT280518– KT280533), McrA ( KT280534– KT280581), Archaeal 16S rRNA gene ( KT280582– KT280632), Bacterial 16S rRNA gene ( KT280633– KT280909), SoxB ( KT280910– KT280928). Gene trees were computed with representative sequences using PhyML 3.0 (Guindon et al., 2010) online execution with defaults on the South of France Bioinformatics platform.

Statistical analysis

Weighted UniFrac (Lozupone & Knight, 2005), Metastats (White, Nagarajan & Pop, 2009), and linear discriminant analysis (LDA) effect size (LEfSe) (Segata et al., 2011) analyses were computed in mothur as outlined in the mothur SOP. Co-occurrence statistical analyses were run using the table of 135 unique OTUs in the format of sequence counts of each OTU per sample. The program SparCC was used to determine significant correlations (Friedman & Alm, 2012). This analysis was run 100 times with default settings, except 10 iterations were used instead of 20. OTUs with SparCC correlations above an absolute value of 0.6 with p-values below 0.01 were considered significant. Resulting associations that occurred in at least 50 out of 100 network iterations are provided in Table S5. Cytoscape (Shannon et al., 2003) was used to display associations in Fig. 1.

CARD-FISH microscopy

A triple CARD-FISH hybridization was performed with bacterial probes listed in Table 1, ANME-1_350 and MSMX-Eel_932. The sample preparation and CARD reaction was performed as per Green-Saxena et al. (2014). After the three CARD reactions, samples were post-stained with DAPI (25 ng/µl). CARD signal within any part of a physically attached group of cells larger than 10 µm was counted as a positive identification. For example, a large EPS matrix that contained many smaller separate ANME-1 and ANME-2 aggregates would count as one positive identification for each clade. This was done to simulate groups that would have been isolated together in a Magneto-FISH enrichment. Since the MSMX-Eel_932 probe also targets the ANME-1 population, only cells with MSMX-Eel_932 signal and no ANME-1_350 signal were recorded as an ANME-2 positive identification to comprehensively target ANME-1, -2, and a bacterial partner in a triple CARD-FISH hybridization set. ANME-3 were not recovered in the iTag dataset and were not considered as potential contributors to MSMX-Eel_932 signal.

Relative sequence abundance of seep microbiome taxa in 16S rRNA gene iTag and libraries

Relative sequence abundances of the methane seep microbiome characteristic taxa, ANME archaea, Deltaproteobacteria, Hyd24-12, and Atribacteria (Ruff et al., 2015), were compared two ways: (1) between iTag and gene library 16S rRNA gene samples to determine how relative sequence abundances differed between sequencing methodologies, and (2) between Magneto-FISH enrichment and bulk sediment to determine taxa-specific relative sequence abundance for each probe (Table 1).

Mock community analysis showed that ANME-2 were always underrepresented in iTag data (0.32–0.81 fold of what was expected), whereas the Deltaproteobacteria and ANME-1b were more faithfully represented (Table S2). ANME-1a was consistently overamplified. By normalizing the relative sequence abundance of ANME-2c, -2a/b, and -1a to the abundance of ANME-1b, the most faithfully amplified archaea in the mock community data (Table S2), we could compute a ratio between the average relative sequence abundance in fixed bulk sediment samples between iTag and the archaeal 16S rRNA gene library. ANME-2c (0.04 iTag:clone ratio), ANME-2a/b (0.12), and ANME-1a (0.40) were all less abundant in iTag sequences as compared to the archaeal gene library (calculated from values in Table 2). Similarly comparing SEEP-SRB1 to Desulfobulbus between the two methods in fixed bulk sediment returns a ratio of 0.41 iTag:clone. Since the iTag methodology recovers far more diversity (e.g., Desulfobacula, Desulfocapsa, Desulfoluna, Atribacteria, and Hyd24-12 were not recovered in the bacterial 16S rRNA gene bulk sediment library), it is expected that the relative sequence abundances of each individual taxon computed from iTag data would be less than the domain targeted 16S rRNA gene libraries. However the ANME-2c abundance ratio was an order of magnitude less than ANME-1a and SEEP-SRB1 ratios, and appears to be an extreme case of underestimation in iTag data. There was also variation between Magneto-FISH enrichment replicates, as indicated by the high standard deviations of Magneto-FISH samples as compared to bulk sediment samples. The degree of variation (average standard deviation across all taxa listed) correlated with the specificity of the probe; where Delta_495a had the lowest average standard deviation and Seep-1a_1441 had the highest average standard deviation.

The high relative sequence abundance taxa (>1.5 fold relative sequence abundance increase over fixed bulk sediment; Table 2) in the averaged Seep-1a_1441 iTag Magneto-FISH enrichments were Desulfoluna (2.20), SEEP-SRB1 (2.36), Hyd24-12 (3.44) and Atribacteria (1.51) (Table 2). The DSS_658 enrichment had fewer high relative sequence abundance taxa with only Desulfoluna (4.62), Spirochaeta (4.36), and Atribacteria (4.80). The Delta_495a enrichment also had three high relative sequence abundance taxa with Desulfobulbus (2.52), Spirochaetae-uncultured (3.70), and Atribacteria (3.02). The MSMX-Eel_932 enrichment had six high relative sequence abundance taxa with Desulfococcus (1.85), Desulfoluna (8.47), SEEP-SRB1 (1.67), Spirochaeta (1.63), Hyd24-12 (1.73), and Atribacteria (7.18). Gene library results showed high relative sequence abundance (>1.5) in both ANME and Deltaproteobacteria with DSS_658 and MSMX-Eel_932 enrichments (Table 2). Similar to the bulk sediment, Desulfobacula, Desulfocapsa, Desulfoluna, Atribacteria and Hyd24-12 were not recovered in the bacterial 16S rRNA gene Magneto-FISH libraries. MSMX-Eel_932 enriched for SEEP-SRB1 (2.73), SEEP-SRB4 (3.28), Desulfococcus (3.82), Spirochaeta (1.64), and ANME-2a/b (2.51) in 16S rRNA gene libraries. There was also a slight enrichment of ANME-2c (1.28). The DSS_658 enrichment had high relative sequence abundance for SEEP-SRB1 (1.74), SEEP-SRB2 (2.78), ANME-2c (1.54), and ANME-2a/b (2.24) with iTag, but these same taxa did not have high relative sequence abundance in the gene library. Spirochaeta and SEEP-SRB1 had high relative sequence abundance in both iTag and gene libraries for MSMX-Eel_932 enrichments. Relative sequence abundances for all non-core methane seep taxa in iTag samples are included in Table 3, and where Magneto-FISH enrichments of these additional taxa support network co-occurrences they are discussed in network results.

Statistical evaluation of Magneto-FISH enrichment

To statistically compare enrichment microbial communities, we used a suite of statistical tests including: non-parametric T-tests (White, Nagarajan & Pop, 2009), LEfSe (Segata et al., 2011), and UniFrac (Lozupone & Knight, 2005). Using the T-test comparison, ten OTUs were significantly (p < 0.001) different between the bulk sediment and Magneto-FISH samples (when only including OTUs with sequences present in both groups). The taxonomic assignments for these ten OTUs were: WCHB1-69, Desulfobulbus, Thaumarcheota, ANME-1a, Bacteroidetes (VC2.1), ANME-2c, Caldithrix, SEEP-SRB1, Candidate Division TA06, and Gammaproteobacteria (CS-B046). LEfSe was then used to determine which OTUs were significantly different between Magneto-FISH enrichments and bulk sediment. We found three OTUs were significantly (p-value < 0.05) higher in relative sequence abundance in Magneto-FISH samples over bulk sediment with the taxonomies: SEEP-SRB1, Desulfobulbus, and Planctomycetes (SHA-43).

Weighted UniFrac analysis was used to compare the community composition between Magneto-FISH iTag enrichments. The UniFrac metric represents the fraction of the branch length that is unique to each sample, or unshared between samples, such that a higher ratio means less similar samples. The Deltaproteobacteria probe enrichment communities were more similar to each other than any of the Deltaproteobacteria probes compared with the MSMX-Eel_932 probe (Table 4). The most distinct communities were MSMX-Eel_932 enrichment and Delta_495a enrichment, with the highest proportion of unshared branch length (0.97; p-value < 0.001). MSMX-Eel_932 enrichment and DSS_658 enrichment had less unshared branch length at 0.88 (<0.001), suggesting MSMX-Eel_932 and DSS_658 probes enrich for a more similar community than MSMX-Eel_932 and Delta_495a probes. Comparison of the MSMX-Eel_932 enrichment and SEEP-1a_ 1441 enrichment communities was not significant at the <0.001 cutoff. Within the Deltaproteobacteria probes, SEEP-1a_1441 enrichment and DSS_658 enrichment had the lowest proportion of unshared community (0.77, <0.001); the most similar community structures were recovered with these two probes. The next lowest proportion of unshared community is between DSS_658 enrichment and Delta_495a enrichment (0.81). SEEP-1a_1441 enrichment and Delta_495a enrichment are least similar, at 0.85. All of these values are highly significant (<0.001). This is consistent with the expectation that the overlap between the target microbial population of the SEEP-1a_1441 probe would be most similar to the target microbial population of the DSS_658 probe, while the Delta_495a enrichment would recover more total Deltaproteobacteria diversity.

Assessing community structure with co-occurrence network analysis

After determination of statistically significant differences between iTag Magneto-FISH and bulk sediment samples, we computed co-occurrence networks to observe which of the 135 most abundant OTUs were correlated in the methane seep microbial community. By combining the results from 100 separate microbial association calculations, we were able to assign confidence to each microbial association and determine the most robust associations. Significant associations are reported in Table S5 and depicted as a network in Fig. 1.

Focusing first on the common ANME syntrophic Deltaproteobacteria partner, SEEP-SRB1, this taxon had the most associations in the network including nine positive associations and one negative association (Fig. 1). There are two separate sets of SEEP-SRB1 & Planctomycetes (AKAU3564 sediment group) positive associations that are both well supported. SEEP-SRB1 is also associated with three other heterotrophic taxa (Candidate Phylum Atribacteria, Spirochaeta, and Bacteroidetes (VC2.1_ Bac22)) and one sulfur-oxidizing taxa (Sulfurovum). SEEP-SRB1 was also associated with Candidate Division Hyd24-12, which has a currently unknown ecophysiology, but could be a heterotroph if the topology of heterotrophic taxa being in the center of the network holds true. Hyd24-12 and Atribacteria are also both associated with the second most associated taxa, Candidate Division OD1, but there was no direct association between SEEP-SRB1 and OD1. SEEP-SRB2 has two of the same associations as SEEP-SRB1 (VC2.1_Bac22 and Atribacteria), but is the only Deltaproteobacteria associated with MBG-B, Anaerolineaceae, and Desulfoluna (another Deltaproteobacteria). SEEP-SRB4 is associated with Desulfobulbus, and the only Deltaproteobacteria associated with and ANME (2a/b), WS3, and Actibacter. WS3 had high relative sequence abundance in both DSS_658 and MSMX-Eel_932 enrichments (Table 3). Desulfobulbus is associated with Desulfococcus, the only Deltaproteobacteria associated with BD2-2, and SAR406. SAR406 had high relative sequence abundance in Seep1a_1441 and Delta_495a enrichments (Table 3). The heterotroph Spirochaeta is also included in the core methane seep microbiome and was associated with Clostridia and WS3, in addition to Hyd24-12 and SEEP-SRB1.

In examination of additional OTUs associated with sulfur metabolisms, we found Sulfurovum and Sulfurimonas (Epsilonproteobacteria) were not associated with each other, but are both associated with Deltaproteobacteria. Sulfurimonas is associated with Desulfocapsa and Sulfurovum is associated with SEEP-SRB1 and Desulfobulbus. Sulfurovum had high relative sequence abundance in MSMX-Eel_932 enrichments and Sulfurimonas had high relative sequence abundance in Seep-1a_1441, DSS_658, and Delta_495a enrichments (Table 3). The Gammaproteobacteria, Thiohalobacter, is only associated with Anaerolineaceae and was not elevated in any of the Magneto-FISH enrichments.

Heterotrophs are the most dominant metabolic guild in the network, and similar to sulfate-reducers, have some of the most connected taxa. The heterotroph OD1 has seven positive correlations, in addition to Atribacteria and Hyd24-12 listed above: Bacteroidetes (BD2-2), Actinobacteria (OM1), Pelobacter, ANME-1b, Chloroflexi (Anaerolineaceae), and Desulfocapsa. Anaerolineaceae and Bacteroidetes (BD2-2) both had seven associations, but with different connectivity. BD2-2 was interconnected with other heterotrophs, sulfate-reducers, and archaeal methanotrophs in the main portion of the network, whereas Anaerolineaceae was connected to three taxa that share no other connections (two heterotrophs and one Gammaproteobacteria sulfur oxidizer). The one other ANME taxa in the network, ANME-1b, is only positively associated with heterotrophs and no known sulfate reducing groups.

Assessing ANME-bacterial partnerships by CARD-FISH

To assess ANME and DSS relative cell abundance, 100 aggregate clusters from the same sediment incubation (see ‘Materials & Methods’) were analyzed with CARD-FISH and the DSS_658/ANME1-350/MSMX-Eel_932 probe combination. Epsi_404, Gam_42a, SEEP-1a_1441, and CF_319A/B probes were also used with the archaeal probe combination to examine non-DSS bacterial diversity recovered in the network analysis ANME associations. All probes, target populations, and references are listed in Table 1.

30% of aggregates contained an ANME-2 signal (see ‘Materials & Methods’; Table 5) and 39% of aggregates had an ANME-1 signal. ANME-1 and ANME-2 identified cells were also consistent with expected morphologies. Multiple clusters of mixed-type ANME/DSS, DSS-only, ANME-only, DSS/non-ANME, and non-DSS/non-ANME aggregates were observed with the ANME-1_350, MSMX-Eel_932, and DSS_658 probe combination (Fig. 2A). There were no clear examples of aggregates with ANME/non-DSS hybridized cells, though we found many instances where both ANME and non-DSS cells were as part of a larger aggregate cluster with other cell types. ANME-1 cells often occurred in the matrix surrounding tightly clustered ANME-2 aggregates. The SEEP-1a_1441 probe, targeting a subgroup of DSS, was observed to hybridize with aggregate clusters that contained ANME-1 and ANME-2 cells, but usually with SEEP-SRB1/ANME-2 in tight association and ANME-1 in more peripheral association. Five of the SEEP-SRB1/ANME-2 aggregate clusters did not have ANME-1 cells (10%) and three of the SEEP-SRB1/ANME-1 aggregate clusters did not have ANME-2 cells.

Ten percent of aggregates (n = 50 counted) hybridized with the Epsi_404 probe, broadly targeting members of the Epsilonproteobacteria. These Epsilonproteobacteria were mostly found in association with other bacteria and occasionally, loosely associated with some ANME. Epsi_404 hybridized cells were generally ovoid and scattered throughout an EPS matrix of cells, as depicted in Fig. 2D. There was no apparent preference for Epsilonproteobacteria association with ANME-1 or ANME-2 aggregate clusters (Table 4). A higher percentage of aggregates had Gammaproteobacteria cells (24% of 50) than Epsilonproteobacteria cells, and there was a slightly higher co-occurrence with ANME-2 (18%) than ANME-1 (12%) hybridized cells. The dominant Gammaproteobacteria morphology observed was a cluster or chain of large (∼1 µm) ovoid cells. Gam_42a hybridizing cell clusters and chains were found both separately and associated with other bacteria, as in Fig. 2B, where they are predominately an unidentified cluster stained by DAPI with a sub-aggregate of ANME-2 cells. CF319A and CF319B were used to target Cytophaga, Bacteroidetes, Flavobacterium, and Sphingobacterium. Eight percent (n = 50 counted) of aggregates contained cells positively hybridizing with the CFB probe, generally observed as clustered filaments or rods (Fig. 2C). Half of these aggregates also had ANME-2 hybridized cells. No CFB cells were observed to co-associate with ANME-1.

Evaluation of Magneto-FISH with iTag

Challenges accompanying downstream analysis of Magneto-FISH enrichments are primarily associated with low DNA yield and poor DNA quality from aldehyde fixation (for further discussion of fixation effects see Trembath-Reichert, Green-Saxena & Orphan, 2013). Low template concentration exacerbates amplification of contaminating sequences since target and non-target templates can approach parity in a PCR reaction. Low template concentration has also been shown to create random variation in amplification products in dilution experiments (Chandler, Fredrickson & Brockman, 1997), which could explain the high variation seen in Magneto-FISH enrichment relative sequence abundances compared to bulk sediment samples. Despite these challenges, the DNA recovered from Magneto-FISH enrichments has been shown to increase the sequence abundance of target organisms relative to the bulk sediment by 16S rRNA gene sequencing and metagenomics on various Next Generation sequencing platforms (Pernthaler et al., 2008; Trembath-Reichert, Green-Saxena & Orphan, 2013). In this study, conventional cloning and sequencing of full-length bacterial and archaeal 16S rRNA genes had fewer contamination issues as compared to iTag sequencing with universal primers. Our Magneto-FISH experiments were designed to mitigate as many sampling and iTag sequencing biases between samples as possible, by concurrently extracting, amplifying, and sequencing all Magneto-FISH samples in parallel, including biological and technical replicates. The relative ratio of contaminant reads to environmental OTU’s were higher in Magneto-FISH enrichments than in bulk sediment samples, but bulk sediment could be used to separate indigenous community members from putative contaminants in the Magneto-FISH samples (see ‘Materials & Methods’). This provided a conservative Magneto-FISH dataset for statistical analyses and demonstrated the importance of parallel processing sequencing of bulk and separated samples.

In addition to issues with contaminating sequences, we also observed bias against some core methane seep microbiome taxa, where these taxa were consistently underrepresented by iTag when compared to gene libraries and CARD-FISH. ANME-2 was the most underrepresented taxon in iTag sequencing of the bulk sediment and mock communities, with much greater relative sequence and relative cell abundance in gene library sequencing and CARD-FISH analysis, respectively. It is most likely that iTag sequencing bias with the EMP primer set is the reason ANME-2c was not enriched in the Magneto-FISH samples and absent from microbial community network analysis. Members of the ANME-2a/b were also, to a lesser extent, underrepresented with iTag. In addition to our gene libraries and CARD-FISH analysis, independent assays using FISH with mono labeled oligonucleotide probes from this sediment incubation further confirmed the abundance of ANME-2 aggregates; 25% of aggregates were ANME-2c and 17% of aggregates were ANME-2b, with about half of ANME-2 aggregates associating with a bacterial partner other than SEEP-SRB1 (Supplement McGlynn et al., 2015). We conclude that while expected ANME-2 associations were not recovered, they can be explained by EMP iTag bias and therefore do not reduce the validity of other non-ANME-2 associations recovered in the co-occurrence analysis (see Tables S2 and S3 captions for further discussion of ANME-2c bias). Although ANME-1a was not underrepresented in the iTag data, it still does not appear in the co-occurrence network. In other co-occurrence network studies dominant OTUs were not associated with the majority of the microbial community, which was thought to be due to a high degree of functional redundancy (Mu & Moreau, 2015). Possible functional redundancy with other archaeal groups, or simply non-specific, loose spatial association with many taxa, as suggested by CARD-FISH analysis, could explain why ANME-1a was not recovered in our network analysis.

Despite this unanticipated methodological bias, iTag sequencing is a valid and valuable tool when combined with Magneto-FISH enrichment techniques for microbial association hypothesis development and testing. For example, we saw more bacterial OTUs, especially among Deltaproteobacteria, in the iTag samples compared with conventional gene libraries and the core methane seep taxon Hyd24-12 was not even observed among gene library sequences.

Magneto-FISH enrichment

This study provides a novel combination of nested Magneto-FISH enrichments and microbial community network analysis methods to develop hypotheses regarding specific lineage associations and, by inference, discusses the potential for additional metabolic interactions relating to sulfur cycling in methane seep sediments. Notwithstanding the low recovery of ANME-2 OTUs, there was statistical support for Magneto-FISH enrichments increasing the relative iTag sequence abundance of target organisms. Statistical analyses demonstrated SEEP-SRB1 and Desulfobulbus OTUs were significantly different in Magneto-FISH samples (t-tests), and these OTUs were significantly more enriched in Magneto-FISH samples using linear discriminant analysis (LDA) effect size (LEfSe). Additionally, weighted UniFrac analysis showed the highest percentage of shared phylogeny was between the clade-specific SEEP-1a_1441 probe and the family-specific Desulfobacteraceae DSS_658 probe enrichments. Therefore these Magneto-FISH samples contain microbial community overlap consistent with probe target specificity, even when some dominant community members are not represented at expected relative sequence abundance in the iTag analysis (ANME-2).

Magneto-FISH enrichment relative sequence abundance followed expected trends for Deltaproteobacteria (Table 2). SEEP-SRB1 had the highest relative sequence abundance in Seep-1a_1441 and MSMX-Eel_932 enrichments, which should target this group. Desulfobulbus had the highest relative sequence abundance in the Delta_495a enrichment, which was the only Magneto-FISH probe that should hybridize to this group (though Desulfobulbus could also be retrieved via association with other target organisms). OTUs affiliated with Desulfoluna (within the Desulfobacteraceae) had the highest relative sequence abundance of all Deltaproteobacteria in the DSS_658 enrichment and are also targeted by the DSS_658 probe. Desulfoluna were not specifically targeted by MSMX-Eel_932 or Seep-1a_1441 probes, but had high relative sequence abundane in these samples and may have a potential association with ANME/DSS consortia. Also, Atribacteria (JS1) was recovered in all iTag sequencing of Magneto-FISH enrichments, suggesting they may associate with either DSS/ANME or DSB/ANME consortia. Members of the Hyd24-12 were only recovered in Seep1a_1441 and MSMX-Eel_932 enrichments and may preferentially associate with SEEP-SRB1a/ANME consortia.

Evaluating our iTag relative sequence abundance data with co-occurrence analysis, we developed hypotheses that were not subject to the variation between Magneto-FISH enrichment replicates; associated taxa should always co-vary, even when they are less abundant than expected. Within the core methane seep taxa, high relative sequence abundances of Atribacteria and Hyd24-12 with SEEP-SRB1 targeting Magneto-FISH enrichments were upheld by the network. Hyd24-12 is highly associated with SEEP-SRB1, whereas Atribacteria is highly associated with both SEEP-SRB1 (DSS) and SEEP-SRB2 (DSB). While Atribacteria have not been cultured, metagenomic sequencing suggests they are likely heterotrophic anaerobes involved in fermentation (Nobu et al., 2015). Hyd24-12 was first cloned from Hydrate Ridge (Knittel et al., 2003) and has been cited as a core methane seep microbial taxon (Ruff et al., 2015), but nothing is known about its physiology. The Hyd24-12/SEEP-SRB1 association was also one of the four unique associations that were recovered in all the network computations (n = 100). These results may aid in determining a role for these enigmatic candidate phyla of the methane seep microbiome.

Methanomicrobia and Deltaproteobacteria only had one co-occurrence in our network. The one statistically supported network ANME/SRB association was between ANME-2a/b and SEEP-SRB4. SEEP-SRB4, belonging to the Desulfobulbaceae (Knittel et al., 2003), and ANME-2a/b both had high relative sequence abundance in the ANME-targeting MSMX-Eel_932 enrichment bacterial 16S rRNA gene library. There have been FISH-confirmed physical associations between ANME-2/ANME-3 and Desulfobulbaceae (Green-Saxena et al., 2014; Löesekann et al., 2007; Pernthaler et al., 2008) in AOM systems. SEEP-SRB4 was also strongly associated with the candidate phyla WS3 in the network, and WS3 was enriched in both DSS_658 and MSMX-Eel_932 enrichments. Both SEEP-SRB4 associations with ANME-2a/b and WS3 warrant future study.

While expected ANME-2/Deltaproteobacteria associations were not recovered (see Evaluation of Magneto-FISH with iTag), network analysis did recover many Deltaprotobacteria co-occurring with bacterial groups. Almost half of all positive associations contained a Deltaproteobacteria OTU (30/61), suggesting a dominant role for the sulfur cycle metabolisms. Of those, 21 associations were with a non-Proteobacteria OTU including a number of candidate organisms as described above. The association between SEEP-SRB1 and ‘AKAU3564,’ a Planctomycetes-affiliated heterotrophic sediment group, was observed twice with two separate OTU associations in this clade that were both strongly supported (occurring 100/100 and 93/100 times, respectively, that the network analysis was run, Table S5). This Planctomycete group was first described in methane hydrate bearing deep marine sediments of the Peru Margin (Inagaki et al., 2006). Planctomycetes-associated sequences were previously recovered in association with ANME-2c Magneto-FISH samples from the Eel River Basin, where the preferred partner was observed to be the SEEP-SRB1 group (Pernthaler et al., 2008). It follows that SEEP-SRB1 may also co-occur with Planctomycetes, if these organisms are affiliated (either directly or indirectly) with ANME-2 consortia. By similar logic, although it did not have high relative sequence abundance in the Seep1a_1441 enrichment, this could explain the high relative sequence abundance of this group in the MSMX-Eel_932 enrichment (Table 3). Planctomycetes targeted CARD-FISH hybridization using the general Planctomycetes probe Pla_886 was attempted; however, many cells with a morphology similar to ANME-1 were hybridized and the results were deemed inconclusive. This ambiguity could be due to the probe’s single base pair mismatch to 97% of ANME-1a, 94% of ANME-1b, and 25% of ANME-2b, even if this mismatch was centrally located (SILVA TestProbe online tool, Greuter et al., 2016). Spirochaeta was also associated with SEEP-SRB1, in addition to Hyd24-12 and WS3, and had high relative sequence abundance in both the DSS_658 and MSMX-Eel_932 enrichments (Table 2). In addition to being core methane seep microbial taxa, some members of the Spirochaetes have sulfide-oxidizing capabilities in mats with sulfidogenic bacteria (Dubinina, Grabovich & Chernyshova, 2004) and its possible that these organisms may be utilizing sulfide produced in seep systems as well.

Epsilonproteobacteria and Deltaproteobacteria were the most common intra-Proteobacteria association in the network and have been shown to co-occur in many sulfidic habitats (Campbell et al., 2006; Omoregie et al., 2008), where Epsilonproteobacteria oxidize sulfur and Deltaproteobacteria disproportionate or reduce sulfur species (Pjevac et al., 2014). In the network, Sulfurovum was associated with both SEEP-SRB1 and Desulfobulbus, and this was also seen in the relative sequence abundance data where Sulfurovum had high relative sequence abundance in all of the Deltaproteobacteria Magneto-FISH enrichments. Epsilonproteobacteria have been shown to oxidize sulfide to S° or HS⁻ to sulfate in microbial mats (Pjevac et al., 2014), allowing some sulfur substrate differentiation between these Epsilonproteobacteria groups in this system. Sulfurimonas was not strongly associated with any Deltaproteobacteria in the network analysis and only had high relative sequence abundance in the MSMX-Eel_932 enrichment (16S rRNA gene iTag, 16S rRNA gene bacterial, and soxB gene libraries; see Fig. S2 for further discussion of metabolic genes). CARD-FISH analysis using probe Epsi_404 confirmed the presence of Epsilonproteobacteria cells within some ANME and other non-hybridized cell-containing loose aggregates, but did not appear to be in the tight physical association characteristic of ANME/SRB consortia. While cultured representatives of these Epsilonproteobacteria have optimum growth with some oxygen present (Inagaki et al., 2003; Inagaki et al., 2004), it is possible that these uncultured methane seep Epsilonproteobacteria may be able to use other oxidants such as nitrate or intermediate sulfur species while in anaerobic incubation conditions.

In comparison to Delta- and Epsilonproteobacteria, there was only one Gammaproteobacteria OTU in the network (Thiohalobacter, with one Anaerolineaceae association). Cultured representatives of Thiohalobacter have diverse sulfur capabilities, including thiocyanate metabolism, but are not known to form associations with other sulfur cycling organisms (Sorokin et al., 2010). This differentiation between Gamma- and Epsilon-/Deltaproteobacteria has been seen in other systems such as sulfidic cave biofilms (Macalady et al., 2008) or in microbial mats on marine sediments (Pjevac et al., 2014). Gam_42a hybridizing cells (Gammaproteobacteria) were observed to form aggregates with non-ANME and non-Desulfobulbaceae (DSS) cells in our CARD-FISH analysis, but the identity of these organisms was not determined. While not recovered in the network, the majority of the Gammaproteobacteria OTUs observed by iTag from the both the bulk sediment and MSMX-Eel_932 Magneto-FISH 16S rRNA gene (Table 1) and aprA gene libraries (see Fig. S2 for further discussion of metabolic genes) were from the SILVA taxonomy endosymbiont clade. This endosymbiont clade houses organisms with a carbon-fixation/sulfur-oxidation metabolism (Duperron et al., 2012; Goffredi, 2010) and is predicted to be an important member of the sulfur and carbon cycles in marine sediments outside of an endosymbiotic lifestyle (Lenk et al., 2011).

There were also three unique, positive Deltaproteobacteria-Deltaproteobacteria associations observed in the network (Desulfobulbus/Desulfococcus, Desulfobulbus/SEEP-SRB4, Desulfoluna/SEEP-SRB2). These multiple intra-Deltaproteobacteria associations suggests there may be further nuances to be explored in the Deltaproteobacteria community structure, perhaps akin to the nitrate based partitioning observed between DSB and DSS in seep sediments (Green-Saxena et al., 2014). Desulfobulbus was also associated with SAR406, and SAR406 had high relative sequence abundance in the Delta495a enrichments. SAR406 (Marine Group A) fosmids contained polysulfide reductase genes that may be used for dissimilatory polysulfide reduction (Wright et al., 2014) Desulfobulbus can also use polysulfide, in addition to a range of other sulfur sources (Fuseler & Cypionka, 1995), potentially linking these two taxa.