Original Research ARTICLE
Position-Specific Metabolic Probing and Metagenomics of Microbial Communities Reveal Conserved Central Carbon Metabolic Network Activities at High Temperatures
- 1School of Life Sciences, University of Nevada, Las Vegas, NV, United States
- 2Department of Biology, California State University, San Bernardino, CA, United States
- 3Department of Energy Joint Genome Institute, Joint Genome Institute, Walnut Creek, CA, United States
- 4Department of Biological Sciences, Center of Ecosystem Science and Society, Northern Arizona University, Flagstaff, AZ, United States
- 5Nevada Institute of Personalized Medicine, University of Nevada, Las Vegas, NV, United States
Temperature is a primary driver of microbial community composition and taxonomic diversity; however, it is unclear to what extent temperature affects characteristics of central carbon metabolic pathways (CCMPs) at the community level. In this study, 16S rRNA gene amplicon and metagenome sequencing were combined with 13C-labeled metabolite probing of the CCMPs to assess community carbon metabolism along a temperature gradient (60–95°C) in Great Boiling Spring, NV. 16S rRNA gene amplicon diversity was inversely proportional to temperature, and Archaea were dominant at higher temperatures. KO richness and diversity were also inversely proportional to temperature, yet CCMP genes were similarly represented across the temperature gradient and many individual metagenome-assembled genomes had complete pathways. In contrast, genes encoding cellulosomes and many genes involved in plant matter degradation and photosynthesis were absent at higher temperatures. In situ 13C-CO2 production from labeled isotopomer pairs of glucose, pyruvate, and acetate suggested lower relative oxidative pentose phosphate pathway activity and/or fermentation at 60°C, and a stable or decreased maintenance energy demand at higher temperatures. Catabolism of 13C-labeled citrate, succinate, L-alanine, L-serine, and L-cysteine was observed at 85°C, demonstrating broad heterotrophic activity and confirming functioning of the TCA cycle. Together, these results suggest that temperature-driven losses in biodiversity and gene content in geothermal systems may not alter CCMP function or maintenance energy demands at a community level.
Temperature is a primary driver of taxonomic diversity of microbial communities (Cole et al., 2013; Sharp et al., 2014; Sunagawa et al., 2015; Power et al., 2018), with an observed maximum at 25°C and decreasing diversity as temperature becomes more extreme (Sharp et al., 2014). This temperature-driven diversity gradient provides an opportunity to assess the consequences of decreasing taxonomic diversity on ecosystem function (Swingley et al., 2012). On the one hand, several studies have suggested that taxonomic and functional diversity are not closely linked. For example, the Human Microbiome Project Consortium (2012) found core functional diversity to be shared by microbial communities in different locations on the human body, despite vast differences in taxonomic composition. Similarly, Sunagawa et al. (2015) showed that core functional diversity was conserved in the global ocean microbiome, despite spatial variability in taxonomic diversity, and that >73% of the ocean and human microbiome core functional gene content was shared. The existence of an ecosystem-independent core functional gene assemblage suggests conservation of ecosystem function across diverse microbial communities and physiochemical conditions (Sunagawa et al., 2015). On the other hand, the upper-temperature limit for photosynthesis is ~73°C (Brock, 1967; Cox et al., 2011; Boyd et al., 2012) and the upper-temperature limit for chemolithotrophic nitrite oxidation may be ~65°C (Lebedeva et al., 2005; Sorokin et al., 2012; Edwards et al., 2013), suggesting temperature can be a strong driver of not just taxonomic diversity but also functional diversity.
Contained in the global ocean and human microbiome core gene assemblages are genes necessary for glycolysis, the pentose phosphate pathway (PPP), and the tricarboxylic acid cycle (TCA) (Sunagawa et al., 2015). These “central carbon metabolic pathways” (CCMPs) are intimately linked to biological energy-generating reactions (catabolism) and to the synthesis of biomass precursors (anabolism) in heterotrophic organisms. If there is an ecosystem-independent core functional assemblage for microbial communities, then these pathways should be present in microbial communities inhabiting extreme environments as well. While metagenomic studies have demonstrated the conservation of ecosystem function at the genomic level in diverse ecosystems, these studies failed to address the activities directly.
The activity of CCMPs in microbial communities can be assessed using position-specific 13C-labeled substrates and modeling metabolic flux using 13CO2 production rates (Dijkstra et al., 2011a,b,c, 2015; Hagerty et al., 2014; Figure 1). The use of position-specific 13C-labeled compounds (i.e., isotopomers) provides greater information than uniformly 13C-labeled substrates, allowing inferences to be made about the flux of carbon through CCMPs and the dynamics of catabolic and anabolic reactions (Dijkstra et al., 2011a,b; Leighty and Antoniewicz, 2013). For studying microbial community activity, this approach is superior to the classic 13C metabolic flux analysis based on amino acid labeling patterns (Wiechert et al., 2001), largely due to the difficulties of obtaining an isotopic steady state for all members of the community (Zamboni et al., 2009).
Figure 1. Simplified demonstration of decarboxylation events with associated position-specific isotopomers used in this study. Locations of decarboxylation events in central carbon pathways (A), and the associated C position in the 13C-labeled compounds used in these experiments (B). C1 from glucose is lost as CO2 in the pentose phosphate pathway with the synthesis of ribulose 5-phosphate by 6-phosphogluconate dehydrogenase or after multiple rounds in the TCA cycle. C1 from pyruvate is lost as CO2 in the formation of acetyl-CoA by pyruvate dehydrogenase. C2 and C3 of pyruvate are lost in multiple cycles of the TCA cycle. C1 from acetate is lost as CO2 in the 2nd cycle of the TCA cycle, while half of C2 is lost in each subsequent cycle. All CO2 producing reactions can produce 13CO2 from uniformly labeled 13C glucose. This example assumes a glycolytic flux through the central carbon pathways. TCA, tricarboxylic acid cycle; PPP, pentose phosphate pathway.
The principle governing using 13C-labeled isotopomers is that if a substrate is used solely for catabolic purposes, then carbon (C) from all positions in the substrate will be oxidized to CO2 in the same ratio they are found in the substrate (i.e., 1-13C: 2,3-13C (13C1:13C2, 3) pyruvate ratio equals 0.5). However, if anabolic reactions are taking place, then this ratio deviates from the hypothetical “catabolic only” ratio. For example, in the complete oxidation of pyruvate through CCMPs, C1 is released as CO2 as a result of pyruvate dehydrogenase or pyruvate-ferredoxin oxidoreductase activity, while C2 and C3 are oxidized to CO2 in the TCA cycle (Figure 1). Similarly, C1 of glucose is oxidized to CO2 during the formation of ribulose-5-phosphate via 6-phosphogluconolactonase in the PPP, while the other five carbons are oxidized by pyruvate dehydrogenase and the TCA cycle. If glucose is metabolized through glycolysis, then C1 of glucose can also be oxidized in the TCA cycle, but only after at least the third cycle, providing a greater opportunity for C1 to be sequestered in biomass. Thus, the relative rate of 13CO2 production from C1 of glucose [relative to uniformly 13C-labled (13CU) glucose] would be high when PPP is active and low during conditions of reduced PPP activity, assuming metabolites are being used for biomass production (Dijkstra et al., 2015). By comparing rates of 13CO2 production from specific isotopomers (e.g., 13C1 vs. 13C2, 3 pyruvate), information is gained on the relative activity of CCMPs and the catabolic and anabolic use of these compounds.
In the present study, we focused on a single terrestrial geothermal system, the Great Boiling Spring (GBS) system, located in northern Nevada, U.S.A. GBS is a circumneutral pH, NaCl-type geothermal system that harbors well-studied microbial communities rich in novel and yet-to-be-cultivated microorganisms (Dodsworth et al., 2012, 2013; Rinke et al., 2013; Becraft et al., 2016, 2017; Eloe-Fadrosh et al., 2016). Within this system, taxonomic diversity has been shown to be inversely proportional to temperature, concomitant with a shift from Bacteria to Archaea (Cole et al., 2013). Native communities in GBS have been shown to utilize a variety of organic electron donors for aerobic respiration (Murphy et al., 2013), and denitrification in GBS appears to be heterotrophic (Dodsworth et al., 2013). These studies, along with stable isotope natural abundance data from Yellowstone National Park (Schubotz et al., 2015; Jennings et al., 2017), point to the importance of both heterotrophy and autotrophy for high-temperature microbial communities. Yet, to our knowledge, the activity of CCMPs within terrestrial geothermal systems has not been addressed directly. This system provides an opportunity to study the functional consequences of temperature-driven losses in taxonomic diversity while minimizing other variables, such as water chemistry, photosynthetically active radiation, and geographic distance.
To examine the relationships between temperature-driven taxonomic diversity and the functioning of CCMPs, we combined 13C-based metabolic probing with analysis of 16S rRNA gene amplicons, shotgun metagenomes, and metagenome-assembled genomes (MAGs) derived from the same samples. This system was then used to test several hypotheses about relationships between temperature-driven losses in taxonomic diversity and ecosystem function: (1) genes associated with photosynthesis and degradation of complex organic matter will decrease at high temperature; (2) PPP gene abundance and ecosystem function will decrease at higher temperature, due to a relative increase in Archaea (Bräsen et al., 2014); (3), genes for other CCMPs and catabolic carbon use will be unaffected by temperature (Valentine, 2007; Sabath et al., 2013; Bräsen et al., 2014); (4) catabolic carbon use will increase with temperature, due to an increase in maintenance energy demand (Tijhuis et al., 1993; Price and Sowers, 2004); and (5) diverse TCA cycle intermediates will be used catabolically at high temperature.
Materials and Methods
Three sediment sample locations within GBS (N40°39′41″ W119°21′58″) and one location from GBS 19 (~30 m from GBS) were selected to encompass a temperature gradient from 60–95°C. Sites were designated GBS 95, GBS 85, GBS 70, and GBS 60 (Figure S1; Table S1), and roughly correspond to previous study sites (Figure S1; Cole et al., 2013). GBS and GBS19 share a single deep subterranean source and have similar water chemistry (Table S1; Anderson, 1978).
Field Measurements and Water Chemistry Analysis
Prior to water and sediment collection, temperature, conductivity, and pH were measured using a LaMotte pH5 Series pH/temperature meter (LaMotte, Chestertown, MD, USA). For major cation and anion analysis, water was filtered through 0.2 μm polyethersulfone (PES) filters (Pall Corp., USA) into 15 mL polypropylene tubes and frozen on dry ice. For trace element analysis, water was filtered through a 0.2 μm polyethersulfone (PES) filters (Pall Corp., USA) into acid-washed 150 mL Nalgene® bottles spiked with 400 μL 10% HNO3 (OmniTrace Ultra). Trace element samples were shaken and degassed three times, then stored at room-temperature until analysis. Field blanks, consisting of sterile Nanopure water processed in tandem with environmental samples, were collected at the time of sampling.
Ion chromatography (IC; Dionex DX-500 chromatography, Dionex Co., USA), direct current plasma optical emission spectrometry (DCP-OES; Beckman, USA), and inductively coupled plasma mass spectrometry (ICP-MS; Varian 820 quadrupole, Agilent Technologies Inc., USA) analysis for major anions, cations, and trace elements were performed as described by Huang et al. (2013).
DNA Extraction and Illumina 16S rDNA Sequencing
The modified FastDNA Spin Kit for Soil (MP Biomedicals, Solon, OH, USA) (Dodsworth et al., 2011) was used for sediment DNA extraction (collected as described below) for both 16S rRNA gene and metagenomic sequencing. DNA was shipped to Micro-Seq Enterprises (Las Vegas, NV, USA), where the amplification and sequencing of the V4 region of the 16S rRNA gene for Bacteria and Archaea was performed using the Illumina MiSeq platform (San Diego, CA, USA), as described in Kozich et al. (2013) and a forward primer designed to increase coverage of Archaea (Hou et al., 2013). Sequences have been deposited to the NCBI Sequence Read Archive under project SRP059341 (GBS 60, SRX1055762; GBS 70, SRX1055763; GBS 85, SRX1055764; GBS 95, and SRX1055765).
Sequences were trimmed, assembled, ambiguous bases were removed and classified using a Bayesian classifier using Mothur according to the standard operating procedure provided by Kozich et al. (2013). Chimeras were identified in the aligned sequences using UCHIME with the Silva SEED database (v119). Following removal of chimeras, sequences were binned into operational taxonomic units (OTUs) above sequence identity of 97.0%. Taxonomic assignment was according to SILVA (v119) (Pruesse et al., 2012; Quast et al., 2013; Yilmaz et al., 2014). Manual curation was performed on several unclassified sequences using BLASTn, Greengenes, EzBioCloud (Yoon et al., 2017), and sequence cut-offs suggested by Yarza et al. (2014) for identification and naming of sequences belonging to candidate phyla.
Community statistics were calculated using rarefaction tool kit (RTK) (Saary et al., 2017) and phyloseq (McMurdie and Holmes, 2013). The OTU abundance table and taxonomy files from mothur were imported into phyloseq, alpha diversity measurements were carried out using “estimate_richness” function. For RTK, samples were rarefied based on GBS 95, which had the lowest number of reads (17,181). OTUs used in the SILVA assignments and the alpha diversity analyses were grouped at 97% sequence identity.
Metagenomic Illumina Sequencing
Detailed descriptions for the processing, sequencing, and assembly of the metagenomes are publicly available from the Integrated Microbial Genomes and Microbiomes platform (IMG/M; https://img.jgi.doe.gov/m/) (Chen et al., 2019), under Gold Analysis Project IDs Ga0197142, Ga0197143, and Ga0197144 (GBS 60, GBS 70, and GBS 85, respectively). Briefly, isolated DNA for metagenomic sequencing was sequenced by the Joint Genome Institute (JGI) using shotgun Illumina HiSeq 2500-1TB, assembled using SPAdes version: 3.10.1, and annotated using the IMG Annotation Pipeline v.4.15.1 (Huntemann et al., 2016). JGI-identified 16S rRNA gene sequences contained in each metagenome were classified using Silva (v128). Metagenome sequencing from GBS 95 was not possible due to the low yield of DNA.
Metagenome-Assembled Genomes (MAGs)
Genome binning was performed using MetaBAT (v0.32.4; Kang et al., 2015) with a 3,000 bp minimum contig cutoff, contig coverage information, and parameter “-superspecific” for maximum specificity.
Bin completeness and contamination were estimated using CheckM (Parks et al., 2015). Bins were deemed to be of high-quality if they were >90% complete, <5% contaminated, and contained a full rRNA operon and >18 tRNAs (Bowers et al., 2017). High-quality bins were analyzed using the Genome Taxonomy DataBase Toolkit and tools contained therein (GTDB-Tk, v0.1.1) for gene calling, protein sequence prediction, and taxonomic identification (Hyatt et al., 2010; Matsen et al., 2010; Eddy, 2011; Jain et al., 2018; Parks et al., 2018). 16S rRNA gene amplicons were matched with metagenome bins using BLAST+ (Camacho et al., 2009) in the web-based Galaxy platform (Cock et al., 2015). 16S rRNA gene amplicons were queried against each metagenomic bin using Megablast (Zhang et al., 2000; Morgulis et al., 2008).
Assessing Metabolic Potential and Pathway Completeness
Protein sequences for high-quality bins and metagenomes were submitted individually to MAPLE (Metabolic and Physiological potential Evaluator, v2.3.1) (Takami et al., 2012b, 2016; Arai et al., 2018) for the evaluation of Kyoto Encyclopedia of Genes and Genomes (KEGG) functional module completion ratios (Takami et al., 2012b). Within MAPLE, the NCBI BLAST search engine was used with the single-direction best hit annotation method for KEGG Orthology (KO) assignment, using all organisms in the KEGG database. Whole-community module completion ratios were considered biologically feasible at Q < 0.5 (https://maple.jamstec.go.jp/maple/maple-2.3.1/help.html).
NOISeq Differential Abundance and KEGG Mapping
KOs identified using the JGI IMG/M pipeline were analyzed using NOISeq-sim (Tarazona et al., 2011). To determine relative abundance of KOs, protein sequences with KEGG Orthology (KO) annotations were tabulated using the JGI-provided estimated gene copy number, normalized to unassembled metagenome size (GBS60, 23,572,233,708 bp; GBS70, 18,720,210,806 bp; and GBS85, 16,850,842,886 bp) using counts per billion bases, and log2-transformed using NOISeq v2.22.0 (Tarazona et al., 2011, 2015) in R v3.4.3 (R Core Team, 2017). To determine KO allele richness, genes with KO annotations were tabulated based on the JGI-provided gene count, normalized to assembled metagenome size (GBS 60, 461,589,717 bp; GBS 70, 297,862,735 bp; and GBS 85, 80,294,486 bp) using counts per billion bases, and log2-transformed as above. Differential abundance and richness analysis was performed using the NOISeq-sim function with the developer-recommended parameters for analyses of simulated replications with no pseudocount added to zero-count hits (Tarazona et al., 2011, 2015). The log2-transformed normalized abundance and normalized richness values were then mapped to relevant KEGG pathway diagrams (Kanehisa and Goto, 2000) and significant differences in KO abundance were noted, using a custom algorithm that combines the outputs of the KEGG reconstruct pathway and KEGG color pathway tools (Kanehisa and Goto, 2000; Kanehisa et al., 2012, 2017). The NOISeq-sim recommended probability value of 0.9 was used to indicate significant differences between any pairwise comparison (Tarazona et al., 2011, 2015). All possible pairwise comparisons were made.
Carbohydrate-active enzyme (CAZyme) (Lombard et al., 2014) annotation was performed on JGI IMG/M protein sequences using hmmer version 3.1b2 (Eddy, 1998) with the dbCAN database (Yin et al., 2012) as a reference. To determine relative abundance of CAZymes, CAZyme HMM hits for each of the metagenomes were tabulated using scaffold average coverage and normalized to unassembled metagenome size using counts per billion bases. To determine CAZyme allele richness, CAZyme HMM hits were tabulated without considering scaffold read-depth and then normalized to assembled metagenome size. Differential abundance and richness of the normalized HMM hits were analyzed as above using the NOISeq-sim function (Tarazona et al., 2011, 2015). Raw reads and assembled metagenomic contigs, as well as functional predictions, are available on the JGI IMG/M platform (Chen et al., 2019) under taxon object IDs 3300020153, 3300020139, and 3300020145 (GBS60, GBS70, and GBS85, respectively).
13CO2 Production Rate Experiments
Hot spring sediment slurry microcosms (top ~1 cm) were incubated in situ to monitor the production rate of 13CO2 from position-specific 13C-labeled isotopomer pairs and uniformly 13C-labeled compounds over an 8-h period. Subsamples of the sediment slurry were aseptically transferred to 2 mL centrifuge tubes and immediately frozen on dry ice for later DNA extraction for 16S rRNA and metagenome work, or were dispersed for the 13CO2 production rate experiments, as detailed below. Slurries were kept near in situ temperatures by placing the bottles in the hot spring during preparation and incubation of samples.
To prepare microcosms, ~20 mL of sediment slurry was aseptically dispersed into sterile, pre-warmed 165 mL serum bottles using a modified 10 mL Oxford BenchMate™ pipette (St. Louis, MO, USA). After slurry addition, serum bottles were capped and sealed and then 1 mL of 13C-labeled substrate was added using a 1 mL syringe and 25 G needle. Each 13C substrate incubation was performed in triplicate except for “killed,” mercury- or glutaraldehyde-treated controls (described below), which were performed in duplicate to limit sample numbers. Microcosms were secured in a wire cage, covered with foil, and completely submerged in the hot spring. Temperature was monitored throughout the experiment using a Traceable® Ultra Waterproof Thermometer (VWR, Radnor, PA, USA) fixed to the wire cage. At the end of the experiment, microcosms were kept on wet ice until sediments could be dried and weighed in the lab. ~10 mL headspace samples were taken ~2, 4, and 8 h after 13C-labeled substrate addition. Atmosphere samples (n = 3) were collected onsite to determine a time-0 13CO2 delta value. Samples were stored in collection syringes until analysis, as described below.
Isotopomer pairs consisted of filter-sterilized solutions (21.4 μmol substrate-C mL−1) of sodium pyruvate (1-13C and 2,3-13C), sodium acetate (1-13C and 2-13C), and glucose (1-13C and uniformly (U) 13C-labeled) (99 atom fraction %; Cambridge Isotope Laboratories, Andover, Massachusetts). For the 85°C site, uniformly 13C-labeled citrate, L-serine, L-cysteine, L-alanine, and succinate (99 atom fraction %; Cambridge Isotope Laboratories, Andover, Massachusetts) were also used at 0.05 mg/mL to evaluate the diversity of organic compounds catabolized by the community. Mercuric chloride and glutaraldehyde were used as poisons [1.0 mM and 2.0% (vol/vol), respectively], individually and in combination, for monitoring abiotic decarboxylation of uniformly 13C-labeled substrates. Mercury was tested at all four sites, while glutaraldehyde was tested at 60 and 95°C sites only.
Headspace samples were analyzed for 13CO2 on a Picarro 2101-i CO2 and CH4 isotope spectrometer (Picarro Inc, Sunnyvale, CA, USA) by first spiking samples with a high concentration of CO2 and then diluting with CO2-free air. This step was necessary to ensure that CO2 concentrations were within the working range of the Picarro (300–2,000 μmol mol−1) and to provide enough gas volume for a 10-min. sample analysis. Spikes were performed by filling a Mason jar (473 mL) with atmosphere and injecting 10 mL of pure CO2 through a rubber stopper inserted in the Mason jar lid. For each sample, a volume of gas from the CO2-spiked Mason jar was added at a volume matching the microcosm gas sample volume, resulting in a 1:1 ratio of sample and CO2-spiked atmosphere. A maximum of four CO2-spiked aliquots were taken from the Mason jar before the jar was opened, vented, and refilled with CO2-spiked atmosphere. Both CO2-spiked atmosphere and samples were injected into a Tedlar® bag (0.5 L) and CO2-free gas was used to increase the mixed sample volume. Data from the Picarro were recorded as 30-second averages of δ13CO2 over a period of near-constant delta readings. Data were expressed as atom percent excess (APE) per hour (relative to time zero measurements).
Field Measurements and Water Chemistry
pH (7.17–7.52), conductance (7.47–7.68 mS), major ions, and trace elements were similar at all four sites along the temperature gradient (Table 1, Table S1). In situ temperatures, measured by tracing thermometers during the experiments, varied by ~4.0°C around the target sample temperatures used as site identifiers (i.e., GBS 60, GBS 70, GBS 85, and GBS 95; Table 1).
Analysis of filtered 16S rRNA gene amplicons showed that observed richness ranged from 113 to 491 species-level OTUs (3%), while Chao1 estimated richness values ranged from 291.75 to 833.16 (Table 2). Simpson, inverse Simpson, and Fisher's alpha diversity ranged from 0.69 to 0.96, 3.27 to 26.05, and 16.01 to 87.53, respectively (Table 2). Shannon diversity ranged from 1.61 to 4.01 (Table 2). All diversity metrics decreased as temperature increased, although the GBS 85 site consistently demonstrated the lowest values (Table 2; Table S2).
Both 16S rRNA gene amplicons and metagenomes showed that Archaea were more abundant at high temperatures, and phototrophs were not detected at high-temperature sites (Figure 2; Table 3). The highest proportion of rare sequences (<1% relative abundance) was found at low-temperature sites (Figure 2; Table 3). Many abundant phyla contained limited numbers of OTUs (Table 3; Table S4), and many OTUs were unique to a single temperature (Table 3; Table S3). At all temperatures, many of the abundant OTUs represented high-level taxa (i.e., phyla, classes, and orders) with no cultivated members.
Figure 2. The relative abundance of phylum-level groups for each sample based on 16S rDNA Illumina Tags (iTags) of the V4 region compared with metagenome-derived 16S rDNA. Colored rows show the relative abundance of phyla present at >1% and correspond to the site label on the left. Unclassified sequences are represented as a single group; see Table S3 for information on individual unclassified groups and their abundance. Gray dashed bar indicates GBS 19 is a separate geothermal pool from the main pool. Metagenomic data from GBS 19 is not available due to low DNA yield.
Chloroflexi (~22%) and Cyanobacteria (~16%) were the most abundant groups at the GBS 60 site (Figure 2; Table 3; Table S3), including the genera Roseiflexus (OTU0014), Chloroflexus (OTU0023), and Caldilinea (OTU0049).
The GBS 70 site had four phylum-level groups >10% relative abundance: Chloroflexi (~18%), “Aigarchaeota” (~16%), “Acetothermia” (~11%), and Thaumarchaeota (~11%) (Figure 2; Table 3; Table S3). Known cultivated or uncultivated genera included the chemoheterotroph Thermoflexus hugenholtzii (OTU0016) (Dodsworth et al., 2014), the “Aigarchaeota” genera Candidatus Calditenuis aerorheumensis (OTU0001) (Beam et al., 2016) and Candidatus Caldiarchaeum subterraneum (OTU0005) (Nunoura et al., 2011), and the chemolithotrophic ammonia-oxidizing archaeon Candidatus Nitrosocaldus yellowstonii (Table S3; De La Torre et al., 2008).
The GBS 85 site was dominated by only two major groups, “Aigarchaeota” (~53%) and bacterial candidate phylum “Fervidibacteria” (~24%) (Rinke et al., 2013). Candidatus Calditenuis aerorheumensis (OTU0001) accounted for 48% of total abundance. “Fervidibacteria” OTU0002 was nearly identical (99% ID) to Candidatus Fervidibacter sacchari (Rinke et al., 2013). OTU0024 represented Candidatus Calescibacterium nevadense (Rinke et al., 2013; Becraft et al., 2016).
The GBS 95 site was dominated by Crenarchaeota (~42%) and “Aigarchaeota” (~25%). Most Crenarchaeota belonged to the order Desulfurococcales. “Aigarchaeota” consisted of a single OTU (0003) that was 99.4% identical to sequences found within Group 4 “Aigarchaeota” (Hedlund et al., 2015), and a MAG retrieved from Jinze Spring in Tengchong, China (Hua et al., 2018).
Metabolic Potential in Metagenomes
KEGG module completion ratios (MCR) showed that the metagenomes had complete or near-complete CCMPs, with most modules being biologically feasible (Q < 0.5) (Figure 3). The only modules in central carbohydrate metabolism that were not biologically feasible were the Entner-Doudoroff pathway (KEGG module M00008) and the semi-phosphorylative Entner-Doudoroff pathway (M00633) in the GBS 60 and GBS 85 metagenomes (75% MCR, Q > 0.5), despite having 100% MCR (Q-value, 0) in the GBS 70 metagenome (Figure 4; Table S5). This was due to the absence of phosphogluconate dehydratase (K01690, EC 220.127.116.11). An alternative module for the semi-phosphorylative Entner-Doudoroff pathway (M00308) was biologically feasible in all metagenomes (Table S5).
Figure 3. Selected MAPLE-derived module completion ratios for assembled metagenomes and GTDB-Tk taxonomically defined, high-quality MAGs. Heat map indicates the degree of a module's completion, ranging from 0 to 100%, for metagenomes or MAGs. Module completion ratios have been rounded to the nearest whole number. *The associated Q-value is below 0.5 and is deemed biologically feasible. The first two-digit number in the Bin ID column indicates the metagenome the bin is from. Letters in parentheses for taxonomic classification indicate p, phylum; c, class; o, order; f, family; g, genus; s, species. Numbers in parentheses for individual modules indicate the number of components included in the module. See Table S4 for full GTDB-Tk classification of MAGs. See Table S5 for all modules, KEGG pathway IDs, and number of components comprising a single module.
Figure 4. Integration of metagenome data and isotopomer experiments, focusing on key decarboxylation reactions. (A) Differential abundance and richness plots for KOs in GBS 60 and GBS 85 metagenomes. Large dots indicate KOs involved in decarboxylation reactions involving the C1 position of pyruvate, acetate, or glucose. Small dark dots indicate KOs that were significantly differentially abundant or rich. All data are shown in Table S6. (B–D) KOs involved in the production of CO2 from the C1 position of substrates. The left side shows simplified KEGG pathways with log2 abundance and richness heatmaps for GBS 60, GBS 70, and GBS 85. The middle shows 13CO2 production rates from isotopomer pairs for the GBS 60 and GBS 85 sites. APE, atom percent excess. The right demonstrates the catabolic and anabolic use of isotopomer pairs. Graphs indicate the ratio of the slopes (n = 9) for APE of 13C of CO2 from the isotopomer microcosms. Dashed line, ratio for complete oxidation; error bar, standard error of the mean; *0.001; **0.0001; ***0.00001 significantly different from complete catabolic use (p < 0.05; one-sample t-test); shared letters within an isotopomer group indicate that they are not distinguishable (ANOVA, Tukey HSD, and Welch t-tests).
The number of complete KEGG modules for carbon fixation pathways decreased as site temperature increased, with GBS 60 having 10 pathways with 100% MCR, GBS 70 having seven, and GBS 85 having six, out of a total of 15 modules (Figure 3; Table S5). The reductive citrate cycle (Arnon-Buchanan cycle) and dicarboxylate-hydroxybutyrate cycle had 100% MCR in all metagenomes. The reductive acetyl-CoA pathway (Wood-Ljungdahl pathway) and the 3-hydroxypropionate bi-cycle both had 100% MCR in GBS 60 and GBS 70 but were not biologically feasible in GBS 85. The reductive pentose phosphate cycle (Calvin cycle) and photosynthetic complexes and pathways had 100% MCR in GBS 60, but were not biologically feasible in GBS 70 (90.9%) and GBS 85 (90.9%). The hydroxypropionate-hydroxybutyrate cycle was not biologically feasible in any of the metagenomes.
Dissimilatory sulfate reduction, thiosulfate oxidation by the SOX complex, denitrification, and nitrification had 100% MCR in all three metagenomes, while no components for nitrogen fixation were found (0% MCR) in any of the metagenomes. MCRs for structural complexes associated with sulfate (M00185) and sulfonate (M00436) transport, and a urea transport system (M00323), all decreased as temperature increased (Table S5). The potential ability to use methanol as a substrate for methanogenesis (M00356) decreased with increasing site temperature (Table S5). Pectin degradation (M00081) had 100% MCR in GBS 60 but 0% MCR in GBS 70 and GBS 85 (Table S5). Modules involved with the biosynthesis of secondary metabolites were largely absent from all metagenomes (Table S5).
Metabolic Potential in MAGs
A total of 167 MAGs were obtained from the three metagenomes of which 33 were considered high quality (Bowers et al., 2017). High-quality bins were taxonomically identified with the GTDB-Tk: 13 from GBS 60, 15 from GBS 70, and 5 from GBS 85 (Figure 3; Table S4). Most contained a 16S rRNA gene sequence identical to a 16S rRNA gene from the amplicon analysis (Table 3; Table S4). High-quality MAGs spanned eight bacterial phyla, one candidate bacterial phylum, and two archaeal phyla, according to the Genome Taxonomy Database (Parks et al., 2018): Acidobacteriota (4 bins), Armatimonadota (1), Bacteroidota (6), Chloroflexota (6), Cyanobacteriota (1), Desulfobacterota (3), Firmicutes (1), Nitrospirota (1), WOR-3 (2) (Baker et al., 2015), Crenarchaeota (6) (including “Aigarchaeota” and Thaumarchaeota), and Hadesarchaeota (1).
Bin 40 from GBS 60 (60_40), belonging to the genus Gloeomargarita, was the only bin to show 100% MCR for photosystems I and II as well as NAD(P)H:quinone oxidoreductase (M00145) (Figure 3; Table S5). Bin 60_40 also had 100% MCR for the reductive pentose phosphate cycle (Calvin cycle). Bin 12 from GBS 60 (Chloroflexus sp.) was the only bin to show 100% MCR for anoxygenic photosystem II. No components were found for anoxygenic photosystem I in any of the bins (0% MCR). All other bins had 0% MCR for all photosystems.
Glycolysis, pyruvate oxidation to acetyl-CoA, the TCA cycle, and the 1st and 2nd carbon TCA cycle oxidations were biologically feasible in most bins (Figure 3). In bins where pyruvate oxidation to acetyl-CoA was not found, gluconeogenesis was biologically feasible (Figure 3). The full PPP (M00004), including both oxidative (M00006) and non-oxidative phases (M00007), were biologically feasible in only five bins (all modules, 100% MCR) (60_12, 60_15, 60_26, 60_34, 85_2), while the non-oxidative phase (M00007) had 100% MCR in most bins (24 bins) (Figure 3). The archaeal PPP (M00580) was biologically feasible in six bins (60_51, 70_34, 70_36, 85_6, 85_8, 85_9), with three belonging to the order Caldiarchaeales [Genus 2 Hedlund et al., 2015]. The Entner-Doudoroff pathway (M00008) was not biologically feasible in any of the bins (Figure 3). The non-phosphorylative Entner-Doudoroff pathway (M00309) was biologically feasible in only one bin (85_6) (Table S5), belonging to the Crenarchaeota family Ignisphaeraceae.
KO and CAZyme Abundance and Richness in Metagenomes
A total of 5,874 KOs and 358 CAZymes were identified in the three metagenomes from the 60, 70, and 85°C sites (Tables S6, S7). GBS 60 contained 5,442 KOs and 329 CAZymes, with total estimated gene copies of ~9,001,975 and ~776,603, respectively. GBS 70 contained 4,852 KOs and 299 CAZymes, with total estimated gene copies of ~7,677,420 and ~629,856, respectively. GBS 85 contained 2,995 KOs and 177 CAZymes, with total estimated gene copies of ~8,403,495 and ~540,463, respectively. Thus, the abundance and richness (i.e., the total number of different genes in each KO) of both KOs and CAZymes both decreased as a function of increasing temperature, in step with the pattern of taxonomic richness (Figure 4A; Tables S6, S7).
After normalization, NOISeq-sim analyses yielded 1,294 total KOs and 62 CAZymes that were differentially abundant for any pairwise comparison (Tables S6, S7), and 2,835 KOs and 134 CAZymes that were differentially rich (Tables S8, S9), with a probability of 0.9 or greater of being in any pairwise comparison. Of these 2,835 KOs, 1,047 were significant in both richness and abundance, 247 were unique to abundance, and 1,788 were unique to richness (Figure S2; Table S10). Of the 134 CAZymes, 48 were significant in both richness and abundance, 15 CAZymes were unique to abundance, and 86 CAZymes were unique to richness (Figure S2; Table S10). Figure 4A shows an example of the difference in abundance and richness of KOs for the GBS 60 to GBS 85 comparison. The plot shows that KOs were more abundant and richer in the GBS 60 metagenome (865 and 1,748), compared with the GBS 85 metagenome (77 and 686).
Many KOs significantly more abundant in GBS 70 and GBS 85 than in GBS 60 were archaeal KOs involved in information processing (e.g., transcription, translation, genome replication, recombination, and DNA repair) (Table S6). For example, archaeal KOs associated with tRNA synthesis, DNA repair, DNA and RNA polymerases, ribonucleases, and ribosomal proteins were more abundant at higher temperature sites. KOs associated with archaeal CRISPR proteins were also more abundant at higher temperature sites (e.g., K19074 and K07725).
Richness for KOs involved in information processing generally followed the same pattern as abundance, with archaeal-associated KOs showing greater richness in metagenomes from high- temperature sites (Tables S6, S8). For example, K03170, reverse gyrase, was more rich at GBS 70 and GBS 85 than GBS 60 (Table S8). In contrast, K02470, DNA gyrase, was less rich in GBS 85 than GBS 60 (Table S8).
KOs associated with starch and glycogen biosynthesis and utilization [e.g., K01176, K05343, K01208, K16147, and K16149 (89.8% probability)] were significantly less abundant at GBS 85 than GBS 60. KOs associated with archaeal flagella increased in abundance as site temperature increased and were more abundant in GBS 85 than GBS 70 (e.g., K07332, K07333), while KOs associated with pili were significantly less abundant at GBS 85 than GBS 60 [e.g., K02658, K02660, and K02662 (89.8% probability)]. Genes associated with the biosynthesis of heat-stable lipids were more rich at GBS 85 than GBS 60: K13787, a type-I geranylgeranyl diphosphate synthase involved in archaeal and bacterial ether-linked lipid biosynthesis (Vandermoten et al., 2009) (abundance significance 60-85 88.8%); K17105, a geranylgeranylglycerol-phosphate geranylgeranyltransferase also involved in archaeal ether-linked lipid biosynthesis (abundance significance 60–85, 89.7%); and K17104, an archaeal phosphoglycerol geranylgeranyltransferase (abundance significance 60–85, 89.8%). KO K13789, a type-I geranylgeranyl diphosphate synthase involved in bacterial and plant ether-linked lipid biosynthesis (Vandermoten et al., 2009), was less rich in GBS 85 than in GBS 60 or GBS 70. KOs involved in lipopolysaccharide (LPS) transport and biosynthesis were less abundant in GBS 85 than in GBS 60 [e.g., K09695, K09694, K19428, and K04744 (89.7% probability)]. Multiple ABC transporters were significantly more abundant at GBS 60 than GBS 85 (Table S6). The archaeal cytoskeletal element crenactin, KO K18641 (Ettema et al., 2011), was more rich in GBS 70 and GBS 85 than in GBS 60 (abundance 60–85, 89.8%).
KOs associated with photosynthesis were significantly more abundant and richer in GBS 60 than in GBS 70 or GBS 85 (>40 KOs). In many cases (>40 KOs), GBS 70 and GBS 85 did not contain any copies of the genes. Many of these KOs were associated with photosystems I and II, chlorophyll biosynthesis, and photosynthetic reaction centers. KOs K01601 and K01602, the large and small chains of ribulose-bisphosphate carboxylase, were more abundant in GBS 60 and GBS 70 than in GBS 85. KO K14138, acetyl-CoA synthase, involved in the Wood-Ljungdahl pathway, was more abundant in GBS 60 and GBS 70 than in GBS 85. In contrast, many KOs associated with the dicarboxylate-hydroxybutyrate cycle and hydroxypropionate-hydroxybutyrate cycle were significantly less abundant at GBS 60 than GBS 85 (e.g., K14465, K14467, K15016, K15020, and K15038).
Few KOs involved in CCMPs were differentially abundant between metagenomes. For example, Figure 4A shows that enzymes carrying out key decarboxylation reactions in the CCMPs were similarly abundant and rich in GBS 60 and GBS 85 (Figures 4A–D). In cases where differences were found, functionally redundant KOs were identified. For example, K00163 is not necessary for pyruvate decarboxylation when K00161 and K00162 are present (Figures 4A,B). Entire CCMP KEGG maps for CCMP functions are shown in Figure S2 and showed similar trends. For example, K01810 is a glucose-6-phosphate isomerase assigned to EC 18.104.22.168 and was significantly less abundant in GBS 85 than in GBS 60 and GBS 70 (Table S6). Yet, KOs K06859, K13810, and K15916, which are also assigned to EC. 22.214.171.124, were found across all three metagenomes with non-significant differences (Table S6). K01622, a fructose 1,6-bisphosphate aldolase/phosphatase thought to be an ancestral gluconeogenic enzyme responsible for removing heat-labile triosephosphates in Archaea and some Bacteria (Say and Fuchs, 2010; Takami et al., 2012a; Bräsen et al., 2014) was more abundant in GBS 85 than in GBS 60 (Table S6), with richness being significantly higher in GBS 85 than in GBS 60 or GBS 70 (Table S8).
The 62 CAZymes determined to be significantly differentially abundant (Table S7) in any pairwise comparison were classified into CAZy families: carbohydrate binding modules (CBM) (9), glycosyltransferases (GT) (6), glycoside hydrolases (GH) (38), polysaccharide lyases (PL) (8), and dockerins (1). The majority decreased in abundance with increasing site temperature. CAZymes associated with peptidoglycan degradation were more abundant in GBS 60 (e.g., GH 73, GH102, and GH104) than GBS 85, and PL9_4 was significantly more abundant at GBS 85 than GBS 70. Family GT8, involved in LPS biosynthesis, was less abundant in GBS 85 than in GBS 60. Some families involved in cellulose degradation decreased in abundance with increasing site temperature and were less abundant in GBS 85 than GBS 60 (e.g., CBM8, CBM63, GH5_46, GH94, and dockerin). Families and sub-families involved in starch degradation and synthesis tended to be significantly more abundant in GBS 60 than GBS 85 (e.g., CBM20, CBM21, GH13_2, GH13_3, GH13_10, GH13_21, GH13_26, GH13_36, GH65, GH88, and GH128), including sub-family GH13_7 which contains characterized α-amylases from Archaea (Pyrococcus sp., Thermococcus sp.). Only six families and subfamilies were significantly more abundant at GBS 85 than GBS 60 (CBM29, GH52, GH96, GH122, PL14_3, and PL22_2). GH122 has a single characterized protein from Pyrococcus furiosus, described as an α-glucosidase with broad substrate specificity for α-glucan carbohydrates (Comfort et al., 2008).
An expanded visualization and discussion of differentially abundant KOs can be found in Figure S2.
13C-Based Metabolic Activity
Three different poisoned controls were tested in an effort to best differentiate biotic and abiotic decarboxylation activity: glutaraldehyde, mercury, and both glutaraldehyde and mercury. The controls varied in terms of suppressing 13CO2 production. Delta 13CO2 values were consistently higher in mercury-poisoned controls using uniformly labeled substrates than in glutaraldehyde-containing controls, with many mercury-poisoned controls producing equal or greater delta 13C values than the non-poisoned microcosms (Table S11). At the GBS 60 and GBS 95 sites, the only sites where glutaraldehyde was evaluated, glutaraldehyde-containing controls consistently produced lower delta 13C values than their substrate-paired samples and poisoned controls containing only mercury. Although a small amount of 13CO2 production was observed in glutaraldehyde-containing controls, they provided confidence that the majority of 13CO2 production from experimental microcosms could be attributed to biological activity.
Figure 4 summarizes key decarboxylation reactions that are diagnostic of glycolysis (Figure 4B), the TCA cycle (Figure 4C), and the oxidative PPP (Figure 4D), with integrated data from the metagenomes (Figure 4A, left of Figures 4B–D) and the isotopomer experiments (figure middle and right of Figures 4B–D). At all sites, added 13C-labeled compounds resulted in a near-linear return of 13CO2 with low variability between replicates for the duration of the experiment (Figures 4B–D middle; Table S11). Within a site, all microcosms contained a similar quantity of sediment (Table S12). All paired isotopomer ratios showed significant differences from complete catabolic use at all sites (p < 0.05; one-sample t-test) (Figures 4B–D right, Table S13), but in most cases the similar ratios of the different decarboxylation reactions regardless of temperature demonstrated that the relative activities of the pathways were unaffected by temperature. An exception was the significantly greater decarboxylation ratio of CU: C1 glucose at GBS 60, with a general trend of a decreasing ratio as site temperature increased, resulting in GBS 95 being significantly lower than GBS 60 and GBS 70 (one-way ANOVA, Figure 4D right, Table S13). This indicates an increasing trend for the C1 glucose atom to be sequestered in biomass or fermentation products at the lower end of the temperature gradient. Pyruvate isotopomer ratios were significantly larger at GBS 95 than at the other sites, indicating an increased sequestration of C2, 3 pyruvate atoms in biomass or fermentation products. Broad heterotrophic activity was observed at the GBS 85 site, as indicated by strong, linear returns from all isotopomer pairs, and uniformly labeled amino acids and TCA cycle intermediates. 13C-labeled amino acids and TCA cycle intermediates were only administered at GBS 85 (Figure 5; Table S11).
Figure 5. 13CO2 production from uniformly labeled 13C citrate, serine, cysteine, alanine, and succinate. Symbols represent the average APE (n = 3) for a given compound at the given time. R2-values ranged from 0.961 to 1.0. The standard error is not shown for simplicity. See Table S11 for slope and linearity data for all uniformly 13C-labeled compounds and replicates for each site, and raw δ13C of CO2 values for the time course. Compounds were tested at the 85°C site only.
Both amplicon-based and metagenome-derived 16S rRNA gene surveys demonstrated a decrease in diversity concomitant with an increase in temperature (Figure 2; Table 2), agreeing with previous diversity studies in GBS and globally (Cole et al., 2013; Sharp et al., 2014; Sunagawa et al., 2015; Power et al., 2018). As the temperature increased, a diverse, Bacteria-dominated community harboring phototrophs was replaced by an Archaea-dominated community consisting of a few thermophilic specialists (Figure 2), following well-known trends in geothermal systems (Hedlund et al., 2016). All sites had similar water chemistry and light exposure, suggesting temperature was the major driver in taxonomic shifts (Figure S1; Table S1).
All sites hosted a high proportion of candidate microbial taxa and unclassified sequences, conforming with previous studies at GBS (Cole et al., 2013; Rinke et al., 2013; Hedlund et al., 2015). In many cases, single OTUs from uncultivated groups accounted for large proportions of the community, highlighting the low diversity and phylogenetic novelty of this system. Additionally, several OTUs were only abundant in a single sample, suggesting distinct temperature niches and exemplifying the large taxonomic differences between samples. The proportion of “rare” 16S rRNA gene sequences (<1% abundance) was higher at low-temperature sites. High-temperature sites might not support a large rare-biosphere due to the energy demands of growth, maintenance, and survival at the physical extremes of life (Price and Sowers, 2004).
KEGG module completion ratios demonstrated the conservation of CCMPs along the temperature gradient for communities and individual MAGs. The fact that many MAGs contained complete CCMPs (Table 3; Table S4), suggests that C flux through these pathways can be facilitated by individual populations and are not reliant on metabolite swapping within the community.
All metagenomes contained a diversity of KOs and CAZymes. KOs associated with information processing were more abundant at higher temperatures, likely due to genome streamlining at higher temperatures (Sabath et al., 2013). Most of these were archaeal KOs, reflecting a shift to Archaea-dominated communities at high temperature. The large number of KOs and CAZymes that were found to be significantly different in richness only may reflect taxonomic differences, where the number of unique copies of a gene would be expected to decrease as taxonomic diversity decreases. Alternatively, the relative abundance of a gene may not decrease as taxonomic diversity decreases, if it provides for an essential function.
As expected, significant decreases in the abundance of genes associated with photosynthesis at high-temperature sites demonstrated the loss of a major energy source for autochthonous carbon fixation (Figure 3; Figure S2; Table S6). Differences in carbon fixation pathways were also observed with respect to temperature, with a lower number of modules for carbon fixation pathways predicted to be biologically feasible as site temperature increased (Figure 3; Table S5). Taken together, the reduction of taxonomic diversity concomitant with increasing site temperatures correlated with a loss in metabolic diversity with respect to primary production.
Genes associated with starch utilization, including an archaeal enzyme, decreased in abundance as temperature increased, suggesting this is not due to domain-level shifts in taxonomy. Genes associated with cellulosomes decreased in abundance with increasing temperature (e.g., dockerin). Cellulosomes are large extracellular protein complexes that may be too costly to produce and maintain in high-temperature environments. In fact, large membrane-bound or surface-associated enzyme complexes are not known to occur in hyperthermophiles (Hedlund et al., 2016). Although GBS 60 had the highest diversity of genes coding for ABC transporters, ABC transporters for various oligo- and monosaccharides, branched-chain amino acids and oligopeptides, and lipopolysaccharides and lipoproteins were present at relatively high abundance at all sites, suggesting heterotrophic flexibility and a conservation of function (Figure S2). Taken together, it appears that starch utilization is not prominent at high temperature, yet other allochthonous biomass from plants and insects may be important substrates for microbial communities in terrestrial geothermal springs (Schubotz et al., 2013).
Despite the temperature-driven decrease in taxonomic diversity, absence of photosynthesis, and shifts in genes associated with polymer degradation and information processing due to domain-level shifts in community composition, the functional potential of CCMPs remained consistent across sites (Figures 3, 4; Figure S2; Table S5), suggesting a conservation of CCMPs. The PPP was well-represented at the metagenomic level (Figures 3, 4D; Figure S2), despite the high relative abundance of Archaea at GBS 85, and representative MAGs (Table 3; Table S4) showing an incomplete PPP (Figure 3). Genes for the glyoxylate cycle were abundant at all sites, suggesting metabolic flexibility to adjust to the utilization of small carbon compounds (e.g., acetate; Figure S2). Abundance- and richness-populated KEGG maps for CCMPs (Figure S2) showed that the relative abundance of many KOs changed with temperature, yet the basic biochemistry was maintained (i.e., pathways were complete, demonstrated by populated EC numbers). This may be due to thermophilic Archaea having distinct enzymes for reactions in the CCMPs, with respect to Bacteria and Eukarya. This disparity may be due to pathway optimization for reduced protein cost and minimization of C and energy loss from thermolabile metabolites being favored over ATP yield, yet the basic biochemistry carried out by CCMPs remains similar (Flamholz et al., 2013; Bräsen et al., 2014). The conservation and diversity of CCMPs combined with shifts in genes associated with thermophily (e.g., DNA gyrases, lipid biosynthesis, cytoskeletal elements, flagella, reduced protein cost, enzymes facilitating reduced thermolabile metabolite pools) suggest that physiological adaptations are necessary to facilitate life at high temperatures, while the basic biochemistry of life is conserved.
In the present study, oxygen consumption rates likely would not deplete the microcosm headspace of O2 in the timeframe of the experiments (inferred from Murphy et al., 2013), yet, there almost certainly was an O2 gradient in the settled sediments of the microcosms, resulting in anaerobic habitats. Conceptually, anaerobic and aerobic respiration would be interpreted similarly with respect to substrate C atom partitioning and C flux through CCMPs. But fermentation, in addition to anabolism, is a likely contributor to the observed isotopomer mineralization ratios in the present study. In a community setting, fermentation products that are being produced faster than they are being metabolized by other community members would provide an overestimation of the contribution of anabolic reactions to experimentally determined isotopomer ratios. Despite this, C being transformed into fermentation products would not alter interpretations of organic carbon mineralization ratios, which could be used to infer C source/sink dynamics.
Substantial 13CO2 production in mercury-poisoned controls suggest that glutaraldehyde is more effective at suppressing heterotrophic activity in the GBS geothermal system (Table S11). GBS sediments consist of fine-grained particles and clay (Costa et al., 2009), potentially providing the opportunity for mercury to quickly sorb to sediment particles (Gabriel and Williamson, 2004). This may render the poison no longer bioavailable in concentrations necessary to inhibit microbial activity. This may have the effect of increasing relative rates of carbon oxidation due to increased maintenance energy demand and futile cycling from exposure to sub-lethal concentrations. Using position-specific 13C-labeled compounds, instead of uniformly 13C-labeled compounds, would allow a test of this hypothesis. Because poisoned controls containing both mercury and glutaraldehyde showed similar results to controls containing only glutaraldehyde, it is unlikely that mercury was reactive with the 13C-labeled compounds, suggesting high 13CO2 production rates in mercury-poisoned controls were due to biological activity. It should be noted that mercury is a better poison for inhibiting autotrophy in planktonic communities of the GBS system and little sediment is found in the water column (Boyd E. S., personal communication). We suggest that care should be taken when choosing poisons for evaluating biological activity and that combinations of poisons may be the most effective at inhibiting the biological activity of whole communities. Although glutaraldehyde was only used at the GBS 60 and GBS 95 sites, the 13CO2 analysis clearly indicated a drastic reduction in the decarboxylation of 13C-labeled compounds at both sites, suggesting 13CO2 from non-poisoned microcosms primarily reflected microbial activity.
Strong linear accumulation of 13CO2 was observed for all substrates at all sites (Figures 4; 5, Table S11), suggesting that communities were in a metabolic steady state (Zamboni et al., 2009). Additionally, this shows that the microbial communities at each site are adapted to the given substrates and are likely actively expressing these pathways, since carbon flux was immediate and constant for the duration of the experiment. These assumptions are supported by previous microrespirometry experiments conducted at GBS showing linear O2 consumption, without lag, upon the addition of organic substrates (Murphy et al., 2013).
A small amount of organic C was used in these experiments to minimize the effects of substrate additions on community metabolism (Dijkstra et al., 2011b). Yet, it has been shown that the communities in GBS sediments are electron-donor limited (Murphy et al., 2013). Thus, organic C additions may have stimulated the metabolic activity of community members and/or potentially altered the community composition over the duration of the experiment. We acknowledge that returns of 13C in the form of 13CO2 is likely due to the metabolism of a subset of microbes within these communities; we have shown that these pathways are intact within many, but not all, individuals of these communities (Figure 3). Attempts to use a metabolic model (Dijkstra et al., 2011b) integrating the different metabolites utilized in this study failed (data not shown). This suggests that pyruvate, glucose, and/or acetate are not being utilized by the same populations; therefore, we interpret our results in light of populations being defined by their use of a substrate (i.e., the glucose-utilizing population, pyruvate-utilizing population, or the acetate-utilizing population). We also acknowledge that several different metabolic strategies may be employed within a substrate-using population, and that we are observing the sum of these metabolic decarboxylation events. Nevertheless, the fast and steady metabolic activity and the care taken to maintain temperatures throughout the experiment provide confidence that experimental observations are ecologically relevant.
All paired isotopomer ratios showed significant differences from complete catabolic use at all sites (p < 0.05; one-sample t-test; Figures 4B–D), indicating biomass production and/or fermentation. Previous studies using position-specific 13C-labeled compounds are few and have focused on metabolism in aerobic soils (Dijkstra et al., 2011a,b,c, 2015; Van Groenigen et al., 2013; Hagerty et al., 2014). Except for the glucose isotopomer ratio at GBS 60, isotopomer ratios presented in this study are similar to those reported for soils. This indicates a similar retention of microbial substrate-derived carbon, relative to carbon uptake, in terrestrial geothermal settings and moderate-temperature soils. In the case of greater allocation of carbon toward anabolism and consequently a lowered relative organic C mineralization rate, a higher carbon use efficiency (CUE) (i.e., the ratio of organic C uptake allocated to biomass; mol C to biomass/ mol C uptake) is achieved. A similar CUE, despite the increased physiological stressor of high-temperature, could be a result of molecular adaptation to chronic energy stress associated with extreme environments (Valentine, 2007; Sabath et al., 2013), where high maintenance energy demands (Price and Sowers, 2004) have provided a selection pressure for efficient energy conservation. Regardless of high CUE (e.g., biomass production) or fermentation, the retention of organic carbon in the ecosystem may promote a more diverse community by providing a diversity of organic products to support diverse heterotrophs. Indeed, extreme ecosystems expected to increase maintenance energy demands (increasing organic carbon mineralization) and potentially lower CUE, such as low/high pH or high salinity, exhibit low diversity (Oren, 2001, 2011; Hou et al., 2013; Power et al., 2018).
Pyruvate isotopomer ratios suggest anabolism and/or fermentation are taking place due to ratios being significantly different from complete catabolic use (Figure 4B). Carbons from position 2 and 3 of pyruvate can be shunted from the TCA cycle for anabolic use or released in fermentation products (Figure 1). GBS 95 had a significantly higher pyruvate ratio indicating a greater proportion of pyruvate-derived carbon is sequestered in biomass and/or fermentation products than at lower temperature sites (Figure 4B). In the case of pyruvate fermentation to acetate, C2, 3 from pyruvate will be retained in acetate, while larger fermentation products generally result in the sequestration of all three C atoms in pyruvate (e.g., lactic acid). In the case of biomass production, a higher CUE would be expected at GBS 95 than at lower temperature sites, suggesting maintenance energy demand (i.e., catabolic carbon use) did not increase with temperature.
Acetate results are also significantly different from complete catabolic use and provide support for TCA cycle compounds being sequestered in biomass (Figure 4C). This is demonstrated by an acetate isotopomer ratio >1, suggesting the C2 position is being retained in biomass. In addition, acetate isotopomer ratios did not significantly differ across the sites. This indicates, again, that no change in the balance of anabolic and catabolic activities occurred, despite a nearly 35°C temperature range, suggesting a conservation of CUE and similar maintenance energy demands across temperatures.
In addition to isotopomer pairs, uniformly labeled citrate, serine, cysteine, alanine, and succinate, administered at the 85°C site only, showed at least partial catabolic utilization (Figure 5, Table S11). The use of these compounds provided additional support for a functioning TCA cycle, as the catabolism of these compounds is intimately linked to the TCA cycle and formation of pyruvate (Figure S2). These results, together with metagenomic data, provide support not just for an ecosystem-independent core functional gene assemblage for microbial communities (Sunagawa et al., 2015) but also for a conservation of function, particularly with respect to carbon flux through metabolic pathways and potential maintenance energy demands in high-temperature environments.
The 60°C site is well-below the upper-temperature limit for photosynthesis (Table 1) and several phototrophs are abundant at this site (Figure 2, green), yet incubations were carried out in the dark. This suggests this study captured the metabolic activity of a phototrophic mat during non-daylight conditions. GBS 60 showed a significant increase in the relative sequestration of C1 of glucose due to anabolism or fermentation, with respect to the other sites (Figure 4D). These results may be due to a reduced or inactive oxidative PPP (Dijkstra et al., 2011c, 2015). Without photosynthetic activity, anabolic reactions involving NADPH and PPP-derived biomass precursors (e.g., ribulose-5-phosphate), and the need to allocate NADPH to quenching reactive oxygen species is decreased, which would result in a low return from C1 of glucose. A reduced PPP is indicative of a lower demand for NADPH, which is used in many anabolic reactions directly linked to DNA replication and for scavenging of reactive oxygen species (Cox and Nelson, 2008). In addition, in the absence of light-driven primary production, a reduction in cellular replication would be expected. This is consistent with lower PPP activity due to a lower demand for the PPP-derived biomass precursor ribose-5-phosphate and NADPH for biosynthesis reactions. This interpretation fits with previous results for phototrophic mat communities in terrestrial geothermal systems that suggest that maximal rates of RNA, DNA, and protein production occur during morning hours with the lowest O2 concentrations found at night, reducing exposure to reactive oxygen species (Steunou et al., 2006; Kim et al., 2015).
In addition to a reduction in anabolism during dark periods, it has been shown that some members of photosynthetic mat communities inhabiting terrestrial geothermal springs switch to fermentative metabolisms during dark cycles (Anderson et al., 1987; Nold and Ward, 1996; Steunou et al., 2006; Kim et al., 2015). Furthermore, decreased PPP activity is a typical response to sub-oxic or anaerobic conditions (Godon et al., 1998; Gombert and Moreira, 2001; Le Goffe et al., 2002; Fredlund et al., 2004; Ralser et al., 2007; Williams et al., 2008; Chen et al., 2011; Dijkstra et al., 2011a). Mats in Mushroom Spring (YNP, U.S.A.) showed acetate and propionate accumulation (Kim et al., 2015), increased transcript levels for genes involved in fermentation, and decreased transcript levels for respiratory proteins at night (Steunou et al., 2006). Our observations of the catabolic use of acetate does not contradict the accumulation of acetate observed in Mushroom Spring (Kim et al., 2015), due to potential acetate production rates for some community members being greater than uptake rates of acetate-utilizing community members.
There is evidence that polyhydroxyalkanoic acid (PHA) and bacteriochlorophyll (BChl) biosynthesis is upregulated in the dark in photosynthetic mat communities, resulting in the majority of C from glucose being retained in biomass (Klatt et al., 2013). Under this scenario, with reduced PPP activity, C flux through biosynthetic pathways would show a strong 13CO2 signal from C3 and C4 of glucose (Dijkstra et al., 2015), from the formation of acetyl-CoA from glycolysis-derived pyruvate, providing a strong 13CO2 return from uniformly labeled glucose and a low return from C1 of glucose. The increased anabolism linked to PHA and BChl at night (Klatt et al., 2013) does not necessarily conflict with an overall reduction in anabolism (Kim et al., 2015). PHAs can be used as storage molecules that can be degraded for mixotrophic use in the light. It has been suggested that BChl biosynthesis in these communities is inhibited by O2 (Klatt et al., 2013). Taken together, the glucose-utilizing populations in GBS 60 may be replicating during the day (requiring total biomass replication and high PPP metabolite and NADPH demand) and building storage molecules and photosynthetic apparatus at night with minimal PPP activity (PHA formation and BChl-specific biomass precursor demand), resulting in greater retention of C1 of glucose, with respect to C1 oxidation in the PPP (Dijkstra et al., 2015). Future studies may achieve greater resolution of C flux through metabolic pathways from the use of a greater diversity of position-specific arrangements (Leighty and Antoniewicz, 2013; Dijkstra et al., 2015) and evaluating carbon flux during diel cycles.
In contrast to the glucose CU: C1 results at 60°C, higher temperature sites exhibited a decreased ratio, with GBS 95 being significantly lower than GBS 60 and GBS 70 (Figure 4D). This is assumed to be the result of 13CO2 being produced from the C1 position in the PPP (Figure 1) and/or reduced fermentation. This may indicate the production of nucleotides and NADPH, both indicative of anabolism directly related to cellular replication, in contrast to the potential sequestration of C1 from anabolism of PAHs and BChl suggested at the 60°C site. On the other hand, an increase in glycolytic C flux through the Entner–Doudoroff (ED) pathway would result in an increase in C1 oxidation during acetate formation. This is due to glucose conversion to pyruvate resulting in one molecule of pyruvate retaining the C1 of glucose in the C1 position of pyruvate, and subsequent pyruvate oxidation in the synthesis of acetyl-coA (i.e., 13CO2 production from C1 of pyruvate). In fact, Archaea are known to have incomplete PPP pathways, which is shown in representative MAGs, while many utilize modified versions of the ED pathway (Bräsen et al., 2014). Populated KEGG maps and MCRs produced in this study contrast with respect to the completeness of the ED pathway. While both confirm the typical ED pathway is only complete in the GBS 70 metagenome, MCR analysis suggests that the non-phosphorylative ED pathway is biologically feasible in all metagenomes and that the semi-phosphorylative ED pathway is not (Table S5). In contrast, populated KEGG maps (Kanehisa and Goto, 2000; Kanehisa et al., 2017, 2012) suggest that the non-phosphorylative ED pathway is not possible in GBS 85, due to the absence of EC 126.96.36.199 (K11395), which has a low relative abundance in GBS 70, while the semi-phosphorylative ED pathway appears possible (Figure S2, Table S6). The ED pathways requires seven-fold less enzymatic protein than glycolysis, the use of which may be a result of environmental selection pressures for reduced protein load at the cost of ATP gain (1 vs. 2 ATP, with respect to glycolysis) as a consequence of increased maintenance energy costs associated with protein repair and turnover at high temperatures (Price and Sowers, 2004; Flamholz et al., 2013).
Future studies should strive to obtain net CO2 production rates and 13C label consumption rates to allow extrapolations to absolute rates of C flux, community respiration rates, and calculations of carbon mass balance. Additional studies would benefit from pairing with additional “omics” techniques (i.e., transcriptomics and proteomics) and organism-level stable isotope techniques [e.g., quantitative stable isotope probing Hungate et al., 2015; Morrissey et al., 2016 FISH-NanoSIMS Dekas and Orphan, 2010; Musat et al., 2012, 2016]. As more microbial communities are probed by isotopomer ratio studies, cross-community comparisons can be made to better inform relationships between diversity and microbially-mediated C biogeochemistry and maintenance energy demands (CUE), and progress can be made toward a quantitative understanding of microbial ecosystems.
In summary, this study demonstrated that, despite domain-level shifts in community composition and a decrease in diversity with increasing temperature, CCMP genes and biochemical efficiency were conserved across a broad temperature gradient in GBS. Evidence for a stable or decreased maintenance energy demand with increasing habitat temperature was found: pyruvate and acetate utilization across all sites, and glucose from 70 to 95°C, suggest a high CUE and retention of organic carbon at all temperatures. In contrast, significant deviations in C flux through CCMPs with respect to glucose were observed in the presence of phototrophs.
The datasets generated for this study can be found in NCBI Sequence Read Archive; Integrated Microbial Genomes and Microbiomes platform, SRP059341, SRX1055762, SRX1055763, SRX1055764, and SRX1055765; Ga0197142, Ga0197143, and Ga0197144.
ST, BH, KT, and PD contributed to experimental design. ST and KT completed field work. ST, KT, and JD completed laboratory work. ST, JD, CS, SM, and DL contributed to bioinformatics work. EE-F oversaw metagenomic sequencing, assembly, and MAG generation. ST analyzed the data and wrote the manuscript with input from all authors.
The work conducted by the U.S. Department of Energy Joint Genome Institute, a DOE Office of Science User Facility, is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. This work was supported by the National Science Foundation (DEB 1557042 and OISE 0968421). BH acknowledges generous support from Greg Fullmer through the UNLV Foundation.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank David and Sandy Jamieson for gracious support and access to GBS. We thank the lab of Dr. Hailiang Dong at Miami University for analyzing samples for major anions, cations, and trace elements. We also thank Dr. Philippos Tsourkas at UNLV for valuable assistance with statistical approach.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2019.01427/full#supplementary-material
Arai, W., Taniguchi, T., Goto, S., Moriya, Y., Uehara, H., Takemoto, K., et al. (2018). MAPLE 2. 3. 0: an improved system for evaluating the functionomes of genomes and metagenomes. Biosci. Biotechnol. Biochem. 82, 1515–1517. doi: 10.1080/09168451.2018.1476122
Baker, B. J., Lazar, C. S., Teske, A. P., and Dick, G. J. (2015). Genomic resolution of linkages in carbon, nitrogen, and sulfur cycling among widespread estuary sediment bacteria. Microbiome 3, 1–12. doi: 10.1186/s40168-015-0077-6
Beam, J. P., Jay, Z. J., Schmid, M. C., Rusch, D. B., Romine, M. F., Jennings Rde, M., et al. (2016). Ecophysiology of an uncultivated lineage of Aigarchaeota from an oxic, hot spring filamentous ‘streamer' community. ISME J. 10, 210–224. doi: 10.1038/ismej.2015.83
Becraft, E. D., Dodsworth, J. A., Murugapiran, S. K., Ohlsson, I., Briggs, B. R., Kanbar, J., et al. (2016). Single-cell-genomics-facilitated read binning of candidate phylum EM19 genomes from geothermal spring metagenomes. Appl. Environ. Microbiol. 84, 992–1003. doi: 10.1128/AEM.03140-15
Becraft, E. D., Dodsworth, J. A., Murugapiran, S. K., Thomas, S. C., Ohlsson, J. I., Stepanauskas, R., et al. (2017). Genomic comparison of two family-level groups of the uncultivated NAG1 archaeal lineage from chemically and geographically disparate hot springs. Front. Microbiol. 8:2082. doi: 10.3389/fmicb.2017.02082
Bowers, R. M., Kyrpides, N. C., Stepanauskas, R., Harmon-smith, M., Doud, D., Reddy, T. B. K., et al. (2017). perspective Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat. Biotechnol. 35, 725–731. doi: 10.1038/nbt.3893
Boyd, E. S., Fecteau, K. M., Havig, J. R., Shock, E. L., and Peters, J. W. (2012). Modeling the habitat range of phototrophs in Yellowstone National Park: toward the development of a comprehensive fitness landscape. Front. Microbiol. 3:221. doi: 10.3389/fmicb.2012.00221
Bräsen, C., Esser, D., Rauch, B., and Siebers, B. (2014). Carbohydrate metabolism in archaea: current insights into unusual enzymes and pathways and their regulation. Microbiol. Mol. Biol. Rev. 78, 89–175. doi: 10.1128/MMBR.00041-13
Brock, T. D. (1967). Life at High Temperatures: Evolutionary, ecological, and biochemical significance of organisms living in hot springs is discussed. Science 158, 1012–1019. doi: 10.1126/science.158.3804.1012
Chen, I. A., Chu, K., Palaniappan, K., Pillay, M., Ratner, A., Huang, J., et al. (2019). IMG / M v. 5. 0: an integrated data management and comparative analysis system for microbial genomes and microbiomes. Nucleic Acids Res. 47, 666–677. doi: 10.1093/nar/gky901
Chen, X., Alonso, A. P., Allen, D. K., Reed, J. L., and Shachar-Hill, Y. (2011). Synergy between 13C-metabolic flux analysis and flux balance analysis for understanding metabolic adaption to anaerobiosis in E. coli. Metab. Eng. 13, 38–48. doi: 10.1016/j.ymben.2010.11.004
Cole, J. K., Peacock, J. P., Dodsworth, J. A., Williams, A. J., Thompson, D. B., Dong, H., et al. (2013). Sediment microbial communities in Great Boiling Spring are controlled by temperature and distinct from water communities. ISME J. 7, 718–729. doi: 10.1038/ismej.2012.157
Comfort, D. A., Chou, C., Conners, S. B., Vanfossen, A. L., and Kelly, R. M. (2008). Functional-genomics-based identification and characterization of open reading frames encoding α-glucoside-processing enzymes in the hyperthermophilic Archaeon Pyrococcus furiosus. Appl. Environ. Microbiol. 74, 1281–1283. doi: 10.1128/AEM.01920-07
Costa, K. C., Navarro, J. B., Shock, E. L., Zhang, C. L., Soukup, D., and Hedlund, B. P. (2009). Microbiology and geochemistry of great boiling and mud hot springs in the United States Great Basin. Extremophiles 13, 447–459. doi: 10.1007/s00792-009-0230-x
Cox, M. M., and Nelson, D. L. (2008). “Chapter 19: oxidative phosphorylation,” in Lehninger Principles of Biochemistry, ed K. Ahr (New York, NY: W.H. Freeman and Company), 707–762. Available online at: https://www.amazon.com/Lehninger-Principles-Biochemistry-David-Nelson/dp/071677108X/ref=pd_lpo_sbs_14_t_0?_encoding=UTF8&psc=1&refRID=BDSRXRGN4ZQE4P3EG2C0
De La Torre, J. R., Walker, C. B., Ingalls, A. E., Könneke, M., and Stahl, D. A. (2008). Cultivation of a thermophilic ammonia oxidizing archaeon synthesizing crenarchaeol. Environ. Microbiol. 10, 810–818. doi: 10.1111/j.1462-2920.2007.01506.x
Dekas, A. E., and Orphan, V. J. (2010). Identification of Diazotrophic Microorganisms in Marine Sediment via Fluorescence in situ Hybridization Coupled to Nanoscale Secondary Ion Mass Spectrometry (FISH-NanoSIMS), 1st Edn. San Diego, CA: Elsevier Inc. doi: 10.1016/B978-0-12-381294-0.00012-2
Dijkstra, P., Blankinship, J. C., Selmants, P. C., Hart, S. C., Koch, G. W., Schwartz, E., et al. (2011a). Probing carbon flux patterns through soil microbial metabolic networks using parallel position-specific tracer labeling. Soil Biol. Biochem. 43, 126–132. doi: 10.1016/j.soilbio.2010.09.022
Dijkstra, P., Dalder, J. J., Selmants, P. C., Hart, S. C., Koch, G. W., Schwartz, E., et al. (2011b). Modeling soil metabolic processes using isotopologue pairs of position-specific 13C-labeled glucose and pyruvate. Soil Biol. Biochem. 43, 1848–1857. doi: 10.1016/j.soilbio.2011.05.001
Dijkstra, P., Salpas, E., Fairbanks, D., Miller, E. B., Hagerty, S. B., van Groenigen, K. J., et al. (2015). High carbon use efficiency in soil microbial communities is related to balanced growth, not storage compound synthesis. Soil Biol. Biochem. 89, 35–43. doi: 10.1016/j.soilbio.2015.06.021
Dijkstra, P., Thomas, S. C., Heinrich, P. L., Koch, G. W., Schwartz, E., and Hungate, B. A. (2011c). Effect of temperature on metabolic activity of intact microbial communities: evidence for altered metabolic pathway activity but not for increased maintenance respiration and reduced carbon use efficiency. Soil Biol. Biochem. 43, 2023–2031. doi: 10.1016/j.soilbio.2011.05.018
Dodsworth, J. A., Blainey, P. C., Murugapiran, S. K., Swingley, W. D., Ross, C. A., Tringe, S. G., et al. (2013). Single-cell and metagenomic analyses indicate a fermentative and saccharolytic lifestyle for members of the OP9 lineage. Nat. Commun. 4:1854. doi: 10.1038/ncomms2884
Dodsworth, J. A., Gevorkian, J., Despujos, F., Cole, J. K., Murugapiran, S. K., Ming, H., et al. (2014). Thermoflexus hugenholtzii gen. nov., sp. nov., a thermophilic, microaerophilic, filamentous bacterium representing a novel class in the Chloroflexi, Thermoflexia classis nov., and description of Thermoflexaceae fam. nov. and Thermoflexales ord. nov. Int. J. Syst. Evol. Microbiol. 64, 2119–2127. doi: 10.1099/ijs.0.055855-0
Dodsworth, J. A., Hungate, B. A., and Hedlund, B. P. (2011). Ammonia oxidation, denitrification and dissimilatory nitrate reduction to ammonium in two US Great Basin hot springs with abundant ammonia-oxidizing archaea. Environ. Microbiol. 13, 2371–2386. doi: 10.1111/j.1462-2920.2011.02508.x
Dodsworth, J. A., Mcdonald, A. I., and Hedlund, B. P. (2012). Calculation of total free energy yield as an alternative approach for predicting the importance of potential chemolithotrophic reactions in geothermal springs. FEMS Microbiol. Ecol. 81, 446–454. doi: 10.1111/j.1574-6941.2012.01369.x
Edwards, T. A., Calica, N. A., Huang, D. A., Manoharan, N., Hou, W., Huang, L., et al. (2013). Cultivation and characterization of thermophilic Nitrospira species from geothermal springs in the US Great Basin, China, and Armenia. FEMS Microbiol. Ecol. 85, 283–292. doi: 10.1111/1574-6941.12117
Eloe-Fadrosh, E. A., Paez-Espino, D., Jarett, J., Dunfield, P. F., Hedlund, B. P., Dekas, A. E., et al. (2016). Global metagenomic survey reveals a new bacterial candidate phylum in geothermal springs. Nat. Commun. 7:10476. doi: 10.1038/ncomms10476
Flamholz, A., Noor, E., Bar-Even, A., Liebermeister, W., and Milo, R. (2013). Glycolytic strategy as a tradeoff between energy yield and protein cost. Proc. Natl. Acad. Sci. U.S.A. 110, 10039–10044. doi: 10.1073/pnas.1215283110
Fredlund, E., Blank, L. M., Schnürer, J., Schnu, U., and Passoth, V. (2004). Oxygen- and glucose-dependent regulation of central carbon metabolism in Pichia anomala. Appl. Environ. Microbiol. 70, 5905–5911. doi: 10.1128/AEM.70.10.5905-5911.2004
Gabriel, M. C., and Williamson, D. G. (2004). Principal biogeochemical factors affecting the speciation and transport of mercury through the terrestrial environment. Environ Geochem Health. 26, 421–434. doi: 10.1007/s10653-004-1308-0
Gombert, A. K., and Moreira, M. (2001). Network identification and flux quantification in the central metabolism of Saccharomyces cerevisiae under different conditions of glucose repression network identification and flux quantification in the central metabolism of Saccharomyces cerev. J. Bacteriol. 183, 1441–1451. doi: 10.1128/JB.183.4.1441-1451.2001
Hagerty, S. B., Van Groenigen, K. J., Allison, S. D., Hungate, B. A., Schwartz, E., Koch, G. W., et al. (2014). Accelerated microbial turnover but constant growth efficiency with warming in soil. Nat. Clim. Chang. 4, 903–906. doi: 10.1038/nclimate2361
Hedlund, B. P., Murugapiran, S. K., Alba, T. W., Levy, A., Dodsworth, J. A., Goertz, G. B., et al. (2015). Uncultivated thermophiles: current status and spotlight on “Aigarchaeota.” Curr. Opin. Microbiol. 25, 136–145. doi: 10.1016/j.mib.2015.06.008
Hedlund, B. P., Thomas, S. C., Dodsworth, J. A., and Zhang, C. L. (2016). “Life in high-temperature environments,” in Manual of Environmental Microbiology, 4th Edn, eds M. V. Yates, C. H. Nakatsu, R. V. Miller, and S. D. Pillai (Washington, DC: ASM Press), 863–887.
Hou, W., Wang, S., Dong, H., Jiang, H., Briggs, B. R., Peacock, J. P., et al. (2013). A comprehensive census of microbial diversity in hot springs of Tengchong, Yunnan Province China using 16S rRNA gene pyrosequencing. PLoS ONE 8:e53350. doi: 10.1371/journal.pone.0053350
Hua, Z. S., Qu, Y. N., Zhu, Q., Zhou, E. M., Qi, Y. L., Yin, Y. R., et al. (2018). Genomic inference of the metabolism and evolution of the archaeal phylum Aigarchaeota. Nat Commun. 9:2832. doi: 10.1038/s41467-018-05284-4
Huang, Q., Jiang, H., Briggs, B. R., Wang, S., Hou, W., Li, G., et al. (2013). Archaeal and bacterial diversity in acidic to circumneutral hot springs in the Philippines. FEMS Microbiol. Ecol. 85, 452–464. doi: 10.1111/1574-6941.12134
Hungate, B., Mau, R. L., Schwartz, E., Caporaso, J. G., Dijkstra, P., van Gestel, N., et al. (2015). Quantitative microbial ecology through stable isotope probing. Soil Biol. Biochem. 81, 7570–7581. doi: 10.1128/AEM.02280-15
Huntemann, M., Ivanova, N. N., Mavromatis, K., Tripp, H. J., Paez-Espino, D., Tennessen, K., et al. (2016). The standard operating procedure of the DOE-JGI Metagenome Annotation Pipeline (MAP v.4). Stand. Genomic Sci. 11:17. doi: 10.1186/s40793-016-0138-x
Hyatt, D., Chen, G., Locascio, P. F., Land, M. L., Larimer, F. W., and Hauser, L. J. (2010). Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinform. 11:119. doi: 10.1186/1471-2105-11-119
Jain, C., Rodriguez, R. L. M., Phillippy, A. M., Konstantinidis, K. T., and Aluru, S. (2018). High-throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9:5114. doi: 10.1038/s41467-018-07641-9
Jennings, R. M., Moran, J. J., Jay, Z. J., Beam, J. P., Whitmore, L. M., Kozubal, M. A., et al. (2017). Integration of metagenomic and stable carbon isotope evidence reveals the extent and mechanisms of carbon dioxide fixation in high-temperature microbial communities. Front. Microbiol. 8:88. doi: 10.3389/fmicb.2017.00088
Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., and Morishima, K. (2017). KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45, D353–D361. doi: 10.1093/nar/gkw1092
Kanehisa, M., Goto, S., Sato, Y., Furumichi, M., and Tanabe, M. (2012). KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Res. 40, D109–D114. doi: 10.1093/nar/gkr988
Kang, D. D., Froula, J., Egan, R., and Wang, Z. (2015). MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities. PeerJ 3:e1165. doi: 10.7717/peerj.1165
Kim, Y. M., Nowack, S., Olsen, M., Becraft, E. D., Wood, J. M., Thiel, V., et al. (2015). Diel metabolomics analysis of a hot spring chlorophototrophic microbial mat leads to new hypotheses of community member metabolisms. Front. Microbiol. 6:209. doi: 10.3389/fmicb.2015.00209
Klatt, C. G., Liu, Z., Ludwig, M., Kühl, M., Jensen, S. I., Bryant, D. A., et al. (2013). Temporal metatranscriptomic patterning in phototrophic Chloroflexi inhabiting a microbial mat in a geothermal spring. ISME J. 7, 1775–1789. doi: 10.1038/ismej.2013.52
Kozich, J. J., Westcott, S. L., Baxter, N. T., Highlander, S. K., and Schloss, P. D. (2013). Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the miseq illumina sequencing platform. Appl. Environ. Microbiol. 79, 5112–5120. doi: 10.1128/AEM.01043-13
Le Goffe, C., Vallette, G., Charrier, L., Candelon, T., Bou-Hanna, C., Bouhours, J. F., et al. (2002). Metabolic control of resistance of human epithelial cells to H2O2 and NO stresses. Biochem. J. 364, 349–359. doi: 10.1042/bj20011856
Lebedeva, E. V., Alawi, M., Fiencke, C., Namsaraev, B., Bock, E., and Spieck, E. (2005). Moderately thermophilic nitrifying bacteria from a hot spring of the Baikal rift zone. FEMS Microbiol. Ecol. 54, 297–306. doi: 10.1016/j.femsec.2005.04.010
Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M., and Henrissat, B. (2014). The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res. 42, 490–495. doi: 10.1093/nar/gkt1178
Matsen, F. A., Kodner, R. B., and Armbrust, E. V. (2010). pplacer: linear time maximum-likelihood and Bayesian phylogenetic placement of sequences onto a fixed reference tree. BMC Bioinform. 11:538. doi: 10.1186/1471-2105-11-538
Morgulis, A., Coulouris, G., Raytselis, Y., Madden, T. L., Agarwala, R., and Schäffer, A. A. (2008). Database indexing for production MegaBLAST searches. Bioinformatics 24, 1757–1764. doi: 10.1093/bioinformatics/btn322
Morrissey, E. M., Mau, R. L., Schwartz, E., Caporaso, J. G., DIjkstra, P., Van Gestel, N., et al. (2016). Phylogenetic organization of bacterial activity. ISME J. 10, 2336–2340. doi: 10.1038/ismej.2016.28
Murphy, C. N., Dodsworth, J. A., Babbitt, A. B., and Hedlund, B. P. (2013). Community microrespirometry and molecular analyses reveal a diverse energy economy in great boiling spring and sandy's spring west in the U.S. great basin. Appl. Environ. Microbiol. 79, 3306–3310. doi: 10.1128/AEM.00139-13
Musat, N., Foster, R., Vagner, T., Adam, B., and Kuypers, M. M. (2012). Detecting metabolic activities in single cells, with emphasis on nanoSIMS. FEMS Microbiol. Rev. 36, 486–511. doi: 10.1111/j.1574-6976.2011.00303.x
Nunoura, T., Takaki, Y., Kakuta, J., Nishi, S., Sugahara, J., Kazama, H., et al. (2011). Insights into the evolution of Archaea and eukaryotic protein modifier systems revealed by the genome of a novel archaeal group. Nucleic Acids Res. 39, 3204–3223. doi: 10.1093/nar/gkq1228
Oren, A. (2001). The bioenergetic basis for the decrease in metabolic diversity at increasing salt concentrations: implications for the functioning of salt lake ecosystems. Hydrobiologia 466, 61–72. doi: 10.1023/A:1014557116838
Parks, D. H., Chuvochina, M., Waite, D. W., Rinke, C., Skarshewski, A., Chaumeil, P., et al. (2018). A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat. Biotechnol. 36, 996–1004. doi: 10.1038/nbt.4229
Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P., and Tyson, G. W. (2015). CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055. doi: 10.1101/gr.186072.114
Power, J. F., Carere, C. R., Lee, C. K., Wakerley, G. L. J., Evans, D. W., Button, M., et al. (2018). Microbial biogeography of 925 geothermal springs in New Zealand. Nat. Commun. 9:2876. doi: 10.1038/s41467-018-05020-y
Price, P. B., and Sowers, T. (2004). Temperature dependence of metabolic rates for microbial growth, maintenance, and survival. Proc. Natl. Acad. Sci. U.S.A. 101, 4631–4636. doi: 10.1073/pnas.0400522101
Pruesse, E., Peplies, J., and Glöckner, F. O. (2012). SINA: accurate high-throughput multiple sequence alignment of ribosomal RNA genes. Bioinformatics 28, 1823–1829. doi: 10.1093/bioinformatics/bts252
Quast, C., Pruesse, E., Yilmaz, P., Gerken, J., Schweer, T., Yarza, P., et al. (2013). The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, 590–596. doi: 10.1093/nar/gks,1219
R Core Team (2017). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Available online at: https://www.R-project.org/
Ralser, M., Wamelink, M. M., Kowald, A., Gerisch, B., Heeren, G., Struys, E. A., et al. (2007). Dynamic rerouting of the carbohydrate flux is key to counteracting oxidative stress. J. Biol. 6:10. doi: 10.1186/jbiol61
Rinke, C., Schwientek, P., Sczyrba, A., Ivanova, N. N., Anderson, I. J., Cheng, J. F., et al. (2013). Insights into the phylogeny and coding potential of microbial dark matter. Nature 499, 431–437. doi: 10.1038/nature12352
Sabath, N., Ferrada, E., Barve, A., and Wagner, A. (2013). Growth temperature and genome size in bacteria are negatively correlated, suggesting genomic streamlining during thermal adaptation. Genome Biol. Evol. 5, 966–977. doi: 10.1093/gbe/evt050
Schubotz, F., Hays, L. E., Meyer-Dombard, D. R., Gillespie, A., Shock, E. L., and Summons, R. E. (2015). Stable isotope labeling confirms mixotrophic nature of streamer biofilm communities at alkaline hot springs. Front. Microbiol. 6:42. doi: 10.3389/fmicb.2015.00042
Schubotz, F., Meyer-Dombard, D. R., Bradley, A. S., Fredricks, H. F., Hinrichs, K. U., Shock, E. L., et al. (2013). Spatial and temporal variability of biomarkers and microbial diversity reveal metabolic and community flexibility in streamer biofilm communities in the Lower Geyser Basin, Yellowstone National Park. Geobiology 11, 549–569. doi: 10.1111/gbi.12051
Sharp, C. E., Brady, A. L., Sharp, G. H., Grasby, S. E., Stott, M. B., and Dunfield, P. F. (2014). Humboldt's spa: microbial diversity is controlled by temperature in geothermal environments. ISME J. 8, 1166–1174. doi: 10.1038/ismej.2013.237
Sorokin, D. Y., Lücker, S., Vejmelkova, D., Kostrikina, N. A., Kleerebezem, R., Rijpstra, W. I., et al. (2012). Nitrification expanded: discovery, physiology and genomics of a nitrite-oxidizing bacterium from the phylum Chloroflexi. ISME J. 6, 2245–2256. doi: 10.1038/ismej.2012.70
Steunou, A., Bhaya, D., Bateson, M. M., Melendrez, M. C., Ward, D. M., Brecht, E., et al. (2006). In situ analysis of nitrogen fixation and metabolic switching in unicellular thermophilic cyanobacteria inhabiting hot spring microbial mats. Proc. Natl. Acad. Sci. U.S.A. 130, 2398–2403. doi: 10.1073/pnas.0507513103
Sunagawa, S., Coelho, L. P., Chaffron, S., Kultima, J. R., Labadie, K., Salazar, G., et al. (2015). Ocean plankton. Structure and function of the global ocean microbiome. Science 348:1261359. doi: 10.1126/science.1261359
Swingley, W. D., Meyer-Dombard, D. R., Shock, E. L., Alsop, E. B., Falenski, H. D., Havig, J. R., et al. (2012). Coordinating environmental genomics and geochemistry reveals metabolic transitions in a hot spring ecosystem. PLoS ONE 7:e38108. doi: 10.1371/journal.pone.0038108
Takami, H., Noguchi, H., Takaki, Y., Uchiyama, I., Toyoda, A., Nishi, S., et al. (2012a). A deeply branching thermophilic bacterium with an ancient Acetyl-CoA pathway dominates a subsurface ecosystem. PLoS ONE 7:e30559. doi: 10.1371/journal.pone.0030559
Takami, H., Taniguchi, T., Arai, W., Takemoto, K., Moriya, Y., and Goto, S. (2016). An automated system for evaluation of the potential functionome: MAPLE version 2. 1. 0. DNA Res. 23, 467–475. doi: 10.1093/dnares/dsw030
Takami, H., Taniguchi, T., Moriya, Y., Kuwahara, T., Kanehisa, M., and Goto, S. (2012b). Evaluation method for the potential functionome harbored in the genome and metagenome. BMC Genomics 13:699. doi: 10.1186/1471-2164-13-699
Tarazona, S., Furió-Tarí, P., Turrà, D., Pietro, A. D., Nueda, M. J., Ferrer, A., et al. (2015). Data quality aware analysis of differential expression in RNA-seq with NOISeq R/Bioc package. Nucleic Acids Res. 43:e140. doi: 10.1093/nar/gkv711
Tijhuis, L., Van Loosdrecht, M. C., and Heijnen, J. J. (1993). A thermodynamically based correlation for maintenance gibbs energy requirements in aerobic and anaerobic chemotrophic growth. Biotechnol Bioeng. 42, 509–519. doi: 10.1002/bit.260420415
Van Groenigen, K. J., Forristal, D., Jones, M., Smyth, N., Schwartz, E., Hungate, B., et al. (2013). Using metabolic tracer techniques to assess the impact of tillage and straw management on microbial carbon use efficiency in soil. Soil Biol. Biochem. 66, 139–145. doi: 10.1016/j.soilbio.2013.07.002
Vandermoten, S., Haubruge, E., and Cusson, M. (2009). New insights into short-chain prenyltransferases: structural features, evolutionary history and potential for selective inhibition. Cell. Mol. Life Sci. 66, 3685–3695. doi: 10.1007/s00018-009-0100-9
Williams, T. C., Miguet, L., Masakapalli, S. K., Kruger, N. J., Sweetlove, L. J., and Ratcliffe, R. G. (2008). Metabolic network fluxes in heterotrophic Arabidopsis cells: stability of the flux distribution under different oxygenation conditions. Plant Physiol. 148, 704–718. doi: 10.1104/pp.108.125195
Yarza, P., Yilmaz, P., Pruesse, E., Glöckner, F. O., Ludwig, W., Schleifer, K. H., et al. (2014). Uniting the classification of cultured and uncultured bacteria and archaea using 16S rRNA gene sequences. Nat. Rev. Microbiol. 12, 635–645. doi: 10.1038/nrmicro3330
Yilmaz, P., Parfrey, L. W., Yarza, P., Gerken, J., Pruesse, E., Quast, C., et al. (2014). The SILVA and “all-species Living Tree Project (LTP)” taxonomic frameworks. Nucleic Acids Res. 42, 1–6. doi: 10.1093/nar/gkt1209
Yoon, S. H., Ha, S. M., Kwon, S., Lim, J., Kim, Y., Seo, H., et al. (2017). Introducing EzBioCloud: a taxonomically united database of 16S rRNA gene sequences and whole-genome assemblies. Int. J. Syst. Evol. Microbiol. 67, 1613–1617. doi: 10.1099/ijsem.0.001755
Keywords: 13C stable isotope, isotopomers, heterotrophy, catabolic, anabolic, thermophile, hyperthermophile, diversity
Citation: Thomas SC, Tamadonfar KO, Seymour CO, Lai D, Dodsworth JA, Murugapiran SK, Eloe-Fadrosh EA, Dijkstra P and Hedlund BP (2019) Position-Specific Metabolic Probing and Metagenomics of Microbial Communities Reveal Conserved Central Carbon Metabolic Network Activities at High Temperatures. Front. Microbiol. 10:1427. doi: 10.3389/fmicb.2019.01427
Received: 23 February 2019; Accepted: 05 June 2019;
Published: 05 July 2019.
Edited by:Don A. Cowan, University of Pretoria, South Africa
Reviewed by:Huiluo Cao, the University of Hong Kong, Hong Kong
Daniel Lundin, Linnaeus University, Sweden
Copyright © 2019 Thomas, Tamadonfar, Seymour, Lai, Dodsworth, Murugapiran, Eloe-Fadrosh, Dijkstra and Hedlund. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Brian P. Hedlund, email@example.com
†Present Address: Kevin O. Tamadonfar, Department of Molecular Microbiology, Washington University, St. Louis, MO, United States Senthil K. Murugapiran, Department of Plant and Microbial Biology, University of Minnesota, St. Paul, MN, United States