About 1 g of solid sample (fumarolic deposits and mud) was added to 20 ml of MilliQ water and extracted as described in Lezcano et al. (2019). The concentration of water-soluble anions (F⁻, Cl⁻, NO₂⁻, Br⁻, NO₃⁻, PO₄³⁻, and SO₄²⁻) in extracts and water samples was measured in triplicate using an 881 compact IC pro ion chromatography system (Metrohm, Switzerland) with a Metrosep A Supp 5-250/4.0 column. The concentration of NH₄⁺ in extracts and water samples was measured using the alkaline hypochlorite/phenol nitroprusside method, after adding sodium citrate to prevent the precipitation of calcium and magnesium salts (Solorzano, 1969). Before measuring concentrations of NH₄⁺ in water samples, samples were diluted 10-fold. Ammonium sulfate (NH₄)₂SO₄ solutions (0, 660, 1,320, 3,300, 6,600, and 13,200 ppm) were prepared and used as standard. The absorbance of each sample was measured in triplicate using an Evolution 60S UV-Vis Spectrophotometer at 640 nm wavelength.

X-ray fluorescence (XRF) was used to examine the chemical composition of mud and fumarolic deposits samples. Samples were analyzed in triplicate using a Horiba XGT-1000WR X-ray Fluorescence with an Rh tube X-ray source and elemental wt.% was determined using the XGT-1000WR software’s quantification (Supplementary Method).

Carbon, nitrogen, and sulfur stable isotope and elemental concentration analyses were performed at the Earth and Planets Laboratory, Carnegie Institution for Science. For organic carbon (TOC), and organic δ¹³C measurements, samples were weighed into silver boats and fumed with 12 N HCl for 12–14 h. δ¹³C and δ¹⁵N isotopes were measured using a Thermo Scientific Delta VPlus isotope ratio gas-source mass spectrometer connected to a Carlo Erba (NA 2500) elemental analyzer (EA/IRMS) via a Conflo III interface. δ³⁴S was analyzed by the same gas-source mass spectrometer but connected to an Elementar Americas vario Micro CUBE elemental analyzer (EA-IRMS) via a Conflo III interface. Stable isotope values are reported in standard delta notation as ‰ variations relative to: Pee Dee Belemnite (PDB) for δ¹³C, atmospheric N₂ gas (AIR) for δ¹⁵N, Vienna Canyon Diablo Troctolite (V-CDT) for δ³⁴S, and with an analytical error of ±0.1‰. In-house δ¹³C and δ¹⁵N standards were also used and calibrated against international and certified standards as well. Additional δ³⁴S standards include the International Atomic Energy Agency reference materials IAEA-S-1 (δ³⁴S = −0.3‰), and IAEA S3 (δ³⁴S = −32.3‰), as well as NBS-123 (δ³⁴S = −17.09‰), NBS-127 (δ³⁴S = −21.17‰), and USGS-42 (δ³⁴S = −7.84‰). A subset of the samples was analyzed in 2021 at the EDGE Stable Isotope Laboratory at the University of California Riverside to confirm the very negative nitrogen isotope values in some of these samples. USGS25 (δ¹⁵N = −30.41‰) and USGS40 (δ¹⁵N = −4.52‰) were used as two of our calibrating standards and confirmed the values measured earlier at the Earth and Planets Laboratory.

About 0.25 g of mud and fumarolic deposits was used for DNA extraction. Before DNA extraction, fumarolic deposits were powdered using a sterilized agate mortar and pestle. Water samples were filtered in the laboratory using a 0.2 μm VWR black polycarbonate filter. Filters were cut into small pieces with a sterile scalpel and used for DNA extraction. DNA was extracted using DNeasy PowerSoil kit (Qiagen Inc., Valencia, CA, United States) with modifications to the manufacturer’s instructions (Supplementary Method) and were stored at −80°C until further processing. Sequencing libraries were prepared using the Accel-NGS 2S Plus DNA Library Kit (Swift Biosciences Inc., Ann Arbor, MI, United States) according to the manufacturer’s protocols. The prepared library was sequenced on the Illumina Novaseq platform at the Oklahoma Medical Research Foundation. Metagenome sequence quality was assessed using FastQC v0.11.5¹ and reads with a mean quality score less than 25 were discarded. Sequences were quality trimmed and Illumina sequencing adapters removed using Trimmomatic v0.36 (Bolger et al., 2014) with the following parameters: SLIDINGWINDOW:4:15, LEADING:3, TRAILING:3, and MINLEN:36. After quality control, taxonomic annotation of metagenomes was performed using Kaiju v1.8.2 (Menzel et al., 2016) against the NCBI nr + euk database (accessed on February 16, 2022) with the following parameters (run mode: greedy, minimum match length: 11, minimum match score: 75, allowed mismatches: 5). To compare the effects of taxonomic profiling using metagenome vs. 16S rRNA gene sequences, the trimmed metagenome sequences were uploaded to the MG-RAST online server (Meyer et al., 2008) and then passed through the MG-RAST QC pipeline. Following the quality control on MG-RAST, 16S rRNA gene sequences present in the metagenome (metagenome-derived 16S rRNA gene sequences) were identified and taxonomy assigned. Taxonomic classification of the metagenome-derived 16S rRNA gene sequences was done using the Greengenes rRNA database (DeSantis et al., 2006) hosted in MG-RAST with a minimum cut-off identity of 60% and e-value of 5. Since contamination leads to an overestimation of diversity, especially in low biomass samples like that of Solfatara and Pisciarelli (Karstens et al., 2019), after assigning taxonomy, we removed molecular biology kit and laboratory contaminants (which hereinafter are referred to as “kitome”), low abundance taxa (<0.01%) and sequences that were unclassified at the domain level from the taxonomy table before performing downstream analyses (Supplementary Method).

Bins were generated from 10 Solfatara and 17 Pisciarelli metagenomes. To maximize genome recovery, we co-assembled several samples based on their taxonomic composition. All 10 Solfatara metagenomes were co-assembled and all 17 Pisciarelli metagenomes were co-assembled using MEGAHIT v1.2.9 (Li et al., 2015). Contigs longer than 1,000 bp were then binned using the metaWRAP binning module (Uritskiy et al., 2018) which incorporates three binning methods: CONCOCT v1.1.0 (Alneberg et al., 2014), MaxBin2 v2.2.6 (Wu et al., 2016), and metaBAT2 v2.12.1 (Kang et al., 2019). The metaWRAP refinement module (Uritskiy et al., 2018) was used to merge results from the three binning methods using the -c 50 and -x 10 options to obtain bins with over 50% completeness and less than 10% contamination according to the CheckM tool v1.0.12 (Parks et al., 2015). Bins with >50% completeness and <10% contamination were then reassembled with SPAdes v3.13.0 (Bankevich et al., 2012) to improve the assembly quality. The contamination and completeness of resulting Metagenome Assembled Genomes (MAGs) were reassessed with CheckM (Parks et al., 2015).

Taxonomy was assigned to MAGs using GTDB-TK v1.7.0 (Chaumeil et al., 2020). To assess whether MAGs belong to the same species (species ANI ≥ 95%), average nucleotide identity (ANI) was calculated for each possible pair of MAGs using FastANI v1.33 (Jain et al., 2018). MAGs with taxonomic assignment similar to previously identified kitome (e.g., Corynebacterium and Staphylococcus; Sheik et al., 2018; Weyrich et al., 2019) were not included in downstream analyses. Functional annotation of MAGs, including gene prediction was done using Prokka v1.14.5 (Seemann, 2014). Predicted genes were compared against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using BlastKOALA server (Kanehisa et al., 2016) to obtain KEGG Orthology (KO) annotations.

To assess the abundance of recovered MAGs, quality trimmed sequences from each metagenome were mapped against each MAG using BBMap v38.96. Sequence counts were normalized as the number of sequences recruited per kilobase of MAG and gigabase of metagenome (RPKG). The normalized sequence counts allowed for direct comparison of genome abundance between metagenome of different depths.

All statistical analyses were performed in R v4.0.2 (R Core Team, 2020). Principal Component Analysis (PCA) was performed with environmental parameters using the “prcomp” function in the stat package with scaling enabled. PCA results were graphed using the “ggbiplot” function from the ggbiplot v0.55 package. Before using the geochemical data for statistical analyses, variables below detection limit in all samples were excluded and the dataset transformed using z-score. After removing kitome and all potential contaminants, we calculated alpha diversity using the species abundance table obtained from metagenome-derived 16S rRNA gene (Greengenes) and metagenome (NCBI nr + euk) taxonomic assignment. Alpha diversity (Shannon, observed richness and Pielou’s evenness) was calculated using the “estimate_richness” function from the phyloseq v1.32.0 package in R. Nonpareil v3.0 (Rodriguez-R et al., 2018) was then used to estimate coverage and also calculate diversity with kmer kernel and default parameters. Kruskal-Wallis test was used to investigate whether alpha diversity varied significantly between locations and substrates. The “lm()” function in R was used to perform linear regression to evaluate potential effects of environmental variables on alpha diversity. Bray–Curtis dissimilarity based on genus and species abundance was used to determine differences in microbial community composition. The Bray-Curtis distance calculated was visualized using a Non-metric Multidimensional Scaling (NMDS) plot. Analysis of Similarity (ANOSIM) was performed on the Bray–Curtis dissimilarity matrix using the vegan v2.5-6 package to evaluate the significance of microbial compositional differences between the two locations and different types of substrates. Similarity Percentage (SIMPER) analysis was performed using the vegan v2.5-6 package to identify genera that contributed to the dissimilarity between two locations and different types of substrates. Mantel test using Spearman’s correlation coefficient with 999 permutations was performed to evaluate the significance of correlation between community composition (Bray–Curtis dissimilarity) and environmental parameters (Euclidean distance).

At Solfatara-Pisciarelli CHSs (Figure 1), the crater structure that forms the rock substrate (alkaline potassic tephra and lava ranging from trachybasalt to phonotrachyte) for both locations was produced and shaped by the same volcanic events (Piochi et al., 2015). Solid substrates from both locations have similar bulk elemental composition [Supplementary Figure 1; Supplementary Table 1 and Piochi et al. (2014), Piochi et al. (2015)]. In addition, the same magma chamber provides a heat source for the fumaroles and mud pools (Chiodini et al., 2001; Valentino and Stanzione, 2003) of both systems. However, tectonic features diverge the fluids that feed Pisciarelli mud pools to pass through old marine sediments with strong organic imprints where the fluids become enriched in NH₄⁺ (Poichi et al., 2019).

PCA based on temperature, pH, Eh, and water-soluble nutrients showed further distinction between Solfatara and Pisciarelli (Figure 2; Supplementary Table 2). Pisciarelli an extremely acidic (pH: 1.5–2.5) and reducing (Eh: −628 to −485) environment has higher temperatures (74°C to 95°C) and gas concentrations including H₂S and CH₄ (Table 1) compared to Solfatara (pH: ~1; temperature: 42°C to 88.7°C; Eh: 198 and 330 mV), making Pisciarelli a more active and extreme hydrothermal environment than Solfatara. These higher temperatures, which inhibit the solubility of oxygen (Boyer et al., 2020), together with higher concentrations of reducing gases may contribute to the more reducing conditions we observed at Pisciarelli. The H₂S gas released in CHSs may be oxidized abiotically by oxygenated meteoritic waters (Nordstrom et al., 2005, 2009) or biotically by chemosynthetic microorganisms leading to the higher concentration of SO₄²⁻ we observed in mud pools (Pisciarelli: 3,326–6,208 ppm; Solfatara: 206–2,888 ppm) compared to fumarolic deposits (Pisciarelli: 118–439 ppm; Solfatara: 567–610 ppm). The concentrations of SO₄²⁻ in Pisciarelli mud pools were two times that of Solfatara mud pool possibly due to the higher concentration of H₂S gas measured at Pisciarelli. The concentration of NH₄⁺ was also higher in mud pools, particularly in Pisciarelli mud pools (Pisciarelli: 1,130–1,998 ppm; Solfatara: 1.1–36 ppm) than in fumarolic deposits (Pisciarelli: 0.8 to 38 ppm; Solfatara: 0.3 ppm) which corresponds with other studies that report higher concentrations of NH₄⁺ in Pisciarelli mud pools (Martini et al., 1991; Valentino and Stanzione, 2003; Glamoclija et al., 2004; Poichi et al., 2019). Temperature and geochemical measurements from Pisciarelli large mud pool revealed an environmental gradient; values were lower at the discharge channel (temperature: 79°C; SO₄²⁻: 134.4–193.6 ppm; NH₄⁺: 46.5–67.8 ppm) compared to the main bubbling pool (temperature: 84°C; SO₄²⁻: 3326–6,208 ppm; NH₄⁺: 1,130–1,998 ppm). The low concentrations of SO₄²⁻ and NH₄⁺ along the discharge channel may result from the removal of H₂S and NH₄⁺ possibly due to oxidation and volatilization as fluid flows into the channel.

In this study, a total of 7 Solfatara and 10 Pisciarelli samples were analyzed. Solfatara samples had between 54,932,932 to 210,855,916 (average = 126,261,781) quality filtered metagenome sequences whereas Pisciarelli had between 119,199,732 to 230,316,988 (average = 179,177,105) quality filtered metagenome sequences. We compared taxonomic profiles generated from metagenome against metagenome-derived 16S rRNA gene sequences and observed that metagenome sequences identified more taxa than metagenome-derived 16S rRNA gene sequences, which corresponds to other studies that detected an increased number of taxa with whole genome shotgun sequencing compared to the 16S amplicon method (Ranjan et al., 2016; Brumfield et al., 2020). For example, Ranjan et al. (2016) reported that with the same number of sequences, whole genome shotgun sequencing identified twice as many species as the 16S method. We observed that metagenome sequences identified 56 phyla while metagenome-derived 16S rRNA gene sequences identified 26 phyla. Only 3 out of the 26 phyla detected by metagenome-derived 16S rRNA gene profiling were not detected by metagenomic profiling. Metagenomic profiling, on the other hand, identified 35 phyla not identified by metagenome-derived 16S rRNA gene profiling with 12 of them being viral and eukaryotic phyla. Overall, the dominant bacterial and archaeal phyla detected were similar irrespective of profiling method. Predominant bacterial and archaeal phyla detected by both profiling methods include Crenarchaeota, Proteobacteria, Firmicutes, and Actinobacteria (Figures 3A,B). At the genus level, metagenomic profiling also identified more genera than metagenome-derived 16S rRNA gene profiling (Figures 4A,B). Of the 245 genera identified across all samples by metagenomic profiling methods only 48 genera (20%) were identified with both profiling methods. Metagenomic profiling identified 197 genera not identified by metagenome-derived 16S rRNA gene profiling including Acidilobus, Acidibacillus, Thermogymnomonas, Ampullavirus, and Bicaudavirus while genera including Acetobacterium and Caldococcus were unique to metagenome-derived 16S rRNA gene profiling. We observed that both methods identified dominant genera including Acidianus, Pyrobaculum, and Sulfobacillus.

Next, we analyzed the effect of both taxonomic profiling methods on diversity (Shannon diversity, observed richness, Pielou’s evenness). Before calculating the diversity, we identified and removed kitome (Supplementary Figure 3; Supplementary Tables 4, 5), low abundance taxa and sequences that were unclassified at the domain level from the taxonomy table to avoid an overestimation of diversity and a misrepresentation of community composition (Karstens et al., 2019). Our results showed that Shannon diversity and observed richness were consistently higher when the taxonomy table from metagenomic profiling was used as input compared to when we used the taxonomy table from metagenome-derived 16S rRNA gene profiling (Figures 5A,B). Overall, our results reflect that the 16S amplicon approach, which has been the most employed approach for studying CHSs’ microbiome, identifies a significantly lower number of bacterial species and completely excludes viruses and fungi which is similar to reports from other studies (Ranjan et al., 2016; Brumfield et al., 2020). Although metagenome sequencing permits the identification of more taxa, the choice between 16S rRNA gene and metagenome sequencing ultimately depends on the ecological questions and objectives of any given study. Since metagenome sequences encompass all members of the microbiome and provide high-resolution diversity analysis, which agrees with our objective to understand the full range of microbial ecology in available substrates from Solfatara-Pisciarelli CHSs, the further analyses and results hereinafter are based on the taxonomy table generated from metagenomic profiling.

Shannon diversity in Solfatara-Pisciarelli CHSs ranged from 1.20 to 4.26 (Figure 5B; Supplementary Table 8). We observed significant variations in Shannon diversity (Kruskal–Wallis: p < 0.05) between Solfatara and Pisciarelli which were supported by Nonpareil diversity values (Supplementary Figure 4; Supplementary Table 8); on average, Shannon diversity was lower at Pisciarelli (Shannon: 1.20–3.52) than at Solfatara (Shannon: 3.55–4.26). Additionally, microbial communities were more uneven (Kruskal–Wallis: p < 0.05) at Pisciarelli (Pielou’s evenness: 0.26–0.74) than at Solfatara (Pielou’s evenness: 0.67–0.84). However, there was no significant difference in observed richness between these locations (Solfatara: 114–201 species; Pisciarelli: 102–166 species; Kruskal–Wallis: p = 0.24). The decreased Shannon diversity at Pisciarelli may have resulted from the combined effects of reducing conditions, extremely low pH, and high temperature. The addition of highly reducing conditions at Pisciarelli to an already extreme environment (high temperature and low pH) may cause further selection against certain microbial groups by hindering their growth and metabolism, thus leading to the lower Shannon diversity observed (Kimbrough et al., 2006; Seo and DeLaune, 2010; Husson, 2013; Zhang et al., 2015). Alpha diversity in substrates showed some variation (e.g., on average, mud and fumarolic deposits had higher Shannon diversity than water, and richness showed less variability in water and mud); however, these differences were not statistically significant (Kruskal–Wallis, Shannon: p = 0.97; Pielou’s evenness: p = 0.83; observed genus richness: p = 0.09) even when substrates were analyzed separately for each location which implies that species alpha diversity does not differ across substrates in the Solfatara-Pisciarelli CHSs. This result contrast with other studies comparing microbial diversity in water and solid substrates from hot springs, which showed differences in richness (Colman et al., 2016) and evenness (Cole et al., 2013) across substrates.

A comparison of diversity and physical parameters showed pH (Shannon: R² = 0.73; Nonpareil: R² = 0.43, p < 0.05) to be the primary driver of microbial diversity which is consistent with several studies that report pH as the major predictor of microbial diversity in some hydrothermal systems (Oliverio et al., 2018; Power et al., 2018). Results also showed that Eh had a significant correlation (Shannon: R² = 0.62; Nonpareil: R² = 0.43, p < 0.05) with diversity. Temperature on the other hand had minimal effect on diversity (Shannon: R² = 0.12, p = 0.09; Nonpareil: R² = 0.02, p = 0.25). This is in contrast to previous studies that reported a strong relationship between temperature and microbial diversity in hot springs (Miller et al., 2009; Sharp et al., 2014). Comparably, minimal effect of temperature on alpha diversity has been reported for hot spring communities in Yellowstone, United States, Fludir, Iceland and Tibetan Plateau, China (Wang et al., 2013; Podar et al., 2020).

After removing kitome and all potential contaminants, a total of 490 species (245 genera) were identified across all samples corresponding to 11 archaeal, 31 bacterial, 10 eukaryotic, and four viral phyla (Supplementary Table 6). Our results showed that overall, 33% (133) of the species detected were shared across the different types of substrates but each substrate contained unique species among which water had the least number of unique taxa and mud contained the highest (Figure 6B). NMDS plot (Figure 7) generated based on Bray–Curtis dissimilarity matrix revealed that the microbial community structure of Pisciarelli was distinct from Solfatara (ANOSIM: R = 0.63 for species; R = 0.67 for genus, p < 0.05). On average, the most abundant phylum at Pisciarelli was Crenarchaeota (46%) followed by Proteobacteria (5%) and Actinobacteria (2%), whereas Solfatara was dominated by Firmicutes (23%), Proteobacteria (21%), and Actinobacteria (12%; Figure 3B). In addition, Pisciarelli substrates had a higher number of viral sequences (4–49%) than Solfatara (1–4%; Supplementary Figure 2). Interestingly, we detected Cyanobacteria (0.5 to 4%) especially at Solfatara even in samples (SF-D_and SF-D_b) with a temperature (89°C) that exceeds the known temperature limit for photosynthesis in acidic environments (Brock, 1985; Cox et al., 2011; Hamilton et al., 2012). The Cyanobacteria we detected may have been introduced into the hydrothermal system from surrounding cooler environments (by wind or fluid circulations). However, we cannot decipher with certainty if they were environmentally introduced or if they are true members of the Solfatara-Pisciarelli hydrothermal community. We observed that at the species level, both locations shared 171 species (35%), however, there were over 100 species unique to each location (Figure 6A). Species unique to Pisciarelli include thermophiles such as Pyrolobus fumarii, Pyrobaculum aerophilum and Metallosphaera sedula. Unique species identified at Solfatara include Thermaerobacter sp. FW80, Methylacidiphilum kamchatkense and Ilumatobacter coccineus.

Our results showed that the microbial community structure in the high temperature environments (74–95°C) of Pisciarelli was different. The source of Pisciarelli large mud pool (84°C) was dominated by Acidianus (21–68%), Pyrobaculum (11–21%) and Bicaudavirus (11–43%) but as temperature decreased away from the source (79°C), we detected a decrease in the abundance of Acidianus (20–23%) and Bicaudavirus (15–17%) but an increase in the abundance of Pyrobaculum (30–43%) in the discharge channel. Water from Pisciarelli small mud pool (89°C) was also dominated by Acidianus (48–59%), but Pyrobaculum accounted for <1% of total sequences in this mud pool. Interestingly, the viral genus, Ampullavirus accounted for over 80% of sequences in Pisciarelli epilithic microbial layer (74°C) and for 26% of sequences in both Pisciarelli fumarolic deposit (93°C) and mud from Pisciarelli small mud pool (88°C). Most of the viruses detected in Pisciarelli fumarolic deposits, epilithic microbial layer and mud from the small mud pool were known archaeal viruses (e.g., Ampullavirus, Bicaudavirus, and Betalipothrixvirus) even though archaea accounted for <3% of the total sequences in these samples. This difference in relative abundance between archaeal viruses and their known host suggests that each archaeal species may host more than one virus type or the archaeal viruses present in Solfatara-Pisciarelli CHSs may have a broader host range that includes bacteria or/and eukaryotes (Munson-McGee et al., 2018). When substrates were compared at each location, we observed that Pisciarelli substrates harbored significantly different (ANOSIM species: R = 0.64; genus: R = 0.62, p < 0.05) microbial communities. SIMPER analysis showed that the sulfur-oxidizing archaeal genus Acidianus was significantly (p < 0.05) more abundant in Pisciarelli water communities, the archaeal genus Pyrobaculum showed a higher relative abundance in mud communities, and Ampullavirus showed a higher relative abundance in fumarole deposits communities (Supplementary Table 9; Figure 4B). Overall, Pisciarelli water and mud communities were similar but distinct from communities in Pisciarelli fumarolic deposits. The distinction between communities in the mud pools and fumarolic deposits at Pisciarelli may be linked to fluids and/or may reflect the availability of oxygen. Although the entire Pisciarelli environment has high temperatures and is highly reducing, the vigorous mixing of hydrothermal gas with air may create pockets of atmospheric O₂ in pore spaces of fumarolic deposits which support aerobic communities that may be absent or in low abundance in mud pools. In contrast, high temperatures of mud pools inhibit the dissolution of oxygen thereby creating microaerobic environments that may favor thermophilic facultative and obligate anaerobes such as Acidianus and Pyrobaculum. Furthermore, the presence of Pyrobaculum specifically Pyrobaculum arsenaticum which is a strict anaerobe (Huber et al., 2000) and members of the genus Sulfuracidifex who are all obligate aerobes in Pisciarelli large pool may be an indication that conditions in this pool fluctuate between microaerobic and completely anaerobic, possibly in congruence with the intensity of hydrothermal activity.

At Solfatara, we observed variability in genus abundance across samples even between replicates (Figure 4B). For example, Acidibacillus accounted for 11% of sequences in SMP-M_a1 but was <1% in replicate mud sample, SMP-M_b. The difference between replicates may derive from the heterogeneous nature of Solfatara-Pisciarelli CHSs; however, seeing that genus abundance also varied in our sequencing replicates (SMP-M_a1 and SMP-M_a1.2), the variability in genus abundance may likely be bias introduced from sequencing. Results also showed that a large percentage of Solfatara sequences (20–36%) were unclassified at the genus level, suggesting a considerable amount of potentially novel genera at this environment. Interestingly, at Solfatara, we detected several genera that are not known thermophiles and not commonly found at pH < 2 (e.g., Rhizopus; 1.7–14%, Actinomadura; 1–9%, Terasakiella; 0.7–5.6%, and Clostridioides; 0.3–2.5%; Wang et al., 2014; Strazzulli et al., 2020). These genera were present in Solfatara water (43°C), mud (69°C), and fumarolic deposits (89°C). Furthermore, we detected these genera in most Pisciarelli samples but to a lesser degree. These microorganisms may not be native to the hydrothermal systems but likely introduced into the habitat from extensive human activity around the crater. In addition, the fluids in mud pools are a mixture of hydrothermal fluids and meteoric water which could deliver genetic material from non-native microorganisms to the hydrothermal systems. This possible introduction of non-native microorganisms into CHSs may be a potential problem for microbial ecology studies since the very nature of the CHSs makes it difficult to completely exclude all environmental contaminants. In contrast to Solfatara water and fumarolic deposits communities, Solfatara mud communities was dominated by Sulfobacillus (10–34%) and Acidibacillus (4–20%). However, unlike Pisciarelli, the differences across Solfatara substrates were not significant (ANOSIM species: R = 0.03, p = 0.50; genus: R = 0.03, p = 0.43) suggesting that variations in microbial community structure across substrates is not a general phenomenon but is specific to hydrothermal systems and possibly linked to the physical and geochemical conditions of individual systems (Wang et al., 2014). The microbial community structure observed in our study is in contrast with studies that investigated the microbial ecology of Solfatara-Pisciarelli CHSs and identified Acidithiobacillus and Metallosphaera as the dominant genus in Solfatara and Pisciarelli, respectively (Iacono et al., 2020; Strazzulli et al., 2020). Furthermore, the dominance of archaea which we observed in Pisciarelli mud pools is in contrast with a study that reported a low archaea/bacteria ratio in high temperature waters (>85°C) of Pisciarelli (Crognale et al., 2022). This contrast between our study and other Solfatara-Pisciarelli studies is likely the consequence of a microbial diversity that reflects differences in sampling year/periods from active CHSs with variable geochemistry. Additionally, the dissimilarity may be related to differences in sampling points (10 cm depth vs. 1 cm depth), sample handling, DNA extraction protocol, bias introduced by 16S rRNA gene primers or sequencing technology. An assessment of the effects of environmental variables on microbial community composition showed that pH and Eh (Mantel: ρ = 0.53 and ρ = 0.47, respectively, p < 0.05) had the strongest correlation to beta diversity which is consistent with reports from other studies that pH and redox potential shape microbial community composition (Alsop et al., 2014; Power et al., 2018). The concentrations of NH₄⁺ and SO₄²⁻ (Mantel: ρ = 0.38 and ρ = 0.23, respectively, p < 0.05) also had significant correlations to beta diversity. However, temperature had no significant (Mantel: ρ = 0.1, p = 0.14) effect on the microbial community composition of Solfatara-Pisciarelli CHSs, which coincides with temperature trends observed for alpha diversity.

To identify the metabolic potential of members of the microbial community and their connection to geochemistry, we generated MAGs. After co-assembly and binning, a total of 22 MAGs (10 from Solfatara and 12 from Pisciarelli) with >50% completeness and <10% contamination were recovered. Out of the 22 MAGs recovered, 16 had completeness >90% and contamination <5%. Taxonomic classification revealed that MAGs belonged to three bacterial (Firmicutes, Proteobacteria, and Aquificota) and 2 archaeal (Thermoproteota also known as Crenarchaeota and Thermoplasmatota) phyla (Supplementary Tables 10, 11), with over 40% (9 out of 22) of unclassified genomes at the genus level suggesting that our MAGs represent novel genera. MAGs belonging to the phylum Thermoproteota were most represented in Pisciarelli water and mud while Firmicutes were more abundant in all Solfatara substrates and Pisciarelli fumarolic deposits (Figure 8A), which is consistent with our results from metagenome fragments. On average, the most abundant MAGs in Pisciarelli water and mud belong to species of the Acidianus and Pyrobaculum genera and an unknown species of the family Sulfolobaceae whereas the most abundant MAGs in fumarolic deposits belong to 2 unknown species: one species of the Alicyclobacillaceae family and the other species of the Alicyclobacillales order and a third species belonging to the Acidianus genus (Supplementary Table 11). In contrast, the most abundant MAG in all Solfatara samples belong to species of the Acidianus genus. When MAGs were dereplicated based on FastANI of ≥95% we found that only 2 species (Acidianus infernus; and an unknown species of the Alicyclobacillaceae family) were shared between locations further confirming that although an overlap exists between both locations; they have different microbial community structures. Since acidic hydrothermal environments are dominated by chemolithoautotrophic communities that derive energy from inorganic compounds, we analyzed MAGs for marker genes involved in carbon fixation, sulfur, nitrogen, and methane metabolism.