Subsurface rock core samples were obtained from different depths (1,679–2,908 mbs) of the Koyna pilot borehole (17°17″57.27” N, 73°44″19.07″ E) of the Deccan Traps, India. Samples were collected following standard techniques and stored anaerobically in sterile containers at 0–4°C and aseptically transported to the laboratory. The detailed processes of handling and storage procedures of the samples were mentioned in the previous study from our group (Sahu et al., 2022). For microbiological investigations, rocks were powdered in N₂-filled condition, using electrical coring machine, and the interior rock pieces/powder were obtained (Sahu et al., 2022). Details of the six rock samples (designated as C1, C2, C3, C4, C6 and C7), have been furnished in Supplementary Table 1, and were used for further enrichment studies.

Reactivation of indigenous microorganisms from the six rock core samples was performed by setting up microcosm-based enrichment setups. The rock powder was checked for the presence of contamination by estimating the sodium fluorescein concentration (used as a marker of potential contamination of drilling fluid) as per the protocol described earlier (Dutta et al., 2018). Rock powders devoid of any sodium fluorescein were used for further analysis. 0.05 g of rock powder samples were incubated anaerobically (in crimp sealed serum vials) at high temperature (50°C) using a basal medium [formulated based on the overall geochemistry of our rock samples as reported in (Sahu et al., 2022) and named as the Koyna medium] and other supplements. A Don Whitley anaerobic workstation (Model A-35) operating with N₂/H₂ + CO₂ + N₂ (10:10:80) gas environment was used to setup the enrichments. The detailed composition of the medium and incubation conditions, etc. have been provided in Supplementary Table 2. All medium components were prepared separately, sterilized either by autoclaving or filtration, (as appropriate) and used. The pH of the medium was set at pH 9.0. For the present study, each of the six samples (C1, C2, C3, C4, C6 and C7) were incubated with CO₂ (along with H₂ as electron donor) or 2 mM HCO₃⁻ (along with NH₄ as electron donor) as sole C source. CO₂ was added from a cylinder containing high purity gas mixture of CO₂ and H₂ in a premixed concentration ratio of 80:20. Ten ml of the gas mixture was added using the liquid displacement technique. A total of 12 experimental setups [marked as C1–C7 HC (for H₂ + CO₂ enriched sets) and C1–C7 BC (for HCO₃⁻ enriched sets)] and two controls (setups devoid of rock samples) were incubated for 180 days at 50°C and considered for further study post-incubation.

Total protein concentration of the enrichment cultures was measured to estimate the microbial biomass content. Protein concentration was measured through a fluorometric method (Qubit, Thermo-Fisher Scientific), using the Qubit Protein Assay Kit. In brief, enrichment samples were sonicated with 30% amplitude and pulse of 15 s on and 15 s off, for 10 min in order to obtain the total protein from the enrichment cultures. After sonication, the solution was centrifuged at 10,000 rpm for 5 min and the supernatant was collected and used as test samples for the estimation. The working solution as well as assay tubes were prepared and the estimation was done in triplicate using the manufacturers protocol (https://research.fredhutch.org/content/dam/stripe/hahn/methods/mol_biol/Qubit_Protein_Assay_QR.pdf). The fluorometer was calibrated against the standards and the test samples were measured in μg ml⁻¹.

A Luminometer (Glomax 20/20 System, Promega with BacTiter-Glo™ Microbial Cell Viability Assay kit) was used for estimating the total ATP of the living cells present in the enrichment samples (Adhikari and Kallmeyer, 2010). The details of the protocol were same as provided by the manufacturer, i.e the BacTiter-Glo™ Microbial Cell Viability Assay Technical Bulletin #TB337. Light emission was measured (as relative light units per second, RLU s⁻¹) for three, five seconds intervals with a five seconds delay before each interval, which was then quantified as mM with reference to standard ATP curve. Sterile tips were used in transferring all solutions and samples to exclude ATP contamination of pipettes and solutions. All procedures were performed in triplicates in the dark and all plastic material, solutions, and pipettes were stored in the dark as light could hamper the results causing delayed fluorescence.

The utilization/reduction of electron acceptors might act as a significant sign of growth of microorganisms. In this study, the initial and final concentrations of terminal electron acceptors like SO₄²⁻ and NO₃⁻, were measured to obtain information about the change in their concentration. Required amounts of enrichment cultures were obtained, centrifuged at 10,000 rpm for 30 min and the supernatant was used for further analysis. SO₄²⁻ estimation was performed through the BaCl₂ based turbidimetric spectroscopy method (Chesnin and Yien, 1951). For measurement of NO₃⁻ concentration, the samples were treated with salicylic acid-H₂SO₄ reagent, followed by NaOH, and estimated via colorimetric procedures (Cataldo et al., 1975). All analyses were performed in triplicates.

Formation of different metabolites as intermediates or end-products of various pathways could provide an indication of successful utilization of the provided enrichment substrates. Metabolites such as lactate, acetate, and formate in the enrichment setups were measured through Ion Chromatography (Eco IC, Metrohm) with PRP-X300 (Hamilton) column following the manufacturer’s protocol. The standard solutions were prepared using highly purified grade sodium salts. 15 ml culture was aseptically and anaerobically withdrawn from each of the serum vials, centrifuged (10,000 rpm for 30 min) and the supernatant was used for IC analysis. The samples were acidified using 0.1 N HCl, filtered using 0.22 μm filter (Merck Milipore, Millex GP) and transferred to the autosampler for further processing.

The entire culture volume of the serum vials (~100 ml) from each enrichment was centrifuged at 10,000 rpm for 30 min and the pellet (containing the rock powder and microorganisms) was used for extraction of total DNA in multiple replicates (n = 4) using QIAGEN power soil DNA isolation kit with certain modifications (30 min incubation after adding C2 and C3 solutions each). The total volume of eluted DNA from each sample was pooled, concentrated and quantified through a fluorometric method (Qubit, Thermo-Fisher Scientific). The V4 region of 16S rRNA gene was amplified from the isolated total DNA with V4-specific barcoded primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACVSGGGTATCTAAT-3′, Bates et al., 2011). The V4 amplicon libraries were purified, processed and sequenced through an in-house Ion S5 (Thermo-Fisher Scientific). Details of the sequencing procedures are described elsewhere (Dutta et al., 2018). Total DNA was extracted from sequencing reagent and sequenced using the same procedure, as control. The sequence reads were submitted to the Short Read Archive (NCBI) under BioProject ID PRJNA642293 (for H₂ + CO₂ set) and PRJNA643233 (for HCO₃⁻ set).

Bacterial abundance in the enrichment cultures was quantified by estimating the copy number of bacterial specific 16S rRNA genes (V3 region). All the amplification reactions were set up in triplicate. Quantitative PCR was performed in QuantStudio 5 using Power SYBR green PCR master mix (Invitrogen), with primer concentration of five pM and the following amplification conditions: 95°C for 10 min, 40 cycles of 95°C for 15 s, 55°C for 30 s and 72°C for 30 s. Melting curve analysis was run after each assay to check PCR specificity. The primers used were P1 (5’-CCTACGGGAGGCAGCAG-3′) and P2 (5’-ATTACCGCGGCTGCTGG-3′). Other details are same as described earlier (Dutta et al., 2018). Samples were run in at least three dilutions to check for PCR inhibitions. Bacterial 16S rRNA gene copy numbers were determined in each sample by comparing the amplification result to a standard dilution series ranging from 10² to 10¹⁰ of plasmid DNA containing the V3 region of the 16S rRNA gene of Achromobacter sp. MTCC 12117. We could not perform the quantification (as we did not obtain any amplification) of archaeal populations, because of lower DNA concentrations.

Raw sequence reads obtained after sequencing were quality filtered (removal of primers, homopolymers run of less than six bp, read length beyond the range of 250–300 bp, mean quality score below minimum of 25, with minimum six primer mismatches) and analysed in QIIME 1.9.1 pipeline (Caporaso et al., 2010). Denovo-based clustering (97% sequence similarity) of filtered reads to OTUs was performed using UCLUST under QIIME workflow. Taxonomic assignment of representative reads from each OTU was performed using SILVA 132 database.¹ The taxonomic assignments could not be performed in newer versions of SILVA database due to unavailability of required computational facilities. OTUs detected in reagent control, singleton OTUs (the cumulative read for which is one amongst all sets, Reitmeier et al., 2021) and OTUs affiliated to potential contaminants of deep subsurface, i.e., members of taxa Staphylococcaceae, Streptococcaceae, Propionibacteriaceae, Lactobacillaceae, E. coli-Shigella, Corynebacterium (Hubalek et al., 2016; Sheik et al., 2018), were removed from the OTU table before further analysis. Samples were normalized to 14,800 reads by random sampling for alpha diversity analysis.

Principal Component Analysis (PCA) and Canonical Correspondence Analysis (CCA) were performed using PAST software version 4.06 (Hammer et al., 2001) on the basis of (i) growth parameters and (ii) concentration of selected geochemical parameters along with abundant taxa (top 50% relative abundance) in both enrichment sets, respectively. Non-metric multidimensional scaling analysis (NMDS, Dixon, 2003; Doi and Okamura, 2011) was performed for comparing the similarity/dissimilarity of the microbial communities at taxonomic level by using metaMDS function of ‘vegan’ package in R software 4.2.1 (Kuhnert et al., 2005; Van der Loo, 2012) which helped to depict the effect of different substrates and depth on the microbial communities. Similarity percentage (SIMPER) analysis was done by using Simper function in ecological vegan package of R software 4.2.1 (R Development Core Team, 2017) to identify the microbial taxa responsible for the dissimilarity between different depths in both the substrate types (HC and BC). OTU overlap among the enrichments from different substrates and among different rocks in the same enrichment were elucidated using InteractiVenn (Heberle et al., 2015). Phyla belonging to individual enrichments were selected for correlation analysis and the correlation heatmap was constructed on the basis of Spearman correlation using METAGENassist (Arndt et al., 2012).

16S rRNA sequence reads of various OTUs were used to construct a phylogenetic tree using a distance matrix-based Neighbor-joining method in MEGA 11 software with 1,000 bootstrap re-iterations (Tamura et al., 2021). Initially, top 15 abundant OTUs from each enrichment set and the OTUs affiliated to Cyanobacteria from the enrichments were analysed by homology search for the closest sequence relatives (matches) using the National Centre for Biotechnology Information (NCBI) BLASTn nucleotide database. Multiple sequence alignments were performed with CLUSTALW package of MEGA 11. The sequences were represented with their accession number followed by the organism’s identity and its location. Finally, phylogenetic tree reconstruction and validation were performed using Neighbor-Joining (NJ) algorithm (Saitou and Nei, 1987) with 1,000 bootstrap re-iterations using Jukes-Cantor distance model (Jukes and Cantor, 1969).

Co-occurrence network analysis was performed at two different levels, considering (i) OTUs enriched in at least four samples out of total six at the lowest taxa level and (ii) association of taxa strongly correlated with depth of the rock samples and with relevant geochemical parameters of the rocks. Correlation among different taxa was calculated using otu.association command in mothur (Schloss et al., 2009). Correlation among different taxa and geochemical parameters was calculated using the PAST software version 4.06 (Hammer et al., 2001). The Spearman correlation values r ≥ ±0.9 and r ≥ ±0.3 were used for construction of co-occurrence network amongst taxa and amongst taxa and parameters, respectively. Visualization of the interaction network was done using Cytoscape 3.7.1 (Shannon et al., 2003).

Microbial metabolic pathways were estimated based on the 16S rRNA gene data from the closed OTU picking method using the PICRUSt v2 software (Douglas et al., 2020). PICRUSt v2 compares 16S rRNA marker gene data with the nearest match of the known genome sequence to allocate the functional attributes. A weighted nearest sequenced taxon index (NSTI) was calculated to evaluate the accuracy of prediction. The predicted KO numbers were plotted on KEGG pathway maps² separately. Abundances of genes related to general metabolic pathways, different C fixation pathways, lithotrophy related metabolisms, secondary metabolite production and xenobiotic degradation were selected for detailed analysis.

Microbial enrichment and response in anaerobic high temperature enrichment microcosms was assessed by estimating the cell biomass in terms of total protein, ATP, as a metabolic marker, utilization of SO₄²⁻ and NO₃⁻ (provided as terminal electron acceptor) and bacterial 16S rRNA gene copy numbers (Figures 1A–E). Substantial amounts of protein and ATP were obtained in both setups indicating a positive response of the rock communities to the substrate provided. Consumption of ~50% of the initial SO₄²⁻ concentration or ~ 94% of the initial NO₃⁻ concentration could be observed. The total amount of bacterial gene copy number in each enrichment culture as estimated by 16S rRNA gene targeted qPCR analysis indicated 2.9 × 10¹⁰ ml⁻¹ to 6.5 × 10¹⁰ ml⁻¹ copies in HC enrichment while in BC enrichments it ranged from 3.7 × 10¹⁰ ml⁻¹ to 5.8 × 10¹⁰ ml⁻¹. However, archaeal 16S rRNA genes could not be amplified during qPCR analysis. A PCA biplot on the basis of observed growth parameters reflected the distribution of the rock core enrichment samples in two halves of the plot (Figures 2A,B).

Identification of different metabolites (small organic acids), via ion chromatography was done. Approximately 77 mgl⁻¹ and 1.0 mgl⁻¹ formate as well as 191.1 mgl⁻¹ and 6.6 mgl⁻¹ lactate in HC and BC microcosms, respectively were estimated. However, ~9.5 mgl⁻¹ acetate was detected only in the BC sets (Supplementary Table 3).

Details of the sequence information have been summarized in Table 1. A total of 389,051 (HC) or 86,100 (BC) 16S rRNA reads was obtained. In HC setup, a maximum of 141,124 and a minimum of 18,136 quality filtered reads were obtained from C1HC and C4HC sets, respectively. Whereas in case of BC enrichments, a maximum of 23,686 and a minimum of 4,400 quality filtered reads were obtained from C7BC and C3BC sets, respectively. On an average more than 99% of the quality filtered reads could be classified as bacteria (368,807 in HC and 85,124 in BC sets) and archaea (16,989 in HC and 286 in BC sets). Following the removal of control and potential erroneous reads, a total of 18,258 OTUs (HC) and 10,482 OTUs (BC) were obtained (Supplementary Figure 1). The Venn diagram illustrated in Supplementary Figures 2A,B depicted the sharing of OTUs amongst the six samples in each enrichment. It was observed that the percentage of unique OTUs and their cumulative abundance decreased as the depth of the rock (used as sample for enrichment) increased. Only 34 and 9 OTUs were common amongst all six samples in HC and BC enrichments, respectively.

Details of the alpha diversity indices were summarized in Table 2. Chao1 estimator, which estimated the minimum total number of OTUs, ranged from ~4176 (C2HC) to ~7629 (C6HC) in HC sets and 770 (C3BC) to 7236 (C7BC) in the BC microcosms. The Shannon diversity index, which was calculated to evaluate species abundance and evenness in each sample, showed values between 8.4 (C2HC) to 10.1 (C4HC) in HC enrichments, whereas in case of BC enrichments, the values ranged between 7.8 (C3BC) to 10.3 (C7BC). The Shannon index values were slightly higher in the deeper horizon samples (C4–C7) under both enrichment conditions. Simpson index, an indicator of species dominance in a community, ranged between 0.98–0.99 across all enrichment samples.

The details of sequence affiliations in each taxonomic level were presented in Table 1. Assigned reads from HC enrichments (total number: 385796) were assigned to 576 genera belonging to 41 phyla (37 bacteria and 4 archaea). Compared to HC, number of taxa detected in BC setups were less. Assigned reads from the BC enrichments (85410) were affiliated to 435 genera belonging to 32 phyla (29 bacteria and 3 archaea).

Members of the phylum Proteobacteria constituted the most dominant group (upto ~63%) across the enrichment conditions (Supplementary Figure 3). The other abundant phyla were Actinobacteria (11%,), Firmicutes (8%), Bacteriodetes (8%), Thaumarchaeota (1.3%), Verrucomicrobia (1.13%), Chloroflexi (2%), Acidobacteria (1.2%) and Cyanobacteria (0.7%). Interestingly, it was observed that for both the enrichment conditions relative abundance of Gammaproteobacteria* (excluding Betaproteobacteriales) and Betaproteobacteriales members were higher in the deeper rock (C4–C7) based enrichments (28 and 53%) than in relatively shallower rock (C1–C3) enrichments (13.5 and 10%). However, the distribution of Actinobacteria, Firmicutes, Bacteroidetes, Cyanobacteria were found to be more in the upper horizon samples (C1–C3HC; C1–C3BC). Members of the phylum Chloroflexi were found to be one of the dominant groups with ~13% abundance in C3BC, whereas Acidobacteria members were found to be relatively abundant in sample C3, in both the HC and BC enrichments.

At the lowest taxonomic level, the most abundant microbial groups in HC enrichments were found to constitute >45% of the total community in each setup (Figure 3A). Comamonas, was found to be the major taxon in HC enrichments. Other major taxa included the members of Burkholderia-Caballeronia-Praburkholderia, Klebsiella, Ralstonia, Enterobacter, unclassified Enterobacteriaceae and Burkholderiaceae. Interestingly, these taxa were more abundant in the rock enrichments from the deeper horizons (C4–C7). On the other hand, mixotrophic Pseudomonas, Rhodanobacter, Methyloversatilis, members of Group 1.1c Thaumarchaeota and Conexibacter members were more abundant in enrichments of shallower rocks with HC (C1–C3). Noticeably, presence of metabolically versatile Oxyphotobacteria (belonging to Cyanobacteria) was observed in these enrichments with an average relative abundance of 0.7%, with higher abundance in the shallower samples (1.25%). In BC enrichments the top taxa contributed around ~45% in average (Figure 3B). The major groups belonged to similar taxa as in HC enrichment of the deeper horizon rocks. Comamonas, Burkholderia-Caballeronia-Paraburkholderia, Klebsiella, unclassified Enterobacteriaceae, Burkholderiaceae, Enterobacter and Ralstonia were the most dominant taxa with higher distribution in enrichments from deeper rocks. Nocardiodes, Sphingomonas, Aeromonas dominated the BC enrichments of shallower rock samples. Other less abundant taxa (relative abundance <0.1%) in HC enrichment microcosms showed the presence of important hydrogen oxidizing and/or CO₂ utilizing taxa like Chloroflexi, Curvibacter, Hydrogenophaga, Hydrogenispora, Methanobacterium, Methanolinea, Acidithiobacillus and Thiobacillus as well as desiccation and radiation resistant, endolithic, extremophilic Chroococcidiopsis. Syntrophic, chemolithoautotrophic, sulfur oxidizing groups belonging to Syntrophobacter, facultatively chemolithoautotrophic Thermincola, Pelobacter and autotrophic Curvibacter, aromatic compound degrader Azoarcus were observed as minor groups in BC enrichment setups. Supplementary Figures 4A,B illustrates the relative abundances of 50 representative less abundant taxa in HC and BC enrichments, respectively.

Rank abundance of observed OTUs was analysed to get an idea of the dominant reactivated community members under different enrichment conditions. In HC enrichments the top 50 OTUs cumulatively contributed to 30% abundance of the total. Rank abundance histograms and taxonomic distribution showed that an OTU, belonging to Burkholderia-Caballeronia-Paraburkholderia was the most abundant one with an average relative abundance of 4.8% (Figure 4A). It was followed by the OTUs related to denitrifying Rhodanobacter (2.1%), Thaumarchaeota (1.8%), Methyloversatilis (1.7%), Pseudarcicella (1.2%), Cloacibacterium (1.1%), taxa capable of hydrogen oxidation like Conexibacter (1%) and Pelomonas (1%). Interestingly, the presence of nine of these top 10 OTUs mentioned (except Conexibacter) were observed to have been enriched with bicarbonate, as enrichment substrate, as well. The other less abundant group contained OTUs belonging to potential H₂ and CO₂ utilizing autotrophic organisms like Ralstonia, Chloroflexi_KD4-96, Acidobacteraceae (Subgroup 1), Comamonas, Luteolibacter, Methylobacterium, Enterobacteraceae and facultative anaerobic Sporichthyaceae.

In BC enrichment microcosms, rank abundance analysis (Figure 4B) of top 50 OTUs, cumulative abundance 24%, revealed that the most abundant OTU to be the same (affiliated to facultative chemolithotrophic Burkholderia-Caballeronia-Paraburkholderia) as in the HC enrichment. OTUs related to Aeromonas (1.6%), Kurthia (1.1%), Nocardiodes (1.1%) and Burkholderia-Caballeronia-Paraburkholderia (1.01%) followed the dominant one. OTUs related to Aeromonas, Pseudomonas, Burkholderiaceae, Aquabacterium, Rhodanobacteraceae Enterobacter, Klebsiella, Nakamurella, Lautropia, constituted the less abundant OTUs. Supplementary Figures 5A,B illustrates the rank abundances of 50 representative less abundant (relative abundance <0.1%) OTUs in HC and BC enrichments, respectively.

The Non-metric Multidimensional scaling (NMDS) ordination was performed to understand the collective role of enrichment substrates and rock sample’s depths on the distribution of microbial taxa. The output reflected that deeper (C4–C7) rock hosted enriched communities occupied a distinct ordination space separated from shallower (C1–C3) rock enriched communities for each substrate (Figure 5). It was observed that microbial communities were clustered based on depth rather than the enrichment substrates (HC or BC). Stress was found to be ≤0.06 which indicated good fit of ordination (Dexter et al., 2018; Mullin et al., 2020). The partitioning was observed clearly in between the two halves of the plot. A canonical correspondence analysis (CCA) based on selected geochemical parameters along with enriched microbial taxa was drawn (Figures 6A,B), which revealed strong significant correlation of enriched microbial taxa with depth and parameters. Separate guilds were formed which were found to be guided by depth and geochemical parameters. Investigation of the individual enrichment revealed positive correlation of Mn with CCA axes 1 and 2 in HC enrichment sets (Supplementary Table 4A). Parameters like CO₂, H₂, and CH₄ correlated positively with CCA axis 1 and negatively with CCA axis 2. Positive correlation of axis 2 was observed with TIC, Fe, NO₃ ^- and SO₄²⁻, which on the other hand correlated negatively with axis 1. For BC enriched samples, CO₂, H₂, and CH₄ correlated positively with both the axes (Supplementary Table 4B). Mn correlated positively with CCA axis 1 and negatively with CCA axis 2. Negative correlation effect of TIC and SO₄²⁻ was observed for both the axes, whereas Fe and NO₃ ^- correlated positively with CCA axis 2. The distribution of the microorganisms in the deeper rock (C4–C7) enrichments were strongly constrained by Mn, CO₂, H₂, and CH₄ whereas the composition of the organisms in the shallower depth enrichment samples (C1–C3) could be guided by TIC, Fe, NO₃ ^- and SO₄²⁻ in both the enrichments. SIMPER analysis was done to determine the microbial taxa responsible for the differences between different depths in both enrichment conditions (Tables 3, 4). The topmost discriminating taxa were Comamonas, Burkholderia-Caballeronia-Paraburkholderia, Klebsiella, Enterobacter, Burkholderiaceae, Enterobacteriaceae and Ralstonia.

The Spearman correlation was used to identify the population at the phylum level that shared strong correlation amongst each other. In HC amended sets (Supplementary Figure 6), enriched members of the bacterial phyla Bacteroidetes, Spirochaetes, Cyanobacteria, Chloroflexi, WPS-2, and Patescibacteria as well as archaeal Euryarchaeota and Thaumarchaeota members constituted a strongly correlated clade. Armatimonadetes, Chlamydiae, Actinobacteria, Fusobacteria, Planctomycetes and Firmicutes populations formed the second guild. Coprothermobactereota, Altiarchaeota, Crenarchaeota, Caldiserica, Thermotogae, Verrucomicrobia constituted the third group. Lentisphaerae, Omnitraphechaeota, Synergistis and members belonging to FCPU426, Latescibacteria and Zixibacteria formed smaller mutually correlated clade. Proteobacteria, in spite of being the most abundant phylum, showed negative correlation with other enriched phyla. In the case of BC enrichments (Supplementary Figure 7) Chloroflexi, Gemmatimonadetes, Verrucomicrobia, Firmicutes, Actinobacteria and Bacteroidetes formed a strongly correlated group. The other two groups with strong correlation were formed by Cyanobacteria, Euryarchaeota, Deinicoccus Thermus and Crenarchaeota, Fibrobacteres, Tenericutes, respectively. Broadly, chemo/litho/autotrophic members of Cyanobacteria, Choloroflexi, Verrucomicrobia, Altiarchaeota were observed to correlate with members of Coprothermobactereota, Tenericutes and Gemmatimonadetes that are known to follow hetero/organotrophic lifestyle mostly with fermentative type of metabolism.

The co-occurrence pattern of the H₂ + CO₂ or HCO₃⁻ enriched microorganisms was analysed using network inferences based on strong and significant correlations, at the lowest taxonomic levels for OTUs present in at least four samples (Figures 7A,B). In HC enriched sets, the resulting microbial network consisted of 51 nodes and 187 edges (Figure 7A). The network consisted of two major hotspots (regions with highly connected taxa), with Bradyrhizobium, Brachybacterium, Streptomyces, Acinetobacter, Corynebacterium 1, being the taxa with the highest number of connections (16 each), followed by Sphingomonadaceae (15 connections). These taxa were found to be connected to Conexibacter, Nocardiodes, Pseudonocardia, Methylobacterium, Marinococcus, Stenotrophomonas, Actinobacillus, Curvibacter, Comamonas and Thermoplasmatales members. The second hotspot was formed by members belonging to Thaumarchaeota_Group 1.1c, Methyloversatilis, Rhodanobacter, Pseudomonas, Microtrichales, Smithella, Ralstonia and Cupriavidus. The most dominant taxa of the HC enrichments like Comamonas, Burkholderia-Caballeronia-Paraburkholderia, Enterobacteriaceae members (maximally dominant in C4–C7 enrichments) showed very few connections (mostly negative) with the other groups, hinting upon the fact that depth of the samples also played a role in determining the correlation amongst the enriched members. Members belonging to Burkholderiaceae, Chloroflexi, Alicycliphylus, Brachymonas, Citrobacter, Burkholderia-Caballeronia-Paraburkholderia, Comamonas, Cupriavidus, Ralstonia and Hydrogenophaga showed strong positive correlation with depth, and other geochemical parameters like CH₄, Mn, H₂, CO₂ as depicted in Supplementary Figure 8A. Bacterial taxa like Anoxybacillus, Conexibacter, Sideroxydans, Anaerococcus, Anaerobranca, Chloroflexi_KD4-96, Pseudomonas as well as archaeal taxa belonging to Methanomassillicoccalles and Methanosaeta were observed to correlate positively with Fe, SO₄²⁻ and NO₃⁻, and negatively with depth and CH₄, Mn, H₂ (Supplementary Figure 8B).

Co-occurrence network for BC enriched samples were formed of three interconnected hotspots comprising of 41 nodes and 56 edges (Figure 7B). Reyranella was found to be the taxa with highest degree connected to seven other taxa nodes including Exiguobacterium, Bacillus, Streptomyces, Actinomycetospora and Pseudomonas. Ralstonia was found to act as the nucleus node of the second hotspot connecting to six other taxa namely, Enterobacteriaceae, Burkholderiaceae, Klebsiella, Comamonas, Enterobacter and Limnohabitans. Members belonging to Thaumarchaeota, Aerococcus, Sphingomonas and Chitinophagaceae formed the third hotspot. Methanolinea, Chroococcidiopsis, Xanthobacter, Ralstonia, Cupriavidus, Curvibacter, Methylophilus, Syntrophaceae and Burkholderiaceae members were found to be correlated positively with depth and parameters including CO₂, H₂, CH₄ and Mn (Supplementary Figure 9A), whereas members of Thaumarchaeota_Group1.1c, Thermomicrobiales, Methylophilaceae, Thermoactinomyces, Thiomonas, Sporosarcina had strong positive correlation with SO₄²⁻, Fe and NO₃ ^- but was found to correlate negatively with depth, Mn, CH₄ and H₂ (Supplementary Figure 9B).

To determine the evolutionary history and gain insight about the relationships of the reactivated organisms, Neighbour Joining (NJ) phylogenetic analysis of 15 most abundant OTUs and cyanobacterial OTUs from both inorganic enrichments (HC and BC, making a total of 36 query OTUs) were performed by constructing a bootstrap tree (Supplementary Figure 10). The presence and relative abundance of these OTUs were depicted in the form of heat map beside each OTU. OTUs observed under both conditions were designated in green, the ones enriched only in BC sets were designated in blue and the OTU that belongs only to HC enrichment set was marked in red. OTUs enriched under both conditions, belonging to bacterial members like Burkholderia-Caballeronia-Paraburkholderia, Klebsiella, Rhodanobacter, Methyloversatilis, Pelomonas, Comamonas, Aeromonas, Kurthia, Nocardiodes, Corynebacterium 1, Dongia, Gemella, Deinococcus, etc. and archaeal group Thaumarchaeota were observed to be phylogenetically connected to sequences retrieved from deep terrestrial subsurface sediments from USA, subsurface clay rock, sediment from Eyreville drill corehole in Chesapeake bay, deep sea, deep hydrothermal vent in Nankai Trough, Hadal Deep Biosphere Mariana Trench and hydrothermal aquifer in Naica Mines. OTU affiliated to Conexibacter (unique to HC enrichment) was found to be phylogenetically similar (with 1,000 bootstrap) to an endolithic bacterial sequence from serpentinite-hosted ophiolite massif of the Alpine-Apennine chain in Italy. Sequences belonging to Fluviicola, Anaerolinea and Prolixibacteraceae were enriched only under BC enrichment condition and were found to be similar to sequences from deep subsurfaces of sedimentary rock Japan, Paris basin aquifer and Mariana Trench. Cyanobacterial sequences enriched from both sets showed strong matches with other cyanobacterial sequences retrieved from granitic subsurface environment in Tokyo, Clay rock borehole water of Switzerland, deep subsurface of South Africa and Mariana Trench.

The bacterial community analysis performed in this study suggested the presence of diverse bacterial populations with varied predicted metabolic functions in deep terrestrial subsurface. PICRUSt v2 based predictive metabolic profiling was done to gain insights into the metabolic potential of such microorganisms that help them to sustain in a subsurface environment. The relative abundances of genes related to different metabolic traits, that were predicted from our enrichments have been presented in Figure 8. Analysis revealed that the genes known to be involved in carbon cycling pathways, mainly in carbon fixation viz. acetyl-CoA decarbonylase (acsABCD), formylmethanofuran dehydrogenase (fmdABCDE; genes for CO₂ fixing WL pathway in acetogens and methanogens respectively), acetyl-CoA carboxylase, (accABCD; for HCO₃⁻ fixing 3-HP bicycle), phosphoenolpyruvate carboxylase (ppc), 4-hydroxybutyryl-CoA dehydratase, (abfD; for DC/4HB cycle, CO₂ and HCO₃⁻ fixing), ribulose-bisphosphate carboxylase (rbcS/L), phosphoribulokinase (prk; for CO₂ fixing CBB cycle) and Incomplete reductive TCA cycle, were predicted in all enrichment samples.

Lithotrophic metabolisms other than carbon dioxide and bicarbonate fixation play an important role in driving the energy cycle in the deep (Lovley and Chapelle, 1995; Osburn et al., 2014). The relative abundance of genes related to various lithotrophy related metabolism (as predicted) were presented in Figure 8, Energy metabolism.

Methane metabolism: Genes related to methane production (methyl-coenzyme M reductase; mcrA) and oxidation (methane monooxygenase; pmo, mmo) were predicted in most of the enrichment samples except in C7HC and C4BC. Overall relative abundance of genes mcr, pmo or mmo was found to be considerably lower than the abundance of other metabolic genes related to lithotrophy.

Nitrogen metabolism: The predicted presence and abundance of nitrogen metabolism genes in all enrichment sets hinted upon the significance of such metabolism. Genes related to nitrite oxidation (nitrate reductase / nitrite oxidoreductase; nxr), nitrate reduction (nitrate reductase; nar/nap, nitrite reductase; nir, nitric oxide reductase; nor and nitrous-oxide reductase; nos) and ammonia oxidation (ammonia monooxygenase; amo, hydroxylamine dehydrogenase; hao) were predicted in all enrichment sets.

Sulfur metabolism: Genes known to be involved in sulfide (anaerobic sulfite reductase; asr, cys) and thiosulfate (soxABCXYZ) oxidation, sulfate reduction (adenylylsulfate reductase; apr, dissimilatory sulfite reductase; dsr, sulfite reductase; cys), were also predicted from the enriched samples. However, the predicted relative abundance of dissimilatory sulfate reduction genes was lower than genes involved in assimilatory sulfate reduction.

Hydrogen oxidizing genes (hydrogenases; hnd, hyd, hox) were predicted to be more abundant in the deeper rock sample enrichments. Genes for production of acetate, lactate and pyruvate (ldh, pox, pfl, ack, ACSS) were predicted in all the enrichment samples.

Presence of genes related to synthesis of secondary metabolites like polyhydroxyalkanoate, polyhydroxybutyrate (pha/phb), antibiotic (vca, ntd, bac, dcs, pac, mbt, str), as well as degradation of chloroalkane (dha), chlorocyclohexane (lin), naphthalene (nah) and polycyclic aromatic hydrocarbon (nah, nid) could also be predicted from both the enrichment sets (Figure 8; Secondary metabolites and Xenobiotic degradation).