Cold Seeps on the Passive Northern U.S. Atlantic Margin Host Globally Representative Members of the Seep Microbiome with Locally Dominant Strains of Archaea

The four cold seep sites investigated on the northern USAM were located at water depths ranging between 335 and 1,545 mbsl along the continental slope (Fig. 1). Seeps were identified by signs of active methane seepage on the sediment surface, namely, white, filamentous microbial mats, live mussels, and/or visible methane gas bubbles (images available in reference 40). At Shallop Canyon East, the shallowest site, methane concentrations averaged 220 μM in a Seep core (C5) sampled from a bacterial mat (Table S1). At New England Seep, pore water methane concentrations were highest (Fig. 2), averaging 780 μM across two separate Seep cores (C1 and C9) sampled from bacterial mats (Table S1). At Veatch Canyon, the deepest site, methane concentrations were lowest (Fig. 2), averaging only 16 μM. Veatch Canyon seeps lacked thick bacterial mats overlying sediment, so Seep cores (C5 and paired cores C2 and C8) were retrieved from sediments covered with thin mats or thick mussel beds (Table S1). Methane concentrations were not measured at Shallop Canyon West (C7). At all sites, methane concentrations were higher in Seep samples than in Near-Seep or Background samples (Fig. 2).

Sulfate and sulfide were negatively correlated in all sediment cores. Higher concentrations of sulfide were observed in cores with higher methane concentrations, consistent with sulfate-dependent AOM. Only New England Seep cores demonstrated a complete conversion from sulfate to sulfide; in all other cores, sulfate was more abundant (Fig. 2). Acetate concentrations decreased with sediment depth in most cores (Fig. 2) but did not show a trend between sites or sediment types (Seep, Near-Seep, or Background). The highest concentrations of acetate (around 100 μM) were observed in a Seep core from Veatch Canyon (C5) and were nearly double those in all other Seep and non-Seep (Near-Seep and Background) samples.

Concentrations of ammonium, nitrate, and nitrite were generally low throughout all samples (Fig. S1). Inconsistent with trends at previously described seeps (41, 42), ammonium was depleted (<15 μM) throughout Seep cores from Shallop Canyon and New England Seep relative to non-Seep cores at those sites, where concentrations were <0.5 μM at the surface but increased to >50 μM at depth. In comparison, ammonium concentrations in Veatch Canyon Seep samples were more elevated (19.2 μM to 61.3 μM) and were within the range observed in surrounding non-Seep sediment at that site (<0.5 μM to ~100 μM). Nitrate and nitrite typically decreased with depth and remained below concentrations of 40 μM and 2 μM, respectively. Both nitrate and nitrite were low at New England Seep, with the exception of the uppermost horizon in C9, where concentrations were 6.14 and 0.92 μM, respectively (Fig. S1).

Concentrations of phosphate (0.5 and 12 μM), bromide (~0.6 to 1 mM), and chloride (~500 to 600 mM) were also measured and demonstrated no trends or variable trends with sediment depth (Data Set S1).

In 26 Seep samples (7 cores) from the four USAM sites, 9,280 unique 16S rRNA gene/16S rRNA amplicon sequence variants (ASVs) were recovered (7,118 in DNA extracts and 4,679 in RNA extracts). The most abundant class-level taxa across all sites (Fig. 3A) were members of the previously documented “seep microbiome” (2, 15, 36), and included Methanomicrobia (which encompasses the ANME archaea), Delta- and Gammaproteobacteria, and Bacteroidia (Fig. 3A). The JS1 lineage was also present at high relative abundances (around 10 to 20%) at the two deeper seep sites, although it made up less than 0.5% of the community at Shallop Canyon. Trends with sediment depth could be observed for some taxa, including a decrease with depth in Gammaproteobacteria (containing sulfide-oxidizing families Beggiatoaceae and Thiotrichaceae and aerobic methane-oxidizing Methylomonaceae) (Fig. S2) and an increase in Methanomicrobia and JS1 (Fig. 3A).

Despite the consistent presence of certain class-level taxa across seep sites (Fig. 3A), heterogeneity was observed at finer taxonomic scales (e.g., family or genus). ANME subgroups within Methanomicrobia displayed particularly different relative abundances across sites (Fig. 4A). At Shallop Canyon East and West, ANME-2 subgroups (specifically, ANME-2a and ANME-2c) made up >90% of Methanomicrobia (Fig. 4A). Relative abundances of ANME-2a decreased with sediment depth, while those of ANME-2c increased. At New England Seep, ANME-1 rather than ANME-2 made up the largest percentage of Methanomicrobia: ANME-1a in C1 and ANME-1b in C9. Veatch Canyon was the most heterogeneous site in terms of ANME composition; all major ANME groups were present, with no sediment depth-related trends. ANME-3 was seen at low abundance in three of four sites (all but Shallop Canyon East) and was mostly detected in samples where the total number of reads assigned to Methanomicrobia exceeded 1,000 (Fig. 4A).

The different SRB clades commonly associated with ANME archaea in AOM consortia also varied in relative abundance across the four seep sites (Fig. 4A). The Seep-SRB2 clade of Desulfobulbaceae, which has been observed in association with ANME-1 (43) and ANME-2c (29), was particularly abundant (accounting for 10 to 30% of the Deltaproteobacteria) at New England Seep, as well as in the shallow, paired cores (<10 cm apart) C2 and C8 at Veatch Canyon. The Seep-SRB1 group within the Desulfobacteraceae, which is highly adapted to symbiotic relationships with ANME-2a/c subgroups (30, 34, 41, 44), as well as ANME-3 (44), appeared in high abundances in all Seep cores (comprising 20 to 60% of the Deltaproteobacteria) and was highest in C1 at New England Seep. This group is a potentially obligate symbiont apparently dependent on the presence of methane for survival (44–47). Seep-SRB4 was more abundant in the upper sediment horizons than the lower and has been previously pinpointed as the partner of ANME-2a/b (48).

In addition to 16S rRNA gene sequences (DNA) (Fig. 3A and 4A), we assessed the relative abundance of 16S rRNA sequences (RNA) (Fig. 3B; Fig. S3a) and found that the class-level trends were extremely similar in both data sets. However, sequences of all major classes of Archaea, with the exception of Methanomicrobia in a Near-Seep core at Shallop Canyon West (C3), were relatively less abundant in the RNA data set than the DNA data set. Relative abundance of rRNA is a not an accurate proxy for relative activity between taxa but can be useful for the more conservative goal of indicating which taxa are likely translationally active (49). The detection of rRNA attributed to all major groups (e.g., ANME subgroups, families of Delta- and Gammaproteobacteria, and JS1) for which 16S rRNA genes were detected suggests that most microorganisms in USAM Seep and non-Seep sediments are active on some level. However, we refer to the detection of rRNA as reflecting potential activity to emphasize the limitations of rRNA to indicate true activity.

In 38 non-Seep samples (8 cores), 23,477 unique 16S rRNA gene/16S rRNA ASVs were recovered (17,256 in DNA extracts and 10,618 in RNA extracts). This was more than twice as many unique ASVs as were recovered in Seep samples. As was observed in Seep samples, microbial communities across sites were similar to one another at the class level. Delta- and Gammaproteobacteria were most abundant, while Alphaproteobacteria, Phycisphaerae, Planctomycetacia, and Nitrososphaeria had high relative abundances in comparison to Seep cores (between 5 and 20 times greater) (Fig. 3A). Across all sites, there was no significant difference between Near-Seep and Background microbial communities at the class level (analysis of similarity [ANOSIM]; P = 0.089) or at the ASV level (ANOSIM; P = 0.077). In contrast, there was a significant difference between microbial communities in Seep and non-Seep samples at each site, both at the class level (ANOSIM; R = 0.351, P = 0.001), and at the ASV level (ANOSIM; R = 0. 409, P = 0.001).

To more specifically target methane-cycling community members, genes and transcripts of the alpha subunit of methyl coenzyme M reductase (mcrA) were also amplified and sequenced from the seep sediments at all four sites. The mcrA analysis generated ~15,000 ANME reads per Seep sample (average) versus ~1,500 in the 16S RNA gene/16S rRNA analysis. Neither mcrA genes nor transcripts were detected in Near-Seep or Background samples. Relative abundance patterns in the mcrA data set (Fig. 4A; Fig. S3b) broadly matched those in the 16S rRNA data set (Fig. 4A; Fig. S3a); ANME-1 was most abundant at New England Seep and in Veatch Canyon paired cores C2 and C8, and ANME-2c was most abundant at Shallop Canyon and in the third Veatch Canyon core (C5). However, ANME-2a and -2b were present in only low numbers across the entire mcrA data set (Fig. 4B), despite constituting ~50% of the Shallop Canyon community in the 16S rRNA data set (Fig. 4A). In contrast, ANME-3 was more relatively abundant in the mcrA data set than in the 16S rRNA gene/16S rRNA analyses, particularly at Shallop Canyon. It is possible that the mcrA primers used do not capture the full extent of the ANME-2a/b sequence space and that this led to an underestimation; other studies using these primers have detected no ANME-2a/b (50, 51) or only a few clones in a clone library (47). An underestimation of ANME-2a/b in Shallop Canyon sediments may therefore result in an overestimation of ANME-3.

The mcrA gene analysis also reveals that although 1,093 unique mcrA ASVs were recovered across the data set, 94% of which are ANME, the ANME population in a given sample is typically dominated by a very small number of mcrA ASVs (Fig. 5; Fig. S4). At Shallop Canyon East and West, for example, only two ANME-2c ASVs make up ~50% of the total mcrA genes (Fig. 5; Fig. S5) and ~30% of the transcripts (Fig. S4 and S6). New England Seep and Veatch Canyon are slightly less uneven, but the 3 or 4 most abundant ASVs still make up ~25% of the total mcrA genes (Fig. 5) and ~40% of the mcrA transcripts (Fig. S4). Furthermore, the most abundant mcrA genes (Fig. S5) and transcripts (Fig. S6) are the same throughout each core. This is in contrast with previous reports of ANME subgroup separation by depth elsewhere (37, 38) and is potentially a result of the relatively small changes in geochemistry with depth observed here (Fig. 2). With the exception of the cores from Shallop Canyon East and West (which share many ASVs), very few of the 10 most abundant mcrA gene or transcript ASVs in each core were also within the 10 most abundant reads in other cores (Fig. S7a), even other cores taken within the same study site. However, those ASVs were occasionally present at other sites and cores, albeit in much lower abundances (Fig. S7b).

Our deep-sequencing approach allowed us to investigate potential associations between ANME and Deltaproteobacteria at USAM seeps using a correlation analysis of 16S rRNA gene and 16S rRNA ASVs across all Seep samples. We used maximal information-based nonparametric exploration (MINE) (52) to create a network diagram of significant correlations (Benjamini-Hochberg-corrected P values < 0.01) (Fig. 6). Correlations between groups are potentially indicative of symbiotic relationships but may also be the result of independent responses to physicochemical drivers. Still, correlation analyses can be effectively used to identify potential symbionts and to generate hypotheses that can be tested with further analyses (for example, see reference 53).

The largest maximal information coefficient in our results (0.928)—and therefore the strongest correlation—was between an ANME-1b ASV and a Seep-SRB1 ASV (Data Set S2). However, the greatest number of correlations existed between ANME-2c and Seep-SRB1 ASVs, consistent with the previously demonstrated symbiosis between these groups. Interestingly, ANME ASVs at USAM seeps occupied central locations in our network, suggesting that they correlated with numerous SRB ASVs of different taxonomic groups. In contrast, SRB ASVs were more often correlated with just one or two ANME ASVs (Fig. 6). One exception was the newly identified Seep-SRB1g (53). While its most significant correlation occurred with ANME-2b, consistent with the recently described symbiosis, other correlations were observed between it and ANME-1a, ANME-2a, and ANME-2c. Our results, while correlative rather than direct observations of symbioses, suggest that ANME may be more promiscuous than SRB at the ASV level in forming syntrophic relationships to support the anaerobic oxidation of methane.

Nonmetric multidimensional scaling (NMDS) ordinations of 16S rRNA genes (Fig. 7A) and 16S rRNA (Fig. 7B) were well represented in two dimensions, with stress values of 0.150 and 0.105, respectively. Microbial communities in each sample were primarily differentiated by sediment type (Seep or non-Seep), which was found to be statistically significant by ANOSIM (R = 0.313; P = 0.001). Furthermore, while Seep samples were significantly differentiated by site (ANOSIM; R = 0.136, P = 0.044), non-Seep samples were not (ANOSIM; R = 0.091, P = 0.124). USAM seep sites therefore contain microbial communities that are more distinct from one another than those of equally spaced background sediments.

Principal-component analysis (PCA) was used to examine whether the greater site-specific community differentiation in Seep samples was a result of greater geochemical variation within Seep samples than within non-Seep samples (Fig. S8). While Seep samples encompassed a greater range in concentration in some geochemical parameters (e.g., methane and sulfide concentrations; PCA axis 1), non-Seep samples included more heterogeneity than Seep samples in other parameters (e.g., nitrogen species, acetate, bromide, and chloride; PCA axis 2). Therefore, greater community differentiation in Seep samples than in non-Seep samples is likely due to the variation in specific key parameters, rather than to heterogeneity overall.

To understand the influence of environmental parameters on the microbial composition of USAM seeps and specifically identify the variables that best explain community variation, we carried out distance-based redundancy analyses (dbRDA) using stepwise model selection. Redundant variables were discarded from the analyses. Early models, which included only samples with a full suite of corresponding geochemical measurements (23 samples of 64), demonstrated that methane was a significant variable in explaining microbial community variation in 16S rRNA genes but not in 16S rRNA (Fig. S9). To nearly double the number of samples and cores in the analysis, 15 additional samples, which were missing only measurements of methane concentration, were added to the final model. With this expanded data set (where methane was not analyzed), the best distance-based linear model of the DNA community indicated that nitrate, sulfide, acetate, and ammonium concentrations, as well as sediment depth and water depth, were significant predictive variables (Benjamini-Hochberg-corrected P value < 0.05), accounting for 53.37% of the total variation in taxonomic composition (Fig. 7C) (sediment depth = 14.7%, nitrate = 14.6%, sulfide = 11.6%, acetate = 4.6%, water depth = 4.3%, ammonium = 3.5%). The primary axis (explaining 25.33% of the variation) separated shallower (0 to 6 cm below the seafloor [cmbsf]). Seep sediment samples from deeper Seep sediment samples and non-Seep samples, while the secondary axis (explaining 13.21% of the variation) separated Seep samples from non-Seep samples. The best distance-based model of the RNA community was highly similar, though predictive variables did not include water depth, and accounted for 39.22% of the total variation in taxonomic composition (sediment depth = 11.7%, sulfide = 11.4%, nitrate = 7.1%, acetate = 5.0%, ammonium = 4.1%). The primary axis (explaining 18.27% of the variation) separated sulfidic New England Seep samples from the others, while the secondary axes (explaining 8.61% of the variation) separated Seep samples from non-Seep samples (Fig. 7D).

To estimate the relative contribution of these geochemical parameters (and other deterministic processes) in the assembly of USAM communities, as opposed to random, stochastic processes, we applied a null model-based framework called the normalized stochasticity ratio (NST) (54). Stochastic, neutral processes (such as random birth/death events, probabilistic dispersal, and ecological drift) generate communities by random chance, while deterministic processes (environmental factors or species interactions) shape communities by fitness. We calculated the phylogenetic NST (pNST) between communities in the DNA data set and found that deterministic processes dominated assembly in Seep sediments, as indicated by a pNST of <50% (pNST = 38.17). Stochastic processes dominated assembly in non-Seep sediments, as indicated by a pNST of >50% (pNST = 59.09). The difference between the pNSTs of the two groups was determined to be statistically significant (permutational multivariate analysis of variance [PERMANOVA]; P < 0.001).

To determine whether the distinct geologic setting of the USAM or its distance from other characterized seeps influenced the composition of its microbial communities, USAM samples were compared to those from other globally distributed sites. While cold seeps have been studied around the world, the microbial communities of surprisingly few have been investigated with deep 16S rRNA gene sequencing (≫100 sequences per sample). We included all published studies with overlapping priming regions and Illumina amplicon sequencing data in our comparison (n = 6). These studies were performed on samples from Eel River Basin, Hydrate Ridge, Juan de Fuca, Astoria Canyon, Coquille, Klamath Knoll, Barkley Canyon, Heceta, and Clayoquot Slope along the Cascadia Margin (36, 55, 56), Oki Trough in the Sea of Japan (57), Haima in the South China Sea (8), and Storfjordrenna in the Barents Sea (58) (Fig. 8A). Most of these cold seep sites are influenced by tectonic activity, including along the Cascadia Margin (59, 60), within the Sea of Japan (61), and at Storfjordrenna (62, 63). Haima cold seep is located on the northwestern continental margin of the South China Sea and is the only other passive margin seep in our data set.

By comparing the total microbial community composition in all samples, we found that the USAM seeps did not contain a distinct subset of the global seep microbiome (Fig. 8B). The average weighted UniFrac distance between USAM seep samples (0.307) was not significantly different from that of equivalently sized collections of samples randomly chosen from the global data set, indicating that the USAM contains a globally representative collection of seep microorganisms. However, we performed the same comparison on the Archaea only (Fig. 8C) and found that the USAM seeps do contain a distinct subset of global seep Archaea. The average weighted UniFrac distance between USAM seep archaeal communities (0.166) was significantly smaller than that of equivalently sized collections of samples randomly chosen from the global data set. A hierarchical clustering dendrogram further revealed the relative similarity between USAM Archaea (Fig. S10).

To determine whether differences in microbial community composition between global sites are related to seep proximity, we compared the weighted UniFrac distance between pairs of samples with their geographic separation (Fig. 8d). Samples from the same latitude and longitude coordinates were removed. In a logarithm-transformed space, the distance-decay relationship was significant but very weak (ANOVA; P < 0.001, r² = 0.037). To test whether such a relationship existed within the USAM alone, the analysis was repeated by including only pairwise comparisons between USAM seep samples (Fig. 8D, inset). The relationship between geographic separation and phylogenetic relatedness at this smaller scale was not statistically significant (ANOVA; P = 0.722). To confirm that our paired cores (Veatch Canyon C2 and C8, 10 cm apart) had not skewed the results, we repeated the analyses without this sample pair (Fig. S11) and obtained the same results.

Most cold seeps with well-characterized microbial communities, such as those in the northeastern Pacific, are influenced by tectonic activity. Along active margins, subduction of oceanic plates creates fractures in the seabed, allowing deeply sourced, methane-rich fluids to advect upward (64). Upward flow of methane along passive margins instead depends on gravitational pressures (65), a constraint that, until recently, was thought to limit passive margin seepage to sedimentary basins, submerged aquifers, or locations with rapid sediment accumulation (66). Along the passive, northern USAM, cold seeps and methane gas plumes had not been expected. Nonetheless, exploration and sampling of the area revealed the widespread existence of seeps (3) and their characteristic surficial communities, including mussels and sulfide-oxidizing microbial mats (40). Whether these seeps hosted subsurface microbial communities similar to those described previously was unknown, and it has implications for methane oxidation at passive margin seeps worldwide.

Our investigation of the sediment-dwelling microbial communities at four USAM cold seep sites revealed that the seep microbiome—the collection of class level taxa common at cold seep environments (e.g., ANME archaea, Delta- and Gammaproteobacteria, and JS1 bacteria)—is present and potentially active. Additionally, we found evidence in our correlational analysis for many of the microbial associations described previously at other seeps, including that between ANME-2c and Seep-SRB-1, suggesting that methane is oxidized via known symbioses. Furthermore, our global analysis that shows the USAM seeps do not form a distinct phylogenetic subset of the global seep microbiome. In fact, the phylogenetic variation observed between seeps at the USAM, all within 50 km of one another, spanned that of the entire global data set (Fig. 8b). Our combined analysis demonstrates that the microbial community composition and ecology of cold seeps at the USAM, and potentially that of seeps on passive margins more generally, is not systematically different from that of seeps on active margins. While valuable in itself, the finding that these seep communities are not systematically distinct also provides an opportunity to identify potentially broadly applicable trends in seep microbiology.

While the microbial communities at USAM cold seeps were similar to one another and those at globally distributed seeps at the class and order level, they were distinct at the genus and ASV levels (Fig. 4A and 7A and B). Within the ANME archaea specifically, there were noticeable differences in the presence and abundance of taxonomic subgroups by site. Each core contained a completely different ANME profile at both the subgroup (approximately genus) and ASV levels (Fig. 4; Fig. S3). Although heterogeneity in genus-level ANME distribution has been observed previously (see below), the microdiversity and distribution within these broader groups are not well known. Interestingly, despite the presence of 257 mcrA ASVs per core on average, around 2 to 4 individual mcrA ASVs in each Seep core were responsible for up to half of the mcrA community (Fig. S5) and potential activity (Fig. S4 and S6) in that core. This suggests that despite considerable strain-level diversity within each core, a small number of ANME phylotypes may be responsible for a large percentage of the methane oxidation.

We hypothesize that there could be a relationship between seep age and ANME community evenness at the ASV level. Although the age of these seeps is not known, bottom waters along the USAM are known to have warmed in the last 5,000 years by 8°C (4, 7), which is currently triggering the destabilization of buried hydrates and resulting in new or increased methane seepage (3, 7). The two Shallop Canyon sites may be younger than the others due to their shallower water depth, on the edge of the hydrate stability zone. This is supported by the much less well developed surficial indicators of seepage, including a lack of carbonate mounds, clam/mussel fields, and thick bacterial mats. Instead, the Shallop sites are characterized by patchy and thin bacterial mats and only sparse macrofauna (40). A connection between community age and evenness has been observed across a variety of species and environments, including epilithic biofilms (67), plant communities after glacial retreat (68), and microbial communities at whale falls (69). In a rare observation of a newly emerged seep, ANME communities in McMurdo Sound, Antarctica, took up to 5 years to emerge after seep development (70), and ANME communities interpreted from 16S rRNA gene sequencing were predominantly made up of three ASVs from one subgroup, ANME-1. The three ASVs detected were more than 97% similar. As more ASV-level studies are performed at additional seep sites (on the USAM and elsewhere), it will become clear whether the local strain dominance observed here is a function of seep age or rather a common feature of ANME populations.

As the climate changes, new sources of methane are expected to emerge on passive margins worldwide, and seep microorganisms will serve as the first barrier to reduce these new methane emissions to the atmosphere. How and how quickly these new seeps are colonized could affect the efficiency of methane oxidation and therefore the effectiveness of this barrier. To better understand this process, we leveraged our deep-sequencing data sets to statistically infer whether deterministic or stochastic processes drove the development of seep microbial communities on the USAM based on the genetic similarity of communities within each sediment core. Our results indicate that USAM cold seep communities are structured primarily by deterministic factors (e.g., environmental selection), while stochastic processes (e.g., random birth/death events, probabilistic dispersal, and ecological drift) dominate assembly within nonseep sediments. While deterministic factors can drive communities to become more or less distinct than the null expectation, microbial communities in USAM seep cores were generally found to be more different from one another than surrounding background sediments (Fig. 7A and B and 8B and C). Our analysis therefore suggests that site-specific differentiation at seeps is a result of niche-driven assembly, likely as a result of specific geochemical variations.

In USAM sediments, the relative abundance of ANME subgroups varied by location (Fig. 4), which has been clearly documented at other cold seep sites (37, 38, 71, 72) and is potentially indicative of separation into ecological niches. Previous work has suggested that dominance of ANME organisms, specifically, could be related to the concentrations of methane or sulfate/sulfide in cores. In the Nyegga cold seep sediments off the Norwegian coast, ANME-2a/b abundances increased in areas with low sulfide concentrations (37). Similar patterns were observed in Black Sea microbial mats (72), consistent with studies suggesting ANME-2 is sensitive to hydrogen sulfide (73). In the Japan Sea, ANME-2 was found in high-sulfate low-sulfide sediments and ANME-1 in low-sulfate high-sulfide sediments (38). In both Guaymas Basin and Nyegga cold seeps, ANME-2 decreased with methane flux, while ANME-1 increased (37, 74). While we did not detect a significant trend in ANME distribution with methane concentrations, this could be due to our relatively small suite of seep samples with methane measurements (only ~40% of the total samples; n = 26) and/or the error associated with measuring gas concentrations in depressurized cores. We did, however, detect significant correlations between the entire seep community and several physicochemical parameters, including nitrate, sulfide, acetate, and ammonium concentrations (Fig. S6c and d). In particular, the importance of sulfide (11.7% DNA; 11.4% RNA), which is a product of sulfate reduction that varies inversely with sulfate, nitrate (14.6% DNA; 7.1% RNA), and sediment depth (14.7% DNA; 11.7% RNA) suggests that the availability of specific electron acceptors and the redox state of the sediment are primary drivers of ASV-level selection at seeps and potentially therefore the differentiation between seeps.

To directly test the role of local dispersion in seep assembly, we explored distance-decay relationships in our data, as well as in globally available seep data. Consistent with the results of our ecological stochasticity analysis, our observations suggest that the composition of neighboring seeps (and therefore local dispersion mechanisms) was not an important determiner of microbial community composition. While paired seep cores (<10 cm apart) are remarkably similar (Fig. 8D), samples in seep cores collected ~425 m apart (our next closest sample pairs) were not. Although a significant distance-decay relationship was observed across the global data set, the low r² value (0.037) suggests that community similarity at the ASV level is only minimally related to geographic separation. In addition, this relationship could in fact be an artifact reflecting the increased likelihood that sample pairs with a small map distance were collected and processed synchronously, with the same methods. The lack of a strong distance-decay relationship suggests that seep communities are colonized primarily by global rather than local dispersal mechanisms, with strong site-specific selection.

In conclusion, we provide spatially resolved geochemical and molecular data to characterize the microbiology at four recently discovered seeps along the USAM. Our data, while demonstrating a lack of systematic differences between microbial communities at seeps on passive versus active margins, provides insight into the strain-level microdiversity and distribution of ANME archaea, as well as the overarching assembly processes and specific environmental parameters shaping seep community structure. To further constrain the spatial scales at which distance decay relationships might operate at cold seeps, a data set including deep sequencing from more seeps, and specifically from seeps separated by distances of 10 cm to 500 m, is necessary. We recommend that location data be reported as specifically and as accurately as possible when sequencing data are made publicly available to facilitate such meta-analyses. Given the likelihood of continued hydrate dissociation on the upper continental slope of passive margins worldwide, and the demonstrated diversity of methanotrophic taxa across this passive margin seep system, understanding the provenance of cold seep communities will be an important avenue for future research.

The four USAM cold seep sites investigated in this study were Shallop Canyon East (335 to 366 mbsl), Shallop Canyon West (390 to 395 mbsl), New England Seep (1,130 to 1,252 mbsl), and Veatch Canyon (1,407 to 1,545 mbsl). The shallowest sites, Shallop Canyon East and West, were located on the upper rim of Shallop Canyon, in gently sloping seafloor. Backscatter data had previously indicated the presence of roughly 15 gas plumes in this area (40). New England Seep was located 13 km southwest of Shallop Canyon on a steeply sloping, northwest-to-southeast-trending ridge, where 5 gas plumes had been previously detected (40). Veatch Canyon, the deepest site sampled, was located 27 km southwest of New England Seep on a steeply sloping, northeast-to-southwest-trending ridge at the base of Veatch Canyon, where 13 gas seeps had been inferred from backscatter data (40).

Sediment was collected from all four sites in July and August of 2016 aboard the R/V Atlantis (cruise AT36). At each site, sediment push cores up to 20 cm long were collected via the HOV Alvin submersible, either within a cold seep or from sediment several meters away. These were categorized as Seep cores and Near-Seep cores, respectively (Table S1). Seeps were delineated by the presence of white, filamentous microbial mats, live mussels, and/or visible methane gas bubbles. The filamentous mats were later confirmed to be rich in bacterial taxa known to oxidize sulfide (40). In addition, multicores were collected using an MC-400 multicoring device equipped with a Nikon D3300 24.2 MP DSLR to collect background sediment from the area surrounding each site; these were categorized as Background cores (Table S1).

Onboard, cores were kept at 4°C until extruded from push core liners and sectioned into 3-cm horizons within 24 h of collection. Several 1-mL subsamples of each depth horizon were immediately flash frozen in liquid nitrogen and preserved at −80°C for later DNA and RNA analyses. For methane measurements, 3-mL subsamples were transferred into 25-mL butyl-rubber-sealed vials filled with 5 mL of 5 M sodium hydroxide solution. Paired, unextruded cores were processed for pore water geochemistry at the same 3-cm horizons, using Rhizon samplers inserted through predrilled holes. Pore water (0.5 mL) was fixed with 0.5 mL of 0.5 M zinc acetate and stored at 4°C; the remaining pore water was stored at −20°C.

Headspace methane concentrations were measured from butyl-rubber-sealed vials (see above) using gas chromatography coupled with flame ionization detection (Shimadzu GC-2014; Boston University, Boston, MA, USA) and back-calculated to pore water concentrations using porosity values determined via weight loss after dehydration (and assuming a sediment density of 2.65 g cm⁻³). Sulfide concentrations were measured in triplicate and determined colorimetrically from the zinc acetate-preserved pore water samples using the methylene blue method (75), with a detection limit of 0.02 mM. All other assays were performed without replication due to limitations on pore water volume. Ammonium concentrations were determined colorimetrically using the indophenol blue method (76), with a detection limit of 0.5 μM. Nitrite and nitrate (NO_x − nitrite) concentrations were determined colorimetrically using the vanadium(III) chloride method (77), with detection limits of 0.1 μM and 1 μM, respectively. Concentrations of sulfate, phosphate, chloride, bromide, and acetate were determined by ion chromatography on frozen pore water samples (Dionex DX-500 ion chromatograph; Stanford University, CA, USA). Detection limits were 100, 0.05, 200,000, 100, and 10 μM, respectively.

RNA and DNA were extracted from flash-frozen sediments using the RNeasy PowerSoil total RNA isolation kit and the RNA PowerSoil DNA elution accessory kit (MoBio Laboratories, Carlsbad, CA, USA) from sediments flash frozen on board and stored at −80°C. RNA extracts were cleaned with the Ambion Turbo DNA-free kit (Thermo Fisher Scientific, Waltham, MA, USA), and reverse transcription of RNA to cDNA was completed using Superscript III first-strand synthesis Supermix (Thermo Fisher Scientific, Waltham, MA, USA).

DNA and cDNA were concentration normalized and amplified using a two-step PCR plan for Illumina amplicon sequencing. In the first step, universal primers 515F-Y and 926R (78) were used to target the V4-V5 region of the 16S rRNA gene sequence. In addition, two separate primer pairs, mcrA_F/mcrA_R (47, 79) and mcraMF/mcraMR, were used to target mcrA genes and transcripts. The latter primer set was designed in house and modified from reference 80. All three sets of primers included an extension complementary to the primers used in the second PCR. The gene-targeting region of the primer sequences are listed in Table 1. The mcrA_F/mcrA_R primer pair amplified the DNA and cDNA from more samples than mcraMF/mcraMR, likely due to increased efficiency with reduced degeneracies (Table S2). PCRs (25 μL) were performed with mixtures containing 0.5 μL of forward and 0.5 μL of reverse primers (10 μM concentration), 10 μL 5PRIME HotMasterMix (2.5×; Quanta-Bio, Beverly, MA, USA), 13 μL DNase-free water, and 1 μL DNA or cDNA template. The thermal cycling conditions were as follows: initial denaturing at 95°C for 180 s; 25 cycles of 95°C for 45 s, 50°C for 45 s, and 68°C for 90 s; a final elongation step at 68°C for 300 s; and refrigeration at 4°C until removal and storage.

In the second step, Illumina adaptors, barcodes, and indices were added to the amplicons. The same PCR mix was used with custom primers targeting the primer extension in the first PCR. The thermal cycling conditions were as follows: initial denaturing at 95°C for 180 s; 8 cycles of 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s; a final elongation step at 72°C for 300 s; and refrigeration at 4°C until removal and storage. Amplicons were cleaned with 0.7× AMPure XP magnetic beads (Beckman-Coulter, Brea, CA, USA), pooled, and quantified before being sent to the UC Davis DNA Technologies Core Facility (Davis, CA, USA) for Illumina MiSeq 2 × 250 bp (16S rRNA) or 2 × 300 bp (mcrA) sequencing. Nine 16S rRNA and five mcrA samples were randomly chosen for duplicate amplification. The average weighted UniFrac distance between duplicate 16S rRNA samples was 0.080. Negative (molecular-grade water) and positive (mock communities of known composition) controls were processed and sequenced in parallel with the samples. Lack of DNA contamination in the RNA extracts was confirmed by processing RNA extracts (without reverse transcription) in parallel and seeing no visible amplification on a gel after the second PCR.

Returned, demultiplexed sequences were trimmed with cutadapt (v. 2.10) (81) and then filtered and processed using the R (v. 4.0.2) package DADA2 (v. 1.16.0) (82). Reads were trimmed to 216 (16S rRNA) or 260 and 230 (mcrA; forward and reverse reads, respectively) base pairs, with those containing more than 2 expected sequencing errors removed. Amplicon sequence variants (ASVs) were then inferred from filtered reads. No reads were recovered from blank samples (9 16S rRNA blanks; 5 mcrA blanks). Phylogenetic classification of 16S rRNA ASVs was based on the SILVA SSU database (v. 132) (83); on average, 1.96 × 10⁴ ± 1.47 × 10⁴ 16S rRNA reads were recovered per sample. Classification of mcrA ASVs was determined manually based on a phylogenetic tree, created as described below.

Sequences generated with the mcrA_F/mcrA_R primers rather than the mcraMF/mcraMR primers were used for further analysis, since they amplified more samples (Table S2). Additionally, across all samples amplified with both primer sets, 272 more unique ASVs were inferred when the mcrA_F/mcrA_R primers were used. ASVs generated with the mcrA_F/mcrA_R primers were clustered in Geneious Prime (v. 2021.1.1) with a minimum overlap identity of 97% and assigned to 122 clusters. The highest-scoring NCBI BLAST hit was identified for the most abundant ASV in each cluster. These sequences (98 unique), as well as several sequences from cultured methanogens, were aligned with MAFFT (v. 7.450) (84) and incorporated into a maximum-likelihood PhyML (v. 3.3) (85) reference tree with 100 bootstraps. The percent identity between individual ASVs inferred from mcrA_F/mcrA_R samples and the NCBI BLAST hit it was assigned to is shown in Fig. S12.

16S rRNA gene sequences for comparison to the USAM data set were downloaded from the Sequence Read Archive (SRA), project numbers PRJNA275905 (55), PRJNA265122 (56), PRJNA381353 (8), PRJNA386387 (36), PRJNA506542 (58), and PRJDB8874 (57). All projects in database searches for “methane seep,” “MiSeq,” and “16S rRNA” and those that used primers targeting the V4 and/or V5 region were included; samples that came from enrichment cultures or nonsediments were manually eliminated. A full list of sample names, as well as primer and geographic information, is provided in Data Set S3.

Sequences from all studies were trimmed using cutadapt with 515F-Y (78) and 806R (86) primers if applicable. (Samples from PRJNA275905, PRJNA265122, PRJNA381353, and PRJNA386387 were pretrimmed before submission to the SRA and thus not trimmed in house.) Sequences were then processed as described above using the R package DADA2, with the exception that reverse reads were trimmed to 100 bp, rather than 216 bp. Additionally, sequences assigned to bacteria were removed from samples in PRJNA381353 and PRJNA506542, as they were originally sequenced with archaeal primers rather than universal primers.

The correlation analysis between ANME archaea and SRB was completed using the Java package MINE (v. 2.0.1) (52). A table of read abundances was generated for the 4,782 16S rRNA gene and 16S rRNA ASVs within the classes Methanomicrobia and Deltaproteobacteria across all 52 USAM Seep samples. ASVs which did not appear in at least 10% of samples were removed. The maximal information coefficient was calculated for each pair of remaining ASVs, and P values were corrected for multiple comparisons in R by the Benjamini-Hochberg method (87). Of the 849 significant (P < 0.01) correlations, 212 were found to involve an ANME archaeon and a member of the Deltaproteobacteria and were included in our network analysis.

Maximum-likelihood trees of 16S rRNA gene and 16S rRNA ASVs from the USAM and the global seep data set were inferred with FastTree (v. 2.1.11) (88) under a general time reversible (GTR) model of evolution. Nonmetric multidimensional scaling (NMDS) was carried out based on the weighted UniFrac distance metric (89) using the R package vegan (v. 2.5-7) (90). NMDS stress values below 0.2 reliably represent the underlying data. Distance-based redundancy analyses (dbRDA) were also carried out on based on the weighted UniFrac distance metric. The environmental data matrix was z-score transformed, and variance inflation factors of >10, representing covariability, were discarded before stepwise model selection using the ordistep function in the R package vegan (v. 2.5-7) (90). Principal-component analysis (PCA) was carried out using singular value decomposition via the prcomp function in the R package stats (v 4.0.2). Geochemical data were centered and scaled.

Analysis of similarity (ANOSIM) was used to determine the significance of microbial community differences between groups of samples, also based on the weighted UniFrac distance metric. Analysis of variance (ANOVA) was used to test the significance of dbRDA models. One-sided Student’s t tests were used to compare the average weighted UniFrac distance between USAM seep samples with that of 10 equivalently sized collections of samples randomly chosen from the global seep data set. (This was conducted for the total microbial community as well as for the Archaea only).

The phylogenetic normalized stochasticity ratio (pNST) was calculated for Seep sediments and for non-Seep sediments using the R package NST (v. 3.0.6) (54), developed by Ning et al. (54). The index is based on the weighted UniFrac similarity between communities of the same treatment and the randomly expected similarity after randomization of the entire metacommunity for that treatment. For these analyses, ASV-level community data for samples within the same sediment core were pooled, and cores were categorized as “Seep” or “non-Seep” cores (giving an n value of 7 in the Seep treatment and 8 in the non-Seep treatment).

Nucleotide sequences from this study were deposited in the European Nucleotide Archive, project number PRJEB45337. Sample metadata are contained in Data Set S1.