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

Methane seep systems along continental margins host diverse and dynamic microbial assemblages, sustained in large part through the microbially mediated process of sulfate-coupled Anaerobic Oxidation of Methane (AOM). This methanotrophic metabolism has been linked to consortia of anaerobic methane-oxidizing archaea (ANME) and sulfate-reducing bacteria (SRB). These two groups are the focus of numerous studies; however, less is known about the wide diversity of other seep associated microorganisms. We selected a hierarchical set of FISH probes targeting a range of Deltaproteobacteria diversity. Using the Magneto-FISH enrichment technique, we then magnetically captured CARD-FISH hybridized cells and their physically associated microorganisms from a methane seep sediment incubation. DNA from nested Magneto-FISH experiments was analyzed using Illumina tag 16S rRNA gene sequencing (iTag). Enrichment success and potential bias with iTag was evaluated in the context of full-length 16S rRNA gene clone libraries, CARD-FISH, functional gene clone libraries, and iTag mock communities. We determined commonly used Earth Microbiome Project (EMP) iTAG primers introduced bias in some common methane seep microbial taxa that reduced the ability to directly compare OTU relative abundances within a sample, but comparison of relative abundances between samples (in nearly all cases) and whole community-based analyses were robust. The iTag dataset was subjected to statistical co-occurrence measures of the most abundant OTUs to determine which taxa in this dataset were most correlated across all samples. Many non-canonical microbial partnerships were statistically significant in our co-occurrence network analysis, most of which were not recovered with conventional clone library sequencing, demonstrating the utility of combining Magneto-FISH and iTag sequencing methods for hypothesis generation of associations within complex microbial communities. Network analysis pointed to many co-occurrences containing putatively heterotrophic, candidate phyla such as OD1, Atribacteria, MBG-B, and Hyd24-12 and the potential for complex sulfur cycling involving Epsilon-, Delta-, and Gammaproteobacteria in methane seep ecosystems.


INTRODUCTION
A central goal in microbial ecology is identifying and understanding microbial interactions in the environment. This goal can be addressed at many scales from statistical analyses of entire ecosystems (Barberán et al., 2012;Malfatti & Azam, 2010;Ruff et al., 2015;Steele et al., 2011;Sunagawa et al., 2015) to high resolution image analysis of specific symbioses (Malfatti & Azam, 2010;McGlynn et al., 2015;Orphan, 2009;Orphan et al., 2001b;Wegener et al., 2015). Previous studies have shown that complex datasets can be distilled to determine primary ecosystem drivers, such as temperature, as main predictors of community variability (Sunagawa et al., 2015). In addition to correlating microbial patterns to environmental factors, interspecies interactions can be evaluated with methods such as co-occurrence analysis (Friedman & Alm, 2012). Statistical significance of co-occurrence can be assessed at scales ranging from the entire genome to the operational taxonomic unit (OTU) (Barberán et al., 2012;Chaffron et al., 2010).
Here we combine Magneto-FISH and Illumina Tag (iTag) sequencing utilizing the Earth Microbiome Project (EMP) universal primer set (Caporaso et al., 2012). The Magneto-FISH method was originally developed to enrich for and characterize multi-species microbial associations in environmental samples (Pernthaler et al., 2008). This method consists of a liquid CARD (CAtalyzed Reporter Deposition)-FISH reaction as a 16S rRNA gene identity-based selection mechanism followed by an immunomagnetic sediment matrix separation mechanism to target specific phylogenetic groups in conjunction with their physically associated microbial partners. By combining this method for phylogenetically targeted physical separation with high throughput amplicon sequencing, we can compare an array of associated microbial communities in parallel, with replicates. This provides statistical power in deriving microbial associations from complex sediment community assemblages, and thereby improving hypothesis development.
Anaerobic methane-oxidizing (ANME) archaea and sulfate-reducing Deltaproteobacteria (SRB) are the predominant community members discussed in methane seep literature and form syntrophic partnerships in physical associations, termed ''aggregates'' or consortia (Boetius et al., 2000;Green-Saxena et al., 2014;Knittel et al., 2003;Orphan et al., 2001a;Schreiber et al., 2010). Since physical association appears to be an important element for consortia activity (McGlynn et al., 2015;Wegener et al., 2015), methods like Magneto-FISH are ideal for probing this system because target organisms are separated from the sediment matrix along with their physically associated partners. A hierarchical probe set was chosen targeting Deltaproteobacteria and their ANME partners to create nested Magneto-FISH enrichments from methane seep sediment incubations under methane headspace. This method allows us to examine potential physical associations between ANME and SRB taxa and other microorganisms using co-occurrence statistical methods applied to iTag sequences from nested Magneto-FISH enrichments. ANME have been broadly divided into three separate groups, which can be further subdivided into 1b,2a,2b,2c,and 2d, and 3. ANME-1 archaea are a unique order-level lineage within the Euryarchaeota, between the Methanomicrobiales and the Methanosarcinales, known to associate with sulfate-reducing bacteria, but obligately associated lineages have yet to be defined. ANME-2 archaea, within the order Methanosarcinales, commonly form associations with Desulfosarcina/Desulfococcus-related (DSS) sulfate-reducing Deltaproteobacteria (Boetius et al., 2000;Orphan et al., 2001a;Schreiber et al., 2010). They have also been found in association with Desulfobulbus-related (DSB) Deltaproteobacteria in the same environments, where geochemical factors have been suggested as a possible explanation for partner differentiation (Green-Saxena et al., 2014). ANME-2a/b and ANME-2c both predominately associate with a subgroup of DSS, SEEP-SRB1 (Schreiber et al., 2010), but also form consortia with DSB (Green-Saxena et al., 2014;Pernthaler et al., 2008). ANME-3 has been found in association with Desulfobulbus-related Deltaproteobacteria (Niemann et al., 2006) and SEEP-SRB1 (Schreiber et al., 2010). These ANME groups have also been observed in the environment without bacterial partners (House et al., 2009;Orphan et al., 2002;Schreiber et al., 2010;Treude et al., 2007). In addition to ANME archaea, other uncultured archaeal lineages commonly recovered from methane seeps include Marine Benthic Group-D (Thermoplasmatales), Deep Sea Archaeal Group/Marine Benthic Group-B (Ruff et al., 2015;Yanagawa et al., 2011), and sometimes methanogens (Orphan et al., 2001a;Ruff et al., 2015;Takano et al., 2013;Vigneron et al., 2015).
Despite the wealth of bacterial and archaeal diversity in methane seep sediments, little is known about potential associations with ANME/SRB, or associations that do not involve ANME or SRB. Our study utilizes the novel combination of targeted Magneto-FISH enrichment of specific microbial taxonomic groups and iTag sequencing to develop statistically supported co-occurrence microbial networks to address knowledge gaps in our understanding of methane seep microbial communities. Network analysis revealed many novel associations between methane seep Proteobacteria taxa and Candidate phyla. The significant co-occurrence observed by these OTUs suggests new avenues for future studies on microbial interactions involved in carbon and sulfur cycling in methane seep systems.

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.02 N 125 • 6.00 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 2 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) andTrembath-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 Table 2 16S rRNA gene iTag and clone library relative sequence abundance for seep microbiome OTUs. Relative sequence abundances were computed for the top 135 OTUs in the iTag dataset. These OTUs correspond to ∼55% of the total sequences in the unfixed bulk sediment. Bacterial and archaeal 16S rRNA gene libraries are included for the core methane seep taxa, with the total number of clones for each library indicated above. Core methane seep taxa were based on Ruff et al. (2015) and include: candidate Phylum Atribacteria, Candidate Division Hyd24-12, Methanomicrobia, Caldilineales, Desulfobacterales, and Spirochaetales. While we did recover other Chloroflexi, no Caldilineales were recovered in iTag or gene library sequencing so they are not included below. Fixed bulk sediment was chosen for baseline comparison (rather than unfixed) since it includes the potential loss of cells due to fixation and wash steps, thereby processed more similarly to the Magneto-FISH samples. An average and standard deviation for relative sequence abundance among replicates was calculated for each sample set. A ratio of the average relative sequence abundance of Magneto-FISH enrichments compared to the fixed bulk sediment value is reported (Rel. fixed). Ratios over 1.5 are underlined. 16S rRNA gene bacteria and archaea clone libraries for two Magneto-FISH enrichments and fixed bulk sediment are also included for comparison to recovered iTag diversity.
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)

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 Nodes represent the taxonomy of the OTUs in the network and edges are the connections between OTUs. Node size is scaled by number of connecting OTUs and colored by simplified, putative metabolic guild (sulfate reducer: blue, small dash; sulfur oxidizer: yellow, medium dash; archaeal methanotroph: magneta, large dash; mixotroph: green, no dash; heterotroph: brown, no outline). Edge thickness is scaled by number of occurrences of this association (from 50 to 100 times) and number of occurrences also included along edge. Negative associations are denoted by hashed lines. The combined network is displayed using Cytoscape, with the average correlation coefficient across all runs determining the distance between nodes and the number of occurrences in 100 network iterations determining edge width.
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.
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 Table 3 iTag relative abundance of remaining 'non-core' methane seep microbiome OTUs. Relative sequence abundances were computed for the top 135 OTUs in the iTag dataset that were not included in the core methane seep microbiome. An average and standard deviation for relative sequence abundance among replicates was calculated for each sample set. A ratio of the average relative sequence abundance of Magneto-FISH enrichments compared to the fixed bulk sediment value is reported (Rel. fixed). Ratios over 1.5 are underlined. 16S rRNA gene bacteria and archaea clone libraries for two Magneto-FISH enrichments and fixed bulk sediment are also included for comparison to iTag enrichment.   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 sulfuroxidizing 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.  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.
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.

CONCLUSIONS
Our findings support the utilization of paired Magneto-FISH and iTag sequencing in developing and testing hypotheses to interrogate complex interactions in microbial communities. Contaminants and amplification bias can be identified and mitigated with diversity assessment by multiple means (i.e., multiple iTag primer sets, FISH surveys, or non-16S rRNA gene surveys) and parallel processing of control samples (bulk sediment and no-template) along with Magneto-FISH enrichments. Since it may not always be known a priori which taxa are in an environmental sample, sequencing of a defined mock community may not be an option for assessing bias. However in our case, prior knowledge of major seep taxa enabled assessment of amplification bias by iTag. It should also be noted that the degree of bias was more pronounced in the environmental samples than our mock samples, therefore mock community samples may not fully capture the degree of bias, but can be useful in identifying which taxa may be the most biased. We found the bulk sediment 16S rRNA gene libraries to be the most useful for determining which of the most abundant taxa were affected by amplification bias. Future studies may benefit more from bulk sediment analysis by a range of iTag primer sets or gene libraries to assess potential sequencing biases in a new microbial community.
Multiple statistical methods supported differences between Magneto-FISH enrichments and the bulk sediment. We also found variation between SparCC network computations. Therefore, we added confidence to network associations by reporting the number of times an association was recovered out of 100 co-occurrence iterations along with correlation and p-value.
Our resultant microbial community network had many statistically significant methane seep taxa correlations beyond the common ANME/SRB association. The downplay of anaerobic methanotrophs in our iTag sequencing may have had the beneficial effect of bringing fermenters to the forefront, highlighting their complex role in methane seep microbial communities. Within the core methane seep microbiome taxa, there were strong associations between Atribacteria and Hyd24-12 and Deltaproteobacteria, but no direct association between Atribacteria and Hyd24-12. This may indicate a different niche for these two currently uncultured groups in methane seep systems. Sulfurovum and Sulfurimonas were differentiated as either Deltaproteobacteria-associated or archaea-associated, respectively. There were statistically significant associations between Deltaproteobacteria and non-Proteobacteria, such as the Planctomycetes sediment group 'AKAU3564,' and groups that contained neither SRB nor ANME but had high statistical significance, such as MBG-B and OM1. Future development and application of more specific FISH probes will assist in further hypotheses development and testing of these associations in Hydrate Ridge methane seeps.
Some groups, such as Gammaproteobacteria, appeared to have associations with other microbes based on broad FISH surveys and Magneto-FISH relative sequence abundance data, but were not recovered in the network analysis. Determination of the specific Gammaproteobacteria involved in associations via FISH probe development or other means (Hatzenpichler et al., in review) will also aid in refining why associations might be missed in the microbial network analysis based on DNA taxa co-occurrence. In summary, a continual feedback loop between microbial identification and isolation techniques and gene based statistical analyses is required to tease apart interactions within complex microbial systems. The combination of Magneto-FISH and high throughput, parallel iTag sequencing provides an effective bridge between these two modes.