Exploration of deep terrestrial subsurface microbiome in Late Cretaceous Deccan traps and underlying Archean basement, India

Scientific deep drilling at Koyna, western India provides a unique opportunity to explore microbial life within deep biosphere hosted by ~65 Myr old Deccan basalt and Archaean granitic basement. Characteristic low organic carbon content, mafic/felsic nature but distinct trend in sulfate and nitrate concentrations demarcates the basaltic and granitic zones as distinct ecological habitats. Quantitative PCR indicates a depth independent distribution of microorganisms predominated by bacteria. Abundance of dsrB and mcrA genes are relatively higher (at least one order of magnitude) in basalt compared to granite. Bacterial communities are dominated by Alpha-, Beta-, Gammaproteobacteria, Actinobacteria and Firmicutes, whereas Euryarchaeota is the major archaeal group. Strong correlation among the abundance of autotrophic and heterotrophic taxa is noted. Bacteria known for nitrite, sulfur and hydrogen oxidation represent the autotrophs. Fermentative, nitrate/sulfate reducing and methane metabolising microorganisms represent the heterotrophs. Lack of shared operational taxonomic units and distinct clustering of major taxa indicate possible community isolation. Shotgun metagenomics corroborate that chemolithoautotrophic assimilation of carbon coupled with fermentation and anaerobic respiration drive this deep biosphere. This first report on the geomicrobiology of the subsurface of Deccan traps provides an unprecedented opportunity to understand microbial composition and function in the terrestrial, igneous rock-hosted, deep biosphere.

are typically low in organic matter, porosity and interconnectivity, and have been subjected to high temperature and/or pressure in the geologic past or at present 8,9 . Microbe mediated mineral precipitation may lead to loss of porosity, although secondary increase in porosity at greater depths via pressure dissolution of minerals may also occur 9 . Endolithic microorganisms in such extreme habitats typically occupy the fractures through available interconnections, with nutrients made available either from the rock itself and/or through transportation (via available interconnections) in the form of dissolved gases, solutes, or colloids 3,8,10 . While the lack of photosynthetically derived organic carbon led 'autotrophy to outstrip heterotrophy' 9 , geogenic supply of CO 2 and CH 4 (as carbon source), H 2 , CH 4 , NH 4 + and HS − (as electron donor) and Fe 3+ , Mn 4+ and CO 2 (as electron acceptors) have potential role in providing minimal nutrient and energy supply 11 . Hydrogen-driven subsurface lithoautotrophic microbial ecosystems (SLiME) wherein microbial populations can thrive on low concentrations of H 2 and CO 2 were initially considered to fuel the subsurface environment [12][13][14][15] . Although the dominance of such lithotrophic metabolism has been reported from several terrestrial subsurface including Columbia River basalt 12 , Fennoscandian shield 16 and Witwatersrand Basin 17 , the role of heterotrophic organisms, capable of utilizing the metabolic products of autotrophs in carbon assimilation and driving the overall community function is highlighted in later studies [18][19][20][21] . Overall, a combination of metabolic processes encompassing both autotrophy (including chemolithoautotrophy) and heterotrophy is implicated in deep terrestrial ecosystems 18,20,21 .
Endolithic microbial populations in deep igneous provinces have been investigated in basalt hosted oceanic crusts [22][23][24][25][26][27] , subsurface crustal environments with high temperature basalts 28 , subsurface gabbros and mantle type rock 29,30 and crustal fluids 31 . Most of these studies reveal the dominance of bacterial members affiliated to Proteobacteria (Alpha-, Gamma-and Delta-sub divisions), Actinobacteria, Firmicutes, Bacteroidetes and Chloroflexi known for sulfate reduction, nitrate reduction, fermentation, and chemolithotrophic metabolism 25,26,29,[32][33][34] . Members of archaea, though present in much lower abundance mostly belong to the Miscellaneous Crenarchaeota group (MCG), Methanosarcinales, and Methanobacteriales taxa. Compared to oceanic system, terrestrial subsurface crustal environments within igneous rocks (basalts and granites) are microbiologically less explored. A few additional studies have investigated the terrestrial subsurface igneous rocks directly 35,36 , while most such studies in the terrestrial realm remain restricted to groundwater samples recovered from basaltic/ granitic aquifers 4,[35][36][37][38] . In contrast to our understanding of structure and function of deep subsurface microbial community through all these studies, the microbiome of Deccan subsurface remains elusive. In this work, we try to demonstrate the nature and function of microbial communities within the deep terrestrial subsurface of Deccan traps using the rock cores recovered through scientific deep drilling 39,40 . Exploration of deep biosphere of granitic-basaltic crustal system is highly imperative to answer the fundamental questions on microbial life within deep, dark, oligotrophic crust, their mode of interaction and role in Earth's biogeochemical processes, origin of life on Earth and even on extra-planetary locations (e.g., Mars). Compared to a reasonable level of understanding of deep life within marine subsurface, geomicrobiology of subterranean igneous rocks at greater depths remains almost elusive.
Scientific deep drilling in the Koyna-Warna region of Deccan traps, western India provides a unique opportunity to investigate microbial life within deep terrestrial igneous rocks (basalts and granite-gneiss basement) 39,40 . The Deccan traps represents a continental flood basalt province covering an area of over 0.5 million km 2 of lava flows (~65 Ma) with a total thickness of over 2000 m near the eruptive centre in western India [41][42][43][44] . This thick pile of Deccan basalts rests on ~2.5 Ga Archean crystalline basement 45 . The Koyna-Warna region has gained further importance owing to the reservoir triggered seismicity ever since the impoundment of the Koyna dam in 1962 39 . Exploratory core drilling has been carried out at multiple sites in the Koyna seismogenic zone to ascertain the subsurface geology and structure in the region and carry out seismological investigations. The boreholes penetrated the entire thickness of Deccan traps (412-1251 m) and passed through a few hundred meters of the underlying granitic basement 39,46 . Study of drill cores also reveal that the transition of Deccan basalts to granitic basement rock occurs over a relatively short span of few tens of centimetres to few meters 47,51 . In these sections, rocks appear to be intermixed with both granitic and basaltic compositions.
In the present study, rock cores obtained from three exploratory boreholes in the Koyna-Warna region are used to decipher microbial community structure and function. Samples are collected from three horizons, i.e., Deccan basalt (BS) (~65 Ma), the underlying granitic basement rock (GR) (~2500 Ma), and the short section between the two formations representing the weathered surface of granitic basement affected by first lava flows (referred here as a transition zone, TZ). To the best of our knowledge, this is the first report describing subsurface microbiology of the late Cretaceous basaltic lava flows over Archean granitic basement.

Results
Geochemical characteristics of the samples. Thirteen rock core samples covering BS, TZ and GR horizons are collected from three bore holes. Core samples from different depths portray geochemical characteristics of the subsurface system ( Fig. 1, Supplementary Tables S1 and S2). Based on the tested parameters we find distinct and characteristic geochemical nature of the granitic and basaltic horizons. Samples from BS zone are mainly composed of feldspar, clinopyroxene (pyroxene mineral), augite (pyroxene mineral), anorthite (calcic plagioclase feldspar), andesine (plagioclase feldspar), hematite (iron bearing mineral), enstatite (pyroxene mineral), diopside (pyroxene mineral), magnetite (iron bearing mineral), magnesioferrite (iron bearing mineral) and albite (plagioclase feldspar). Samples from GR zone show the prevalence of quartz, albite, anorthite, microcline (alkali feldspar), anorthoclase (alkali feldspar) and andesine (plagioclase feldspar). Samples from TZ zone show a mixed nature. Samples from GR zone are more alkaline (pH range 7. 8 The distribution pattern of carbon, nitrogen and phosphorous has been analysed with respect to depth in order to assess the general nutritional condition of the subsurface crustal system. Total organic carbon content shows a strong negative correlation with depth followed by total carbon and total inorganic carbon. Nitrogen (NO 3 − and NO 2 − ) is generally low whereas phosphate is below detection limit except for samples obtained from the deeper horizon (PV8, P1 and U11).
Sequencing data. Sequence information and diversity indices are summarized in Table 1 . With respect to distribution of estimated bacterial and archaeal cells, higher relative abundance (average of 6.09%; SD = 0.91%) of archaeal cells are observed in samples from the transition zone followed by the basaltic horizon (4.81%; SD = 2.78%). Samples from granitic horizon exhibit on average 1.32% (SD = 1.38%) archaea of the total estimated number of cells. The observed abundance of archaea in samples from TZ and BS zone samples corroborates with higher level of total organic carbon and near constant level of total inorganic carbon in these samples compared to the GR horizon. Abundance of two major biomarker genes (dsrB involved in sulfate reduction and mcrA involved in methanogenesis) of sulfur and methane metabolism is studied by quantitative PCR across the different rock samples. Presence of dsrB gene is noted in all but one sample with copy numbers varying from 7.6 × 10 2 to 2.7 × 10 5 per gram of rock with an average of 7.7 × 10 4 (SD = 6.88 × 10 4 ). Gene fragment of mcrA is observed in 10 out of 13 samples and its copy number per gram of rock varies from 4.1 × 10 1 to 2.5 × 10 3 with an average of 7.7 × 10 2 (SD = 7.54 × 10 2 ). Interestingly, the overall abundances of both the genes are higher by at least one order of magnitude in BS samples when compared with those from TZ and GR.
Microbial community analysis. Sequences assigned to bacteria cover 54 phyla, 127 classes, 212 orders, 377 families, and 1049 genera whereas sequences assigned to archaea cover 12 phyla, 12 classes, 8 orders, 15 families, and 23 genera. On the basis of average relative abundance analysis at phylum level, Proteobacteria, Actinobacteria and Firmicutes constitute the three major populations of the deep biosphere of Deccan trap region with average relative abundance >11% (Fig. 2). Among those, Proteobacteria is the most dominant (20-77%) followed by Actinobacteria (8-73%) and Firmicutes (2-19%) (Fig. 2). Bacteroidetes, Chloroflexi, Cyanobacteria, Planctomycetes, Acidobacteria and Deinococcus-Thermus are the other phyla constituting the communities, with >1% average relative abundance across the samples, but often much higher levels in specific samples/ type of rocks. Euryarchaeota is the most dominant archaeal phylum detected across all the samples followed by Thaumarchaeota whose presence is observed in all but PV8 and U7 samples. On the basis of cumulative abundance, Gammaproteobacteria, Actinobacteria, Betaproteobacteria, Alphaproteobacteria and Bacilli are identified as the five major bacterial classes. Average relative abundance of Gammaproteobacteria in the BS cores is lower (23%; SD = 10%) than the GR (28%; SD = 12%) and TZ (30%; SD = 7%) cores. On the other hand, average abundance of Alpha-and Beta-Proteobacteria are higher in the BS cores (12% and 15% respectively; SD of 5% and 8% respectively) than that of GR (~7% for both; SD of 3% and 2% respectively) and TZ (9% and 7% respectively; SD of 5% and 1% respectively). Out of 12 assigned archaeal classes, Thermoplasmata exhibits highest average abundance followed by Methanobacteria, Methanomicrobia and Halobacteria.
Taxonomic affiliation at the lowest level (genus) is possible only for 89% of total classified reads (91% of the OTUs). Acinetobacter followed by Corynebacterium, Pseudomonas and Staphylococcus represent the most abundantly detected genera. Other bacterial genera which are detected in significant amount across all the samples are Enhydrobacter, Halomonas, Micromonospora, Rhizobium, Agromyces, Bacillus, Sphingomonas, Shewanella and Deinococcus. On the other hand, at genus level, Thermoplasma, Methanobacterium and Methanosaeta are found to be the most abundant archaeal genera. Similarity among the bacterial communities. Relationships between microbial community structures at different horizons (BS, TZ, GR) and environmental parameters have been established. Distance based redundancy analysis (db-RDA) separate the samples into two distinct abundance-weighted categories. Regardless of the borehole sites, subsurface rock microbiomes partition into (i) BS zone communities and (ii) TZ-GR zone communities (Fig. 4). Average taxonomic distribution patterns demonstrate that compared to both TZ and GR or only GR, BS microbiomes have relatively greater relative abundance of Alpha-and Betaproteobacteria and lesser abundance of Gammaproteobacteria, Actinobacteria, Chloroflexi and Euryarchaeota. GR microbiome in this reference show elevated dominance of Actinobacteria compared to its BS and TZ counterparts. SIMPER analysis is used to determine the microbial taxa (class level), responsible for the differences observed between the BS, TZ and GR microbiomes (Table 2). Additionally, this data also demonstrates the groups responsible for the overlap of transition zone and granitic horizon microbiomes. The top most 'discriminating' taxa are Gammaproteobacteria, Actinobacteria, Beta-and Alphaproteobacteria followed by Bacilli, Thermoplasmata, etc.
Endemic Deccan subsurface microbiome and OTU level analysis. Core OTUs, i.e., those present in all samples of a particular rock type from each of the three geologic horizons BS, TZ and GR are analyzed. Overall, 131, 337 and 143 OTUs constitute the core OTUs of BS, TZ and GR microbiomes respectively ( Supplementary Fig. S1). Among those OTUs, a considerable proportion (74%, 84% and 66% respectively) is unique to each geologic horizon and only 15 OTUs are found to be shared among all the 13 samples ( Supplementary Fig. S2). Taxonomic distribution of the core OTUs present a considerable difference among the three zones ( Supplementary Fig. S3). Compared to GR microbiome, BS exhibits relatively higher abundance of Actinobacteria (Corynebacteriales and Propionobacteriales), Alphaproteobacteria (Rhizobiales, Sphingobacteriales, Rhodobacteriales) and Betaproteobacteria (Burkholderiales), but lower abundance of Bacilli and Gammaproteobacteria (Pseudomonadales, Enterobacteriales, Oceanospirillales and others). TZ microbiome core community is more enriched. Actinobacteria, Alphaproteobacteria and Bacilli are mainly present with elevated abundance. However, lower abundance of Gammaproteobacteria, and unique presence of Clostridia, Bacteroidia,      Table S3). Metagenomic data is analysed to gain a first insight into the factors that allow sustenance of observed microbial communities, particularly the interactions among community members and their contribution in carbon, nitrogen and sulfur cycles in the nutrient-limited igneous subsurface of Deccan traps. Genes involved in carbon cycling pathways mainly in carbon fixation viz. rbcL, acsB, accC and ACLY, methane production (mcrA), methane utilization (pmoA) and different fermentative processes (ACSS, ldhA, poxB, pflD, ackA and others) have been detected ( Fig. 6a and Supplementary Fig. S5a). Key genes coding for CO 2 fixation by Calvin-Benson-Bassham cycle (rbcL), 3-hydroxypropionate cycle (accC) and reverse TCA cycle (ACLY) are observed in all the samples with varying abundance whereas the gene responsible for CO 2 fixation in Wood-Ljungdahl pathway (acsB) is observed in U7 (BS) and U11 (GR) only. Gene involved in methane production, mcrA is detected in all the samples whereas gene involved in methane utilization, pmoA is detected in BS and TZ zone samples only. Overall abundance of affiliated reads to mcr and pmo are considerably lower than the reads affiliated to marker genes of different carbon fixation pathways. 3-hydroxypropionate pathway is found to be the most dominant carbon fixation pathway in the Deccan subsurface (Fig. 6a). Presence of different marker genes taking part in various nitrogen cycling pathways is detected in all the three samples. Evidence of denitrification (nar, nap, nir and nor) is found. Presence of both cytoplasmic (narABDG) and periplasmic (napABC) nitrate reductases and all other enzymes that perform steps in the complete denitrification process suggests that dissimilatory nitrogen-transforming metabolisms are common ( Fig. 6b and Supplementary Fig. S5b). Furthermore, presence of nitrogen fixation gene (nifH) is observed in all the metagenomes whereas genes aiding in the process of nitrification (amoA and hao) are observed in U7 (BS) and U6 (TZ). Sulfur cycle is one of the important processes that take place in the subsurface environment where microorganisms are a significant contributor. Presence of genes taking part in dissimilatory sulfur reduction (dsr and apsAB), assimilatory sulfur reduction (cysC, cysH, cysJI and sir) and sulfur oxidation (sqr, sox and SUOX) are detected in the subsurface metagenomes ( Fig. 6c and Supplementary Fig. S5c). However, in Deccan subsurface, the number of reads affiliated to dissimilatory sulfur reduction genes is lower than genes involved in assimilatory sulfur reduction and sulfur oxidation (Fig. 6c). Presence of genes responsible for hydrogenase activity (hypABC-DEF) is also observed in low and varied abundances across the three samples.

Discussion
The  to minerals present, alkalinity, and low organic carbon along with distinct pattern of abundance for oxides. Nevertheless, variations as evident in a few samples could be attributed to geological phenomena. The variations are mostly in terms of relative abundance of oxides, in particular, higher weight percentage of Fe 2 O 3 and lower weight percentage of SiO 2 than the standard basalts and granites 53 . Further investigation on rock mineralogy by our laboratory (data not shown) indicates that apart from the characteristic minerals of igneous rocks (like pyroxenes, plagioclase feldspar and alkali feldspar) 54 , some samples contain metasomatic minerals (like berlinite, dravite, pyrochlore and phlogopite). Presence of these metasomatic minerals indicates secondary alterations and possible movement of fluids through different faults and fractures. These results are consistent with a previous study which depicts structural deformations and several instances of water movement through different fractures and faults associated with seismicity in the Koyna-Warna region 51 . Microbial life within the nearly impermeable, highly oligotrophic basalts and granitic bedrock is presumed to be constrained by prevailing physical and chemical factors. The igneous rocks in the deep subsurface have very low porosity (<0.2 ± 0.2% for granitoids) 50 and moderate-to-high temperature and pressure (25 °C increase per km in basalt and 15 °C increase per km in granite; 26.7 MPa increase in lithostatic pressure per km) 55 . In comparison to several other subsurface igneous provinces, the BS, TZ and GR samples of the present study show significantly low (over one order of magnitude) organic carbon (TOC) [56][57][58] . Deep subsurface igneous provinces are often devoid of photosynthetically derived organic carbon. Nevertheless, the geogenic carbon (in the form of dissolved inorganic carbon, organic acid, abiogenic CH 4 , C 2-4 hydrocarbons), electron donors (CO 2 , CH 4 , H 2 , reduced -S and -N) and electron acceptors (NO 3 − , NO 2 − , SO 4 2− , Fe 3+ , Mn 4+ , etc.) can act as essential metabolic resources and facilitate microbial life 59 . From the overall geochemical data, we could elucidate the nutritional status of this igneous realm which suggests its potential to support microorganisms that prefer chemolithotrophic/ chemolithoautotrophic lifestyle. The chemolithotrophic/chemolithoautotrophic mode could be linked to heterotrophic metabolism through metabolic products of lithotrophic organisms, thus driving the community's overall function.
Our qPCR results on 16S rRNA gene abundance indicate the presence of nearly 10 5 cells per gram of rock. Bacteria outnumber archaea in all samples. Estimation of cell numbers from qPCR based gene copies is done considering an average 16S gene copy per genome and this estimate may be prone to bias 27 . Additional caveat may be due to primer bias during 16S rRNA gene amplicon preparation. Nevertheless, the 16S rRNA gene amplicon data, qPCR based analysis and whole metagenome based data corroborate among each other with respect to microbial community composition. Cell number observed in our study is consistent with the "rocky deep biosphere trend" 9 . Current estimated cell number in various deep subsurface igneous provinces range between 10 4 and 10 8 cells per gram of rock 27,60-64 . Lack of correlation of cell number with depth within deep subsurface of Deccan traps is consistent with previous reports on several other terrestrial deep subsurface studies 33,34,65 . Relative consistency of cell numbers across depth in Deccan traps could be attributed to a common set of metabolic or ecophysiological constraints to the inhabiting cells. Considering the geologic age of formations and very low porosity (<0.2 ± 0.2%) of the rocks 50 , the origin of the observed subsurface basaltic and granitic crustal communities remains enigmatic. The only plausible explanation is dispersion of cells through water movement in various fractures and faults developed due to tectonic perturbations. A recent study suggests movement of water through fissures and fractures in the Deccan subsurface of Koyna seismic zone which is evident from secondary mineralization in basalts and basement granitoids 51 . 16S rRNA gene based community composition corroborates with our qPCR based estimation suggesting a more diverse bacterial population when compared to archaea. Considerably higher diversity of bacteria as observed in this study is in the same order of magnitude or even higher than those reported from other deep terrestrial habitats from a similar depth range 16,[35][36][37] . Sequencing depth attained in our work is significant leading to the detection of more OTUs and diversity. In addition to this methodological aspect (sequencing depth), previous studies have linked higher diversity in igneous rocks with the rock age and thus its alteration state. Younger basalts are reported to host less diverse communities while older basalts (>20000 years) generally display richer communities adapted to local alterations of the rocks 24 .
Phylum level composition of microbial communities of deep Deccan subsurface corroborates with recent surveys on microbial life within igneous rocks (basalts, granites and volcanic glasses) indicating predominance of Gamma-and Alpha/Betaproteobacteria, Actinobacteria, Chloroflexi, Firmicutes, Nitrospirae, Deltaporteobacteria, and Epsilonproteobacteria 22,38,61,[66][67][68][69] . The global occurrence of most of these taxa in almost all basaltic-granitic environments could possibly be considered as an indication for igneous rock hosted core group 27 . Members of Alpha-and Betaproteobacteria are relatively more abundant in BS horizon. Most of these proteobacterial taxa are known to be of diverse metabolic types especially H 2 /Fe 2+ /S oxidizing autotrophs; C1 compound-, hydrocarbonmetabolizing, microaerophilic, fermentative and strict anaerobic organotrophs 70 . Actinobacteria, though found in all three zones, is relatively more abundant in GR. Among the actinobacterial classes, presence of obligate acidophilic Acidimicrobiales members commonly found in Fe-rich ecosystem (including sites with volcanic rocks) and with Fe 2+ oxidizing or Fe 3+ reducing abilities justifies their niche specific colonization within the deep subsurface of Deccan traps 71 . A small proportion (average abundance ~2%; SD = 1.7%) of the communities is represented by Cyanobacteria. Despite that the 'rocks at this site have not seen sunlight' for several thousands to millions of years, presence of Cyanobacteria in the deep subsurface of Deccan traps environment corroborates well with several earlier deep biosphere studies, indicating their potential role in supporting sulfate reducing bacteria by providing dark fermentation derived H 2 72-75 . Cyanobacteria are considered to be ubiquitous in aquatic environments and it is possible that a cyanobacterial bloom at or near surface was trapped in the aquifer long time ago and/or the water might have intruded into the deeper layers of rock strata through faults/fractures. Whether these phototrophs survived because of their ability to switch to a fermentative lifestyle 76 or they are less prone to degradation cannot be determined with the data in hand and remains an open question to answer through further research.
Correlation among the microbial classes based on their abundance highlights the possible interplay among the groups. On the basis of known metabolic properties of these interrelated taxa, we hypothesize a possible metabolic roadmap followed by the Deccan deep subsurface microbiome which allows utilization of available resources and provides survival benefits. Close association of several taxa previously known from diverse deep, oligotrophic igneous subsurface and capable of autotrophic carbon fixation (Anaerolinea, Cyanobacteria, Epsilonproteobacteria), lithotrophic nutrition with oxidation of inorganic nitrogen (NO 2 − ), sulfur (S 0 ) and iron (Fe 2+ ) (Spirochitae, Nitrospira, Acidimicrobia) and organotrophic nutrition with fermentative, denitrifying and sulfate reducing activities (Firmicutes, Actinobacteria, Deltaproteobacteria) can be noted 21,[77][78][79][80][81][82][83][84] . Presence of methanogenic archaea (Methanomicrobia, Methanobacteria) along with Deltaproteobacteria is a clear indication of interspecies metabolite transfer through syntrophic activities. Close metabolic tie-up between autotrophic and heterotrophic members of the communities is previously reported from several oligotrophic deep subsurface ecosystems 20,21,75,77 . Although present in low abundance, correlation among Nitrospira, Spirochaetes, Chloroflexi, Cyanobacteria and Epsilonproteobacteria indicates that they could be the key players in driving the redox reactions and carbon flux through oxidation of inorganic substrates (e.g., reduced sulfur, nitrite, hydrogen) and dark fixation of carbon in the nutrient deprived, mineral rich, deep, dark Deccan subsurface [82][83][84] . Assimilated carbon and its fermentation products (e.g. H 2 , acetate, formate, etc.) drive other lithotrophic and heterotrophic metabolism including anaerobic nitrate/sulfate reduction and methanogenesis. Lower level of correlation among more abundant proteobacterial groups as evident from our study has been previously noted in the Fennoscandian shield, and could support the hypothesis that close interrelations among these taxa are possibly more relevant on or near surface rather than in deep terrestrial subsurface 16 .
It is interesting to note that a high percentage of OTUs are unique to BS, TZ and GR microbiomes (as sampled from various depths) and only a handful of OTUs (15/107721) covering 2.8% of the total sequences are shared among all the 13 samples. This possibly indicates community isolation between different horizons (lava flows in basalt, transition zone and granitic province). This particular observation presents a sharp contrast to the high degree of overlap in OTU distribution pattern in deeply buried oceanic crust and sedimentary habitat above 26 . Interestingly, the observed partitioning of BS and GR microbiomes (as observed in db-RDA) corroborates well with the previous reports of a basalt specific microbiome 85,86 . Based on these results, we propose that although taxonomically similar organisms have occupied different igneous rocks at varied depths, environment specific taxon recruitment based on local availability of electron donors, acceptors and carbon sources, and limited opportunities of migration led to the distinct taxonomic assemblages. Our db-RDA data clearly indicate such partitioning of interrelated guilds and these distinct microbial assemblages of GR and BS horizons are controlled by iron oxide, total organic carbon, nitrate, nitrite and sulfate. The BS microbiome guilds could potentially be regulated by higher nitrate but relatively lower nitrite and indiscriminately present (2 out of 6 samples) low sulfate levels. GR microbiome on the other hand is strongly affected by sulfate. Since we could not measure hydrogen and sulphide from the rock samples, their potential in influencing the microbial assemblages could not be ascertained in this analysis. Differential abundance of genes related to dissimilatory sulfate reduction (dsrB) and methanogenesis (mcrA) within BS and GR horizons possibly imply the disparity in recruitment of such processes within the subsurface provinces. These results suggest that microbial assemblages are affected by different geochemical parameters across the subsurface horizons of Deccan traps. Microbial co-occurrence analysis indicates the presence of distinct and well defined networks across the three horizons, which corroborates the correlation based observation. Families affiliated to the first hotspot have previously been reported to be involved in diverse, yet deep biosphere relevant metabolic processes, viz., heterotrophic, anaerobic/microaerophilic fermentation, nitrate and sulfate reduction [87][88][89][90][91][92] . Nitrate and sulfate reducers of this group may be fuelled by hydrogen generated by the fermentative members. Members of the second hotspot are mainly anaerobes/facultative anaerobes, syntrophs, autotrophs, nitrate reducers, hydrogen-and iron oxidizers known to withstand different environmental stresses 21,[93][94][95][96][97] . This hotspot represents microbial groups which are generally observed in extreme and oligotrophic environments and further highlight the link between autotrophic-fermentative/organotrophic-syntrophic organisms that could drive communities through anoxygenic carbon fixation, hydrogen metabolism, and nitrate/sulfate reduction. Nitrospiracea, Cyanobacteria and the other members of micro-network 1 play important roles in subsurface nitrogen and carbon cycle [98][99][100][101] , while Thermoplasmataceae and Terrestrial Miscellaneous Group archaea are known thermophiles, previously observed in hot springs and other igneous provinces 102 .
Microbial life in oligotrophic deep terrestrial subsurface hosted by igneous rocks are highly constrained by metabolic resources. Due to the extremely low cell turnover rates (hundreds to thousands of years), it is difficult to describe microbial processes by direct measurements 20 . Nevertheless, whole genome metagenome based studies, supported by 16S rRNA gene based analysis allow us to gain a better understanding. The metagenomic inventory supports our observation that there is a syntrophic relationship between autotrophic and heterotrophic organisms which could exchange fixed carbons, CO 2 , H 2 , etc. for their sustenance. Though there are presence of both denitrification and dissimilatory nitrate and sulfate reduction pathways for harvesting energy in the Deccan subsurface, nitrate appears to be the preferred energy source. Reduced sulfur, nitrite, iron and hydrogen could be potent electron donors to fuel the community metabolism. This inference is in line with previous reports on thermodynamics of microbial processes within deep biosphere as well as their metagenomic analyses 20,103 . The presence of different genes related to assimilatory as well as dissimilatory nitrogen metabolism indicates that there is likely a dynamic nitrogen cycling occurring within the Deccan subsurface. Overall abundance of nirB in all the samples (ammonia assimilation) is higher when compared to denitrification genes (nirK, nirS, norB and nosZ). Results from previous studies indicate that ATP synthesis in denitrification is far lower than that expected from the free energy changes and even lower than in nitrate ammonification 104 . This justifies the fact that ammonification could be more common pathway when compared to denitrification in the Deccan subsurface not only for generating more ATP but also the product of ammonification (NH 3 ) is more bioavailable than the product of denitrification, i.e., N 2 . Relative abundances of genes related to carbon cycle give a sense that CO 2 could be the driving force in subsurface igneous provinces of Deccan traps. 3-hydroxypropionate pathway (3-HP) is found to be the most abundant carbon fixation pathway. In addition to CO 2 fixation, this pathway enables synthesis of necessary metabolic intermediates and fatty acids from HCO 3 − 105 . Ample presence of HCO 3 − in the igneous rocks including the Deccan traps, which contains highest HCO 3 − concentration among other similar rock types and basaltic watersheds from around the world 106 , justifies the fact that 3-HP could be an important and common carbon fixation pathway. High abundance of HCO 3 − in the Deccan traps selectively allows both chemoautotrophic and chemoorganotrophic organisms harbouring accC to thrive in the Deccan subsurface. Similar result with ubiquity of accC gene is obtained in deep Fennoscandian bed rock fluid samples 18 .
In view of the geochemical environment, taxonomic composition and metagenomic inventory, we infer that autotrophic carbon fixation, sulfur/hydrogen oxidation, fermentation and nitrification-denitrification are likely to be the key reactions driving the subsurface ecosystem linking the carbon, nitrogen and sulfur cycles. Products of fermentation and other organic molecules (generated through accC activities) are used as potent substrates in the nutrient deprived subsurface provinces of Deccan traps. Low levels of organic molecules may be utilized by chemoorganotrophs for gaining energy or may be assimilated as a carbon source. CO 2 is a useful fuel for the autotrophs which further produce nutrients for other heterotrophic community members. The coupling of redox transformation of sulfur and nitrogen substrates to carbon flux is clearly revealed from metagenomic data. Previous reports suggest that though organotrophic denitrification is more spontaneous than lithotrophic denitrification in presence of reduced sulfur (on the basis of ΔG values), presence of acetate may accelerate the process of denitrification and sulfur oxidation simultaneously 107,108 . Higher abundance of acetyl CoA synthase and acetate kinase among other genes related to fermentation in all the three metagenomes gives us a sense that the presence of acetate is prevalent in the Deccan subsurface. Though metagenomic binning would yield better results and provide improved understanding of the coupled metabolic processes, from the present scenario it can be hypothesized that sulfur oxidation coupled with denitrification in presence of acetate may be a common process in Deccan subsurface, considering the subsurface anaerobic environment of Deccan traps. Similar results have been reported in a 1.34 km deep fault zone of Witwatersrand Basin, South Africa 17 . Overall metagenomic data are in accordance with the taxonomic outcomes where close association of chemolithoautotrophs and heterotrophs are evident.
In conclusion, this is the first study that provides a comprehensive overview of the microbial community structure and their probable ecological role in the subterranean province of seismically active Koyna-Warna region in the Deccan traps. Our results reveal that the microbial communities are predominated by Proteobacteria, Actinobacteria and Firmicutes. Although less abundant, members of Nitrospira, Spirochaetes, Chloroflexi, Epsilonproteobacteria and Cyanobacteria could play key roles in driving the community function. Partitioning of interrelated microbial guilds on the basis of rock geochemistry suggests environment specific taxon recruitment. Overall, 16S rRNA gene based community composition supported by metagenomic inventory indicate that chemoautotrophic assimilation of carbon with oxidation of sulfur, hydrogen and nitrite coupled with fermentation and anaerobic respiration using nitrate as preferred terminal electron acceptor fuel the community. Limited  (Fig. 1). Samples were checked for the presence of sodium fluorescein (500 mg m −3 ) used during drilling according to the protocol of Nyyssönen et al. 16 .
Only interior pieces of rock were selected for microbiological study to avoid potential drilling mud contamination as recommended elsewhere 109 . Samples were collected following aseptic techniques and stored in sterile containers at 4 °C for shipment. In laboratory, the samples were stored at −20 °C to limit the chance of microbial contamination, till further analysis.
Rock processing. Thirteen rock core samples, four each from Panchgani and Phansavale, and five from Ukhalu were collected and processed for further study. Rock cores were washed thoroughly under sterile condition with autoclaved, DNA free water (Thermo-Fisher Scientific) and sub-coring of the rock was done to get rock samples devoid of any possible contamination that might have occurred during drilling. At first, the uneven ends of the cores were removed using granite cutter fitted with sterilized cutting blade. Subsequently, with either a mechanical drilling instrument fitted with sterilized drill bit or with sterilized chisel and hammer, rock powder was obtained from the core interior part. The rock powder thus obtained was stored in sterilized DNA free containers at −20 °C for further geochemical and microbial analysis.
Geochemical analysis. Rock powder was sieved through 2 mm mesh. Acid digestion of the rock powders was carried out in a microwave digester (Milestone SK12) followed by the measurement of major elements by using ICP-MS (iCAP-Q, Thermo Scientific). Rock powders were ultra-sonicated in presence of water for major anions (Cl − , NO 2 − , SO 4 2− , NO 3 − ,and PO 4 3− ) analysis by using Dionex ICS-2100 (Thermo Scientific). Total organic carbon, total inorganic carbon and total carbon were measured on OI analytical TOC analyzer. Before analysing TOC, rock powders were fumigated with 1 N HCl for 48 hours under fume hood. Quantification of the elemental oxides was done using PANalytical Epsilon3 XRF instrument. All the measurements for major elements, anions, carbon content and oxides were carried out in triplicates and the standard error was <1% from the mean for all the analyses. Mineralogical analysis of the powdered rocks was conducted using PANalytical XRD X'pert Powder and on the basis of 2θ and d-spacing, mineral composition was evaluated using X'pert HighScore Plus software. Alkalinity of the rocks was measured using USEPA method 310.2. pH measurement was done by incubating 1 g of rock powder in 0.1 M CaCl 2 at 1:10 ratio (w/v), and highly sensitive probes (Orion) fitted with an Orion multimeter (Thermo electron corporation, Beverly, MA) were used 110 .
Extraction of metagenomic DNA. Metagenomic DNA was extracted from 13 samples using MoBio power soil DNA isolation Kit (MoBio) following the manufacturer protocol. Total DNA from reagent control was also extracted using the same procedure (as mentioned before) and were used subsequently to check any possible contamination. Quality of the extracted metagenomic DNA and its concentration was measured using NanoDrop 2000 spectrophotometer, followed by fluorometric quantitation using Qubit (Thermo-Fisher Scientific).
Quantitative polymerase chain reaction (qPCR). The bacterial and archaeal populations in a rock sample were quantified by estimating the copy number of bacterial and archaeal specific 16S rRNA gene. Copy numbers of two functional genes (dsrB and mcrA) were also quantified using qPCR-based technique to estimate the sulfate reducing and methanogenic populations. Details of qPCR primers used in the study are provided in Supplementary Table S4. 2 µl of metagenomic DNA was added to the PCR master-mix with a total volume of 10 µl. All the reactions were set up in triplicate. Quantitative PCR was performed in QuantStudio 5 with Power SYBR green PCR mastermix (Invitrogen), with primer concentration of 5 pM and the following amplification conditions: 95 °C for 10 minutes, 40 cycles of 95 °C for 15 sec, 55 °C for 30 sec and 72 °C for 30 sec. Melting curve analysis was run after each assay to check PCR specificity. Quantitative PCR was performed in triplicate and some of the samples having higher DNA yield were run in at least two 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 2 to 10 10 of plasmid DNA containing the 16S rRNA gene of Achromobacter sp. MTCC 12117. Genes encoding archaeal 16S rRNA, dsrB, and mcrA were PCR amplified from the metagenome, cloned in TA cloning vector and plasmid DNAs for each with copy numbers 10 2 to 10 10 were used as standards for quantitation purpose.
16S rRNA gene amplification and sequencing. Following the quality check, metagenomes were subjected to amplification of V4 region of 16S rRNA gene using primers 515 F/806R 111 (Supplementary Table S4). Each forward primer was tagged with sample specific 10-12 bp barcode for multiplexing during sequencing run. Duplicate PCR reactions were performed for each samples using AmpliTaq Gold 360 Mastermix (Thermo-Fisher Scientific) and the PCR products were gel purified using QIAquick Gel Extraction Kit (Qiagen). Concentrations of the purified PCR products were measured through Qubit and the purified amplicons from each sample were pooled with equimolar concentration and sequenced on Ion S5 platform (Thermo-Fisher Scientific). The sequence reads were submitted to short read archive under BioProject ID PRJNA389253. Quality filtering, operational taxonomic unit (OTU) picking and annotation. Raw data obtained from Ion S5 were preprocessed and quality filtered using the Quantitative Insights Into Microbial Ecology (QIIME 1.9.1) 112 bioinformatics pipeline in which sequences having lengths outside bounds of 200 and 310, mean quality score below minimum of 25, maximum homopolymer run exceeding limit of 6 and number of mismatches in primer exceeding limit of 6 were filtered out. After filtering and removing potential erroneous reads, OTUs were picked for further analysis. De novo based clustering of reads to form OTU was performed using UCLUST under QIIME workflow. Sequences having greater than 97% similarity were assigned to same OTU. Representative reads from each OTU were assigned taxonomy using UCLUST trained SILVA 128 database 113 . OTUs present in the reagent control were removed from the OTU pool of the samples. Alpha diversity parameters were calculated using alpha_diversity.py under QIIME workflow. For comparable alpha diversity analysis, the data set was normalized by random sampling of 115812 reads per sample (i.e., lowest number of filtered reads detected across all the samples) and control OTUs (i.e., OTUs present in reagent control) were removed and alpha diversity parameters were calculated as mentioned before. Average abundance of each phylum was calculated considering values obtained for samples from each horizon.
Statistical analysis. Principal component analysis (PCA) was performed using PAST software 114 on the basis of geochemical parameters. Geochemical parameters were normalized on the basis of feature scaling before PCA was plotted. Statistical significance of similarity of geochemical parameters across different horizons was evaluated on Euclidean matrix in PAST software using analysis of similarities (ANOSIM). Average abundances of each phylum across all the samples were calculated and compared. Standard deviations (SD) were calculated during average calculations and mentioned in brackets. Classes having average abundance of 0.1% across all the samples were selected for correlation analysis and correlation heatmap was constructed on the basis of Spearman correlation using METAGENassist 115 . Distance based redundancy analysis (db-RDA) of Bray-Curtis distance on the basis of all the microbial classes present across 11 samples, and 10 different physicochemical parameters was performed using Vegan package in R statistical software 116 . Similarity percentage (SIMPER) analysis was also done using PAST software for 11 samples to understand the microbial classes responsible for the grouping of different rock types. Two samples, PV2 and PV8 were not considered for db-RDA and SIMPER analysis due to over-representation of a single OTU in these samples. OTU overlap among the samples from similar rock type and among different rock types were elucidated using InteractiVenn 117 and Venny 2.1 118 .
Network analysis. Co-occurrence network analysis was performed considering the microbial families detected in the samples. Microbial families having greater than 0.1% average abundance across all the samples was used for the network construction. Correlation among different families was calculated using otu.association command in mothur 119 . Based on pairwise Spearman correlations with correlation values >0.75 and <−0.75 were used for construction of co-occurrence network using Cytoscape 3.4.0 120 .