Quantifying genome-specific carbon fixation in a 750-meter deep subsurface hydrothermal microbial community

Abstract

Dissolved inorganic carbon has been hypothesized to stimulate microbial chemoautotrophic activity as a biological sink in the carbon cycle of deep subsurface environments. Here, we tested this hypothesis using quantitative DNA stable isotope probing of metagenome-assembled genomes (MAGs) at multiple ¹³C-labeled bicarbonate concentrations in hydrothermal fluids from a 750-m deep subsurface aquifer in the Biga Peninsula (Turkey). The diversity of microbial populations assimilating ¹³C-labeled bicarbonate was significantly different at higher bicarbonate concentrations, and could be linked to four separate carbon-fixation pathways encoded within ¹³C-labeled MAGs. Microbial populations encoding the Calvin–Benson–Bassham cycle had the highest contribution to carbon fixation across all bicarbonate concentrations tested, spanning 1–10 mM. However, out of all the active carbon-fixation pathways detected, MAGs affiliated with the phylum Aquificae encoding the reverse tricarboxylic acid (rTCA) pathway were the only microbial populations that exhibited an increased ¹³C-bicarbonate assimilation under increasing bicarbonate concentrations. Our study provides the first experimental data supporting predictions that increased bicarbonate concentrations may promote chemoautotrophy via the rTCA cycle and its biological sink for deep subsurface inorganic carbon.

deep biosphere, DNA stable isotope probing, microbial carbon fixation, qSIP

Introduction

Deep subsurface hydrothermal environments house enigmatic microbial communities in the absence of sunlight, which are understood to be largely sustained by chemolithoautotrophic metabolism (Stevens 1997, Magnabosco et al. 2016, Power et al. 2018). Since photosynthetically derived organic carbon is often a limiting carbon source for microbial communities in the deep subsurface (Pedersen 2000), in situ production of dissolved organic carbon (DOC) by chemolithoautotrophs helps to sustain subsurface microbial life subsisting under extreme energy limitation (Lau et al. 2016, Magnabosco et al. 2016, Momper et al. 2017, Probst et al. 2020). Subsurface environments are the largest ecosystems for prokaryotes with a biomass estimated to be 2–6 × 10²⁹ cells, accounting for 23–31 petagrams of carbon (PgC) (Magnabosco et al. 2018), and act as a subsurface biological sink for deep sources of inorganic carbon and CO₂ (Magnabosco et al. 2018, Merino et al. 2019, Fullerton et al. 2021).

CO₂ deriving from deep geological sources is a carbon source for chemolithoautotrophs that assimilate a fraction of released dissolved inorganic carbon (DIC) into biomass (Barry et al. 2019, Fullerton et al. 2021). Different microbial carbon-fixation pathways might be present in deep subsurface environments depending on the energy sources and chemical composition of the fluids (Lau et al. 2016, Momper et al. 2017, Probst et al. 2017, Fullerton et al. 2021). Predictions from a metagenomic study indicated that the reverse tricarboxylic acid (rTCA) cycle encoding chemolithoautotrophic microbes act as a sink for inorganic carbon in deep subsurface hydrothermal fluids and abundances of rTCA encoding genes positively correlate with increasing DIC and DOC concentrations (Fullerton et al. 2021), suggesting that higher DIC concentrations stimulate chemolithoautotrophically derived in situ DOC production. In pure cultures of sulfate reducing bacteria, autotrophic growth conditions (Mall et al. 2018, Nunoura et al. 2018) and high partial pressure of CO₂ (Steffens et al. 2021) have been found to stimulate carbon fixation via the recently discovered reverse oxidative tricarboxylic acid (roTCA) cycle pathway, indicating the potential for DIC concentrations to increase carbon fixation at low pH.

In neutral and slightly alkaline hydrothermal fluids with pH 7–9, the majority of the DIC is not in the form of CO₂ but rather bicarbonate (HCO₃⁻) (Amend and Shock 2001). In many deep subsurface settings, the fluids are neutral-alkaline and, therefore the DIC available to carbon fixing microbes is in the form of bicarbonate, not CO₂ (Amend and Shock 2001). However, experimental studies investigating the relationship between DIC concentrations and microbial carbon fixation in alkaline subsurface hydrothermal fluids are lacking.

In order to provide insights into whether higher DIC concentrations could stimulate primary producers in alkaline hydrothermal subsurface fluids, we investigated the effect of DIC concentration on microbial carbon fixation in alkaline hydrothermal fluids from a 750-m deep borehole. A positive or negative effect of DIC concentrations on biological carbon fixation could influence the role of microbial communities to act as a subsurface biological carbon sink: the microbial conversion of DIC into organic carbon (either particulate organic carbon or DOC, or both) via carbon-fixation pathways. To this end, we used ¹³C-labeled bicarbonate in quantitative stable isotope probing (qSIP) experiments. The results showed that microbial taxa encoding the Calvin–Benson–Bassham (CBB) cycle contributed most to carbon fixation across bicarbonate concentrations spanning 1–10 mM, with slightly reduced carbon fixation at higher DIC concentrations. In contrast, chemolithoautotrophic bacteria affiliated with the Aquificae encoding the rTCA cycle exhibited increased carbon fixation at higher bicarbonate concentrations, whereas all other microbes encoding different carbon-fixation pathways had either reduced, or no change, in carbon fixation at higher bicarbonate concentrations. These experimental results from 750 m deep alkaline hydrothermal fluids support predictions (Barry et al. 2019, Fullerton et al. 2021) of increased carbon fixation by rTCA encoding chemolithoautotrophic bacteria at high inorganic carbon concentrations in the deep subsurface.

Materials and methods

Study site information and fieldwork sampling for microbial biomass and geochemical analyses

The Biga Peninsula is in the northwest of Turkey and is surrounded by the sea of Marmara in the north and Aegean Sea in the west and south (Fig. S1, Supporting Information) (Okay et al. 1991). Samples were obtained from the Kazdagi Thermal Resort and SPA hotel, which is located in Bardakçılar village in the Kaz Mountain, Biga Peninsula, Turkey (Mount Ida) (Fig. S1, Supporting Information). The hotel facility sits on the Middle Miocene İlyasbaşı formation (Fig. S1, Supporting Information). The reservoir bedrock of the Bardakçılar geothermal field consists of carbonate rock from the Permo-Triassic Kazdag massif. The fluids rise through Upper Oligocene–Lower Miocene Hallaçlar volcanic rocks (Toh), which consists of andesite, volcanic tuff, dacite, rhyodacite, and anglomerate (Akkuş et al. 2005). Granodioritic intrusions due to Oligocene–Miocene volcanic activity are outcropped in the studied area (Sarp et al. 1998). North–east–south–west-trended faults in Bardakçılar geothermal field as well as in nearby geothermal fields in Can and Etili are the main source of the heat transfer (Fig. S1, Supporting Information). The geological units in this area of drilling consist of andesite and volcanic tuffs (Sarp et al. 1998).

Alkaline hydrothermal fluids were collected from a 750-m well (Fig. S1, Supporting Information). The hotel uses two different pumps located 75 m below the ground and the temperature of the deeply sourced fluids was measured at the time of sampling to be 62°C. Electrical conductivity (Ec) (WTW, Cond 3110 Set 1, Weilheim, Germany), temperature, salinity (WTW, Cond 3110 Set 1), pH (WTW, pH 3110 Set 2), and dissolved oxygen (WTW Oxi 3310) were measured in the field using handheld probes, in fluids filled into glass bottles. A volume of 10–20 l of water directly from borehole were filtered through 0.22 µm hydrophilic polyethersulfone (PES) filters (Millipore Express, Merck, Darmstadt, Germany) for 20 min at a medium speed using a peristaltic pump (Masterflex E/S 07571–05, Cole Parmer, USA). Prior to sampling, the fluids from the 750-m well were flushed for 20 min until the fluid temperature reached a steady temperature of 62°C. This effort was made since according to the United States Geological Survey (USGS) standards, it is recommended to flush the well before sampling. Four replicates of PES filters were collected for the microbial analyses. Filters were stored in 15 ml sterile falcon tubes and immediately placed on dry ice, shipped back to Germany, and stored at −80°C until DNA extractions were performed at the Geomicrobiology Lab at the University of Munich (LMU Munich). Geochemical analysis of the fluids was performed as described previously (Einsiedl et al. 2020, Coskun et al. 2023). Briefly, samples for laboratory-based measurements of major anions and cations concentrations and water isotopes (δ²H and δ¹⁸O) were collected in 1.5 ml glass flasks after filtering with 0.2 µm PES filters. The pH of the samples collected for the major anions was adjusted to pH 10 while samples collected for the major cations were adjusted to pH 2. These samples were kept at 4°C. Samples for analysis of DOC concentrations were collected in 50 ml glass bottles after filtering through 0.45 µm polyvinylidene difluoride (PVDF) filters. Samples for DOC analysis were immediately acidified (pH = 2) to avoid changes in organic matter due to biological activity and to remove the inorganic carbon. Samples were stored at 4°C. Samples for the concentrations and isotope analysis of methane (δ¹³C) were transferred into 200 ml glass vials without headspace and sealed with crimped butyl stoppers and stored at 4°C. Laboratory measurements were done in the hydrogeology department of the Technical University of Munich (TUM) and the General Directorate of Mineral Research and Exploration of Turkey (MTA).

Experimental setup for SIP incubations and DNA extraction

For incubation experiments, 5 l borosilicate bottles (VWR, Germany) were filled with 62°C geothermal water and sealed gas-tight with rubber stoppers to prevent intrusions of atmospheric oxygen (O₂). Three 5 l bottles containing hydrothermal fluid were amended with either ¹³C-labeled 1, 5, or 10 mM sodium bicarbonate (99% ¹³C-content, Sigma-Aldrich, Darmstadt, Germany). As a controls, three additional 5 l bottles were filled with hydrothermal fluid and amended with unlabeled bicarbonate at the same concentrations (six 5 l bottles in total). The large (5 l) volumes were necessary in order to be able to collect enough microbial biomass for DNA extraction and qSIP, considering the ultra-low biomass in the starting fluids (Fig. 1B). Attempts at qSIP using smaller volumes failed due to insufficient microbial biomass being collected on the filters. The larger volumes (5 l), however required relatively large bottles that resulted in insufficient space in the hydrothermal heating pool for replicate bottle incubations (biological replicates). For this reason, three technical replicates of density gradient centrifugation and fractionation were performed for DNA extracted from each bottle (see below). While we acknowledge that biological replicates are ideal for any experiment, we also note that technical variation in ¹³C-EAF from ultracentrifugation/density fractionation is often similar to, or greater than, biological ¹³C-EAF variation between incubations in replicate bottles (Coskun et al. 2018, 2019).

Figure 1.

Organic and inorganic carbon concentrations and 16S rRNA gene concentrations (A) A positive correlation between DOC and DIC in hydrothermal fluids at multiple sites in Biga Peninsula (see Fig. S1, Supporting Information, for map). The 750-m deep borehold (Bardakçılar geothermal field) chosen for the qSIP experiment is indicated. (B) Concentration of 16S rRNA genes at the beginning and end of the qSIP incubations, using hydrothermal fluids obtained from the 750-m deep borehole at Bardakçılar geothermal field.

The bottles containing the fluids (and added ¹³C-bicarbonate) were incubated in a hot water tank that was constantly filled with hydrothermal fluids (62°C) from the borehole to maintain a consistent temperature, which was measured throughout the incubation and determined to be steady at 59°C. The gas and water tight sealed bottles were completely submerged under the geothermal water, and had no contact with the atmosphere over the entire course of the incubation. The experiments were terminated after 135 h (5 days and 15 h).

The pH, Ec, salinity, and dissolved O₂ measurements were conducted before and after the incubation (Supplemental Data S1). The initial O₂ concentration of the hydrothermal fluids was in the low dysoxic range (<60 µM) (Supplemental Data S1). We measured dissolved O₂ and redox at the end of incubation and determined the conditions to be slightly reducing (+4 mV) and in the dysoxic to suboxic range (33.9 +/− 2.6 µM O₂, Supplemental Data S1), conditions that promote anaerobic respiration of nitrate as terminal electron acceptor (Wright et al. 2012). Sulfide could be smelled in the fluids, a further indication of reducing conditions.

No acidification treatment was done to convert the initially existing bicarbonate (1.5 mM) to CO_2. Therefore, considering the dilution from the naturally existing (unlabeled) 1.5 mM bicarbonate concentration in the hydrothermal fluids, the final percentage of ¹³C-bicarbonate concentrations should have been roughly 40%, 77%, and 87% in the 1 mM, 5 mM, and 10 mM incubations, respectively. No additional terminal electron acceptors or electron donors were added to the incubations. Therefore, the only energy yielding substrates available to the carbon fixing microbes were those naturally occurring in the fluids. At the end of the incubation, the fluids were immediately filtered onto 0.22 µm filters using a peristaltic pump (Masterflex E/S 07571–05) and frozen immediately after filtering using dry ice in the field.

DNA was extracted from (t₀) and SIP experiments as described previously (Orsi et al. 2022). In brief, for each 15 ml falcon tube holding the filters, the contents of four 2 ml Lysing Matrix E tubes (MP Biomedicals, Solon, OH, USA) were added. A volume of 4 ml of a sterile-filtered sucrose ethylene–diaminetetraacetic lysis buffer (0.75 M sucrose, 0.05 M Tris-Base, 0.02 M ethylenediaminetetraacetic, 0.4 M NaCl, 4 ml 10% sodium dodecyl and pH 9.0) was added to the tubes and then, beat beating was performed for 40 s using a Fast-Prep 5 G homogenizer (MP Biomedicals) at a speed of 4 m/s. Samples were subsequently heated for 2 min at 99°C. After heating, 25 ml of 20 mg/ml proteinase K was added, and tubes were incubated at 55°C overnight. DNA was extracted and purified from the lysate using the DNeasy Blood and Tissue Kit (Qiagen). DNA was suspended with 200 µl diethyl pyrocarbonate (DEPC)-treated (sterile, nuclease free) water. Extracted DNA was quantified fluorometrically using Qubit 3.0 fluorometer (Invitrogen, Eugene, OR, USA).

Density-gradient centrifugation and gradient fractionation

DNA samples were prepared for density-gradient centrifugation according to previously defined protocols for qSIP (Hungate et al. 2015, Coskun et al. 2018) with minor modifications. In brief, density gradient centrifugations were carried out in a TLN-100 Optima MAX-TL ultracentrifuge (Beckman Coulter, Brea, CA, USA) near-vertical rotor at 18°C for 72 h at 165 000 × g. A volume of 50 µl of DNA was added to a solution of cesium chloride (CsCl) and gradient buffer (0.1 M Tris, 0.1 M KCl and 1 mM EDTA) in order to achieve a starting density of 1.705–1.715 g/ml in a 3.3-ml polyallomer OptiSeal tubes (Beckman Coulter). For each SIP incubation, we performed three separate ultracentrifuge runs as technical SIP replicates to account for the variation in density gradient fractionation. These three technical SIP replicates were performed using three separate 50 µl aloquots of extracted DNA from each SIP incubation. The variation between the density gradient replicates (controls and experiments, at each bicarbonate concentration) can be observed in Fig. 2. We focused our replication efforts on the technical associated with density gradient fractionation, because the technical variation introduced at this step tends to be comparable with, or greater than biological variation between replicate SIP incubations (Coskun et al. 2018, 2019).

Figure 2.

Increased DNA buoyant density of 16S rRNA genes in ¹³C-bicarbonate incubations relative to unlabeled controls, at 1 mM (A), 5 mM (B), and 10 mM (C) bicarbonate concentrations. Three separate ultracentrifugation and density-gradient fractionation replicates were performed for DNA extracted from each ¹³C (solid lines, filled symbols) and control (dashed lines, open symbols) bottle incubation at the three bicarbonate concentrations. The plots show the normalized relative abundance of 16S rRNA genes measured by qPCR (y-axis) as a function of density (x-axis) in three replicate ultracentrifugation tubes. The distribution of 16S rRNA gene densities in each replicate tube provided an estimate of the average DNA density and 90% CI (n = 3), which is represented as diamonds and error bars above each panel. The shaded area represents the density range for fractions that were selected for high-throughput Illumina sequencing of 16S rRNA genes for qSIP.

After ultracentrifugation, the density gradients were fractionated into 20 equal fractions of 165 µl from the bottom of polyallomer OptiSeal tubes by using a syringe pump and fraction recovery system (Beckman Coulter). The density of these fractions was measured with an AR200 digital refractometer (Reichert Analytical Instruments, Depew, NY, USA). DNA was precipitated from the fractions using two volumes of polyethylene glycol with 2 µl (10 mg/ml) glycogen and precipitated overnight at room temperature. DNA was pelleted by centrifugation (13 000 × g; 40 min), washed with 70% ethanol, and resuspended with 30 µl Milli-Q ultrapure water (Smart2Pure, Thermo Scientific). DNA was quantified fluorometrically using a Qubit 3.0 fluorometer (Thermo Scientific).

qPCR of 16S rRNA gene copies

Universal primers targeting the V4 hypervariable region of 16S rRNA genes were used in quantitative PCR (qPCR) to determine density shifts in the peak DNA of buoyant density for each incubation. We used a version of the 515F primer with a single-base change (in bold) to increase the coverage of certain taxonomic groups (515F-Y, 5′-GTGYCAGCMGCCGCGGTAA) (Parada et al. 2016). qPCR reactions were carried out as described previously (Coskun et al. 2019). The quantified 16S rRNA genes from each density fraction were plotted against their measured densities for each of the three density gradient fractionation replicates (n = 3 per sample). From each replicate, the DNA from 10 to 12 density fractions was selected for 16S rRNA gene sequencing for qSIP (Fig. 2, shaded areas). Two 16S PCR amplicons from each density fraction (technical replicates to reduce PCR bias) were pooled and subjected to dual-indexed barcoded sequencing of 16S rRNA gene amplicons on the Illumina MiniSeq as described previously (Pichler et al. 2018).

Bioinformatic and qSIP analysis

The sequence reads were processed as previously described (Pichler et al. 2018). In brief, reads were quality trimmed and assembled using USEARCH version 11.0.667 with the default parameters (Edgar 2010) resulting in 28.4 million quality checked 16S rRNA V4 reads. Reads were then de novo clustered at 97% identity using UPARSE; OTUs represented by a single sequence were discarded (Edgar 2013). Taxonomic assignments were generated by QIIME 1.9.1 (Caporaso et al. 2010) using the implemented BLAST method against the SILVA rRNA gene database release 132 (Quast et al. 2013). The level of contamination in each density fraction for qSIP analysis was determined using DNA sequences of dust samples collected from three different laboratories in our building as well as DNA extraction blanks. The common contaminant genera such as Pseudomonas, Ralstonia, Variovorax, or Streptococcus were also removed as these are common contaminants of molecular reagent kits (Salter et al. 2014). Next, only OTUs greater than 10 sequences in total in each replicate for the control and SIP-labeled fractions were selected for further study (Coskun et al. 2018) since low abundance taxa may cause large variations in qSIP calculations (Morrissey et al. 2016). Absolute abundance of OTUs determined by qPCR normalized relative abundance (fractional 16S rRNA gene sequence abundance) are used to calculate shifts in OTU specific DNA buoyant density across the density gradients. For the EAF calculations, the determining factor for how much ¹³C is assimilated by an OTU is the quantitative increase in DNA buoyant density, a calculation, i.e. enabled by the qPCR normalization of fractional 16S sequence read abundance across the CsCl gradient. Therefore, biases imposed by the PCR primers (as is the case for all PCR based studies of microbial diversity), the shifts in OTU specific DNA buoyant density are less biased by this because they are quantitatively scaled against the internal standards for qPCR. The excess atom ¹³C fractions (EAF) values were calculated for the 16S rRNA genes corresponding to OTUs according to previously described study (Hungate et al. 2015) using a qSIP workflow embedded in the HTS-SIP R package (Youngblut et al. 2018). An OTU was considered as an ¹³C-incorporator if the lower CI was greater than zero (Hungate et al. 2015). Statistical analyses and plots were performed using RStudio Version 3.3.0 (Team 2015). The nonmetric multidimensional scaling (NMDS) using Bray–Curtis matrices was generated in vegan package (Oksanen 2010) and pairwise Tukey’s Honest Significant Difference (HSD) post hoc test was performed to distinguish dissimilarities among the samples. Furthermore, the studied geothermal field was compared to other hydrothermal locations in the vicinity using already published data (Coskun et al. 2023). The 16S sequence data were entered in the NCBI Short Read Archive under BioProject ID PRJNA837050.

Metagenomic analysis of ¹³C-labeled fractions

Based on the shift in DNA buoyant density of 16S rRNA genes between the experimental treatments (which indicate ¹³C-labeling of DNA), ¹³C-labeled SIP fractions resulting from the density gradient ultracentrifugation were selected for metagenomic shotgun sequencing (Fig. S5, Supporting Information). Metagenomic libraries were prepared using Nextera XT DNA Library Prep Kit (Illumina) and following the manual provided by the manufacturer. Quality control and quantification of the metagenomics libraries were done on an Agilent 2100 Bioanalyzer System using high sensitivity DNA reagents and DNA chips (Agilent Genomics). Metagenomic libraries were diluted to 1 nM and pooled together to sequence on Illumina MiniSeq platform. SqueezeMeta (Tamames and Puente-Sánchez 2018) and Anvi’o snakemake workflow (Köster and Rahmann 2012, Eren et al. 2015) were used for downstream analysis using coassembly mode with default settings. In short, SqueezeMeta workflow deployed Trimmomatic for adapter removing, trimming and quality filtering by setting the parameters: leading = 8, trailing = 8, slidingwindow = 10:15, and minimum length = 30 (Bolger et al. 2014). Contigs were assembled using Megahit assembler using minimum length 200 nucleotides (Li et al. 2015). Open reading frames (genes and rRNAs; ORFs) were called using Prodigal (Hyatt et al. 2010), rRNAs genes were determined using barrnap (github.com/tseemann/barrnap). Diamond software (Buchfink et al. 2015) was deployed to search for gene homologies against the databases Genbank nr database for taxonomic assignment, eggNOG v4.5 (Huerta-Cepas et al. 2016) and KEGG database (Kanehisa and Goto 2000). Cutoff values for assigning hits to specific taxa were performed at e value as e⁻³, and a minimum amino acid similarity of 40 for taxa and 30 for functional assignment, which were default settings of SqueezeMeta. Bowtie2 (Langmead and Salzberg 2012) was used to map the read onto contigs and genes. Anvi’o snakemake workflow (anvi’o v7) was used to bin and to refine the MAGs. MaxBin (Wu et al. 2016), concoct (Alneberg et al. 2014), and Metabat2 (Kang et al. 2019) were deployed for the binning and DAS Tool (Sieber et al. 2018) was deployed to choose the best bin for each population. The completeness and contamination of the bins were checked in anvi’o v7. Accordingly, 8 high-quality, 13 medium-quality, and 3 low-quality MAGs were constructed and used for the analysis (Supplemental Data 3). KEGGdecoder (Graham et al. 2018), anvi’o and COG annotations were used to estimate the completeness of the predicted metabolisms. MAGSIPTR2 and MAGSIPTR39 were published previously with MAG names Candidatus Bipolaricaulota CK101 and CK84 (Coskun et al. 2023). The metagenomic sequences were stored in the NCBI Short Read Archive under BioProject ID PRJNA837050. Metagenomic dataset and intermediate files to produce qSIP results were deposited under https://figshare.com/authors/_mer_Coskun/9725927.

Phylogenetic tree construction

16S rRNA gene phylogenetic trees (Fig. S8, Supporting Information) were constructed using GTR algorithm with 100 bootstraps in Seaview (Gouy et al. 2010). In brief, rRNA and single copy genes were extracted from the MAGs using hidden Markov model (HMM) sources of barnapp and “Bacteria_71” using the “anvi-run-hmms” command in anvi'o (Eren et al. 2015). Ribosomal phylogenetic tree (Fig. S8C, Supporting Information) was constructed using Seaview (Gouy et al. 2010) with five concatenated ribosomal proteins (L18, L3, L4, S11, and S3) extracted from MAGs and closest complete or incomplete genomes to our MAGs determined by BLAST searches against NCBI-nr. For the taxonomic assignment of the MAGs that contained a 16S rRNA gene, we used the 16S rRNA gene phylogeny for taxonomic assignment. For those MAGs that did not contain a 16S rRNA gene, we used the anvi’o SCG classification that was based on an alignment of five concatenated ribosomal proteins and the corresponding phylogenetic tree (Fig. S8, Supporting Information).

Linking ¹³C-labeled OTUs to MAGs

We linked the ¹³C-labeled OTUs to MAGs using phylogenetic analyses (Fig. S8, Supporting Information). MAGs were grouped into two categories (1) those encoding a binned 16S rRNA gene, and (2) those not encoding a 16S rRNA gene. For those MAGs encoding a binned 16S rRNA, we were able to directly link the identity of the MAG to its corresponding 16S rRNA gene within the qSIP dataset using 16S rRNA gene phylogenies (Fig. S8A, Supporting Information). For all MAGs that contained a 16S rRNA gene (n = 11), a ¹³C-labeled 16S OTU detected in qSIP could be directly linked to the MAG. For those MAGs that did not contain a binned 16S rRNA gene (n = 12), we indirectly linked the identity of the MAG to its corresponding ¹³C-labeled 16S rRNA gene OTU. This was done by comparing the MAG ribosomal protein phylogeny to a 16S rRNA gene phylogeny of the same affiliated named species and “Candidatus” organisms (Fig. S8B, Supporting Information). This allowed us to indirectly link MAGs lacking a 16S rRNA gene with their most likely ¹³C-labeled OTUs in the qSIP dataset (Fig. S8B, Supporting Information). The confidence in the assignment between 16S and MAG taxonomy comes from (1) the relatively short phylogenetic distance to the same organism (or closely related organism) and (2) bootstrap support values on the nodes connecting the environmental sequences to named taxa at a relatively short distance, most of which were >90% (Fig. S8B, Supporting Information). Once all MAGs (n = 24) were linked (either directly or indirectly) to the labeled OTUs, we were able to compare the gene content of the MAGs with their ¹³C-assimilation across the three tested bicarbonate concentrations (Fig. 2).

Calculating absolute ¹³C-labeled carbon of MAGs linked with ¹³C-labeled OTUs

We grouped ¹³C-incorporating MAGs into autotrophs, mixotrophs (or facultative autotrophs), and heterotrophs based on our own criteria in order to trace carbon flows between these functional groups. Potential autotrophs were grouped as ¹³C-assimilating MAGs with a carbon-fixation pathway >80% complete and <2 ABC-type transporters for utilization of organic matter (n = 4). Potential mixotrophic or facultative autotrophic taxa (n = 11) were defined as ¹³C-assimilating MAGs with a carbon-fixation pathway >80% completeness and >1 ABC-type transporters for utilization of organic matter. Heterotrophic taxa were identified as ¹³C-assimilating MAGs that did not encode a carbon-fixation pathway, or encoded <20% of a carbon-fixation pathway (n = 9). We could group all ¹³C-assimilating MAGs into one of these three categories (potential autotroph, potential mixotroph, and heterotroph) according to these criteria. Criteria for defining heterotrophy, mixotrophy, and autotrophy may differ from study to study and because the criteria we use are specific to our study they are arbitrary. However, we believe that our definition is relatively robust because it is based on (1) the genomic potential for utilization of external organic substrates and CO₂ fixation, and (2) quantified taxon-specific CO₂ fixation activity (from qSIP).

The absolute ¹³C-labeling in MAGs was calculated using the following steps. First, each MAG was linked to a ¹³C-labeled OTU (Fig. S8, Supporting Information) that has a corresponding ¹³C-EAF value (Fig. S7, Supporting Information). Second, the MAG-specific ¹³C-EAF values were normalized against the absolute abundance of the MAG-specific 16S rRNA gene copies (from qPCR). Third, the absolute abundances of MAG-specific ¹³C-EAF was summed for each of the four carbon-fixation pathway groups (rTCA, WLP, CBB, and 3-HP) to determine the total amount of ¹³C-fixed by microbes containing each carbon-fixation pathway. For a detailed summary of how these calculations were made, please consult Supplemental Data S4. For MAGs categorized as heterotrophs (lacking a CO₂-fixation pathway), absolute ¹³C-assimilation was calculated using the same steps. Technical variation associated with the MAG-specific ¹³C-EAF were assessed using the lower and upper bound of 90% confidence intervals of ¹³C-EAF measurements from the qSIP results of the corresponding 16S rRNA gene OTUs (Fig. S6, Supporting Information).

Results and discussion

Correlation between DOC and DIC in hydrothermal fluids

We sampled fluids for geochemical measurements from numerous hydrothermal aquifers ranging from surface hot spring pools to fluids from 1350-m subsurface deep boreholes, across the Biga Peninsula (Coskun et al. 2023) (Fig. S1, Supporting Information). This revealed a positive relationship between DIC and DOC in the greater subsurface aquifer system (Fig. 1A) and is consistent with prior observations in deep hydrothermal fluids from the Costa Rican convergent margin (Barry et al. 2019, Fullerton et al. 2021) indicating that DIC is potentially related to increased primary production in the hydrothermal fluids and autochthonous DOC production. However, some of the DOC could be derived from mixing with allochthonous shallow surface sources, and it is possible that high DIC could be the result of heterotrophic activity under high DOC conditions. Moreover, the aquifer is associated with carbonate bedrock, which is likely a source of higher DIC in the fluids sampled at some locations (Yalcin 2007) (Fig. S1, Supporting Information). The relation between DIC, DOC, and primary production in hydrothermal fluids is poorly understood, as are the rates of primary production and DOC production in hydrothermal systems (Le Bris et al. 2019, Lai et al. 2023). Thus, we investigated the effect of DIC concentration on microbial carbon fixation in alkaline subsurface hydrothermal fluids from deep subsurface (750-m) fluids from one of the sampling locations in the Bardakçılar geothermal field (Fig. 1A; Fig. S1, Supporting Information).

Our target location for the qSIP experiments was a 750-m deep borehole with alkaline (pH = 8.27), hydrothermal (62°C) fluids containing a DIC concentration of 1.5 mM. This DIC concentration was at the lower end of the bicarbonate concentrations that exists within the greater aquifer system that we sampled (from 1.5 to 11 mM; Fig. 1A) (Yalcin 2007). By incubating the fluids with ¹³C-bicarbonate at three different concentrations characteristic of the natural state (1, 5, and 10 mM), we were able to investigate effects of higher DIC concentrations on primary producers sampled from subsurface fluids with a lower DIC concentration. An increased carbon fixation by particular microbial primary producers at higher DIC would help explain the DIC and DOC correlation observed across the aquifer system (Fig. 1A). If DIC concentrations influence biological carbon fixation it would highlight a potentially important factor controlling how microbial communities act as a subsurface biological carbon sink by converting DIC into particulate and DOC.

Biogeochemical characterization of the deep alkaline fluids

Prior to performing the qSIP experiments, geochemical and microbiological samples were obtained over a 3-year period (2019, 2020, and 2021) from the 750-m deep borehole (Fig. S1, Supporting Information). Untreated, pristine alkaline hydrothermal fluids were collected directly from the borehole (see the section “Materials and methods”) and measured to have an alkaline pH (pH = 8.27) with a temperature of 62°C and concentration of 1.5 mM of DIC (Fig. 1A). At this alkaline pH, the majority of DIC is in the form of bicarbonate (Amend and Shock 2001). The study site bedrock is metamorphically altered carbonate rocks (marble) embedded in Kazdağ massif (Yalcin 2007) (Fig. S1, Supporting Information). This carbonate aquifer becomes sodium and sulfate-enriched while rising through the Upper Oligocene–Lower Miocene volcanic rocks (Akkuş et al. 2005).

The fluids were depleted in nitrate and nitrite (Supplemental Data S1), in line with previous reports (Yalcin 2007). The concentration of DOC was relatively low (0.8 mg/l) and reflects a characteristic of deep aquifer systems that are generally oligotrophic (Lopez-Fernandez et al. 2018, Barry et al. 2019, Probst et al. 2020, Fullerton et al. 2021), where biological production is relatively limited. Indeed, the microbial 16S ribosomal RNA (rRNA) gene concentrations in the alkaline hydrothermal fluids from the 750-m deep aquifer were extremely low, with 1.8 × 10² (± 33) copies per ml (Fig. 1B). Our results are in line with previous continental subsurface studies from a depth of 500–800 m where similarly low 16S rRNA gene copy concentrations (10¹–10⁴ copies per ml) were detected (Lavalleur and Colwell 2013, Bomberg et al. 2016, Magnabosco et al. 2018). The 16S rRNA gene concentrations measured were consistently higher than our detection limit of 10² total 16S rRNA gene copies per contamination control (DNA extraction blanks and quantitative PCR no template controls) (Coskun et al. 2022), showing that we detected the natural microbial community emanating from the deep hydrothermal fluids with a minimal influence of contamination.

The low microbial biomass can be explained by the ca. 10-fold lower concentration of DOC (Supplemental Data S1) compared to other surface-influenced continental subsurface aquifers (Lopez-Fernandez et al. 2018), because DOC availability is known to shape microbial communities in continental subsurface crustal ecosystems (Lopez-Fernandez et al. 2018, Westmeijer et al. 2022). Here, we also acknowledge that our filter size (0.22 µm) might have discarded the low-nucleic acid microbial cells smaller than our filter size (in other words, ultramicrocells; Duda et al. 2012), leading to potential underestimation of microbial biomass such as members of Candidate Phyla Radiation (Luef et al. 2015).

Over the three sampling years, the initial microbial communities (t₀) were dominated by different taxonomic groups (Fig. S2, Supporting Information) and were markedly different from one another (Fig. S3A, Supporting Information). In 2019, the microbial community was dominated by operational taxonomic units (OTUs) affiliated with Sulfurihydrogenibium (Aquificae; 70% relative abundance) and Deinococcus–Thermus (23%) whereas samples obtained in 2020 were comprised of OTUs belonging to Deinoccoccus–Thermus (31%), Chloroflexi (31%), and Hadesarchaeaeota (7%) (Fig. S2, Supporting Information). In the samples collected in 2021, Gammaproteobacteria (53%), Nitrospirae (20%), and Aquificae (12%) predominated the initial microbial community of the subsurface alkaline hydrothermal fluids, which were used for the stable isotope probing (SIP) incubations. This overall microbial community structure is similar to thermophilic microbial communities found in other alkaline hot springs such as those in the Tibetan Plateau (Wang et al. 2013), Yellowstone National Park (Schubotz et al. 2015), Costa Rican convergent margin (Fullerton et al. 2021), deep-subsurface aquifers (Spanevello and Patel 2004), and hydrothermal fluids in Bardakçılar geothermal field (Coskun et al. 2023).

The oxidation–reduction potential (ORP) of the fluids was +4 mV, which falls with the ORP range that promotes denitrification via anaerobic nitrate respiration (Fuerhacker et al. 2000). The initial dissolved O₂ concentration was in the subhypoxic range (51.3 µM) (Supplemental Data S1), and a strong odor of hydrogen sulfide (H₂S) was evident during sampling indicating the presence of H₂S in the fluids. Methane was detectable in the fluids, albeit at relatively low concentration of 0.079 (mg/l), and the δ¹³CH₄ was −10.4 ‰ (Supplemental Data S1) indicating that the methane was likely derived from thermogenic pyrolysis of organic matter (Schoell 1988). The δ²H and δ¹⁸O composition of the fluids suggested a meteoric origin of the water as the isotopic values ranged between both Global and Eastern Mediterranean meteoric water lines (Craig 1961) (Fig. S4, Supporting Information). Moreover, the δ¹⁸O values were within the range of annual average δ¹⁸O of precipitation and shallow groundwater (−6‰ to −9‰) from nearby sites in Turkey (Jasechko 2019). This suggests that the origin of the subsurface hydrothermal fluids in our tested ecosystem is mostly modern meteoric waters and not affected by saltwater intrusion, consistent with prior studies (Yalcin 2007). A meteoric water origin also suggests the possibility that sites of high DOC in the fluids (Fig. 1A) could be due to allochthonous sources of organic matter.

¹³C-labeling of OTUs in qSIP analyses varies across substrate concentrations

One of the major challenges in microbial ecology is to link the identities of microorganisms with their quantitative ecological roles (Orsi 2023). One way to do this is to use quantitative DNA stable-isotope probing methods (Radajewski et al. 2000, Neufeld et al. 2007). Specifically, quantitative DNA stable isotope probing (qSIP) uses labeled carbon (¹³C) substrates to quantify the amount of carbon assimilated by all detectable microbial taxa in an environmental sample (Hungate et al. 2015). The qSIP methodology with ¹³C-labeled substrates has been applied to broad range of environments such as soil (Hungate et al. 2015, Morrissey et al. 2016, Starr et al. 2021), marine sediments (Orsi et al. 2020), seawater (Orsi et al. 2021), ocean crust (Coskun et al. 2022), and geothermal terrestrial springs (Lai et al. 2023). Here, we performed qSIP experiments by adding the H¹³CO₃⁻ at three different concentrations (1, 5, and 10 mM) to test the effect of naturally available range of bicarbonate concentrations that the microbial communities are exposed to as they are transported through the subsurface aquifer (from 1 to 11 mM; Fig. 1A). No additional energy yielding terminal electron acceptors or electron donors were added to the incubations. Therefore, the only energy yielding substrates available were those naturally occurring in the fluids.

The final concentration of ¹³C-bicarbonate in the incubations was affected by the naturally occurring concentration of bicarbonate in the fluids (1.5 mM, >99% naturally abundant ¹²C-bicarbonate, see the section “Methods”) (Supplemental Data S1). We chose not to acidify the samples in order to remove the naturally abundant starting DIC concentrations, because the acidic pH change would affect the microbial community composition living in the alkaline fluids. However, the addition of bicarbonate was associated with the modest changes in pH which decreased from pH 8.27 to 7.6 (± 0.3) (Supplemental Data S1). The equilibrium between the carbonate species could be slightly altered by this pH reduction, but the available DIC species for the microorganisms should have been mostly bicarbonate (Amend and Shock 2001, Dettori and Donadio 2020). The overall redox state of the qSIP incubations was relatively low (+4 mV), and O₂ concentrations at the end of the incubations were in the subhypoxic range (33.9 +/− 2.6 µM) (Supplemental Data S1). We acknowledge that our incubation conditions of the hydrothermal fluids do not reflect the exact physicochemistry of the in situ deep aquifer environment, and that the results obtained are a function of the incubation conditions and may not completely reflect the exact ecosystem of the subsurface aquifer that can vary in terms of physicochemistry.

After 5 days, we observed that the amount of 16S rRNA gene copies in all incubations increased roughly 10-fold (Fig. 1B). Moreover, increased DNA buoyant density of 16S rRNA genes was observed in all three ultracentrifugation replicates for each ¹³C-bicarbonate concentration tested (Fig. 2) demonstrating assimilation of bicarbonate by the microbial communities. The microbial community structures in the incubations were found to have shifted from the initial microbial composition from 2021, but this change was not significant (Tukey’s HSD, P = .80–.99) (Fig. S3A, Supporting Information), indicating that physicochemical conditions in our SIP incubations did not have drastic effects on the microbial community composition. The qSIP analysis revealed which OTUs were labeled in the presence of ¹³C-bicarbonate at the different concentrations (Fig. S6, Supporting Information). The assimilation of ¹³C-bicarbonate was measured in the qSIP method via EAF ¹³C in the 16S rRNA genes of OTUs, whereby ¹³C-EAF values were reported on a scale between 0.0 and 1.0 with 1.0 representing 100% labeling of carbon atoms in the 16S rRNA gene (Hungate et al. 2015).

The mean EAF values of all the ¹³C-labeled OTUs in the incubations showed that the 1 mM incubation exhibited the highest degree of ¹³C-labeling (mean EAF: 0.17 ± 0.01), followed by the 10 mM (mean EAF: 0.16 ± 0.01), and 5 mM (mean EAF: 0.13 ± 0.01) incubations (Fig. S6 and Table S1, Supporting Information; Supplemental Data S2). The 1 mM ¹³C-bicarbonate incubations also had the highest proportion of labeled OTUs (38 OTUs out of total 65 OTUs; 58%), followed by 5 mM (40 OTUs out of total 91 OTUs; 44%) and 10 mM ¹³C-bicarbonate incubations (25 OTUs out of 59 OTUs; 42%) (Table S1, Supporting Information; Supplemental Data S2). A total of 60 OTUs (out of 100 total) incorporated ¹³C-bicarbonate with at least one ¹³C-bicarbonate concentration. Among these ¹³C-labeled OTUs, 24 OTUs assimilated more ¹³C-bicarbonate with 1 mM bicarbonate than the incubations with higher bicarbonate concentrations, whereas 26 and 10 OTUs had higher ¹³C-labeling at the incubations amended with 5 mM and 10 mM bicarbonate, respectively (Fig. S7, Supporting Information).

Among labeled OTUs in the 1 mM bicarbonate incubation, the OTU showing the highest degree of labeling was affiliated with the family Hydrogenophilaceae (Gammaproteobacteria, also assigned to Hydrogenophilalia; Boden et al. 2017) (0.37 mean EAF) (Fig. S6, Supporting Information), i.e. mainly comprised of taxa that are thermophilic chemolithoautotrophs that can gain energy from H₂ oxidation and use DIC for biomass production via the CBB cycle (Orlygsson and Kristjansson 2014). Taxa affiliated with Hadesarchaeaeota, Hydrothermae, and candidate phyla Bipolaricaulota (formerly known as Acetothermia) also assimilated ¹³C-bicarbonate at 1 mM albeit to a relatively lower degree (Figs S6 and S7, Supporting Information). Despite the 1 mM H¹³CO₃⁻ concentration being diluted more than 50% by the naturally abundant ¹²C bicarbonate (1.5 mM; Fig. 1A), the carbon fixation at this diluted concentration was similar compared to the 5 mM ¹³C-bicarbonate incubation, and greater compared to the 10 mM ¹³C-bicarbonate incubation (Fig. 2), where the labeled substrates were available to the microbes at higher concentrations. This shows that most of the microbes in the sampled hydrothermal fluids were adapted to fix carbon at concentrations closer to the in situ bicarbonate concentration. We considered the potential for enzymatic discrimination (due to kinetic isotope fraction effects) against an over-abundance of “heavy” (¹³C) isotopes to influence the results at the 5 and 10 mM concentration. However, biological isotope fractionation occurs in the per-mille range (Kendall and Caldwell 1998) whereas ¹³C-EAF values are in the % range (Hungate et al. 2015) and, therefore an order of magnitude higher. Moreover, ¹³C assimilation was comparable between 1 and 5 mM ¹³C -bicarbonate incubations (Fig. 2), and indicating enzymatic isotope discrimination against higher ¹³C concentrations was minor compared to the overall ¹³C assimilation.

In the 5 mM ¹³C-bicarbonate incubation, the mean OTU-specific EAF values were lower compared to 1 mM (Table S1, Supporting Information). The dominant ¹³C-bicarbonate assimilating OTUs at 5 mM bicarbonate were affiliated with Tepidomonas (Gammaproterobacteria), Thermus (Deinococcus–Thermus), family Hydrogenophilaceae, and Hydrogenobacter (Aquificae) (Fig. S7, Supporting Information; Supplemental Data S2). OTUs affiliated with candidate Phyla such as Sumerlaeota, Hadesarchaeaeota, Candidate class Bathyarchaeia also assimilated ¹³C-bicarbonate at 5 mM (Fig. S7, Supporting Information) showing that Archaea were also contributing to carbon fixation in the subsurface alkaline hydrothermal fluids. An OTU affiliated with the Firmicute genus Moorella, a group of thermophilic acetogens, had one of the highest amounts of ¹³C-bicarbonate incorporation with 5 mM bicarbonate (Fig. S7, Supporting Information).

In the incubations amended with 10 mM bicarbonate, ca. 40% fewer number of ¹³C-bicarbonate assimilating OTUs were observed compared to 1 and 5 mM (Table S1, Supporting Information; Supplemental Data S2). This was consistent with a relatively lower increase in microbial DNA buoyant density at 10 mM ¹³C-bicarbonate, compared to 1 and 5 mM ¹³C-bicarbonate (Fig. 2). Many phyla that were important for carbon fixation at 1 mM bicarbonate concentration (e.g. Acetothermia, Hydrothermae, Elusimicrobia, Crenarchaeota, and Hadesarchaeaeota) were no longer found to assimilate ¹³C-bicarbonate at the 10 mM concentration (Fig. S7, Supporting Information). However, the highest degree of ¹³C-bicarbonate assimilation at 10 mM was found in bacterial taxa affiliated with Aquificae, Thermodesulfovibrio (Nitrospirae), Hydrogenobacter, and Hydrogenophilaceae (Fig. S7, Supporting Information). Many of these taxa including the Aquificae are often found in deep subsurface or hydrothermal habitats (Henry et al. 1994, Takai et al. 2003, Hügler et al. 2007, Schubotz et al. 2015, Fullerton et al. 2021, Lai et al. 2023). The carbon-fixing community at 10 mM bicarbonate concentration was significantly different compared to those in 1 and 5 mM [analysis of similarity (ANOSIM) of 1 and 10 mM: R² = 0.676, P = .0045; ANOSIM of 5 and 10 mM: R² = 0.906, P = .0021; Fig. S7, Supporting Information; Supplemental Data S2), indicating that the concentration of bicarbonate is an important factor that selects for distinct carbon-fixing taxa in the alkaline hydrothermal fluids from this subsurface ecosystem. The activity of the carbon-fixing taxa during the incubations measured with qSIP, is reflected in the net growth of the microbial community by roughly an order of magnitude (Fig. 1B).

Carbon-fixation pathways and energy metabolism in ¹³C-labeled metagenome-assembled genomes

Linking the ¹³C-labeled OTUs to metagenome-assembled genomes (MAGs) using phylogenetic analyses (see the section “Materials and methods”; Fig. S8, Supporting Information) enabled us to compare the gene content of all MAGs (n = 24) with their ¹³C-assimilation across the three tested bicarbonate concentrations (Fig. 3; Fig. S9, Supporting Information). A total of 15 of the DIC-incorporating MAGs were found to encode carbon-fixation pathways, indicating their potential to contribute to primary production (Fig. 3). This included the CBB cycle, Wood–Ljungdahl Pathway (WLP), rTCA cycle, and the 3-hydroxypropionate bicycle (3-HP bicycle) (Fig. 3). The remaining DIC-incorporating MAGs lacked carbon-fixation pathways and therefore displayed the potential for a heterotrophic lifestyle.

Figure 3.

Linking ¹³C-DIC-incorporating MAGs with the encoded carbon-fixation pathways and mixotrophic or facultative autotrophic potential. From top to bottom: ¹³C-labeled MAGs were grouped based on the presence or absence (heterotrophs) of a carbon-fixation pathway. MAGs encoding a carbon-fixation pathway (left side) were sorted based on the increase in mixotrophic or facultative autotrophic potential defined as encoded ABC-transporters for external organic substrates (bottom heatmap). Quantitative assimilation of ¹³C-bicarbonate by each MAG at the three substrate concentrations is shown in the middle plot. Individual points represent EAF values of specific OTUs phylogenetically linked to the MAGs (Fig. S9, Supporting Information). Error bars: 90% CI based on three separate CsCl density gradients (see Fig. 2). The lower heatmaps shows the metabolic pathway (KEGG) completeness of each MAG.

The energy metabolism of many MAGs displays the ability for anaerobic dissimilatory nitrate reduction and utilization of sulfur compounds (S⁰, H₂S) as an electron donor (Fig. S9, Supporting Information), similar to other deep subsurface and hydrothermal environments that often favor a sulfur-oxidizing, denitrifying microbial metabolism (Amend and Shock 2001, Shock et al. 2010, Osburn et al. 2014, Lau et al. 2016). Especially, sulfide oxidation was found to be a dominant potential energy conservation mechanism for most of the detected ¹³C-labeled MAGs (Fig. S9, Supporting Information). This included two rTCA-encoding MAGs affiliated with Aquificae, five CBB cycle encoding MAGs affiliated with Thermus (Deinococcus–Thermus), Meiothermus, Sulfuritortus (Gammaproteobacteria), and Tepidomonas (Gammaproteobacteria) mostly through the presence of sulfide: quinine reductase (sqr) and/or flavocytochrome c (fccAB) (Fig. S9, Supporting Information). Sulfur oxidation through the Sox system and thiosulfate oxidation using tdsA were also detected in several of the MAGs (Fig. S9, Supporting Information; Supplemental Data S5). These data fit with prior studies showing that thermophilic CBB-encoding Thermus and Sulfuritortus bacteria obtain energy by sulfur oxidation and reducing O₂ or NO₃⁻/NO₂⁻ as electron acceptors (Bjornsdottir et al. 2009, Kojima et al. 2017). The ability to use nitrate as terminal electron acceptor in anaerobic respiration is a potential reason why the carbon-fixation activity (¹³C-bicarbonate assimilation) of these MAGs was favored under the suboxic conditions in the fluids and SIP incubations (23–51 µM O₂, Supplemental Data S1), because anaerobic respiration via dissimilatory nitrate reduction is a favored terminal electron accepting process in this oxygen concentration range (Hügler and Sievert 2011, Wright et al. 2012). Nitrate reduction as a terminal electron accepting process is further consistent with the depletion of dissolved nitrate in the hydrothermal fluids (Supplemental Data S1).

We do not know what the primary electron donor is for the microbial communities in the subsurface alkaline hydrothermal fluids, which could range from organic matter, H₂, reduced iron, or sulfur compounds as previously described for hydrothermal systems (Amend and Shock 2001, Shock et al. 2010, Osburn et al. 2014). However, at slightly alkaline pH (7–8) there is an affinity for the reaction of nitrate reduction coupled to sulfur oxidation (Shock et al. 2010). Considering the pH of the studied fluids falls close to this range, the use of sulfur as an electron donor would explain the widespread genomic potential for sulfur oxidation in many MAGs (Fig. S9, Supporting Information).

Assessing the potential for autotrophy, mixotrophy, and heterotrophy in ¹³C-labeled MAGs

Potentially mixotrophic or facultatively autotrophic DIC-incorporating MAGs were identified based on (a) the presence or absence of CO₂ fixation pathways, and (b) an enriched (or lack of) potential for uptake of external organic substrates (e.g. sugars and or amino acids) via genes encoding ABC-transporters and or TonB-dependent transport proteins (Fig. 3). Our genomic and qSIP experimental findings do not discern between mixotrophy (an ability to utilize organics and CO₂ simultaneously) and facultative autotrophy (an ability to switch between the autotrophy and heterotrophy) but rather highlight those ¹³C-labeled MAGs from the studied alkaline hydrothermal aquifer that have the genomic capacity to perform one of these lifestyles. Out of 15 MAGs that encoded partial/complete carbon-fixation pathways, 11 of them also contained relatively higher numbers of ABC-type transporters for classes of organic molecules (Fig. 3), suggesting a mixotrophic or facultative autotrophic lifestyle whereby they can use CO₂ and organic carbon from the environment as major carbon sources for either heterotrophic and/or autotrophic growth. All of the MAGs encoding the CBB cycle, most of the MAGs encoding the WLP cycle and one MAG encoding the 3-HP bicycle displayed a genomic potential for mixotrophic or facultative autotrophic lifestyle (Fig. 3). The WLP encoding Ca. Bipolaricaulota MAGs assimilated more ¹³C-DIC in the incubation amended with 1 mM bicarbonate than that in 5 mM and did not get labeled in 10 mM incubation (Fig. 3), suggesting that their carbon fixation was reduced under the increased DIC concentrations of our experimental conditions. Our results may reflect an adaptation of Ca. Bipolaricaulota to lower DIC concentrations, or an inhibition of carbon fixation at higher DIC concentrations. We are unaware of any other studies that measured carbon-fixation rates of Ca. Bipolaricaulota at different DIC concentrations, so we are not able to compare our results to studies making similar measurements. Nevertheless, our results provide a first indication that some bacteria of the Ca. Bipolaricaulota exhibit more carbon fixation under a relatively lower DIC concentration of 1 mM.

Notably, ¹³C-labeled MAGs encoding the essential genes for rTCA cycle (exclusively Aquificae) had a relatively reduced suite of ABC-transporters for organic molecules suggesting a stricter and possibly obligate autotrophic metabolism (Fig. 3). The only ABC transporter for organic substrates encoded by the Aquificae MAGs was an incomplete ABC-type phospholipid transporter system (MlaABCDEF). This transporter has an important role in maintaining outer membrane lipid asymmetry (Chong et al. 2015) and phospholipid-trafficking between outer and inner membrane (Hughes et al. 2019), suggesting a role in maintaining cell membrane integrity as opposed to transport of external organic carbon sources. Together with the rTCA carbon-fixation cycle metabolism (Fig. S9, Supporting Information), these results indicate Aquificae MAGs were fixing ¹³C-bicarbonate via autotrophy and the rTCA cycle, which is a wide-spread trait in the thermophilic Aquificae (Hügler et al. 2007, Hügler and Sievert 2011, Gupta 2014). Some Aquificae have exhibited mixotrophic metabolism in SIP experiments from hydrothermal springs, via ¹³C-labeling from organic sources such as glucose (Schubotz et al. 2015) and amino acids (Lai et al. 2023). However, the lack of genes encoding transport proteins for sugars or amino acids in the Aquificae MAGs reported here indicate that the Aquificae taxa inhabiting the alkaline hydrothermal fluids of the 750-m borehole lack this capacity.

Total microbial carbon assimilation with different carbon-fixation pathways

We grouped ¹³C-incorporating MAGs into potential autotrophs, mixotrophs (or facultative autotrophs), and heterotrophs based on criteria (see the section “Materials and methods”) in order to trace carbon flows between these functional groups. Our qSIP analysis only captures the microorganisms that were able to assimilate the ¹³C-bicarbonate, or ¹³C-organic matter secondarily deriving from primary production, into their biomass. We considered that anaplerotic reactions can contribute to inorganic carbon assimilation in heterotrophs accounting for 2%–8% of cell carbon (Spona-Friedl et al. 2020). However, many of MAGs defined as heterotrophs here reach ¹³C-assimilation values >8% (Fig. 3). Therefore, anaplerotic carbon fixation alone is insufficient to explain the ¹³C-labeling of these heterotrophic bacteria. We, therefore interpret labeling of heterotrophic MAGs with an EAF >8% as a mixture of anaplerotic reactions and the secondary production of organic ¹³C, that was first created by primary producers (encoding a carbon-fixation pathway) from ¹³C-bicarbonate (Fig. 3).

We estimated the total contribution of each of the four groups of carbon-fixation pathway encoding MAGs (rTCA, WLP, CBB, and 3-HP) to primary production. This was done using qPCR normalized absolute abundances of 16S rRNA genes from ¹³C-labeled MAGs, to convert MAG-EAFs to total ¹³C-assimilation. The total estimated carbon-fixation contributions by each of the four detected carbon-fixation pathway groups (rTCA, WLP, CBB, and 3-HP) was determined as the sum of qPCR-normalized ¹³C-assimilation of all MAGs corresponding to each group (Supplemental Data S4). The same estimation of total carbon assimilation was applied to heterotrophic MAGs (¹³C-labeled but no carbon-fixation pathway), that serves as an estimate for bacterial secondary production from organic ¹³C that was originally produced by facultative or strict autotrophs and or mixotrophs encoding one of the carbon-fixation pathways (Fig. 3).

The CBB cycle dominates carbon fixation in the geothermal fluids

Across all three concentrations of bicarbonate tested, CBB-encoding MAGs had a significantly higher contribution to total carbon fixation by roughly an order of magnitude compared to taxa encoding the rTCA cycle (Two side t-test P = .02) (Fig. 4). The total amount of carbon flux attributed to mixotrophic bacteria encoding CBB cycle was at least 8.8 times higher compared to rTCA encoding taxa and more than an order of magnitude higher compared to WLP and 3HP encoding MAGs (Fig. 4). These differences can be explained by the higher absolute abundances of the mixotrophic bacteria encoding CBB cycle (64%–70% of total 16S rRNA gene copies), compared to MAGs with different carbon-fixation pathways. Therefore, under our tested conditions the CBB-encoding MAGs were the group of carbon fixing microbes that had the greatest contribution to primary production and carbon turnover. CBB cycle in mixotrophic bacteria was shown to be an important trait, which enhances the microbial survivability in oligotrophic groundwater ecosystems (Taubert et al. 2022) and CO₂-rich cold subsurface aquifer system (Probst et al. 2017). Bicarbonate concentrations did not positively correlate with total carbon incorporation of CBB-encoding MAGs, but rather coincided with a decrease in their ¹³C-EAF at 10 mM bicarbonate (Fig. 3, Supplemental Data S4).

Figure 4.

Total carbon assimilation by ¹³C-labeled MAGs encoding different carbon-fixation pathways. The y-axis is the sum of MAG ¹³C-EAF values per group using MAG absolute abundance based on qPCR. Error bars: 90% CI based on three density-gradient fractionation replicates per treatment (Fig. 2). The total ¹³C-assimilation by rTCA-encoding bacteria was positively correlated with bicarbonate concentration (R² = 0.94), no positive correlations were observed for any of the other groups of organisms encoding different carbon-fixation pathways (CBB, WLP, and 3-HP).

Chemolithoautotrophic Aquificae encoding the rTCA cycle increased carbon fixation at higher bicarbonate concentrations

The only carbon-fixation pathway that was associated with a positive correlation between DIC concentration and quantitative assimilation of DI¹³C was the rTCA cycle (Fig. 4). At 10 mM DIC, there was a 3- to 10-fold higher carbon assimilation by rTCA encoding organisms, compared to the 1 mM and 5 mM DIC concentrations (Fig. 4). The technical variation associated with the quantitative ¹³C assimilation by rTCA encoding organisms at 10 mM does not overlap with the variation for estimated carbon fixed at 1 and 5 mM (Fig. 4). Technical variation tends to be comparable to, if not greater than, biological variation between replicate SIP incubations using our qSIP protocol (Coskun et al. 2018, 2019), and therefore differences in quantitative ¹³C labeling of rTCA MAGs between the different DIC concentrations is unlikely to be explained by random variation. We, thus interpret these results as indicating an increased carbon fixation at higher DIC concentration by the rTCA cycle.

The rTCA cycle encoded by two carbon-fixing Aquificaceae MAGs affiliated with Sulfurihydrogenibium azorense (family Hydrogenothermaceae) and Aquificaceae family (most likely Thermocrinis or Hydrogenobacter; Supplemental Data S3 and S5) encode the key genes of two alternative versions of the rTCA-cycle (Hügler et al. 2007) that contain either citryl-CoA lyase/ citryl-CoA synthetase (CCL/CCS) or the ATP citrate lyase (aclAB) (Supplemental Data S5). The MAG affiliated with S. azorense encodes ATP citrate lyase (aclAB), which is responsible for converting citrate into acetyl-CoA and oxaloacetate using ATP and CoA. This MAG also encodes the citryl-CoA lyase (CCL), but we could not identify a homologue of citryl-CoA synthetase (CCS). CCS is required for the alternative version of the rTCA cycle (Hügler et al. 2007). Hence, because the CCS is not present, the MAG associated with S. azorense may use ATP citrate lyase (aclAB) in its rTCA cycle.

In contrast, the second carbon-fixing MAG affiliated with Aquificaceae (Fig. 3) encodes an alternative rTCA cycle exhibited by featuring both CCS and CCL enzymes (Supplemental Data S5). This MAG also exhibited an increased total amount of carbon fixation with increasing bicarbonate concentrations (R² = 0.94, Fig. 3), contributing at least 30 times more DIC assimilation at 10 mM bicarbonate compared to the S. azorense MAG that encoded the canonical rTCA cycle (Supplemental Data S4). The increased carbon fixation at higher bicarbonate concentrations of this particular Aquificaceae taxon in the alkaline hydrothermal fluids may be related to its use of the alternative CCS/CCL-bearing rTCA cycle, potentially in tandem with other as-of-yet unknown metabolic features. We can only speculate as to the physiological explanation for CCS/CCL-bearing rTCA cycle encoding Aquificae exhibiting higher carbon fixation at increased bicarbonate concentrations. However, a positive effect of DIC concentration on anaerobic carbon fixation has been observed recently at low pH (high CO₂ partial pressures) by the roTCA cycle in pure cultures of sulfate reducing bacteria (Steffens et al. 2021). Our findings indicate that higher bicarbonate concentrations at alkaline pH may also result in increased carbon fixation by particular rTCA encoding chemolithoautotrophic bacteria.

The quantitative measurements shown here, of carbon fixation by chemolithoautotrophic bacteria in a natural microbial community at multiple bicarbonate concentrations in a deep subsurface hydrothermal environment support the models of subsurface carbon cycling that have proposed DIC concentrations influence chemoautotrophic production (Barry et al. 2019, Fullerton et al. 2021). Our findings show that the stimulated carbon fixation by rTCA encoding primary producers resulted in organic carbon that moved relatively quickly through the microbial food web and was taken up by heterotrophic bacteria for secondary production. The microbial carbon fixation, therefore, supported the greater subsurface microbial community, and should act as a biological sink for deep inorganic carbon by converting it to particulate and DOC. Our findings may be relevant for global warming mitigation strategies (Ali et al. 2022) that plan to pump atmospheric CO₂ into the deep subsurface for potential carbon sequestration.

Acknowledgments

We are grateful for the support of the management and the staff of Kazdagi Thermal Resort & Spa. We also thank to Mustafa Coskun, Halil Gültekin, and Safa Erden, who participated a part of field campaign in 2019 and 2020.

Author contributions

Ömer K. Coskun (Data curation, Formal analysis, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing), Gonzalo V. Gomez-Saez (Resources, Supervision, Writing – original draft, Writing – review & editing), Murat Beren (Methodology), Doğacan Özcan (Methodology), Suna D. Günay (Investigation, Methodology), Viktor Elkin (Methodology), Hakan Hoşgörmez (Methodology), Florian Einsiedl (Methodology, Writing – review & editing), Wolfgang Eisenreich (Methodology), and William D. Orsi (Data curation, Methodology, Resources, Supervision, Writing – original draft, Writing – review & editing)

Conflict of interest

We declare no conflict of interest.

Funding

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 364653263—TRR to W.D.O and W.E., under Germany’s Excellence Strategy—EXC 2077-390741603, and by the DFG—Project-ID GO 3267/2-1—to G.V.G.-S.

Data Availability

The 16S rRNA gene and metagenomic sequence data were entered in the NCBI Short Read Archive under BioProject ID PRJNA837050. Metagenomic dataset and intermediate files to produce qSIP results were deposited under https://figshare.com/authors/_mer_Coskun/9725927. All remaining data needed to evaluate this article are available in the main text and/or the supplementary materials.