This study focuses on six selected hot springs, located along the northeastern foothills of the Sierra Madre Oriental (SMO) in the Mexican states of Nuevo León, Tamaulipas, and San Luis Potosí, Mexico (Figure 1). The springs examined include: Potrero del Prieto (PP), El Bañito (EB), Ojo Caliente (OC), Taninul (TA), Balneario El Bañito (BEB), and Mainero Azufroso (MA). Notably, Mainero Azufroso, although exhibiting a cooler temperature similar to that of the adjacent stream, releases hydrogen sulfide vapors that contribute to the formation of unique aquatic ecosystems characterized by the presence of “green mats” and localized precipitation of sulfur, gypsum, calcite, or halite, warranting its inclusion in this study.

The PP hot spring is situated near the Prieta Linda waterfall and the town of El Potrero del Prieto de Arriba, nestled between the Iturbide anticline and the El Mezquital syncline in the Sierra El Baño, on the bed of the Cabezones River, at an elevation of 1,229 meters above sea level (masl). This spring is recognized as the highest hot spring in the SMO. The closest meteorological station -19073 Galeana- (Servicio Meteorológico Nacional (Mexico), 2021), reports an annual precipitation of 361.8 mm, and the region is classified as having a dry semi-warm climate (http://www.inegi.org.mx/temas/climatologia/). This hot spring emerges from the Lower Tamaulipas Formation of the Lower Cretaceous and currently has no designated use.

The EB hot spring is situated at the central part of the anticline of the sierra Cerro de La Silla, near the Rodriguez Gómez dam and the La Chueca creek at an elevation of 405 masl. The closest meteorological station -19069 La Boca- (SMN, 2021), reports an annual precipitation of 1001 mm, and the region is classified as having a semi-warm climate with summer rain (www.inegi.org.mx/temas/climatologia/). This hot spring emerges from the La Casita Formation of the Upper Jurassic and currently is used for recreational purposes.

The OC hot spring is located at an altitude of 364 masl, at the base of the El Platanillo mountain range in the Sierra El Filo. The nearest meteorological station, -28218 La Boca- (Servicio Meteorológico Nacional (Mexico), 2021), documents an average annual precipitation of 743 mm, with the area experiencing a temperate subhumid climate (http://www.inegi.org.mx/temas/climatologia/). This spring originates in the San Felipe Formation of the Upper Cretaceous and is situated on private property, currently utilized for livestock.

The TA hot spring is well-regarded locally for its medicinal and recreational applications, possibly dating back to pre-Hispanic times. It is part of a hotel complex and is located in the Sierra El Abra-Tanchipa, at the foot of the El Abra mountain range, emerging from the Cretaceous El Abra Formation at an elevation of 64 masl. The nearest meteorological station -3145 El Choy- (Servicio Meteorológico Nacional (Mexico), 2021) indicates an annual precipitation of 1165.4 mm, with the region characterized by a warm subhumid climate (http://www.inegi.org.mx/temas/climatologia/).

The BEB hot spring is currently used for recreational purposes within the municipality of Ciudad Valles. It is situated on a gently sloping hilltop, at 55 masl. According to the Ciudad Valles -24012- meteorological station (Servicio Meteorológico Nacional (Mexico), 2021) annual precipitation is reported at 1241.2 mm. This spring emerges from the San Felipe Formation, existing within a warm subhumid climate (http://www.inegi.org.mx/temas/climatologia/).

The MA spring rises at 715 masl, along the bed of an intermittent stream in the Sierra La Guitarra, at the base of the San Manuel Mountain range. The Villa Mainero −3735- meteorological station (Servicio Meteorológico Nacional (Mexico), 2021), records an annual precipitation of 993.8 mm, and the area is characterized by a temperate subhumid climate (http://www.inegi.org.mx/temas/climatologia/). This spring originates from the Taraises Formation of the Lower Cretaceous and currently lacks a specific use.

At each spring location, a comprehensive assessment was conducted involving the measurement of 11 physicochemical variables: temperature, dissolved oxygen, electrical conductivity, salinity, oxidation–reduction potential, pH, turbidity, alkalinity, OH⁻, CO₃²⁻, and HCO₃⁻. The methodologies employed for the precise measurement of each variable are detailed in Pantoja-Irys et al. (2025).

Water samples were systematically collected directly from the spring source to minimize any influence from external water flows and to maintain the integrity of the native microbial communities. At each designated sampling site, five 1-liter replicates were collected using sterile plastic bottles. Upon filling, the bottles were promptly placed into a cooler to sustain a low temperature and inhibit microbial growth prior to the filtration process.

Filtration was conducted to effectively isolate and concentrate the microorganisms present in the water samples. This was achieved using sterile cellulose ester membranes with pore sizes of 0.45 μm, which are capable of retaining bacteria and other microbial cells. The filtration equipment consisted of a filtration funnel, a Kitazato flask connected to a vacuum pump, and sterilized membrane filters. During the assembly of the equipment, the membranes were meticulously positioned in the funnel using sterile gloves and dissecting forceps to prevent any risk of contamination.

The funnel was securely attached to the Kitazato flask, which was connected to a vacuum pump to establish negative pressure. Subsequently, the membranes were transferred to sterile 10 mL Falcon tubes, appropriately labeled and stored at 4 °C to preserve the DNA until subsequent analysis.

The filtered samples were shipped to MR DNA Laboratory (Shallowater, Texas, United States) for comprehensive microbial community profiling utilizing 16S rRNA gene amplicon sequencing. DNA was extracted using proprietary MR DNA protocols specifically optimized for environmental samples. The V4 region of the 16S rRNA gene was amplified using the universal primer pair 515F (GTGYCAGCMGCCGCGGTAA) and 806R (GGACTACNVGGGTWTCTAAT). PCR amplification was conductedperformed in a single-step reaction using the HotStarTaq Plus Master Mix Kit (Qiagen, USA) withemploying the following thermocycling conditions: an initial denaturation at 95 °C for 5 min; 30 cycles of denaturation at 95 °C for 30 s, annealing at 53 °C for 40 s, and extension at 72 °C for 1 min; followed by a final extension at 72 °C for 10 min. The PCR products were visualized on 2% agarose gels, pooled in equimolar ratios based on concentration and molecular weight, and purified using calibrated SPRI (Solid Phase Reversible Immobilization) beads. Sequencing was performed on the Illumina NovaSeq 6000 platform using paired-end chemistry (2 × 250 bp) according to the manufacturer’s protocols.

The sequence data were processed the MR DNA proprietary bioinformatics pipeline and QIIME2 v2023.2 (Bolyen et al., 2019). Paired-end reads were joined, and sequences shorter than 150 bp or containing ambiguous base calls were removed. Primer sequences were trimmed using Cutadapt, and reads were quality-filtered using a maximum expected error threshold of 1.0. Unique sequences were dereplicated and denoised using the UNOISE3 algorithm to generate amplicon sequence variants (ASVs), also referred to as zero-radius operational taxonomic units (zOTUs). Chimeras were detected and removed with UCHIME (Edgar et al., 2011). Taxonomic classification of ASVs was performed using BLASTn against a curated version of the NCBI nucleotide database (Edgar, 2010). The final outputs included absolute abundance tables and relative abundance matrices at various taxonomic ranks from phylum to species, as well as zOTU-to-sample mapping files (Edgar, 2016).

In order to assess species diversity within six hot spring communities, species richness and diversity were estimated using Hill numbers of order q = 0, 1 and 2. Hill numbers offer a comprehensive framework for quantifying diversity based on the effective number of species, with variations contingent on the parameter q. This parameter allows for an adjustment in the weighting of species abundances, effectively reflecting diverse aspects of community structure and enabling comparisons across samples with distinct dominance patterns (Jost, 2006; Moreno et al., 2011). Specifically, species richness (q = 0) quantifies the total number of species, without regard to their abundances. Shannon diversity (q = 1) provides a balanced estimate that incorporates both richness and evenness, thereby moderately weighting species according to their relative abundances. In contrast, Simpson diversity (q = 2) places greater emphasis on the most abundant species, consequently reducing the influence of rare species and effectively highlighting patterns of dominance within the community.

Inventory completeness was standardized through the use of sample coverage (Ĉn) facilitating a meaningful comparison of spring diversity, across various communities and ensuring that all were analyzed at consistent level of sampling completeness (Chao and Jost, 2012; Chao and Hsieh, 2016). Effective diversity estimates and sample coverage were derived utilizing the “iNEXT” function from the iNEXT package in R (https://www.r-project.org/). Statistical comparisons were conducted on the 95% confidence intervals of the Hill numbers, where significant differences were inferred if the 95% confidence intervals did not overlap.

Additionally, spring communities were classified based on species composition using a cluster dendrogram constructed using the Bray–Curtis dissimilarity index and Ward’s agglomeration method. Given the sensitivity of the Bray–Curtis index to species abundances, bacterial species abundances were log-transformed (x + 1) prior to analysis to mitigate this influence and to achieve a balanced representation of both common and rare species. The analysis was performed using the “hclust” function from the vegan package in R.

A Principal Component Analysis (PCA) was performed to identify the physicochemical variables associated with the variance among community groups observed in the studied springs. Prior to conducting the analysis, we examined the correlations among variables, identifying pairs with high correlation coefficients (> 0.85). Representative variables were selected and subsequently log-transformed, with the exception of pH, which is inherently expressed on a logarithmic scale. The PCA facilitated a visualization of the influence of environmental variables on the spring groups, allowing us to and to discern key environmental gradients. The analysis was executed using the “rda” function from the vegan package in R, resulting in a two-dimensional ordination plot. In this plot, community groups are represented as points, while physicochemical variables are denoted as vectors (arrows). The length and direction of each vector indicate the magnitude and direction of influence exerted by the respective variable. Groups situated near the terminus of the vectors exhibit a strong association with corresponding environmental gradients, whereas those positioned closer to the origin exhibit reduced influence from the measured variables.

Additionally, we estimated the indicator values of species within the community groups to identify the most robust and ecologically relevant indicator species for each group. To enhance the reliability of this analysis, we filtered the dataset used in the clustering analysis to include only 168 species, which collectively represented 95% of the total abundance. This approach allowed us to concentrate on species that significantly contribute to community structure. The indicator value for each species was calculated using the method proposed by Dufrêne and Legendre (1997) known as IndVal. This method quantifies both the specificity and fidelity of each species to a particular habitat, where specificity denotes the exclusivity of a species to a given group, and fidelity refers to the frequency of occurrence within that group. Indicator values were calculated independently for each taxon and expressed as percentages, yielding a robust metric for assessing species associations with various spring types. The analysis was conducted using the statistical software PAST version 4.17 (https://www.nhm.uio.no/english/research/resources/past/).

The TA hot spring exhibited the highest temperature and turbidity among the studied springs, while the PP spring demonstrated the highest salinity levels. Although no significant thermal anomaly was identified in the MA spring, it recorded the highest concentration of dissolved oxygen (Table 1). All springs maintained a neutral pH, and generally exhibited high alkalinity, predominantly in the form of bicarbonate, along with a negative redox potential, with the exception of the OC spring (Table 1).

The microbiome of the six hot springs comprised 425 species, with the majority (409) attributed to the domain Bacteria, alongside 13 to Eukarya, and 3 to Archaea. The bacterial dataset encompassed 177,136 counts across the six springs, distributed among 31 phyla, 43 classes, 73 orders, 138 families, 245 genera, and 409 species. Notably, the highest bacterial counts were recorded in the TA, BEB, and PP springs, with 32,718, 31,580, and 30,921 reads, respectively. Conversely, the MA, the OC, and the EB springs exhibited the lowest counts, with 26,301, 27,803, and 27,813 reads, respectively.

The bacterial communities within the MA and the OC springs displayed the highest taxonomic richness across all taxonomic levels, from phylum to species. In contrast, the TA and the BEB springs exhibited the lowest overall taxonomic representation. While the BEB and the PP springs contained a greater number of sequence-assigned species compared to the TA spring, both communities exhibited lower richness at the phylum level, indicating a more constrained phylogenetic breadth. It is noteworthy that despite a similar number of species in the MA and OC springs, the MA spring was characterized by a higher richness at the phylum level, suggesting a broader phylogenetic diversity (Table 2).

The phylum Pseudomonadota was found to be the most abundant across all springs with the exception of the BEB and the TA springs, where Campylobacterota and Chlorobiota, respectively, dominated. Additionally, Bacteroidota ranked as the second most abundant phylum within spring EB (Figure 2).

The taxonomic diversity observed in the investigated ecosystems revealed several dominant classes, notably Gammaproteobacteria, Betaproteobacteria, Epsilonproteobacteria, Chlorobiia, Flavobacteriia, and Alphaproteobacteria. Epsilonproteobacteria was particularly prevalent in the BEB spring, representing 80.2% of the total microbial abundance. Similarly, Gammaproteobacteria showed significant abundance in the MA spring, accounting for 71.5% of the total, while Flavobacteriia was the predominant group in the EB spring, contributing 40% to the overall abundance. It is noteworthy that Chlorobiia was exclusively identified in the TA spring, where it constituted 36.3% of the total microbial population.

At the order level, key representative groups included Campylobacterales in both the PP and BEB springs; Chromatiales and Chlorobiales in the TA spring; Flavobacteriales in the EB spring; Chromatiales again in the MA spring. In the OC spring, the dominant orders were and Burkholderiales, Moraxellales, and Rhodobacterales in the OC spring.

A detailed analysis of representative families and genera across the springs reveals that in the PP spring, Arcobacteraceae and Chromobacteriaceae were predominant, with the genera Halarcobacter and Vogesella; the BEB spring was characterized by Sulfurovaceae and Thiovulaceae, with Sulfurovum and Sulfuricurvum identified; the TA spring exhibited Halothiobacillaceae and Chlorobiaceae, featuring Thiofaba and Chlorobaculum; in the EB spring, Flavobacteriaceae, particularly Flavobacterium was prevalent; while the MA spring showcased, Halothiobacillaceae, represented by Thiofaba; finally, in the OC spring was marked by Comamonadaceae and Moraxellaceae, including Limnohabitans and Acinetobacter.

The sample coverage across all six springs surpassed 99%, with values ranging from 0.9985 to 0.9997; (Table 3). Notably, species richness (q = 0) demonstrated significant variation among springs. Springs OC and MA displayed the highest levels of bacterial richness compared to the other sites. While no statistically significant differences in species richness were detected between the PP and the EB springs, both sites exhibited greater richness compared to the BEB and the TA springs. The TA spring recorded the lowest level of richness (Table 3; Figure 3).

In terms of Shannon diversity (q = 1), the OC spring exhibited the highest value among all studied springs, with a q1 = 31.36. This was followed, in descending order, by the MA, the PP, and the EB springs (Table 3; Figure 3). In contrast, the BEB and the TA springs, which demonstrated similar diversity values, recorded the lowest levels of diversity, indicating a reduced evenness in species abundance. For inverse Simpson diversity (q = 2), the OC spring again represented the apex of diversity, achieving a q2 = 11.13, while the BEB spring marked the lowest at q2 = 2.83, suggesting a pronounced dominance of a select few species. The PP, the EB, and the MA springs displayed moderately high diversity values, reflecting varied degrees of evenness in species distribution (Table 3). Notably, despite the TA spring’s relatively low species richness, its abundances were distributed more evenly than those observed in the BEB spring (Table 3).

The classification of springs based on species composition revealed five distinct groups (Figure 4). Group I was exclusively comprised of the TA spring, while the BEB, EB, and PP springs formed separate clusters, designated as Groups II, III, and IV, respectively. The OC and MA springs were classified together in Group V, exhibiting the highest similarity in bacterial community composition (51%).

The Principal Component Analysis (PCA) revealed that the first two components accounted for 92.2% of the total variation in the dataset (Figure 5). The first principal component (PC1) accounted for 73.9% of the variability and was significantly influenced by electrical conductivity (EC) and bicarbonates (HCO₃⁻). The second principal component (PC2) accounted for 18.3% of the variation, with dissolved oxygen (DO) and pH contributing to its high values. Groundwater temperature (T°gw) had a moderate influence on PC2, although to a lesser extent compared to the other key other variables (Figure 5).

Group I demonstrated a strong association with elevated levels of HCO₃⁻ and EC along PC1. Group II was positioned near the center of the plot, indicating intermediate values across the majority of physicochemical variables. Conversely, Group III exhibited a pronounced correlation with pH and DO, suggesting conditions characterized by higher oxygen availability and a neutral to basic pH. Group IV displayed some relationship with T°gw, although no other physicochemical variable appeared to exert a dominant influence. Lastly, Group V was located near the origin, indicating more balanced or less differentiated physicochemical characteristics across the measured variables (Figure 5).

A total of 98 species were identified as having indicator potential (IndVal ≥ 50%), with 30 species demonstrating statistical significance as indicators for the various groups (Table 4). Notably, species classified as indicators in each spring exhibited over 97% homology with existing sequences.

In Group I, 13 species showed indicator potential; however, only two species, Enterococcus dispar (IndVal = 45.81%, p = 0.0492) and Lactococcus garvieae (IndVal = 63.97%, p = 0.0176) reached statistical significance.

Group II comprised four species with indicator potential, three of which were statistically significant (Table 4). Among these, Dyella ginsengisoli (IndVal = 67.93%, p = 0.0358) and Serratia marcescens (IndVal = 59.65%, p = 0.018) were particularly noteworthy.

Group III exhibited the highest number of significant indicator species, with 21 species showing indicator potential, of which 16 were statistically significant (Table 4). Strongly associated species in this group included Luteolibacter yonseiensis (IndVal = 79.06%, p = 0.0179), Paracidovorax avenae (IndVal = 74.2%, p = 0.0179), and Imtechium assamiensis (IndVal = 72.74%, p = 0.0179).

Group IV revealed a considerable number of potential indicators, totaling 20 species, among which nine were statistically significant. Highlighted species included Lysobacter daecheongensis (IndVal = 61.53%, p = 0.0178), Erythrobacter donghaensis (IndVal = 61.27%, p = 0.0255), and Luteimonas lutimaris (IndVal = 55.75%, p = 0.035).

Lastly, Group V comprised 40 species with indicator potential; however, no species were identified as statistically significant indicators.