Geochemical context of hot springs
MV2 (also referred to as “Figure 8”; 44°36′36.14″, −110°26′21.95″) is an acidic spring in the Greater Obsidian Pool Area of the Mud Volcano geyser basin of eastern-central YNP (Supplementary Fig. S1). At the time of initial sampling for metagenomics analysis in 2014, the temperature of MV2 was 62.0 °C and the pH was 3.8. Over the next 4 years of repeated sampling, the geochemistry of MV2 fluctuated (Supplementary Table S1). For example, the pH of MV2 waters ranged from 3.0 to 3.8 and the temperature ranged from 59 to 69 °C (Supplementary Table S1). Published data suggests even wider geochemical variance, with pH ranging as low as 3.0 and as high as 4.8 [25,26,27]. The relatively high levels of SO42− and Cl− within MV2 waters are indicative of an acid-sulfate chloride spring (Supplementary Fig. S3), sourced by a mixture of a deep hydrothermal water end member enriched in Cl− and a vapor-phase fluid end member enriched in SO42− via oxidation of sulfide (Supplementary Fig. S3) [43].
SJ3 spring is in the Smokejumper geyser basin within southwestern YNP (44°24′57.42″, −110°57′20.76″, Supplementary Fig. S1) [28]. At the time of sampling for metagenomics analysis in 2014, SJ3 waters exhibited a moderately acidic pH of 5.4 and a temperature of 61.9 °C. The geochemical composition of SJ3 (moderate SO42− levels and extremely low Cl−) is consistent with mixing of a vapor-phase fluid end member with surface-derived meteoric waters [29, 43] (Supplementary Table S1, Supplementary Fig. S3).
Recovery of euryarchaeal MAGs from YNP hot spring metagenomes
Metagenomic characterization of the MV2 and SJ3 sediment communities from 2014 yielded five (Supplementary Text; Supplementary Table S2; Supplementary Fig. S4) and 82 [28] moderate to high-quality MAGs, respectively. A high-quality MAG estimated to be nearly complete (98% estimated completeness, Supplementary Table S3) was recovered from the MV2 community metagenome and contained a near full-length 16S rRNA gene (length = 1425 nt) that was 100% identical to those from the “Thermoplasmatales A10” group (Euryarchaeota phylum) that have been previously identified in MV2 sediments [25, 26] (Supplementary Figs. S5 and S6). This MAG, hereafter referred to as MV2-Eury, comprised 1.30 Mbp and 1384 predicted PCGs. An additional MAG recovered from the SJ3 hot spring community metagenome (hereafter SJ3-Eury) was represented by an unbinned 16S rRNA gene (length = 994 nt) that was 100% identical to that of MV2-Eury (Supplementary Figs. S5 and S6). Further, pairwise average amino acid and nucleotide identities (AAI/ANI, respectively) were nearly identical among the two Euryarchaeal MAGs (two-way AAI: 99.36% ± 1.50%, n = 1201 encoded proteins; two-way ANI: 99.08% ± 0.78%, n = 5599 alignments), suggesting a very high degree of genomic similarity, although unique subsets of proteins were present in each MAG. In addition, the two euryarchaeal MAGs also exhibited highly similar G + C contents (see Supplementary Table S3 for additional MAG statistics). The SJ3-Eury MAG was also estimated to be nearly complete (98% estimated completeness) but was smaller than that of MV2, comprising a length of 1.19 Mbp and 1254 predicted PCGs (Supplementary Table S3). The MAGs were estimated to have <1.0% redundancy (i.e., contamination) in single-copy marker genes (Supplementary Table S2), thus meeting the criteria of high-quality MAG drafts [44]. Lastly, a second MV2-Eury MAG was generated from sediments collected in 2018 that was effectively identical to the MV2-Eury MAG (two-way ANI: 99.85% ± 0.28%; unbinned 16S rRNA gene 100% identical to the 2014 MV2-Eury MAG), and is thus not discussed further.
Phylogenomic analyses distinguished the MV2/SJ3 MAGs from other higher-order clades within the recently proposed euryarchaeal superclass Diaforarchaea [45] (previously referred to as Thermoplasmata; Fig. 1a, b) and better resolved the branching of nodes within the Diaforarchaea compared to the 16S rRNA gene analyses (Supplementary Figs. S5 and S6). Pairwise comparisons among other members of the major Diaforarchaea groups indicated that the MV2-/SJ3-Eury MAGs shared <60% mean AAI with all members of the Diaforarchaea and <50% AAI with members of the most closely-related phylogenetic clade (Fig. 1b), the Aciduliprofundum genus (previously known as the Deep Sea Hydrothermal Vent Euryarchaeota group 2 [DHVE2] Archaea). Such values are consistent with a class or order-level distinction from the Aciduliprofundum based on previously published suggestions [46]. Consequently, the analyses indicated that the MV2- and SJ3-Eury MAGs represented members of a previously uncharacterized group of Diaforarchaea that are distantly related to the Aciduliprofundum spp., which together comprise an outgroup to acidophilic Thermoplasmatales-related Archaea (Fig. 1a, b). The moderately acidophilic, S0/ferric iron reducing DHVE2 Archaea were first isolated from deep sea hydrothermal vents [47], and DHVE2 cultivars or MAGs are thus far only known from such environments. This biogeographic distinction suggests that the MV2/SJ3-Eury clade represents a terrestrial hydrothermal lineage that is ecologically distinct from currently characterized DHVE2. In support of this assertion, all publicly available 16S rRNA gene phylotypes associated with the MV2/SJ3-Eury phylotypes were recovered from terrestrial hydrothermal springs/pools (Supplementary Figs. S5 and S6), with no close representatives from marine vents. Specifically, 16S rRNA gene phylotypes closely related to this lineage have also been detected in acidic/moderately acidic and moderate temperature hot springs in geothermal fields in YNP, Japan, Kamchatka, Taiwan, and elsewhere (Supplementary Fig. S6), suggesting that the group is widely distributed among terrestrial hydrothermal systems with these characteristics. The clade is thus referred to hereafter as the YNP terrestrial euryarchaeotal group (YNP-TEG).
a Maximum-likelihood (ML) phylogenomic reconstruction of the MV2/SJ3-Eury MAGs and other Diaforarchaea. The analysis comprised 103 housekeeping protein homologs conserved across Archaea. Branch length is relative to the scale provided at the top indicating the expected number of substitutions per site. Bootstrap values >90 for nodes are indicated by black circles while those between 50 and 90 are denoted by white circles. Values are not depicted for nodes with bootstrap support <50 (out of 100 ML bootstraps). Archaeoglobus fulgidus and Methanobacterium formicicum were used as the outgroups (not shown). Clade-level triangles indicate the phylogenetic diversity within each group via side lengths that are proportional to the distances between the clade’s closest and furthest taxa. b Heatmap showing the pairwise average amino acid identities (AAI) between members of the major groups of Diaforarchaea that were used in the phylogenomic analysis in a. Groups are indicated by colored lines based on the legend at the bottom right. Pairwise AAI values are colored according to the scale on the right. The dendrogram was constructed based on hierarchical clustering of the pairwise average AAI values and the branch length is proportional to the scale at the bottom of the dendrogram.
A new taxonomic clade of SRO
Reconstruction of metabolic pathways involved in energy conservation for the YNP-TEG group indicated that the organisms are putatively capable of SO42−/SO32− reduction coupled to hydrogen and/or organic acid (e.g., lactate, acetate, and formate) oxidation (Fig. 2, Supplementary Text, Supplementary Dataset 1). Both MV2 and SJ3 MAGs contained adjacent dsrA and dsrB protein encoding genes that were distantly related to other DsrAB. Sequence alignment of DsrA confirmed the presence of characteristic siroheme binding (CX5CXnCX3C; Supplementary Fig. S7) and [Fe4S4] cluster binding (CX2CX2C; Supplementary Fig. S7) motifs [48]. DsrAB are catalytically active in both the forward (reductive) and reverse (oxidative) directions [49]. The former is used during SO42−/SO32− reduction in the six bacterial and archaeal lineages with characterized SRO: the Archaeoglobales, Deltaproteobacteria, Firmicutes, Crenarchaeota, Nitrospirae, and Thermodesulfobacteria [50]. In contrast, the latter is used by sulfur oxidizing organisms primarily affiliated with anoxygenic phototrophs (i.e., Chlorobi) and the Proteobacteria [50] and may also be used for sulfur disproportionation in D. alkaliphilus [51].
Question marks indicate unknown metabolic transformations, while gray lines indicate hypothesized metabolic pathways. Abbreviations are as follows: Dsr dissimilatory (bi)sulfite reductase, MQ menaquinol, Nuo archaeal-like ferredoxin-dependent Nuo complex, Hdr heterodisulfide reductase, Mvh F420 non-reducing [NiFe] hydrogenase (Group 3c type), LacZ beta-galactosidase, EM Embden–Meyerhoff glycolysis, Dld/Lld D/L-lactate dehydrogenase, Fdh formate dehydrogenase, FtfL formate-tetrahydrofolate ligase, FolD methylenetetrahydrofolate dehydrogenase, MTHFD methylenetetrahydrofolate dehydrogenase, GlyA glycine hydroxymethyltransferase, Ser D-3-phosphoglycerate dehydrogenase, Apg 2,3-bisphosphoglycerate-independent phosphoglycerate mutase, Fd ferredoxin, Sdh succinate dehydrogenase, Por pyruvate ferredoxin oxidorectase, Acs acetyl-CoA synthase, Scs succinyl-CoA synthetase, Aor aldehyde ferredoxin oxidoreductase, Ior indolepyruvate ferredoxin oxidoreductase, Vor 2-ketoisovalerate ferredoxin reductase, Kgor 2-ketoglutarate ferredoxin oxidoreductase, Fdo formate dehydrogenase, Hyh NADP(H)-coupled bidirectional [NiFe] hydrogenase (Group 3b type).
Three primary phylogenetically distinct clusters of DsrAB have been previously identified, and these coherently group proteins with physiological biases for DsrAB directionality. These include the (1) bacterial (and Archaeoglobales) reductive DsrAB, (2) bacterial oxidative DsrAB, and (3) crenarchaeal reductive DsrAB (nomenclature following [49]). It should, however, be noted that a recent characterization of the elemental sulfur (S0)/polysulfide disproportionating D. alkaliphilus deltaproteobacterium revealed the capacity for S0/polysulfide oxidation or disproportionation that is inferred to be catalyzed by a bacterial reductive type DsrAB, suggesting that the directionality of DsrAB based on phylogenetic clustering may not always be coherent [51]. Nevertheless, the majority of characterized SRO exhibit phylogenetic placements of DsrAB homologs that are consistent with their directionality, and it is not yet clear if the above example is unique to D. alkaliphilus or S0/polysulfide disproportionators, or if it is more widespread. In addition, the recent suggestion of putative SRO or S0 oxidizing lineages from uncultured candidate division MAGs based on bioinformatics data could potentially further confound the inference of DsrAB catalysis from sequence-based homology [50]. Physiological characterization of these uncultured taxa and their capacity for SO42–/SO32− reduction or S0 oxidation will be necessary to fully discern whether they validate the use of DsrAB phylogenies to support catalytic directionality.
An unrooted phylogenetic analysis of concatenated YNP-TEG DsrAB homologs (n = 1218, comprising homologs from previously described datasets [49, 50]) indicated that they form an outgroup to reductive crenarchaeal DsrAB (Fig. 3a, Supplementary Fig. S8). Little is known of crenarchaeal DsrAB and crenarchaeal SRO, in general, with only four demonstrated SO42/SO32− reducing genera in the phylum, all within the Thermoproteales order. SO42− reduction has been suggested or observed in the Vulcanisaeta, Caldivirga, and Thermoproteus [52,53,54], while SO32− reduction (but not SO42− reduction) has been observed in Pyrobaculum [55]. Further, SO42− reduction has been suggested in Thaumarchaeota MAGs recovered from YNP springs, wherein the thaumarchaeal DsrAB show very high homology to those from Vulcanisaeta [32]. Thus, phylogenetic relatedness of the DsrAB homologs among these strains and the YNP-TEG MAGs (albeit with amino acid identities <45%) suggests that the latter are involved in SO42− or SO32− reduction. The unrooted phylogenetic analyses also suggested an early DsrAB evolutionary trifurcation, with one branch comprising the YNP-TEG/crenarchaeal DsrAB homologs. The second branch comprised duplicate DsrAB copies within Moorella spp. genomes and uncultivated organisms within the candidate division (c.d.) “Rokubacteria”, c.d. “Hydrothermarchaeota”, and Verrucomicrobia (Fig. 3a; Supplementary Fig. S8, Supplementary Dataset 2), which have all previously been shown to branch deeply in DsrAB phylogenies [50]. The third branch comprised the canonical reductive and oxidative type DsrAB, in addition to deeply-branching archaeal homologs in the “Aigarchaeota” and “Hydrothermarchaeota” candidate divisions (c.d.).
a Unrooted phylogenetic reconstruction of a concatenated DsrAB alignment block (alignment length = 1230 positions, n = 1218 homolog pairs). The tree is displayed with an artificial midpoint-rooting visualization. Branches are identified according to major taxonomic/functional groups (following the classification scheme of [49]). Bootstrap values >90 for clade-level bifurcations are indicated by black circles while those between 50 and 90 are denoted by white circles. Values are not depicted for nodes with bootstrap support <50. Branch length is relative to the scale provided at the top indicating the expected number of substitutions per site. The clade-level triangles indicate the phylogenetic diversity within each group via side lengths that are proportional to the distances between the clade’s closest and furthest taxa. b Paralogous rooting of the DsrA (top) and DsrB (bottom) subunits using a subset of taxa (n = 62) from the major lineages shown in the more expansive phylogenetic analysis in a. Bootstrap values are depicted the same as in panel a. Branch length is relative to the scale provided at the left indicating the expected number of substitutions per site. The root is labeled where the DsrA/DsrB subunits are presumed to have originated by a gene duplication prior to expansion and radiation. The clade-level triangles indicate the phylogenetic diversity within each group via side lengths that are proportional to the distances between the clade’s closest and furthest taxa.
DsrA and DsrB share significant homology and are likely to have arisen from gene duplication prior to the radiation of SRO [48, 49]. Thus, to better discern the phylogenetic placement among early-branching DsrAB, paralogous rooting was conducted for DsrA and DsrB (Fig. 3b). These analyses indicated that the YNP-TEG archaeal DsrAB lineage is among the earliest evolving DsrA and DsrB groups, which is consistent with the basal branching of crenarchaeal DsrA and DsrB in other analyses [49, 50]. Identical basal-branching topologies were recovered in phylogenetic reconstructions using multiple Maximum-Likelihood (ML) algorithms (RAxML and IQ-TREE) and/or with the use of site-heterogeneous amino acid frequency mixture models that have been shown to mitigate long branch attraction of rapidly evolving lineages [56]. Specifically, the placement of the branches shown in Fig. 3b were invariant to ML algorithm or substitution scheme. Given the deep-branching of the YNP-TEG DsrAB, monophyly with Crenarchaeota/Thaumarchaeota DsrAB, and the general concordance of Thermoproteales DsrAB with taxonomic phylogenies suggesting vertical transfer (Supplementary Fig. S9), one of several evolutionary scenarios could explain the presence of SRO in the Diaforarchaea. This includes (1) DsrAB were horizontally transferred into an ancestor of the YNP-TEG group from an organism harboring primitive DsrAB homologs, followed by a subsequent transfer to the ancestor of the Thermoproteales order comprising Vulcanisaeta, Caldivirga, Thermoproteus, and Pyrobaculum (and a later transfer to Thaumarchaeota; [32]), (2) DsrAB were present in the most recent common ancestor of the YNP-TEG and the Thermoproteales and were subsequently lost among other euryarchaeal and crenarchaeal lineages, or (3) DsrAB were acquired relatively recently via a horizontal transfer into the YNP-TEG lineage, and separately into the Thermoproteales. The second evolutionary scenario is highly implausible considering the breadth of phylogenetic distance separating the Diaforarchaea within the Euryarchaeota and the Thermoproteales within the Crenarchaeaota (Supplementary Fig. S9). In addition, if the third scenario was likely, the DsrAB homologs of the YNP-TEG and Thermoproteales lineages would be phylogenetically nested within other ancestral donor lineages, as is well-documented for DsrAB in the euryarchaeal order Archaeoglobales [49].
Consequently, the combined phylogenetic analyses suggest that the YNP-TEG group diverged from the S0/Fe(III) reducing DHVE2 group after the divergence of the combined YNP-TEG/DHVE2/Thermoplasmatales group from other Diaforarchaea. This taxonomic divergence was likely associated with either the maintenance of ancestral SO32− reduction PCGs in the YNP-TEG lineage or acquisition of these PCGs from an unidentified or extinct lineage prior to the early divergence of all DsrAB into the three contemporary higher-order groups described above. Nevertheless, identification of additional deep-branching archaeal DsrAB, if these indeed exist, would improve the resolution of these potential evolutionary scenarios.
Core metabolism of the YNP-TEG group
The YNP-TEG MAGs encode a number of additional proteins commonly involved in the dissimilatory SO32− reduction pathway (Supplementary Dataset 1) including DsrC (with conserved C-terminal cysteine residues; Supplementary Fig. S7) that acts as a mediator of SO32− reduction via the formation of a trisulfide bond [57]. It has been proposed that the membrane bound DsrMK(JOP) complex couples the reduction of the trisulfide bond to menaquinol oxidation [57]. This reaction would result in the downstream transport of protons (e.g., via Nuo) across the plasma membrane that could then be used in ATP synthesis (Fig. 2).
Neither of the YNP-TEG assemblies encoded identifiable homologs of sulfate adenylyl transferase (Sat) or adenosine-5′-phosposulfate reductase AB subunits (AprAB) that catalyze the activation of SO42− to adenosine-5′-phosphosulfate (APS) and APS reduction to SO32−, respectively. HMM-based searches of Sat proteins did not yield any positive matches, while those with AprA models yielded only potential orthologs that upon further inspection did not exhibit conserved biochemical features of characterized AprA of SRO (discussed further below). Moreover, the absence of these proteins was confirmed using tBLASTn searches of homologs from the most closely-related taxa with DsrAB homologs, Vulcanisaeta distributa, Caldivirga maquilingensis, and Pyrobaculum islandicum. While AprA-like proteins were detected in both YNP-TEG assemblies, as indicated above, phylogenetic analyses and the lack of conserved AprA-typical residues indicated that they did not share homology with AprA from SRO that genomically co-localize with SO42− reduction genes (qmoABC and/or sat), but rather belonged to a divergent family of paralogous proteins (Supplementary Fig. S10, Supplementary Text). The lack of identifiable Sat and AprAB homologs in MV2-Eury and SJ3-Eury MAGs, despite being nearly (98%) complete and from springs separated by >100 km and from different years suggests that these organisms lack the ability to activate SO42− and reduce APS. Moreover, a second MV2-Eury MAG produced from sediments collected from MV2 in November 2018 (estimated completeness: 98%; ~100% ANI to the 2014 MV2-Eury MAG) confirmed the lack of Sat and Apr homologs. This is unlikely to be an artifact of genome binning, since a survey of all unbinned contigs and other MAGs from the 2014 MV2 and SJ3 metagenome did not yield homologs that could be attributed to the two organisms (Supplementary Dataset 3) combined with three nearly identical, high-quality, effectively complete MAGs produced from different springs and from multiple years. Lastly, a screen of the MV2-Eury and SJ3-Eury MAGs indicated the absence of other previously published marker gene homologs for S0/sulfide oxidation (Sqr; Sox; Fcc; DsrEFH) and SO42− reduction (DsrD) (Supplementary Text).
Consistent with the conclusion that YNP-TEG are unlikely to be capable of SO42− reduction, enrichment cultures from MV2 spring sediments sampled in May 2018 that were supplied with HSO3− (most abundant form at pH 4–5) as the oxidant and peptone as a source of carbon and reductant (in addition to H2) yielded considerable sulfide (final concentration of 629 µM) after incubation for several weeks at 60 °C. However, despite numerous attempts, successful transfer of initial enrichments have not been achieved, disallowing more complete characterization. In contrast, enrichment cultures provided with SO42− did not produce sulfide over this time interval, including those amended with H2/CO2, lactate, peptone, or H2/CO2 plus peptone or H2/CO2 plus lactate. Importantly, samples are acidified in the methylene blue assay used to detect sulfide, which allows for detection of metal sulfides that might have been generated by SRO in sediments included as inoculum. Thus, binding of sulfide in mineral matrices would not likely affect the sulfide production values observed here. These results are thus consistent with predictions from genomic data and suggest the capacity for HSO3− (but not SO42−) reduction in MV2 populations. We attribute the activity observed in the enrichments from MV2 to MV2-Eury populations, since the only other putative SRO identified in the MV2 sediment community was a Thermoproteales-associated organism with potential capacity for SO42− reduction via Sat and AprAB (Supplementary Table S3, Supplementary Dataset 3). Nevertheless, because the inoculation sample was not taken at the same time as the metagenomic analysis sample, it is possible that other low-abundance SO42− reducers could have been present when sampling for cultivation. However, given that the same sediment samples were used to inoculate both the SO42− and SO32− cultures, it is unlikely that an SRO would be active in SO32− reduction, but not SO42− reduction, unless it was an obligate SO32− reducer, as is inferred for the MV2-Eury populations. Subsequent ongoing mixed-population enrichments from samples collected in Fall 2019 have also indicated the presence of SO32− reduction activity after months of incubation. Moreover, the presence of the MV2-Eury populations in these cultures was confirmed based on PCR-based detection of MV2-Eury-specific dsrA using population-specific primers described in the materials and methods.
Possible reductants capable of coupling with SO32− reduction in the YNP-TEG organisms are hydrogen (H2) or organic acids (e.g., formate, acetate, and lactate; Fig. 2, Supplementary Text). In addition, the MAGs encoded numerous peptidases and aminotransferases, in addition to dipeptide and single amino acid transporters. This indicates a general ability to metabolize cellular and extracellular proteins and amino acids via heterotrophic metabolism (Fig. 2) and is potentially consistent with the aforementioned enrichment data. Further, a hypothesized CO2 fixation pathway was present in the YNP-TEG MAGs (Supplementary Text), but it is unclear if the pathway is active and whether these organisms can grow autotrophically or mixotrophically.
Ecological distribution of the YNP-TEG
MV2-Eury is a consistent member of the sediment community of MV2, as evidenced by the detection of 16S rRNA genes with homology to this MAG in samples taken between 2010 and 2018 [25, 26] despite fluctuating geochemistry over these time intervals (Supplementary Table S1). Intriguingly, MV2-Eury populations were present at multiple sampling events between 2016 and 2018 over several seasons as shown by detection of DsrA homologs in DNA extracts (data not shown). However, active transcription of DsrA was only detected in sediments when MV2 exhibited a more moderately acidic pH (>3.5) in June–July 2018 than at any other sampling period over the previous two years (e.g., when spring pH < 3.5; Fig. 4). Given the high abundances of MV2-Eury populations in samples taken when the pH of MV2 was >3.5 [25, 26], in addition to the presence of a nearly identical MAG in SJ3 spring that exhibits a higher pH of 5.4, it is likely that the YNP-TEG generally occupy a higher pH niche that is variably present in the geochemically dynamic MV2 spring.
The mean values of triplicate qPCR assays are shown and error bars indicate standard errors. Sediment samples are arranged by MV2 spring pH with corresponding spring temperatures and sampling dates below each. n.d. indicates that no DsrA expression was detected in the sediment community RNA.
A survey of MV2-Eury-like 16S rRNA gene phylotypes (>97% nucleotide identity) among previously published archaeal community 16S rRNA gene datasets from 48 YNP springs [25, 26] revealed a patchy distribution in thermal springs with moderate temperatures and moderately acidic pH (Fig. 5a). Specifically, MV2-Eury phylotypes were only present in >1.0% relative abundance in five of 48 YNP springs with available data (two of which included MV2 from two different years). The pH of these spring waters ranged from 1.8 to 5.4 and the temperatures ranged from 42.5 to 62.0 °C. Likewise, a survey of MV2-Eury homologs (DsrA and RpoB) among 33 chemosynthetic metagenomes we generated from YNP springs as part of this and other projects ([28, 32, 58]; Supplementary Table S4), revealed a similar distribution (pH: 3.1–5.4; Temperature: 62–76 °C), albeit with lower estimated population abundances than the 16S rRNA gene BLAST searches (Fig. 5b). Relative abundances estimated from 16S rRNA gene abundances are known to be problematic for several reasons including overestimation due to multiple 16S rRNA gene copies within genomes [59], PCR primer bias [60], and the inability of certain PCR primers to amplify certain taxa, an issue that is particularly exacerbated for Archaea, and specifically Crenarchaea [61], that are often abundant in thermal springs. Thus, the estimated abundances of the MV2-Eury-like phylotypes in the metagenomic datasets are likely to be more accurate estimates of in situ abundances. Nevertheless, the pH and temperature ranges where YNP-TEG-like populations were identified were broadly consistent in both datasets, despite being derived from different thermal springs. As discussed below, the stability and chemical speciation of SO32− shifts considerably in the pH range of ~3–6, which may ultimately bound the niche space occupied by SO32− reducing YNP-TEG within hydrothermal systems.
a 16S rRNA gene distribution revealed by comparison against 32 YNP hot spring communities described in Colman [26] and 15 communities evaluated in Colman et al. [25]. Positive detection was defined as when 16S rRNA gene identities >97% were identified relative to the MV2-Eury 16S rRNA gene. Points are colored by estimated relative abundance, as given by the scale on the right. Gray circles indicate the lack of detection in >0.01% relative abundance. b Inferred distribution of YNP-TEG-like phylotypes using an in-house database of 33 chemosynthetic hot spring community metagenomes that we have previously generated from YNP springs (Supplementary Table S4). Relative abundances were estimated based on DsrA/RpoB amino acid identities >95% to the MV2-Eury phylotype in MAG-binned contigs, followed by relative abundance estimates calculated by read-mapping as described in the materials and methods. The distribution of MV2-Eury-like DsrA and RpoB homologs were identical, and thus only one scatter plot is shown. Of the four metagenomes with positive YNP-TEG detection, two were MV2 from different years and another was SJ3 (Supplementary Table S4). Note that the samples in a and b derive from different datasets and do not necessarily represent the same springs. Likewise, the temperature and relative abundance scales differ between the two panels to emphasize within-dataset variation.
(Bi)sulfite prevalence in hydrothermal environments
SO32− is extremely unstable in modern oxic environments and rapidly oxidizes abiotically to SO42− [16, 62]. However, SO32− (primarily HSO3− in aqueous solutions with pH 2.0–5.5; Supplementary Fig. S11; [63]) can be measured in appreciable concentrations in some hydrothermal springs [64]. At temperatures <150 °C [15], volcanic SO2 can hydrate and ionize to form SO32−. Further, thiosulfate (S2O3−) produced during oxidation of sulfide [62] is unstable at pH < 4.0 and degrades to S0 and SO32− [62]. Thus, several abiotic mechanisms exist to supply SO32− in volcanically-influenced hydrothermal environments. Indeed, SJ3 and other springs within the Smokejumper geyser basin and MV2 within the Mud Volcano geyser basin exhibit some of the highest volcanically-sourced gas concentrations in YNP [29]. Thus, volcanic gas is especially abundant in these two areas of YNP, which may deliver a consistent supply of SO32− to these springs via SO2 or H2S processing.
The primary mechanism leading to SO32− oxidation to SO42− in contemporary hydrothermal environments is via O2− or Fe3+ ions as oxidants [16], although atmospheric photolysis of SO2 could potentially contribute SO32− to near surface waters [2]. Thus, in the absence of either O2 or Fe3+ ions, SO32− should be relatively stable, in particular in subsurface environments or within spring sediments where MV2-Eury is predominantly found [25]. Indeed, laboratory experiments conducted here indicate that SO32− solutions that were prepared anoxically maintained >80% of the original SO32− after 24 h incubations (Fig. 6a), while those incubated in the presence of O2 had very little SO32− remaining after incubation over the same period (Fig. 6b). SO32− stability was more pronounced in solutions with moderately acidic pH, with ~70–80% of SO32− remaining after anoxic incubations for 24 h at pH between 4 and 6. In addition, increased temperatures (i.e., ≥60 °C) resulted in increased SO32− oxidation under oxic conditions, but not under anoxic conditions. Thus, hot springs that harbor moderately acidic pH, are of moderate temperature, are anoxic or suboxic, and that feature high inputs of volcanic gas favor SO32− stability that, in turn, increases its bioavailability for organisms such as YNP-TEG. Importantly, these are distinguishing characteristics of the sulfur-rich MV2 and SJ3 springs, but such springs are relatively uncommon in the YNP geothermal system [32], which may explain the relatively limited ecological distribution of YNP-TEG organisms among YNP springs (Fig. 5a, b).
a (Bi)sulfite stability under anoxic conditions and at pH values that encompass the range observed in MV2 and SJ3 springs. The mean percentage of SO32− remaining 24 h after adding Na2SO3 to a concentration of ~50 µM is shown for triplicate microcosms prepared at each pH interval and incubated at 40 °C, while error bars show the standard deviation of the assays. b The effect of temperature and oxygen (20% O2, 80% N2) exposure on the stability of (bi)sulfite at pH 4. Microcosms were prepared and measured as in Fig. 6a, but at varying temperatures and with or without preparation under anoxic conditions.
Implications for the origins of SRO
The identification of deeply-branching DsrAB homologs in organisms that characteristically inhabit moderately acidic hot springs prompts the question as to whether such environments might have existed at the time that SRO are thought to have originated. SRO are thought to have been present as early as 3.47 Ga, based on sulfur isotopic data from sulfides in barite deposits from the Dresser Formation [2]. Recent data suggest that these deposits formed in an ancient volcanic caldera and several of these, in particular the barites, were likely formed by circulating hydrothermal fluids in a subaerial hot spring setting [4]. Moreover, the detection of minerals, specifically kaolinite and illite minerals, suggest that acid-sulfate conditions were likely present within this caldera-like environment, at the very minimum, during the final stages as this high-sulfidation volcanic system evolved [4]. Like moderately acidic springs in YNP, such as SJ3 and MV2, the acidity and sulfate in springs inferred to have existed in this ancient volcanic caldera in the Dresser Formation were likely generated by condensation of volcanic gases containing SO2 and H2S with groundwater, which can produce acidic, sulfur-rich conditions [62]. When combined with evidence suggesting that volcanic activity was more prevalent during the Archean eon (4.0–2.5 Ga) relative to today [11], these data lead to the conclusion that moderately acidic sulfur-rich springs capable of supporting ancestors or analogs of YNP-TEG likely existed on early Earth and were potentially prevalent.
The phylogenetic evidence described above, which is broadly consistent with previous analyses [49, 50], suggests that cultivated thermophiles or inferred thermophiles (Moorella spp.; the YNP-TEG; Crenarchaeota genera ubiquitous among hot springs including Vulcanisaeta, Pyrobaculum, Thermoproteus, and Caldivirga; the c.d. “Hydrothermarchaeota”; and c.d. “Aigarchaeota”) from hydrothermal settings harbor the earliest diverging DsrAB homologs currently known, and are thus reflective of the earliest SRO (Fig. 3b). While two other enzyme complexes, Asr and Mcc, can also reduce SO32− in dissimilatory metabolisms [65, 66], both exhibit patchy taxonomic distributions and are generally found in derived lineages [50], observations that together suggest them to be recently evolved. The only significant exceptions to inferred thermophily among the early-branching DsrAB homologs are from MAGs representing the c.d. “Rokubacteria” from subsurface environments and a Verrucomicrobia MAG recovered from peat soils (Fig. 3b), both of which are suggested to have acquired their putative sulfur oxidizing or SO32− reducing/sulfur oxidizing functions, respectively, from one or multiple horizontal transfers of DsrAB and other companion proteins [50]. Among the inferred metabolic coupling of early-diverging DsrAB: Moorella spp. do not reduce either SO42− or SO32− (but generally reduce S2O3−) rendering the function of both of their DsrAB copies unclear. In contrast, the YNP-TEG are only potentially capable of reducing SO32− while the Crenarchaeota SRO comprise taxa reported to reduce SO32− and SO42−. The Verrucomicrobia MAG SbV1 [67] lacks Sat and AprA (confirmed with HMM searches) rendering SO42− reduction unlikely while only one of the two “Hydrothermarchaeota” MAGs with DsrAB (JdFR-18) contains Sat and AprAB-like proteins (JdFR-17 does not). However, the Sat and AprAB-like proteins in JdFR-18 are not genomically co-localized with either of two copies of DsrAB but are co-localized with a homolog of F420-linked sulfite reductases of methanogens, and the AprA-like protein lacks the characteristic residues strictly conserved among AprA homologs (Supplementary Fig. S10) This rendors the function of both DsrAB copies to be unknown and suggests they are not likely be coupled to SO42− reduction. Likewise, the c.d. “Aigarchaeota” MAG that contains DsrAB does not contain AprAB or Sat. Finally, the c.d. “Rokubacteria” MAGs likely acquired sulfur oxidation potential through horizontal transfer of various Dsr proteins from various sources [50], including DsrAB potentially transferred from organisms harboring deeply-rooted DsrAB homologs like the Crenarchaea [50]. Thus, the potential for the earliest diverging DsrAB homologs to be coupled to SO42− reduction appears to be highly limited. A more parsimonious explanation of the above data is that the early-diverging DsrAB homologs are generally linked to SO32− reduction rather than SO42− reduction, as exemplified by the YNP-TEG MAGs.
Although SO42− was likely to be generally less abundant in Archean environments when SRO are thought to have originated when compared with contemporary conditions [8, 9], the possibility that SO42− utilizing SRO emerged in a thermal environment that hosted elevated SO42− is an equally plausible scenario to explain the isotopically depleted sulfides in the Dresser Formation barites. At temperatures greater than ~150 °C, volcanic SO2 can disproportionate to form S0/H2S and SO42− [15] allowing for localized enrichment of SO42− at concentrations capable of supporting SRO. Alternatively, SO42− derived from photochemical oxidation of atmospheric SO2 and subsequent rain-out could have also contributed SO42− to nonthermal surface waters [12]. Further, localized evaporation of these waters could generate still higher SO42− concentrations [14]. Intriguingly, it was recently shown using purified DsrAB that the primary enzymatic influence on sulfur isotope fractionation during SO42− reduction is at the level of reduction of SO32− to sulfide (via DsrAB) [68], although a fractionation effect by enzymatic SO42– activation and reduction (i.e., by APS reductase; AprAB) may also occur [69]. Thus, it is unclear whether SO32− or SO42− reducing organisms (or both) were responsible for the measured sulfur isotope fractionations between sulfides produced from such activities and barites in the ~3.5 Ga Dresser Formation hydrothermal deposits. An alternative explanation, and one that is consistent with the discovery of deeply diverging SO32− reducing (but not SO42− reducing) thermophilic populations described here and elsewhere, is that SO32− reducers were present in the Archean in volcanically-influenced environments either as predecessors or as contemporaries of SO42− reducers. If true, then the ultimate global dominance of SO42− reducers would have occurred with increasing input of SO42− into oceans via oxidative weathering of continental sulfide minerals. Such an explanation is consistent with previous suppositions regarding the primitive nature of SO32− reducers and the later dominance by SO42− reducers based on several lines of inference [21,22,23]. The buildup of atmospheric O2 would have allowed for the global dominance of SO42− reducers and also would have had the opposite effect on SO32− and the organisms that solely use SO32− as an oxidant since it is unstable in the presence of O2. One notable exception is anoxic or suboxic, moderately acidic and moderate thermal environments that are infused by SO2, such as volcanically-influenced hydrothermal environments that likely have persisted in a relatively similar state since early in Earth history.
In support of this hypothesis, SO32− reduction represents a simpler metabolic pathway, involves fewer enzymes, and, following the ideas of the evolution of metabolic pathways proposed by Granick [70], may therefore have originated prior to the more complex SO42− reduction pathway. In addition to the above evidence for the evolutionary primacy of SO32− reducers, the majority of the free energy yield from SO42− reduction arises from SO32− reduction to sulfide, while the first SO42− activation step with ATP is endergonic, and the free energy yield from APS reduction to SO32− is minimal [24]. The retention of the ability of many SO42− reducers to grow on SO32− combined with the considerably higher growth yields and efficiencies achieved by widely-studied SO42− reducers including Desulfotomaculum orientis [71] and Desulfovibrio vulgaris [72] when grown on SO32− relative to SO42−, may provide additional evidence for the primitive nature of the SO32− reduction pathway. Consequently, the YNP-TEG archaeal group may represent descendants or analogs of primitive SO32− reducing organisms and important models to understand the early evolution and emergence of SO32− and/or SO42− reduction pathways. While preliminary data indicate that MV2 sediment populations can reduce SO32− (but not SO42−), further efforts are needed to identify suitable conditions to maintain and propagate these cells for use in physiological, biochemical, and isotopic characterization.
Source: Ecology - nature.com
