Phylogenetic and Functional Analysis of Metagenome Sequence from High-Temperature Archaeal Habitats Demonstrate Linkages between Metabolic Potential and Geochemistry

Abstract

Geothermal habitats in Yellowstone National Park (YNP) provide an unparalleled opportunity to understand the environmental factors that control the distribution of archaea in thermal habitats. Here we describe, analyze, and synthesize metagenomic and geochemical data collected from seven high-temperature sites that contain microbial communities dominated by archaea relative to bacteria. The specific objectives of the study were to use metagenome sequencing to determine the structure and functional capacity of thermophilic archaeal-dominated microbial communities across a pH range from 2.5 to 6.4 and to discuss specific examples where the metabolic potential correlated with measured environmental parameters and geochemical processes occurring in situ. Random shotgun metagenome sequence (∼40-45 Mb Sanger sequencing per site) was obtained from environmental DNA extracted from high-temperature sediments and/or microbial mats and subjected to numerous phylogenetic and functional analyses. Analysis of individual sequences (e.g., MEGAN and G + C content) and assemblies from each habitat type revealed the presence of dominant archaeal populations in all environments, 10 of whose genomes were largely reconstructed from the sequence data. Analysis of protein family occurrence, particularly of those involved in energy conservation, electron transport, and autotrophic metabolism, revealed significant differences in metabolic strategies across sites consistent with differences in major geochemical attributes (e.g., sulfide, oxygen, pH). These observations provide an ecological basis for understanding the distribution of indigenous archaeal lineages across high-temperature systems of YNP.

Keywords: archaea; functional genomics; geochemistry; phylogeny; thermophilic archaea and bacteria.

Figure 1

Site photographs of high-temperature chemotrophic systems in Yellowstone National Park (YNP) selected for metagenome sequencing and described in the current study. The sites cover a range in geochemical conditions including sulfur-rich sediments from pH 2.5 (CH_1) to 6.4 (WS_18) and higher oxygen environments where the oxidation of ferrous-iron results in the formation of Fe(III)-oxide microbial mats (OSP_8) [pH and temperature measured on site; yellow arrows indicate sampling locations; all site locations referenced in Table 1].

Figure 2

Frequency plot of the G + C content (%) of random shotgun sequence reads (Sanger) obtained from archaeal habitats in Yellowstone National Park (YNP). Phylogenetic classification of each sequence (∼800 bp) was performed using MEGAN (“blastx”), which shows the predominant populations that contribute to metagenome sequence in these environments (dark-gray = total reads, pink = domain Archaea (shown in sites OSP_8 and JCHS_4 only), yellow = Sulfolobaceae (identity to Sulfolobus sp. also shown), red = Metallosphaera sp., green = Aeropyrum pernix, violet = Euryarchaeota, light-gray = Thermoproteaceae, dark-blue = Caldivirga sp., light-blue = Pyrobaculum sp.). Site WS_18 is shown in Figure A1 in Appendix.

Figure 3

Nucleotide word frequency PCA plots of metagenome assemblies (>5 kb contigs) from high-temperature, archaeal-dominated geothermal environments in YNP. (A) Sites: Crater Hills (CH_1) = gold; Nymph Lake (NL_2) = yellow; Monarch Geyser (MG_3) = green; Joseph’s Coat Hot Springs (JCHS_4) = dark-blue; Cistern Spring (CIS_19) = light-blue; One Hundred Spring Plain “streamer” (OSP_14) = light-red; One Hundred Spring Plain Fe-oxide mat (OSP_8) = red; (B) Order Level: The identical PCA orientation was maintained as presented in A, but with phylogenetic analysis of contigs to the closest reference genome (yellow = Sulfolobales; green = Desulfurococcales; blue = Thermoproteales; purple = Nanoarchaeum; black = unassigned); and (C) Genus Level (yellow = Sulfolobus; red = Metallosphaera; green = Aeropyrum; light-blue = Pyrobaculum; dark-blue = Vulcanisaeta; gray-blue = Caldivirga; purple = Nanoarchaeum; black = unassigned).

Figure 4

Phylogenetic tree (16S rRNA gene) showing major phylotypes identified in assembled metagenome sequence across high-temperature sites dominated by archaea. Entries from metagenome sequence are labeled and colored by site (% nucleotide identity to closest cultivated relative; scaffold ID; dashed boxes indicate the dominant sequence types characterized in this study; bootstrap values based on neighbor joining tree with 1,000 replications).

Figure 5

Nucleotide word frequency PCA scatter plots of non-viral (black) versus viral (red) scaffolds identified within the metagenome assemblies of archaeal-dominated sites (Table 2 provides additional details on the characteristics of viral scaffolds).

Figure 6

Principal components analysis (PCA) of relative gene abundances across seven high-temperature archaeal-dominated chemotrophic communities. The three panels show pairwise plots of the first three principal components (PC1 and PC2 account for 76% of the variation across sites, while PC3 only represents ∼6%). (A) All TIGRFAMs grouped into functional categories, and (B). Only those TIGRFAMs associated with the role category “Electron Transport.” Sites are colored as before: Crater Hills (CH_1) = gold; Nymph Lake (NL_2) = yellow; Monarch Geyser (MG_3) = green; Joseph’s Coat Hot Springs (JCHS_4) = dark-blue; Cistern Spring (CIS_19) = light-blue; One Hundred Spring Plain (OSP_8) = red; Washburn Spring (WS_18) = open.

Figure 7

Hierarchical cluster analysis of relative gene abundances across seven archaeal-dominated communities using all TIGRFAMs grouped into functional categories. Broad TIGRFAM categories include all cellular processes such as regulatory functions, energy metabolism, central C metabolism, mobile elements, transcription, cofactors, and transporters. Data was standardized by functional category before clustering to avoid biasing analysis by a few categories with high gene abundance. Pearson correlation was used as the distance measure for average linkage agglomerative clustering. TIGRFAMs from WS_18 and OSP_8 form separate functional clades consistent with the phylogenetic uniqueness of these sites.

Figure A1

G + C content (%) frequency plot of random shotgun sequence reads (Sanger) obtained from Washburn Springs (WS_18) sediments rich in elemental sulfur and pyrite. Phylogenetic classification of each sequence read (∼800 bp) was performed using MEGAN (“blastx” analysis) (gray = total reads; orange = domain Bacteria; light-pink = domain Archaea; dark-green = Aquificales (primarily Sulfurihydrogenibium-like); light-green/yellow = class Thermodesulfobacteria; bronze = phylum Korarchaeota; white = class Thermoprotei; dark-blue = Pyrobaculum-like; light-blue = Thermofilum-like; other minor assignments with <100 total reads are not shown).

Figure A2

Hierarchical cluster analysis of relative gene abundances in the TIGRFAM role category “Electron Transport” across seven archaeal-dominated geothermal communities. TIGRFAM gene families with low variation across the sites were removed before the clustering to retain ∼50 of the most variable families. Subunits of the protein complexes were only represented by one representative TIGRFAM family. Pearson correlation was used as the distance measure for average linkage agglomerative clustering.

Figure A3

Conservation of open reading frames found in the heterodisulfide reductase (HDR) gene complex of several acidophilic bacteria and archaea known to be utilized during the oxidation of reduced sulfur. Highly syntenous components of the HDR-gene complex are conserved in several of the thermophilic archaea found in YNP, as well as bacteria within the order Aquificales (Takacs-Vesbach et al., 2013). Percent amino acid similarities of the deduced HDR proteins are shown for M. sedula and M. yellowstonensis relative to S. tokodaii. The Sulfolobales Type I and Hydrogenobaculum sp. Y04AAS1 have two copies of drsE [gene names: hdrB1/hdrB2 = heterodisulfide reductase, subunit B (COG2048); hdrC1/hdrC2 = heterodisulfide reductase, subunit C (COG1150); orf2 = hypothetical conserved protein; hdrA = heterodisulfide reductase, subunit A (COG1148); dsrE = peroxiredoxin family protein (COG2210); tusA = sirA family of regulatory proteins; rhd = rhodanese-related sulfurtransferase (COG0607)].