Community ecology and functional potential of bacteria, archaea, eukarya and viruses in Guerrero Negro microbial mat

1 Vol.:(0123456789) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports Community ecology and functional potential of bacteria, archaea, eukarya and viruses in Guerrero Negro microbial mat P. Maza‑Márquez 1,2*, M. D. Lee 1,3 & B. M. Bebout 1 In this study, the microbial ecology, potential environmental adaptive mechanisms, and the potential evolutionary interlinking of genes between bacterial, archaeal and viral lineages in Guerrero Negro (GN) microbial mat were investigated using metagenomic sequencing across a vertical transect at millimeter scale. The community composition based on unique genes comprised bacteria (98.01%), archaea (1.81%), eukarya (0.07%) and viruses (0.11%). A gene‑focused analysis of bacteria archaea, eukarya and viruses showed a vertical partition of the community. The greatest coverages of genes of bacteria and eukarya were detected in first layers, while the highest coverages of genes of archaea and viruses were found in deeper layers. Many genes potentially related to adaptation to the local environment were detected, such as UV radiation, multidrug resistance, oxidative stress, heavy metals, salinity and desiccation. Those genes were found in bacterial, archaeal and viral lineages with 6477, 44, and 1 genes, respectively. The evolutionary histories of those genes were studied using phylogenetic analysis, showing an interlinking between domains in GN mat. Microbial mats are one of the most ancient ecosystems known, having persisted through around 85% of the Earth’s history1 and played a key role in the evolution of Earth’s atmosphere2,3. Today’s mats are modern analogues of these first ecosystems on the Earth. As microbial mats were likely the locations at which oxygen was first produced in an otherwise anoxic world, they may offer an ecological model to understand both the evolution of biochemical cycles and microbial adaptation to drastic environmental changes, including the advent of an oxygenated atmosphere. Microbes found in microbial mats have been shown to exhibit a number of adaptive responses to extreme environmental conditions. In one of the few metagenomics studies of microbial mats, genes involved in adaptive responses and resilience against high-UV irradiation, elevated salinity conditions, oxidative stress and heavy metal resistance were described from microbial mats from Shark Bay4– one of the most extensive marine microbial mat systems in the world. Guerrero Negro microbial mat is one of the best studied microbial mat ecosystems, and the use of molecular analysis based on sequencing of small-subunit rRNA genes and 454 sequencing have revealed the high complexity of this stratified ecosystem5,6. Mats located in salterns managed for the production of salt are permanently submerged by a hypersaline water column which serves to protect the vertical structure of the mats—relative to mats found in intertidal environments subjected to more profound environmental disturbance7. Earlier work detected a vertical stratification of oxygen and sulfide as well as a surprising occurrence of sulfate reduction in the oxic zone of the mat8. A recent study revealed a vertical patterns of nitrogen cycling genes with respect to vertical variations in oxygen concentration9. The vertical organization of other functions in these communities has been less well studied. GN microbial mats are characterized by extreme vertical chemical gradients at micrometer to millimeter spatial scales, and these vertical gradients have been shown to correspond to changes in microbial community composition6. Bacteria and archaea have been well documented in GN, as well as fungi in a recent study (within the eukarya domain)10. However, the composition and variability of viruses and their functions at fine scale resolutions in these communities is less known. Through metagenomics, the current study provides a taxonomic description of the vertical taxonomic organization as well as a functional organization delineated between bacteria, archaea, eukarya and viruses in a GN microbial mat—revealing new insights into the ecology of these communities. The goals of the present study were to characterize the community structure and their functional potential in this microbial mat through: (I) analyzing the functional annotations and taxonomic OPEN 1Exobiology Branch, NASA Ames Research Center, Moffett Field, CA, USA. 2University of Granada, Granada, Spain. 3Blue Marble Space Institute of Science, Seattle, WA, USA. *email: [email protected] 2 Vol:.(1234567890) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ classification of assembled genes spanning bacterial, archaea, eukarya and viruses; (II) examining potential genetic mechanisms of adaptation in microbial mat; and (III) exploring the interlinkages between genes present in bacteria, archaea, eukarya, and viruses through gene-level phylogenetic analysis. Material and methods Microbial mat sampling The sampled microbial mats are located in salterns managed by the world’s largest salt-producing company (Exportadora de Sal SA ESSA), situated on the Pacific Ocean side of the Baja peninsula. Exportadora de Sal has minimum human activity and consists of a series of 13 concentration areas with an approximate extension of 28,184 hectares. The concentration areas were constructed approximately 2m above sea level in the low lands adjacent to the ponds “Ojo de Liebre” and “Guerrero Negro” (Fig.S1). The microbial mats were collected in June 2019 from hypersaline ponds concentration area 4 (Fig.S1), as previously described10. At the time of mat collection, salinity of brine was 125 ppt, temperature 24.4 C, pH 8.3, ammonium concentration 0.12μM, dissolved oxygen 7mg/L, and nitrate concentration was below the limit of detection (< 0.5μM). Samples from the microbial mats were collected with a stainless-steel corer of 1cm diameter as previously described10. Three replicate cores were placed into sterile centrifuge tubes (Falcon®, Corning, Corning, NY, USA), capped, and immediately frozen in liquid nitrogen. To get the vertical layers at one-millimeter intervals of the first four layers (0–1, 1–2, 2–3 and 3–4mm from the top of the mat), the mat was sectioned using sterile scalpels. Three replicates were pooled for metagenomic analysis, resulting in a single pooled metagenome for each depth for library preparation and sequencing. DNA extraction Total DNA extraction (from approximately 0.20g per sample) was performed from each microbial mat layer, using a DNeasy Power Biofilm Kit, (Qiagen, Venlo, The Netherlands), according to manufacturer’s instructions. A nanophotometer (Implen GmbH, München, Germany) was used for checking the quality (A260/A280) and quantity (A260) of extracted genomic DNA. Library preparation and metagenomic sequencing were performed at Molecular Research (MR DNA, Texas, USA,http:// www. mrdna. org/ conta c t. html). Libraries were prepared using the Nextera DNA Flex library kit (Illumina) following the manufacturer’s instructions, and sequencing was performed on the NovaSeq 6000 platform (2 × 150 nucleotides). Metagenomic data processing Metagenomics processing with annotated code is documented at our open-Science Foundation site, https:// osf. io/ 9kwn3/ wiki 11. Conda (2020; www. anaco nda. com (accessed on 10 March 2021)) was utilized for program installation and environment management. Read quality was scanned with FastQC v0.11.912 and reads trimmed/filtered with trimmomatic v0.3913. A co-assembly of all 4 depths was performed with SPAdes v3.14.014, the assembly was filtered and summarized with bit v1.8.1615 (see TableS1 for assembly summary statistics), and each individual samples’ reads were mapped to the filtered co-assembly with bowtie2 v2.3.5.116 and sorted and indexed with samtools v1.917. Metagenomic sequence data from the 4 depths are available through NCBI’s Sequence Read Archive at BioProject PRJNA688760. NCBI accession numbers: SRX9761389, SRX9761388, SRX9761387 and SRX9761386. Taxonomic classification and functional annotation Co-assembly and read-mapping files were integrated into anvi’o v6.218 for annotation with the KEGG database19 and parsing and extraction of the gene-level coverage and detection data (with “detection” being the proportion of a gene that recruited any reads to it). Gene-level taxonomic classification was performed with CAT v5.1.220 against the NCBI nr database. Normalization and analyses were performed with R v3.6.3 (Core Team 2017) in Rstudio v1.1.456 (www. rstud io. com). To mitigate non-specific read-recruitment, gene-level coverage information was filtered based on detection (proportion of a gene that recruited any reads to it), such that those with a detection less than 50% had their coverage set to 0. This had a net effect of removing less than 3% of the prefiltered total coverage. Information on KEGG’s Metabolism pathway was accessed with KEGGREST v1.26.0, defining the KEGG Orthology (KO) terms. The gene table coverages were normalized across the 4 samples by dividing each value by its sample’s total coverage and multiplying by 1 million to generate values of Coverage per Million (herein referred to as CPM). Entrez Direct v13.9, www. ncbi. nlm. nih. gov/ books/ NBK17 9288/) was utilized to search and retrieve reference sequences from NCBI. Heatmap and cluster analyses based on Bray Curtis similarities were calculated for all genes with a mean coverage of >= 9 in when summed across the 4 samples. Tree construction For making phylogenetic trees, sequences were aligned with Muscle v3.8.155121, and then trimmed using trimal v1.4.1. Trees were generated using the open-source software IQ-TREE v2.2.022. We inferred the maximum-likelihood tree with auto-model selection via the built-in ModelFinder (option `-m MFP`) using 1000 bootstraps23. The trees were edited through the Interactive Tree of Life web-interface24. Accurately rooting is essential for the correct interpretation of the genetic changes between sequences since the root gives the directionality of evolution within the tree25. However, due to the span of diversity we were considering (across all domains), and to have a consistent method of rooting across all the generated gene trees, we utilized mid-point rooting26. All trees in Newick format are archived with https:// doi. org/ 10. 6084/ m9. figsh are. 25018 007 and available here: https:// doi. org/https:// doi. org/ 10. 6084/ m9. figsh are. 25018 007. 3 Vol.:(0123456789) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ Viral analysis For viral specific identification, contigs from the initial co-assembly were screened for viral sequences using VIRSorter2 v3, then checkV was passed for quality control of the VirSorter2 results, following previous protocols27. Viral read of each sample were mapped using bowtie2 v2.3.5.116 and sorted and indexed with samtools v1.917. Anvi’o v6.218 was used to annotate viral gene and calculate coverage profiles. For virus gene level taxonomic classification, two databases were used: RefSeq8428 and IMGVR29. To avoid false positive viral gene identification, the genes were filtered based on two criteria: genes taxonomic classified in the RefSeq84 and/or IMGVR databases; and functional viral protein annotated in KEGG database. Results A gene‑focused view of bacteria, archaea, eukarya, and viruses. As described in Methods, metagenomes from all 4 depths were co-assembled together, genes were identified, and normalized coverage values were attained by recruiting the individual sample reads to the assembled contigs. Mean coverage values for genes were extracted, and for exploratory purposes these mean coverage values for each sample were normalized to be out of 1 million (coverage per million (CPM)). Genes identified were taxonomically classified and functionally annotated, and here we break down those results. We also include information on read-based classification, as a means to potentially identify if any large biases might have been introduced through the assembly-to-gene-calling process, but the results did not largely vary in any cases. A gene‑focused view of Bacteria Taxonomically classified within the bacterial domain, 764,307 unique genes were assembled and identified, with 461,332, 679,559, 656,264 and 677,165 of them having a coverage greater than 0 in Layer 1, Layer 2, Layer 3, and Layer 4, respectively (TablesS2 and S3; Fig.1A1). However, the total normalized coverages of the bacterial genes decreased with depth (Fig.1A1). Layer 1 contained the greatest coverage of bacterial genes, with 901,325.3 CPM, while Layers 3 and Layer 4 contained the lowest, 864,618.79 and 861,536.72 CPM, respectively. Heatmap and cluster analyses based on Bray Curtis resemblance of the genes are shown in Fig.S2A1. The global similarity was > 40%, with Layer 1 being the least similar and clustering out separately; the similarity between Layer 2, Layer 3, and Layer 4 was > 80%. SIMPROF analysis based on Bray Curtis resemblance at 5% of significance level detected a significant difference between Layers 3 and Layer 4 (p = 0.001) and Layer 3-Layer 4, and Layer 2 (p = 0.001), and between Layer 3-Layer 4Layer 2 and Layer 1 (p = 0.001) 32.8% of bacterial genes were successfully functionally annotated with KO terms, comprising 357,984 ± 13,375 CPM across the 4 layers (mean ± 1SD; TableS2). According to KEGG’s groupings, genes with >= 9 average CPM across the 4 layers were ascribed to 42 metabolic pathways, with the most abundant being related to genetic information processing, signaling and cellular processes, carbohydrate metabolism, and energy metabolism (TableS2). Gene-level taxonomic classification within the Bacteria detected 44 phyla, 106 families, and 430 species (Fig.1A2, S3, TableS3). The dominant phyla were Cyanobacteria, followed by Chloroflexi, Proteobacteria, Firmicutes and Bacteroidetes. Cyanobacteria were dominant in Layer 1, while Chloroflexi was dominant in Layer 4. At family level, Microcoleaceae was the most prevalent family, followed by Coleofasciculaceae and Oscillatori‑ aceae. At species level, the dominant species were Chloroflexi bacterium, Anaerolineae bacterium, Anaerolineales bacterium, Geitlerinema sp. PCC 9228 and Cyanobacteria bacterium J055. Chloroflexi bacterium, Geitlerinema sp. PCC 9228, Cyanobacteria bacterium J055 and Coleofasciculus chthonoplastes were the predominant species in Layer 1, while Chloroflexi bacterium, Anaerolineae bacterium, Anaerolineales bacterium and Thermoflexia bacterium were dominant species in Layer 4. A gene‑focused view of Archaea For the archaeal domain, 14,148 unique genes were identified with 2,762, 8,338, 13,250, 13,793 having coverages > 0 for Layer 1, Layer 2, Layer 3, and Layer 4, respectively. The highest coverages of genes were detected in Layer 3 and 4 (Fig.1B1), which were > 7 times more abundant than in Layer 1. The dendrogram in Fig.S2A2 shows the heatmap and cluster analyses of samples based on Bray Curtis similarities of the archaeal composition of the community. Layers 3 and 4 were more similar overall, at > 70% similarity, while Layer 1 was the least similar (Fig.S2A2). According to SIMPROF analysis, there was a significant difference in the archaeal gene coverages (based on an alpha value of 0.05) between Layer 3 and Layer 4 (p = 0.022), Layers 3Layer 4 and Layer 2 (p = 0.001), and Layers 3Layer 4Layer 2 and Layer 1 (p = 0.001). 25.7% of archaeal genes were successfully functionally annotated with KO terms, with those unannotated accounting for 7,792.4 ± 5,215.8 CPM across the depths (TableS4). The 14,148 genes were classified into 26 KEGG metabolic pathways (TableS4). Genetic information processing was the most abundant pathway, while the lowest abundant pathways were related to amino acids, cellular processes, cell motility, and sulfur metabolisms (TableS4). Gene-level taxonomic classification of archaea is summarized in Figs.1B2, S4 and TableS5. A total of 20 phyla, 36 families and 496 species were identified. The dominant phyla were Euryarchaeota, followed by Candidatus Lokiarchaeota, Candidatus Micrarchaeota, Candidatus Woesearchaeota, Candidatus Thorarchaeota. At family level, Methanosarcinaceae were the dominant family, followed by Methanotrichaceae, Methanobacteriaceae and Candidatus Methanoperedenaceae. At species level, Candidatus Lokiarchaeota archaeon were the dominant species, followed by Thermoplasmata archaeon, Candidatus Micrarchaeota archaeon and Candidatus Woesearchaeota archaeon. A gene‑focused view of Eukarya For the eukaryal domain, 547 genes were detected, with 403, 396, 298, 318 of them having coverages > 0 for Layer 1, Layer 2, Layer 3 and Layer 4, respectively. The greatest coverages were identified in Layer 1, with 1,696.08 CPM 4 Vol:.(1234567890) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ (Fig.1C1). The coverages in Layers 2–4 were > 6 times lower than in Layer1, with 241.6, 188.6, 257.7 CPM for Layer 2, 3 and 4, respectively. The heatmap and cluster analysis representing the similarity of the 547 genes with depth is show in Fig.S2A3. The global similarity of the eukaryal community was low, < 5%. Layer 1 represented Figure1. Summaries of normalized gene-level coverages broken down by domain; note the x-axes vary between A1-D1. (A1) Bar plots of the genes identified in bacteria across the uppers 4 layers examined [Layer 1 (0-1mm from surface), Layer 2(1–2mm from surface), Layer3 (2–3mm from surface), and Layer 4 (3-4mm from surface)]. (A2) Heatmap showing the read-based taxonomic classification of bacteria at species level with coverages >= 9 across the uppers 4 layers examined. (B1) Bar plots of the CPM of the 14,184 genes identified in archaea across the uppers 4 layers examined [Layer 1 (0–1mm from surface), Layer 2(1–2mm from surface), Layer 3 (2–3mm from surface), and Layer 4 (3–4mm from surface)]. (B2) Heatmap showing the read-based taxonomic classification of archaea at species level across the uppers 4 layers examined. (C1) Bar plots of the CPM of the 547 genes identified in eucaryote across the uppers 4 layers examined [Layer 1 (0–1mm from surface), Layer 2(1–2mm from surface), Layer 3 (2–3mm from surface), and Layer 4 (3–4mm from surface)]. (C2) Heatmap showing the read-based taxonomic classification of eucaryotes at species level across the uppers 4 layers examined. (D1) Bar plots of the CPM of the 394 genes identified in viruses across the uppers 4 layers examined [Layer 1 (0–1mm from surface), Layer 2(1–2mm from surface), Layer 3 (2–3mm from surface), and Layer 4 (3–4mm from surface)]. (D2) Heatmap showing the read-based taxonomic classification of viruses at species level across the uppers 4 layers examined. 5 Vol.:(0123456789) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ the least similar one, while the similarity between Layers 3 and 4 was > 80%. According to SIMPROF analysis, there was a significant difference in the coverage (based on an alpha value of 0.05) between Layer 1 and Layer 2-Layer 3-Layer 4 (p = 0.001), and between Layer 2-Layer 3 and Layer 4 (p = 0.001), but not between Layer 3 and Layer 4 (p = 1), where the similarity was > 80%. 26.7% of eukaryal genes were successfully functionally annotated with KO terms, with those unannotated accounting for 377.6 ± 352.6 CPM across the depths (TableS6). According to KEGG orthology classification, the 547 genes were grouped in 11 metabolic pathways (TableS6). The most abundant pathways were genetic information processing and photosynthesis. Overwise, the minor coverages of genes were related to the pathways sulfur relay system and quorum sensing (TableS6). Gene-level taxonomic classification of eukarya identified 19 phyla, 119 families and 15 species (Fig.1C2, S5, TableS7). Bacillariophyta were the dominant phyla, followed by Streptophyta and Ascomycota. At family level, Thalassiosiraceae were the most predominant family, followed by Symbiodiniaceae and Bacillariaceae. At species level, Thalassiosira pseudonana was the predominant species, followed by Symbiodinium microadriaticum and Fistulifera solaris. A gene‑focused view of Viruses 835 genes were assembled and taxonomically classified as viral, with 451, 474, 791 and 759 of them having a coverage greater than 0 in Layers 1, 2, 3, and 4, respectively. The coverages of viral genes increased with depth, with Layers 3 and 4 having coverage values > 3 times higher than in Layer 1 (Fig.1D1). Figure S2A4 shows the heatmap and cluster analysis based on the coverages of the 835 genes identified. The global similarity between the layers was < 20%, with Layer 1 again being the least similar. A SIMPROF test at 5% of significance level revealed no significance difference in the coverage of the genes between Layer 3 and Layer 4 (p = 0.35), but detected significant differences between Layers 3-Layer 4 and Layer 2 (p = 0.001), and between Layer 2-Layer 3Layer 4 and Layer 1 (p = 0.001). 17.9% of identified viral genes were successfully functionally annotated with KO terms, with those unannotated accounting for 1531.0 ± 1095.9 CPM across the depths (TableS8). KEGG grouped these genes into 10 pathways (TableS8), the highest coverages were related to genetic information processing and phage terminase large subunit/ phage replication initiation protein. FigureS6 shows the read-based classification of the identified virome at phyla and family levels: 3 phyla, 16 families and 149 species were identified (TableS9). The dominant phyla were Uroviricota, followed by Nucleocytoviricota and Hofneiviricota. At species level, the most abundant species was related to uncultured Caudoviri‑ cetes phage, following by Prokaryotic dsDNA virus sp., uncultured Mediterranean phage uvMED, uncultured Mediterranean phage, Streptomyces phage Brock, Cellulophage phage phi17:1, Podoviridae sp. cty5g4, Marseillevi‑ rus LCMAC101 (Fig.1D2). Living in Guerrero Negro microbial mat: bacteria, archaea, and virus genes related to poten‑ tial adaptation mechanisms Here we focus on recovered genes related to antibiotic and multidrug resistance, heavy metal toxicity, oxidative damage genes, cold, heat and phage shock proteins, UV-radiation stress genes, salinity and desiccation stress conditions. We recovered 6477 unique genes annotated with those functions that were classified within the bacterial domain and 44 within the archaeal domain (Fig.2, TablesS10, S11, S12, S13, S14, S15, S16 and S17). Moreover, one gene (related to UV-DNA damage endonuclease TableS18) was recovered that was classified as originating from a virus. According to KEGG, those genes were classified within membrane transport (76), signaling and cellular processes (488), antimicrobial resistance genes (245), signal transduction (394), genetic information processing (4220), metabolism (758), transport and catabolism (24), replication and repair (10) and energy metabolism (301) (Fig.2). The presence of these genes may have implications for adaptability and resilience to stress conditions, as has been previously described in another microbial mat4. Genes related to potential adaptation mechanism in Bacteria For Bacteria, 6477 genes related to potential environmental adaptation were detected (Fig.2A and supplementary TablesS10–S16). According to KEGG orthology classification, the 6477 genes were classified within membrane transport (76), signaling and cellular processes (480), antimicrobial resistance genes (245), signal transduction (393), genetic information processing (3834), metabolism (749), transport and catabolism (24), replication and repair (370) and energy metabolism (301): 564 unique genes were related to multidrug resistance protein/multidrug efflux pump (TableS10); 86 unique genes were detected for multiple antibiotic resistance protein; 280 genes were detected for fluoroquinolone, catechol, tetracycline, fosmidomycin, quaternary compound, tetracycline and vancomycin resistance protein; 92 unique genes were identified to resistance to arsenical, copper, mercuric, tellurite and zinc (TableS11); 25 genes were related to heavy metal (TableS12); 410 involved in dealing with oxidative stress (TableS13); 237 genes were related to cold shock proteins (TableS14); 283 genes were identified as heat shock proteins (TableS14); 281 genes for phage shock protein (TableS14); 997 genes were related to UV damage: excinuclease ABC subunit ABC, DNA helicase II /ATP-dependent DNA helicase PcrA, UV DNA damage endonuclease, and ATP-dependent DNA helicase UvsW (TableS15).; and 3222 genes were associated with desiccation and salinity conditions (TableS16): 1329 genes were annotated as RNA polymerase sigma70 factor, ECF subfamily (rpoE), 301genes for F-type H + /Na + -transporting ATPase subunit beta (ATPF1B, atpD), 281 genes for chaperonin GroEL (groEL, HSPD1), 251 genes for chaperonin GroES (groES, HSPE1), 259 genes for molecular chaperone GrpE (GRPE), 121 genes for molecular chaperone DnaJ (dnaJ), 116 genes for molecular chaperone DnaK (dnaK, HSPA9), 15 genes for trehalose 6-phosphate synthase (otsA), 35 genes for trehalose 6-phosphate phosphatase (otsB), 81 genes for osmoprotectant transport system substrate-binding protein (opuC), 113 genes for glycine betaine/proline transport system ATP-binding protein (proV), 106 genes for glycine betaine/proline transport system permease protein (proW), 206 genes for glycine betaine/proline 6 Vol:.(1234567890) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ transport system substrate-binding protein (proX) and 8 genes for L-ectoine synthase (ectC). Summary of these genes and their normalized coverages are showed in Fig.2A. Overall, layer 1 contained the highest coverages of genes involved in antibiotic and multidrug resistance (pumps ATP-binding cassette, subfamily B, multidrug efflux pump and membrane fusion protein, multidrug efflux system), resistance to metals (arsenic, copper resistance protein, tellurite, tellurium and tetracycline genes), heavy metal, oxidative stress genes (superoxide dismutase, Cu–Zn family, superoxide dismutase Fe, Mn family and superoxide oxidase) and phage and heat shock proteins. On the other hand, deeper layers contained the greatest coverages of genes related to multidrug resistance genes (MFS transporter, ACDE family, multidrug resistance protein, multidrug resistance protein, MATE family and outer membrane protein, multidrug efflux system), zinc resistance gene, superoxide reductase gene, cold shock protein and genes associated with UV-resistance/repair (uvrA, uvrB, uvrC and uvrD genes). Genes involved with desiccation and salinity stress conditions had similar Figure2. (A) Heatmap showing the coverages of potentially adaptation-relevant genes present in bacteria across the uppers 4 layers examined. (B) Heatmap showing the coverages of potentially adaptation-relevant genes present in archaea across the uppers 4 layers examined. (C) Heatmap showing the UV gene present in virus across the uppers 4 layers examined. 7 Vol.:(0123456789) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ coverages across layers, detecting the greatest coverages for genes to encoded the production of exopolysaccharides (rpoE) and molecular chaperones (GroES, GroEL, DnaJ and DnaK). Deeper discussion to explain their relevance in the context of potential ecological adaptations of microbial mats is detailed in Section“Discussion”. Genes related to potential adaptation mechanism in Archaea In the archaeal domain, 44 genes identified in the metagenome data were related to stress conditions: 2 for multidrug resistance, 4 for multiple antibiotic resistance, 8 for interactions with oxidative genes, 16 shock proteins and 14 for UV genes (Fig.2B, TableS17). Layer 3 and 4 contained the greatest coverages of multidrug resistance genes, antibiotic, oxidative genes (superoxide reductase and superoxide dismutase), cold and phage shock proteins and UV genes. On the other hand, antibiotic resistance genes were not recovered in Layer 1 or 2; while heat shock protein hsIJ was higher than htpX and the greatest coverages were found in Layer 2. According to KEGG orthology classification, the 44 genes were classified within signaling and cellular processes (8), signal transduction (1), genetic information processing (18), metabolism (9), replication and repair (8). Deeper discussion to explain their relevance in the context of potential ecological adaptations of microbial mats is detailed in Section“Discussion”. Interlinking between the genes present in bacteria, archaea and viruses using phylogenetic analysis To further examine the distribution of genes related to potential adaptative mechanisms in terms of their evolutionary relatedness, we built phylogenetic trees. A total of 12 trees were built for the genes: MFS transporter, ACDE family, multidrug resistance protein (KO ID K08221), small multidrug resistance pump (KO ID K03297), superoxide dismutase, Fe–Mn family (KO ID K04564), superoxide reductase (KO ID K05919), heat shock protein (KO ID K03799), cold shock protein (KO ID K03704), phage shock protein (KO ID K03973), excinuclease ABC subunit A (KO ID K03701), excinuclease ABC subunit B (KO ID K03702), excinuclease ABC subunit C (KO ID K03703), DNA helicase II ATP dependent DNA helicase (KO ID K03657) and UV DNA damage endonuclease (KO ID K13281). These genes were present in bacteria, archaea and viruses domains as follow: for KO ID K08221: 14 sequences were present in bacteria and 1 in archaea; for KO ID K03297: 21 sequences were detected in bacteria and 1 in archaea; for KO ID K04564: 171 sequences were present in bacteria and 1 in archaea; for KO ID K05919: 107 sequences were found in bacteria and 7 in archaea; for KO ID K03799: 106 sequences were present in bacteria and 4 in archaea; for KO ID K03704: 221 sequences were found in bacteria and 7 in archaea; for KO ID K03973: 114 sequences were present in bacteria and 2 in archaea; for KO ID K03701: 362 sequences were present in bacteria and 5 in archaea; for KO ID K03702: 193 sequences were present in bacteria and 1 in archaea; for KO ID K03703: 203 sequences were found in bacteria and 3 in archaea; for KO ID K03657: 214 were identified in bacteria and 1 in archaea; for KO ID K13281: 17 sequences were detected in bacteria, 2 in archaea and 1 in viruses. Figures3, 4, 5 and 6 show the phylogenetic relationships between domains; the complete trees are in the supplementary material Figs.S7–S15. And see Table1 and the corresponding supplemental tables for gene IDs and full amino-acid sequences. Figure3 shows the phylogenetic tree of MFS transporter, ACDE family, multidrug resistance protein (yitG, ymfD, yfmO; 3A) and small multidrug resistance pump (emrE, qac, mmr, smr; 3B). For ACDE family, multidrug resistance protein (yitG, ymfD, yfmO) (Fig.3A), the gene ID 57669 present in archaea belonging to Candidatus Methanolliviera hydrocarbonicum was between clades holding genes classified as sourced from the following bacterial lineages: Spirochaetes bacterium (gene ID 391435), Candidate division Zixibacteria bacterium (gene ID 13873), Bacteroidetes bacterium (genes IDs 219193, 511031), Tangfeifania diversioriginum (gene ID 127482) and the genes IDs 158770, 316573, 83331 and 254423. With respect small multidrug resistance pump (emrE, qac, mmr, smr) (Fig.3B), the gene ID 419991 present in Archaeoglobales archaeon formed a clade with 4 genes present in bacteria: the genes IDs 424272, 737623, gene ID 683720 present in Desulfonemais himotonii and gene ID 212872 present in Olavius sp. associated proteobacterium Delta1. Figure4 and Fig.S7–S8 show the phylogenetic trees for superoxide dismutase, Fe–Mn family (SOD2) (Fig.4A and S7) and superoxide reductase (dfx) (Fig.4B and S8). For SOD2, the gene ID 550693 present in archaea formed a clade with a gene ID 201256 present in bacteria. For dfx gene, the gene ID 612886 present in archaea formed a clade with 3 genes present in bacteria lineages: two present in Deltaproteobacteria (genes IDs 56860, 26015) and one present in bacteria gene ID 612887. The gene ID 571069 present in Candidatus Micrarchaeota was closely related to genes present in Deltaproteobacteria (gene ID 59597) and Desulfobacteraceae (gene ID 241743). Finally, 5 genes present in archaea with genes IDs 323416, 494295, 151464, 732097, 733106 presents in Thermoplasmatales and Thermoplasmata formed a clade. This clade was closely related to a clade formed by the bacterial lineages Candidatus Cloacimonas sp. (411781), Peptoclostridium litorale (685247), Clostridia (847882), Planctomycetes bacterium (315242) and the gene IDs 40055, 525801, 693963, 411780, 6209. Figure5 and FiguresS9–S11 show the phylogenetic trees corresponding to heat, cold and phage shock protein. For heat shock protein HtpX (htpx) (Fig.5A and S9), the gene IDs 443088 and 222562 from the archaea lineages of Thermoplasmata and Candidatus Micrarchaeota formed a clade with a gene recovered from a Deltaproteobac‑ teria (gene ID 17526). The genes IDs 350980 and 60713 present in archaea deeply branched in between clades holding bacterial genes from Planctomycetes (563474, 523502, 102322), Phycisphaerales (754540), candidate division Zixibacteria (73164) and Desulfofustis sp (446773). Regarding to cold shock protein (cspA) (Fig.5B and S10), the genes ID 806860, 245087, 649029, 199423 and 224539 detected from the archaeal lineages of Thermo‑ plasmata, Candidatus Aenigmarchaeota archaeon and Thermoplasmatales formed a clade with two genes present in Chloroflexi bacterium (genes IDs 641093, 589131, 623971, 155751). The genes ID 454595 and 191917 presents in Candidatus Woesearchaeota and Euryarchaeota were in a deeply branching clade with genes from Planctomycetota (gene ID 76796) and Bacteroidetes (genes IDs 20000, 628969). For phage shock protein C (pcpc) (Fig.5C and S11), the gene ID 259327 present in archaea formed a clade with genes from Bacteroidetes bacterium (gene 8 Vol:.(1234567890) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ ID 585054). The gene ID 254336 present in Methanomas siliicoccalesa formed a clade with bacterial genes from Alphaproteobacteria bacterium (gene ID 125844) and Bacteroidetes bacterium (gene ID 585054). The phylogenetic trees representing the UV genes (uvrA, uvrB, uvrC, uvrD, pcrA, uvsE, UVE1) are shown in Fig.6 and in the supplementary Figs.S12–S15. For uvrA gene (Fig.6A and S12), three genes classified as coming from the archaeal lineages of Euryarchaeota and Thermoplasmata (genes IDs 848457, 361714, 636321) formed a clade with a bacteria gene ID 519469 from Candidatus Buchananbacteria bacterium. One gene classified as coming from Methanothermobacter (gene ID 727168) was within a deep clade with bacterial genes from Phycisphaerae (gene ID 701638) and Brachyspira (gene ID 577760). For uvrB gene (Fig.S13), the gene ID 653100 taxonomically classified as Thermoplasmata branched near genes classified as Planctomycetota (genes IDs 728796, 54766). For uvrC gene (Fig.6B and S14), the genes IDs 699307, 905297, 647629 corresponding to archaea formed a clade with 6 genes present in bacterial lineages: Pseudomonadota (gene ID 873573), Spirochaetes (gene ID 664616), Spirochaetia (gene ID 391558), Spirochaetia (gene ID 399873), Spirochaetia (gene ID 179402) and Gemmatimonadetes bacterium (gene ID 82391). For uvrD, pcrA gene (Fig.S15), the gene ID 181004 present in Candidatus Bathyarchaeota formed a clade with a gene ID 297131 present in Bacteroidales bacterium. For uvsE, UVE1 gene (Fig.6C), two genes present in Thermoplasmata (genes IDs 807563 and 510363) and one gene present in uncultured Caudoviricetes phage (gene ID 667218) formed a clade with 3 genes present in bacteria (genes IDs 572443, 553365 and 74783, representing a Chloroflexi bacterium, an unclassified bacterium, and a Deltaproteobacteria bacterium, respectively). Discussion Guerrero Negro microbial mat is one of the best studied microbial mat ecosystems; however, the vertical functional organization has been less well studied. In this study, 922,765 unique gene-copies were recovered (meaning assembled and predicted), with 84.51% of those being classified to at least the domain level, leaving 15.49% unclassified (TableS19). The greatest coverages of bacteria and eukarya genes were detected in Layer 1, while the highest coverages of viruses and archaea genes were found in Layers 3 and 4 (Fig.1). The upper one-millimeters 254423 Tangfeifaniadiversioriginum(127482) Bacteroidetesbacterium (511031) 83331 316573 158770 Bacteroidetesbacterium (219193) CandidatusMethanolliviera hydrocarbonicum(57669) Spirochaetesbacterium(391435) candidatedivision Zixibacteriabacterium(13873) Gemmatimonadetesbacterium(161363) 6839 625332 130037 240999 Treescale: 0.1 135285 Halieaceae bacterium (692358) 333319 Xanthomonadales bacterium (165873) 15297 360116 654474 498441 Rhodovulum sp.12E13 (803754) Roseicyclus(115830) Rhodobacteraceae bacterium HLUCCO18 (404976) 789305 209865 392516 572799 26137 Archaeoglobales archaeon (419991) 424272 Desulfonemaisishimotonii (683720) 737623 Olaviussp.associated proteobacteriumDelta1(212872) Roseospira visakhapatnamensis (440177) Treescale:0.1 A. MFS transporter, ACDE family, multidrug resistance protein (yitG, ymfD, yfmO) B. Small multidrug resistance pump (emrE, qac, mmr, smr) Bacteroidota Candidate division Zixibacteria Gemmatimonadota Pseudomonadota Spirochaetota A A r r c c h h a a e e a a p p h h y y l l a a Bacteria phyla Euryarchaeota 75 14 47 57 55 54 62 97 47 17 100 33 83 99 32 26 29 90 68 35 98 72 59 82 45 100 26 63 77 70 87 Figure3. (A) Mid-point rooted phylogenetic tree based on amino-acid sequences detected in our study from MFS transporter, ACDE family, multidrug resistance protein (yitG, ymfD, yfmO). (B) Phylogenetic tree based on amino-acid sequences from small multidrug resistance pump (emrE, qac, mmr, smr). Amino-acid sequences from bacteria (black color), amino acid sequences from archaea (purple color). For built the tree, 37 sequences were included in the phylogenetic analysis (15 sequences for yitG, ymfD, yfmO genes and 22 sequences for emrE, qac, mmr, smr genes). The sequences were aligned with Muscle and the tree was generated with IQTREE2 with 1000 bootstraps, with auto-model selection via the built-in ModelFinder. 9 Vol.:(0123456789) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ Figure4. (A) Mid-point rooted phylogenetic tree based on amino-acid sequences detected in our study from superoxide dismutase, Fe–Mn family (SOD2). (B) Phylogenetic tree based on amino-acid sequences from superoxide reductase (dfx). Amino-acid sequences from bacteria (black color), amino acid sequences from archaea (purple color). For built the tree, 286 sequences were included in the phylogenetic analysis (172 sequences for SOD2 gene and 114 sequences for dfx gene). The sequences were aligned with Muscle and the tree was generated with IQTREE2 with 1000 bootstraps, with auto-model selection via the built-in ModelFinder. 16 Vol:.(1234567890) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ 56. Kurth, D. et al. Arsenic metabolism in high altitude modern stromatolites revealed by metagenomic analysis. Sci. Rep. 7, 017–00896. https:// doi. org/ 10. 1038/ s4159801700896-0 (2017). 57. Tam, H. K., Wong, C. M. V. L., Yong, S. T., Blamey, J. & Gonzalez, M. Multiple-antibiotic-resistant bacteria from the maritime antarctic. Polar Biol. 38, 1129–1141. https:// doi. org/ 10. 1007/ s003000151671-6 (2015). 58. Wei, S., Higgins, C., Adriaenssens, E., Cowan, D. & Pointing, S. Genetic signatures indicate widespread antibiotic resistance and phage infection in microbial communities of the McMurdo Dry Valleys, East Antarctica. Polar Biol. 38, 919–925. https:// doi. org/ 10. 1007/ s003000151649-4 (2015). 59. Zaikova, E. et al. Antarctic relic microbial mat community revealed by metagenomics and metatranscriptomics. Front. Ecol. Evol. https:// doi. org/ 10. 3389/ fevo. 2019. 00001 (2019). 60. Breitbart, M. et al. Metagenomic and stable isotopic analyses of modern freshwater microbialites in Cuatro Ciénegas. Mexico. Environ. Microbiol. 11, 16–34. https:// doi. org/ 10. 1111/j. 14622920. 2008. 01725.x (2009). 61. Giaquinto, L. et al. Structure and function of cold shock proteins in archaea. J. Bacteriol. 189(15), 5738–5748. https:// doi. org/ 10. 1128/ JB. 0039507 (2007). 62. Zhang, B. et al. Conserved TRAM domain functions as an archaeal cold shock protein via RNA chaperone activity. Front. Microbiol. 8, 1597. https:// doi. org/ 10. 3389/ fmicb. 2017. 01597 (2017). 63. Chanda, P. K., Mondal, R., Sau, K. & Sau, S. Antibiotics, arsenate and H2O2 induce the promoter of Staphylococcus aureus cspC gene more strongly than cold. J. Basic Microbiol. 49(2), 205–211. https:// doi. org/ 10. 1002/ jobm. 20080 0065 (2009). 64. Matz, J. M., Blake, M. J., Tatelman, H. M., Lavoi, K. P. & Holbrook, N. J. Characterization and regulation of cold-induced heat shock protein expression in mouse brown adipose tissue. Am. J. Physiol. 269(1), R38–R47. https:// doi. org/ 10. 1152/ ajpre gu. 1995. 269.1. R38 (1995). 65. Cao, Y. et al. TGF-beta1 mediates 70-kDa heat shock protein induction due to ultraviolet irradiation in human skin fibroblasts. Pflugers Arch. 438(3), 239–244. https:// doi. org/ 10. 1007/ s0042 40050 905 (1999). 66. Luo, Z. H. et al. Diversity and genomic characterization of a novel parvarchaeota family in acid mine drainage sediments. Front. Microbiol. https:// doi. org/ 10. 3389/ fmicb. 2020. 612257 (2020). 67. Lebre, P., De Maayer, P. & Cowan, D. Xerotolerant bacteria: Surviving through a dry spell. Nat. Rev. Microbiol. 15, 285–296. https:// doi. org/ 10. 1038/ nrmic ro. 2017. 16 (2017). 68. Assaha, D. V. M., Ueda, A., Saneoka, H., Al-Yahyai, R. & Yaish, M. W. The role of Na+ and K+ transporters in salt stress adaptation in glycophytes. Front. Physiol. 8, 509. https:// doi. org/ 10. 3389/ fphys. 2017. 00509 (2017). 69. Reno, M. L., Held, N. L., Fields, C. J., Burke, P. V. & Whitaker, R. J. Biogeography of the Sulfolobus islandicus pan-genome. Proc. Natl. Acad. Sci. U. S. A. 106(21), 8605–8610. https:// doi. org/ 10. 1073/ pnas. 08089 45106 (2009). 70. Shu, W. S. & Huang, L. N. Microbial diversity in extreme environments. Nat. Rev. Microbiol. 20, 219–235. https:// doi. org/ 10. 1038/ s4157902100648-y (2002). 71. Decho, A. W., Norman, R. S. & Visscher, P. T. Quorum sensing in natural environments: Emerging views from microbial mats. Trends Microbiol. 18, 73–80. https:// doi. org/ 10. 1016/j. tim. 2009. 12. 008 (2010). 72. Gallagher, K. L., Kading, T. J., Braissant, O., Dupraz, C. & Visscher, P. T. Inside the alkalinity engine: The role of electron donors in the organomineralization potential of sulfate-reducing bacteria. Geobiology. 10, 518–530. https:// doi. org/ 10. 1111/j. 14724669. 2012. 00342.x (2012). 73. Wong, H. L. et al. Dynamics of archaea at fine spatial scales in Shark Bay mat microbiomes. Sci. Rep. 8, 46160. https:// doi. org/ 10. 1038/ srep4 6160 (2017). 74. Daffonchio, D. et al. Stratified prokaryote network in the oxicanoxic transition of a deep-sea halocline. Nature 440, 203–207. https:// doi. org/ 10. 1038/ natur e04418 (2006). Acknowledgements This study was supported by a grant from NASA’s Exobiology Program. The authors would like to thank the Bay Area Environmental Research (BAER) Institute for managing the postdoctoral fellowships awarded to P.M.M and M.D.L. We also thank representatives of Exportadora de Sal of Guerrero Negro, Baja California Sur. México for access to sites and assistance. We also thank our colleagues from México: José Q. García Maldonado, Jacob Alberto Valdivieso Ojeda, Santiago Cadena Rodriguez, Alejandro López Cortés and Hever Latisnere Barragán for field support and logistical assistance. Author contributions B.M.B., P.M.M. and M.D.L. conceptualized the study and/or supported the data analysis. P.M.M. performed the nucleic acid extractions. M.D.L. performed the bioinformatic analysis and P.M.M. drafted the manuscript. P.M.M., M.D.L. and B.M.B. all participated in writing, reviewing and editing the draft. All authors have read and approved the final manuscript. Funding This study was funded by NASA’s Exobiology Program. Ref. 17-EXO17-2-0134. Competing interests The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. Additional information Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s4159802452626-y. Correspondence and requests for materials should be addressed to P.M.-M. Reprints and permissions information is available at www.nature.com/reprints. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 17 Vol.:(0123456789) Scientific Reports | (2024) 14:2561 | https://doi.org/10.1038/s41598-024-52626-y www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. © The Author(s) 2024