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⁺, Mg²⁺, and Ca²⁺) and anions (SO₄^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.

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 Q₂₀ 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.

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:

Diversity at phylum level was inferred from MG-RAST (maximum e-value, 1 × 10⁵ 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.

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:

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:

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.

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:

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.

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).

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₄^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.

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).

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 unclassifed 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; 5.9%). Other genera with less than 6% abundances in all three habitats were also detected as shown in Figure 2B.

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: J 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, 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-sulfosulfanyl-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.

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 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 cysJ 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 genes for assimilatory reduction pathway are diversifying under moderate selection pressures.

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 (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).