https://orcid.org/0000-0002-4249-0488
Abstract
The Cuatro Ciénegas Basin is a biodiversity hotspot known for its unique biodiversity. However, this ecosystem is facing severe anthropogenic threats that are drying its aquatic systems. We investigated microbial communities at three sites with different physicochemical and environmental characteristics (Pozas Rojas, Archean Domes, and the Churince system) within the basin to explore potential connections to deep aquifers and determine if the sites shared microorganisms. Utilizing 16S rRNA gene data, we identified a core microbiota between Pozas Rojas (PR) and Archean Domes (AD). Sulfur reduction appears to shape the microbial connectivity among sites, since sulfur-reducing bacteria has the highest prevalence between samples from PR and AD: Halanaerobium sp. (88.46%) and Desulfovermiculus halophilus (65%); and between the Churince system and AD: Halanaerobium sp. (63%) and D. halophilus (60%). Furthermore, metagenome-assembled genomes from Ectothiorhodospira genus were found in both AD and Churince, suggesting microbial dispersal. An important finding is that microbial diversity in the AD system declined, from 2016 to 2023 the ecosystem lost 29 microbial phyla. If this trend continues, the basin will lose most of its water, resulting in the loss of various prokaryotic lineages and potential biotechnological solutions, such as enzymes or novel antibiotics. Our findings highlighting the need for water extraction regulations to preserve the basin’s biodiversity.
Introduction
The Cuatro Ciénegas Basin (CCB), located in the Chihuahuan Desert of Coahuila, Mexico, is a highly diverse oasis. Renowned for its unique and fascinating biodiversity, the CCB harbors diverse animals, plants, and microbes, including stromatolites and microbial mats. The Basin features a complex system of ponds, lakes, and rivers interconnected by underground water sources (Wolaver et al. 2012, Souza et al. 2018), creating a “lost world” of microbial diversity (Souza et al. 2006, 2018, 2023).
Sadly, the CCB is now endangered (Souza et al. 2023), as evidenced by the complete desiccation of the Churince system in 2011 (Souza et al. 2018) and the loss of all surface water in the Archean Domes (AD) system by 2023 (Souza et al. 2023).
The CCB, harbors interesting—and apparently ancient and unique—microbial lineages (Moreno-Letelier et al. 2011, Souza et al. 2018). It has been suggested that CCB microbial diversity survived within the deep aquifer of the Sierra de San Marcos y Pinos (see Fig. 1A and B), creating a “lost world” of microbial diversity dating back to the late Precambrian (Souza et al. 2018). One remarkable species that has only been found in CCB is Bacillus coahuilensis (Alcaraz et al. 2008, Torres et al. 2018, Moreno-Letelier et al. 2012), an ancient lineage that split from an ancestor of most Bacilli ca. 1000 million years ago (Moreno-Letelier et al. 2011), which has the capability of synthesizing membrane sulfolipids, a trait probably acquired through horizontal gene transfer from cyanobacteria.
(A) Map of Coahuila, northern Mexico, with a marker indicating the location of the CCB. (B) Aerial view of the CCB. (C) Distribution of the Pozas Rojas ponds within the CCB. Images obtained from Google Earth (2024). Source: Google Earth, Maxar Technologies.
Microbial beta diversity within CCB aquatic systems is high, even at small spatial scales (Espinosa-Asuar et al. 2022), suggesting that each sample represents a larger “seed bank” community in the deep aquifer (i.e. a preserved microbial community within the groundwater potentially dispersing throughout the entire CCB; Souza et al. 2018). Despite water overexploitation (Pérez-Ortega 2020, Leal Nares et al. 2022), remnants of the deep aquifer may still safeguard microbial communities from extinction (Souza et al. 2008), enabling their recovery through appropriate water management policies for agriculture and the closure of canals.
In ecosystems lacking light, such as groundwater, microbial communities are expected to remain stable over time yet differ between sites due to environmental heterogeneity (Griebler and Lueders 2009). Groundwater can support high microbial diversity, shaped by ecosystem characteristics (e.g. rock type), microbial interactions, biogeochemical cycles, and recently by human disturbance (Fillinger et al. 2023). It has been shown that microbial dispersion between groundwater and surface ponds is possible (see Ji et al. 2022). However, the potential for microbial dispersion across freshwater ponds (whether oligotrophic or not) via a shared deep aquifer remains largely unexplored in the literature.
Among the numerous aquatic systems and ponds in CCB that have been extensively studied, the AD system (Cisneros-Martínez et al. 2023, Espinosa-Asuar et al. 2022, Madrigal-Trejo et al. 2023, Rodriguez-Cruz et al. 2024), Churince (de Anda et al. 2018a, Souza et al. 2018, 2023), and Pozas Rojas (PR) (Peimbert et al. 2012, Garcia-Ulloa et al. 2022), are the best understood.
AD is an extreme environment characterized by high methane levels and elevated salinity, that presents dome-shaped flexible microbial mats. During the rainy season, the pH and salinity at the AD site were 9.8% and 5.28%, respectively, whereas in the dry season, the pH dropped to 5, while salinity reached saturation. A C:N:P imbalance of 122:42:1 was observed during the rainy season (Medina-Chávez et al. 2023). AD is the only CCB site studied so far that harbors high Archaea diversity (Medina-Chavez et al. 2023, Rodriguez-Cruz et al. 2024), with most taxa being poorly known and representing phylogenetically basal lineages (Rodriguez-Cruz et al. 2024). In contrast, the Archaea domain is underrepresented at the other studied CCB sites. For instance, previous metagenomic diversity profiles of microbial mats at CCB, such as the red mat located in the Los Hundidos region (near PR site, at 26° 52′ 17″ N, 102° 01′ 11.3″ W) and the green mat, at Pozas Azules Ranch (near AD site, at 26° 49′ 24.4″ N, 102° 00′ 53.2″ W) (Bonilla-Rosso et al. 2012), revealed that Archaea represented only 1.78% of the red mat, and 2.06% of the green mat (bacteria, accounting for the rest).
Similarly, a study of the Churince system (Souza et al. 2018) using 16S rRNA gene tags, detected 5167 OTUs (Operational Taxonomic Units) across multiple sampling sites, representing 60 different phyla, of which only three were Archaea. This evidence underscores the significance of the high Archaea diversity found at AD, compared with their underrepresentation in other CCB sites.
The Churince system, located in the western region of the CCB, was extensively studied from 2000 until 2011, when it dried out (Souza et al. 2018, 2023). It was considered an interesting system due to its strong gradients in salinity, temperature, pH, and dissolved oxygen (Cerritos et al. 2011). Additionally, it was surrounded by large and pure gypsum dunes (Johannesson et al. 2004). The environment at Churince was extremely poor in phosphorus (concentrations of PO−4 lower than 0.1 μmol), but rich in sulfate and magnesium (Elser et al. 2005).
On the other hand, the PR system, located on the eastern side of the CCB, is considered the best-preserved system (García-Ulloa et al. 2022). It consists of several small ponds surrounding a deeper lagoon and is characterized by an extremely unbalanced stoichiometry (C:N:P ratio of 820:157:1) (Bonilla-Rosso et al. 2012, Peimbert et al. 2012). However, its extreme oligotrophy was disrupted by Hurricane Alex in April 2010, which introduced relatively high concentrations of nitrogen and phosphorous (García-Ulloa et al. 2022). This site retains dynamic interactions between the deep aquifer and surface water, with microbial mats remaining intact. In soil microbial communities, natural sites harbor greater compositional diversity than anthropized sites, where microbial diversity tends to be lower (Mo et al. 2024). Similarly, in mangrove sediments, conserved mangroves exhibit higher microbial diversity and distinct signature taxa (Esguerra-Rodriguez et al. 2024).
We observed a progressive decline in the wetland at AD over the 8 years of our study, likely driven by poor management of the deep aquifer in the area, where alfalfa is the main crop (Pérez Ortega 2020). This decline is comparable to the total desiccation of the Churince system by 2011, as mentioned above. The reduction in water availability appears to have led to a significant loss of microbiological diversity, including unique microorganisms, and what is likely one of the most diverse archaeal communities in the world (see Rodríguez-Cruz et al. 2024).
We hypothesized that the natural water circulation between the deep aquifer and the surface water in the CCB, as described by Wolaver et al. (2012) (see also Souza 2012, 2018, Souza and Eguiarte 2018), promotes the dispersal of microbiological diversity across the basin. As such, we predicted the presence of core microbial taxa shared among different ponds, despite their different surface environmental conditions and geographic distance.
To date, no evidence supports the transport of microbial diversity between ponds in the CCB mediated by animals such as insects or birds (Souza et al. 2006). Current understanding suggests that shared environmental factors, such as the potential influence of a deep aquifer, likely explain the observed microbial connectivity. This aquifer may also serve as a common source facilitating dispersal across the CCB, which has preserved microbial diversity with marine affinity (Escalante et al. 2008, Alcaraz et al. 2008, Minckley and Jackson 2008, Souza et al. 2006), despite geological evidence indicating that marine waters retreated from the valley at the end of the proto-Gulf of Mexico formation 90 million years ago (Souza et al. 2008). Thus, the CCB possesses two key ingredients for hyperdiverse microbial endemism: isolation and long-term continuity. However, the role of animals or other vectors in microbial dispersion remains an intriguing hypothesis that warrants further investigation.
In this study, we found that several core genera are prevalent across sites, with Halanaerobium sp. and Desulfovermiculus sp. being particularly common. This suggests at least some degree of microbial dispersal between sites, possibly mediated by the deep aquifer, as previously proposed by Souza et al. (2006) in PR and Churince. Additionally, we report an important loss of microbial diversity in the CCB, as evidenced by the decline in the number of metagenome-assembled genomes (MAGs) at the AD site—from 125 MAGs in 2016 to only 22 in 2023. Our findings highlight the potential role of the deep aquifer as a microbial reservoir and underscore the urgent need for regulating water extraction to preserve what remains of the CCB’s unique biodiversity.
Methodology
Sampling
From 2020 to 2023, samples were collected from the AD site, in both the dome-shaped microbial mats (D) and orange circle-forming zones (C) (26° 49′ 41.7″ N 102° 01′ 28.7″ E), at depths ranging from the surface to ∼50 cm. A total of 20 samples were collected: 13 from the dome zone and 7 from the orange circle zone.
Additionally, in 2020, 2021, and 2023, six ponds within the PR system (26° 52′ 00.0″ N 102° 01′ 00.0″ W) were sampled (refer to Table S1 for detailed information on the distances between the sampled ponds within this system). The average distance between ponds was 0.1163 km, ranging from 0.026 km (between PR4 and PR5) to 0.026 km (between PR2 and PR6). For more information on the metadata associated with each sample refer to Table S2.
Sampling was conducted under collection permit SGPA/DGVS/03 188/20, issued by the Subsecretaría de Gestión para la Protección Ambiental, Dirección General de Vida Silvestre.
For Churince (26° 50′ 55.5″ N 102° 08′ 43.9″ W), we used data published by De Anda et al. (2018b), which consists of 12 metagenomes from water samples collected at the Churince site between 2012 and 2014.
Following collection, the samples were immediately transferred to liquid nitrogen and stored until total DNA extraction.
Total DNA extraction and amplicons sequencing
The DNA extraction followed a column-based protocol with a Fast DNA Spin Kit for Soil (MP Biomedical; Santa Ana, CA, USA). DNA was sequenced at the Laboratorio de Servicios Genómicos, LANGEBIO (http://langebio.cinvestav.mx/labsergen/) using the Illumina MiSeq 2 × 300 bp platform for the V4 region of the 16S rDNA gene. Universal primers for prokaryotes were employed for Polymerase Chain Reaction (PCR) amplification: the 515F primer 5′-GTGYCAGCMGCCGCGGTAA-3′ and the 806R primer 5′-GGACTACNVGGGTWTCTAAT′ (Walters et al. 2015).
Sequencing data preprocessing
The 16S rRNA gene sequence data underwent initial quality control procedures. FastQC (v0.11.8) (Andrews 2010) was employed to assess raw data quality. Subsequently, the cutadapt tool (v3.2) (Martin 2011), integrated directly within the R script, was used for trimming and filtering low-quality bases. The impact of this quality trimming on the distribution of read lengths was then evaluated to ensure the suitability of the data for further analysis. Following this, the required R libraries, including DADA2 (V.1.32.0) (Callahan et al. 2016), were loaded, and the trimmed and filtered sequencing data were imported into R as a DADA2-compatible object.
Denoising procedures were carried out using the DADA2 package. We used the filterAndTrim function for additional filtering and trimming based on quality. The learnErrors function was then used to estimate and model error rates in the data. Denoising was accomplished using the dada function, allowing for the inference of exact sequence variants (amplicon sequence variants, ASVs hereafter). This step aimed to improve the accuracy of the representation of biological sequences by addressing sequencing errors. Following denoising, chimeric sequences were identified and removed using the removeBimeraDenovo function. Subsequently, a nonchimeric sequence table was generated for taxonomic assignment and additional analyses.
Taxonomic assignments for the resulting ASVs were performed using the IdTaxa function from the DECIPHER package (v2.0) (Wright et al. 2012) in R. The SILVA_SSU_r138_2019 database was used as a reference. Nonchimeric sequences and their associated taxonomies were directly combined into a final table.
Alpha and beta diversity analysis of the 16S rRNA sequence data
For the 16S rRNA sequence data, phyloseq (v1.48.0) (McMurdie and Holmes 2013) in R was used to assess alpha diversity with Chao1 and Shannon indices, which provide information on both the richness and diversity of microbial communities, respectively. The analysis was conducted using the plot_richness function from the phyloseq package, which allows for the visualization of diversity indices based on the specified variables. In this analysis, the object ASV_phyloseq, which contains the ASV data, was used as input for the function. Beta diversity analysis was performed using the DESeq2 package (v1.44) in R (Love et al. 2014). Raw sequencing reads were transformed into a count table and metadata detailing sample grouping was incorporated. A DESeqDataSet object was created, specifying the experimental design, and subsequent normalization and variance-stabilizing transformations were applied to address variations in library sizes.
A dissimilarity matrix, representing pairwise distances between samples, was calculated based on the transformed count data using Euclidean distance metrics with the hclust function (Guénard and Legendre 2022) in R Hierarchical clustering, using Ward’s method to group samples according to their beta diversity profiles.
To further investigate beta diversity, we employed a two-step approach. First, count data were transformed into proportions using the transform_sample_counts function from the phyloseq package. This step ensures that subsequent analyses are not biased by differences in sequencing depth between samples, a common practice when working with compositional microbiome data.
Second, nonmetric multidimensional scaling (NMDS) was performed on the Bray–Curtis dissimilarity matrix, to explore the compositional dissimilarities between samples. The plot_ordination function, also from the phyloseq package, facilitated the visualization of the NMDS results.
To investigate the effects of depth (0–10 cm, 10–20 cm, 20–30 cm, 30–40 cm, and 40–50 cm), location (PR, D, and C), and sampling year on microbial community composition, a Bray–Curtis dissimilarity matrix was constructed using the phyloseq package. Subsequently, a permutational multivariate analysis of variance (PERMANOVA) was performed using the adonis2 function from the vegan package (v2.6–8). The analysis included 1 × 106 permutations to ensure robust results. To identify significant differences among groups for each variable pairwise comparisons were conducted using the pairwise.adonis function (v0.4) (Martinez-Arbizu 2020), with a Bonferroni adjustment for multiple testing.
Change in relative abundance at different detection thresholds between PR and AD
To visualize the core microbiome, we used the phyloseq package in R. First, the detection of each taxon was calculated, defined as the proportion of samples in which the relative abundance exceeded 1%. The prevalence of each taxon was then calculated, and a minimum threshold of 10% was applied to select the ASVs of interest. Then the data were normalized, by transforming them into proportions. Finally, a heatmap of the core microbiome was generated using the plot_core function, specifying prevalence thresholds (0%–100%) and detection thresholds (0.001, 0.003, and 0.009).
Comparison of microbial communities between Churince and AD using whole genome metagenomic sequencing
To further assess and compare diversity within and among sites in CCB, we analyzed shotgun metagenomic sequencing data from two sites: AD [including both the dome-shaped microbial mats (D) and orange circle-forming zones (C)] and Churince. For AD, we used 22 previously published metagenomes by Rodríguez-Cruz et al. (2024) along with three unpublished metagenomes obtained in 2023. Details on DNA extraction and sequencing can be found in Rodríguez-Cruz et al. (2024).
As mentioned above, for Churince we used data from 12 metagenomes derived from water samples collected between 2012 and 2014 (De Anda et al. 2018a). For detailed information on each sample refer to Table S3.
A heatmap of the core microbiome was generated using the plot_core function, from the microbiome package in R, specifying prevalence thresholds (0%–100%) and detection thresholds (0.001, 0.003, and 0.009).
The make_network function from the phyloseq package was used with a maximum distance value (max.dist) of 0.7 to define the connections between samples. The Bray–Curtis distance metric was used to quantify dissimilarity between samples, ranging from 0 (identical composition) to 1 (completely different composition), while accounting for the relative abundances of taxa.
MAGs from AD and Churince sites
We analyzed changes in the number of MAGs from 2016 to 2022 at the AD site (Rodríguez-Cruz et al. 2024), including unpublished MAGs from a 2023 sampling season, resulting in a total of 25 metagenomes analyzed.
The sampling period at AD consisted of four stages within the microbial mat-forming dome area (D). The first stage (2016–2019) included collections in April (coded as D1M01 in Table S3) and October 2016 (D1M02), representing dry and wet seasons, respectively. Additional annual samples were collected in February 2017, October 2018, and March and September 2019 (D1M01, D1M02, D1M03, D1M04, D1M05, and D1M06). The second stage, conducted in October 2020, involved samples collected at 10 cm intervals, covering depths from 0 to 50 cm, in both the dome area (D) and the orange-circle-forming area (C). This resulted in six samples (C1M08, C4M09, C5M10, D1M07, D4M11, and D5M12). The third stage (2021–2022) included 12 additional samples (D1M13, D2M14, D3M15, D4M16, D5M17, D1M18, D2M19, D3M20, D4M21, D5M22, D1M23, and C1M24), following the same sampling scheme as the second stage.
The fourth sampling period was in 2023, from which three metagenomes were obtained (D1M23, C1M24 y, and C1M25).
MAGs were assembled from the metagenomic data obtained in 2023, using the binning tools MaxBin2 (v2.2.7) (Wu et al. 2015), with the options run_MaxBin.pl -thread 24 -contig file_scaffolds.fasta -out directory_maxbin -abund file_scaffolds.bam.abundance.txt, Metabat2 (v2:2.15) (Kang et al. 2019) was executed with runMetaBat.sh -m 2500 file_scaffolds.fasta file_scaffolds.bam, and MetaCOAG (v1.0) (Mallawaarachchi and Lin 2022) was used with the options MetaCoAG –assembler spades –graph file_assembly_graph_with_scaffolds.gfa –contigs file_contigs.fasta –paths file_contigs.paths –abundance file_bam.abundance.txt –output directory_metacoag_out. The completeness and contamination of each MAG were assessed using the CheckM (v1.1.3) (Parks et al. 2015) software.
The same procedure was applied to the Churince data (De Anda et al. 2018a), which covers the period from 2012 to 2014 and includes 12 metagenomes. It is important to note that this dataset is the only one available for comparison with the shotgun sequencing data obtained from the AD site, as most of the data for the Churince system have been obtained from sequencing fragments of the 16S rRNA gene using different primers than those employed in this study.
Among the MAGs reconstructed from the AD and El Churince sites, two high-quality MAGs classified within the genus Ectothiorodospira were identified. One MAG originated from the AD dataset and the other from the El Churince dataset. Both MAGs met the MIMAG (Minimum Information about a MAG) criteria proposed by Bowers et al. (2017), which defines high-quality MAGs as having >90% completeness, <5% contamination, the presence of 23S, 16S, and 5S rRNA genes, and at least 18 tRNAs. To evaluate the nucleotide identity between these MAGs from the same genus, we used the software fastANI (Jain et al. 2018).
Given the number of MAGs analyzed in this study for the AD site (818 MAGs), we only performed functional annotation for the Churince MAGs. For this, we used freely available hidden Markov model (HMM) databases for microbial metabolic genes of environmental or biogeochemical relevance. For example, we used the metabolic-hmms database (available at https://github.com/banfieldlab/metabolic-hmms), FOAM (Functional Ontology Assignments for Metagenomes) (Prestat et al. 2014), TIGRFAMS (Haft 2003), and Pfam (Mistry et al. 2021) (V36.0). This annotation was then mapped to KEGG orthologs and normalized to the total number of coding sequences per genome.
A phylogenetic perspective of the MAGs obtained from the AD site (samples collected between 2016 and 2023) and Churince (samples collected between 2012 and 2014) was inferred using the tree generated by GTDB-tk (Chaumeil et al. 2020) (v2.3.1) and visualized with iTOL (Letunic and Bork 2021).
Furthermore, to determine whether the observed loss of diversity in terms of MAG number was due to variations in sequencing depth, we analyzed the average sequencing depth across years and sites using the Kruskal–Wallis test and Dunn’s post hoc test. Normality was assessed with the Shapiro–Wilk test. We also evaluated the correlation between the proportion of MAGs by taxonomical class and average sequencing depth using Pearson’s correlation. The same procedure was applied to the Churince data.
To assess microbial diversity across the three study sites (AD: D and C, PR, and Churince), we predicted 16S rRNA genes from metagenomic assemblies at Churince, due to the absence of direct 16S rRNA gene sequencing data. Genes were identified using Barrnap (v0.9) (Seemann 2013), classified taxonomically with Kaiju (Menzel et al. 2016), and visualized using ggplot2 in R.
Results
Alpha and beta diversity analysis for 16S rRNA data from AD and PR
We analyzed 16S rRNA tag data from two Eastern sites—PR, and AD in both the dome-shaped microbial mats (D) and orange circle-forming zones (C). In total, we processed 26 sequencing samples, each starting with a variable number of initial reads (refer to Tables S2, S4, and S5 for more information on each sample). Filtering and denoising significantly reduced the number of reads, with postfiltering counts ranging from 29 152 in sample D6104 (September 2021 at 0–10 cm depth) to 445 898 in sample D7110 (September 2021 at 10–20 cm depth). Denoised reads showed a gradual decrease compared to initial forward and reverse reads, highlighting the effectiveness of the cleaning steps.
Taxonomic analysis revealed important diversity at the phylum level. Halobacterota was the predominant phylum in multiple samples, peaking in sample C123 (March 2023 at 0–10 cm depth) at 0.698. Bacteroidota was also highly represented, particularly in D6104 (March 2022 at 30–40 cm depth) (0.495) and D8106 (March 2022 at 20–30 cm depth) (0.333). Nanoarchaeota was present in D6109 (September 2021 at 0–10 cm depth) (0.264) and PR6 (0.158). Halanaerobiaeota was prominent in PR5 with an abundance of 0.212, while Firmicutes, although less dominant, had a notable presence in D10113 (September 2021 at 40–50 cm depth) (0.152) and D8106 (March 2022 at 20–30) (0.140). Desulfobacterota also showed relevant abundance in D852 (October 2020 at 10–20 cm depth) (0.174) and D953 (October 2020 at 20–30 cm depth) (0.134).
The samples exhibited a wide range of species richness (see Tables S2, S4, and S5 for more information on each sample). Summary statistics, including read counts, and the number of unique ASVs are shown in Fig. 3A (see Fig. S1 for rarefaction curves).
For the dome-shaped microbial mats (D), observed ASVs ranged from 39 (sample D7105, September 2021, 0–10 cm depth) to 2668 (sample D7110, September 2021, 10–20 cm depth), highlighting a broad range of species richness in both Archaea and bacteria. Similarly the Shannon Index varied from 2.38 (D7105) to 6.54 (D7110), indicating a high variability in microbial diversity within this site.
For the orange circle-forming zones (C), observed ASVs ranged from 58 (sample C123, 2023, 0–10 cm depth) to 851 (sample C348, 2020, 20–30 cm depth). The Shannon index followed a similar trend, spanning from 3.52 (C123) to 5.60 (C348).
At the PR site, species richness was the highest among all sites, with observed ASVs ranging from 306 (sample PR3, 2020, 0–10 cm depth) to 2005 (sample PR6, 2020, 0–10 cm depth). The Shannon Index values for this site ranged from 4.69 (PR3) to 6.47 (PR2).
For the dome-shaped microbial mats (D site), archaeal richness ranged from 0 to 585 ASVs, with the lowest value in D7105 and the highest in D9112. The Shannon diversity index varied from 0.00 to 5.03, with the lowest in D7105 and the highest in D7110.
For the orange circle-forming zones (C), archaeal richness ranged from 32 to 224 ASVs, with the lowest value in C123 and the highest in C348. The Shannon diversity index varied from 2.88 to 4.49, with the lowest in C123 and the highest in C348.
For the PR site, archaeal richness ranged from 74 to 529 ASVs, with the lowest value in PR3 and the highest in PR2. The Shannon diversity index varied from 3.08 to 5.33, with the lowest in PR5 and the highest in PR2.
For the dome-shaped microbial mats (D site), bacterial richness ranged from 42 to 2407 ASVs, with the lowest value in D7105 and the highest in D7110. The Shannon diversity index varied from 2.49 to 6.40, with the lowest in D7105 and the highest in D7110.
For the orange circle-forming zones (C), bacterial richness ranged from 27 to 688 ASVs, with the lowest value in C123 and the highest in C348. The Shannon diversity index varied from 2.94 to 5.28, with the lowest in C123 and the highest in C348.
For the PR site, bacterial richness ranged from 246 to 1735 ASVs, with the lowest value in PR3 and the highest in PR6. The Shannon diversity index varied from 4.43 to 6.21, with the lowest in PR3 and the highest in PR2.
In the NMDS analysis (Fig. 2), samples PR1, PR2, PR4, and PR6 exhibit a distinct microbial profile, setting them apart from other groups. Their positive position on the MDS2 axis suggests that these microbial profiles are characterized by unique features not found in samples from orange circle-forming zones (C) and dome-shaped microbial mats (D) within the AD site.
NMDS ordination plot based on the Bray–Curtis dissimilarity matrix, illustrating compositional dissimilarities among microbial communities across samples. The data were obtained from 16S rRNA gene sequencing. Each point represents a unique microbial community, with distances between points reflecting the degree of dissimilarity. The ordination analysis yielded a final stress value of 0.1309, indicating a good two-dimensional representation of the data (stress <0.1 suggests a strong ordination with no real risk of false inferences, while stress <0.2 can still provide a usable visualization, though caution is needed for higher values) (Clarke 1993). Samples are color-coded by location to highlight spatial patterns in beta diversity, corresponding to sampling depths: A = 0–10 cm, B = 10–20 cm, C = 20–30 cm, D = 30–40 cm, and E = 40–50 cm.
The total dispersion within the PR group indicates substantial variability in microbiological composition. Samples PR3 and PR5 show relative proximity to some samples from group D (D953 and C449, taken at a depth of 30–40 cm) in the NMDS space (Fig. 2), implying potential similarities in their microbial profiles. However, the differences between these groups remain significant enough to clearly distinguish between them. This may indicate that while groups PR and D share some microbial characteristics, distinct environmental factors likely drive their divergence.
The results from the PERMANOVA analysis indicated that all considered variables were relevant concerning microbial composition. Sample depth had the most significant effect on microbial composition. (Pseudo-F = 7.189, R2 = 0.2348, P < .001), followed by location (Pseudo-F = 5.973, R2 = 0.1853, P < .001) and year (Pseudo-F = 3.210, R2 = 0.1021, P = .004). Post hoc pairwise comparisons for Depth indicated significant differences between groups, particularly between A (0–10 cm) and D (30–40 cm) (P < .001), and A (0–10 cm) and E (40–50 cm) (P = .009). For Location, significant differences were observed between C and D (P = .005), C and PR (P = .003), and D and PR (P = .001). Regarding Year, significant differences were noted between 2021 and 2023 (P = .004) and 2021 and 2022 (P = .001).
Core, shared, and unique ASVs obtained using 16S rRNA Sequence Data in AD and PR
The distribution and overlap of ASVs between different sampling sites are depicted in Fig. 3(A). We found a global core ASVs microbiota (i.e. shared taxa in the PR and AD sites), comprised by 237 ASVs (see Fig. S2 for a detailed distribution of the shared ASVs between PR and AD). Of these, 215 ASVs were classified at the order level. Among these, the most predominant orders were Halobacterales, with 27 ASVs, representing ~11.39% of the classified ASVs; followed by Halanaerobiales, with 16 ASVs (6.75%), and Woesearchaeales, with 15 ASVs (6.33%).
Core Microbiome Heatmap and Upset Plot Analysis from 16S rRNA gene sequencing samples. (A) Upset plot showing the shared and unique ASVs from the dome-shaped microbial mats (D) and orange circle-forming zones (C), both located within the AD site. “PR” corresponds to samples from the “Pozas Rojas” site. (B) Heatmap depicting the core microbiome from AD and PR sites, with prevalence thresholds set at 0%, 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 90%, and 100%, and detection thresholds at 0.001, 0.003, and 0.009.
Additional orders with important representation in the global core ASVs included Bacteroidales with 12 ASVs (5.06%), and Phycisphaerales and Spirochaetales, each with 10 ASVs (4.22%). Together, these five orders account for 33.97% of the classified ASVs. It is noteworthy that 22 ASVs were not classified at the order level, accounting for 9.28% of the total ASVs. Orders such as Actinomarinales, Alicyclobacillales, Anaerolineales, and Caulobacterales, among others, are represented by a single ASV each, contributing only 0.42% to the total classified ASVs of the global set. The shared ASVs among the three sites suggest the presence of a core microbial community, likely resulting from connectivity (dispersion) among sites, involving microbes capable of adapting and surviving in different physicochemical conditions.
The ASVs data also highlight the notable number of unique ASVs in PR, D, and C sites (1844, 1767, and 173 ASVs, respectively)
Core, shared, and unique microbial taxa between AD sites and Churince using whole genome metagenome data
Overall, alpha diversity values, measured using the Shannon index from whole genome metagenomic data collected at the AD site between 2016 and 2023 at depths of 0–50 cm, exhibited considerable variability. The dataset spans 8 years with samples collected at various depths: 0–10 cm (depth A), 10–20 cm (depth B), 20–30 cm (depth C), 30–40 cm (depth D), and 40–50 cm (depth E).
The Shannon diversity index values for both Archaea and Bacteria varied considerably across the sampled zones (refer to Table S5 for additional information on the Shannon index values for each sample). In the orange circle-forming zones (C), the Shannon indices for Archaea and Bacteria ranged from 5.56 to 6.15. In the dome-shaped microbial mats (D), these diversity values exhibited greater variability compared to the other zones, with indices ranging from a minimum of 3.90 for D1M13 (2021, Depth D) to a maximum of 7.19 for D2M14 (2021, Depth D).
In the Churince area, the Shannon indices for Archaea and Bacteria varied from 5.72 to 7.42, indicating a diverse microbial community.
Overall, the Shannon indices for both Archaea and Bacteria indicate that the dome-shaped microbial mats not only exhibited the greatest variability but also included both some of the highest and lowest among all zones.
The Shannon diversity index values for Archaea varied significantly across different zones. In the orange circle-forming zones (C), the indices remained consistently high, ranging from 4.10 to 4.79. In contrast, the dome-shaped microbial mats (D) exhibited the greatest variation among samples, with a minimum of 4.05 for D1M01 (2016, Depth D) and a maximum of 4.85 for D5M12 (2020, Depth D). In the Churince area (CH), the Shannon indices for Archaea ranged from 3.39 to 4.44, reflecting a diverse microbial community.
The Shannon diversity index values for Bacteria varied significantly across the sampled zones. In the orange circle-forming zones (C), the indices ranged from 6.36 to 7.09. In the dome-shaped microbial mats (D), Bacterial diversity exhibited greater variability than in the other zones, with indices spanning from a 5.16 for D1M13 (2021, Depth D) to of 7.79 for D2M14 (2021, Depth D).
In the Churince area (CH), the Shannon indices for Bacteria ranged from 6.48 to 8.06. Sample S2 displayed the lowest diversity at 6.48, while sample S10 had the highest diversity at 8.06.
Overall, the Shannon indices indicate that the dome-shaped microbial mats (D) not only exhibited the greatest variability but also included some of the highest and lowest values among all zones.
The prevalence analysis of taxa in the AD and Churince samples, using a minimum detection threshold of 0.1% relative abundance, revealed a set of taxa consistently detected across a large fraction of the samples (Fig. 4A). This suggests that these organisms may be capable disperse among the different ponds in the CCB.
Heatmap and network analysis of the microbiome from the AD site and Churince, based on whole metagenome sequencing data. At the AD site, (D) represents dome-shaped microbial mats, and (C) represents "orange circle-forming zones". Samples from Churince are labeled as CH (S1 to S12). (A) Heatmap of the core microbiome with prevalence thresholds set at 0%, 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80%, 90%, and 100%, and detection thresholds at 0.001, 0.003, and 0.009. (B) Network analysis illustrating relationships among samples based on relative abundance data. Networks were constructed using the make_network function in phyloseq with a maximum distance (max.dist) of 0.7. Bray–Curtis dissimilarity was used to quantify sample composition (ranging from 0 for identical to 1 for fully distinct). Sampling depths are color-coded as follows: A = 0–10 cm, B = 10–20 cm, C = 20–30 cm, D = 30–40 cm, and E =40–50 cm.
Shotgun metagenomic sequencing data indicate that ~21.2% of shared microbial taxa belong to Archaea, including Halodesulfurarchaeum formicicum, Halorubrum, Halorhodospira halophila, Haloterrigena spp., Halovivax asiaticus, Salinarchaeum sp. Harcht-Bsk1, and Thermoplasmata archaeon. In contrast, Bacteria account for 78.8% of the shared microbiome, encompassing diverse groups such as Desulfobacteraceae (e.g. Desulfosalsimonas propionicica, Desulfatitalea tepidiphila, Desulfonatronovibrio magnus, Desulfonema ishimotonii, Desulfohalobium retbaense), Deltaproteobacteria (e.g. Deltaproteobacteria bacterium), and Spirochaetaceae (e.g. Spirochaetaceae bacterium). Additional bacterial taxa include Actinobacteria, Acidobacteria, Gemmatimonadetes, Halanaerobium, Halomicrobium, Halomonas, Halorussus, Mariniphaga, Nitrospirae, Phycisphaeraceae, Planctomycetes, Puniceicoccaceae, Rhodobacteraceae, Rhodosalinus, and Candidatus Melainabacteria.
Taxa with high prevalence, such as Halanaerobium sp. (63%), Desulfovermiculus halophilus (60%), and Spirochaetaceae bacterium, 57%) (see Fig. 4A), may indicate their ability to disperse among different aquatic systems. Their presence suggests a significant role in microbial communities, contributing to environmental stability and facilitating communication between communities sharing the same deep aquifer within the basin.
To explore the relationships between samples based on taxon composition, a network analysis was performed using metagenomic data transformed to relative abundance. Notably, we identified a connection between sample D4M11 (AD depth 30–40 cm), and sample S8 (Churince), as well as between metagenomes D5M17 (AD depth 40–50 cm) and S8 (Churince) (Fig. 4B). These findings suggest a possible link—via the deep aquifer—between these two distant sites in the CCB, despite the samplings being taken a decade apart. This connection becomes more evident at greater depths, as samples from 30 to 40 cm at AD may be more strongly influenced by the deep aquifer, facilitating the dispersal of microorganisms also found in Churince.
MAGs from 2016 to 2023 show a decrease of microbial diversity in AD site
From the 25 metagenomes generated for AD site between 2016 and 2023, we assembled a total of 818 MAGs (171 from the domain Archaea and 647 from bacteria; see supplementary file 1 for more information on each MAG) with completeness values >20% (see Fig. 5A and B for a phylogenetic perspective). High-quality MAGs must meet the MIMAG criteria (Bowers et al. 2017): completeness above 90%, contamination below 5%, presence of the 23S, 16S, and 5S rRNA genes, and at least 18 tRNAs. In our study, bacterial MAGs were classified into 70 taxonomic classes for (see Fig. S3) and used to reconstruct a phylogeny (Fig. 5A and B). Similarly, for the archaeal MAGs, were assigned to 13 taxonomic classes (see Fig. S3).
Phylogenetic analysis of MAGs from the AD and Churince sites (2016–2023). The phylogenetic tree was constructed using GTDB-tk v2.3.1 and visualized in iTOL. Panel A shows bacterial MAGs, while panel B shows archaeal MAGs.
We observed a declining trend in the number of MAGs recovered as the samples became more recent (see Fig. 6 and Fig. S3B). Specifically, the oldest sample, D1M01, had 125 MAGs, while the most recent sample, D1M23, contained only 22 MAGs. This decrease is evident in several intermediate samples: D1M02 (34 MAGs, 2016), D1M03 (46 MAGs, 2017), D1M04 (47 MAGs, 2018), D1M05 (91 MAGs, 2019), D1M06 (92 MAGs, 2019), D1M07 (89 MAGs, 2020), D1M13 (39 MAGs, 2021), D1M18 (8 MAGs, 2022), and D1M23 (22 MAGs, 2023), as shown in Fig. 6 and Fig. S3B. Examining these data series, we note that after sample D1M01, which has the highest number of MAGs, there is an important drop in sample D1M02, which has only 34 MAGs. This trend remains relatively constant, with minor fluctuations until reaching sample D1M18, which has the lowest number of MAGs at just 8. Although there is a slight increase in sample D1M23 (22 MAGs) compared to D1M18, it is still significantly lower than in the older samples.
Heatmap comparing taxonomic classes present or absent at the AD sight from 2016 to 2023. A taxonomic class is considered present if at least one MAG from that class was detected in the first year of sampling (2016). Rows represent taxonomic classes, and columns represent the sampling years.
A Kruskal–Wallis test found no significant differences in sequencing depth across years (Kruskal–Wallis chi-squared = 18.006, df = 18, P = .4552. Normality was assessed with the Shapiro–Wilk test, W = 0.95134, P = .2688). However, Dunn’s post hoc test revealed a significant difference between 2022 and 2023 (P = .018), while no other years showed significant differences. Additionally, Pearson correlation coefficient indicated no significant linear relationship between the proportion of MAGs by taxonomical class, and average sequencing depth (r = −0.112, P > .05).
We attempted to assemble MAGs from 12 metagenomes collected from Churince between 2012 and 2014 (see Fig. S5 for relative abundance estimates based on predicted 16S rRNA genes). However, we recovered only 5 MAGs with more than 20% completeness, and less than 10% contamination, all belonging to the Bacteria domain (listed in Table 1). Notably, MAGs S2_1 and S5_1, both classified within the family Ectothiorhodospiraceae, meet the high-quality MAG criteria established by Bowers et al. (2017) under the MIMAG (Minimum Information about a MAG) guidelines.
MAGs from the Churince site. This table lists MAGs with completeness values >20% and contamination values under 10%.
| Phylum | Class | Order | Family | Completeness | Contamination | |
|---|---|---|---|---|---|---|
| S2_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 96.44 | 0 |
| S5_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 96.15 | 0.85 |
| S2_2 | Bacteroidota | Rhodothermia | Balneolales | Cyclonatronaceae | 63.8 | 1.09 |
| S9_1 | Cyanobacteriota | Cyanobacteriia | Thermostichales | Thermostichaceae | 55.26 | 0.88 |
| S8_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 36.21 | 0 |
| Phylum | Class | Order | Family | Completeness | Contamination | |
|---|---|---|---|---|---|---|
| S2_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 96.44 | 0 |
| S5_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 96.15 | 0.85 |
| S2_2 | Bacteroidota | Rhodothermia | Balneolales | Cyclonatronaceae | 63.8 | 1.09 |
| S9_1 | Cyanobacteriota | Cyanobacteriia | Thermostichales | Thermostichaceae | 55.26 | 0.88 |
| S8_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 36.21 | 0 |
MAGs from the Churince site. This table lists MAGs with completeness values >20% and contamination values under 10%.
| Phylum | Class | Order | Family | Completeness | Contamination | |
|---|---|---|---|---|---|---|
| S2_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 96.44 | 0 |
| S5_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 96.15 | 0.85 |
| S2_2 | Bacteroidota | Rhodothermia | Balneolales | Cyclonatronaceae | 63.8 | 1.09 |
| S9_1 | Cyanobacteriota | Cyanobacteriia | Thermostichales | Thermostichaceae | 55.26 | 0.88 |
| S8_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 36.21 | 0 |
| Phylum | Class | Order | Family | Completeness | Contamination | |
|---|---|---|---|---|---|---|
| S2_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 96.44 | 0 |
| S5_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 96.15 | 0.85 |
| S2_2 | Bacteroidota | Rhodothermia | Balneolales | Cyclonatronaceae | 63.8 | 1.09 |
| S9_1 | Cyanobacteriota | Cyanobacteriia | Thermostichales | Thermostichaceae | 55.26 | 0.88 |
| S8_1 | Pseudomonadota | Gammaproteobacteria | Ectothiorhodospirales | Ectothiorhodospiraceae | 36.21 | 0 |
Discussion
Alpha and beta diversity analysis for 16S rRNA data from AD and PR
Taxonomic analysis of microbial diversity at sites AD and PR revealed a wide range of phyla. Notably, the dominance of Halobacterota in several samples, such as C123 (0.698), aligns with this phylum’s known prevalence in hypersaline environments (Oren 2003). Members of the archaeal class Halobacteria represent the most successful microbial group in hypersaline conditions and thus have been used as key model organisms in astrobiological experiments (Stan-Lotter et al. 2002, Wu et al. 2022). Their resilience and ability to withstand multiple environmental extremes (polyextremophiles), including high salinity, radiation, vacuum, extreme cold, and growth in the presence of caustic materials like perchlorate, highlight their relevance in astrobiological contexts (Stan-Lotter et al. 2002, Wu et al. 2022). The importance of CCB as a key site for astrobiological research, as previously reported (Souza et al. 2012, 2020), is underscored by these findings.
The detection of Nanoarchaeota in samples D6109 (0.264) and PR6 (0.158) is noteworthy given this phylum’s typical association with extreme environments (St. John and Reysenbach 2019). As previously reported for site AD, Nanoarchaea constitute ~4.27% of the total MAGs (35 out of 818) microbial community (Rodriguez-Cruz et al. 2024). The metagenomes and MAGs reported in this study provide a unique opportunity to explore the role of Nanoarchaea in microbial communities, potentially shedding light on their interactions with diverse archaeal and bacterial hosts (Medina-Chávez and Travisano 2022), particularly given their limited metabolic capabilities (Huber and Kreuter 2014, St. John and Reysenbach 2019).
Alpha diversity analysis for both bacteria and archaea between samples at different depths revealed substantial variability in microbial structure. Surface samples from PR, such as PR1, PR2, and PR6, exhibited high Shannon diversity indices (5.376, 5.793, and 5.802, respectively), suggesting local factors shape surface microbial diversity. In contrast, intermediate and deep samples at AD, such as D7110 (10–20 cm) and D9112 (30–40 cm), displayed higher diversity (5.84 in both samples), possibly due to the more stable environments that promote diverse microbial communities (Lozupone and Knight 2007, Thompson et al. 2017).
The NMDS analysis indicated that samples from PR exhibited distinct microbial profiles compared to samples from AD, both the dome-shaped microbial mats (D) and orange circle-forming zones (C). This separation in NMDS space suggests significant differences in taxonomic composition, likely driven by variations in the physicochemical conditions of these environments. The distinct profiles highlight variations in both the structure and function of microbial communities across these diverse habitats (Stahl et al. 2013).
According to the PERMANOVA test, sampling depth emerged as the most important variable, explaining 23.48% of the variation in microbial composition, indicating that changes in environmental conditions associated with different depths can have a major impact on microbial communities. Location also exhibited a significant effect, explaining 18.53% of the variation. Sampling year was significant, explaining 10.21% of the variation, but its impact is lower compared to depth and location.
Core, shared, and unique ASVs obtained using 16S rRNA Sequence Data in AD and PR
Halanaerobium ASV12 (see Fig. 3B) exhibited a notably high prevalence, being detected in 88.46% of the samples, suggesting its widespread distribution and potential ecological significance in these environments. Members of this genus have been found in high abundance in subterranean systems characterized by high salinity, such as the Barnett Shale geological formation located in the Bend Arch-Fort Worth basin (USA), where their relative abundances range from 17% to 33% (Liang et al. 2016). Halanaerobium spp. are recognized for their roles in carbohydrate fermentation and sulfur production via thiosulfate reduction (Ravot et al. 1997, 2005). Therefore, this bacterium likely plays a crucial role in organic matter biodegradation and sulfur production in the CCB.
Another highly prevalent ASV among CCB ponds is Desulfovermiculus ASV48 (see Fig. 3B), detected in 65.38% of the samples. Desulfovermiculus spp. (Belyakova et al. 2006) has been identified as a halophilic sulfate reducer, suggesting that sulfate reduction may contribute to sulfide production in high-salinity subterranean systems. The high prevalence of this genus across ponds indicates its potential role as a key component of the shared microbial community.
Notably, both Halanaerobium sp. and Desulfovermiculus sp. have been found in high-salinity environments and at shallow depths, such as the Barnett Shale geological formation (Liang et al. 2016). This suggests these taxa may be transported across different ponds via a common deep aquifer. However, further sampling at greater depths in CCB is needed to better understand the microbial community dynamics in the region.
ASVs such as Spirochaeta_2 ASV8 and Halocella ASV30, with prevalences of 69.23%, also showed a widespread distribution among the samples, reinforcing the hypothesis of a set of common and predominant taxa in these sites. Desulfosalsimonas ASV46, has been described as a halophilic Gram-negative sulfate-reducing bacterium (Kjeldsen et al. 2010). Consequently, it is a common member of extreme hypersaline sediments, where salinities range from 120 to 270 g of NaCl per liter (Kjeldsen et al. 2010). In turn, Desulfonatronovibrio ASV14 has been described as an alkaliphilic, sulfate-reducing bacterium (Zhilina et al. 1997) and has been isolated from a soda-depositing lake (Zhilina et al. 1997). Desulfonatronovibrio comprises motile vibrio bacteria that utilize only hydrogen and formate as electron donors, while using sulfate, sulfite, and thiosulfate—but not sulfur—as electron acceptors (Zhilina et al. 1997). Desulfonatronovibrio ASV14 was detected in 69.23% of the samples, suggesting its potential role in specific biogeochemical processes, particularly in the sulfur cycle within the ponds analyzed in this study.
Other ASVs exhibited lower prevalence, such as Kocuria ASV28, which was detected in only 19.23% of the samples, suggesting that they may occupy more restricted ecological niches or that their populations are more vulnerable to environmental fluctuations. Alternatively, it could indicate a limited dispersal capacity, hindering their ability to spread across the different ponds within CCB.
Core, shared, and unique microbial taxa between AD sites and Churince using whole genome metagenome data
Samples from El Churince exhibited the highest Shannon diversity indices across all sites when considering both Archaea and Bacteria. For example, S10 (depth A: 0–10 cm) had a combined index of 7.42, with Bacteria contributing 8.06 and Archaea contributing 4.21. In all samples from El Churince, the bacterial Shannon diversity index consistently exceeded the archaeal index. For instance, in S1 (depth A), the bacterial index was 7.93, compared to 3.90 for Archaea.
Samples from orange circle-forming zones displayed intermediate Shannon diversity index values when both Archaea and Bacteria were considered. The indices ranged from 5.56 in C1M08 (depth A) to 6.15 in C5M10 (depth E). A consistent pattern emerged when analyzing the domains separately: the bacterial Shannon diversity index consistently surpassed the archaeal index. For instance, in C1M24 (depth A), the bacterial index was 7.07, compared to 4.79 for Archaea, highlighting bacterial predominance in this environment.
Samples from dome-shaped microbial mats (D) exhibited a range of Shannon diversity index values for Archaea and Bacteria, showing lower diversity compared to other sites. For instance, D1M13 (depth A) had the lowest combined index (3.90), while D2M14 (depth B) exhibited a higher index (7.19). When analyzing each domain separately, the bacterial Shannon diversity index consistently exceeded the archaeal index across all samples. For example, in D5M12 (depth E), the bacterial index reached 7.00, while the archaeal index was 4.85, underscoring the dominant role of bacteria in these mats.
Among the samples, D4M11 (depth D) stood out with a total diversity index of 6.83, potentially due to hydrological connectivity with the deep aquifer, which may facilitate microbial dispersal and mixing.
Previous isotopic and other studies have shown that groundwater from the deep aquifer is the primary water source for the different CCB aquatic systems (Felstead et al. 2015, Johannesson et al. 2004, Wolaver et al. 2013). We propose that this deep aquifer has preserved bacterial lineages by maintaining conditions similar to those of the ancient ocean (Espinosa-Asuar et al. 2020, Gomez-Lunar et al. 2018, Souza et al. 2018). These stable conditions have kept ancient microbial lineages isolated from their marine relatives for millions of years (Moreno-Letelier et al. 2011, 2012), allowing sufficient time for the emergence of high microbial diversity.
To date, no evidence supports microbial transport between ponds in the CCB via by animals such as insects or birds (Souza et al. 2006). Instead, current understanding suggests that shared environmental factors, particularly the influence of a deep aquifer, explain microbial connectivity. This aquifer may also serve as a common source of microbial diversity through dispersion within the CCB, preserving microbial lineages with marine affinities (Alcaraz et al. 2008, Escalante et al. 2008, Minckley and Jackson 2008, Souza et al. 2006).
Additionally, extensive phylogenetic evidence indicates that multiple microorganisms from the CCB, such as members of the genus Halothece sp., cluster monophyletically. This suggests that these organisms are more closely related to each other than to external entities, implying shared unique characteristics absent in their distant ancestors (Rodriguez-Cruz et al. 2024). A similar monophyletic clustering pattern has been previously reported in culturable microorganisms from the phylum Actinobacteria in the CCB (Arocha-Garza et al. 2017).
Nevertheless, we believe that additional geophysical studies, along with further analyses incorporating new metagenomic data from deeper samples across different sites in CCB and nearby areas, will be necessary. However, in some systems, this may not be possible, as the Churince area and nearby sites in the west wing of the CCB have now dried up. Given this limitation, we explicitly compared the number of MAGs lost over just 8 years of studying the AD site.
MAGs from 2016 to 2023 indicate a decline of microbial diversity at the AD site
Our data suggest a progressive loss of microbiological diversity over time (Fig. 6). Fluctuations observed in intermediate samples (D1M05, D1M06, and D1M07) may reflect specific events that temporarily impacted microbiological diversity. However, the overall trend points to a decline. This loss of diversity in the more recent AD samples may be linked to changes in dispersal, possibly influenced by the connectivity between the deep aquifer and surface layers, driven by environmental or anthropogenic factors.
While global climatic change may have influenced climate patterns and water availability in the CCB, evidence suggests that anthropogenic disturbances are key drivers of microbial diversity loss throughout the CCB. For instance, between 2000 and 2009, natural vegetation cover declined by 604.91 hectares, primarily due to agricultural expansion and plowing activities (González Hernández et al. 2003, FMCN–CONANP 2009).
Additionally, the Cuatro Ciénegas Protected Natural Area (84 347.47 hectares), has lost 43.02 hm³ of water volume between 2000 and 2024, as documented by the Mexican federal government (Secretaría del Medio Ambiente y Recursos Naturales and Instituto Mexicano de Tecnología del Agua 2024) and discussed by Leal Nares et al. (2022).
If these trends persist, microbial diversity could continue to decline, potentially leading to the complete disappearance of the aquatic system, similar to what occurred with the Churince system, which dried up entirely in 2011. Furthermore, the PR system remains at risk.
A total of 5 MAGs were recovered from Churince (Table 1), three of which belong to the genus Ectothiorhodospira. Interestingly, four MAGs from the same genus were also assembled from AD. When comparing the average nucleotide identity between two MAGs of Ectothiorhodospira from each site (AD and Churince; see Fig. S4A), focusing on those with the highest completeness and lowest contamination values (93.9% completeness and 1.42% contamination for MAG C1M0812 from the AD site), we observed an average nucleotide identity of 90.49%. According to the commonly accepted threshold of 95%–96% similarity to consider two genomes as belonging to the same species (Richter and Rosselló-Móra 2009), we find that MAGs C1M08_12 from orange circle at AD and S2_1 from Churince may represent different species within the Ectothiorhodospira genus. However, since we are not comparing complete genomes, as C1M08_12 has 93.9% completeness and S2_1 a completeness of 96.44%, we cannot entirely rule out the possibility that these two MAGs belong to the same species.
Furthermore, when examining the metabolic capabilities of MAGs belonging to the genus Ectothiorodospira, C1M08_12 and S2_1. (see Fig. S4B), we found no notable differences between their metabolisms (see Fig. S4B for details): both C1M08_12 and S2_1 have genes associated with assimilatory sulfate reduction, including PAPSS (3′-phosphoadenosine 5′-phosphosulfate synthase) (Peng and Verma 1995), sat (sulfate adenylyltransferase) (Carlson et al. 2021), CysNC (bifunctional enzyme CysN/CysC), CysH (phosphoadenosine phosphosulfate reductase), and cysJ (sulfite reductase [NADPH] flavoprotein alpha-component) (Hummerjohann et al. 1998). Additionally, we identified the soxB, soxC, and soxY genes as part of the sulfur oxidation (Sox) multienzyme complex (Ghosh et al. 2009), which is crucial in the oxidation of reduced sulfur compounds, such as thiosulfate to sulfate. This complex is believed to play a key role in sulfur oxidation, an essential process for energy extraction and sulfur cycling in sulfur-rich environments, as previously reported in the CCB (De Anda et al. 2018b).
We also detected genes in the C1M08_12 and S2_1 MAGs from the genus Ectothiorodospira, involved in dissimilatory sulfate reduction, such as aprAB (adenylyl sulfate reductase) and dsrA (dissimilatory sulfite reductase alpha subunit) (Simon and Kroneck 2013). In addition to sulfur metabolism, we detected in the C1M08_12 and S2_1 MAGs, genes linked to nitrogen metabolism, including nitrogen fixation genes nifD and nifK, encoding nitrogenase molybdenum–iron proteins involved in the conversion of atmospheric nitrogen into ammonia (Fani et al. 2000).
Genes like nirS and norB (Black et al. 2016) were also found in the C1M08_12 and S2_1 MAGs, involved in denitrification, by reducing nitrite to nitric oxide and subsequently to nitrogen gas.
Furthermore, the detection of nosZ (Black et al. 2016) in the C1M08_12 and S2_1 MAGs, which is responsible for reducing nitrous oxide to nitrogen, and amoA (Li et al. 2015), involved in ammonia oxidation during nitrification, suggests a comprehensive capacity for nitrogen cycling within these MAGs.
Conclusions
In this new study, we combined results from three aquatic systems (AD, PR, and Churince) and identified common microorganisms across the different sites. For instance, Desulfovibrio spp. was found in the three aquatic systems. Additionally, in the comparison of MAGs that persisted over the 8 years of sampling, Desulfovibrio was the only taxonomic class of MAGs that remained prevalent throughout the entire period. These sulfate-reducing microorganisms have been reported in granitic groundwater sampled taken from depths of 50–600 m, where they are frequently found in deep granitic rock aquifers (Pedersen 1997, 2000). This finding may support our hypothesis that microorganisms are transported from the deep aquifer to various surface sites in CCB. Moreover, the detection of Ectothiorhodospira sp. at both AD and Churince sites could suggest microbial transport from the deep aquifer to different ponds within CCB. Alternatively, it may indicate that the physicochemical conditions at these sites are suitable for sustaining these microbes.
Acknowledgments
Ulises Erick Rodriguez Cruz is a doctoral student from the the Programa de Doctorado en Ciencias Biomédicas, Universidad Nacional Autónoma de México (UNAM) and received SECIHTI fellowship 857544. Manuel Ochoa-Sánchez is a doctoral student from Posgrado en Ciencias Biológicas, Universidad Nacional Autonóma de México (UNAM) and acknowledges the SECIHTI fellowship 917392. This paper was written during the sabbatical leaves of VS and LEE 2024–2025, and they wish to thank the UNAM for allowing them to have this sabbatical. We extend our gratitude to Dra Rosalinda Tapia-Lopez, Dra Erika Aguirre-Planter, Silvia Barrientos, and Manuel Rosas from the Instituto de Ecología, Universidad Nacional Autónoma de México, for their technical and field assistance. We also thank PRONATURA Noreste for granting access to the Pozas Azules ranch, and to Dra Eria Rebollar-Caudillo and Dr Arturo Becerra-Bracho for their valuable feedback on the manuscript. Additionally, we are grateful to Rodrigo Garcia-Herrera, for facilitating the use of the high-performance computing cluster “Patung” located at the Laboratorio Nacional de Ciencias de la Sostenibilidad (LANCIS), Instituto de Ecología (UNAM).
Author contributions
Ulises E. Rodriguez-Cruz (Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Writing – original draft, Writing – review & editing), Manuel Ochoa-Sánchez (Conceptualization, Methodology, Visualization, Writing – original draft, Writing – review & editing), Luis E. Eguiarte (Conceptualization, Data curation, Funding acquisition, Investigation, Project administration, Resources, Software, Supervision, Validation, Writing – review & editing), and Valeria Souza (Funding acquisition, Investigation, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing)
Conflict of interest
All authors declare that they have no conflicts of interest.
Funding
SECIHTI supported U.E.R.-C. with the doctoral scholarship CVU: 857544 and M.O.-S. with the doctoral scholarship CVU: 917392. This research was funded by PAPIIT-DGAPA, UNAM IG200319 grants (awarded to V.S. and L.E.E.), PAPIIT-DGAPA, UNAM IN204822 (awarded to V.S.), and the operating budget of the Instituto de Ecología, Universidad Nacional Autónoma de México provided to V.S. and L.E.E.
