Guaymas Basin sites were visited and sampled with R/V Atlantis, human-occupied vehicle (HOV) Alvin, and autonomous underwater vehicle (AUV) Sentry during cruise AT37-06 (December 6–29, 2016). Alvin dives targeted previously explored sampling areas in the southern axial valley of Guaymas Basin (Teske et al., 2016), or newly identified sites found by AUV Sentry. Hydrothermal sediment samples (Alvin push cores) used in this study were obtained during Alvin Dive 4872 (Dec. 24, 2016) in the Cathedral Hill area (27°00.684 N/ 111°24.266 W) at 2014 m depth. The sample set covers a hydrothermal gradient from hot hydrothermally charged sediments covered with orange and white Beggiatoaceae mats to temperate bare sediment, within ~1 m distance (Figure 1). Core 4872–01 represents bare sediment that was not covered by a Beggiatoaceae mat, whereas Core 4872–06 was covered by a white Beggiatoaceae mat, and Core 4872–14 was covered by an orange Beggiatoaceae mat. Once the sediment cores were returned to the surface, porewater profiles were collected using rhizons (Rhizosphere Research Products, Wageningen, The Netherlands). Then, these cores were divided into 3 cm sections and methane samples were collected by injecting 2 mL sediment plugs into 30 mL serum vials containing 2 mL of 1 M sodium hydroxide solution each. The serum vials were immediately sealed with thick rubber stoppers, crimped with aluminum seals, and stored at 4°C. The remaining sediment samples were frozen at −80°C by the shipboard science crew, and later used for DNA extraction shoreside.

Thermal profiles were measured in surficial sediments using Alvin’s 50 cm heat flow probe.¹ The 50 cm probe contains thermal sensors every 10 cm, starting 5 cm below the attached plastic disk (the “puck”) that limits probe penetration and rests on the seafloor once the probe was inserted. After 5 to 10 min of temperature reading stabilization, temperature readings were recorded.

For sulfide analysis, 1 mL porewater subsamples were fixed with 0.1 mL of 0.1 M zinc acetate solution to preserve sulfide as zinc sulfide until analysis using the methylene blue method (Cline, 1969). For sulfate measurements, porewater samples were diluted in water (1:50) and analyzed using ion chromatography (Metrohm 930 Compact IC flex with Metrosep A PCC HC/4.0 preconcentration column, and Metrosep A Supp 5 Guard/4.0 chromatography column). Sulfate for the white and orange cores was measured by Wetland Biogeochemistry Analytical Services (WBAS) using ion chromatography on a Dionex ICS-1000 using a RFIC IonPac AS22 4 × 250 mm column and a RFIC IonPac AG22 Guard column 4 × 50 mm. The IC was calibrated using a 5-point standard curve and a minimum R² value of 0.999. The detection limit measured at 0.10 ppm sulfate. A check standard was made from a Dionex 7 anion standard (1,000 ppm) to a concentration of 15 ppm. An external Standard from HACH (100 ppm sulfate) was used to make a 20 ppm quality check standard. In general, quality checks were analyzed for 10% of the samples.

Porewaters from the Guaymas Basin were analyzed at WBAS for nitrate/nitrite (NO_x) and ammonium (NH₄⁺) concentrations colorimetrically, using a Flow Solutions IV segmented flow Auto Analyzer from O.I Analytical, College Station, TX. The HCl-acidified samples were neutralized with sodium hydroxide (NaOH) before analysis. NO_x was determined using the cadmium reduction method and NH₄⁺ was determined using the phenate method. Both nutrients were diluted to get their concentrations within the linear range of the auto analyzer. Quality control standards from certified stock standards (Environmental Research Associates), were analyzed every 15–20 samples.

For combined concentration and δ¹³C analysis of methane, 2 mL sediment subsamples were added to 30 mL serum vials containing 2 mL of 1 M NaOH solution, sealed with thick butyl rubber stoppers, crimped with aluminum seals, and stored at 4°C. Since cores were retrieved unpressurized, outgassing may have particularly impacted the measurements of methane concentrations near and above saturation. After the cruise, the methane samples were analyzed by headspace gas chromatography-flame ionization detection (GC-FID) at Florida State University (Magen et al., 2014). Additionally, the gas samples were analyzed for δ¹³CH₄ by injecting 0.1 to 0.5 mL of sample into a gas chromatograph interfaced to a Finnigan MAT Delta S Isotope Ratio Mass Spectrometer inlet system as previously described (Chanton and Liptay 2000). Values are reported in the per mil (‰) notation relative to Vienna Pee Dee Belemnite (VPDB).

DNA was extracted from Guaymas Basin sediment samples using a published protocol for soil DNA (Zhou et al., 1996). The extraction method was modified to optimize DNA recovery from deep-sea sediments. The detailed protocol can be found in the Supplementary Data Sheet 1.

Sequencing libraries were prepared from extracted genomic DNA with the commercially available LoopSeq kits from Loop Genomics (Supplementary Data Sheet 2).² Synthetic long reads (SLRs) were constructed from the short-read sequencing reads using the standard Loop Genomics informatics pipeline and processed as previously described (Callahan et al., 2021). The method involves attaching two DNA tags: one Unique Molecular Identifier (UMI) to each unique “parent” DNA molecule and one sample-specific tag (i.e., a sample index) equally to all molecules in the same sample. Barcoded molecules are amplified, multiplexed, and each UMI is distributed intramolecularly to a random position within each parent molecule. Molecules are then fragmented into smaller units at the position of each UMI, creating a library of UMI-tagged fragments with an average length of 400 bp compatible with an Illumina sequencing platform run in PE150 mode (Callahan et al., 2021).

This protocol allows for the assembly of continuous 16S rRNA gene long reads from individual DNA molecules, providing sequencing data like classical long read sequencing methods but with much higher throughput. The terminal primers determine the extent of the 16S rRNA gene sequence yielded by this approach (27F, 5’-AGAGTTTGATCMTGGCTCAG-3′; 1492R, 5’-TACCTTGTTACGACTT-3′). LoopSeq can be performed on existing DNA sequencing instruments, and due to the chemistry and the tagging of individual molecules the method is semi-quantitative providing a major advantage toward other amplification-based methods.

We note that LoopSeq reconstructs a single long read from several short reads. Therefore, we define the resulting amplicon sequence variants (ASVs) as LoopSeq variants (L-ASVs) to differentiate them from “traditional” ASVs that are essentially PCR amplicons.

The statistical programming language R 4.2.1 (R Core Team, 2022) and the DADA2 (v. 1.20.0) bioinformatics package (Callahan et al., 2016), along with a LoopSeq suitable DADA2 workflow (Callahan et al., 2021) were used to analyze 16S rRNA gene sequences obtained by LoopSeq sequencing. This workflow is available at Github.³ DADA2 creates ASVs (or LoopSeq variants (L-ASVs) in our case) by grouping together amplicon reads of the same sequence.

The function “removePrimers” was used to select for 16S rRNA amplicon sequences with both the forward and reverse primers and allowing a maximum mismatch of two base pairs, as reads lacking either primer do not represent the complete gene. The following parameters were used to filter the remaining sequences: a minimum length of 1,000 bp (minLen = 1,000), a maximum length of 2000 bp (maxLen = 2000), a minimum quality score of 3 (minQ = 3), a maximum expected error rate of 1 (maxEE = 1), and a maximum Ns of 0 (maxN = 0). These parameters may be used to filter reads without selecting for reads with a forward and reverse primer to capture long, but not complete genes. We found that 20.4% of reads that fulfilled the previously stated parameters lacked the forward/reverse primer. Only full-length gene sequences were kept in our analysis. The function “dereqFastq” was used to dereplicate reads. The function “learnErrors” was used with a band size of 32 (BAND_SIZE = 32) to learn error rates. The error estimate function was set to the PacBio Error Function. DADA2 was ran with a p-value cutoff of 1⁻¹⁰ (OMEGA_A = 1e-10), single detection (DETECT_SINGLETONS = TRUE), and a band size of 32 (BAND_SIZE = 32). The function “removeBimeraDenovo” (method = consensus) was used to remove chimeras. An overview of raw read processing is depicted in Supplementary Table 1.

The function “assignTaxonomy” was used to assign taxonomy using the Silva reference database (v138) (Quast et al., 2013). A total of 12,494 L-ASVs were identified. Figures depicting taxonomic distribution, alpha diversity, and beta diversity (Figure 2; Supplementary Figures 2, 3) were created using the R Phyloseq package (McMurdie and Holmes, 2013). All sequence reads were submitted to the NCBI Sequence Read Archive (BioProject: PRJNA1105367).

We compared our L-ASV sequences against reference sequences obtained from the NCBI GenBank database. All sequences were aligned using Multiple Sequence Comparison by Log-Expectation (MUSCLE) in MEGA11 (Stecher et al., 2020). From these alignments, MEGA11 was used to construct minimum-evolution phylogenetic trees. MEGA11 provides the option to use the Kimura 2-parameter (K2P) model to calculate gene distances (Kimura, 1980). The K2P model is well-suited to clustering similar sequences, and stays close to empirically observed 16S rRNA gene distances (esp. for low sequence divergence <5%, common in our datasets) that form the basis for taxonomic recommendations (Stackebrandt and Goebel, 1994; Qin et al., 2014; Konstantinidis et al., 2017). The resulting phylogenies were checked by the aggregate of 1,000 bootstrap replicates (Rzhetsky and Nei, 1995; Felsenstein, 1985), in addition to cross-checking with published phylogenies for consistency. Bootstrap values greater than 70% are shown and NCBI GenBank accession numbers are provided on tip labels. Phylogenies were checked for consistency against independently derived phylogenies in relevant papers (Knittel et al., 2003; Salman et al., 2011; Schutte et al., 2018; Jochum et al., 2018; Krukenberg et al., 2016).

Distance matrices corresponding to each phylogenetic tree are shown in Supplementary Data Sheet 3. The number of base substitutions between sequences are shown as a percentage (i.e., per 100 sites). Analyses were conducted using the Maximum Composite Likelihood model (Tamura et al., 2004). Each analysis involved all 16S rRNA gene sequences displayed in the corresponding phylogenetic tree. All ambiguous positions were removed for each sequence pair (pairwise deletion option). Evolutionary analyses were conducted in MEGA11 (Tamura et al., 2021; Stecher et al., 2020).

The color of Guaymas Basin Beggiatoaceae mats indicated the relative temperature of the underlying sediment as measured using the Alvin heat flow probe; orange mats covered hotter sediments and white mats covered comparatively cooler sediments (McKay et al., 2012). Temperature varied greatly with increasing sediment depth at the three sample sites (Figures 1, 3A). Temperatures at 40 centimeters below the seafloor (cmbsf) exceeded 100°C below the orange mat, reached 75°C below the white mat, and remained near 20°C in sediments without a Beggiatoaceae mat (Figure 3A).

Temperature considerably impacts microbial population density, and therefore the amount of DNA available for extraction. Measurable concentrations of DNA persisted the deepest (40 cmbsf) in Core 1, from bare sediment neighboring the warm mats. In contrast, DNA recovery is limited to surficial sediments (~ 20 cmbsf) in Core 14 from the hot center of the orange mat (Figure 3B.). Interestingly, the highest DNA yield (ca. 2.5 to 7 ug DNA/cm³ sediment) are found within the upper ~6 cm of sediment in all three cores, indicating that the impact of hydrothermal activity extends beyond the Beggiatoaceae mats to the bare sediments of Core 1 (Figure 3B). Previous research in Guaymas Basin consistently reveals a very rapid decline in cell counts and 16S rRNA gene abundances with increasing depth and temperature in hot hydrothermal sediments (Meyer et al., 2013; Møller et al., 2018; Lagostina et al., 2021).

Micromolar concentrations of sulfide are observed throughout the bare sediment core, whereas the white and orange mat cores reach concentrations of sulfide up to 2 millimolar (Figure 3C). At the surface, sulfide concentrations are near zero for the bare sediment and white mat cores. In contrast, the orange mat core has elevated (~1 millimolar) sulfide levels at the surface, consistent with high hydrothermal flow in the mat center. Sulfate is abundant in all three cores, with depletion only occurring at ~32 cm and below in the orange mat core, indicating that hydrothermal circulation contributes to the replenishment of sulfate in surficial sediments (Figure 3D).

Ammonium is observed throughout our white and orange mat sediment cores (Figure 3E) Consistent with previous research, ammonium concentrations increase as temperature increases downcore (Ramírez et al., 2021). Ammonium in Guaymas Basin hydrothermal sediments is primarily derived from the thermal degradation of organic matter (Von Damm et al., 1985). Deep in the orange mat core, ammonium concentrations reach ~8 mM (Figure 3E), close to the ~10–15 mM concentrations observed in Guaymas Basin hydrothermal fluids (Von Damm et al., 1985).

Elevated nitrate/nitrite concentrations are observed in surficial mat covered sediments (Figure 3F). Bottom water concentrations of nitrate in Guaymas Basin are up to 50 μM (Winkel et al., 2014). Elevated nitrate levels (up to nearly 0.25 mM; Figure 3F) in mat-covered surficial sediments are at least in part due to damaged Beggiatoaceae mats leaking nitrate stored in their vacuoles into surficial sediments (Schutte et al., 2018).

High concentrations of methane are observed in all three cores (Figure 3G). A gradient of peak methane concentrations (~1 mM in the bare sediment core, ~2 mM in the white mat core, and ~ 8 mM in the orange mat core) illustrates the impact of hydrothermal fluids that are supersaturated in methane. High methane concentrations (multiple millimolar) are a characteristic feature of Guaymas Basin hydrothermal sediments, similar to those found at cold seep sites (Welhan, 1988; McKay et al., 2016; Song et al., 2021).

Hydrothermal methane in Guaymas Basin unmodified by microbial activity has a δ¹³CH₄ value of approximately −42 ‰ (McKay et al., 2016; Song et al., 2021). The preferential oxidation of ¹²CH₄ by methane-oxidizing archaea leads to the orange and white mat cores containing methane with a heavier carbon isotopic signature that is relatively enriched in ¹³C, with values ranging from −30 to −20 ‰ (Figure 3H). The lighter carbon isotopic signature of methane throughout the bare sediment core and near the surface of the orange and white mat cores, indicates some contribution due to microbial methanogenesis (Figure 3H). The methanogens in these sediments include methylotrophs, which thrive in sulfate-replete sediments (Zhuang et al., 2018); numerous additional types of methanogens have been detected as well (Dhillon et al., 2005; Lever and Teske, 2015).

LoopSeq recovered nearly full-length sequences (Supplementary Figure 1, in Data Sheet 4) that affiliated with several phyla of heterotrophic bacteria, such as most Proteobacteria, Desulfobacterota, Bacteroidota, and Chloroflexi (Figure 2). These phyla have been detected by previous Guaymas Basin sequencing surveys using clone libraries of near full-length 16S rRNA genes (Teske et al., 2002; Dowell et al., 2016; McKay et al., 2016).

We analyzed the alpha diversity of each mat site using five measures (Observed, Chao1, Shannon, Simpson, and Inverse Simpson). We expected hydrothermal stress to influence alpha diversity based on previous research showing that alpha diversity decreases with increasing hydrothermal fluid temperature (Barosa et al., 2023) and with increasing sediment depth and temperature (Ramírez et al., 2021). However, we found no significant difference in alpha diversity between the three coring sites as indicated by the Kruskal-Wallis rank sum test (Supplementary Figure 2, Data Sheet 4). This unexpected result may reflect the previously observed high degree of habitat connectivity in Guaymas Basin sediments (Meyer et al., 2013), emerging from the coexistence of different electron acceptors and donors in non-consolidated sediments that undergo rapid hydrothermal pulses (McKay et al., 2016) and allow the coexistence of physiologically distinct microbial groups in surficial sediments (Engelen et al., 2021).

To explore potential differences in beta diversity, we created non-metric multi-dimensional scaling (NMDS) and dendrogram plots based on the Bray–Curtis dissimilarity method. Based on the dendrogram we highlighted a cluster of deeper samples (≥ 21 cmbsf) from the bare sediment site (Supplementary Figure 3, Data Sheet 4).

Sulfur-cycling lineages that frequently interface with hydrocarbon metabolism are present throughout surficial hydrothermal sediments in the Guaymas Basin (Figure 2). Numerous sulfate-reducing lineages (SEEP-SRB lineages, Candidatus Desulfofervidus) exist in syntrophic assemblages with methane- and alkane-oxidizing archaea (Knittel et al., 2003; Krukenberg et al., 2016), while other sulfate reducers (Desulfatiglans) perform terminal oxidation of hydrocarbons as free-living organisms (Teske et al., 2019). Sulfur-oxidizing Beggiatoaceae are known to assimilate methane-derived carbon (McKay et al., 2012). We performed phylogenetic analyses on these sulfur-cycling lineages to determine their intragroup diversity. Of 12,494 L-ASVs, 799 were identified as Candidatus Desulfofervidus, 585 as Desulfatiglans, 158 as SEEP-SRB2, 135 as SEEP-SRB4, and 29 as Beggiatoaceae.

The detection of Beggiatoaceae is an interesting test of LoopSeq capabilities. Previous efforts to generate full-length Beggiatoaceae sequences required collecting large amounts of individually prepared and cleaned filaments (McKay et al., 2012). Due to the size of individual Beggiatoaceae cells, they are outnumbered by other much smaller marine bacteria, complicating sequencing efforts and requiring high-throughput methods to overcome the ubiquitous presence of other marine bacterial groups. In our survey, all LoopSeq variants (L-ASV) identified as Beggiatoaceae matched reference sequences of the white or orange phylotype, with no additional phylotypes identified (Figure 4; Supplementary Figure 5, in Data Sheet 4). We note that light microscopy observations have revealed substantial diversity in the filament sizes of Beggiatoaceae (Nelson et al., 1989, Salman et al., 2013), indicative of further sequence diversity that remains to be detected. There is no evidence for introns within the 16S rRNA genes of Guaymas Basin Beggiatoaceae, although they have been observed in non-Guaymas Basin Beggiatoaceae (Salman et al., 2011; Salman et al., 2012). The LoopSeq primers coincide positionally with PCR primers previously used to amplify Beggiatoaceae sequences (Teske et al., 2002; McKay et al., 2012), and in addition accommodate ambiguous nucleotide positions; therefore, primer mismatches can be ruled out. Recovered L-ASVs were generally found to be associated with surficial sediments underlying white and orange Beggiatoaceae mats. Some L-ASVs were found below the surface, resulting from the coring device pulling some of the long filaments down into the sediment column (Teske, 2009; Supplementary Figure 6 in Data Sheet 4).

In contrast to Beggiatoaceae, Candidatus Desulfofervidus sequences are often recovered in Guaymas Basin sequencing surveys (McKay et al., 2016; Dowell et al., 2016). Ca. Desulfofervidus is a sulfate-reducing syntrophic partner of methane- and alkane-oxidizing archaea (Laso-Pérez et al., 2016) with a thermal optimum of 50°C (Holler et al., 2011) to 60°C (Kniemeyer et al., 2007; Wegener et al., 2015). Although often observed as a syntroph, Ca. Desulfofervidus can also survive as a free-living hydrogenotroph without methane- or alkane-oxidizing partners (Krukenberg et al., 2016). Ca. Desulfofervidus is a phylogenetically narrow clade, as it is classified as its own family, order, and class-level lineage within the phylum Desulfobacterota (Waite et al., 2020). However, previous research has not examined the phylogenetic diversity within Ca. Desulfofervidus.

Our phylogeny identified several phylotypes of near-identical 16S rRNA gene sequences within Candidatus Desulfofervidus, each consisting of reference sequences from 50°C (Holler et al., 2011) or 60°C (Kniemeyer et al., 2007; Wegener et al., 2015) enrichments. All reference sequences were from Guaymas Basin, with no close matches from other sites. The close alignment of our sequences with established Guaymas Basin reference sequences attests to their quality. Our phylogeny supports the previous assertion of Krukenberg et al. (2016) that all Ca. Desulfofervidus lineages are closely related, since they form a well-defined 16S rRNA gene sequence cluster of approximately 98.5% sequence identity (Supplementary Data File 1, Sheet 2). The repeated detection of near-identical phylotypes within the Ca. Desulfofervidus cluster (Figure 5) implies the presence of unrecognized species diversity within Ca. Desulfofervidus.

We found sequences in GenBank that formed two deeply branching clusters with strong bootstrap support that are distinct from the major group of Ca. Desulfofervidus (Figure 5). These clusters were separated from the main group by sequence divergence up to 5.8% (Supplementary Data File 1, Sheet 2). The recovery of these sequences across independent studies at different vent sites suggests unrecognized microbial diversity related to Ca. Desulfofervidus. A previous metagenomic survey of Guaymas Basin hydrothermal sediments detected metagenome-assembles genomes (MAGs) branching between Ca. Desulfofervidus and Thermodesulfobacteriaceae (Dombrowski et al., 2018). Detailed analyses of these MAGs (labeled as B74_G16 and B4_G9) would reveal the genomic content and physiological potential of Desulfofervidus-affiliated lineages.

Like Ca. Desulfofervidus, SEEP-SRB2 sequences are often recovered in Guaymas Basin sequencing surveys (Dowell et al., 2016). Our survey recovered sequences from slightly below the surface to approximately 15 cmbsf. Unlike Ca. Desulfofervidus, SEEP-SRB2 sequences are recovered from diverse non-hydrothermal environments; cold seeps, salt lakes, coastal sediments, and estuaries (Kleindienst et al., 2014; Ruff et al., 2015). Despite its ecological importance as a syntrophic partner of methane and alkane oxidizers at cool and warm temperatures (enriched at 20°C and 37°C; Krukenberg et al., 2018), the internal phylogenetic structure of this group has not been examined in any detail.

Our SEEP-SRB2 phylogeny is largely comprised of a well-supported major cluster with ~97.5% sequence identity (Supplementary Data File 1, Sheet 3, equivalent to Genus level; Schloss and Handelsman, 2005). The closest cultured relatives of SEEP-SRB2, Dissulfuribacter thermophilus and Dissulfurimicrobium hydrothermale, are thermophilic sulfur disproportionators isolated from hydrothermal environments (Slobodkin et al., 2013; Slobodkin et al., 2016). These bacteria constitute separate lineages but share a deep branching point with SEEP-SRB2 (Figure 6). The distinctiveness of these lineages is further supported by sequence distances of up to 11% (Supplementary Data File 1, Sheet 3). Our 16S rRNA gene sequences indicate that SEEP-SRB2 would form a stand-alone group within Dissulfuribacteraceae (Waite et al., 2020). Further phylogenetic diversity that would fall between the pure culture examples (Dissulfuribacter thermophilus and Dissulfurimicrobium hydrothermale) and SEEP-SRB2 remains to be explored. Previously, the creation of a taxonomic family to comprise Dissulfuribacter thermophilus and Dissulfurimicrobium hydrothermale has been suggested (Slobodkin et al., 2013). Our phylogeny shows that such a group would also include SEEP-SRB2 (Figure 6).

In contrast to the syntrophic symbionts SEEP-SRB1 and SEEP-SRB2, SEEP-SRB4 has largely been overlooked since it was proposed 20 years ago (Knittel et al., 2003). FISH surveys indicate that SEEP-SRB4 is not a syntroph but occurs as individual free-living cells in cold seep sediments (Knittel et al., 2003). The closest cultured relatives of SEEP-SRB4 are non-thermophilic sulfur-disproportionating genera such as Desulforhopalus (Isaksen and Teske, 1996), Desulfocapsa (Finster et al., 1998), Desulfofustis (Friedrich et al., 1996), and the sulfate-reducing psychrophile Desulfotalea (Knoblauch et al., 1999). While SEEP-SRB2 has been enriched from Guaymas Basin sediments at 37°C (Krukenberg et al., 2018), there is currently no evidence that SEEP-SRB4 is moderately thermotolerant. We seek to test this assumption in our survey, checking the distribution of SEEP-SRB4 in surficial hydrothermal sediments. As with SEEP-SRB2, we seek to investigate the previously unrecognized phylogenetic diversity of SEEP-SRB4.

Our SEEP-SRB4 phylogeny depicts one major cluster with strong bootstrap support (100%) and > 98.9% sequence identity (Supplementary Data File 1, Sheet 4). This cluster contains several well-supported, phylogenetically narrow clusters (Figure 7). While the phylogeny shows a potential Guaymas Basin-specific subcluster of four sequences, the other subclusters contain globally distributed sequences from diverse cold seep sites and from Guaymas Basin.

The majority of L-ASVs for SEEP-SRB4 are retrieved from the top ~6 cm of the sediment across the entire lateral hydrothermal gradient (Supplementary Table 3). The preference for surficial sediments may be due to the abundance of elemental sulfur and polysulfides located there (Teske et al., 2016). A similar preference for surficial sediments with coexisting sulfate and sulfide was noted in the first study of SEEP-SRB4 (Knittel et al., 2003). SEEP-SRB4 bacteria seem to be able to colonize surficial hydrothermal sediments as long as in-situ temperatures remain relatively cool, similar to the strategy utilized by Guaymas Basin Beggiatoaceae (McKay et al., 2012).

Desulfatiglandales is an order-level lineage within the phylum Desulfobacterota that contains polyaromatic degrading organisms, both isolated and uncultured species. The order and family level taxonomy are named after Desulfatiglans annilini, the first cultured representative (Suzuki et al., 2014). Sequences related to Desulfatiglans are often recovered in Guaymas Basin sequencing surveys (Ramírez et al., 2021; Dowell et al., 2016), including the shallow subsurface (Edgcomb et al., 2022). In contrast to the tightly clustered groups of SEEP-SRB2 and SEEP-SRB4, Desulfatiglans exhibits a greater degree of phylogenetic diversity.

Our Guaymas Basin sequences and almost all GenBank sequences were not closely related to cultured isolates that degrade polyaromatic and substituted-aromatic compounds (Figure 8). Even finding uncultured sequences that were closely related (>98% sequence identity) to some of our Guaymas Basin ASVs proved to be difficult. This could imply an undersampling of the natural diversity of Desulfatiglans, or their ongoing diversification and speciation. Previously, Desulfatiglans has been classified into two groups distinguished by genomic and physiological differences, Group 1 consisting of sulfate reducers and Group 2 consisting of dehalogenating subsurface bacteria (Jochum et al., 2018). Most of our L-ASVs belong to Group 1 and a few belong to Group 2.

Outside of these two groups, we found two well-supported hydrothermal clusters (Figure 8). The first cluster contained L-ASV 190 (predominately recovered from the white mat site), sequences from 37°C and 50°C enrichments of Guaymas Basin sediment (Kellermann et al., 2012; Wegener et al., 2015) and an environmental sequence from Guaymas Basin hydrothermal sediment (McKay et al., 2016). The second cluster contained L-ASV 207 (predominately recovered from the orange mat site), a sequence from a Desulfatiglans MAG that was prevalent in a 70°C enrichment of Guaymas Basin sediment (Zehnle et al., 2023), and environmental sequences from diverse hydrothermal sediments (Li et al., 2014; Li et al., 2020). We examined these sequences via distance matrix analyses and find that these hydrothermal clusters are generally separated from other sequences by 11–12% sequence divergence (Supplementary Data File 1, Sheet 5).

In contrast to other sulfate reducers that occur in warm and hot Guaymas Basin sediments, we observed most of our Desulfatiglans L-ASVs in cool sediments, suggesting that they do not depend on methane and light hydrocarbons derived from hydrothermal fluids, and instead utilize polyaromatics (Schnell and Schink, 1991; Zehnle et al., 2023). Polyaromatics from hydrothermal oils are prevalent throughout Guaymas Basin sediments of diverse thermal regimes (Kawka and Simoneit, 1990), and Desulfatiglans appears to have evolved to utilize polyaromatics across the thermal spectrum, as suggested by its hydrothermal clusters.

DADA2 detected chimeras in approximately 4% of our reads from a natural hydrothermal community. Previous research found structural errors (chimeras and introgressions) in approximately 2% of reads from a Zymo mock community (Callahan et al., 2021). Specifically, we found that LoopSeq produced chimeras in the majority of Candidatus Desulfofervidus sequences we examined. Segments of these sequences were entered into BLASTn, and we discovered that the chimeras consistently occurred in the first ~200 bp. BLASTn searches of the chimeric segments revealed their relation to other marine sediment bacteria such as Chloroflexi. The chimera breaks are located within the highly variable Helix 11 of the 16S rRNA gene (Van de Peer et al., 1996), and therefore appear to be secondary structure dependent. These issues with Ca. Desulfofervidus sequences were observed in samples across the DNA concentration range, from very high (4,240 ng/cm³) to very low (0.5 ng/cm³). Automated chimera identifying programs can have sensitivity issues when used with uncultured lineages that have poorly resolved phylogenies. We found that trimming off the first ~200 bp remedied the chimera issue, resulting in sequences with high alignment quality and > 97% identity to previously archived Ca. Desulfofervidus sequences available in GenBank.

LoopSeq recovered eight sequences identified by the DADA2 annotation pipeline as “Thermodesulfobacterium.” Initial analysis of theses sequences showed that they were plagued with chimera and artifact issues, a potential consequence of low DNA yields from deep, hot samples. Two sequences were salvaged; one did not have chimera issues (L-ASV 12240), while the other was trimmed to only include the first ~340 bp (L-ASV 11484). These two sequences came from samples with low DNA concentrations (Sample 4872-06-08, 50 ng/cm³ and Sample 4872-06-09, 12 ng/cm³). These two sequences formed a Guaymas Basin specific cluster (Supplementary Figure 7) with the thermophilic ANME-1 syntroph Thermodesulfobacterium torris (Benito Merino et al., 2022) and a nearly identical hydrothermal sediment clone (McKay et al., 2016).