Details of the sampling campaign, DNA/RNA extraction, and sequencing were previously described (Ayala-Muñoz et al., 2020). Briefly, 5 L of water were collected in triplicate at CM at 35 m depth using a Van Dorn limnological ‘horizontal’ sampling bottle (KC Denmark A/S, Silkeborg, Denmark). For each replicate water sample, 3 L were used for collection of microbial biomass for DNA extractions, and 1 L was filtered for RNA extractions. Water was prefiltered through a 2 µm pore size glass fiber filter to reduce clogging, followed by a 0.2 µm pore size polyethersulfone (PES) sterivex filter. As in previous studies, the data pertained to planktonic cells and particles retained on a 0.2 µm filter but not retained on a 2 µm filter (Ayala-Muñoz et al., 2020; Ayala-Muñoz et al., 2022). This protocol was necessary to exclude abundant mineral particles that clogged 0.2 µm filters before sufficient microbial biomass could be obtained. We are aware that there might be microbial populations associated with mineral precipitates or debris bigger than 2 μm, or microbial populations smaller than 0.2 µm, but we expect most of the biomass to fall within 0.2–2 microns size range, a follow-up study will be necessary to characterize microbial populations outside this range.

Sterivex filters were placed immediately on ice and cryo-shipped to the Pennsylvania State University, Pennsylvania, United States, where they were stored at -80C until DNA/RNA extraction. Extractions were conducted with Qiagen DNAeasy Powerwater Kit for DNA and Qiagen RNAeasy PowerMicrobiome Kit (Qiagen, Venlo, Netherlands) for RNA. Quality and quantity of the extracts was verified using Qubit^® 2.0 Fluorometer (Invitrogen, Carlsbad, CA, United States) and Bioanalyzer 2100 RNA 6000 Pico Assay (Agilent, Santa Clara, CA, United States). Only the best two replicates were used for nucleic acid sequencing. Metagenome library preparation was performed using Illumina’s NexteraXT library preparation kit (Illumina, San Diego, CA, United States). Double stranded cDNA synthesis and metatranscriptome library preparation was performed using the Tecan RNA Trio library preparation kit (Tecan, Mannedorf, Switzerland). Sequencing was conducted on an Illumina Hiseq 4000 platform (Illumina, San Diego, CA, United States) using 150 bp paired end chemistry.

Details of the processing of metagenomes and metatranscriptomes were also previously described (Ayala-Muñoz et al., 2020). Two metagenomic and two metatranscriptomic datasets were obtained. After quality-filtering and trimming, metagenomic short-reads were de-novo assembled with Megahit v1.1.2 (Li et al., 2016). We assembled the two replicate metagenomes separately and also produced a co-assembly. The metatranscriptomes were subjected to in silico removal of ribosomal RNA sequences using sortmeRNA v2.1 (Kopylova et al., 2012). The remaining RNA reads (mostly mRNAs) were mapped to the metagenomic co-assembly with BBMap (JGI, 2017) (26) (JGI, 2017) (26) (JGI, 2017) (26) (JGI, 2017) (26) (JGI, 2017) (26) (Kopylova et al., 2012) (25) (Kopylova et al., 2012) (24) (Ayala-Muñoz et al., 2022) (min_id = 0.95 and slow mode). Approximately 50% of the mRNA reads from only one of the two metatranscriptomes (CM35_1) mapped to the respective metagenome contigs. The other metatranscriptomic dataset (CM35_2) with less than 1% of mapped mRNA reads to its respective metagenomic co-assembly was not used for downstream analysis. We speculate that the metatranscriptome from CM35-2 failed to map to its corresponding metagenome because of a contamination issue.

We applied a combination of bioinformatic approaches to reconstruct metagenome-assembled genomes (MAGs) from metagenomes. Two de-novo individual assemblies (one per replicate) and one co-assembly were obtained using Megahit v1.1.2 with default parameters (Li et al., 2016). A subset of randomly selected reads representing 100% coverage of the most abundant taxa in each metagenome was also assembled with metaSPAdes employing the default parameters (Nurk et al., 2017). The latter approach was conducted to maximize the binning of MAGs from the most abundant taxa (i.e., Thermoplasmatales). We mapped reads to each assembly with Bowtie v2.2.4 using default parameters (Langmead et al., 2009). Maxbin 2.0 (-min_contig_length = 1500, -prob_threshold 0.95, and -markerset 40 for archaeal dominated samples) (Wu et al., 2016) and MetaBat v2.12.1 (Kang et al., 2019) were used to bin the two individual assemblies, the co-assembly, and the read-subset assemblies. A total of 361 MAGs were obtained using this branched pipeline. DASTool (Sieber et al., 2018) and dRep v2.3.2 (Olm et al., 2017) with default options were used to select the best non-redundant MAGs (a total of 43). CheckM was used to analyze the quality of the MAGs. Finally, we used FastAni v1.3 (Jain et al., 2018) to select MAGs with an average nucleotide identity (ANI) of <96.5% as representatives of potentially different species. MAGs with ANI >96.5% were grouped together and only the best quality MAG from each group was further analyzed (Supplementary Table S1), resulting in a total of 18 MAGs (Supplementary Table S2).

Taxonomic annotation of the 18 MAGs was conducted with GTDB-Tk v1.3 (Genome Taxonomy Database) (Chaumeil et al., 2019) and the microbial genome atlas webserver (MiGA, NCBI-Prok Database) (Rodriguez-R et al., 2018). In parallel, genes from each MAG were predicted with Prodigal v2.6.3 (Hyatt et al., 2010). Genes and contigs from each MAG were taxonomically annotated with the bins annotation tool (BAT), using the CAT_prepare_20200618 database files and flag --top 11 (von Meijenfeldt et al., 2019). Contigs from a MAG corresponding to a different phylum from that assigned by GTDB-Tk were manually removed prior to downstream analyses. CheckM and GTDB-Tk were run again over the curated MAGs. A multiple sequence alignment obtained from GTDB-Tk and based on 122 archaeal- and 120 bacterial-specific marker proteins was used to construct a concatenated phylogenetic tree. Maximum likelihood (ML) phylogenies for archaeal and bacterial genomes were inferred with raxML (-x 12345 -p 2352890 -# 1000 -m PROTGAMMAAUTO) (Stamatakis, 2014). The tree was customized with Itol (Letunic and Bork, 2016). DNA reads were mapped to the contigs of each MAG using BBmap.sh (min_id = 0.95, slow mode), and relative abundances of taxa were calculated based on total DNA reads mapped to the metagenomic coassemblies.

The curated versions of the 18 MAGs were functionally annotated with METABOLIC v4.0 (Zhou et al., 2020). This bioinformatic pipeline predicts genes with Prodigal V2.6.3. Annotation was conducted with hmmsearch (HMMER 3.1b2) (Eddy, 2011) using databases including KOFAM (Aramaki et al., 2019), TIGRFAM (Haft et al., 2003), Pfam (Finn et al., 2014), CAZy (Lombard et al., 2013) and a database of 143 custom HMM profiles for genes involved in sulfur metabolism. The quality-filtered DNA reads (metagenome) and mRNA reads (metatranscriptome) from CM35_1 (replicate 1), were mapped to the predicted genes of each MAG with BBmap (min_id = 0.95, slow mode) (JGI, 2017; Miller et al., 2011) (JGI, 2017; Miller et al., 2011). TPM values summed per functional group were used to yield normalized values for gene or transcript abundance considering the relative abundance of each MAG. Normalized gene and transcript abundances were calculated as TPM values (number of DNA or mRNA reads mapped to the predicted gene (*10⁶) divided by length/sum of DNA or mRNA reads mapped to all genes divided by length per MAG). The gene and transcript profiles were converted to functional profiles by summing the normalized abundance of the annotated genes in each functional group. Plots were generated in R and customized in Inkscape. The KEGG annotations of each MAG coupled with KEGG Mapper (Kanehisa and Sato, 2020) were used for generating conceptual representations of MAGs with predicted in-situ activity for biosulfidogenesis. The representations were constructed in BioRender.com.

The 18 MAGs retrieved in this study represented the majority of phyla previously described in the deep layer of Cueva de la Mora. Total relative abundances were based on 16S rRNA amplicon sequencing and 16S rRNA reconstruction from metagenomes using EMIRGE (Miller et al., 2011). The most abundant archaeal phylum (Euryarchaeota; 46.0% by EMIRGE, 46.3% by amplicon sequencing) was represented by four MAGs (EUR_01, 06, 14 and 15) with a total relative abundance of 23.8% (Figure 1 lower panel). The most abundant bacterial phylum (Parcubacteria; 18.7% by EMIRGE, 17.6% by amplicon sequencing) was represented by only one MAG (PAT_13) with a relative abundance of 0.2%. The other phyla were represented by one or two MAGs with total relative abundances comparable to what was previously described based on 16S rRNA amplicon sequencing (Ayala-Muñoz et al., 2022) and 16S rRNA reconstruction from metagenomes using EMIRGE (Ayala-Muñoz et al., 2020). DOR_08 represented the Candidatus Dormibacteria phylum. This phylum, formerly known as candidate division AD3, has been reported in soil or sediment samples of extreme environments including acid mine drainage sites (Mesa et al., 2017). In the Silva taxonomic database v.132 (Quast et al., 2013), AD3 was a class of the phylum Chloroflexi, which is how it was presented previously (Ayala-Muñoz et al., under revision). No MAGs from the Firmicutes or Thaumarchaeaota phyla were reconstructed.

None of the eighteen MAGs had high ANI (>95%) to genomes included in the Genome Database Taxonomy or the NCBI prokaryotic database (Supplementary Table S2). Eleven MAGs belonged to the Bacteria domain (Figure 1). ACI_09 (Acidobacteria) formed a separate branch from the clade that includes Granulicella, Terracidiphilus and Acidipila. NIT_04 (Nitrospirae) is closely affiliated with the uncultured Thermodesulfovibrionales JdFR-88 and separate from the clade formed by cultured Thermodesulfovibrio strains. PRO_16 (Proteobacteria) formed a clade with Candidatus Acidulodesulfobacter (formerly Sva0485) strains. PRO_05 clustered with Desulfomonile tiedjie in a clade shared with other Desulfobacteraceae. ACT_02 and ACT_11 (Actinobacteria) fell into a clade with Thermoleophilia strains. PAT_13 (Parcubacteria) was loosely affiliated with Colwellbacteraceae. Finally, DOR_08, CHL_20, CHL_03 are part of the phyla Candidatus Dormibacteria (formerly candidate division AD3 represented as a class division of Chloroflexi in the Silva taxonomic database v.132) and Chloroflexi. DOR_08 formed a clade with Dormibacteria strains, while CHL_20 and CHL_03 fell into a clade with Dehalococoidia strains. CHL_20 and CHL_03 potentially belong to different families. CHL_20 clustered with Dehaloccoidia CF_162 while CHL_02 clustered with Dehalococcoidia SZUA-161.

Seven MAGs belonged to the Archaea domain (Figure 1). EUR_01, 06, 14, and 15 (Euryarchaeota) clustered within the order Thermoplasmatales. EUR_01, 15 and 14 are affiliated with the uncultivated Thermoplasmatales sp01856825. EUR_06 is part of a different clade where cultured Thermoplasmatales were found. NAN_12 (Nanoarchaeota) is affiliated with the Woesearchaeles. MIC_10 (Micharchaeota) is affiliated with the Micharchaeaceae order. CRE_19 and CRE_07 (formerly Crenarchaeota, now Thermoprotei) are part of the Nitrosphaerales clade. CRE_19 is affiliated with Nitrosospaherales UBA141 and CRE_07 is affiliated with Nitrososphaerales UBA183.

While MAGs are in silico representations of environmental species populations, for brevity, we will refer to MAGs containing and expressing genes. The majority of MAGs contained and expressed genes involved in organic carbon oxidation (Figure 2, Supplementary Table S3). We have chosen to use dot plots in Figures 2–4, and Figure 7 primarily to present both the presence/absence of functional potential and gene expression, and to distinguish gene complexes with missing genes (filled vs. unfilled dots). We do not evaluate or discuss our results based on quantitative TPM values because mapping genomic and transcriptomic reads back to MAGs could produce artifacts. As far as we know, there is no standard way yet to perform this type of analysis. Based on recent studies, we decided to implement an approach like the one described by Salazar et al. (2019) (Salazar et al., 2019).

Bacterial MAGs, with the exception of PAT_13, have complete (100%) or close to complete pathways (>50%) for the tricarboxylic acid (TCA) cycle. Among the archaeal MAGs, only CRE_19 and EUR_16 had a TCA cycle with >50% completeness. With respect to the pentose phosphate pathway (PPP), most of the bacterial MAGs had almost complete PPP, oxidative (ox) and non-oxidative (non-ox) branches. Archaeal MAGs only had the PPP non-oxidative branch with CRE_19 and all EUR_MAGs having the most complete pathways. Most MAGs contained >50% complete pathways for gluconeogenesis and glycolysis. The bacterial MAGs ACI_09, CHL_03, CHL_20, DOR_08, PRO_05, PRO_16 had >50% complete pathways for beta oxidation of fatty acids.

With respect to one-carbon metabolisms, most of the MAGs lacked genes involved in acetoclastic methanogenesis or had incomplete pathways (<50%, Figure 2, Supplementary Table S3). We also searched for genes involved in carbon monoxide (CO) oxidation (coxSML), methanol oxidation (mdh), and formate oxidation (fdo, fdh). Strikingly, all four EUR_ MAGs and CRE_19 had genes involved in CO oxidation, some of which were actively expressed. Two bacterial MAGs (ACI_09 and PRO_05) contained at least one gene involved in CO oxidation whereas CRE_19 had all the Cox subunits (coxSML). ACI_09, EUR_01, and EUR_06 only had coxS, and EUR_14 and EUR_15 only had coxSM. Carbon monoxide can be a product of peroxidation of membrane lipids due to oxidative stress from metal toxicity (Wolff and Bidlack, 1976). Almost all archaeal MAGs (CRE_19 and EUR_01, 06, 14, 15), contained, and in some cases expressed, genes involved in methanol oxidation, whereas few bacterial MAGs (CHL_03, CHL_20, and PRO_05) contained those genes. Methanol can be a subproduct of dead phytoplankton sinking from the upper layers (Mincer and Aicher, 2016). Finally, with respect to formate oxidation, eight MAGs contained genes encoding for one of the subunits in either FdoGH or FdhAB while three MAGs (ACT_02, NIT_04, and EUR_14) expressed these genes.

All of the bacterial MAGs and three archaeal MAGs contained and expressed genes involved in carbon fixation (Figure 2, Supplementary Table S3). The ACT_MAGs, CHL_03, NIT_04, and PRO_05 had nearly complete Wood-Ljungdahl pathways (>70%) and all of these MAGs (except for ACT_11) expressed at least one gene involved in the Wood-Ljungdahl pathway. Some MAGs had >50% complete reductive citrate cycle (rTCA) pathways but none of them contained the rTCA autotrophy marker gene coding for ATP-citrate lyase (aclAB). Three MAGs (CHL_03, PRO_05, and PRO_16) had a >50% complete dicarboxylate hydroxybutyrate cycle (2HC-4HB) and also had the marker gene abdfD. Finally, despite that most MAGs had a >50% complete pathway for the Calvin-Benson-Bassham cycle (CBB), only CRE_07, MIC_10, NAN_12 had the rubisco genes. CRE_07 and MIC_10 expressed the rubisco genes. It should be noted that Rubisco genes in Archaea typically encode the form III and II/III Rubisco that perform light independent incorporation of CO₂ into sugars derived from ribonucleotides (Wrighton et al., 2016).

Both bacterial and archaeal MAGs contained and expressed genes involved in fermentation (Figure 2, Supplementary Table S3). Genes involved in transformation of Acetyl-CoA to acetate (reversible) are found in almost all MAGs with the exception of PAT_13 and EUR_06, and seven MAGs expressed gene(s) related to this reaction. Twelve MAGs contained genes to ferment acetate to pyruvate (porA), but only ACT_02, NIT_04, EUR_01 and EUR_14 expressed them.

The majority of MAGs contained genes for degrading oligo- and polysaccharides but not all expressed these genes (Supplementary Figure S1). Twelve MAGs expressed at least one gene in a family of glycosyl hydrolases involved in degradation of oligo- and polysaccharides. MAGs ACT_02, CHL_03, NIT_04, PRO_05, EUR_01, and EUR_14 expressed multiple (≥2) gene families. From these MAGs, all except PRO_05 had a relative abundance of ≥3%. The most abundant MAG (EUR_01) had transcripts in gene families GH 15, 57 and 133. GH 15 includes an anomer-inverting glucosidase, trehalase (EC 3.2.1.28 found in Thermoplasma acidophilium and T. volcanium), and glucoamylase (EC 3.2.1.3 found in T. acidophilium). GH57 includes glycosyl hydrolases that act on starch, glycogen, and related oligo- and polysaccharides. Finally, GH133 includes an amylo-alpha-1,6-glucosidase (EC 3.2.1.33).

The bacterial MAGs mainly expressed genes involved in dissimilatory sulfate reduction (Figure 3, Supplementary Table S4). The first step in this pathway is the reduction of sulfate to sulfite and catalyzed by sat and aprAB. The sat gene encodes for sulfate adenyltransferase that catalyzes the formation of adenosine 5′-phosphosulfate (APS) from inorganic sulfate and ATP. The aprAB gene encodes for an adenosine 5′-phosphosulfate reductase that catalyzes the reduction of APS to sulfite while interacting with the membrane QmoABC complex in the electron transport chain (Anantharaman et al., 2018). Among the eighteen MAGs, thirteen contained sat but only four expressed this gene (bacterial MAGs ACT_02, NIT_04, PRO_05, and the archaeal MAG EUR_14). Ten MAGs (including the archaeal CRE_17, CRE_19, and EUR_06) contained genes for AprAB. CRE_19 and all the bacterial MAGs with the exception of CHL_20 expressed aprAB. Gene parts of the QmoABC complex were found in almost all bacterial MAGs but in only one archaeal MAG (CRE_19). The genes qmoABC were expressed in ACT_02, ACT_11, NIT_04, and PRO_16.

The second step in dissimilatory sulfate reduction is the reduction of sulfite to sulfide (Figure 3, Supplementary Table S4). Genes involved in this step are dsrAB, dsrC, dsrD, dsrTMK, and dsrJOP. The dsrAB genes encode for a dissimilatory sulfite reductase which forms S(II) and S (0) intermediates. DsrAB directly interacts with DsrC for the conversion of the S (0) intermediate to sulfide. The dsrMKJOP genes encode for a transmembrane complex involved in electron transport to restore DsrC to its reduced state (Venceslau et al., 2014). DsrMK are the minimum subunits required for electron transfer for sulfite reduction (Pereira et al., 2011). The dsrD and dsrT genes are likely involved in a regulatory function (Venceslau et al., 2014; Anantharaman et al., 2018). All bacterial MAGs with the exception of CHL_20 contained the genes dsrAB but only six MAGs expressed these genes. Among the archaeal MAGs, CRE_19 and EUR_06 contained dsrAB but only CRE_19 expressed these genes. Seven MAGs contained dsrC but only ACT_02, CHL_03, NIT_04 and CRE_19 expressed this gene. Only bacterial MAGs contained dsrTMK and/or dsrJOP genes. Finally, dsrD was found in all bacterial MAGs with the exception of CHL_20 and PRO_05. The MAGs ACT_02, ACT_11, CHL_03 and NIT_04 expressed dsrTMK and dsrD.

Other genes involved in reductive S processes are ttrABC, phsA, and hydBG (Figure 3, Supplementary Table S4). TtrABC are involved in tetrathionate reduction to thiosulfate. PhsA is involved in disproportionation but also in thiosulfate reduction to sulfide. HydBG are involved in sulfur (0) reduction to sulfide (Vigneron et al., 2021). Five MAGs contained and expressed ttrABC including the archaeal MAG CRE_19. Three MAGs contained phsA and ACT_11 and PRO_05 expressed this gene. Four of the bacterial MAGs contained the hydBG genes but only ACT_02 expressed them.

Both bacterial and archaeal MAGs contained genes and/or transcripts involved in sulfide and/or S (0) oxidation (Figure 3, Supplementary Table S4). In this pathway, Sqr oxidizes sulfide to polysulfide, which can be further oxidized to sulfite by the complex DsrABCEFH (Gregersen et al., 2011) or by a periplasmic sulfur dioxygenase (Sdo) (Rohwerder and Sand, 2003). DsrEFH is central in this pathway since it donates sulfur to DsrC (Stockdreher et al., 2012). DsrC then transfers the sulfur to DsrAB acting in the oxidative direction to produce sulfite (Venceslau et al., 2014). On the other hand, sulfur oxidation can also be catalyzed by Sdo, which could also act on perthiols formed by the spontaneous interaction between thiols and elemental sulfur and sulfide (Rohwerder and Sand, 2003). The generated sulfite could be further converted into sulfate by the Sat/AprAB complex (Anantharaman et al., 2018). Nine MAGs had sqr and six expressed this gene. Among the ones expressing sqr are three of the Euryarchaeota MAGs, including the second most abundant EUR_14. Four MAGs had sdo but only the Chloroflexi MAG (CHL_03) expressed this gene. Finally, the Actinobacteria MAGs contained and expressed genes for the complete complex DsrEFH.

A phylogenetic analysis of the dsrAB genes found in the MAGs showed that the majority of them clustered with either the bacterial or archaeal reductive type (Supplementary Figure S2). However, the Actinobacteria MAGs each had two copies of the dsrAB genes. Copy 1 clustered with the reductive bacterial type while copy 2 clustered with the oxidative bacterial type. From this genomic context analysis, the dsrAB copy 1 from ACT_MAGs were closer to genes coding for DsrMKCBD while copy 2 is closer to the genes coding for DsrEFH.

Few MAGs contained genes involved in the oxidation of thiosulfate (Figure 3, Supplementary Table S4). Three bacterial MAGs contained genes for thiosulfate oxidation to tetrathionate (tsdA and doxD) or sulfate (soxXAYZB). The Nitrospirae MAG (NIT_04) is the only MAG that expressed soxXAYZB. Rhodanase-like sulfur-transferase (Tst) is another enzyme known for its involvement in thiosulfate oxidation (Vigneron et al., 2021). However, it was also found under sulfur-reducing conditions (Florentino et al., 2019). We found tst in three MAGs, two of which (CHL_20 and NIT_04) expressed this gene.

Most of the bacterial MAGs expressed at least one gene involved in dimethyl sulfoxide (DMSO) metabolism (Figure 3, Supplementary Table S4). Ddh specifically catalyzes the transformation of DMSO to dimethyl sulfide (DMS) (Vigneron et al., 2021). Nine MAGs contained genes from the dimethyl sulfoxide dehydrogenase complex (ddhABC) and were expressed in five bacterial MAGs (ACT_02, CHL_03, CHL_20, NIT_04 and PRO_05). Seven MAGs contained dmoB, which encodes for a dimethyl-sulfide monooxygenase subunit B, but only NIT_04 expressed this gene. Only one MAG (CHL_03) contained dmsAB, encoding for an anaerobic dimethyl sulfoxide reductase, but these genes were not expressed.

Bacterial MAGs contained genes involved in oxidative phosphorylation, oxygen reduction, hydrogen metabolism, nitrate/nitrite reduction, and Fe(III) reduction (Figure 4, Supplementary Table S5). All MAGs except for PAT_13 contained genes involved in the complex I (nuoABC) and complex V (atpAD F-type) for oxidative phosphorylation, and ACT_02, CHL_03, and NIT_04 expressed these genes. Most bacterial MAGs contained a gene involved in oxygen reduction (cydAB, cyoABCD, ccoNOP, and coxAB). ACT_02 expressed ccoNOP and NIT_04 expressed cydAB. Five of the bacterial MAGs contained genes involved in hydrogen metabolism (nife-group1). NIT_04 and PRO_05 expressed nife-group-1 involved in hydrogen oxidation. Three bacterial MAGs contained genes involved in nitrate reduction (narGH) and NIT_04 expressed narGH. The gene nrfA involved in nitrite reduction to ammonia was expressed by PRO_16. Finally, the gene mtrC involved in Fe(III) reduction was expressed by NIT_04.

In contrast to the bacterial MAGs, archaeal MAGs only contained genes involved in oxidative phosphorylation (Figure 4, Supplementary Table S5). The Euryarchaeota-MAGs, EUR_01, EUR_06, and EUR_14 contained nuoABC genes. All archaeal MAGs except for NAN_12 contained atpAB V/A-type genes, and these were expressed by CRE_19, EUR_01, and EUR_14. These three MAGs and EUR_15 also expressed a gene from the CydAB complex, a high affinity oxygen reductase associated with oxygen-limited conditions (Borisov et al., 2011).

Microbial strains unrelated to cultured representatives of known sulfate reducing bacteria actively contributed to biosulfidogenesis in the deep layer of CM and were metabolically versatile. Six novel strains represented by bacterial MAGs ACI_09, ACT_02, ACT_11, CHL_03, NIT_04, and PRO_16 contained genes for the reduction of sulfate to sulfide, including the dsrD gene (Figure 3). Populations represented by these six MAGs accounted for roughly 15% of the microbial community (Figure 1). The dsrD gene is an important marker for sulfate reduction given its absence in the genomes of sulfide oxidizing organisms (Rabus et al., 2015; Anantharaman et al., 2018). Reduction of sulfate could be coupled to oxidation of organic carbon or to hydrogen and indeed, these six MAGs expressed genes for acetate (Figure 2) or H₂ oxidation (Figure 4). All of these MAGs except for ACI_09 contained a nearly complete Wood-Ljungdhal pathway, suggesting that they could fix inorganic carbon when organic carbon (e.g., organic acids) is not available (Figure 2). This hypothesis is supported by the high concentration of CO₂ (29 mM) available for the microbial populations in the deep layer (Sánchez-España et al., 2020b). Some of these MAGs also contained genes involved in sulfide, S (0), or thiosulfate oxidation (Figure 3), which could be coupled with oxygen, nitrate, or ferric iron reduction (Figure 4). Our results suggest that metabolic versatility is beneficial to these populations, and along with evidence for sulfide oxidation in the abundant Thermoplasmatales populations, points toward transient changes in redox potential and dissolved oxygen concentrations in the CM deep layer.

Remarkably, the most abundant bacterial MAGs (ACT_02, CHL_03, and NIT_04) expressed genes for almost the entire set of enzymes related to biosulfidogenesis (Figure 3), and represented novel taxa. The taxonomic analysis of these MAGs with GTDB-Tk did not report ANI ≥96.5% to any genome found in the Genome Taxonomy Database. ACT_02 (Actinobacteria, class Thermoleophilia (Figure 1, Supplementary Table S1)) populations made up roughly 3% of the total community. The closest genome related to this MAG was Adlercreutzia equolifaciens DSM 19450 with an AAI (average aminoacid identity) of 43% according to the NCBI_Prok database in the Microbial Genome Atlas (Supplementary Table S2). The few acidophilic Actinobacteria described to date include ferrous-iron oxidizers (e.g., Acidimicrobium ferrooxidans), ferric-iron reducers (Aciditerrimonas ferrireducens), and/or sulfur-oxidizers (Acidithiomicrobium spp.) (Clark and Norris, 1996; Davis-Belmar and Norris, 2009; Itoh et al., 2011), and none with the ability for performing biosulfidogenesis, only, the neutrophile Gordonibacter pamelaeae is among the few cultivated Actinobacteria with dsrAB genes (Müller et al., 2015) along with other neutrophilic, non-cultured Actinobacteria (Anantharaman et al., 2018). Within the Thermoleophilia class, the few known representatives, retrieved from either hot springs or soil samples, have no reported contribution to sulfur cycling (Hu et al., 2019). However, the metagenomic evidence found in this work suggests that ACT_02 populations are organoheterotrophs able to reduce sulfate and oxidize sulfide (Figure 5). They have the potential to oxidize pyruvate to acetyl-CoA by a pyruvate ferredoxin oxidoreductase (PorABCD) and a 2-oxoglutarate/2-oxoacid ferredoxin oxidoreductase (KorAB). ACT_02 populations also have an acs involved in Acetyl_CoA generation from acetate, complete glycolysis (Embden-Meyerhof) and non-oxidative pentose phosphate pathways, and all the genes encoding for the methyl and carbonyl branches of the Wood-Ljungdhal pathway except for fdhAB. The lack of the genes encoding for Pta and AcdA suggests that ACT_02 could not grow autotrophically (Schuchmann and Müller, 2014). Energy generation is supported by the NADH quinone dehydrogenase complex (Nuo) and dissimilatory sulfate reduction (Sat, AprAB, Qmo, DsrABCDMK). ACT_02 populations were also capable of sulfide oxidation (Sqr and DsrEFH). They have a second (oxidative) copy of dsrAB located closer to the genes dsrEFH (Supplementary Figure S2), and the genes necessary for O₂ respiration, including CydBD.

CHL_03 (Chloroflexi, class Dehalococcoidia (Figure 1, Supplementary Table S1)) populations made up roughly 3% of the total community. The closest genome related to this MAG was Dehalococcoides mccartyi VS. NC 013552 with an AAI of 44% according to the NCBI_Prok database in the Microbial Genome Atlas (Supplementary Table S2). Previous studies presented evidence of acidophilic Chloroflexi populations in environments including acidic pit lakes (Santofimia et al., 2013; García-Moyano et al., 2015; Mesa et al., 2017), but their metabolisms has not determined. A single-cell genome from this class was retrieved from a marine subsurface sediment and is the only reported genome in the class to date with dsrAB genes (Wasmund et al., 2016). Analysis of the CHL_03 MAG points to a fermentative lifestyle with the additional ability to reduce sulfate (Figure 5). CHL_03 lacked four out of ten subunits of the NADH-ubiquinone oxidoreductase (Nuo) complex. The interaction of NuoEF with its membrane anchor NuoLMN subunits is improbable given the lack of the connecting subunits NuoCDGI. CHL_03 also lacked standard electron transport chain components including cytochrome c reductases (Complex III) and cytochrome c oxidases (Complex IV). It contained acs and porA involved in acetate and pyruvate oxidation to Acetyl-CoA. CHL_03 had also an incomplete TCA cycle with a lack of citrate synthase, citrate lyase, and aconitate hydratase. CHL_03 contained a complete set of genes involved in glycolysis and gluconeogenesis. Such evidence suggests a primarily fermentative metabolism relying on sugars for energy generation. Acetyl-CoA could also be generated by fatty acid degradation and by CO₂ reduction via the Wood-Ljungdahl pathway. Autotrophic growth by acetogenesis is potentially feasible thanks to ATP generated by AcdA. CHL_03 populations could also carry out sulfate reduction (Sat, AprAB, Qmo, DsrABCDTMKJOP) coupled to hydrogen oxidation (nife-1) to generate a proton motive force which fuels ATP formation by an F-type ATPase. CHL_03 MAG was 97.7% complete and, although unlikely, it is possible that the lacking genes were missed due to the incompleteness of the genome.

NIT_04 (Nitrospirae, Figure 1, Supplementary Table S1) populations made up roughly 3% of the total community. The closest genome related to this MAG was Thermodesulfovibrio yellowstonii DSM_11347_NC_011296 with an AAI of 45.97% according to the NCBI_Prok database in the Microbial Genome Atlas (Supplementary Table S2). The most renowned acidophilic isolates from this phylum are Leptospirillum strains known for their low pH iron-oxidation capability (Méndez-García et al., 2015). However, NIT_04 populations were affiliated with the genus Thermodesulfovibrio. To date all Thermodesulfovibrio isolates are neutrophiles or alkaliphiles (Henry et al., 1994; Frank et al., 2016), there are no reported genomes of acidophilic Thermodesulfovibrio strains. However there is molecular evidence of acidophilic Thermodesulfovibrio presence in at least one acidic environment from anoxic hot spring sediments (Willis et al., 2019). Some members of this genus such as Thermodesulfovibrio thiophilus (Sekiguchi et al., 2008) have the reductive bacterial-type DsrAB (Müller et al., 2015). NIT_04 populations are organoheterotrophs that can reduce sulfate and oxidize sulfide and thiosulfate (Figure 5). Acetate and pyruvate could be oxidized to acetyl-CoA thanks to the presence of Acs and PorA. Fatty acid degradation could also be used to produce Acetyl-CoA. NIT_04 also contained complete sets of genes for glycolysis, gluconeogenesis, and the pentose phosphate pathway (oxidative and non-oxidative branches). NIT_04 had an incomplete methyl branch and a complete carbonyl branch of the Wood-Ljungdahl pathway. However, NIT_04 also had genes encoding for MtaA, a methyl transferase involved in the methyl branch of the WL (not shown in Figure 5), and AcdA involved in energy generation, suggesting that NIT_04 could utilize methyl compounds to get energy using the Wood-Ljungdahl pathway. NIT_04 reduced sulfate all the way to sulfide (Sat, AprAB, Qmo, DsrABCDTMKJOP) coupled to hydrogen oxidation (nife-group-1). NIT_04 also expressed sqr, involved in sulfide oxidation, and the Sox complex (soxXAYZB), involved in thiosulfate oxidation. Sulfide or thiosulfate oxidation could be coupled to oxygen respiration (CydBD). Coupling of electron transfer to energy conservation could be mediated in NIT_04 by an H/Na-pumping Rnf complex, an NADH-quinone oxidoreductase (Nuo) and an F-type ATPase synthase. Genes involved in nitrate reduction (narGH) and metal reduction (mtrC) were also found, providing evidence of the potential to use a diverse range of electron acceptors in respiration.

Cueva de la Mora (CM) is one the world’s most studied AMD sites and represents a model system for discovering novel acidophilic and acid-tolerant organisms. While restoration of water quality to pre-mining conditions may be financially prohibitive, stimulation of biosulfidogenesis in the deep layer of permanently stratified pit lakes such as CM may be a cost-effective means to limit environmental risks posed by elevated metal (loid) concentrations. The most abundant populations in the deep layer are Thermoplamatales-group archaea that likely scavenge organic carbon and potentially oxidize sulfide using trace concentrations of oxygen. Therefore, we focus our discussion on MAGs representing abundant populations that contribute to biosulfidogenesis under in situ conditions. Based on our findings, we speculate on how biosulfidogenesis could be enhanced in meromictic acidic pit lakes.

Our results strongly suggest that organic carbon (C_org) additions will be required to stimulate biosulfidogenesis in the CM deep layer, both to fuel sulfide production directly and to consume oxidants (e.g., O₂) which, although present at low concentrations, currently allow for the rapid re-oxidation of reduced sulfur compounds. A stumbling block in historical biosulfidogenesis efforts in AMD systems is the high toxicity of organic acids at low pH (Pinhal et al., 2019). Depending on metabolic networks in the C_org-amended community, this toxicity may become limiting for sulfide production by organic acid-oxidizing sulfate-reducing microbial populations. In this case, the addition of S (0) in combination with C_org at lower loadings (to fuel S (0) reduction) or no C_org addition (to fuel S (0) disproportionation) are alternative strategies that may result in higher sulfide production due to lower organic acid toxicity. These alternative strategies will be tested in future laboratory work with enrichments from the CM deep layer.