# Microbial Ecology of Sulfur Biogeochemical Cycling at a Mesothermal Hot Spring Atop Northern Himalayas, India

Shekhar Nagar, Chandni Talwar, Mikael Motelica-Heino, Hans-Hermann Richnow, Mallikarjun Shakarad, Rup Lal, Ram Krishan Negi

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Shekhar Nagar, Chandni Talwar, Mikael Motelica-Heino, Hans-Hermann Richnow, Mallikarjun Shakarad, et al.. Microbial Ecology of Sulfur Biogeochemical Cycling at a Mesothermal Hot Spring Atop Northern Himalayas, India. Frontiers in Microbiology, 2022, 13, 10.3389/fmicb.2022.848010. insu-03650992

HAL Id: insu-03650992

https://insu.hal.science/insu-03650992

Submitted on 25 Apr 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

Check for
updates

# Microbial Ecology of Sulfur Biogeochemical Cycling at a Mesothermal Hot Spring Atop Northern Himalayas, India

Shekhar Nagar<fnref n="1" />, Chandni Talwar<fnref n="1" />, Mikael Motelica-Heino², Hans-Hermann Richnow³, Mallikarjun Shakarad⁴, Rup Lal⁵ and Ram Krishan Negi<fnref n="1" />*

## OPEN ACCESS

**Edited by:**

Vasvi Chaudhry,
University of Tübingen, Germany

**Reviewed by:**

Digvijay Verma,
Babasaheb Bhimrao Ambedkar
University, India

Puneet Singh Chauhan,
National Botanical Research Institute
(CSIR), India

**\*Correspondence:**

Ram Krishan Negi
rknegi@zoology.du.ac.in;
negigurukul@gmail.com

**Specialty section:**

This article was submitted to
Evolutionary and Genomic
Microbiology,
a section of the journal
Frontiers in Microbiology

Received: 03 January 2022

Accepted: 09 March 2022

Published: 13 April 2022

**Citation:**

Nagar S, Talwar C,
Motelica-Heino M, Richnow H-H,
Shakarad M, Lal R and Negi RK
(2022) Microbial Ecology of Sulfur
Biogeochemical Cycling
at a Mesothermal Hot Spring Atop
Northern Himalayas, India.
Front. Microbiol. 13:848010.
doi: 10.3389/fmicb.2022.848010

Sulfur related prokaryotes residing in hot spring present good opportunity for exploring the limitless possibilities of integral ecosystem processes. Metagenomic analysis further expands the phylogenetic breadth of these extraordinary sulfur (S) metabolizing microorganisms as well as their complex metabolic networks and syntrophic interactions in environmental biosystems. Through this study, we explored and expanded the microbial genetic repertoire with focus on S cycling genes through metagenomic analysis of S contaminated hot spring, located at the Northern Himalayas. The analysis revealed rich diversity of microbial consortia with established roles in S cycling such as *Pseudomonas*, *Thioalkalivibrio*, *Desulfovibrio*, and *Desulfobulbaceae* (*Proteobacteria*). The major gene families inferred to be abundant across microbial mat, sediment, and water were assigned to *Proteobacteria* as reflected from the reads per kilobase (RPKs) categorized into translation and ribosomal structure and biogenesis. An analysis of sequence similarity showed conserved pattern of both *dsrAB* genes (n = 178) retrieved from all metagenomes while other S disproportionation proteins were diverged due to different structural and chemical substrates. The diversity of S oxidizing bacteria (SOB) and sulfate reducing bacteria (SRB) with conserved (r)dsrAB suggests for it to be an important adaptation for microbial fitness at this site. Here, (i) the oxidative and reductive *dsr* evolutionary time-scale phylogeny proved that the earliest (but not the first) *dsrAB* proteins belong to anaerobic *Thiobacillus* with other (rdsr) oxidizers, also we confirm that (ii) SRBs belongs to δ-*Proteobacteria* occurring independent lateral gene transfer (LGT) of *dsr* genes to different and few novel lineages. Further, the structural prediction of unassigned DsrAB proteins confirmed their relatedness with species of *Desulfovibrio* (TM score = 0.86, 0.98, 0.96) and *Archaeoglobus fulgidus* (TM score = 0.97, 0.98). We proposed that the genetic repertoire might provide the basis of studying time-scale evolution and horizontal gene transfer of these genes in biogeochemical S cycling.

**Keywords:** metagenomics, hot spring, biogeochemical cycle, sulfur spring, evolution

<footnote n="1">¹Fish Molecular Biology Laboratory, Department of Zoology, University of Delhi, New Delhi, India, ²UMR 7327, Centre National de la Recherche Scientifique, Institut des Sciences de la Terre D'Orleans (ISTO), Université d'Orleans-Brgm, Orleans, France, ³Department of Isotope Biogeochemistry, Helmholtz Centre for Environmental Research – UFZ, Leipzig, Germany, ⁴Evolutionary Biology Laboratory, Department of Zoology, University of Delhi, New Delhi, India, ⁵NASI Senior Scientist Platinum Jubilee Fellow, The Energy and Resources Institute, New Delhi, India</footnote>
# INTRODUCTION

The untapped sulfur (S) compounds oxidizing microorganisms (SOM) and S compounds reducing microorganisms (SRM) microbial communities residing in extreme and contaminated environmental conditions such as hot water, sulfide contaminated springs offer an intriguing opportunity to explore the unique microbial diversity with uncovered metabolic potential (Ghilamicael et al., 2017). The investigations of such microbiota began with focus on identifying and culturing novel thermostable biocatalysts with huge biotechnological applications (Inskeep et al., 2010; Li et al., 2014; Ayangbenro and Babalola, 2017). However, little progress has been made in exploring the correlation between microbiome and geochemistry of hot spring systems particularly that possess mesothermic hot waters with neutral pH and elemental S and sulfate richness (Ghosh et al., 2012; Roy et al., 2020). Moreover, the survival of microbes in these niches is often supported by community dynamics and interactions. Studies of such ecosystems may provide insights into the microbial evolution of specific pathways for microbial biogeochemical cycling of minerals. However, with about 400 thermal hot water springs located in India, less than 15% have been explored for biogeochemical and taxonomical classification using genomics and metagenomics approaches (Cinti et al., 2009; Saxena et al., 2017). Sulfur springs provide harsh physiochemical conditions to sustain the growth of only meso- and hyper-thermophilic microbes which includes S oxidizers and sulfate reducers (Chan et al., 2015; Gonsior et al., 2018). The survival could also be achieved with "microorganism adaptation" by several resistance mechanism such as activity of biprecipitation, biosorption, extracellular sequestration, and/or chelation (Haferburg and Kothe, 2007). During these changes, the exchange of genetic material by means of horizontal gene transfer (HGT) is prevalent and necessary for the adaptation of microbes through the acquisition of novel genes.

Khirganga, the mesothermal S spring in Northern Himalayas discharging waters rich in sulfate, chlorine, sodium, and magnesium ions has remained uncharted so far (Shirkot and Verma, 2015; Poddar and Das, 2018). High levels of sulfides in the environment accounts for the milky appearance of the hot spring water with white microbial mats predicted to be formed from sulfide reduction by the S-related prokaryotes (SRP) enriched at this site (Sharma et al., 2004; Dong et al., 2019). The microbial S disproportionation one of the oldest (about 3.5 billion years ago; Finster, 2011) biological processes on Earth producing sulfide, sulfite, and sulfate compounds establishes a complex network of pathways in the biogeochemical S cycle. Thus far, it is the very foremost metagenomic investigation of microbial communities in Khirganga (average atmospheric temperature 6.9 ± 0.3) focused on exploring the microbial biogeochemical S cycling with a complex of disproportionation of elemental S conforming intermediary compounds. The current study was carried out via microbial mats, sediments, and hot spring water samples in hot spring to decipher the stabilized and diversified genes involved in S cycle intermediary process in anoxygenic, photolithotrophic and chemolithotrophic S-oxidizing and reducing bacteria (Dahl and Truper, 1994; Hipp

et al., 1997). The work expands the genetic and evolutionary information for S cycling genes and evaluates the biodiversity and applications for screening of the novel thermostable enzymes from microorganisms. Further, understanding these adaptations vis-à-vis the physiological properties and metabolic processes in these springs could be monitored as the engineered SRP consortia could develop into an effective tool in optimizing degradation of sewage waste in industrial processes (Ayangbenro et al., 2018). Also, the sulfate reducing bacteria (SRB) implied to treat various environment contaminants including metals (Mothe et al., 2016; Zhang et al., 2016), metalloids (Battaglia-Brunet et al., 2012; Sahinkaya et al., 2015), various non-methane hydrocarbons (Callaghan et al., 2012), alicyclic hydrocarbons (Jaekel et al., 2015), nitroaromatic compounds (Boopathy, 2014; Mulla et al., 2014), and aromatic hydrocarbons (Stasik et al., 2015; Meckenstock et al., 2016; Kamarisima et al., 2019).

# MATERIALS AND METHODS

## Sample Collection, Physicochemical Analysis, and Helium Ion Microscopy

Samples of microbial mat deposits (250 g), sediment (250 g) and water (5L) were collected from Khirganga hot water spring (31°59'34" N, 77°30'35" E) in February 2017. Sampling was performed in two replicates for each habitat from two closely located primary thermal outlets (31°99'18" N, 77°50'96" E) and secondary outlets (31°99'19" N, 77°50'96" E). The surface temperature and pH of each habitat were recorded on site.

First, microbial mats and sediment were digested in pure nitric acid and water samples were filtrated to 0.1 µm prior to chemical analysis. All samples were subjected to physicochemical analysis for major elements. Concentrations of major cations (Na+, K+, Mg2+, and Ca2+) and anions ($SO_4^{2-}$ and $Cl^-$) were analyzed by ionic chromatography (Dionex ICS-2000, Sunnyvale, CA) using the columns CS16A for measuring cations and AS17 for anions. An elemental analysis of minor and trace elements through inductively coupled plasma mass spectrometry (ICP-MS) Agilent ICP-MS 7,900 with ultra-high matrix introduction (UHMI). The samples of sediments and microbial mats were desiccated overnight followed by ethanolic dehydration and microstructure was studied using a scanning electron microscope at the Center for Chemical Microscopy (ProVIS). Images were captured using a high efficiency detector.

## Metagenomic DNA Extraction, Sequencing, and Assembly

For the extraction of total DNA from microbial mats, 0.25-g samples were processed following a method described by Varin et al. (2010). The total community DNA from 0.25-g sediment samples and 5 L of filtered water (0.45 µm) were extracted using PowerMax Soil DNA isolation kit (MoBio Laboratories Inc., Carlsbad, CA, United States) following the manufacturer's instructions. Sequencing was performed at Beijing Genome Institute (BGI), Hongkong, China using Illumina Hiseq 2,500 platform. Paired end libraries of read length 100 base pairs
(bp) were generated with insert size of 350 bp. The raw sequences were quality filtered using SolexaQA (Cox et al., 2010), and the low-quality sequences below Q20 quality cut-off and artificially duplicated reads (ARDs) were castoff using Illumina-Utils (Eren et al., 2013) and duplicate read inferred sequencing error estimation (DRISEE) (Gomez-Alvarez et al., 2009), respectively. Further, the assembly was integrated in IDBA-UD (Peng et al., 2012) with 50-bp insertion length, minimum k-mer: 31, maximum k-mer: 93 (61 for water) using seed k-mer size for alignment 30 bp and the minimum size of contig as 200 bp while allowing minimum multiplicity for filtering k-mer while building the graph.

## Taxonomic and Functional Assignments

Alpha diversity within each sample was estimated as abundance-weighted average of annotated species from source databases built in MG-RAST v3.0 (Meyer et al., 2008) and expressed as Shannon diversity index transformed based on rarefaction curve using the following formula:

$$ \log_{10} (\alpha - \text{diversity}) / \ln (10) $$

Diversity at phylum level was inferred from MG-RAST (maximum e-value, $1 \times 10^5$ and minimum percentage identity cutoff, 60%). A paired-sample t-test was applied on the phylum determined in any habitat pair to estimate significant similarities based on taxonomic mean abundance using SPSS (SPSS Inc., version 20.0, IBM). Microbial genera were deciphered based on clade-specific markers to identify taxonomy up to species level using MetaPhlAn v2.0 (Truong et al., 2015) and a heatmap was constructed using Bray-Curtis dissimilarity with supporting dendrograms for both species and samples. We used HUMAnN2 (Franzosa et al., 2018) to perform phylum-resolved functional profiling of the communities that maps contigs onto the pangenomes of the known species of the community and quantifies the pathways and UniRef90 gene families database (Suzek et al., 2015). Later, these UniRef90 families were regrouped as clusters of orthologous groups (COGs) annotations based on eggNOG (Huerta-Cepas et al., 2017). Open reading frames (ORFs) of assembled metagenome were predicted using Prodigal v2.6.1 (Hyatt et al., 2010) and annotated at hierarchy levels, namely, subsystems, protein families, and individual enzymes using Prokka v1.12 (Seemann, 2014). The amino acid sequences were mapped against Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2004) and top 50 metabolic pathways in all of the six samples were compared through heatmap constructed using package pheatmap (Kolde and Kolde, 2015) and ggplot2 (Wickham, 2009) in R (R Development Core Team, 2011). Identification of S substrates disproportionation genes was performed by mapping all predicted ORFs on the HMM databases obtained from TIGRfam v10 (Haft et al., 2012) and Pfam (Finn et al., 2014) using hmmscan v3.1b2 (Eddy, 2011). An abundance of each enzyme was plotted as number of copies annotated within each sample. The sequences with more than 150 amino acids were queried against the National Center for Biotechnology Information (NCBI) Microbial proteins from RefSeq nr database

(04 April 2020) using BLASTp (Altschul et al., 1990) to identify the sequences producing significant alignments for taxonomic confirmation.

## Analysis of Diversity of Sulfate Reduction Proteins

For sequence similarity networks (SSN), amino acid sequences of sulfide oxidation and sulfate reduction proteins annotated in all six samples annotated by KEGG Ids were implied over an empirical measurement of diversity. For this, an all-vs.-all BLAST was performed to define the similarities/variations between sequence pairs of diversifying sulfate reduction proteins. A user defined threshold was optimized according to the alignment score and maximum length of BLAST results in diversifying and stabilized protein sequences. Clustering was performed using CD-HIT (Li and Godzik, 2006) on the scores of BLASTp pairwise alignments at a threshold value (e-value of 1e-30). The networks were visualized in Cytoscape v3.7.1 (Shannon et al., 2003). The average number of degree and neighbors for a protein sequence or a node was calculated as:

$$ k = \frac{2k}{N} $$

where, K is denoted with number of edges and N is denoted with total number of nodes. Also, to determine the divergence/similarity among nodes or protein sequences was calculated as:

$$ D = \frac{2k}{N(N-1)} $$

The attributes of node degree distribution, average clustering coefficient, average neighborhood connectivity, and closeness centrality were studied through power law fits to determine their correlation with number of neighbors. Sulfur oxidizing bacteria (SOB) and SRB were identified for the sequences that could be classified up to genus level to study the distribution of S substrates oxidation and reduction genes in the different clusters.

## Sequence Alignment, Phylogeny and Structure Prediction of Putative Unidentified Dsr and Asr Enzymes

To elucidate the phylogeny of key sulfite reductases, the DsrA/B and AsrA/B protein sequences (more than 150 amino acids) were individually aligned using MUSCLE v3.8.31 (Edgar, 2004) and clustered using UPGMB (unweighted pair group method with arithmetic mean). All the alignments were end trimmed manually and maximum likelihood (ML) phylogeny was inferred with 500 bootstrap resampling using RAxML v8.0.26 (Stamatakis, 2014). For this, we used standalone version of RAxML which was called as follows:

$$ \begin{array}{l} \text{raxmlHPC} - \text{PTHREADS} - \text{s input} - \text{N 500} - \text{n result} - \text{f a} \\ \qquad - \text{p 12345} - \text{x 12345} - \text{m PROTGAMMAGTR.} \end{array} $$

The resulting phylogenies were also confirmed using most complex general time-reversible model (CAT-GTR; Tavaré, 1986)
with PhyloBayes v1.7b using CIPRES Science Gateway v 3.3 (Lartillot et al., 2009) that incorporates different rates for every change and different nucleotide frequencies.

For proteins showing similarity with those from uncultured bacteria, we determined the structures using I-TASSER suite (Yang and Zhang, 2015). These predicted structures were then aligned onto their top structural analogs and C-scores, TM-scores, and RMSD were computed and ligand binding sites with conserved residues were identified. The TM-score is to compare two models based on their given residue equivalency (i.e., based on the residue index in the PDB file). It is usually NOT applied to compare two proteins of different sequences. The TM-score predicted from structural alignment of two proteins while comparing them based on residue equivalency such that a score of 0.6 and above denote the two proteins to be fairly aligned (Yang and Zhang, 2015). The TM-align will first find the best equivalent residues of two proteins based on the structure similarity and then output a TM-score.

# RESULTS

## Description of Sampling Site and Microscopic Analysis of Samples

The Khirganga is a natural hot spring setting that lies in the Parvati Valley in the Northern hemisphere of the great Himalayas (31°59'34" N, 77°30'35" E, altitude 2,978 m) at district Kullu, Himachal Pradesh, India (Figures 1A,B). For this study, samples from all three habitats, namely, microbial mats, sediments, and water were collected proximal to the major opening (KgM1, KgS1, KgW1) and from a distance of 10 m (KgM2, KgS2, KgW2) as shown in Figure 1C.

Using an electron microscopy, we dissected the microstructure of the niches and were able to visualize cellular structures on complex sample matrices. The microbial diversity was visualized as numerous filamentous structures in microbial mats and sediments that resembled Cyanobacteria. In addition, rod and cocci shaped cells of varying sizes were also observed in the sample matrices that providing a visual insight into the microbial diversity at this mesothermic site (Figure 1D).

## Physicochemical and Elemental Analysis

The in-situ measures of water temperature were from 59°C at the outlet to 55°C at 10 m distance (Table 1). Microbial mat deposits and sediments had much lower temperature (42-45°C) than water. The pH of the hot spring water was 6.7 while sediments and mats were slightly acidic with pH 6.1 and 6.3, respectively. Thus, all three habitats were recorded to be mesothermic. The physicochemical composition of the hot spring is dominated by anions of chloride (up to 11,024 ug/g) and sulfate (up to 10,079 ug/g) while ions of calcium and potassium were abundant (Table 1 and Supplementary File 1). Importantly, sulfates ($SO_4^{2-}$) concentration in microbial mats and sediments were higher (9,529 ± 313.29 ug/g; 10,079 ± 863.29 ug/g, respectively) and exceeded the limit of 8,000 ug/g standardized by Environment Protection Act (EPA, 2001) and also found to be exceeded the limit of 53 mg/L in surface waters

(79.82 ± 1.85 mg/L) (EPA, 2001; Table 1). The chlorides (1,456.77 ± 367.27 mg/L), manganese, sodium, and silicon constituents in the hot spring waters were surpassing the normal average concentrations of 250, 0.05, 200, and 4 mg/L, respectively (EPA, 2001) in surface water samples (Supplementary File 1). Among others, the predominant elements and minerals in water samples were aluminum, magnesium, copper, zinc, and arsenic.

## Metagenomic DNA Sequencing and Assembly

A large metagenomic dataset was obtained from sequencing having number of reads sized up to ~18 Gb for each sample. We retrieved a total number of reads ranging between 1.1 × 10⁸ to 1.5 × 10⁸ in all samples which were assembled into 180,849-519,194 (more than 200 bp) contigs. After assembly, the metagenomes sizes varied between 329 and 600 Mbp. A summary of characteristics of the datasets and assembled metagenomes is provided in Table 2. The alpha diversity estimated as the Shannon diversity indices ranged between 2.5 and 3 (Table 2 and Supplementary File 2).

## Microbial Consortia and Proportionality of S Oxidizing Bacteria and Sulfate Reducing Bacteria

Bacteria belonging to 15 different phyla dominated the microbial communities. The average percentage relative abundances of major phyla in the three habitats shown in parentheses in the order microbial mat, sediment, and water is as follows: *Proteobacteria* (62.1, 50.5, and 58.7%), *Bacteroidetes*, *Firmicutes*, *Cyanobacteria*, *Planctomycetes*, and *Chloroflexi* (Figure 2A). Species belonging to phylum *Proteobacteria* are found in varied temperature ranges which results in their dominance in various hot springs (Bowen de León et al., 2013; Singh and Subudhi, 2016; Saxena et al., 2017) and disproportionation of S compounds is mainly carried out by SRM of *Proteobacteria* (Finster, 2011), Besides, *Actinobacteria*, *Spirochaetes*, *Verrucomicrobia*, *Acidobacteria*, *Deinococcus-Thermus*, *Deferribacteres*, *Chlorobi*, *Gemmatimonadetes*, and *Nitrospirae* were also detected in all three habitats with relative abundances less than 3%. Among the three habitats, microbial diversity profiles of water and sediments were more similar compared to those of microbial mats.

The highest genus level diversity was revealed in sediment (n = 196) followed by water (n = 132) and mat (n = 63). The top-50 genera in all habitats were plotted (Figure 2B). The microbial mats were dominated by *Pseudomonas* (51.8%) followed by unclassified genera of family *Desulfobulbaceae* (SRB; 5.2%) and *Flavobacterium* (4.3%). Among the abundant genera in the microbial mats were *Thioalkalivibrio* (1.3%, SOB), *Aeromonas* (1.3 %), *Klebsiella* (1%), *Exiguobacterium*, *Enterobacter*, and *Escherichia* (more than 1%) (Figure 2B). On the other hand, sediment habitats were found to be enriched in *Thioalkalivibrio* (SOB; 18.9%), *Desulfobulbaceae* (SRB; 9.9%), *Halothiobacillus* (8.3%), *Burkholderia* (7.1%), and unclassified genera of families *Acetobacteraceae* (6.1%). The hot spring waters with highest diversity of bacterial genera were dominated by *Thioalkalivibrio* (SOB; 20.5%), *Acetobacteracea* (9.6%), and *Desulfovibrio* (SRB; 4.3%).
<figure><img src="image_4.png" /><figcaption>**FIGURE 1** | (A) Geographical location of Himachal Pradesh, Northern Himalayas. (B) Khirganga hot spring as shown in the maps are located in Parvati valley in Kullu district of Himachal Pradesh, India. (C) The different habitats from where the samples were collected are shown: microbial mat (green), sediment (yellow), and water (blue). Samples were collected in replicates from two outlets located 10-m distance apart. (D) The scanning electron micrographs of different habitat samples shown with arrow demarcating the filamentous (cyanobacteria), cocci-shaped and rod-shaped bacteria (SRB) in pink, yellow, and white, respectively.</figcaption></figure>

5.9%). Other genera with less than 6% abundances in all three habitats were also detected as shown in **Figure 2B**.

## Metabolic Functions of the Community and S Disproportionation Genes

The major gene families inferred to be abundant across all three habitats were assigned to *Proteobacteria* followed by *Chloroflexi*, *Firmicutes*, *Bacteroidetes*, and *Spirochaetes* as reflected from the reads per kilobase (RPKs) in the metagenomes. These gene families were then regrouped as COGs and the top functions were determined to be translation, ribosomal structure and biogenesis (COG: J), amino acid transport and metabolism (E), general function prediction (R), energy production and conversion (C), replication, recombination and repair (L), and carbohydrate transport and metabolism (G). These functions in microbial mats were carried out by *Proteobacteria* (COGs: J, E, C, and G) and unclassified bacteria (COGs: R and L); in sediments by unclassified bacteria (COG: J) and *Proteobacteria* (COGs: E, R,

C, L, and G) and in water by *Proteobacteria* (COGs: I and C), *Firmicutes*, and *Chloroflexi* (COGs: E, R, G, and L) (**Figure 2C** and **Supplementary File 3**).

The ORFs that were categorized on the basis of KEGG categories were mapped onto the metabolic functions and the pathways that could be reconstructed with more than 60% completeness were used to define the metabolic potential of the habitats. Based on this criterion, we studied the top-50 functional pathways of each habitat and identified the core functions (n = 37 pathways) of the communities that included the common pathways for metabolism of nucleotides, carbohydrates, and amino acids (**Figure 3A** and **Supplementary File 4**). In addition, we determined differentially abundant pathways in each habitat: microbial mats (n = 7), sediments (n = 7), and water (n = 2). The community functional profiles of sediment and water were more similar compared to those of microbial mats which may be due to the stratified layered organization of the mats which are different in sediment and water. The microbial communities in mats were optimized for metabolism of methane specifically,
<table><caption>**TABLE 1** | Physicochemical and elemental analysis of the niche samples.</caption><thead><tr><th rowspan="2">Environmental data</th><th colspan="3">Value for indicated sampling sites</th></tr><tr><th>Mat</th><th>Sediment</th><th>Water</th></tr></thead><tbody><tr><td>Temperature (°C)</td><td>42.4 ± 0.56</td><td>44.55 ± 1.06</td><td>57 ± 2.2</td></tr><tr><td>pH</td><td>6.2 ± 0.14</td><td>6.5 ± 0.14</td><td>6.7 ± 0</td></tr><tr><th>Major ions/metals</th><th>(µg/g)</th><th>(µg/g)</th><th>(mg/L)</th></tr><tr><td>Cl<sup>-</sup></td><td>11024.88 ± 983.17</td><td>15547.08 ± 204.96</td><td>1456.77 ± 367.27</td></tr><tr><td>SO4<sup>2-</sup></td><td>9529.06 ± 313.29</td><td>10079 ± 863.29</td><td>79.82 ± 1.85</td></tr><tr><td>Na<sup>+</sup></td><td>976.88 ± 64.83</td><td>991.47 ± 7.88</td><td>1338.52 ± 167.97</td></tr><tr><td>Ca<sup>2+</sup></td><td>7529.71 ± 1034.68</td><td>5111.4 ± 588.21</td><td>30.27 ± 7.32</td></tr><tr><td>K<sup>+</sup></td><td>3681.66 ± 886.65</td><td>7743.91 ± 412.12</td><td>41.56 ± 9.03</td></tr><tr><td>Al</td><td>7205.41 ± 1685.37</td><td>22011.1 ± 2312.9</td><td>0.05 ± 0</td></tr><tr><td>Mg</td><td>3528.62 ± 735.04</td><td>8204.28 ± 964.02</td><td>8.83 ± 0.02</td></tr><tr><td>Si</td><td>3.22 ± 0.66</td><td>5.28 ± 0.42</td><td>31.5 ± 1.36</td></tr><tr><td>Mn</td><td>1149.7 ± 178.85</td><td>1219.29 ± 13.93</td><td>0.34 ± 0.03</td></tr><tr><td>Cu</td><td>14.4 ± 2.62</td><td>27.23 ± 1.71</td><td>&lt;0.01</td></tr><tr><td>Zn</td><td>27.44 ± 6.36</td><td>67.35 ± 0.47</td><td>&lt;0.01</td></tr><tr><td>As</td><td>3.9 ± 0.64</td><td>6.87 ± 0.22</td><td>&lt;0.01</td></tr><tr><td>Ag</td><td>5.51 ± 3.7</td><td>5.1 ± 0.03</td><td>&lt;0.01</td></tr><tr><td>Pb</td><td>17.74 ± 2.54</td><td>11.4 ± 1.27</td><td>&lt;0.00</td></tr></tbody></table>

Values for microbial mats and sediments samples are given in (µg/g) and those of water samples are shown in parts per million (ppm).

members of genus *Methanospirillum*, which were abundant in mats (Figure 2B). Other metabolic pathways such as lysine biosynthesis, C5-branched dibasic acid metabolism, thiamine metabolism, pyrimidine metabolism, vitamin B6 metabolism, and other glycan degradation were also found to be abundant in microbial mat.

The S metabolism pathway could be reconstructed within a range of 76.31–84.21% which was maximum in sediment and minimum in microbial mat. A total of 75 genes were responsible for S metabolism present in all the samples with a mean copy number of 1580 ± 249.2. To gain insights into S disproportionation potential across all the habitats, we mapped the TIGRfam and Pfam (Supplementary File 4) ids of the 25 associated genes (mean copy number of 980 ± 222.1) on to the ORFs and copy numbers of these genes involved in S oxidation, sulfide oxidation, and sulfate reduction (as described in Figure 3B) were estimated. Thiosulfate ions are fused to the carrier complex of S-oxidizing proteins *soxYZ* (*soxY* = 145, *soxZ* = 71), while L-cysteine S-thiosulfotransferase (*soxX* = 76, *soxA* = 86) and S-sulfosfanyl-L-cysteine sulfohydrolase (*soxB* = 145) mediate the hydrolytic release of reduced S ions from S-bound *soxYZ*. Besides, the dissimilatory sulfite reductase (DsrAB) encoding genes: *dsrA* = 95, *dsrB* = 83 and anaerobic sulfite reductase subunits *asrA* = 44, *asrB* = 30, *asrC* = 9 reduces sulfite to sulfide (Supplementary File 4). The other 15 genes for reduction of sulfate ions included solute binding protein (*cysP* = 37), ATP binding protein (*cysA* = 168), transport system permease proteins (*cysU* = 158, *cysW* = 154), ATP sulfurylase (*sat* = 254), transferases (*cysC* = 348, bifunctional *cysNC* = 155, *cysN* = 75), phosphoadenosine phosphosulfate reductase (*cysH* = 314), adenylylsulfate reductase subunit A and

B (*aprA* = 88, *aprB* = 104), sulfite reductase flavoprotein alpha-component (*cysJ* = 111), sulfite reductase hemoprotein beta-component (*cysI* = 180), homocysteine desulfhydrase (*mccB* = 6), and cysteine synthase (ATCYSC1 = 4). Figure 3C shows the comparative abundances of these proteins in the three habitats. Here, the results revealed that the microbial S disproportionation occurs largely through the dissimilatory pathway carried out by DsrAB as compared to assimilatory sulfite reductase (AsrABC) mediated reduction. In nature, the dissimilatory pathway is shorter and thus, preferred route of microbial sulfate reduction (Figure 3B; Kushkevych et al., 2020).

We identified the sequences of S disproportionation producing significant alignments from the *nr* database for taxonomic confirmation and assigned each sequence that could be classified up to genus level to either SOB or SRB (Figure 4A and Supplementary File 5). The taxonomy and evolutionary phylogenetic topologies are discussed in detail in the next section.

## Diversification and Evolution of S Disproportionation Proteins

To gain insights in the differentiation of S disproportionation genes, we study the diversity and evolution of the key enzymes of SOB and SRB communities in this environmental biosystem. Therefore, we employed a two-step strategy of comparing similarities of all sequences in a pairwise fashion through SSN analysis and further estimated the measures of the rates of non-synonymous to synonymous substitutions in their orthologous proteins between each pair of habitats. SSN effectively resolves the pairwise similarities of each sequence (node) with every other sequence of an enzyme or a group of enzymes for a pathway such that any two nodes are connected by edges only if they share
<table><caption>TABLE 2 | Characteristics of sequenced datasets generated and assembled metagenomes obtained for each sample from different habitats.</caption><thead><tr><th rowspan="2">Attributes</th><th colspan="2">Microbial mat</th><th colspan="2">Sediment</th><th colspan="2">Water</th></tr><tr><th>KgM1</th><th>KgM2</th><th>KgS1</th><th>KgS2</th><th>KgW1</th><th>KgW2</th></tr></thead><tbody><tr><td>Sample name</td><td>110,861,650 × 100</td><td>152,895,302 × 100</td><td>148,273,248 × 100</td><td>145,473,498 × 100</td><td>144,915,928 × 100</td><td>141,139,772 × 100</td></tr><tr><td>(Clean reads × 100) bases</td><td>180,849</td><td>216,350</td><td>404,684</td><td>398,499</td><td>519,194</td><td>506,942</td></tr><tr><td>No. of contigs</td><td>329</td><td>397</td><td>592</td><td>583</td><td>600</td><td>591</td></tr><tr><td>Size of assembly (Mbp)</td><td>7,668</td><td>8,123</td><td>2,672</td><td>2,664</td><td>3,945</td><td>3,985</td></tr><tr><td>N<sub>50</sub></td><td>55.55 ± 12.97</td><td>54.0 ± 12</td><td>59.70 ± 10.74</td><td>59.61 ± 11</td><td>56.43 ± 12.97</td><td>56.34 ± 13</td></tr><tr><td>GC-content (%)</td><td>2.6</td><td>2.6</td><td>2.7</td><td>2.7</td><td>3.0</td><td>3.03</td></tr><tr><td>Shannon diversity index</td><td></td><td>0.194</td><td></td><td>0.334</td><td></td><td>0.416</td></tr><tr><td>Paired-sample t-test (p-value < 0.05)</td><td>305,932</td><td>369,390</td><td>570,660</td><td>562,428</td><td>545,606</td><td>537,498</td></tr><tr><td>Predicted proteins</td><td>452</td><td>595</td><td>560</td><td>570</td><td>1,105</td><td>1,081</td></tr><tr><td>CRISPRs</td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

sequence homology above a certain cutoff (here e-value of 1e-30). Thus, SSN provides for an accurate placement of a sequence among its putative homologs (Talwar et al., 2020).

Here, we examined the diversity among 19 key S substrates oxidizing and reducing proteins determined from the communities as shown in **Figure 3C**, except membrane permeases (CysPUWA) and genes for cysteine synthesis (MccB, ATCYSC1). In total, we retrieved 2,413 protein sequences (mean; *M* = 254, *S* = 480, *W* = 472) denoted by nodes in SSN. The network was organized into 88 connected components including 46 isolated nodes with an average clustering coefficient of 0.84. The connected components were represented by the homologous and heterologous clusters depending on whether they were constituted by the same gene or a number of different genes involved in a pathway, respectively. The number of connected components formed through SSN analysis of S metabolic proteins distributed into gene clusters and the isolated nodes denoted the diverging and highly diverged sequences, respectively (**Figure 4A**). Hence, we looked into these components to study the diversity of each gene that were distributed as shown in **Table 3**. The proteins CysNC, CysH, CysI, and CysJ catalyze important steps and act as cofactors for the AsrABC which were all found to be diverging with many isolated nodes and loosely connected components. The enzymes for S oxidation (Sox) were also found to be diverging as observed from loosely formed clusters. On the other hand, all sequences of the key enzyme of dissimilatory pathway, DsrAB, formed only one connected component, which suggested that they might be under convergent evolution at this site (**Figure 4A**). Further, we compared the distribution of diverged sequences that could be separated as isolated nodes and found that hot spring sediments harbored a high diversity of these enzymes (*n* = 23) in comparison with microbial mats (*n* = 12) and water (*n* = 11).

The node degree distribution estimated to be decreasing with increasing protein quantity (correlation = 0.52, *r*² = 0.29), average neighborhood connectivity within the networks interpreted as function in *k* was increasing and positively correlated (correlation = 0.90, *r*² = 0.74) (**Figure 4B**). Furthermore, closeness centrality curve that measures closeness between nodes was unable to reach the bench top (correlation = 0.03, *r*² = 0.01), might be due to the maximum number of connected components and less sequence homology. So, we also analyzed each protein cluster individually by using network analysis, 178/2,413 nodes of the network DsrAB protein cluster found to be conserved showed higher clustering coefficient values (0.94), followed by AprA (0.93), AprB, CysNC (0.88), and others (**Table 3**). The evolutionary selection pressures on these genes were studied through estimation of *dN/dS* values calculated for a subset of conserved gene sequences in all three habitats (**Figure 4C**). The number of core genes and the range of *dN/dS* values identified for each gene are shown in **Table 3**. The *cysI* and *cysI* genes were found to be under moderate selection pressures with *dN/dS* values in the range 0.4–0.7 (**Figure 4C** and **Supplementary File 6**). The results supported the observation as these enzymes code for important co-factors for the AsrABC that were found to be diverging in through SSN analysis. Therefore, the microbial
<figure><img src="image_5.png" /></figure>

**FIGURE 2** | Relative abundance of phylum and genera. (A) The stacked bar representation shows the dominating phylum in all three habitats. The dendrograms show hierarchical clustering between species and samples. (B) The abundant common genera *Pseudomonas* (20–60%), *Desulfobulbaceae_unclassifed* (15–20%), *Burkholderia*, *Desulfovibrio*, and *Thioalkalivibrio* (1–20%) are shown here relative abundance × log<sub>2</sub> scale. (C) Taxa based functional profiles demonstrating the major phylum contributing toward COG subsystems and percentage of proteins annotated within each COG category for the habitat sites. The distribution of the RPKs were mapped in accordance to the member abundance in the habitats. Dots were colored according to the phylum.

genes for assimilatory reduction pathway are diversifying under
moderate selection pressures.

Divergence, Phylogeny, and Structural
Relationships of Dissimilatory Sulfite
Reductase and Assimilatory Sulfite
Reductase

The enzymes catalyzing the reductive (Dsr) or oxidative (rDsr)
transformation between sulfite and sulfide appear to be related
with respect to their subunit composition and catalytic properties
(Loy et al., 2009). The dsr genes have been characterized from
bacterial as well as archaeal domains (Chang et al., 2001; Grim
et al., 2011; Colman et al., 2020). However, their evolution
in these domains has long remained a subject of discussion.
Our preliminary results showed that both subunits of dsr
genes (dsrA = 79, dsrB = 72) corresponds to about 70 newly
identified organism for both oxidation and reduction processes
(Supplementary File 8). Through RAxML phylogenetic analysis,
it can be confirmed that the dsrAB genes have been introduced in
most of the newly identified members by a multiple independent
LGT (Anantharaman et al., 2018). Importantly, organisms
from Acidobacteria, Candidate division Zixibacteria, Chloroflexi,
and β-Proteobacteria form completely novel lineages other
than known DsrAB clusters identified through RefSeq nr
protein database with accession numbers (National Center for
Biotechnology Information; Supplementary File 5). Hence,
the dsr from sulfate reducers formed a separate cluster, with
sequences from Desulfarculus, Desulfocarbo, Desulfarcinum,
Thermodesulfobacteria, Syntrophobacter, Desulfomonile,

Desulfovibrio, *Desulfatirhabdium*, and *Desulfobacteriaceae* in both DsrA and DsrB phylogenies and additionally, *Dissulfuribacter* in DsrB (Supplementary File 8). We proposed that these organisms with newly identified lineages of *dsr* genes involved in sulfite/sulfate oxidation and reduction likely serves an important control on S cycling on terrestrial subsurface. The divergence of *dsrAB* between unrelated taxa could be driven through combination of speciation, functional diversification, and LGT. Also, there is equal possibility of non-functionality of the genes in these taxa (Loy *et al.*, 2008; Anantharaman *et al.*, 2014).

Through time-scale evolutionary phylogeny of the sequences
of DsrAB and AsrAB identified from the metagenomes, we
determine the most basal and earliest evolved lineages involved
in dsr and rsdr pathways (Figure 5 and Supplementary File 7).
Our results suggested that both the subunits of the oxidative
type reverse-Dsr evolved much earlier than the reductive type
Dsr subunits (Figure 5). We used the phylogenetic analysis
to further assign taxa to the sequences that showed similarity
with yet uncultured bacteria and predicted their structures to
gain insights into the more common phylogenetic ancestor
of the two Dsr subunits. Interestingly, these sequences of the
oxidative type Dsr (rDsr) formed monophyletic clade with a more
recently identified genus, *Sulfuritortus* in both DsrA and DsrB
phylogenies. Our analysis revealed similar tree topologies with
these unassigned sequences forming clade with *Sulfuritortus*,
*Thiobacillus*, and *Hydrogenophilales bacterium* in both DsrA
(n = 1) and DsrB (n = 3). Although the sequences were similar
to *Thiobacillus* phylogenetically, prediction of their structures
revealed DsrA to be highly similar to that of *Desulfovibrio gigas*
A

<figure><img src="image_6.png" /></figure>

B

<figure><img src="image_7.png" /></figure>

C

<figure><img src="image_8.png" /></figure>

**FIGURE 3** | (A) Reconstruction of top 50 pathways annotated using KEGG automatic annotation server. Heatmap matrix representation and clustering was performed by using "pheatmap" package (Kolde and Kolde, 2015) in R (R Development Core Team, 2011). (B) The sulfate reduction pathway involved a group of reductases, kinases, and transferases with the product chemical structures generated through chemDraw7 and Inkscape v0.9 (Inkscape Project, 2020). (C) The gene copy number of both sulfate reduction and sulfide oxidation pathway that were partitioned in different habitats showed here using ggplot2 in R (R Development Core Team, 2011).

(TM score = 0.75; analog TM = 0.86; C-score = 0.25, 1.6, and 1.58; **Figure 6**) and structures of DsrB proteins aligned most closely with *Archaeoglobus fulgidus* (TM score = 0.94; analog TM = 0.97; C-score = 1.64; **Supplementary File 9**). While the other DsrB subunit (n = 2) formed clades with *Thioalkalivibrio* (*β*-**Proteobacteria**) also showing similarity to *A. fulgidus* (TM score = 0.94; analog TM = 0.98; C-score = 1.64). In the reductive type Dsr clades, two DsrA sequences from uncultured bacteria formed clades with *Nitrospiraceae* and *Chloroflexi* bacteria (**Figure 5**). However, a consistent observation was seen in the similarity of all these structures of DsrA proteins with *Desulfovibrio vulgaris* (TM score = 0.94; analog TM = 0.96; C-score = 1.58).

DISCUSSION

Hot water is continually discharged from a major outlet at Khirganga from where it deposits S upon microbial mats over the sediments along its course (Sharma et al., 2004). Microscopic analysis showed that cyanobacteria are widely distributed in mats and sedimentary deposits of thermal springs (van Gemerden, 1993; Podar et al., 2020). The physio-chemical data signified that the hot spring waters emerges from the confluence of rivers Parvati and Beas have high concentrations of chlorides

and sulfates that are characteristic of majority of other hot springs in the Himalayan ranges (Cinti et al., 2009; Sangwan et al., 2015). The alpha diversity was higher in the water samples than microbial mat and sediments as has been reported previously (Ghilamicael et al., 2017; Nagar et al., 2021). Species richness as rarefaction curves obtained for all samples attained a plateau indicating optimum metagenomic sequencing data and sampling of a reasonable number of species for all metagenomes. The Bray-Curtis index calculated and plotted using non-metric multidimensional scaling (NMDS) demonstrated a significant difference in the beta-diversity of all three habitats at phyla level (PERMANOVA; *p* < 0.01). Abundance of class *γ*-**Proteobacteria** in microbial mats significantly distinguished latter from the other two habitats. Similar results have been reported in previous studies (Selvarajan et al., 2018; Pohlner et al., 2019). Communities in sediment and water samples may be varied from each other majorly due to differences in the abundance of *δ*-**Proteobacteria**. The dominance of aerobic and facultative anaerobic bacteria like *Pseudomonas* in all three habitats could be possible due to mesothermic environment and natural subsurface water hydrodynamics (Nazina et al., 2005). The other taxa that are often found associated with mat deposits are active biofilm producers that use adherence to the surface as a strategy to survive, evolve, and to cope with various abiotic stresses at such extreme habitats (López et al., 2006; Wright et al., 2013).
A

Sequence Similarity Networks

<figure><img src="image_9.png" /></figure>

B

Network Analysis

<figure><img src="image_10.png" /></figure>

C

dN/dS

<figure><img src="image_11.png" /></figure>

**FIGURE 4** | Sequence similarity network analyses. (A) Diversity of sulfate reduction genes of both assimilatory and dissimilatory pathways in microbial mat (diamond-shaped), sediment (square), and water (spherical) habitats visualized in cytoscape v3.7.1. Highlighted only the classified taxa, where color cyan belong to SRBs and yellow to SOBs. The network was set at threshold e-value cutoff of 1e-30 and nodes represent sequences connected through edges if the similarity exceeds the cutoff score. Here, the clusters and isolated nodes were showing the conserved pattern and diversified pattern of the proteins significantly playing an important role in sulfate reduction. (B) Topological properties of the similarity networks: degree distribution, average clustering coefficient, average neighborhood connectivity, and closeness centrality are plotted against the number of neighbors. The power law fit curves are shown within each graph. (C) Habitat vs. habitat dN/dS values of all S cycle genes were estimated and plotted using xy-plot in R (R Development Core Team, 2011).

In contrast, abundance of SOB and SRB in sediment and water was observed (Agostino and Rosenbaum, 2018). These SOB and SRB are usually categorized as lithoautotrophs that play key microbial role in biogeochemical cycling of S in various habitats. In general, hydrogen sulfide account for the S present in the underground geothermal waters originating from pyrites or leaching of other sulfides by deep hypothermal waters (Picard et al., 2016). Sulfide (S<sup>2</sup><sup>-</sup>) is oxidized to sulfate (SO<sub>4</sub><sup>2</sup><sup>-</sup>) as the water rises to the surface and under mild oxidizing conditions,

sulfide is only oxidized to sulfate or S dioxide (Picard et al., 2016; Wu et al., 2021). The results provided pertinent information on the geochemical composition of the three habitats to be correlated with the microbial diversity and community functions. More importantly, high concentrations of sulfate ions in microbial mats and sedimentary deposits supported the hypothesis of a key role of the bacterial S cycle in sustaining the microbial community at the hot water spring. Sulfur oxidizing bacteria oxidize the reduced S compounds such as hydrogen sulfide
<table><caption>TABLE 3 | Attributes of the SS N and dN/dS analysis of S metabolism genes.</caption><thead><tr><th rowspan="2">Network parameters</th><th colspan="10">Network analysis of homologous and heterologous clusters</th></tr><tr><th>CysNC</th><th>Sat</th><th>CysH</th><th>DsrAB</th><th>CysI</th><th>CysJ</th><th>AprA</th><th>AprB</th><th>AsrAB</th><th>AsrC</th><th>SoxA</th><th>SoxB</th><th>SoxX</th><th>SoxY</th><th>SoxZ</th></tr></thead><tbody><tr><td>Sulfate reduction genes</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>No. of nodes</td><td>2,413</td><td>254</td><td>346</td><td>178</td><td>170</td><td>102</td><td>87</td><td>98</td><td>68</td><td>9</td><td>86</td><td>145</td><td>76</td><td>145</td><td>71</td></tr><tr><td>No. of edges</td><td>1,62,389</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Average degree</td><td>134.59</td><td>205.1</td><td>121.2</td><td>77.8</td><td>65.2</td><td>57.1</td><td>72.4</td><td>61.0</td><td>17.8</td><td>3.4</td><td>35.3</td><td>119.6</td><td>52.1</td><td>43.3</td><td>21.7</td></tr><tr><td>Connected components</td><td>88</td><td>17</td><td>14</td><td>1</td><td>11</td><td>8</td><td>2</td><td>3</td><td>6</td><td>4</td><td>4</td><td>2</td><td>2</td><td>7</td><td>7</td></tr><tr><td>Isolated nodes</td><td>46</td><td>12</td><td>8</td><td>0</td><td>10</td><td>5</td><td>1</td><td>0</td><td>4</td><td>2</td><td>0</td><td>0</td><td>0</td><td>2</td><td>2</td></tr><tr><td>Network density</td><td>0.05</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Characteristic path length</td><td>1.89</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Shortest path</td><td>14%</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Network centralization</td><td>0.16</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Clustering coefficient</td><td>0.84</td><td>0.86</td><td>0.80</td><td>0.94</td><td>0.85</td><td>0.86</td><td>0.93</td><td>0.88</td><td>0.79</td><td>0.76</td><td>0.97</td><td>0.91</td><td>0.89</td><td>0.83</td><td>0.83</td></tr><tr><td>Core genes clusters</td><td>-</td><td>14</td><td>6</td><td>5</td><td>3</td><td>2</td><td>1</td><td>2</td><td>1</td><td>0</td><td>2</td><td>4</td><td>3</td><td>6</td><td>2</td></tr><tr><td>Range dN/dS (≤0)</td><td>-</td><td>≥0.4</td><td>≥0.5</td><td>≥0.6</td><td>≥0.5</td><td>≥0.6</td><td>≥0.4</td><td>≥0.1</td><td>0</td><td>-</td><td>0</td><td>≥0.4</td><td>≥0.1</td><td>≥0.3</td><td>≥0.3</td></tr></tbody></table>

(H₂S), elemental S (S⁰), sulfite (SO₃²⁻), thiosulfate (S₂O₃²⁻), and various polythionates (SₙO₆²⁻ or -SₙO⁶⁻) into sulfate (SO₄²⁻). On the contrary, SO₄²⁻ can serve as an electron acceptor of SRB under anaerobic conditions, and they reduce the SO₄²⁻ and other oxidized S compounds (S₂O₃²⁻, SO₃²⁻, and S⁰) into H₂S (Agostino and Rosenbaum, 2018). The abundance of SOB such as *Thioalkalivibrio* and *Burkholderia* as well as SRB such as *Desulfobulbaceae* unclassified and *Desulfovibrio* in spring is not surprising as high levels of sulfate dominate the site and relative abundance of these bacteria provide evidence of an active S cycling mediated by microbial communities. The enriched diversity for Sox and sulfate reduction as well as the geochemical analysis of sulfide rich habitats compelled us to mine the regulatory genes involved in the different pathways of S cycle. In natural system, the S intermediates are reduced by different bacteria through two different reduction processes, namely, dissimilatory and assimilatory reactions (Vermeij and Kertesz, 1999; Zavarzin, 2008; **Figure 3B**). In dissimilatory reduction, SRB utilize three enzymes [(ATP sulfurylase (sat), APS reductase (apr), and sulfite reductase (dsr)] to reduce sulfate and produce toxic hydrogen sulfide (Agostino and Rosenbaum, 2018; Kushkevych et al., 2020). On contrary, sulfate is assimilated into organic compounds under assimilatory process (Kushkevych et al., 2020).

Based on sequence similarity, majority of the AsrABC genes (assimilatory) that were taxonomically related to SRB could be distinguished from the rest of the sequences that were not identified as either SOB or SRB. Thus, assimilatory reduction was diverged among the SRB communities while those for dissimilatory pathway were rather conserved at the site. Previous reports of sequence comparisons have confirmed that DsrAB, Dsr, to be highly conserved enzyme that could serve as marker gene for SRBs (Loy et al., 2009). The DsrAB and AprAB enzymes were contributed by both SOB and SRB with syntrophic interaction which suggests for the presence of reverse dsr (rdsr) mediated oxidation of S substrates in addition to dissimilatory reduction (Kumar et al., 2017). The inherent complexity of S-based metabolic network revealed that there are controlled mutation rates in dsrAB genes in presence or increased selective pressure of contamination and extreme conditions. A relatively high diversity of the other sulfate disproportionation proteins in all three segments unveiled the high nutritional demands and efficiency of the microbes toward uptake of a wide range of structurally and chemically diverse amino acid side chains from environment (Talwar et al., 2020). The syntropy of SOB and SRB prevailing in anoxic and anaerobic conditions governs the dissimilatory S metabolism (oxidation and reduction parallelly) and indirectly promoting the growth of diverse microbes in this natural ecosystem (Bhatnagar et al., 2020).

It was noted that mutation rates were low in all subunits of Dsr proteins which we readily analyzed to trace down the ancestor among SOBs or SRBs in dissimilatory and reverse dissimilatory (rdsr) pathways. The result of the time-scale phylogeny suggested that an increase in substitution rate in both subunits of DSR might have occurred on the branch connecting *δ-Proteobacteria* to all other taxa as
A Diseliminatory reductase subunits (DsrAB)

<figure><img src="image_12.png" /></figure>

B Anisomilary reductase subunits (AsrAB)

<figure><img src="image_13.png" /></figure>

**FIGURE 5** | Divergence estimation over time. Reconstruction of the phylogenetic tree of optimized full length DsrAB and AsrABC subunits in three habitats using PhyloBayes with the CAT-GTR model. The highlighted squares consist of clades with proteins that were remained unclassified through *nr* database. (A) (i) Among, 78 DsrA nodes that showed here the earliest evolution of the rDsr oxidative proteins occurred in *Thiobacillus* sp. (ii) 72 nodes of DsrB proteins with similar results. (B) (i) 38 nodes of AsrA and (ii) 21 nodes of AsrB were also compared as a control for branch length shown here.

A

<figure><img src="image_14.png" /></figure>

B

<figure><img src="image_15.png" /></figure>

**FIGURE 6** | Structural similarity and analogy of unclassified proteins from PDB (Protein data bank) using i-TASSER. (A) DsrA protein subunits: (i) np252939 and (ii); (iii) np500880; np481424 showing Tm-align similarity with PDB ids 3or1 and 2v4j (Desulfovibrio gigas and D. vulgaris), respectively; with SF4 (iron S cluster) ligand binding sites (B) DsrB protein subunits (i) np24977, np24617, np116484 and (ii) np199275, np59533 showing Tm-align similarity with PDB ids 3c7b; with siroheme ligand binding sites.
observed from the branch lengths. The different rates of substitution of the two DSR subunits has this far only been reported in *δ-Proteobacteria* lineage (Wagner et al., 1998). The templates of these *dsrAB* genes have potential to study the genotypic and phenotypic traits in SRPs and the dissimilatory S metabolism processes which will expand the gene-environment interaction mechanism. Also, prior analysis has proved that the evolution of *dsrAB* have been influenced by LGT only among major taxonomic lineages (Klein et al., 2001; Müller et al., 2015) but the findings here provide evidence of independent multiple LGT events distributed throughout the dissimilatory gene clusters. Currently, the time-scale study of this site cannot produce evidence of the progenitor lineages, as the evolutionary history of dissimilatory reduction is complex and yet ascertain. Although it had provided information of the earliest lineage where sulfate/sulfite oxidation and reduction appeared. The genus of SOB currently has few recognized species and is closely related to the members of genus *Thiobacillus* (Kojima et al., 2017). Through structural homology, we predicted that the genus derives its two subunits of Dsr from different ancestors. A plausible explanation for this is observed in the previous reports of high sequence homology between the Dsr of *A. fulgidus* and *D. vulgaris* that suggest a common origin of archaeal and bacterial DSRs or their HGT (Karkhoff-Schweizer et al., 1995). In addition to homology in their sequences, the evolutionary distance separating the enzymes from *A. fulgidus* and *D. vulgaris* was deciphered. For DsrB subunits, the archaeal and bacterial sequences were not particularly distant; such that the branches with structural homology to *A. fulgidus* were approximately the same length as branches leading to bacteria such as *Thiobacillus*, *Thioalkalivibrio*, and *β-Proteobacteria bacterium* (Larsen et al., 1999).

## CONCLUSION

The mesothermic hot spring have been composed of a diverse group of microbes (Bacteria and Archaea) and genotypes (*dsrAB*) that could be screened out as novel thermozymes that cannot be underestimated. From the results, it could be concluded that the microbial community functions were distinguished in microbial mat from water and sediment. Here, the genomic repertoire suggested the ongoing specific adaptations to cope up with extreme values of sulfide content in this ecological setting. The S metabolic pathways are completed where inorganic S compounds being the main source for SRB releasing toxicity in the form of sulfides (S²⁻). The sulfate reduction profiling in all three habitats reveals dissimilatory sulfate reduction process (*dsrAB*) is active than assimilatory sulfate reduction (*asrAB*). Later, the genes involved in S reduction/oxidation were classified and belong mostly to *Proteobacteria* with maximum homologous proteins classified in anoxygenic SOBs. In all S disproportionation proteins, the sulfite reductase DsrAB proteins showed conserved behavior with 0/1 isolated nodes that have been signified as phylogenetic markers for SRBs. The evolutionary phylogenetic

analysis showed that the oxidative rDsr were the earliest than the reductive Dsr which may predict that the condition with more sulfides oxidized in more sulfates, directing SRB to perform dissimilatory reduction later. Phylogenetic clades of DsrAB proteins showed unanimous distribution of taxa except the *δ-Proteobacteria* which could be the reason for occurring LGT to other phyla. On the basis of structural alignments, the lineages with unclassified clades have shown different analogy in both Dsr subunits where DsrB derived from Archaea and DsrA are *δ-Proteobacteria* in origin at this mesothermic niche. The stabilization and evolutionary time-scale phylogeny of DsrAB revealing a positive syntrophic relationship between SOB and SRB. These thermophilic microbial inhabitants are very crucial in expanding the metal toxification, ion exchange, and biogeochemical cycling of the elements.

## DATA AVAILABILITY STATEMENT

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

## AUTHOR CONTRIBUTIONS

SN: conceptualization, methodology, investigation, formal analysis, data curation, methodology, writing—original draft, review, and editing. CT: formal analysis, data curation, methodology, writing—original draft, review, and editing. MM-H and H-HR: formal analysis, writing—review, and editing. MS: writing—review and editing. RL and RN: conceptualization, Writing—review, and editing. All authors contributed to the article and approved the submitted version.

## FUNDING

This work was supported by funds from the Department of Biotechnology (DBT), National Bureau of Agriculturally Important Microorganisms (NBAIM), University Grants Commission-Career Advancement Scheme, and Department of Science and Technology-Purse grant.

## ACKNOWLEDGMENTS

The sequence data were produced by the Beijing Genome Institute, China (BGI Tech. Solutions Co., Ltd., Hongkong) in collaboration with the user community. We thank Matthias Schmidt for providing the electron micrographs at the Center for Chemical Microscopy (ProVIS) scanning at UFZ Leipzig, which was supported by European Regional Development Funds (EFRE-Europe Funds Saxony), and the Helmholtz Association, and The Energy and Resources Institute (TERI). SN and CT thank Council
of Scientific and Industrial Research (CSIR) for providing doctoral fellowships.

## SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2022.848010/full#supplementary-material

**Supplementary File 1** | NCBI accession and *in-situ* measures of temperatures, pH, and elemental analysis of each sample from three studied niches.

**Supplementary File 2** | (A) Rarefaction curves based on measures of alpha diversity species richness calculated for each of the six samples. (B) The interquartile range of alpha diversity measures based on Shannon indices. (C) Non-metric multidimensional scaling plots based on Bray-Curtis indices of community dissimilarities in the six samples types (*p* < 0.05, PERMANOVA).

**Supplementary File 3** | The list of all metabolic gene families in each sample (KgM1, KgM2, KgS1, KgS2, KgW1, and KgW2) and the RPKs abundances of taxa (*Proteobacteria*, *Bacteroidetes*, etc.) summarized in details.

**Supplementary File 4** | Detailed information of percentage reconstruction of metabolic pathways; Common and unique pathways in all six samples; Identification sources for 19 important proteins involved in S disproportionation using TIGRfam and Pfam conserved domain database, x indicates no hits found. Also, the copy number of 25 S disproportionation proteins recorded.

**Supplementary File 5** | Details of taxonomic information of 2,413 nodes of 19 S cycle genes identified from nr NCBI microbial proteins database.

**Supplementary File 6** | Habitat vs. habitat *dN/dS* values of the S cycling genes.

**Supplementary File 7** | Sequence alignment of DsrA, DsrB, AsrA, and AsrB proteins identified in this study.

**Supplementary File 8** | Phylogenetic analysis of DsrAB proteins: showing percentage similarity with taxa on NCBI nr database. Few newly identified lineages were marked in red. In DsrAB, maximum number of LGT events were occurred in Proteobacteria class itself *clade (A2,B2)* while the independent events of LGT toward different phyla were occurred through δ-Proteobacteria depicted in *clade (A1,B1)*.

**Supplementary File 9** | Structural models for other ligand binding sites of unclassified proteins of DsrA and DsrB subunits.

## REFERENCES

Agostino, V., and Rosenbaum, M. A. (2018). Sulfate-reducing electro-autotrophs and their applications in bioelectrochemical systems. *Front. Energy Res.* 6:55. doi: 10.3389/fenrg.2018.00055

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. *J. Mol. Biol.* 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2

Anantharaman, K., Duhaime, M. B., Breier, J. A., Wendt, K. A., Toner, B. M., and Dick, G. J. (2014). Sulfur oxidation genes in diverse deep-sea viruses. *Science* 344, 757–760. doi: 10.1126/science.1252229

Anantharaman, K., Hausmann, B., Jungbluth, S. P., Kantor, R. S., Lavy, A., Warren, L. A., et al. (2018). Expanded diversity of microbial groups that shape the dissimilatory sulfur cycle. *ISME J.* 12, 1715–1728. doi: 10.1038/s41396-018-0078-0

Ayangbenro, A. S., and Babalola, O. O. (2017). A new strategy for heavy metal polluted environments: a review of microbial biosorbents. *Int. J. Environ. Res. Public Health* 14:94. doi: 10.3390/ijerph14010094

Ayangbenro, A. S., Olanrewaju, O. S., and Babalola, O. O. (2018). Sulfate-reducing bacteria as an effective tool for sustainable acid mine bioremediation. *Front. Microbiol.* 9:1986. doi: 10.3389/fmicb.2018.01986

Battaglia-Brunet, F., Crouzet, C., Burnol, A., Coulon, S., Morin, D., and Joulian, C. (2012). Precipitation of arsenic sulphide from acidic water in a fixed-film bioreactor. *Water Res.* 46, 3923–3933. doi: 10.1016/j.watres.2012.04.035

Bhatnagar, S., Cowley, E. S., Kopf, S. H., Pérez Castro, S., Kearney, S., and Dawson, S. C. (2020). Microbial community dynamics and coexistence in a sulfide-driven phototrophic bloom. *Environ. Microb.* 15:3. doi: 10.1186/s40793-019-0078-0

Boopathy, R. (2014). Biodegradation of 2, 4, 6-trinitrotoluene (TNT) under sulfate and nitrate reducing conditions. *Biologia* 69, 1264–1270. doi: 10.2478/s11756-014-0441-1

Bowen de León, K., Gerlach, R., Peyton, B. M., and Fields, M. W. (2013). Archaeal and bacterial communities in three alkaline hot springs in Heart Lake Geyser Basin, Yellowstone National Park. *Front. Microbiol.* 4:330. doi: 10.3389/fmicb.2013.00030

Callaghan, A. V., Morris, B. E., Pereira, I. A., McInerney, M. J., Austin, R. N., Groves, J. T., et al. (2012). The genome sequence of *Desulfatibaculum alkenivorans* AK-01: a blueprint for anaerobic alkane oxidation. *Environ. Microbiol.* 14, 101–113. doi: 10.1111/j.1462-2920.2011.02516.x

Chan, C. S., Chan, K. G., Tay, Y. L., Chua, Y. H., and Goh, K. M. (2015). Diversity of thermophiles in a Malaysian hot spring determined using 16S rRNA and shotgun metagenome sequencing. *Front. Microbiol.* 6:177. doi: 10.3389/fmicb.2015.0017

Chang, Y. J., Peacock, A. D., Long, P. E., Stephen, J. R., McKinley, J. P., MacNaughton, S. J., et al. (2001). Diversity and characterization of sulfate-reducing bacteria in groundwater at a uranium mill tailings site. *Appl. Environ. Microbiol.* 67, 3149–3160. doi: 10.1128/AEM.67.7.3149-3160. 2001

Cinti, D., Pizzino, L., Voltattorni, N., Quattrochi, F., and Walia, V. (2009). Geochemistry of thermal waters along fault segments in the Beas and Parvati valleys (north-west Himalaya, Himachal Pradesh) and in the Sohn a town (Haryana), India. *Geochem. J.* 43, 65–76. doi: 10.2343/geochemj.1.0011

Colman, D. R., Lindsay, M. R., Amenabar, M. J., Fernandes-Martins, M. C., Roden, E. R., and Boyd, E. S. (2020). Phylogenomic analysis of novel *Diaforarchaea* is consistent with sulfite but not sulfate reduction in volcanic environments on early Earth. *ISME J.* 14, 1316–1331. doi: 10.1038/s41396-020-0611-9

Cox, M. P., Peterson, D. A., and Biggs, P. J. (2010). SolexaQA: at-a-glance quality assessment of Illumina second-generation sequencing data. *BMC Bioinformatics* 11:485. doi: 10.1186/1471-2105-11-485

Dahl, C., and Truper, H. G. (1994). Enzymes of dissimilatory sulfide oxidation in phototrophic sulfur bacteria. *Methods Enzymol.* 243, 400–421. doi: 10.1016/0076-6879(94)43030-6

Dong, Y., Sanford, R. A., Inskeep, W. P., Srivastava, V., Bulone, V., Fields, C. J., et al. (2019). Physiology, metabolism, and fossilization of hot-spring filamentous microbial mats. *Astrobiology* 19, 1442–1458. doi: 10.1089/ast.2018.1965

Eddy, S. R. (2011). Accelerated profile HMM searches. *PLoS Comput. Biol.* 7:e1002195. doi: 10.1371/journal.pcbi.1002195

Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. *Nucleic Acids Res.* 32, 1792–1797. doi: 10.1093/nar/gkh340

EPA (2001). *Parameters of Water Quality: Interpretation and Standards*. Wexford: Environmental Protection Agency.

Eren, A. M., Vineis, J. H., Morrison, H. G., and Sogin, M. L. (2013). A filtering method to generate high quality short reads using Illumina paired-end technology. *PLoS One* 8:e66664. doi: 10.1371/journal.pone.0066664

Finn, R. D., Bateman, A., Clements, J., Coggill, P., Eberhardt, R. Y., and Eddy, S. R. (2014). Pfam: the protein families database. *Nucleic Acids Res.* 42, D222–D230. doi: 10.1093/nar/gkt1223

Finster, K. (2011). Microbiological disproportionation of inorganic sulfur compounds. *J. Sulf. Chem.* 29, 281–292. doi: 10.1080/17415990802105770

Franzosa, E. A., McIver, L. J., Rahnavard, G., Thompson, L. R., Schirmer, M., Weingart, G., et al. (2018). Species-level functional profiling of metagenomes and metatranscriptomes. *Nat. Methods* 15, 962–968. doi: 10.1038/s41592-018-0176-y
Ghilamicael, A. M., Budambula, N. L. M., and Anami, S. E. (2017). Evaluation of prokaryotic diversity of five hot springs in Eritrea. BMC Microbiol. 17:203. doi: 10.1186/s12866-017-1113-4

Ghosh, W., Mallick, S., Haldar, P. K., Pal, B., Maikap, S. C., and Gupta, S. K. D. (2012). Molecular and cellular fossils of a mat-like microbial community in geothermal boracic sinters. Geomicrobiol. J. 29, 879–885. doi: 10.1080/01490451. 2011.635761

Gomez-Alvarez, V., Teal, T. K., and Schmidt, T. M. (2009). Systematic artifacts in metagenomes from complex microbial communities. ISME J. 3, 1314–1317. doi: 10.1038/ismej.2009.72

Gonsior, M., Hertkorn, N., and Hinman, N. (2018). Yellowstone hot springs are organic chemodiversity hot spots. Sci. Rep. 8:14155. doi: 10.1038/s41598-018- 32593-x

Grim, F., Franz, B., and Dahl, C. (2011). Regulation of dissimilatory sulfur oxidation in the purple sulfur bacterium *Allochromatium vinosum*. Front. Microbiol. 2:51. doi: 10.3389/fmicb.2011.00051

Haferburg, G., and Kothe, E. (2007). Microbes and metals: interactions in the environment. J. Basic Microbiol. 47, 453–467. doi: 10.1002/jobm.200700275

Haft, D. H., Selengut, J. D., Richter, R. A., Harkins, D., Basu, M. K., and Beck, E. (2012). TIGRFAMs and genome properties in 2013. Nucleic Acids Res. 41, 387–D395. doi: 10.1093/nar/gks1234

Hipp, W. M., Pott, A. S., Thum-Schmitz, N., Faath, I., Dahl, C., and Truper, H. G. (1997). Towards the phylogeny of APS reductases and sirohaem sulfite reductases in sulfate-reducing and sulfur-oxidizing prokaryotes. Microbiology 143, 2891–2902. doi: 10.1099/00221287-143-9-2891

Huerta-Cepas, J., Forslund, K., Coelho, L. P., Szklarczyk, D., Jensen, L. J., and Von Mering, C. (2017). Fast genome-wide functional annotation through orthology assignment by eggNOG-mapper. Mol. Biol. Evol. 34, 2115–2122. doi: 10.1093/ molbev/msx148

Hyatt, D., Chen, G. L., Locascio, P. F., Land, M. L., Larimer, F. W., and Hauser, L. J. (2010). Prodigal: prokaryotic gene recognition and translation initiation site identification BMC Bioinformatics 11:119. doi: 10.1186/1471-2105-1 1-119

Inkscape Project (2020). Available online at: https://inkscape.org/en/

Inskeep, W. P., Rusch, D. B., Jay, Z. J., Herrgard, M. J., Kozubal, M. A., Richardson, T. H., et al. (2010). Metagenomes from high-temperature chemotrophic systems reveal geochemical controls on microbial community structure and function. PLoS One 5:e9973. doi: 10.1371/journal.pone.0009773

Jaekel, U., Zedelius, J., Wilkes, H., and Musat, F. (2015). Anaerobic degradation of cyclohexane by sulfate-reducing bacteria from hydrocarbon-contaminated marine sediments. Front. Microbiol. 6:116. doi: 10.3389/fmicb.2015.00116

Kanehisa, M., Goto, S., Kawashima, S., Okuno, Y., and Hattori, M. (2004). The KEGG resource for deciphering the genome. Nucleic Acids Res. 32, 277–280. doi: 10.1093/nar/gkh063

Karkhoff-Schweizer, R. R., Huber, D. P., and Voordouw, G. (1995). Conservation of the genes for dissimilatory sulfite reductase from Desulfovibrio vulgaris and Archaeoglobus fulgidus allows their detection by PCR. Appl. Environ. Microbiol. 61, 290–296. doi: 10.1128/aem.61.1.290-296.1995

Klein, M., Friedrich, M., Roger, A. J., Hugenholtz, P., Fishbain, S., and Abicht, H. (2001). Multiple lateral transfers of dissimilatory sulfite reductase genes between major lineages of sulfate-reducing prokaryotes. J. Bacteriol. 183, 6028– 6035. doi: 10.1128/JB.183.20.6028-6035.2001

Kojima, H., Watanabe, M., and Fukui, M. (2017). Sulfuritortus calidifontis gen. nov., sp. nov., a sulfur oxidizer isolated from a hot spring microbial mat. Int. J. Syst. Evol. microbiol. 67, 1355–1358. doi: 10.1099/ijsem.0. 001813

Kolde, R., and Kolde, M. R. (2015). Package 'Pheatmap'. R Package 1, 790.

Kumar, S. S., Malyan, S. K., Basu, S., and Bishnoi, N. R. (2017). Syntrophic association and performance of Clostridium, Desulfovibrio, Aeromonas and Tetrathiobacter as anodic biocatalysts for bioelectricity generation in dual chamber microbial fuel cell. Environ. Sci. Pollut. Res. Int. 24, 16019–16030. doi: 10.1007/s11356-017-9112-4

Kushkevych, I., Abdulina, D., Kováë, J., Dordeviæ, D., Vítizov, M., Iutynska, G., et al. (2020). Adenosine-5'-Phosphosulfate-and sulfite reductases activities of sulfate-reducing bacteria from various environments. Biomolecules 10:921. doi: 10.3390/biom10060921

Larsen, O., Lien, T., and Birkeland, N. K. (1999). Dissimilatory sulfite reductase from Archaeoglobus profundus and Desulfotomaculum thermocisternum:

phylogenetic and structural implications from gene sequences. *Extremophiles* 3, 63–70. doi: 10.1007/s007920050100

Lartillot, N., Rodrigue, N., Stubbs, D., and Richer, J. (2009). Phylobayes-MPI 1.7: a Bayesian software for phylogenetic reconstruction using mixture models. *Syst. Biol.* 62, 611–615. doi: 10.1093/sysbio/syt022

Li, S. J., Hua, Z. S., and Huang, L. N. (2014). Microbial communities evolve faster in extreme environments. Sci. Rep. 4:6205. doi: 10.1038/srep06205

Li, W., and Godzik, A. (2006). Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics 22, 1658–1659. doi: 10.1093/bioinformatics/btl158

López, M. A., Zavala-Díaz de la Serna, F. J., Jan-Roblero, J., Romero, J. M., and Hernández-Rodriguez, C. (2006). Phylogenetic analysis of a biofilm bacterial population in a water pipeline in the Gulf of Mexico. FEMS Microbiol. Ecol. 58, 145–154. doi: 10.1111/j.1574-6941.2006.00137

Loy, A., Duller, S., and Wagner, M. (2008). "Evolution and ecology of microbes dissimilating sulfur compounds: insights from siroheme sulfite reductases," in Microbial Sulfur Metabolism, ed. C. Dahl (Berlin: Springer), 46–59. doi: 10.1007/ 978-3-540-72682-1_5

Loy, A., Duller, S., Baranyi, C., Mussmann, M., Ott, J., Sharon, I., et al. (2009). Reverse dissimilatory sulfite reductase as phylogenetic marker for a subgroup of sulfur-oxidizing prokaryotes. Environ. Microbiol. 11, 289–299. doi: 10.1111/ j.1462-2920

Meckenstock, R. U., Boll, M., Mouttaki, H., Koelschbach, J. S., Cunha Tarouco, P., Weyrauch, P., et al. (2016). Anaerobic degradation of benzene and polycyclic aromatic hydrocarbons. J. Mol. Microbiol. Biotechnol. 26, 92–118. doi: 10.1159/ 000441358

Meyer, F., Paarmann, D., D'Souza, M., Olson, R., Glass, E. M., Kubal, M., et al. (2008). The metagenomics RAST server - a public resource for the automatic phylogenetic and functional analysis of metagenomes. BMC Bioinformatics 9:386. doi: 10.1186/1471-2105-9-386

Kamarisima, Miyanaga, K., and Tanji, Y. (2019). The utilization of aromatic hydrocarbon by nitrate- and sulfate-reducing bacteria in single and multiple nitrate injection for souring control. Biochem. Eng. J. 143, 75–80. doi: 10.1016/ j.bej.2018.12.006

Mothe, G. K., Pakshirajan, K., and Das, G. (2016). Heavy metal removal from multicomponent system by sulfate reducing bacteria: mechanism and cell surface characterization. J. Hazard. Mater. 324, 62–70. doi: 10.1016/j.jhazmat. 2015.12.042

Mulla, S. I., Talwar, M. P., and Ninnekar, H. Z. (2014). "Bioremediation of 2, 4, 6-trinitrotoluene explosive residues," in Biological Remediation of Explosive Residues, ed. S. N. Singh (Cham: Springer), 201–233. doi: 10.1007/978-3-319- 01083-0

Müller, A. L., Kjeldsen, K. U., Rattei, T., Pester, M., and Loy, A. (2015). Phylogenetic and environmental diversity of DsrAB-type dissimilatory (bi)sulfite reductases. ISME J. 9, 1152–1165. doi: 10.1038/ismej.2014.208

Nagar, S., Talwar, C., Bharti, M., Yadav, S., Siwach, S., and Negi, R. K. (2021). Metagenome-assembled genomes recovered from the datasets of a high-altitude Himalayan hot spring Khirganga, Himachal Pradesh, India. Data Brief. 39:107551. doi: 10.1016/j.dib.2021.107551

Nazina, T. N., Sokolova, D. Sh, Shestakova, N. M., Grigor'ian, A. A., Mikhailova, E. M., Babich, T. L., et al. (2005). The phylogenetic diversity of aerobic organotrophic bacteria from the Dagan high-temperature oil field. Mikrobiologija 74, 401–409.

Peng, Y., Leung, H. C., Yiu, S. M., and Chin, F. Y. (2012). IDBA-UD: a de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth. Bioinformatics 28, 1420–1428. doi: 10.1093/bioinformatics/bts174

Picard, A., Gartman, A., and Girguis, P. R. (2016). What do we really know about the role of microorganisms in iron sulfide mineral formation? Front. Earth Sci. 4:68. doi: 10.3389/feart.2016.00068

Podar, P. T., Yang, Z., Björnsdóttir, S. H., and Podar, M. (2020). Comparative analysis of microbial diversity across temperature gradients in hot springs from Yellowstone and Iceland. Front. Microbiol. 11:1625. doi: 10.3389/fmicb.2020. 01625

Poddar, A., and Das, S. K. (2018). Microbiological studies of hot springs in India: a review. Arch. Microbiol. 200, 1–18. doi: 10.1007/s00203-017-1429-3

Pohlner, M., Dlugosch, L., Wemheuer, B., Mills, H., Engelen, B., and Reese, B. K. (2019). The majority of active Rhodobacteraceae in marine sediments belong to uncultured genera: a molecular approach to link their distribution
to environmental conditions. *Front. Microbiol.* 10:659. doi: 10.3389/fmicb.2019.00659

R Development Core Team (2011). *R: A Language and Environment for Statistical Computing*. Vienna: R Foundation for Statistical Computing.

Roy, C., Rameez, M. J., and Haldar, P. K. (2020). Microbiome and ecology of a hot spring-microbialite system on the Trans-Himalayan Plateau. *Sci. Rep.* 10:5917. doi: 10.1038/s41598-020-62797-z

Sahinkaya, E., Yurtsever, A., Toker, Y., Elcik, H., Cakmaci, M., and Kaksonen, A. H. (2015). Biotreatment of As-containing simulated acid mine drainage using laboratory scale sulfate reducing upflow anaerobic sludge blanket reactor. *Miner. Eng.* 75, 133-139. doi: 10.1016/j.mineng.2014.08.012

Sangwan, N., Lambert, C., Sharma, A., Gupta, V., Khurana, P., Khurana, J. P., et al. (2015). Arsenic rich Himalayan hot spring metagenomics reveal genetically novel predator-prey genotypes. *Environ. Microbiol. Rep.* 7, 812–823. doi: 10.1111/1758-2229.12297

Saxena, R., Dhakan, D. B., Mittal, P., Waiker, P., Chowdhury, A., Ghatak, A., et al. (2017). Metagenomic analysis of hot springs in Central India reveals hydrocarbon degrading thermophiles and pathways essential for survival in extreme environments. *Front. Microbiol.* 7:2123. doi: 10.3389/fmicb.2016.02123

Seemann, T. (2014). Prokka: rapid prokaryotic genome annotation. *Bioinformatics* 30, 2068–2069. doi: 10.1093/bioinformatics/btu153

Selvarajan, R., Sibanda, T., and Tekere, M. (2018). Thermophilic bacterial communities inhabiting the microbial mats of “indifferent” and chalybeate (iron-rich) thermal springs: diversity and biotechnological analysis. *Microbiol. Open* 7:e00560. doi: 10.1002/mbo.3.560

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. *Genome Res.* 13, 2498–2504. doi: 10.1101/gr.1239303

Sharma, A. K., Sharma, R., and Dandi, H. R. (2004). *Mineral Resources of Himachal Pradesh*. Available online at: https://inkscape.org/en/

Shirkot, P., and Verma, A. (2015). Assessment of thermophilic bacterial diversity of thermal springs of Himachal Pradesh. *ENVIS Bull. Himal. Ecol.* 23, 27–34. doi: 10.1007/s40009-018-0676-4

Singh, A., and Subudhi, E. (2016). Structural insights of microbial community of Deulajhari (India) hot spring using 16s-rRNA based metagenomic sequencing. *Genom. Data* 7, 101–102. doi: 10.1016/j.gdata.2015.12.004

Stamatakis, A. (2014). RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. *Bioinformatics* 30, 1312–1313. doi: 10.1093/bioinformatics/btu033

Stasik, S., Wick, L. Y., and Wendt-Potthoff, K. (2015). Anaerobic BTEX degradation in oil sands tailings ponds: impact of labile organic carbon and sulfate-reducing bacteria. *Chemosphere* 138, 133–139. doi: 10.1016/j.chemosphere.2015.05.068

Suzek, B. E., Wang, Y., Huang, H., McGarvey, P. B., and Wu, C. H. (2015). UniRef clusters: a comprehensive and scalable alternative for improving sequence similarity searches. *Bioinformatics* 31, 926–932. doi: 10.1093/bioinformatics/btu739

Talwar, C., Nagar, S., Kumar, R., Scaria, J., Lal, R., and Negi, R. K. (2020). Defining the environmental adaptations of genus *Devosia*: insights into its expansive short peptide transport system and positively selected genes. *Sci. Rep.* 10:1151. doi: 10.1038/s41598-020-58163-8

Tavaré, S. (1986). Some probabilistic and statistical problems in the analysis of DNA sequences, lectures on mathematics in the life sciences. *Am. Math. Soc.* 17, 57–86.

Truong, D. T., Franzosa, E. A., Tickle, T. L., Scholz, M., Weingart, G., Pasolli, E., et al. (2015). MetaPhlAn2 for enhanced metagenomic taxonomic profiling. *Nat. Methods* 12, 902–903. doi: 10.1038/nmeth.3589

van Gemerden, H. (1993). Microbial mats: a joint venture. *Mar. Geol.* 113, 3–25. doi: 10.1016/0025-3227(93)90146-M

Varin, T., Lovejoy, C., Jungblut, A. D., Vincent, W. F., and Corbeil, J. (2010). Metagenomic profiling of Arctic microbial mat communities as nutrient scavenging and recycling systems. *Limnol. Oceanogr.* 55, 1901–1911. doi: 10.4319/lo.2010.55.5.1901

Vermeij, P., and Kertesz, M. (1999). Pathways of Assimilative Sulfur Metabolism in *Pseudomonas putida*. *J. Bacteriol.* 181, 5833–5837. doi: 10.1128/JB.181.18.5833-5837.1999

Wagner, M., Roger, A. J., Flax, J. L., Brusseau, G. A., and Stahl, D. A. (1998). Phylogeny of dissimilatory sulfite reductases supports an early origin of sulfate respiration. *J. Bacteriol.* 180, 2975–2982. doi: 10.1128/JB.180.11.2975-2982.1998

Wickham, H. (2009). Elegant graphics for data analysis. *Media* 35, 10–1007.

Wright, K. E., Williamson, C., Grasby, S. E., Spear, J. R., and Templeton, A. S. (2013). Metagenomic evidence for sulfur lithotrophy by Epsilonproteobacteria as the major energy source for primary productivity in a sub-aerial arctic glacial deposit, Borup Fiord Pass. *Front. Microbiol.* 4:63. doi: 10.3389/fmicb.2013.00063

Wu, B., Liu, F., Fang, W., Yang, T., Chen, G. H., He, Z., et al. (2021). Microbial sulfur metabolism and environmental implications. *Sci. Total Environ.* 15:778. doi: 10.1016/j.scitotenv.2021.146085

Yang, J., and Zhang, Y. (2015). Protein structure and function prediction using I-TASSER. *Curr. Protoc. Bioinform.* 52, 5.8.1–5.815. doi: 10.1002/0471250953.bi0508s52

Zavarzin, G. A. (2008). *Microbial Cycles, Encyclopedia of Ecology*. Cambridge, MA: Academic Press, 2335–2341. doi: 10.1016/B978-008045405-4.00745-X

Zhang, J., Tian, Y., Zhang, J., Li, N., Kong, L., Yu, M., et al. (2016). Distribution and risk assessment of heavy metals in sewage sludge after ozonation. *Environ. Sci. Pollut. Res.* 24, 5118–5125. doi: 10.1007/s11356-016-6313-1

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

**Publisher's Note:** All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Copyright © 2022 Nagar, Talwar, Motelica-Heino, Richnow, Shakarad, Lal and Negi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.