Introduction

Yellowstone National Park (YNP), Wyoming, USA, hosts the world’s largest active continental hydrothermal system and its surface expression manifests as the >10,000 fumaroles, geysers, mud pots, and hot springs that are distributed across its landscape1,2,3. Numerous studies of hydrothermal features in YNP have been performed since the late 1800s, making it arguably the most well-studied hydrothermal system in the world. These studies have collectively documented extensive spatial variation in the geochemistry and microbial diversity of surface hydrothermal features. However, far less is known of the subsurface hydrothermal system in YNP, including whether it is habitable for microbial life.

In the 1960s, the United States Geological Survey drilled a series of 13 experimental wells in thermal areas of YNP to study groundwater chemistry, heat flow, and mineralogy4. Shortly after drilling, a sample of water from one of the wells (not identified) was filtered to concentrate microbial cells for microscopic examination5. No cells were detected, which was perhaps not unexpected given the temperature of the water was estimated to be ~200 °C. Three decades later, another of these wells, well Y-7 in Biscuit Basin that intersects fracture fluids sourced by the parent hydrothermal aquifer water, was examined for the presence of microbial cells6. In this study, glass slides were suspended from the surface (47 °C) to near the bottom of the well at 27 m (93 °C). Following incubation for several weeks to several months, the slides were fixed for microscopic analyses. Few if any cells were detected on slides incubated at depth while those incubated near the surface exhibited more colonization. Colonization differences were attributed to the lack of available oxidants capable of supporting microbial life. Indeed, chemical analyses of the waters indicated depletion in oxidized sulfur species, including sulfate that would be expected to be present at ~30 mg L−1 in circumneutral to alkaline waters sourced by the parent hydrothermal aquifer2,7. Importantly, the well was cased to a depth of 32 m, which combined with mineral precipitation and self-sealing that took place over the intervening 30 years between when it was initially drilled and then later sampled, effectively limited the exchange of well waters with the local aquifer and led to unnatural hydrologic and geochemical conditions8. Thus, it is unclear how well these samples represent actual subsurface conditions where water would be expected to move more freely along natural fractures and that would allow for mixing with other fluid types. To the authors’ knowledge, only one other study has evaluated the presence of microbial communities in high-temperature hydrothermal field waters that are potentially reflective of subsurface conditions via Icelandic borehole fluids collected at surface wellheads9. Although microbial community analyses of these waters via amplified 16S rRNA gene inference revealed the presence of diverse putative thermophiles, the collection of these fluids at the wellheads renders it unclear how well these communities reflect native subsurface communities. Consequently, the presence and characteristics of subsurface microbial ecosystems in continental geothermal fields in YNP, or elsewhere, remain inconclusive.

Deep hot spring pools offer an alternative mechanism to examine subsurface geochemistry and microbiology under conditions that allow the natural movement of subsurface fluids. Cinder Pool (CP) is a near-boiling spring in the Norris Geyser Basin (NGB) of YNP and is one of the deepest springs in the Yellowstone hydrothermal system1. CP is an acid–sulfate–chloride spring (pH ~ 4.0) with a geochemical composition consistent with it being sourced from the deep parent hydrothermal aquifer7 in addition to minor input of near-surface meteoric water10,11,12,13. Mapping of the physical structure of CP in August 1969 revealed a funnel-like architecture at the surface (89.9 °C) with a 9 m diameter surface circumference that narrows to roughly 1 m at a depth of 11.6 m. This architecture extends to 18.3 m were a “mushy” substance with a temperature of 94.0 °C was recorded10,11. After allowing a thermistor and lead weight to settle to a depth of 19.2 m, the temperature rose to 96.0 °C and by 20.4 m, the temperature increased to 117.5 °C. Retrieval of the lead weight revealed that the mushy substance was molten elemental sulfur (S°) that exhibited a metallic gray color due to dispersed iron sulfide10,11.

The molten S° at depth in CP is unique and has been hypothesized to result from the spring intersecting an ancient solfatara that accumulated substantial S° 11. This S° ‘lid’ allows for the unusual thermal gradient in CP10. The hollow “cinders” found floating on the surface of CP form from rising gas bubbles that buoyantly carry molten sulfur into the cooler overlying water column, allowing them to rapidly cool and trap gas10. The molten S0 pool and the sulfur spherules in CP also control pathways of sulfur cycling that maintain spring pH to ~411, while also giving rise to the highly atypical geochemistry and sulfur cycling dynamics of CP11. In particular, the formation of thiosulfate (S2O32−) via S0 hydrolysis and the oxygen-dependent decomposition of S2O32− to polythionate (SxO62−) catalyzed by trace iron sulfide in cinders are reactions that are uniquely observed in CP relative to other characterized YNP springs11. SxO62− decomposition can occur via reaction with sulfide yielding S2O32− and S0 or by disproportionation to SO42− and S2O32−. The combination of these reactions limits the availabilities of both O2 and sulfide that are needed to drive the formation of sulfuric acid (pKa = 2.5–3.0 between 70 and 100 °C14), the primary buffer of highly acidic springs in YNP15,16,17.

The temperature of CP at the surface has historically ranged from ~77–93 °C10,11,12,18,19,20, and fluctuates seasonally, likely due to the near-surface input of cooler meteoric water13. Regardless of the surface temperature at CP, several studies have shown that temperature increases with depth, with an abrupt increase to the melting point of S0 (~120 °C) at ~18 to 20 m10,12. The unique geochemical conditions of CP and the presence of a molten pool of S° in the subsurface are a consequence of its unique depth, size, and history (intersecting an ancient solfatara), rendering it highly unique among global hot springs and leading to its description as “one of the most unique springs”11 and “one of the most extraordinary features” in YNP10. The unique conditions of CP have apparently persisted since it was first scientifically described in 19271, with numerous subsequent studies indicating relatively little change in its overall characteristics. However, between autumn 2018 and spring 2019, the pH of CP waters decreased to pH ~ 2.5, its sediments became visibly coated in S0, and the floating cinders disappeared. It is not yet clear what caused this pronounced change, although it is possibly related to increased surface deformation, subsurface gas accumulation, and seismic activity in NGB at this time21 and that have together been suggested as responsible for the reawakening and increased eruptive activity of the world’s tallest geyser, Steamboat Geyser21,22. As such, CP provides a unique opportunity to not only examine the potential presence and temporal dynamics of microbial life in the subsurface of YNP, but also to assess the potential mechanisms by which subsurface microbial communities could influence the temporal dynamics of hot spring geochemistry. As mentioned above, while numerous studies have examined spatial distributions of microbial taxonomy and functionality in hot springs as a function of geochemistry, far fewer have examined temporal dynamics in hot springs, and none have evaluated subsurface microbial community dynamics in continental volcanic systems.

Temporal variation in spring geochemistry is related to fluctuations in one or more of a thermal feature’s three main components: water, heat, and the flow path that delivers fluids to springs23. Short-term and minor fluctuations in these components can be caused by changes in the circulation of subsurface waters due to deposition of silica and subsequent sealing of the flow paths that source springs or to alteration of these flow paths due to chemical and/or physical weathering24. For example, a seasonal study of four springs in the NGB that included CP showed coupled shifts in the geochemical and surface water microbial compositions of springs that were associated with near-surface input of oxidized meteoric water. The timing of these coupled changes was largely synchronous with recent precipitation events, suggesting that hydrologic factors can influence spring geochemical and microbial compositions on relatively short time scales. Alternatively, large-scale changes in the three components of hot spring systems can instantaneously occur due to seismic activity24. For example, a study of Evening Primrose25, a moderately acidic spring in the Sylvan Springs area of YNP, showed that the geochemical composition of the spring has been influenced by several periods of increased seismic activity over the past 90 years. While the spring was circumneutral to alkaline between 1927 and 1959, a large (magnitude 7.3) earthquake altered the hydrologic paths sourcing the spring, allowing for increased input of volcanic gases, that when oxidized, contribute to acidity. By early 1960, the spring rapidly became sulfur-rich and highly acidic (pH < 2.0) and remained this way until ~2000–2004 when another period of elevated seismic activity is hypothesized to have again altered the flow path of fluids sourcing the spring leading to increased pH (~5.4) that has remained until at least 201625. Collectively, these studies are consistent with the notion that hot spring geochemistry can vary over a range of time scales3 and changes are driven by dynamic feedbacks between hydrologic, geologic, and biologic factors3,13,25.

Here, we examine the geochemistry and microbial community in CP (1) to examine the potential for a subsurface hydrothermal biosphere in YNP and (2) to investigate coupled subsurface geobiological changes associated with the acidification of this spring. Samples were collected from a vertical depth profile in CP in 2016 and again in 2020 for geochemical and microbial characterization. The results of these studies indicate that the pronounced, rapid, and coupled geobiological change in CP is likely attributed to the disappearance of the molten S0 at a depth that allowed for a shift in biogeochemical sulfur cycling pathways from reducing conditions to those that favored oxidation and formation of sulfuric acid. Results implicate an uncultured sulfur-oxidizing Sulfolobales archaeon in the acidification of the spring waters and suggest that this organism could be endemic to the shallow subsurface of YNP where the majority of sulfur oxidation and water acidification is likely taking place10,26. These findings underscore the potential for dynamic and rapid feedbacks between the geosphere and biosphere in the subsurface of active, volcanic hydrothermal environments.

Results and discussion

Acidification of CP

Historical geochemical data suggest that the water chemistry of Cinder Pool (CP) has been relatively stable from the time of first reported geochemical data in 1947 until autumn 2018, followed by pronounced acidification between winter and spring 2019 (Supplementary Data 1, Fig. 1a, b). Images and documentation dating to even earlier (1927) reveal the presence of cinders covering ~50% of the spring surface at that time, a temperature near boiling (91.5 °C), and a description of having high sulfate and chloride levels (although data was not provided), suggesting that its chemistry has been generally stable since its discovery1. Spring pH ranged between ~3.6 and 4.5 in 22 yearly measurements spanning 71 years (1947–2018; multiple measurements in the same year were averaged to represent each year) (Fig. 1b), while the pH has been subsequently measured after 2018 as low as 2.5 (Fig. 1b). A single pH measurement of 2.5 was also recorded in a 2003 publication27, although other measurements in 2003, 2000, and 2001 were more consistent with the long-term average (i.e., pH 4.2–4.3; Supplementary Data 1). Scrutiny of chemical data accompanying the pH 2.5 measurement in 2003 indicates a SO42− concentration (~48 mg L−1) that is considerably lower than would be expected for CP, even when the pH is much higher (SO42− = 80 mg L−1; pH = 4.2–4.3). Considering that sulfuric acid is the predominant buffer of pH in these systems7,28, the pH 2.5 reading in 2003 is considered questionable. Nevertheless, the 2018 shift in pH towards more acidic conditions was accompanied by a notable change in the appearance of CP. Prior to autumn 2018, the spring waters were cloudy gray with the considerable suspension of kaolinite clay particles20 and black cinders10. However, between autumn 2018 and spring 2019, the spring waters visibly turned blue-green and contained colloidal S° particles that were also deposited along the pool shelves, while the pool also lacked its characteristic black cinders (Fig. 1a). The spring has maintained this appearance since spring 2019 until at least July 2022.

Fig. 1: Historical geochemistry of Cinder Pool (CP).
figure 1

a Top panel shows the visual change in the appearance of CP in 2016 (left) and 2020 (right). Scale bars in the bottom right are ∼1 m. b Measurements of pH (n = 21; black line) and sulfate (SO42−) concentrations (n = 12; red line) in CP waters between 1947 and 2021. Years with multiple measurements were averaged to represent the entire year. c Paired measurements of SO42− and chloride (Cl−) concentrations (n = 12) between 1947 and 2021 in the context of the same measurements for 488 YNP springs derived from previous studies. Paired points for CP are colored based on the year they were recorded (averaged for multiple measurements/year as described above). End member fluid compositions as described in the manuscript text are indicated based on the abbreviations: MO meteoric only, HO hydrothermal only, MG meteoric plus gas, HB hydrothermal plus boiling, HBG hydrothermal plus boiling plus gas. Points for 2016, 2018, 2019, 2020, and 2021 are indicated by “16”, “18”, “19”, “20”, and “21”, respectively.

Full size image

The source of fluids in YNP hot springs can be broadly defined by concentrations of sulfate (SO42−) and chloride (Cl−)2,7. These indicators have been previously used to define the source of YNP springs as either (1) hydrothermal only (HO) waters that have moderate concentrations of SO42− (~30 mg L−1 depending on the depth of boiling; described below) but high concentrations of Cl− (~300 mg L−1), (2) meteoric-only (MO) waters containing lower concentrations of both solutes, or (3) MO waters infused with gas (MG) that have lower Cl− concentrations and higher SO42− concentrations (Fig. 1c). Subsequent boiling and/or evaporation of HO waters can concentrate Cl− and SO42− to higher concentrations (termed hydrothermal plus boiling; HB), while additional gas input into HO or HB waters can lead to particularly high concentrations of both Cl− and SO42− (hydrothermal + boiling + gas; HBG)7 (Fig. 1c). Geochemical data from surveys spanning 1947 to 2018 suggest that CP was largely sourced by hydrothermal (HO) waters that have undergone boiling and/or evaporation (HB) during this time frame (Fig. 1c).

HO and HB waters are typically circumneutral7, while CP (which is also sourced by HB waters) has maintained a moderately acidic pH of ~4 until autumn 2018 (Fig. 1b). Several other low pH HB waters have been previously observed within the NGB7. The moderately acidic pH in CP (prior to 2018) has been attributed to the hydrolysis of molten S° that occurs at depths of >18 m that leads to the formation of S2O32– 11. Oxygen (O2)-dependent oxidation of S2O32−, catalyzed by trace iron sulfide in the cinders, forms SxO62− that can then react with sulfide to yield S2O32− and S° 11. Alternatively, SxO62− can be disproportionated to form S2O32− and SO42− 11. The relative rates of these reactions in CP prior to 2018 are not known although similar concentrations of S2O32− measured between 1995 and 1997 suggest that rates of S° hydrolysis and rates of S2O32− formation have been relatively constant over yearly time scales11. The consumption of O2 by reaction with S2O32− and the consumption of sulfide involving reactions with SxO62− would limit the amount of sulfuric acid that could be formed, thereby maintaining a less acidic pH than other sulfuric acid buffered acidic springs in YNP7.

Between November 2018 and March 2019, the pH of CP markedly decreased to 2.8 in 2019, 2.7 in 2020, and 2.6 in 2021. This coincided with a marked increase in SO42− concentrations of ~3–5 fold above historical ranges (Fig. 1b), while Cl− concentrations fluctuated without clear trends during this time (Supplementary Fig. 1c). Thus, CP transitioned from an HB water type to an HBG water type between autumn 2018 and spring 2019 and has remained this way since (Fig. 1c). This is interpreted to reflect a substantial increase in H2S/S° oxidation that results in the formation of SO42− and H+ (sulfuric acid). Several observations suggest a fundamental restructuring of CP’s unique sulfur cycling due to dramatic physical and chemical changes at this time. As described in more detail below, the molten S° layer was detected at a depth of 18 m in 2016. However, in 2020 and 2021 there was no evidence of molten S° at ~18 to 20 m depth as previously documented, and sampling equipment could be freely dropped to a depth of 22 m (length of the cable) without interruption. In the absence of the molten S° at depth, the S° hydrolysis product S2O32−, and the cinders that catalyze SxO62− formation from S2O32− and H2S, it is possible that such reactions that previously competed for H2S or O2 (i.e., those involving S2O32− and SxO62−) are no longer taking place in CP. This in turn would allow for sulfur compounds (H2S and S°) to now be oxidized, thereby contributing to spring acidification.

Alternative scenarios underlying the dramatic changes in CP waters also warrant consideration, and the three most logical are presented below. First, it is possible that the waters sourcing CP may have shifted either via replacement of the primary source or by altered mixing of multiple water sources. Water isotope values (δ2H and δ18O) can be used to further deconvolute the sources of hydrothermal waters because distinctive isotope values are associated with distinct water sources and the various influences upon them including meteoric water recharge, boiling (and/or evaporation), and water–rock interactions7,29. The water isotope values measured among the measured depths in CP in 2020 were near the range of water isotope values observed in CP across multiple months in 201613 (depth-resolved water isotope measurements were not made in 2016). The 2020 CP water isotope values were slightly right-shifted relative to those of 2016, suggesting a minor increase in the evaporation and concentration of CP water isotopes between 2016 and 20207 (Supplementary Fig. 2). These data thus do not support the hypothesis that the source of waters in CP dramatically shifted between 2016 and 2020, consistent with the SO42− and Cl− measurements indicating that the primary change to CP waters was increased input or availability of H2S for oxidation.

A second alternative explanation is that a change in the water level of CP could potentially alter residence times which could allow for more oxidation of sulfur compounds in the spring and increased acidification. Such a scenario would also likely result in increased evaporation and concentration of solutes. However, the minimal increase in water isotope values (Supplementary Fig. 2) and similar Cl− concentrations (Supplementary Fig. 1c) accompanying a ~3–5 fold increase in SO42− concentration pre- and post-acidification (Fig. 1b) argue that increased residence time was of minimal importance in acidification.

A third possible explanation is that a change in the plumbing system of CP is now delivering more vapor phase gas that contributes H2S and acidity when oxidized. Such a scenario could be consistent with increased surface deformation, subsurface gas accumulation, and seismic activity that has been taking place near NGB just prior to these changes21, and the transition from HB-type to HBG-type waters in CP. Sulfur species isotope analyses would help deconvolute the sources of SO42− in CP, but samples for sulfur isotopic analyses were not collected prior to acidification. Thus, it is unclear if this process may also be contributing to the acidification of CP. Regardless, the disappearance of the molten S° cap either by consumption or displacement would in effect make H2S more available for oxidation, similar to increased vapor phase input. The acidification of hot springs involves the oxidation of H2S by O230. More specifically, partial oxidation of H2S at acidic pH (<4.0) leads to the formation of S2O32− which is unstable and disproportionates to form sulfite and S0 28. In the absence of microbial S0 oxidation or reduction activities and at temperatures <100 °C, S° is stable and often accumulates in acidic and sulfidic springs28. Aerobic S° oxidizing thermoacidophiles are thought to be the primary catalysts of S° oxidation and this involves reactions that contribute acid26,28. At temperatures >80 °C, such as in CP, S° oxidation is primarily catalyzed by members of the archaeal order, Sulfolobales7,28. Taken together, these observations imply that the near-surface acidification of CP waters was enabled by the disappearance of molten S° at depth and the cessation of reactions that S° hydrolysis made possible, including those that consumed H2S and O2. It is also possible that the increased input of vapor phase H2S into CP contributed to its acidification and increased overall sulfur content. In either case, H2S and O2 would have been made available for the oxidation of H2S or S° by thermoacidophiles, possibly related to members of the Sulfolobales7,10,15,16,31.

Subsurface geochemistry and microbial communities in CP

To better understand the existence of a subsurface biosphere in CP, depth-resolved sampling of the water column was conducted in July 2016. The in situ measured temperature of CP water at the surface in 2016 was 85.8 °C and gradually increased to 88.3 °C at 9 m and then abruptly to 93.1 and 95.6 °C at 12 and 15 m, respectively (Fig. 2b; Supplementary Data 2). At ~18 m, the temperature dramatically increased to ~120 °C. Upon retrieval, a metallic gray solidified material covered the sampling apparatus thermistor, weight, and support cable (Supplementary Fig. 3), consistent with previous reports of encountering a similarly viscous material at ~18 to 20 m with a temperature of ~120 °C and that had a metallic appearance upon solidification10. The marked increase in temperature over this short distance is attributed to the molten S0 acting as an insulating lid on the hydrothermal system10. Sampling was consequently restricted to the 0–15 m depth interval. pH did not vary considerably over the entire depth profile (Fig. 2a, Supplementary Data 2). Like temperature, ORP did not markedly vary over the 0–9 m intervals (−78 to −81 mV) but decreased to ~−117 mV at 12 and 15 m (Fig. 2c; Supplementary Data 2). A minimal correlation was observed between SO42− and Fe(II) concentrations with depth, but total sulfide concentrations notably increased with depth (Fig. 2d–f; Supplementary Data 2). These data suggested the presence of two thermally and chemically stratified bodies of water in the water column of CP: one that was more reduced, likely due to the influence of molten S° at depth, and the other that was more oxidized and influenced by contact with the atmosphere.

Fig. 2: Depth-resolved geochemical measurements for Cinder Pool (CP) waters collected in 2016 and 2020.
figure 2

Measurements of a temperature, b pH, c oxidation–reduction potential (ORP), d SO42−, e total sulfide (S2−), and f Fe(II) were taken from waters pumped from 0, 3, 6, 9, 12, and 15 m depth or using an in situ thermocouple (temperature only). Blue circles show measurements taken in 2016 and red circles show measurements taken in 2020.

Full size image

To continue to assess the possibility of a subsurface biosphere in CP, and to assess its potential role in spring acidification between autumn 2018 and spring 2019, depth-resolved sampling was conducted again in August 2020. The same sampling strategy was used to sample CP as in 2016. As observed in 2016, the temperature varied minimally between the surface (86.4 °C) to 9 m (87.3 °C) and slightly increased to ~90 to 91 °C at 12–15 m (Fig. 2b; Supplementary Data 2). As mentioned above, a molten S° layer was not encountered while lowering the sampling equipment to 22 m in 2020 and 2021 (Supplementary Fig. 4) and the temperature increased only to 93 °C, consistent with the absence of an insulating S° lid. Spring pH did not vary substantively across the depth profile in 2020. In contrast to the 2016 sampling profile, ORP did not markedly vary with depth, nor did SO42− and Fe(II) concentrations, although sulfide again increased with depth (Fig. 2c–f; Supplementary Data 2). Notably, measurements of ORP in the 2020 CP waters were substantially higher by ~400 mV (~+270 to +280 mv) than those of the 2016 CP waters (−78 to −117 mV), indicative of a much more oxidized and acidic spring. This is consistent with sulfuric acid buffering of this spring to a pH of 2.5–2.67.

To investigate the abundance, composition, and functional potentials of CP communities in 2016 and 2020, biomass was filtered from 0, 9, and 15 m depths and subjected to DNA extraction and shotgun metagenomic sequencing. Total DNA concentrations, normalized to L of spring water filtered, were much lower in the 2016 water samples than in the 2020 samples (Fig. 3a; Supplementary Data 2), suggesting lower overall planktonic biomass in 2016. The decreased biomass detected in samples from 2016 may be due to increased evidence for viral infections in those cells at that time (discussed below) that, if lytic, would be expected to decrease the number of cells captured on filters. Interestingly, while DNA concentrations considerably declined with depth in the 2016 waters, DNA concentrations increased with depth in the 2020 waters (Fig. 3a; Supplementary Data 2). This may point to the presence of populations in the 2020 samples that were better adapted to growing under subsurface conditions (e.g., lower O2 and increased H2S/S° concentrations relative to the surface). Measurements of the abundance of planktonic cells in waters collected from 2021 at 0, 9, and 21 m confirmed the increase of biomass with depth (Fig. 3b). Moreover, the presence of abundant cells at 21 m in samples from 2021 suggests that cells could be present at even greater depths within the CP subsurface.

Fig. 3: Characterization of Cinder Pool (CP) planktonic communities in depth-resolved samples taken in 2016 and 2020.
figure 3

a Total DNA extracted from biomass filtered from CP waters obtained from 0, 9, and 15 m, normalized to liters of water filtered. b Abundances of cells in waters collected from 0, 9, and 21 m depth of CP in 2021. Cellular abundances were not enumerated in 2016 water samples. Error bars show standard deviations from triplicate measurements. c Composition of CP water communities sampled from 0, 9, and 15 m. The relative abundance of taxonomic orders was evaluated based on the percentage of quality-filtered metagenomic reads mapped to each metagenome-assembled-genome (MAG) shown in the legend on the upper right. Viral and non-cellular reads are those that were assigned as a viral contig or were otherwise not clearly associated with cellular genomes. ‘Others’ comprise the remaining reads. Only the most abundant unbinned contigs that comprised 75% of all reads mapped to unbinned contigs were evaluated for viral/non-cellular assignment (Supplementary Data 3).

Full size image

Shotgun metagenomic sequencing, assembly, and binning were used to reconstruct population-level genome bins from samples collected in 2016 and 2020. Taxonomically, both 2016 and 2020 CP communities comprised only Archaea and were dominated by members of the class Thermoprotei (Table 1 and Fig. 3c), consistent with previous characterizations of acidic, high-temperature hot spring communities in YNP16,32, including CP13,20. Taxonomically, the community structure did not considerably differ with depth in either the 2016 or 2020 water samples (Fig. 3c), suggesting the presence of generally homogenous communities along depths. Metagenome-assembled-genomes (MAGs) belonging to a group of uncharacterized Sulfolobales were dominant members of all communities in 2016 and 2020, followed to a lesser extent by those of uncultured Acidilobus sp. and Vulcanisaeta sp. (Fig. 3c). These compositions are also consistent with 16S rRNA gene community profiles of CP surface waters from across ~5 months of 201613. In addition, a second Vulcanisaeta sp. MAG was also present in low abundance in all three 2020 water communities along with several low abundance Candidatus Nanopusillus-affiliated MAGs (Table 1). Lastly, a Stygiolobus-like population was also present in the 0 m 2020 water communities (Table 1).

Table 1 Cinder Pool archaeal population and genome information.
Full size table

Although the community compositions were generally similar in 2016 and 2020, their structures varied, with the uncultured Sulfolobales more abundant in the 2020 post-acidification waters, regardless of sample depth (Fig. 3c). Differences in estimated abundance were largely due to the presence of much higher proportions of putatively viral (and/or non-cellular) genomic DNA in the 2016 community metagenomes than those in 2020 (Fig. 3c). Although further characterization of these putative viral contigs was not conducted, the similarity of many of these contigs to those from previously characterized archaeal viruses from YNP (Supplementary Data 3) suggests that the viral load in CP waters in 2016 was overall higher than in waters from 2020. These results are potentially consistent with more stressful conditions for thermoacidophiles (i.e., from higher temperatures and more reducing conditions) in CP waters sampled in 2016 compared to those from 2020. The near homogenous composition of the communities sampled in 2016 and 2020 (Fig. 3c), albeit in decreasing absolute abundance with depth (based on total DNA concentrations) in 2016 and increasing in abundance with depth in 2020 (Fig. 3a, b), suggests that strain level functional variation among populations comprising each assemblage may enable higher fitness in surface waters in 2016 and deeper waters in 2020.

Strain-level shifts in CP populations coinciding with acidification

To further evaluate the involvement of CP populations in hot spring acidification and to identify potential adaptations that may enable higher fitness in surface versus subsurface waters, the MAGs were subjected to phylogenetic and metabolic reconstructions. MAGs for the three primary populations (Uncultured Sulfolobales, Acidilobus sp., and Vulcanisaeta 1) that were present among all three water depth communities and in both 2016 and 2020 communities comprised distinct monophyletic groups corresponding to each year (Fig. 4). The 2016 and 2020 Sulfolobales and Acidilobus populations were not most closely related to corresponding groups in CP waters from the other sampling year, but rather to genomes previously sampled from other YNP hot springs in other studies (Fig. 4). In contrast, the Vulcanisaeta 1 populations were also monophyletic within each year’s communities, but the 2016 and 2020 MAGs were most similar to one another, to the exclusion of Vulcanisaeta MAGs recovered from other springs in YNP. Thus, the altered geochemical conditions of CP that occurred between 2016 and 2020 selected for slightly different strain-level variants of the three pre-dominant populations, in addition to other low-abundance Archaea that were only present in 2020 (Vulcanisaeta 2, Stygiolobus, and Ca. Nanopusillus). The Sulfolobales and Acidilobus MAGs from CP in 2016 and 2020 (in addition to others sampled from CP sediments in 2018 and 2019) were phylogenetically separated by MAGs recovered from numerous other hot springs across YNP and from years ranging from 2008–2018, suggesting that the shifting dominance of strain-level variants in CP between 2016 and 2020 was due to the recruitment and/or proliferation of strains from a pool of metapopulation genotypes distributed across YNP.

Fig. 4: Phylogenomic analysis of metagenome-assembled genomes (MAGs) recovered from Cinder Pool (CP) in depth-resolved samples taken in 2016 and 2020.
figure 4

Maximum-Likelihood phylogenomic trees are shown for the three dominant populations present in both the 2016 and 2020 CP water communities: a Uncultured Sulfolobales, b Acidilobus, and c Vulcanisaeta. A MAG corresponding to a second Vulcanisaeta strain was present in the 2020 water samples and is differentiated by those detected in both sampling years by the “Vulcan2” designation. Phylogenies were constructed from 103 universal archaeal housekeeping single-copy genes (e.g., ribosomal proteins, RNA polymerase subunits, etc.), and available reference genomes from other Thermoproteales. MAGs highlighted in orange are those from the 2020 CP water communities, with the depth given next to the name, while those highlighted in blue were recovered from the 2016 water communities. Bolded, black text indicates MAGs recovered from CP sediments in 2018 and 2019 via our other studies. Asterisks indicate MAGs recovered from YNP springs (as part of this study, or other studies). Bootstrap values were evaluated with 1000 replicates and only those <90 are shown. Scale bars next to individual trees show the expected number of substitutions per site.

Full size image

Genomic plasticity of a subsurface archaeon

The uncultured Sulfolobales that dominated CP waters in 2016 and 2020 belonged to a family-level lineage comprising closely related genomes (all with >95% whole-genome amino acid identity to each other) recovered from numerous acidic or slightly acidic springs in YNP (Figs. 4a and 5)26,32,33,34,35, and that is hereafter referred to as the uncultured YNP Sulfolobales (UYS) group. A total of 22 UYS genomes were subjected to further analyses, comprising six from CP waters in 2016 and 2020, three from CP sediments in 2018 and 2019, four previously published MAGs from other hot springs, and nine other MAGs generated from eight other acidic springs across YNP. While the UYS were present in hot springs with pH ranging from 1.6 to 5.5, the spring pH from where they were identified was highly associated with the phylogenetic placement of individual MAGs (Fig. 5). Specifically, the deepest branching MAGs were generally present in the least acidic springs (pH 4.0–5.5; with the exception of one MAG from Moose Pool [pH 1.6]; inclusive of the 2016 CP MAGs), while the most derived MAG populations were only present in highly acidic springs (pH 2.3–3.7; with the exception of one MAG from Cinder Pool sediments in 2018 [pH 4.0]; inclusive of the 2020 CP MAGs) (Fig. 5). These trends are consistent with the hypothesized trajectory of thermoacidophile evolution towards greater acid tolerance/preference16. Moreover, these trends may also be consistent with niche construction models, wherein H2S/S° oxidation by thermoacidophiles helps construct increasingly acidic environments via acid production, which then concurrently drives adaptation to these increasingly acidic environments16.

Fig. 5: Metabolic reconstructions of uncharacterized Sulfolobales metagenome-assembled genomes (MAGs) recovered from Cinder Pool (CP) and other Yellowstone National Park (YNP) hot springs.
figure 5

Uncultured Yellowstone Sulfolobales (UYS) MAGs are shown in a cladogram on the left in the same order as the phylogenomic tree in Fig. 4a. MAGs recovered from CP waters in 2016 are highlighted in blue and those from 2020 are highlighted in orange. MAGs with superscript labels were recovered from YNP hot spring sediments in previously published studies as follows (a33; b35; c34; d32; e26; note the single MAG generated from a water community of ‘Red Bubbler’ spring). All other MAGs were recovered from YNP hot spring sediments as part of our other studies. The name of the spring is given first, followed by the MAG identifier, and the estimated completeness of the MAG is in parentheses. The year of sample collection and spring pH is given in the column in the middle. The heatmap to the right of the labels shows the presence of genes encoding proteins involved in sulfur cycling within each MAG, with protein names and functions listed at the top of the heatmap. White cells indicate the lack of gene presence, dark blue cells indicate gene presence in the MAG, and light blue cells indicate the presence of homologs identified in the metagenome assemblies that are attributable to the UYS (>90% amino acid identity to other homologs from UYS MAGs), but that was only present on unbinned contig sequences. Proteins are grouped based on their functionalities and associations in complexes. TetH (tetrathionate hydrolase), SQO sulfide:quinone oxidoreductase, SOR sulfur oxygenase reductase, SoxABCD Sulfolobus oxidase, SoxM Sulfolobus oxidase, CbsAB cytochrome b 558/566, SoxLN cytochrome ba complex, DoxBCE Desulfurolobus oxidase, DoxAD/TQOab Desulfurolobus oxidase/thiosulfate-quinone oxidoreductase, HdrAB1C1B2C2 (heterodisulfide reductase), DsrE3 DsrE3 sulfurtransferase, Dld dihydrolipoamide dehydrogenase, LplA lipoate-protein ligase A, LbpA lipoate binding protein A/glycine cleavage system H protein, TusA tRNA 2-thiouridine synthesizing protein A, SreABC sulfur reductase, SAOR sulfite:acceptor oxidoreductase, HcaLS [NiFe]-hydrogenase group 1 g. SoxEFGHI and FoxABCDEFGH (ferrous iron oxidation) gene sets were also investigated, but not identified in any of the MAGs and not shown here for brevity. A complete description of the enzymes/proteins found in individual UYS MAGs is provided in Supplementary Data 4.

Full size image

To assess the potential role of the UYS in sulfur biogeochemical cycling, the metabolic functional potentials of these populations were evaluated in greater detail based on their reconstructed genomes (Fig. 5, Supplementary Data 3). The UYS encoded the capacity for autotrophy via full complements of enzymes involved in the 3-hydroxypropionate/4-hydroxybutyrate cycle (3HP-4HB) (Supplementary Data 4), consistent with the general potential for autotrophy in most other Sulfolobales36. Consistently, the SoxM subunit that has been suggested as a marker for (facultatively) heterotrophic growth of Sulfolobales37 was absent in all UYS MAGs (Fig. 5, Supplementary Data 4). Given that all known Acidilobus and Vulcanisaeta spp. are characterized heterotrophs without known autotrophic capacity38,39, the UYS are likely the sole primary producers in the CP surface and subsurface waters, consistent with their considerable dominance in CP water communities over time.

Also consistent with almost all other Sulfolobales36, the UYS universally encode the ability to reduce O2 via terminal cytochrome oxidases, although not via Sulfolobus oxidase (SoxABCD) complexes that are common among many Sulfolobales36 but rather via Desulfurolobus oxidase complexes (DoxBCE) (Fig. 5, Supplementary Data 4). An additional terminal oxidase complex (CbsAB-SoxLN) was encoded in the 2020 CP MAGs along with several other UYS MAGs from other YNP springs, although homologs of CbsAB-SoxLN were not present in the 2016 CP MAGs or several others recovered from sediments of other hot springs (Fig. 5). Thus, a potentially important metabolic difference between the pre- and post-acidification (2016 and 2020, respectively) CP Sulfolobales was the ability to use different terminal cytochrome oxidase compliments for aerobic respiration. The capacity to use multiple terminal oxidases has been suggested as an adaptation to varying oxygen tensions/availabilities37,40 that likely substantively differed between the low ORP 2016 CP waters and the high ORP 2020 CP waters (Fig. 2c). Consequently, these data point to the ecological succession of UYS strains within CP that are, at least in part, related to strain-level differences in aerobic respiration capacities.

A defining feature of most cultured Sulfolobales is the ability to grow chemolithoautotrophically by coupling the oxidation of sulfur compounds (e.g., S0) to aerobic respiration37. The slow kinetics associated with abiotic oxidation of S0 with O2 at temperatures <100 °C combined with the widespread distribution of Sulfolobales in moderately acidic to acidic high-temperature hot springs has been used to support their role as key mediators of sulfur compound oxidation and hot spring acidification7,16,28. The membrane-bound tetrathionate hydrolase (TetH) and intracellular sulfur oxygenase reductase (SOR) proteins that are present among many Sulfolobales and that mediate the production of S° from tetrathionate (S4O62−) and the disproportionation of S° to H2S and SO32−/S2O3−, respectively36,37, were absent in all UYS MAGs (Fig. 5). The absence of SOR has been suggested to preclude the ability of Sulfolobales to use S° as an electron donor37, although Sulfolobales (i.e., Metallosphaera) are known to use S° as electron donor but do not encode SOR41,42. Further, SOR is not necessary for the growth of acidophilic bacteria (e.g., Acidithiobacillus) via S° oxidation and is thus, likely a supplementary enzyme for S° oxidation43.

Recent experimental studies of acidophilic S°-oxidizing bacteria have indicated that the heterodisulfide reductase complex (HdrABCB2C2)/sulfane-sulfur carrier protein TusA/DsrE3 sulfurtransferase/lipoate-binding protein (LbpA) intracellular system is necessary for growth via S° oxidation44,45. These proteins are highly conserved across Sulfolobales36 and are proposed to operate similarly as in bacterial S° oxidation42,44, although the exact mechanism by which S° is transferred to the sulfur-carrier proteins is currently unclear44. These proteins are nevertheless conserved in nearly all UYS MAGs (and in nearly all Sulfolobales36), with the notable exception of the 2016 CP UYS (Fig. 5, Supplementary Data 4). Given the high level of estimated genome completeness for the 2016 CP UYS (95–100%; Supplementary Data 4), their high level of genomic sequencing coverage (~2200–16,700×; Supplementary Data 4), the presence of most intracellular S° oxidation proteins in most UYS genomes, and the presence of at least some of these genes in MAGs or the unbinned contigs of their corresponding metagenomes suggests that the 2016 CP UYS indeed lacked the genes encoding these proteins and intracellular sulfur oxidation capacity in general. The co-localized HdrABCB2C2 gene structure, along with associated lipoate binding/ligase proteins and sulfur-carrier proteins, that are all syntenous within a 30–40 gene locus present among all UYS MAGs with available data for this genomic region, supports this assertion (Supplementary Fig. 5). Moreover, the core syntenic gene structures share considerable similarity to model Sulfolobales shown to upregulate Hdr expression during autotrophic S° oxidizing growth (e.g., Metallosphaera cuprina, Supplementary Fig. 5)46, along with bacterial S° oxidizers, where this complex has been validated44.

All UYS MAGs also encoded the ability to oxidize H2S via sulfide quinone oxidoreductases (SQO), with many MAGs even encoding multiple copies of sqo (Supplementary Data 4). Further, the coupling of H2S oxidation to O2 reduction could be mediated by thiosulfate–quinone oxidoreductase (TQOab) complexes that are encoded by nearly all UYS. Thus, the UYS generally appear to be able to couple the oxidation of S2O3−, S° (via sulfur carrier proteins), and H2S via TQO, Hdr/TusA/DsrE3/lipoate-metabolizing proteins, and SQO, respectively, to the generation of a reduced quinone pool that can be used to reduce O2 via cytochrome oxidases (Fig. 5, Supplementary Data 4).

Aerobic oxidation of S2O3−, S°, and H2S can contribute to the production of H+ (either through direct oxidation to SO42− and H+ or through the production of intermediates that could then be further oxidized to SO42− and H+), and all of these sulfur compounds are critical focal points in the acidification of hot spring waters47. Thus, aerobic oxidation of these sulfur compounds is likely also occurring in the near subsurface of other regions of YNP where similar geochemical conditions exist. The lack of proteins allowing for intracellular S° oxidation by the 2016 CP UYS suggests that they may have lacked the ability to oxidize sulfur species and were limited to the use of H2S via SQO and S2O3− via TQO, potentially generating S° as a by-product. As such, they may have had a decreased potential to contribute to sulfur oxidation and H+ production compared to 2020 CP UYS. Notably, CP was previously considered anomalous due to its high concentrations of S2O3− which were among the highest in YNP up until 201811. This is due to S2O3− being unstable at pH < 428 and thus not generally being found at high concentrations in acidic hot springs, with the exception of CP11. Consequently, the 2016 CP UYS may have been uniquely adapted to the S2O32− replete historical conditions of CP, but this population may have been quickly outcompeted by the 2020 CP UYS strain upon the disappearance of the molten S° and subsequent (or concomitant) acidification of CP (and thus, the disappearance of abundant S2O3−). Consequently, the 2020 CP UYS strain is uniquely poised to utilize H2S and S° compounds that would have been more available in CP following its acidification, and thus the strain would have outcompeted the 2016 strain based on the availability of sulfur compounds as electron donors.

Conclusions

Prior to autumn 2018–spring 2019, the waters sourcing CP were classified as boiled deeply sourced hydrothermal waters (HB) and were buffered to pH ~4.0 via hydrolysis of a molten S° pool at ~18 to 20 m depth. During this time, reactions involving S° hydrolysis products (S2O32− and SxO62−) that were catalyzed by trace pyrite in cinders formed by gas ebullition through the molten S° would have limited the amount of H2S and O2 available to drive further hot spring acidification10,11. However, between autumn 2018 and spring 2019, the molten S° pool appears to have disappeared, the cinders disappeared, the amount of SO42− increased ~3–5 fold above historical ranges, and the pH of the waters decreased from ~4.0 to ~2.5. While the fate of the molten S0 at depth in CP is not clear, the change in the chemistry of CP and the change in its appearance is intriguingly coincident with the “re-awakening” of the world’s tallest geyser, Steamboat Geyser, which is also present within the NGB. Steamboat Geyser began prolifically erupting in 2018 after years of dormancy, including 109 major eruptions between 2018 and 2020 alone22. The increased geyser activity may have been initiated by deformational uplift in the NGB due to a magmatic intrusion in the years preceding 2018 that was also accompanied by a substantial accumulation of volatiles on the near surface of NGB, perhaps even within a few hundred meters of the surface21. It is possible that the same processes that awakened Steamboat Geyser also impacted the plumbing system sourcing CP and/or disrupted the molten S0 layer at depth, allowing for its acidification after at least ~100 years of geochemical stability. Further geophysical work is warranted to identify potential subsurface connections among the near synchronous changes observed for Steamboat Geyser, CP, and possibly other hot springs in the NGB.

Depth-resolved sampling of CP in 2016 and 2020 identified a member of the Sulfolobales within a group referred to here as the UYS, as the most dominant organism throughout the depths of CP. These data represent the first conclusive evidence to our knowledge of a high-temperature subsurface biosphere in Yellowstone, or in other continental volcanic geothermal regions. In 2020, the abundance of this population increased with depth, suggesting that its realized niche may reflect the more reduced, sulfide-rich, and higher temperature conditions of the CP subsurface. This observation indirectly points to the organism as being active in the subsurface of CP. Metabolic inference indicates that the strain is an autotroph that obtains energy for growth through aerobic oxidation of sulfur compounds that generate acidity. Such observations suggest that UYS may have contributed to the acidification of CP. Other studies in the NGB point to the UYS as dominant microbial populations that potentially inhabit shallow subsurface waters, indicating that they may be widely distributed among acidic subsurface high-temperature environments in YNP. Further work is necessary to identify the cardinal growth and physiological properties of this organism and to assess the extent to which it is active in subsurface waters elsewhere in YNP. Collectively, the results presented herein provide evidence for the existence of a subsurface hydrothermal biosphere in CP, its potential impact on subsurface geochemistry, and the unique geobiological feedbacks that drive geochemical and biological variation in the YNP hydrothermal system. Evidence for a potentially extensive, biogeochemically active, subsurface biosphere in YNP warrants further exploration of other high-temperature continental hydrothermal system subsurface systems to describe their extent, characteristics, activities, and contribution to continental hydrothermal field biogeochemistry.

Methods

Collection of hot spring waters

Water samples were collected from Cinder Pool (CP) in the Norris Geyser Basin (NGB; N 44°43′56.8″, W 110°42′35.1″) using a peristaltic pump (Geotech, Grand Rapids, MI) and silicone tubing (3/16″ ID × 3/8″ OD). A weight attached to a vinyl-coated stainless-steel cable (3/16″ diameter) that was held in place using a two-point hoist was used to position the tubing near the center of the spring. An insulated type K Digi-Sense thermocouple and type K/J Digi-Sense meter (Cole-Parmer, Vernon Hills, IL) were attached at the tip of the tubing to monitor in situ temperature at each depth. The tubing inlet was lowered to sampling depths every 3 m, and water was pumped from each depth at a rate of ~0.42 L min−1 through 0.22 µm polycarbonate filters (Pall Corp, Port Washington, NY) that were contained in an inline filter housing at the outlet of the tubing. After the water collection, the tubing inlet was repositioned to the next sampling depth. Prior to filtering water from each depth, the tubing was flushed for ~10 min to ensure that the waters collected were representative of those from that depth. Triplicate filters for each depth were collected, immediately frozen on dry ice, and then stored at −80 °C until subsequent DNA extractions. Filtered waters were also collected for geochemical analyses (described below) in deionized water-washed Nalgene bottles that were stored at 4 °C prior to analysis upon return to the laboratory.

Physical and chemical measurements

The pH, oxidation–redox potential (ORP or redox potential), and conductivity of spring water sampled from each depth interval were measured in the field with a temperature-compensated model WTW 3100 pH meter (WTW, Weilheim, Germany) a waterproof ORPTestr 10 redox meter (Oakton Instruments, Vernon Hills, IL), and a temperature-compensated EC300 conductivity meter (YSI, Yellow Springs, OH), respectively. Waters were allowed to cool to approximately the same temperature (~50 °C) for conductivity and pH measurements. Ferrous iron (Fe2+) and total sulfide (S2−(T): H2S + HS− + certain metal sulfides) were quantified in the field using Hach reagent kits and a Hach DR/890 field portable spectrophotometer (Hach Company, Loveland, CO). Sulfate (SO42−) and chloride (Cl−) concentrations were also measured using Hach reagent kits and a portable spectrophotometer upon return to the laboratory. The detection limits for Fe(II), S2−(T), SO42−, and Cl− are 0.01, 0.01, 2, and 10 mg/L, respectively. To evaluate trends in CP water geochemistry over time, previously published geochemical data were gathered from published reports (Supplementary Data 1). Where multiple measurements were published for a single year from different reports, values were averaged to represent yearly means. These measurements were necessarily made using different techniques, instruments, and research practices across ~80 years of study, and these differences likely affect the precision and accuracy of measurements. Thus, the trends in these datasets are the primary focus of these comparisons, rather than quantitative geochemical analyses.

Isolation of genomic DNA and shotgun metagenomic sequencing

DNA was extracted from triplicate filters containing biomass collected at depths of 0, 9, and 15 m using a modified phenol–chloroform extraction procedure48. Briefly, filters were re-suspended in 150 mM Tris–HCl buffer (pH 7.5) containing 15 mM EDTA (pH 8.0) and 1 mg mL−1 of lysozyme, followed by incubation at 37 °C for 20 min. After incubation, proteinase K (200 µg mL−1 final concentration) and sodium dodecyl sulfate (1% final concentration) were added to the extract and incubated at 50 °C for 20 min. After incubation, a lysis buffer containing guanidine thiocyanate (1.8 M final concentration) and sarcosyl (0.23% final concentration) was added and the mixture was incubated at 50 °C for 20 min. Cellular and filter debris were removed by centrifugation (9000×g for 20 min). DNA was extracted from the clarified phase, precipitated, and resuspended in molecular grade water as previously described48. The concentration of resuspended DNA was determined using a Qubit 2.0 Fluorimeter (Invitrogen, Carlsbad, CA, USA) and a Qubit dsDNA HS Assay kit (Molecular Probes, Eugene, OR, USA).

Equal volumes of triplicate DNA extracts from 0, 9, and 15 m depths from 2016 were pooled and then subjected to shotgun metagenomic sequencing as part of the Deep Carbon Observatory Census of Deep Life sequencing program. ~62–132 × 106 paired-end (2 × 150 bp) reads (total ~9–19 Gbp/library) were generated on the Illumina HiSeq platform for each metagenomic library. Equal volumes of triplicate DNA extracts from 0, 9, and 15 m depths from 2020 samples were pooled and subjected to paired-end 2 × 150 bp sequencing on the Illumina NovaSeq platform (~70–116 × 106 reads; total of ~10–17 Gbp/library). Sequencing adapters and quality filtering of raw reads were conducted using the TrimGalore v.0.6.0 (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) implementation of the Cutadapt pipeline v.2.1.0 with default parameters. Quality-filtered reads from both datasets were first down-sampled to minimize data redundancy and mitigate uneven population sequencing depth coverage using the BBnorm script (https://github.com/BioInfoTools/BBMap), as described elsewhere25. Down-sampled reads were then assembled using the MetaSPAdes (v.3.14.0) assembler49 with a range of k-mer lengths and default parameters. Assembly contigs were then subjected to MAG binning using the MetaWRAP v.1.3.2 pipeline50. Briefly, read-mapping to the assembled contigs was conducted with the BWA aligner (v.0.7.17)51, followed by binning of contigs >2.5 kbp in length with the metaBAT v.252, MaxBin v.253, and CONCOCT v.1.1.054 algorithms. The “bin_refinement” module of MetaWRAP was used to identify the highest quality bins among the three binning strategies based on CheckM v.1.1.355 estimates of completeness and contamination. Only those MAGs exhibiting >50% estimated completeness and <10% estimated contamination are reported here. In the case of one MAG (the 2020 0 m Sulfolobales MAG), an initially relatively incomplete MAG was first obtained (estimated ~75% completeness). Thus, the nearly complete (~98% estimated completeness) MAG from the 9 m metagenomes was used to recruit contigs from the unbinned 0 m metagenome contigs to produce a more complete 0 m representative MAG.

Bioinformatic and phylogenetic analyses

The relative abundances of the final MAG bin dataset for each metagenome were assessed based on the percentage of quality-filtered reads mapped to the MAG assemblies. Given the dominance of all metagenome assemblies by organisms within the same order and with similar estimated genome sizes (discussed below), reads mapped to the MAGs were used to estimate population abundances, without additional correction for MAG sizes. Protein annotations were generated using the PROKKA v.1.14.5 pipeline56.

Phylogenomic reconstructions were generated to evaluate the evolutionary relationships among taxa. Other closely related reference MAGs and isolate genomes from the Thermoproteales order were gathered from public databases based on BLAST searches of conserved housekeeping genes from the MAGs. A total of 103 archaeal-specific single-copy housekeeping phylogenetic marker genes (e.g., those encoding ribosomal proteins, RNA polymerase subunits, etc.) within the Amphora2 dataset57 were used for phylogenomic reconstructions. Alignments for each gene were constructed using Clustal Omega (v.1.2.4)58, and the 103 alignments were concatenated and subjected to phylogenomic analyses in IQ-TREE59 specifying the LG + G amino acid substitution model and evaluation of node support with 1000 ultrafast bootstrap replicates60. Select functional genes were also subject to phylogenetic reconstruction using the methods described above but using the amino acid substitution indicated as optimal based on Bayesian Information Criterion rankings of 546 models implemented within the model finder plus feature of IQ-TREE. Whole-genome amino acid identity comparisons were calculated with the AAI workflow pipeline within CompareM (v.0.0.23; https://github.com/dparks1134/CompareM).

Metabolic reconstructions for Sulfolobales populations were conducted by querying MAGs for homologs of protein-coding genes using BLASTp searches (considering thresholds of >30% AAI and >60% query coverage) using Sulfolobales references considered to be important for energy conservation or central metabolism36,37. Protein encoding genes identified within a previously published uncultured Sulfolobales ‘Acd1’ MAG34 that was estimated to be 100% complete were also used as queries against the MAGs described here, due to the high level of genomic identity of Acd1 to the MAGs of this study. If CP MAGs did not encode genes found in Acd1 or other closely related MAGs generated in this study, homologs were also queried within all unbinned contigs of the corresponding metagenome assemblies (considering positive identifications with thresholds of >90% AAI and query coverage >60% to those from closely related MAGs) to evaluate whether “missing” genes of MAG assemblies were rather due to incomplete MAG binning.

Abundant short genomic fragments were observed in the 2016 metagenomes and these were queried for the existence of viral genomic signatures to putatively assess the presence of viral (or non-cellular) genomes within the metagenomes. VirSorter (v.2.2.2) was used to identify putative “viral” contigs using the default settings and by only querying those contigs of length >1 kbp. A viral hallmark value of ≥1 or a viral predictive score of >0.5 was used to identify putative viral contigs. Putative contigs comprising <20% of predicted viral genes or >0% predicted “host” sequences were then removed from the putative viral contig results. In addition, putative viral (or non-cellular) contigs were identified among abundant unbinned contigs via BLASTp searches of encoded proteins of the contigs against the Genbank nr database. Specifically, the contigs comprising up to 75% of the total mapped reads of all unbinned contigs for a given metagenome were queried. If a protein within the contig looked plausibly cellular (e.g., it was derived from thermoacidophiles with >40% amino acid identity), then the contig was considered “cellular”. The best protein matches to previously characterized viruses were used to identify putative viral contigs. In addition, contigs that predominantly contained genes without hits within the nr database were assigned as putative viral contigs. The latter was assumed to comprise non-cellular contigs since homologs would be expected to be encoded in described taxa that include numerous closely related genomes from YNP hot springs that are available within the nr database.

Cell density measurements in post-acidification waters

To enumerate cell densities in hot spring waters, samples were collected from 0, 9, and 21 m in 2021 and fixed in formaldehyde in the field (2% v/v final concentration). Triplicate 0.5 mL aliquots of fixed cells from each depth were mixed with 0.1 µL of 4′,6-diamidino-2-phenylindole (DAPI stain; 4 µg mL−1 final concentration) in a sterile microcentrifuge tube. Cells were then incubated for 30 min, collected onto black polycarbonate filters (0.22 µm pore size; Millipore, Billerica, MA, USA), and enumerated with an Evos fluorescence microscope (Thermo Fisher Scientific, Waltham, MA, USA).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.