Distribution and Activity of Sulfur-Metabolizing Bacteria along the Temperature Gradient in Phototrophic Mats of the Chilean Hot Spring Porcelana

In terrestrial hot springs, some members of the microbial mat community utilize sulfur chemical species for reduction and oxidization metabolism. In this study, the diversity and activity of sulfur-metabolizing bacteria were evaluated along a temperature gradient (48–69 °C) in non-acidic phototrophic mats of the Porcelana hot spring (Northern Patagonia, Chile) using complementary meta-omic methodologies and specific amplification of the aprA (APS reductase) and soxB (thiosulfohydrolase) genes. Overall, the key players in sulfur metabolism varied mostly in abundance along the temperature gradient, which is relevant for evaluating the possible implications of microorganisms associated with sulfur cycling under the current global climate change scenario. Our results strongly suggest that sulfate reduction occurs throughout the whole temperature gradient, being supported by different taxa depending on temperature. Assimilative sulfate reduction is the most relevant pathway in terms of taxonomic abundance and activity, whereas the sulfur-oxidizing system (Sox) is likely to be more diverse at low rather than at high temperatures. Members of the phylum Chloroflexota showed higher sulfur cycle-related transcriptional activity at 66 °C, with a potential contribution to sulfate reduction and oxidation to thiosulfate. In contrast, at the lowest temperature (48 °C), Burkholderiales and Acetobacterales (both Pseudomonadota, also known as Proteobacteria) showed a higher contribution to dissimilative sulfate reduction/oxidation as well as to thiosulfate metabolism. Cyanobacteriota and Planctomycetota were especially active in assimilatory sulfate reduction. Analysis of the aprA and soxB genes pointed to members of the order Burkholderiales (Gammaproteobacteria) as the most dominant and active along the temperature gradient for these genes. Changes in the diversity and activity of different sulfur-metabolizing bacteria in photoautotrophic microbial mats along a temperature gradient revealed their important role in hot spring environments, especially the main primary producers (Chloroflexota/Cyanobacteriota) and diazotrophs (Cyanobacteriota), showing that carbon, nitrogen, and sulfur cycles are highly linked in these extreme systems.


Introduction
Sulfur is one of the six most abundant chemical elements on planet Earth [1]. Its terrestrial biogeochemical cycle is influenced by volcanic activity, decomposition of organic matter, and fossil fuel combustion [2]. In terrestrial hot springs, solubilization and surface release of metallic sulfides occur due to increased hydrostatic pressure when groundwater is heated in the proximity of magma [2,3]. Subsurface geochemistry highly influences the concentration of sulfur compounds, which can be found at high concentrations in acidic sulfur springs (pH < 3), but at low sulfate concentrations (<100 mg/L) in chloride hot springs with more neutral pH [4]. These sulfur species can be used as electron acceptors or donors by a variety of sulfur-metabolizing microorganisms [5].
Microbial communities of terrestrial hot springs can develop tight mats stratified longitudinally by physicochemical gradients generated in part by their metabolism [6,7]. Temperature and sulfur concentrations are major drivers of the structure of these microbial communities [8][9][10]. In the upper layers of non-acidic hot spring mats, oxygenic phototrophic Cyanobacteriota (former Cyanobacteria) and anoxygenic phototrophic Chloroflexota members dominate, assimilating and utilizing sulfur as part of their metabolism [7,[11][12][13][14][15]. In the deeper layers of the mat, anaerobic microorganisms capable of dissimilative sulfate reduction metabolism thrive [16,17].
Several metabolic pathways are involved in the sulfur cycle [16], with assimilative sulfur reduction (ASR), dissimilative sulfur reduction (DSR), and sulfur oxidation (SOX) being the most studied. ASR is mostly related to the assimilation of sulfur as a substrate to synthesize biological molecules [16], whereas DSR and SOX are mainly associated with energy processes that use sulfur compounds during electron transfer-associated processes [16]. Since the rise of the omics era, next-generation sequencing has been used to study biological diversity related to the sulfur cycle [18,19], helping to understand new aspects of this cycle in the environment.
In hot springs, the diversity of sulfur-reducing bacteria (SRB) and sulfur-oxidizing bacteria (SOB) has previously been studied based on the 16S rRNA gene [20][21][22], as well as the functional marker genes aprA, soxB, and dsrA [9,[23][24][25]. The aprA gene encodes an adenosine-5 -phosphosulfate (APS) reductase that is highly conserved among SRB and is also present in several SOBs that utilize the so-called reverse sulfate reduction pathway [26]. During the sulfate dissimilative reduction process, phosphorylated APS, which is generated by the heterodimeric enzyme ATP sulfurylase, is directly reduced to sulfite by APS reductase. Thus, the aprA gene has been useful for simultaneously determining the diversity of both SRB and SOB microbial groups [27,28]. Indeed, this has allowed the detection of a variety of SRB, such as the genus Thermodesulfobacterium [24]; photolitotrophic sulfur oxidizers and strict anaerobes of the phylum Chlorobiota; and chemolithotrophs of the classes Alpha-, Beta-, and Gamma-proteobacteria found in both terrestrial and submarine hot springs [29,30].
Conversely, the soxB gene, encoding soxB thiosulfohydrolase, is the most conserved component of the sulfur-oxidizing Sox system, and catalyzes dissimilative oxidation of thiosulfate and other sulfur compounds [31]. The Sox system has been described as an alternative oxidation pathway to the reverse sulfate reduction pathway and has been reported in several SOB, including Chlorobiota, Pseudomonadota (Alpha-and Gamma-proteobacteria), and Acidithiobacillia [32]. In environments such as hydrothermal vents and terrestrial hot springs, the soxB gene has been used as a complimentary marker to the aprA gene to determine SOB diversity [28,30,[32][33][34].
However, to date, only a few studies have described the microbial community variability related to sulfur cycling along the temperature gradient of terrestrial hot springs [9,22,35]. In addition, the diversity and transcriptional activity of microbial communities have been overlooked in non-acidic thermal systems. Thus, the objective of this work was to address the knowledge gap related to the variation and activity of sulfur cycle-related organisms along a temperature gradient in non-acidic hot springs. Accordingly, we used a complementary combination of methodologies to analyze the genomic taxonomy, diversity, distribution, and transcriptional activity of the sulfur-metabolizing bacterial community in microbial mats from the Porcelana hot spring (northern Chilean Patagonia). The main results reveal the diverse functions of sulfur-cycling microorganisms in these thermophilic microbial communities, many of which have been previously described for their particular roles in the carbon and nitrogen cycles [36,37].  [37,38]. Temperature, conductivity, dissolved oxygen, and pH were monitored in situ at all sampling sites using a multiparameter instrument (Oakton, Des Plaines, IL, USA; model 35607-85) [37,38]. Water samples (15 mL) were collected in March 2013 to determine the inorganic sulfur species. To inhibit the microbial and inorganic oxidation of hydrogen sulfide, water samples were stored with a 2 M zinc acetate solution (Winkler ™, Santiago de Chile, Chile), and the pH was adjusted to ≥9 using a 6 M NaOH solution.

Determination of Inorganic Sulfur Species in Porcelana Thermal Water
For the quantification of inorganic sulfur compounds, the protocol described by the American Public Health Association was followed [39]. The sulfide calibration curve (five points from 0.05 to 0.5 mg/L S −2 ) was prepared from a serial dilution of a standard 10 mg/mL Na 2 S solution (Sigma-Aldrich™, St. Louis, MO, USA). Sulfide was measured spectrophotometrically at 664 nm using the methylene blue method [40]. The sulfate calibration curve (five points from 5 to 40 mg/L) was prepared using a 1 g/L stock solution of anhydrous Na 2 SO 4 . Sulfate was measured by turbidimetry analysis using spectrophotometry at 420 nm [41]. The thiosulfate calibration curve (five points from 0.5 to 4 mg/L) was prepared using a standard 0.01 M sodium thiosulfate pentahydrate solution. Thiosulfate levels were determined by measuring absorbance at 350 nm [42].  [37,38]. PCR amplification of the aprA and soxB genes was performed using the primer sets aprA-1-FW/aprA-5-RV [43] and soxB-693-F/soxB-1446-R [34], respectively. Primer sequences and final PCR conditions using GoTaq ® DNA Polymerase (Promega, Madison, WI, USA) are provided in Table S1. PCR products were purified using the Wizard ® Clean-up System kit (Promega™) and cloned with the pJET ® 1.2/blunt cloning kit (Thermo Scientific, Waltham, MA, USA). Plasmids were transformed into chemocompetent Escherichia coli DH5α following the manufacturer's protocol. Escherichia coli DH5α colonies were verified by PCR using the same primers as noted above. Verified clones were subjected to a second colony-PCR using the pJetF/pJetR primers (Thermo Scientific). From 100 clones per sample, 20% of the PCR products of the expected size were sent for sequencing (Macrogen Inc., Seoul, South Korea). Identification and assignment for a given operational taxonomic unit (OTU) were performed for the aprA and soxB clones by comparing restriction fragment length polymorphism patterns using NlaIV and SduI restriction enzymes. These enzymes were selected based on a cutting simulation using A Plasmid Editor (v2.0.45) software performed on randomly selected aprA and soxB gene clone sequences, which resulted in specific patterns for both genes. Enzymatic digestion of the PCR products of the clones from each clone library was performed according to the manufacturer's instructions (Thermo Scientific). Eighty percent of the aprA and soxB sequences were identified by restriction patterns.
The protein sequences of the three assembled Porcelana metagenomes [44] were predicted using the software Prodigal [54]. Taxonomy was assigned based on that provided by GTDB-tk (v0.3.2) R95 software [53] for proteins belonging to the MAGs. In addition, the taxonomic identity of unbinned proteins was determined by the closest match against predicted proteins from the GTDB R95 database [51] using DIAMOND software (blastp, e-value = 1 × 10 −6 ) [55]. Functional annotation was performed using eggNOG-mapper software [56,57].
Trimmed reads obtained from the three metagenomes and metatranscriptomes (48, 58, and 66 • C) from samples obtained in 2013 [36] were aligned to the gene set from 67 MAGs or from unbinned contigs coming from the metagenomes using bowtie2 (v2.3.4.1) (--end-toend--sensitive-k) [58]. The per-gene count table was processed in R using the Tidyverse package (https://www.tidyverse.org/, accessed in March 2021). The analysis focused on three main pathways of the sulfur cycle: (i) assimilative sulfate reduction (ASR), which is involved in sulfur assimilation and includes genes of the cys family; (ii) dissimilative sulfate reduction, which uses sulfur species as electron acceptors and includes the dsr and apr genes; and finally, (iii) the Sox system, which comprises the Sox family of proteins involved in thiosulfate metabolism. To prevent overcounting of KEGG orthologs (Kos), the number of reads associated with each gene was divided by the number of Kos associated with the respective gene, then equality was distributed among the respective Kos [59], and finally normalized by quantiles. Kos related to the sulfur cycle were obtained from the KEGG website (pathway ID: KO00920). Normalized sulfur counts per function were summarized by KO per taxa in each sample, while genes were grouped and summed according to the sulfur-cycle step in which they were involved. Heatmaps were generated using the ggplot package in R.

Taxonomic Description and Phylogenetic Reconstruction
Sequences corresponding to the aprA (TIGRFAM accession: TIGR02061) and soxB (TIGRFAM accession: TIGR04486) genes were extracted using the hmmsearch tool of the HMMER (v3.2.1) software (http://hmmer.org/, accessed in November 2018) from the gene set predicted from the MAGs and unbinned contigs assembled from the three Porcelana metagenomes. The obtained sequences were checked against the nr and refseq_protein databases to confirm annotations, avoiding the recruitment of other protein sequences with shared domains. Taxonomic assignment of the soxB and aprA genes from the cloned OTUs produced in this study and the three metagenomic samples was performed by the closest BLASTP match against the GTDB R95 custom protein database [51]. The aprA and soxB gene sequences from clones and metagenomes, as well as their closest references, were then aligned using MUSCLE (v6.0.98) software [60]. Maximum-likelihood trees were generated with the Iqtree (v.1.5.5) software, using the TESTNEW option to choose the best substitution model and the non-parametric ultrafast bootstrap option with 10,000 replicates [61]. Phylogenetic reconstructions, node collapse, and tree rooting were managed with the iTOL web server [62].

Changes in Sulfur Compounds along the Temperature Gradient at Porcelana Hot Spring
The sulfide, sulfate, and thiosulfate levels of samples collected in 2013 at the Porcelana hot spring (pH = 6 to 7; temperature = 48 to 69.7 • C) were measured spectrophotometrically (See Material and Methods). Table S2 summarizes these results and also includes data retrieved from previous studies [37,38]. The low sulfate levels (i.e., <100 mg/L according to [4]) observed in Porcelana (56.02 to 66.24 mg/L) at 48, 58, and 66 • C were similar to those previously reported in other non-acidic or neutral chloride hot springs from the Tibetan Plateau [9], Japan [22,63], Iceland, and Yellowstone National Park, USA [17,64,65]. Sulfate levels decreased with distance from the source, specifically from 66.24 mg/L at 66 • C to 57.82 mg/L at 58 • C and 56.02 mg/L at 48 • C. These levels were similar to those previously reported at Porcelana (between 40.6 and 66.9 mg/L, [37]), suggesting a constant presence of sulfate over time. Sulfide levels remained between 1.13 and 1.51 mg/L, while thiosulfate increased from 4.42 to 8.09 mg/L with decreasing temperature (Table S2), resembling values reported for other hot springs [64,66]. Changes in sulfur species concentrations with temperature, such as lower sulfate levels at lower temperatures, suggest dynamic biological activity related to the sulfur cycle [37,38].

Key Active Sulfur-Metabolizing Organisms and Pathways in the Sulfur Cycle of Porcelana Revealed by Meta-Omics
Since organisms potentially related to the sulfur cycle were found in Porcelana microbial mats [37,38], we analyzed the already available meta-omic data from the temperature gradient at 48, 58, and 66 • C [36,44] to obtain a general but detailed overview of the presence and activity of microorganisms involved in the sulfur cycle.
The sulfur-cycle-related gene abundances were quantified by read alignment to genes predicted from MAGs and unbinned contigs and then normalized and processed as described in Materials and Methods. Individual gene counts were grouped according to their order and phylum and associated pathway (pathways and KO IDs shown in Table S3), revealing the presence of the pathway ( Figure 1) and its transcriptional levels ( Figure 2). ASR and DSR reactions are mentioned in the reductive way for formality, as some of these enzymatic reactions can function also in an oxidative way [67]. The main results show that the ASR pathway was the most abundant and active sulfite-to-sulfide reduction pathway at 48, 58, and 66 • C ( Figures 1A and 2A). This pathway was mainly associated with dominant primary producers in Porcelana, such as members of the order Cyanobacteriales, belonging to the phylum Cyanobacteriota, and the order Chloroflexales, belonging to the phylum Chloroflexota. In contrast, the DSR pathway was only detected as an active sulfite-to-sulfide reduction reaction pathway in the order Burkholderiales from the phylum Pseudomonadota (Proteobacteria, Figure 2B), suggesting lower (but still present) activity compared to ASR.
For sulfate-to-sulfite reactions, both ASR and DSR were present and active in Porcelana ( Figure 1C,D and Figure 2C,D). In terms of ASR abundance and activity, sulfate-tosulfite reactions were more variable than sulfite-to-sulfide reactions ( Figures 1C and 2C). Cyanobacteriota and Chloroflexota, the most important primary producers in the Porcelana hot spring [36], remain the most abundant and active taxa for both sulfate-to-sulfite and sulfite-to-sulfide reactions. Cyanobacteriota were mostly present and active at 48 and 58 • C, while Chloroflexota were present and active at 58 and 66 • C. In addition, some phyla, such as Myxococcota, Pseudomonadota (Proteobacteria), Planctomycetota, and Bacteroidota, were present in both sulfate-to-sulfite and sulfite-to-sulfide reactions along the temperature gradient but were slightly more abundant at 48 • C than at the higher temperatures (Figures 1 and 2), suggesting their contribution in both pathways of the sulfur cycle. However, it is likely that sulfate-to-sulfite reactions are more diverse than sulfite-to-sulfide reactions, as several taxa were absent for the sulfite-to-sulfide reactions, such as the phylum Armatimonadota at lower temperatures, or Thermotogota at higher temperatures. In the case of DSR, most reads were associated with the sat gene (Table S4), which functions in both ASR and DSR, and is required for the synthesis of APS [16,67], an intermediate of the sulfate-to-sulfite reaction. Therefore, these results alone cannot truly discriminate the major taxa in DSR. However, most of the reads belonged to Cyanobacteriota, Chloroflexota, and Pseudomonadota (Proteobacteria), suggesting that these phyla might be involved in the DSR pathway. It has been previously reported that some bacteria belonging to Pseudomonadota, which is one of the main phyla found in the Porcelana hot spring, have the classical sulfur-reducing DsrAB and a sulfur-oxidizing DsrAB variant named rDsrAB [67,68], which may suggest that sulfur-redox reactions are occurring in both directions.  Table S3.  Table S3. * C denote "Candidatus".  Table S3.
For sulfate-to-sulfite reactions, both ASR and DSR were present and active in Porcelana ( Figures 1C,D and 2C,D). In terms of ASR abundance and activity, sulfate-tosulfite reactions were more variable than sulfite-to-sulfide reactions (Figures 1C and 2C). Cyanobacteriota and Chloroflexota, the most important primary producers in the Porcelana  Table S3. * C denote "Candidatus".
The Sox system proved to be the least abundant and active pathway of those studied in Porcelana (Figures 1E and 2E). The most abundant taxa showing Sox ( Figure 1E) were Chloroflexota, Pseudomonadota (Proteobacteria), and Deinococcota. However, Sox activity changed along the temperature gradient, specifically, Pseudomonadota were the most active community members at 48 and 58 • C, while Chloroflexota were the most active community members at 66 • C. Additionally, the activity of Deinococcota was highest at 66 • C ( Figure 2E). These results suggest that thiosulfate metabolic activity in the Porcelana phototrophic mats is carried out by several members of the bacterial community along the temperature gradient, being more taxonomically diverse at the phylum level at lower temperatures. This may be exemplified by the dominance of a Gammaproteobacteria soxB OTU detected at higher temperatures (see below). The lower sulfate levels far from the spring source could be explained by higher SOX activity at lower temperatures, resulting in increased thiosulfate downstream. Although we cannot exclude an effect of net ASR and DSR activity pathways on sulfate concentrations, this is not supported by the low variations in sulfide levels.
Interestingly, some "microbial dark matter" taxa in the hot spring (e.g., Binatota and Verrucomicrobiota) showed high activity at 66 • C, despite their low abundances at the DNA level (Figures 1 and 2). This suggests the potential and still unknown importance of these phyla in the sulfur cycle. In particular, the presence of 16S rRNA gene sequences from the recently described phylum Verrucomicrobiota was reported in a sulfur-active environment [69], while some MAGs belonging to Binatota (also referred to as Candidatus UBP10 [51]) were reported encoding for several genes involved in the metabolism of methylated sulfur compounds [70], suggesting, along with our results, that these phyla could have active sulfur metabolisms. However, further research is needed to confirm the relevance of these less abundant microorganisms in the sulfur cycle.

Complementary Diversity Analysis of Sulfur-Metabolizing Organisms Using aprA and soxB Marker Genes
Available metagenomic data allowed us to obtain an overview of the key microbial players in the sulfur cycle of the Porcelana hot spring. However, the lack of biological replicates (i.e., our analysis was restricted to a single metagenome per sampled temperature) limited us to a higher level of resolution and confidence for applying alpha diversity analyses. For this reason, we applied an additional complementary methodology that employed a combination of two marker genes, such as those targeting the Sox system (soxB gene) and APS reductase (aprA gene). Differential recovery and classification of taxa at different taxonomic levels have been previously reported using these marker genes [30,33].
For the clone libraries of the seven Porcelana samples along the temperature gradient (48-69 • C), the alpha rarefaction curve plateaued for the aprA gene, but not for the soxB gene ( Figure S1A,B). For soxB, two OTUs recovered in all samples were annotated as a hypothetical protein and a malate synthase, respectively, and were not considered in the subsequent analyses. An average of 100 clones per sample were analyzed for both genes using the restriction fragment length polymorphism method (with 20% clone sequencing, see Materials and Methods; Figure S1). A total of six different OTUs were obtained for the aprA gene (Figure 3), where OTU-1 and OTU-5 were ubiquitous, and OTU-1 was the most abundant across the temperature gradient. Alpha diversity (according to the Simpson and Shannon indices) for this gene (Table 1) was higher and had higher evenness indices at 58 and 66 • C compared to all other temperatures. For the soxB gene, only two ubiquitous OTUs were obtained, in which OTU-1 was dominant at all temperatures except 48 • C ( Table 1). For this gene, the Simpson index ranged from 0.3 to 0.5 at all temperatures, except at 61.2 • C, which showed the lowest diversity (0.147). Although similar diversities, as shown by the Shannon and Simpson indices, were obtained for the microbial communities of both markers (aprA and soxB) at all temperatures, a decrease in diversity was observed at temperatures close to 60 • C. However, more data are required to corroborate a correlation between diversity and temperature at the Porcelana hot spring.
Similar to Porcelana, diversity estimates for sulfur-metabolizing microorganisms in other hot springs were lower than those for the entire microbial community, as measured by 16S rRNA gene analysis (e.g., [71]). In Porcelana, the Simpson and Shannon indices for the aprA and soxB OTUs were similar to other hot springs, such as those in Tibet [9]. However, sulfur-metabolizing bacteria are more diverse in other environments, such as in freshwater ecosystems, which have higher aprAB gene richness [72], or in sediment samples from the China Sea, which show higher Shannon index values for the dsrB and soxB genes [73]. This suggests an important role for a few taxa in the generally less diverse and extreme thermal ecosystem.

Phylogenetic Affiliation and Abundance of Sulfur-Metabolizing OTUs and MAGs
The phylogenetic affiliation (Figures 3 and 4) and relative abundances of the aprA and soxB sequences ( Table 2) were estimated for each temperature and year for clone libraries. The two recovered OTUs of the soxB gene were assigned to the phylum Pseudomonadota (Proteobacteria), whereas OTUs of the aprA gene were distributed in four phyla: Pseudomonadota, GTDB Bacillota_B (former Firmicutes), Nitrospirota (Nitrospirae), and Desulfobacterota (Thermodesulfobacteria). Furthermore, proteins predicted from the MAGs and unbinned contigs of the Porcelana metagenomes at 48, 58, and 66 • C [36,44] were analyzed, providing additional data on the organisms involved in the sulfur cycle.
Conversely, OTU-2, also belonging to the order Burkholderiales (Figure 3), was the second most abundant sequence for the aprA gene, but it only dominated at 66 • C. For this OTU, the closest sequence representing a cultured organism was the isolate Thioalkalivibrio sp. HK1 (90% identity), whereas the closest sequence representing an uncultured organism corresponds to a MAG recovered from a groundwater metagenome (95% identity; GCA_004297625.1; [76]), which was suggested to be involved in the sulfur cycle [77,78] and form a new genus and family in the order Burkholderiales. Both OTU-1 and OTU-2 reveal potential SOB in new Burkholderiaceae genera that have not yet been cultivated. Furthermore, the OTU-2 clade suggests polyphyly among aprA sequences within Gammaproteobacteria, being more closely related to the GTDB phylum Nitrospirota than to the OTU-1 clade (Figure 3), which, according to other studies, could be explained by extensive horizontal gene transfer of aprAB genes [71]. . Phylogenetic reconstruction of the soxB marker gene. Maximum-likelihood tree reconstruction for the alignment of 246 sequences that include 189 CDS from the GTDB R95 protein database, 50 CDS from the Porcelana metagenomes (blue circles), and two soxB OTUs from the clone libraries (blue triangles). The tree reconstruction was performed with IQtree software using the LG + F + R7 substitution model and a non-parametric UF-bootstrap support of 1000 replicates. The leaf labels are colored according to GTDB phyla, and specifically by class for sequences from the phylum Pseudomonadota (Proteobacteria). The tree was rooted in the node between a clade of 9 outgroup sequences (trifunctional nucleotide phosphoesterase protein YfkN) and the other 246 sequences. Purple dots represent > 95% bootstrap support for that node. Maximum-likelihood tree reconstruction for the alignment of 246 sequences that include 189 CDS from the GTDB R95 protein database, 50 CDS from the Porcelana metagenomes (blue circles), and two soxB OTUs from the clone libraries (blue triangles). The tree reconstruction was performed with IQtree software using the LG + F + R7 substitution model and a non-parametric UF-bootstrap support of 1000 replicates. The leaf labels are colored according to GTDB phyla, and specifically by class for sequences from the phylum Pseudomonadota (Proteobacteria). The tree was rooted in the node between a clade of 9 outgroup sequences (trifunctional nucleotide phosphoesterase protein YfkN) and the other 246 sequences. Purple dots represent > 95% bootstrap support for that node. 1 Showing taxonomy from order to genus. Complete taxonomy can be found in Table S5. * The GTDB taxonomy of the GCF_001418245.1 genome has changed at family level between database versions being classified in f_Hydrogenophilaceae in the R95 and in f_Rhodocyclaceae in older and newer database versions.
OTU-3 also belongs to the order Burkholderiales, with 87% identity to a new genus of the family Thiobacillaceae according to GTDB R95 (closest sequence: Thiobacillus sp. UBA2186) ( Table 2), as well as with sequences from the family Hydrogenophilaceae. Members of the Thiobacillaceae and Hydrogenophilaceae families are abundant in hot spring microbial mats [20,29], and their aprA sequences [75] have been associated with aerobic and facultative chemolithotrophs isolated from hot springs [79,80].
As for the soxB gene (Figure 4), the two recovered OTUs belong to the phylum Pseudomonadota (Proteobacteria). The less abundant OTU-2 ( Table 2) was affiliated with the order Rizhobiales (Alphaproteobacteria) and showed 80% identity to Filomicrobium insigne CGMCC 1.6497. Neither the species nor the genus have been widely reported in hot springs, but detected at temperatures above 40 • C [86,87]. OTU-1 was affiliated (100% identity) with the hot spring-predominant species Hydrogenophilus thermoluteolus (GCF_003574215.1; Figure 4). Members of the family Hydrogenophilaceae (order Burkholderiales) were also found in Porcelana by aprA gene analysis. The presence of the soxB and aprA genes has been reported in some members of the phylogenetically close genera Sulfuriferula and Thiobacillus [33,88]; however, the presence of the aprA gene in the genus Hydrogenophilus has never been reported before [28,89]. In addition, two MAGs recovered at 48 • C were classified by GTDB R95 as a new genus of the family Burkholderiaceae and order Burkholderiales. These MAGs recruited one aprA and one soxB sequence each (Table S6). Their aprA sequences were associated with OTU-1, whereas the soxB sequences grouped with the Methylomirabilota clade and the genus Cupriavidus (family Burkholderiaceae) (Figure 3), respectively. This could suggest the existence of a new member of the order Burkholderiales along the temperature gradient in Porcelana that could present both functional genes.
Furthermore, compared to the clone libraries, Porcelana metagenomes recovered more diverse sequences for the soxB gene. These were associated with four different orders of Alphaproteobacteria (Table S5), revealing a high number of SOB that remain to be studied and isolated from hot springs. These soxB sequences from Pseudomonadota were associated with Gammaproteobacteria (eight sequences forming two different clades in the family Burkholderiaceae; Figure 4) and Alphaproteobacteria (13 sequences associated with the order Acetobacterales, 2 with Rodobacterales and 2 with Rizhobiales; Table S5). Another 16 sequences grouped with the family Thermaceae of the GTDB phylum Deinococcota, while 8 were closely affiliated with (but not part of the same clade as) sequences of the order Rokubacteriales (phylum Methylomirabilota in the GTDB; phylum Candidatus Rokubacteria in the NCBI database) (Figure 4), for which the potential conversion of various sulfur compounds has been suggested [90].
In summary, our results reveal that several members of the phototrophic community from Porcelana hot spring mats, especially, Pseudomonadota (Proteobacteria), Chloroflexota, and Cyanobacteriota, are more involved in the sulfur cycle than previously believed ( Figure 5). While sulfur was one of the main discriminating factors between different types of microbial communities in different hot springs in Yellowstone National Park, USA [91], and further studies determined a closed sulfur cycle in the undermat of one of the photoautotrophic microbial communities (60 • C) [17], here, we show an important role of this biogeochemical cycle in non-acidic hot spring microbial communities, as has been seen in thermophilic streamer communities and archaeal systems [91] in which the photoautotrophic metabolism is not the main input of carbon and energy.
In this case, the sulfate assimilative reduction pathway was the most represented in terms of taxonomic abundance and activity. As each of the taxa related to the sulfur cycle has different activity and occurrence along the thermal gradient at the Porcelana hot spring, the observed changes in sulfate and thiosulfate concentrations along the gradient could potentially be explained by the presence of these sulfur metabolizers and their differential activities, with the main primary producers being closely associated to the sulfur cycle, which are also reported to inhabit at different temperatures along the thermal gradient [36]; however, the role of the primary producers in the sulfur cycle of Porcelana hot spring and others non-acidic thermal systems requires further confirmation. 5). While sulfur was one of the main discriminating factors between different types of microbial communities in different hot springs in Yellowstone National Park, USA [91], and further studies determined a closed sulfur cycle in the undermat of one of the photoautotrophic microbial communities (60 °C) [17], here, we show an important role of this biogeochemical cycle in non-acidic hot spring microbial communities, as has been seen in thermophilic streamer communities and archaeal systems [91] in which the photoautotrophic metabolism is not the main input of carbon and energy.  (Table S4) for 8 different phyla and sulfur cycle steps. The size of the colored circle denotes abundance, and the color denotes the respective phylum. Temperature is shown on the Y-axis and sulfur cycle steps on the X-axis. The phylum names are written on the Y-axis, represented as a function of temperature and the sulfur cycle step where they appeared most relevant. The colored area denoting the most relevant Porcelana sulfur-associated bacteria (Pseudomonadota, Chloroflexota, and Cyanobacteriota) was hand-drawn behind the dots. In the central panel, concentrations of three sulfur species are shown based on field-measured concentrations (Table S2). In the right panel, the abundance level of photosynthesis-associated genes the of phyla Cyanobacteria, Chloroflexota, and Pseudomonadota. Temperature is shown on the Y-axis and concentration on the X-axis. * Based on data obtained from Alcamán-Arias et al. [36].  (Table S4) for 8 different phyla and sulfur cycle steps. The size of the colored circle denotes abundance, and the color denotes the respective phylum. Temperature is shown on the Y-axis and sulfur cycle steps on the X-axis. The phylum names are written on the Y-axis, represented as a function of temperature and the sulfur cycle step where they appeared most relevant. The colored area denoting the most relevant Porcelana sulfur-associated bacteria (Pseudomonadota, Chloroflexota, and Cyanobacteriota) was hand-drawn behind the dots. In the central panel, concentrations of three sulfur species are shown based on field-measured concentrations (Table S2). In the right panel, the abundance level of photosynthesis-associated genes the of phyla Cyanobacteria, Chloroflexota, and Pseudomonadota. Temperature is shown on the Y-axis and concentration on the X-axis. * Based on data obtained from Alcamán-Arias et al. [36].

Conclusions
In this work, we used a combination of metagenomic, metatranscriptomics, and molecular markers to obtain novel information on the diversity of microorganisms involved in the biogeochemical sulfur cycle in microbial mats along a temperature gradient of the neutral pH Porcelana hot spring. The presence of SOB (including highly abundant Chloroflexota) and the diversity of SRB in these microbial mats explain the complete reduction and oxidation of inorganic sulfur species in this habitat. In addition to the phyla Chloroflexota and Cyanobacteriota (previously determined to be the main carbon and nitrogen fixers in this microbial community), which were the most abundant sulfur-related organisms in the Porcelana hot spring at high and low temperatures, respectively, the phylum Pseudomonadota (Proteobacteria) showed the widest distribution across the thermal gradient. The richness of the SOB and SRB communities remained constant, although their relative abundances varied at different temperatures along the gradient. These differences allow for a decrease in sulfate levels along the gradient from the high-temperature source to the lower-temperature sites downstream, as well as a slight increase in thiosulfate levels along the same gradient. Analysis of diversity indices and phylogenetic composition of sulfur-metabolizing bacteria benefited from the complementary use of different methods and molecular markers, revealing the aprA marker in eight phyla and the soxB marker in four phyla, including a group of microorganisms with no related sequences in the databases. Recruitment of sequences in MAGs recovered from Porcelana and unbinned contigs further revealed the potential of uncultured microorganisms for both processes. In-terestingly, our study in a sulfur-poor non-acidic hot spring such as Porcelana shows similar results to a previous study indicating that metabolic redundancy is relevant in sulfur-rich environments [92], suggesting that sulfur-related metabolic redundancy is widespread in microbial-dominated environments. Taken together, these results are key to further deciphering the role of cultured and uncultured microorganisms in the sulfur cycle, and to better understanding their relationship with primary producers in the complex microbial metabolic networks of hot spring microbial mats.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/microorganisms11071803/s1, Figure S1: Rarefaction curves for aprA (a) and soxB (b) genes; Table S1: PCR primer sequences for aprA and soxB and PCR parameters; Table S2: Porcelana hot spring physicochemical parameters; Table S3: gene names and KO IDs of sulfur cycle-associated pathways; Table S4: gene counts for sulfur-associated KO IDs grouped at order level in metagenomes and metatranscriptomes obtained from Porcelana hot spring; Table S5: aprA and soxB CDS recovered from metagenomes and OTUs from Porcelana hot spring; Table S6: Genomic taxonomy of aprAor soxB-containing MAGs obtained from Porcelana hot spring according GTDB r89.