The Microbial Composition in Circumneutral Thermal Springs from Chignahuapan, Puebla, Mexico Reveals the Presence of Particular Sulfur-Oxidizing Bacterial and Viral Communities

Terrestrial thermal springs are widely distributed globally, and these springs harbor a broad diversity of organisms of biotechnological interest. In Mexico, few studies exploring this kind of environment have been described. In this work, we explore the microbial community in Chignahuapan hot springs, which provides clues to understand these ecosystems’ diversity. We assessed the diversity of the microorganism communities in a hot spring environment with a metagenomic shotgun approach. Besides identifying similarities and differences with other ecosystems, we achieved a systematic comparison against 11 metagenomic samples from diverse localities. The Chignahuapan hot springs show a particular prevalence of sulfur-oxidizing bacteria from the genera Rhodococcus, Thermomonas, Thiomonas, Acinetobacter, Sulfurovum, and Bacillus, highlighting those that are different from other recovered bacterial populations in circumneutral hot springs environments around the world. The co-occurrence analysis of the bacteria and viruses in these environments revealed that within the Rhodococcus, Thiomonas, Thermonas, and Bacillus genera, the Chignahuapan samples have specific species of bacteria with a particular abundance, such as Rhodococcus erytropholis. The viruses in the circumneutral hot springs present bacteriophages within the order Caudovirales (Siphoviridae, Myoviridae, and Podoviridae), but the family of Herelleviridae was the most abundant in Chignahuapan samples. Furthermore, viral auxiliary metabolic genes were identified, many of which contribute mainly to the metabolism of cofactors and vitamins as well as carbohydrate metabolism. Nevertheless, the viruses and bacteria present in the circumneutral environments contribute to the sulfur cycle. This work represents an exhaustive characterization of a community structure in samples collected from hot springs in Mexico and opens opportunities to identify organisms of biotechnological interest.


Introduction
Terrestrial thermal springs are widely distributed throughout the world. They harbor a significant number of microorganisms of biotechnological interest. These ecosystems have been classified in low-temperature (<55 • C) and high-temperature (>55 • C) springs; in terms of pH, the springs are acidic (pH < 4), intermediate (pH~4), circumneutral or neutral (pH~7), or alkaline (pH > 7) [1][2][3]. Additionally, the thermal springs are classified according to their origin in magmatic waters, which are born in volcanic areas and at high temperatures (>50 • C), and telluric waters, which are formed when underground water currents pass along deep hot rocks [4].
Mexico contains a wide diversity of thermal springs, steam vents, geothermally heated soils, boiling mud pools, and geothermal zones [25,27,28]; however, few studies have described the diversity of microbial communities in thermal environments. In particular, in the acidic hot spring "Los Azufres" (pH 3.6 and 65 • C) located in the state of Michoacan, Rhodobacter, Acidithiobacillus, and Lysobacter [25], among other bacteria, have been identified; whereas the Sulfolobales archaeon has been discovered in "Los Azufres" [29]. Finally, the viruses identified correspond to archaeal Fusellovirus, archaeal Rudivirus, and Sulfolobales Archaeon AZ1 [30]. In addition, in a hot spring located in the Araro region, Michoacan, different genera of bacteria were found, such as Bacillus, Aeromonas, and Pseudomonas [31,32], whereas in the "Carrizal" thermal pool hot spring and in "Los Baños" in Veracruz, Mexico, bacteria of the genera Geobacillus, Anoxybacillus, and Aeribacillus have been identified [33].
In this work, we explore the microbial community and functional composition in the thermal spring "Baños Termales de Chignahuapan", which is located in the geothermal region of Tulancingo-Acoculco, Sierra Norte, Puebla. This is a mountainous complex, whose origin dates back to the Pleistocene (1.7-0.9 million years ago), where there is a magmatic hot spring with travertine sediment compositions and rocks of dacites and rhyolites [27,28].
We consider that the hot springs in Chignahuapan present a particular combination of physicochemical characteristics, such as high concentrations of calcium, carbonate, and sulfur, making it an excellent spot to determine the microorganisms and viruses that make up that ecosystem, as well as their functional potential. In addition, we achieved a comparative analysis with 11 circumneutral hot springs, to determine differences in composition and diversity in microorganisms and highlight the influence of the environment in the community structure.

Sample Collection and Processing from Chignahuapan Puebla, Mexico
Two samples of water (20 L) were collected from the recreational center "Baños Termales de Chignahuapan" located at coordinates 19 • 50'30' N 97 • 59'41" W, at an altitude of 2136 m above sea level, during April 2019. The samples were obtained with sterilized tools in 1 L containers, from the water emerged in the mountain before the pools were supplied, and transported at room temperature to the laboratory of the Benemérita Universidad Autónoma de Puebla, where they were filtered through 0.22 µm Millipore filters. The DNA was obtained from the filters and isolated using ZymoBIOMICS DNA kits (MoBio, West Carlsbad, CA, USA). The DNA concentration was determined using a NanoDrop 1000 spectrophotometer (Thermo Scientific), and fluorometry was measured using a Qubit 4 fluorometer (Invitrogen). DNA was sequenced using the Illumina NextSeq 500 platform with the Nextera reagent kit V3.0 for a read length 2 × 75 bp at the Instituto de Biotecnología of Universidad Nacional Autónoma de México.

Taxonomic Annotation of Metagenome
We analyzed two metagenomic samples from Chignahuapan, Puebla, and 11 shotgun metagenomic samples retrieved from the Sequence Read Archive (SRA) database, from thermal environments with physicochemical characteristics of circumneutral hot spring and pH 7 (Table S1). The quality control of sequences and the removal of adapters was performed by using Trimmomatic v.036 [39] with a sliding window of 4 bp, an average quality per base of 30, and a minimum read length of 75 bp. Reads were assembled in contigs using MEGAHIT [40], under default parameters in paired-end mode, and contigs with a minimum length of 1000 bp were considered for further analysis.
Taxonomic assignments were performed with the software Kaiju v1.7.3 [41] against the nonredundant protein database v1.7 sourced from the National Center for Biotechnology Information (NCBI) databases using the maximum exact matches, and 11 as a minimum match length. Finally, the results were displayed with the library Pavian in R [42].
The composition of the prokaryotic communities was evaluated using statistical analyses in R [43]. The nonmetric multidimensional scaling (NMDS) plot was performed in Vegan v2.3-1, with the stress function, to determine the goodness.
The taxonomic profiles at the genus level were used to calculate the diversity indices from all data. Diverse alpha-diversity descriptors were obtained using the Phyloseq function in R [44].
The beta diversity was determined by Bray-Curtis dissimilarity, and the sampling effort was evaluated through the rarefaction curves using the Vegan library implemented in R [45].
The metagenomes were deposited in the Joint Genome Institute (JGI) Integrated Microbial Genomes and Microbiomes database, with accession number: Gs014786.

Co-Occurrence Network Analysis
The co-occurrence analysis was implemented using the igraph library [46] and bipartite library in R [47], implemented under development Virome Network Analysis (ViNA) [48]. In brief, we computed a table of incidences of the relevant bacteria and viruses at the species level. These tables indicate the presence or absence of each taxon in the metagenome. After that, the network displayed the taxon associations and locations, which was built using the Kamada-Kawai algorithm layout [49].

Identification and Annotation of Viral Genomes
For virus classification, two approaches were implemented. In the first approach, Virsorter [50] was used to determine the viral contigs, using the viral hallmark genes annotated as "main capsid protein", "portal", "large subunit of the terminase", "tail", and "envelope", among others. The entire contig was considered viral if more than 80% of predicted genes on a contig had a viral signal. This software finds new viruses with different confidence categories from 1 to 6, with 3 and 6 as the least confidence level.
The categories 1, 2, 4, and 5 were concatenated, and contigs were compared with BLASTn against the nonredundant database (nr), with the following parameters: -num_alignments 20, -num_descriptions 20, e-value 0.0001, -word_size 11. The results were visualized in MEGAN v5.10.6 considering the lowest common ancestor (LCA) method with the following parameters that reduce the rate of false positives and false negatives [51]: minimum support of 2; minimum score of 70; top percent of 10.
In the second approach, the viral contigs were recovered and the auxiliary metabolic genes (AMG) were obtained using the program VIBRANT [52]. This program is a hybrid machine-learning algorithm and similarity comparisons of protein sequences. It annotates the genes supporting metabolism and recovers the metabolic pathways where these genes are involved. We considered a minimum length of 1000 bp, with summary plots on and function virome off. The retrieved contigs were compared with BLASTn against the Viral RefSeq database and analyzed using the software MEGAN v5.10.6 with the same conditions [51].

Functional Annotation
To predict protein-coding genes in the assembled contigs, we used Prodigal v2.6.3 [53] with the metagenomic mode. The functional annotation was achieved using SUPERFOCUS [54], which contains the SEED database with an E-value of 0.0003 and 60% identity. The results were displayed in a heatmap using the ggplot2 library in R [55]. The metabolic pathways were displayed using MG-RAST server with the database KO [56].

Field Sampling and Physicochemical Characterization
The water samples were collected from two thermal springs at Chignahuapan, Puebla, Mexico. The first sample was collected from the thermal spring (Mex_Chig_S1) with a temperature of 49-50 • C, and a pH of 7.02. The second sample was collected from the water that supplies the pool (Mex_Chig_S2). This sample temperature was 45 • C and the pH was 6.66.
The compositions of both samples were compared with respect to the content of different salts, quantified in Table 1. The sulfate was found with a concentration of 25.6 and 30.2 mg L −1 in Mex_Chig_S1 and Mex_Chig_S2, respectively; in the site, there was an intense smell of hydrogen sulfide (H 2 S), indicating that the sulfate was being reduced to hydrogen sulfide by sulfate-reducing microorganisms. These concentrations were in the range of thermal spring waters. In contrast, when Ca 2+ , carbonate (HCO 3 −1 ), and Na 2+ were quantified, these ions were present in high concentrations, suggesting that they are out of range according to water quality standards in Mexico. The presence of high levels of carbonates and calcium in both samples collected could be associated with hydrothermal travertine deposits found in hot springs, as these deposits are mainly composed of CaCO 3 (calcite) [6,19]. Likewise, there are previous reports of the presence of travertine deposits (calcium carbonate), rhyolite, and dacite in the Chignahuapan springs [28], and the presence of these carbonates could be involved in the modification of the microbial structures of the communities. The structure of the microbial community is also modified by the concentration of Na + ; it is known that decreased concentrations of salts lead to a higher diversity of bacteria, while archaea are abundant in a higher concentration of salts [57,58]. In the thermal environments of Puebla, a concentration outside the range for consumption according to the Mexican regulations was also found; thus, the high salt concentrations found in both samples will determine the diversity and structure of microbiomes in thermal springs.

Microbial Community Composition of Chignahuapan Metagenomes
The diversity and abundance of microorganisms in the circumneutral thermal spring samples from Chignahuapan, Puebla, were determined by shotgun metagenomic sequencing. From these metagenomes obtained from two locations, the diversity of the microbiome was determined. To this end, the assembled contigs were classified with Kaiju (Table 2), and the results at the domain level showed that Bacteria represent between 88.4 and 91.8% of the total microorganisms, followed by Archaea (1.3-1.8%), Viruses (0.8-0.9%), and Eukarya (0.7-0.8%) ( Figure S1). The microbial composition in the metagenomes was in concordance with that found in similar environments with a neutral pH, moderate temperature (50 • C), and similar chemical composition, such as the Jordanian hot springs [17,21] (Figure 1). In the sample Mex_Chig_S1, bacteria from the phylum Actinobacteria (64.21%) were the most abundant, followed by Proteobacteria (36.09%). In contrast, in the sample Mex_Chig_S2, Proteobacteria was found as the most abundant (76.81%), followed by Actinobacteria (20.14%) ( Figure 1). This result is consistent with previous descriptions in circumneutral water, sulfur water springs, and volcanic terrain, respectively, where Proteobacteria and Actinobacteria are predominant in sulfur water springs and volcanic terrain, respectively [17,48].  In Mex_Chig_S1, the most abundant genera were Rhodococcus (59%), followed by Acinetobacter (13%), Azotobacter (5%), Halothiobacillus (1.2%), and Bacillus (0.9%). At the species level, Rhodococcus erythropolis was identified as the most abundant ( Figure 1). It is worth mentioning that R. erythropolis, originally isolated from crude oil [59], is a biologically important bacterium because it possesses selective desulfurization activity and the capacity to degrade alkanes (C 8 to C 20 n-alkanes) and methyl benzenes such as toluene [59,60].
Another interesting genus was Halothiobacillus, which is an obligate chemolithoautotroph and sulfur oxidizer. In particular, Halothiobacillus neapolitanus was found in the samples that encode a complete SOX complex, involved in sulfur oxidation [61].
At the species level, Bacillus licheniformis (0.8%) was also found-an interesting bacterium which was also isolated from Jordanian hot springs [17]. Bacillus licheniformis is widely distributed in thermal springs, and it is considered for commercial use as it has been used in the production of enzymes, antibiotics, and detergents, but some species of Bacillus are involved in carbon metabolism [62,63].
Overall, these bacteria are moderate thermophilic organisms isolated from hydrothermal springs with diverse enzymatic activities characterized as amylases, cellulases, and lectinases [65]. T. bhubaneswarensis and T. intermedia are widely distributed in hot springs from India, which are rich in arsenic and contain low levels of organic matter; these bacteria are sulfur and thiosulfate-oxidizing [66].
The sulfur-oxidizing bacteria were highly abundant in both samples, correlating with the concentrations of sulfate in water. These bacteria can oxidize sulfur compounds (thiosulfate, tetrathionate, sulfide, and polysulfide) to produce energy, which has been previously reported in high-temperature sulfidic hot springs [14].
In general, these bacteria belong to Actinobacteria, Gammaproteobacteria, Betaproteobacteria, and Epsilonproteobacteria phyla, chemoheterotrophs or chemolithoautotrophs in the microbial community, with the ability to use electrons from inorganic compounds as an energy source. Overall, many of them are sulfur-oxidizing bacteria.
In the circumneutral thermal spring samples from Chignahuapan, Puebla, archaeal organism abundance was low, similar to findings in other circumneutral hot springs [17,50]. This archaeal composition is similar to those found in Malaysia's circumneutral hot spring, where a low proportion of archaea was found [67,68]. The most abundant classes identified in our samples correspond to Halobacteria within the Euryarchaeota phylum. Within the few species found were Halolamina sediminis and Candidatus altiarchaeum sp. Both are interesting because the first is associated with hypersaline aquatic environments, which may be possible due to the high salt concentrations present in Chignahuapan, and the second is involved with carbon fixation and plays an essential role in biogeochemical cycles [69,70]. These results contrast with previous reports indicating that Crenarchaeota is most abundant in terrestrial thermophilic environments in acid hot springs with high temperatures [63,65].

Comparative Analysis and Ecological Indices
To determine whether the samples of circumneutral hot springs from Chignahuapan, Puebla, share similarities with other hot springs samples, we compared the microbial community with 11 circumneutral metagenomes from the SRA database. These were obtained from different sites in the world, selected based on physicochemical characteristics similar to those of the sample's Mex_Chig_S1 and Mex_Chig_S2 (Table S1). The aquatic thermal terrestrial metagenomes were selected with a nearly neutral pH. The comparison showed that the Puebla hot springs metagenomes did not have a similar bacterial composition to the other metagenomes. To analyze this, we carried out a nonmetric multidimensional scaling (nMDS) analysis. The nMDS analysis showed that the Chignahuapan samples were grouped together, while the other metagenomes were grouped according to their geographical area, with the stress of 0.03; the lower the stress value, the better the goodness of fit ( Figure 2A).   Likewise, differences were found in the composition at the genus level in the bacterial communities. The microbial community structure was particular in the Chignahuapan samples, mainly of the Rhodococcus genus, which was found in great abundance, followed by genera such as Acinetobacter, Thermomonas, and Azotobacter ( Figure 2B).
This particular community of microorganisms is possibly associated with the physicochemical characteristics, including abiotic factors such as Ca, SO 4 , HCO 3 ion levels, that have previously been reported to have an essential role in the microbial composition and have been associated with some bacteria, such as Thermomonas hydrotermalis, Bacillus licheniformis, Bacillus subtilis, and Anoxybacillus kamchatkensis [6]. The composition of the bacteria and minerals reported is similar to those of Chignahuapan; therefore, the water's physicochemical factors very likely influence the structure of the microbial community.
Concerning the other metagenomes, the IN_SRR3050168_50 and IN_SRR8613699_52 from India were clustered near the Mexican samples, and the others, IN_SRR3961733_52, IN_SRR3961734_55, and IN_SRR3961739_43 samples were separated from all others ( Figure 2A). In this regard, the genera Acidovorax, Microbacterium, Pseudomonas, and Caulobacter were associated with India's location. A similar finding was observed with the samples from New Zealand, ZNL_SRR10063242_27.2, ZNL_SRR10063241_30, and ZNL_SRR10063240_34.5, where Streptomyces, Phenylobacterium, and Asticcacaulis were clustered. The hot springs of Canada (CAN_SRR10095342_45) and Japan (JP_SRR7905022_60) were clustered together, whereas the genera, Streptomyces, Roseiflexus, Desulfobacca, and Chloroflexus, were clustered ( Figure 2B and Figure S3).
The results show differences in the recovered genera and that they are grouped by geographic area. However, this grouping can also be driven by the composition of ions, minerals, and elements present in the water. In previous studies, measurements of different elements, compounds, and ions were taken in the water samples. For example, in India's selected samples, which have high concentrations of Co, La, Fe, Hg, and Si, the predominant bacteria were Pseudomonas stutzeri and Acidovorax sp.; these findings are in accordance with the taxonomic assignment made by us [9]. Furthermore, samples from another study from the same country were nearly clustered with them. However, in this other site, different dissolved solids were measured: there were high concentrations of phosphorus and sulfur ( Figure S3), and the genera Microbacterium, Propionibacterium, Caulobacter, and Rhodococcus predominated, among others [71].
Whereas the New Zealand samples had different values in the parameters evaluated, one of the samples had a high concentration of methane (ZNL_SRR10063242_27.2), and ammonia (ZNL_SRR10063240_34.5) and iron were present in all the samples. The genera that predominated in these metagenomes were different than those found in the other metagenomes [72].
On the other hand, in the Japanese metagenome (JP_SRR7905022_60), high levels of iron and dissolved oxygen have been reported, and these findings correlate with the more abundant bacteria genus Chloroflexus, which has photosynthetic activity [73].
However, the differences observed in the microbial community are associated with the temperature presented by thermal environments. For example, a study of the hot springs in Canada and New Zealand showed that some phyla had trends that changed with temperature, where Cyanobacteria, Acidobacteria, Verrucomicrobia, and Planctomycetes were absent at high temperatures, while other phyla did not show changes [74]. In our cluster analysis, we can also observe that there is certain proximity of the metagenomes that share a similar temperature ( Figure 2).
The multivariate approaches performed revealed that diversity patterns changed in each geographical location, where specific genera predominated in each of the metagenomes, and possibly this predominance was a consequence of the different physicochemical compositions and temperatures present in the water [26].
These differences can also be clearly distinguished in the relative abundance analysis. Proteobacteria were predominant, correlating with other moderate temperature circumneutral springs, such as the Jordanian Main springs, where about 50% corresponds to these phyla [17]. However, Puebla's microbial community structure was mainly compared to other terrestrial hot springs ( Figure S3). Overall, these results indicate that the microorganism communities change depending on the geographical area and physicochemical composition ( Figure S3).
The circumneutral hot springs were evaluated through the ecological indices or diversity indices, and species accumulation curves were compared among all the metagenomes. The rarefaction curve associated with the Chignahuapan metagenomes showed a low number of species, probably due to a low yield obtained from the sequencing run; therefore, these samples did not have asymptotic behavior ( Figure 3A). When the alpha and beta biodiversity indices were analyzed, a variation was shown along with the α-index. spe.dis Alpha-diversity values were lower for the Chignahuapan samples compared to the other metagenomes. However, it has been reported that diversity decreases with increasing temperature; these studies showed that the alpha-index is higher in those places with a neutral pH and a temperature of 50 • C [74]. Therefore, it was expected that the samples from Mexico had high diversity indices compared to other samples. In the case of the samples analyzed, the sample from Japan had the highest alpha index, compared to the other metagenomes. Therefore, it can be suggested that while diversity is modified by physicochemical changes, it is also modified by experimental yields and conditions.
The low values obtained for the β-diversity suggest that Chignahuapan metagenomes are less diverse than the other metagenomes ( Figure 3B). In this regard, the clustering using β-diversity showed that the hot spring metagenomes were clustered primarily based on geographic locations and temperature. The Chignahuapan thermal springs clustered independently, indicating that they have a different microorganism community than the other metagenomes ( Figure 3C). Interestingly, in this analysis, the sample from India was grouped with the samples from Mexico. The Indian sample also had low diversity indices indicating less diversity.
In summary, comparative analysis shows that Chignahuapan metagenomes are less diverse than other hydrothermal environments. However, they have a unique microorganism composition, because they do not group with the other samples with similar temperatures. Therefore, the chemical composition of water could determine the microbial structure of Chignahuapan hot springs.

Co-Occurrence Network Analysis
We performed a co-occurrence network analysis to evaluate possible genera interactions between the 13 metagenomes. Each of them evaluated the association between species and metagenomes; the genus Rhodococcus (Actinobacteria) was considered because it is the most abundant bacterial genus in Chignahuapan samples. Another of the abundant phyla was Proteobacteria. However, two genera of sulfur-oxidizing Beta-proteobacterium were considered in this analysis (Thiomonas and Thiobacillus) and Bacillus (Firmicutes). Although these genera were not the most abundant, they are sulfur-oxidizing bacteria, and so it is interesting to observe their interaction; the main network metrics were evaluated, which were connectance, nestedness, modularity index, and weighted closeness.
The Rhodococcus network showed a connectivity value of 0.54, weighted nestedness of 0.73, and modularity index of 0.39. The higher weighted closeness in the network was for Rhodococcus erythropolis (node A), with a value of 0.32. This result suggests that R. erythropolis is present in all circumneutral host springs; however, the value of weighted nestedness was closer for the Ching_S1_Mex and Ching_S2_Mex metagenomes, indicating the main predominance of this bacterium in metagenomes of Chignahuapan in comparison with other metagenomes. These results also confirm that R. erythropolis is the most abundant species ( Figure 4A).
The Thiomonas cluster showed the following network metric values: connectance 0.66; modularity index 0.10; weighted nestedness 0.65 (node B); Thiomonas intermedia and Thiomonas sp. FB-6 had the highest value-weighted closeness with 0.13 and 0.19, respectively. It was interesting that it was closer to most of the nodes in the graph or present in all terrestrial hot springs. In contrast, Thiomonas sp. CB6 (node H), Thiomonas sp. ACO7 (node I), and Thiomonas sp. B1 (node J) had the lowest value for weighted closeness, with 0.0067. Overall, in the Mexico metagenomes, Thiomonas species were less shared with the other metagenomes, indicating the microbial composition is particular ( Figure 4D   The Bacillus cluster showed the following network metric values: connectance, 0.51; modularity index, 0.471; weighted nestedness, 0.47. Furthermore, the value-weighted closeness of most higher values was 0.13 for to Bacillus cytotoxicus (node U), Bacillus subtilis (node A) with 0.09, and Bacillus cereus (node C) with 0.09; these showed higher centrality or connections within the network in all metagenomes. The Chignahuapan samples were closer to Bacillus licheniformis (node B) at 0.06, whereas the lower weighted closeness was B. cereus R309803 (node D) at 0.0009, Bacillus sp. OxB-1 (node G) at 0.003, B. amyloliquefaciens (node M) at 0.006, B. testis (node J) at 0.001, and Bacillus thuringiensis serovar israelensis ATCC 35646 (node P) at 0.0009, which indicates that these species of bacteria are poorly connected in the network and are unique within the metagenomes Chig_S2_Mex (M1) and ZNL_SRR10063240 (M3). Interestingly, Bacillus cereus R309803 is a unique species in the Chignahuapan metagenome ( Figure 4B). The value of modularity suggests that the network of Bacillus has a modular structure in this cluster. Modularity with values above 0.44 indicates that the networks are more connected.

Functional Metagenomics Analysis
A functional analysis was performed with all metagenomes to determine whether the functional activity was similar in all the metagenomes; the amino acid sequences were annotated using SEED subsystems. The SEED subsystem is a classification system that organizes the coding sequences for functional categories into a hierarchy with 5 levels of resolution; in level 1, the families of proteins that share function. The results shown in Figure 5 correspond to the level 2 families.
The metabolism of carbohydrates (~13-22%) (central carbohydrate metabolism, CO 2 fixation, and fermentation) are the hot springs' main functional processes. From these, the metabolic pathways identified in high abundance correspond to the Calvin-Benson cycle, CO 2 uptake carboxysome, Tricarboxylic acid (TCA) cycle, and oxidative phosphorylation pathways involved in CO 2 fixation. Carbon fixation is a process where inorganic carbon (in the form of CO 2 ) is transformed into organic compounds and is an essential process for the production of anabolism precursors. Four metabolic pathways in bacteria have the capacity of Carbon fixation: Calvin-Benson cycle, the reverse TCA cycle, the Wood-Ljungdahl pathway, and the 3-hydroxypropionate (3-HP) bicycle [75,76]. Two carbon fixation pathways (Calvin-Benson cycle and TCA cycle) were identified in the metagenomes, whereas only the gene encoding the ribulose bisphosphate carboxylase (RuBisCO) was found in the samples. The results suggest that the Calvin-Benson cycle is the primary metabolic pathway in the Chignahuapan terrestrial hot springs water ecosystem ( Figure S4).
Carbonate was found in high concentrations in the hot water spring of Chignahuapan; the high abundance of HCO 3 could correlate with the presence of pathways involved in the fixation of CO 2 , and the carbonic anhydrase, involved in the interconversion of CO 2 to carbonate [77]. Additionally, the CO 2 is dissolved in water and can form different compounds such as carbon dioxide, carbonic acid, bicarbonate, and carbonate, and sulfur-oxidizing bacteria can fix it. In this context, Acidithiobacillus ferrooxidans was found in metagenomes and uses CO 2 for growth [78]. Furthermore, other thermophilic bacteria have been reported, including Bacillus schlegelii and P. thermocarboxydovorans growing with CO 2 as a unique carbon and energy source [79,80]. At the first level within the SEED category, the amino acids and derivatives, cofactors, vitamins, prosthetic groups, and pigments were the most abundant categories (~2-14%), followed by lysine, threonine, methionine, and cysteine biosynthesis. Most of these amino acids involved in protein synthesis and as cofactors in many metabolic reactions; however, methionine and cysteine biosynthesis are involved in the sulfur metabolism for biosynthesis of sulfur-containing biomolecules, such as the cysteine sulfinate-dependent pathways to produce sulfite or turine [81].
This metabolic category has also been reported in abundance in the hot springs of the Araró region, located in the Trans-Mexican Volcanic Belt [32], suggesting a central role in the biosynthesis of diverse compounds, for assimilatory sulfate reduction, SOX pathway, and for obtaining energy in thermal water environments ( Figure 5).
Sulfur metabolism has been reported in microbial mats, hydrothermal vents, and YNP hot springs, where microorganisms called sulfur oxidizers live, as well as many chemolithotrophic Proteobacteria [77,82]. Overall, in the metagenomes analyzed, the proportion of sulfur metabolism was present in the metagenomes (0.9-1.9%).
Chignahuapan hot springs contain the complete SOX pathway. This pathway is involved in the oxidation of sulfide (S 2− ) and thiosulfate (S 2 O 3 2− ) to sulfate (SO 4 2− ) ( Figure S5). This pathway has been reported present in alpha and epsilon-proteobacteria [18,59,83], identified in the two metagenomes of Puebla, suggesting that epsilon-proteobacteria are contributing to sulfide and thiosulfate oxidation, as an energy source. It has been reported that some bacteria can remove inorganic sulfur from oil, and it has been investigated that the enzymes to carry out this process are within the sox pathway; this was characterized in Rhodococcus sp. strain IGTS8 [84], the most abundant genera identified in the hot springs. The bacteria H. neapolitanus and Acidithiobacillus caldus have sulfur-oxidizing enzyme systems involved in the SOX pathway, and both bacteria were present in Chignahuapan.
In summary, sulfur oxidation is carried out by bacteria within the phylum Proteobacteria and Actinobacteria in Chignahuapan thermal springs. Where the assimilatory sulfate pathway involved in the reduction in sulfate (SO 4 2− ) to sulfide (S 2− ) was complete ( Figure S5) [85,86]. These results contrast with bacteria involved in sulfur metabolism such as Deltaproteobacteria and Firmicutes found in Yellowstone National Park [87,88]. Nitrogen fixation is associated with carbon fixation in the microbial mat communities, and the nitrogen fixation occurred in many environments [89]. Nitrogen-fixing enzymes could be expected to be present in metagenomes. However, in general, the nitrogen fixation pathway enzymes found to be abundant in the samples of Chignahuapan were the enzymes involved with the pathway assimilatory nitrate reduction nitrate to ammonia (nasA and nirB genes) and ammonia-lyases also involved in ammonia production. As mentioned above, this pathway is associated with nitrogen fixation, and carbonic anhydrase was found here, which is associated with converting carbon dioxide to carbonates ( Figure S6).
The stress oxidative and heat shock category was highly abundant in Chignahuapan of the reactive oxygen species (ROS) and can cause irreversible damage to cells, and indifferent thermophilic bacteria, and it has been reported to present a superoxide dismutase; this enzyme catalyzes the dismutation of the superoxide (O 2 1− ) into oxygen (O 2 ). In the case of the samples, the bacteria that presented the putative enzyme were Thioalkalivibrio spp., and Acinetobacter spp., while Rhodococcus opacus B4 had NrdH, which mediates resistance to oxidative stresses. Previously, bacteria with superoxide dismutases had been reported, such as Aquifex pyrophilus, Hydrogenobacter thermophilus, Thermus thermophilus, Propionibacterium shermanii, and Rhodothermus sp. among others. This system is essential to avoid damage to extreme environments [90]. This abundant category of functions also relates to other thermal environments, such as the mats from Araró Mexico [32].
Similarly, as there are mechanisms to prevent cell damage, there are DNA repair mechanisms; previously in hypersaline waters, we observed that microorganisms have different enzymes involved in DNA repair [58]. In the case of thermophilic microorganisms, there are efficient mechanisms to prevent DNA repair and proteins and the lipid membrane from preventing these damages: for example, the composition of fatty acids changes with increasing temperature [90]. The thermophilic microorganisms of Chignahuapan also have repair systems; however, they are abundantly presented with the enzyme exonuclease SbcC, which is involved in DNA repair when alkylation is damaged.
Overall, the annotation allowed us to predict the functional potential of the thermophilic community of circumneutral metagenomes. Whereas mainly fixation of carbon is the crucial pathway and amino acids and derivatives, many of them contribute to sulfur metabolism and fixation of carbon. Likewise, the pathway to oxidation and reduction in metabolites of sulfur was completed. Many of the thermophilic microorganisms have a mechanism to prevent oxidative stress damage and repair DNA damage.

Viral Community Composition
For the viral community in a circumneutral terrestrial hot spring, the analysis shows that overall, nine viral families were retrieved, and many of them infect bacteria, such as Siphoviridae, Myoviridae, Podoviridae, Corticoviridae, and Herelleviridae. These viruses also infect invertebrates Baculoviridae, eukaryotic algae Phycodnaviridae, and Protozoan Mimiviridae.
Clustering analysis revealed that most abundant families are Myoviridae and Siphoviridae, which are ubiquitous in all metagenomes. These samples also showed that the metagenomes of Mexico can be grouped with samples from India and Canada with similar temperatures (Figure 6). subtilis and B. licheniformis.
VIRSorter and VIBRANT programs were used in order to retrieve viral contigs. This software failed to recover the complete genome from samples of Chignahuapan, only partial sequences, probably because of the low performance in the metagenome sequencing. However, according to the taxonomic classification from viral contig retrieval, the most abundant viruses were Acidithiobacillus phage AcaML1, Bacillus phage SIOphi, Bacillus virus Bobb, Acinetobacter virus R3177, and Bacillus phage Shbh1 in the samples from Mexico, correlating with what was found in bacteria. According to the taxonomic assignment, 200 viruses were recovered with Vibrant and 125 viruses with Virsorter in all samples. In the Japanese metagenome, there was a greater abundance of viruses ( Figure S7).
Overall, the viral communities in moderate thermophilic environments that infected moderate to thermophilic bacteria demonstrated distinct viral community structures among the circumneutral thermal springs, compared with acid or hyperthermophile hot springs such as Yellowstone [92]. These results indicate that the circumneutral thermal springs harbor viral communities with phage double-stranded DNA that infect mainly bacteria, followed by viruses that infect invertebrates or eukaryotic algae.
Through an analysis of occurrence carried out with the recovered viruses from Vibrant, the results from the analysis revealed there are few connections between virus species in the metagenomes, each having many particular species. These results correlate with the classifications of bacteria, where each one observed that each environment has a genus of particular bacteria, and that terrestrial hot springs have a common Acidithiobacillus phage AcaML1, which infects Acidithiobacillus caldus ( Figure S8).
Acidithiobacillus caldus has been reported in the various thermal environments. It is a moderately thermoacidophilic bacteria that contributes to the carbon and sulfur cycles, as it obtains energy from the oxidation of elemental sulfur for carbon dioxide fixation, and it is ubiquitous in sulfide mineral environments [95]. These results contrast with previous reports in thermal or geothermal terrestrial environments; it has been determined that the main families of viruses present are Fuselloviridae, Bicaudaviridae, Turriviridae, Ampullaviridae, Guttaviridae, Lipothrixviridae, Rudiviridae, and Globuloviridae [3,91]. Most of these virus families infect Archaea that live in thermal environments with a temperature above 80 • C, and these viruses have been called Archeoviruses [92,93].
However, it was expected that the viruses recovered in these metagenomes mainly infect bacteria, since the temperature and pH are also moderate in the analyzed metagenomes. Thus, the dominance of viruses that infect prokaryotes found are moderate thermophilic bacteria. Nevertheless, in the hot springs of Chignahuapan, besides Siphoviridae and Myoviridae, the Herelleviriadae family was abundant.
The Herelleviriadae family was recently described: the phylogenetic evidence performed at the genome and proteome level, and by a single gene (capsid, tail protein, DnaB395 like-helicase), showed that these viruses are polyphyletic since they are grouped into five different clades or subgroups. This family contains linear viral genomes, with a length of 125-170 kbp, that infect the genus Bacillus [94].
The abundance of bacteriophages of the Herelleviriadae family correlates with the taxonomic assignment results, where the genus Bacillus was present in both metagenomes, and according to the co-occurrence analysis, the main associated species of the metagenomes of Mexico were mainly B. subtilis and B. licheniformis.
VIRSorter and VIBRANT programs were used in order to retrieve viral contigs. This software failed to recover the complete genome from samples of Chignahuapan, only partial sequences, probably because of the low performance in the metagenome sequencing. However, according to the taxonomic classification from viral contig retrieval, the most abundant viruses were Acidithiobacillus phage AcaML1, Bacillus phage SIOphi, Bacillus virus Bobb, Acinetobacter virus R3177, and Bacillus phage Shbh1 in the samples from Mexico, correlating with what was found in bacteria. According to the taxonomic assignment, 200 viruses were recovered with Vibrant and 125 viruses with Virsorter in all samples. In the Japanese metagenome, there was a greater abundance of viruses ( Figure S7).
Overall, the viral communities in moderate thermophilic environments that infected moderate to thermophilic bacteria demonstrated distinct viral community structures among the circumneutral thermal springs, compared with acid or hyperthermophile hot springs such as Yellowstone [92]. These results indicate that the circumneutral thermal springs harbor viral communities with phage double-stranded DNA that infect mainly bacteria, followed by viruses that infect invertebrates or eukaryotic algae.
Through an analysis of occurrence carried out with the recovered viruses from Vibrant, the results from the analysis revealed there are few connections between virus species in the metagenomes, each having many particular species. These results correlate with the classifications of bacteria, where each one observed that each environment has a genus of particular bacteria, and that terrestrial hot springs have a common Acidithiobacillus phage AcaML1, which infects Acidithiobacillus caldus ( Figure S8).
Acidithiobacillus caldus has been reported in the various thermal environments. It is a moderately thermoacidophilic bacteria that contributes to the carbon and sulfur cycles, as it obtains energy from the oxidation of elemental sulfur for carbon dioxide fixation, and it is ubiquitous in sulfide mineral environments [95].

Auxiliary Metabolic Genes and Whole Viral Genomes
We evaluated the presence of AMGs in the viral contigs recovered from all the metagenomes, but in the case of the Chignahuapan metagenomes, no AMGs were obtained. The viral genes recovered in these correspond mainly to structural parts of the virion. Thus, this comparison was carried out only with those viral contigs that contained AGMs.
The AMGs have been related to an increase in the fitness and altering or complementing of their host's metabolism, facilitating adaptation under adverse conditions [96][97][98], and they have been related to photosynthesis, carbon fixation [99], and sulfur and nitrogen biogeochemical cycles [100].
The most abundant AMGs were classified within the functional categories metabolism of cofactors and vitamins (MCV), carbohydrates metabolism (CM), amino acids metabolism (AAM), metabolism of terpenoids, and polyketides ( Figure 7).
In the MCV category, the virus contributes to the folate biosynthesis pathway (KEGG entries identified: K00287, K01495, K01737, K06920, K09457, K10026). Tetrahydrofolate is a cofactor present in all bacteria, and it is essential to the growth; synthesis of formylmethionyl tRNAfMet, carried out by the enzyme dihydrofolate reductase (DHFR); and in thermophilic bacteria, has been identified to use modified folates [101], and within the folate biosynthesis pathway, the AMG 7-cyano-7-deazaguanine synthase was also found, a hypermodified 7-deazaguanosine, which has been previously reported in viruses; it is proposed that phages have taken the 7-deazaguanine from bacteria to evade the restriction-modification system (RM system) [101,102]. Interestingly, viruses found in thermophilic environments have these AGMs.

Conclusions
The composition of microorganisms in Chignahuapan is driven by chemical composition and geographic location. Since the microbial community's structure was particular where the bacteria of the genera Rhodococcus, Acinetobacter, Thermomonas, Tepidimonas, and Azotobacter predominated, in comparison to other circumneutral environments, these bacteria are sulfur oxidizers, which is consistent with the functional analysis where the sulfur reduction and oxidation pathways were complete. The functional analysis also predicted that the Calvin-Benson cycle metabolic pathways are the main pathways to contribute to carbon fixation. Furthermore, the microorganisms present in circumneutral environments have mechanisms that prevent cellular and DNA damage. Therefore, the microbial community structure in particular in each location is driven by physicochemical properties, but many metabolic pathways were common in circumneutral terrestrial hot springs, Additionally, other metabolic pathways were found within this category MCV such as Thiamine metabolism (K03153, K04487), biotin metabolism (K00059, K09458), pantothenate and coenzyme (CoA) biosynthesis (K13038), and retinol metabolism (K11153), and within nicotinate and nicotinamide metabolism (K00858, K01916, K03462, K13522) viruses have a gene that could modify or complement these metabolic pathways.
The categories least abundant but not less interesting, are those for obtaining energy and degradation of hydrocarbons. Some viruses have been reported in oil reserves, which could have interesting functions with the biotechnological applications [103].
The energy category was found in a low proportion; the pathways involved were sulfur metabolism and methane metabolism. One of the most important enzymes that we found was cysH phosphoadenosine phosphosulfate reductase, which is involved in the reduction in 3 -phosphoadenosine-5 -phosphosulfate (PAPS) to sulfite, an intermediate step of the pathway assimilatory sulfate reduction. Additionally, sufS which encodes cysteine desulfurase was found: this enzyme mobilizes the sulfur from L-cysteine.
In previous work, it has been reported that the viruses inside hydrothermal vents carry out AMG related to dissimilatory sulfite; specifically, the enzyme reverse dissimilatory sulfite reductases encoded by rdsrA and rdsrC in charge of oxidizing elemental sulfur, and some oxidizing sulfur bacteria lack the Sox system and use rdsrA for oxidizing elemental sulfur [104]. Nevertheless, here, the AMGs found were cysH and sufS involved in sulfur metabolism. Therefore, viruses and bacteria present in circumneutral environments contribute to the sulfur cycle.
Finally, there are some bacteria that are capable of carrying out the degradation of hydrocarbons; interestingly, we found viral AMGs in low abundance involved in the degradation of xenobiotics: for example, in the Benzoate degradation pathway (K01055), pcaD genes encoding to 3-oxoadipate enol-lactonase. Toluene degradation and fluorobenzoate degradation pathway (K01061) were found through the hydrolase carboxymethylenebutenolidase. This represents a biotechnologically important finding since these sequences could be obtained for use in the industry. These results reveal that viruses contribute to carbohydrate metabolism, sulfur cycles, and the degradation of aromatic compounds.

Conclusions
The composition of microorganisms in Chignahuapan is driven by chemical composition and geographic location. Since the microbial community's structure was particular where the bacteria of the genera Rhodococcus, Acinetobacter, Thermomonas, Tepidimonas, and Azotobacter predominated, in comparison to other circumneutral environments, these bacteria are sulfur oxidizers, which is consistent with the functional analysis where the sulfur reduction and oxidation pathways were complete. The functional analysis also predicted that the Calvin-Benson cycle metabolic pathways are the main pathways to contribute to carbon fixation. Furthermore, the microorganisms present in circumneutral environments have mechanisms that prevent cellular and DNA damage. Therefore, the microbial community structure in particular in each location is driven by physicochemical properties, but many metabolic pathways were common in circumneutral terrestrial hot springs, which contribute to carbon fixation. In the hot springs' viral community, prokaryotic viruses predominate overall, but in Chignahuapan, the Herelleviridae family was mainly abundant. The analysis of the auxiliary genes revealed that the viruses also contribute to metabolic pathways and the sulfur cycle. This study shows the information on the microbial and viral diversity in Mexican hot springs and compares these microbial communities with other metagenomic samples, thus providing an opportunity to understand the role of the viral AMGs and the structure of the viruses in the adaptation process to circumneutral hot springs. Additionally, this report serves as a reference to viromes in extremophile environments from Mexico.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2076-2607/8/11/1677/s1, Figure S1: Taxonomic profiling and community structure comparison at the domain level. Our profiling showed a dominance of Bacteria, at around 88.4 to 91.8%, followed by Archaea, viruses, and eukaryotes; Figure