Skip to main content
Download PDF

Genome-centric metagenomes unveiling the hidden resistome in an anchialine cave

Environmental Microbiome volume 19, Article number: 67 (2024) Cite this article

Abstract

Background

Antibiotic resistance is a critical global concern, posing significant challenges to human health and medical treatments. Studying antibiotic resistance genes (ARGs) is essential not only in clinical settings but also in diverse environmental contexts. However, ARGs in unique environments such as anchialine caves, which connect both fresh and marine water, remain largely unexplored despite their intriguing ecological characteristics.

Results

We present the first study that comprehensively explores the occurrence and distribution of ARGs and mobile genetic elements (MGEs) within an anchialine cave. Utilizing metagenomic sequencing we uncovered a wide array of ARGs with the bacitracin resistance gene, bacA and multidrug resistance genes, being the most dominant. The cave’s microbial community and associated resistome were significantly influenced by the salinity gradient. The discovery of novel β-lactamase variants revealed the cave’s potential as a reservoir for previously undetected resistance genes. ARGs in the cave demonstrated horizontal transfer potential via plasmids, unveiling ecological implications.

Conclusions

These findings highlight the need for further exploration of the resistome in unique environments like anchialine caves. The interconnected dynamics of ARGs and MGEs within anchialine caves offer valuable insights into potential reservoirs and mechanisms of antibiotic resistance in natural ecosystems. This study not only advances our fundamental understanding but also highlights the need for a comprehensive approach to address antibiotic resistance in diverse ecological settings.

Introduction

Antibiotic resistance poses a significant global health challenge, complicating the treatment of microbial infections and resulting in an estimated 4.95 million deaths in 2019 [1]. Current projections predict this number could escalate to 10 million by 2050 if no proactive measures are taken. It is thus imperative to closely monitor the emergence and prevalence of antibiotic resistance to safeguard public health and preserve the effectiveness of antibiotic [2]. Several factors contribute to the development of antibiotic resistance, such as the misuse and overuse of antibiotics, poor infection control practices, and genetic exchanges between bacteria [3, 4]. Central to this genetic interchange are the antibiotic resistance genes (ARGs), that confer antibiotic resistance traits to the bacteria carrying them. Through a process of horizontal gene transfer (HGT), mobile genetic elements (MGEs) can spread genes encoding for antibiotic resistance within and between bacterial populations, enabling rapid dissemination of resistance [5, 6]. Although there is a widespread recognition that the emergence and spread of ARGs are closely linked to human activities [7], bacteria have been developing resistance mechanisms for millions of years through natural selection [8].

ARGs are detectable across a range of environments such as soil, air, water, and biofilms, with aquatic environments being recognized as major reservoirs of ARGs [9,10,11]. Additionally, aquatic ecosystems play a central role in the cycle of antibiotic resistance, aligning with the One Health concept, which emphasizes the interdependence and mutual influence of human, animal, and environmental health. Anchialine caves, also known as subterranean estuaries, represent unique aquatic ecosystems characterized by a complex mixture of freshwater and seawater, resulting in stratified water bodies with distinct salinity gradients [12]. These gradients create diverse ecological niches that host a wide array of microbial communities including halophilic, halotolerant, and freshwater-adapted species that have partially been analysed in some caves [13,14,15,16]. Due to their isolated nature, anchialine caves are minimally impacted by anthropogenic activities, providing an ideal environment for studying ARGs in a less disturbed context. Anchialine caves are considered extreme habitats due to the varying physico-chemical conditions, including salinity, temperature, dissolved oxygen, nutrient availability, and limited light penetration. These factors can act as selective pressures on the resident organisms, encouraging the evolution of novel genetic adaptations [17]. Although anchialine caves are important habitats, they are still one of the least studied environments [18], with no existing data on the presence of ARGs.

We collected water samples from Sarcophagus Cave, an anchialine cave situated on the Croatian karst near the coast. The cave is characterized by a challenging entrance access and is minimally impacted by anthropogenic activities. Through the utilization of metagenomic sequencing, we present the first study addressing the existing research gap by conducting a comprehensive investigation of the presence and distribution of ARGs and MGEs in anchialine cave ecosystems. We hypothesize that these unique cave ecosystems serve as reservoirs for novel ARGs and MGEs. The main objectives of this study are: (i) to assess the abundance, diversity, and distribution of ARGs and MGEs across the salinity gradient within the anchialine cave, (ii) to identify the ARG-carrying microbial hosts, (iii) to elucidate the potential role of HGT in the dissemination of ARGs within these unique ecosystems. The outcomes of this study will significantly contribute to our understanding of antibiotic resistance in anchialine environments and shed light on the potential presence of novel resistance genes, while also addressing environmental concerns.

Materials and methods

Sampling site, physical and chemical measurements

The sampling site was an anchialine cave (Sarcophagus cave) located in the eastern coastal area of the Adriatic Sea near the Krka River estuary, Croatia (Fig. 1). It has a small vertical entrance shaft located at about 580 m distance from the estuary shore at an elevation of 5 m from the sea level. The entire depth of the cave is roughly 18 m, and the water body’s known depth is approximately 12 m. There are no visible signs of anthropogenic influence in the surrounding area, and the environment is considered minimally impacted.

Fig. 1
figure 1

Study area. Location map of Croatia, Šibenik, and the Sarcophagus Cave in Čapljina. Water samples were collected in August 2021 from the anchialine cave at six different depths (SC1-SC6), ranging from 0 to 12 m, along the salinity gradient. The cave’s entrance is inaccessible to most, and skilled divers are needed to collect the samples. The topology of the anchialine cave was adapted from the plan of Vedran Jalžić

Full size image

A professional cave scuba diver collected water samples (5 L) in a Niskin bottle at six depths along the salinity gradient at 0 (SC1), 2 (SC2), 4 (SC3), 7 (SC4), 10 (SC5), and 12 (SC6) m depths in August, 2021 (Table S1). In-situ measurements of physical parameters, including temperature, salinity, conductivity, dissolved oxygen, pH was carried out along different depths using a diver-operated Multisensor CTD probe (EXO2, YSI, USA). Subsequently, water samples intended for chemical analysis were gathered in acid-cleaned LDPE (Nalgene) containers and stored at -20 °C until analysis, a procedure executed within a month. DOC was analyzed using the HACH QBD1200 analyzer. Samples for the analysis of reduced sulfur species (RSS) were collected without exposure to oxygen and analyzed within 24 h. The total reduced sulfur species (RSStot) were analyzed by electrochemical methods at the hanging mercury electrode (663 VA Stand) connected to the potentiostat (PG STAT 128N, Metrhom, Netherlands) as previously described [19, 20]. More detailed information about the sampling site and analyses can be found in Kajan et al. [21].

Sample processing and DNA extraction

Each water sample (1–4 L) was filtered onto 0.22 µm pore-size polycarbonate filters, Whatman Nuclepore Track-Etch Membrane, diam. 47 mm, and then processed to extract total genomic DNA using the DNeasy PowerWater Kit (Qiagen, Inc., Valencia, CA, USA) following the manufacturer’s instruction. DNA samples were analyzed using 16S rRNA gene amplicon sequencing to define the microbial community and shotgun sequencing to characterize the resistome, which encompasses the total content of ARGs within the microbial community, as well as to investigate the presence of MGEs, and to retrieve the metagenome assembled genomes (MAGs).

Amplicon sequencing

A fraction of the isolated DNA was employed to amplify the hypervariable V4 region of the prokaryotic 16S rRNA gene using primer pair 515F Parada (5′-GTG YCA GCM GCC GCG GTA A-3′) and 806R Apprill (5′ GGA CTA CNV GGG TWT CTA AT-3′) [22, 23], modified with two 16 bp sequences which allow for sample barcoding in a second PCR step as described in detail in [24]. Amplicon sequencing was performed using Illumina MiSeq System in paired-end mode (v3 chemistry, 2 × 300 bp; Illumina, San Diego, CA, United States) at the Joint Microbiome Facility of the Medical University of Vienna and the University of Vienna.

Sequence data were processed in R using DADA2 following the workflow by Callahan et al. [25]. For the microbial community composition, the amplicon sequence variants (ASVs) were inferred across all samples in pooled mode. Further details on sequence trimming and settings for quality filtering are described in Pjevac et al. [24]. Taxonomic classification was performed by mapping 16S V4 ASV sequences against the SILVA SSU Ref NR 99 database (v. 138.1). Sequences identified as prokaryotic taxa in the 16S V4 region were retained for further analysis. Before analysis, chloroplastic and mitochondrial sequences were excluded and the ASV table was rarefied.

Shotgun metagenomic sequencing

Short-reads-based ARG annotation

Shotgun metagenome sequencing was performed using Illumina—NovaSeq 6000 at the Eurofins Genomics Europe Sequencing GmbH, Germany.

To identify potential ARG sequences using short reads as input data, DeepARG [26] was run using the short reads mode and setting the following parameters: –arg-alignment-identity 90.0 –min-prob 0.8 –arg-alignment-evalue 1e−10. The relative abundance of ARGs was calculated by normalizing the absolute abundance of ARGs to the 16S rRNA content obtained from amplicon sequencing in each sample. A graphical overview of the steps involved in the ARG and MGE annotation is presented in Fig. S1.

Assembly–based analysis for ARGs and MGEs annotation

Trimmed reads were obtained through cutadapt (v. 2.10) [27] and assembled using megahit (v. 1.1.2) [28], while contigs shorter than 1 kbp were removed using seqtk (v. 1.3) (github.com/lh3/seqtk). Metagenomic reads were mapped (68 Gbp of mapped data) to the assemblies using minimap2 (v. 2.17) [29]. The calculated read mappings were then converted using samtools (v. 1.11) [30] for binning using metabat2 (v. 2.15) [31]. We co-assembled all generated metagenomes and performed joint binning into combined MAGs following the same method to optimize the binning result. The quality of the genomes was checked using QUAST (v. 5.0.2) [32], CheckM (v. 1.1.1) [33], and classified using GTDBtk (v. 1.5.0) [34]. MAGs with completeness ≥ 50% and contamination ≤ 10% were selected for further analysis (770,953,733 bp bin size). Genome abundance was determined by the proportion of total reads mapping to the contigs in each bin, with counts normalized to the depth of the sequenced metagenomes.

To identify ARGs using MAGs as input data, three different types of software were used to screen contigs from MAGs for the presence of antibiotic resistance genes. First, Resistance Gene Identifier (RGI) software [35] which uses reference data from the Comprehensive Antibiotic Resistance Database (CARD) and detects homolog sequences by DIAMOND was run with default settings, but including loose hits to possibly detect novel ARGs. Second, we choose a deep learning method, DeepARG, with the following options “–model LS”, “–min-prob 0.8” and “–arg-alignment-identity 50” [26]. The amino acid sequences from the DeepARG candidate ARGs were retrieved manually by inspecting the genes predicted by Prodigal [36]. Finally, farGene [37], based on Hidden Markov models (HMMs), was applied to the same dataset to discover other drug classes and to confirm β-lactamase candidate sequences. All obtained candidate ARG-carrying sequences were manually checked for conserved domains using the NCBI’s web conserved domain search (CD-Search) with default settings (Table S6).

For the prediction of plasmid sequences in ARG-carrying contigs we used PlasFlow [38] Galaxy Version 1.0 under default settings (threshold for probability filtering = 0.7). Integron finder (with option—local-max included) was used to screen for the possible presence of complete integrons, ln0 (integron-integrase without gene cassettes) and CALIN elements (if composed of at least 2 attC sites) [39]. The criterion for an ARG to be part of an integron element was its presence within 12 kb ARG-ORF flankings as suggested by previous work [40]. Annotation of insertion sequences in ARG-carrying contigs was done with ISEScan [41]. We treated an insertion sequence as relevant for the mobility of an ARG if it was detected at a maximum 5 kb distance downstream or upstream of an ARG–ORF according to previous studies [42].

We conducted an additional analysis using the mobileOG-db (the curated and homologues set of mobile elements) to expand the search for MGE-related proteins. We performed a BLASTp search with and e-value threshold of 1e − 10, using all translated proteins as queries. Subsequently, we filtered the results based on the best e-value, query coverage (≥ 80%), and percent identity (≥ 80%).

β-lactamase and MAGs phylogenetic tree construction

All candidate β-lactamase sequences detected by three types of software (RGI-CARD software, DeepARG and farGene) were aligned with lactamase representative sequences (downloaded from the KEGG database) using MAFFT [43] implemented in Geneious® 9.1.8. Maximum-likelihood phylogenetic analysis was performed using IqTree v2.0. The best—fit model for alignment was selected and the branch supports were assessed with ultrafast bootstrap approximation (1000 replicates). FigTree was used for the tree visualization. To search for the best hits, we used detected β-lactamase sequences as queries for BLASTp through NCBI GenBank Database and through a local CARD database. For the β-lactamase carrying MAGs and other representative MAGs a phylogenetic tree was constructed using the classify_wf workflow in the Genome Taxonomy Database Toolkit (GTDB-TK) [34].

Statistical analysis

Data tables were imported into R studio v 2022.12.0 [44] and R v 4.2.2 [45] for statistics and visualization using the packages vegan [46], MASS [47], dplyr [48], ggplot2 [49], ggVennDiagram [50], pheatmap [51], circlize [52], cowplot [53], ggdendro [54], and ggpubr [55]. First, we analyzed the microbial community and resistome obtained from read-based analysis, using ASV read abundance and ARG copy number per 16S rDNA gene copy numbers as input data. We calculated the richness, as number of different ASVs/ARGs per sample, and constructed a matrix of dissimilarities using the abundance-based Bray–Curtis index, as an estimator of beta-diversity. The matrix was used to plot samples in a tree after hierarchical linkage analysis. We investigated if the differences among samples were related to the salinity, applying a negative binomial generalized linear model (NB-GML) on richness, since it represented counted data, and a PERMANOVA on beta-diversity. Moreover, we calculated the total normalized ARG abundances per sample and tested them against the salinity gradient, using a liner model (LM). Similarly, we assessed (via LM) if ARG abundances summed per class of resistance were dependent on the salinity. The correlation between microbial community and resistome composition was explored by Mantel test.

For MAGs, we quantified the abundance of ARGs and MGEs based on the cumulative abundance of a particular ARG-containing MAG at a specific sampling depth. Prior to heatmap visualization using the pheatmap package in R, the abundance of ARGs was log-transformed and the clustering method “complete” by rows was used.

Pearson’s correlation was used to measure the strength of the linear relationship between resistome abundances and environmental parameters using function ggscatter (package ggpubr). Presented are only those correlation graphs that yielded statistically significant results. Additionally, variation partitioning analyses (VPA) was performed to quantify the relative contributions of temperature, salinity, dissolved oxygen, and sulfide concentration to microbial community with ARGs (UpSetVP package in R).

Results

Microbial community composition

We retrieved 4782 ASVs, covering 65 microbial phyla, 147 classes, 203 orders, 232 families and 247 genera (Table S2). The most abundant phylum, across all the samples, was Nanoarchaeota, followed by Campylobacterota, Pseudomonadota, Bacteroidota, Verrucomicrobiota and Crenarchaeota (total abundance > 5000 reads) (Table S2). Regarding the genera, Sulfurimonas had the highest abundance, followed by Candidatus Omnitrophus, Sulfurovum and Nitrosarchaeum (each one with a total abundance > 4000 reads; with 99,660 ASVs not identified at genus level) (Table S2). Richness of microbial community showed a growing trend, and it was significantly higher with the increase of salinity (NB-GLM: p = 8.26 × 10−5) (Fig. S2A). The salinity was also a relevant driver in shaping the microbial community composition and explained 49.6% of sample variance (PERMANOVA: p = 0.0028). Indeed, in the tree depicting sample composition, SC1 isolated in a branch (Fig. S2B). Then, SC2 separated from the other samples, with SC3 clustered with SC4 and SC5 with SC6 (Fig. S2B).

Resistome composition

ARGs were annotated in all samples. We identified 130 ARGs, belonging to 20 different resistance classes (Fig. 2A, Table S3). The resistance against the bacitracin was the most represented resistance class, across all the samples, followed by multidrug, peptide, glycopeptide, and rifamycin resistance (Fig. 2A, Table S3). Considering the single genes, bacA (encoding resistance against bacitracin) had the highest abundance, followed by rpoB2 (a multidrug resistance, MDR, gene), ugd (encoding resistance against peptides), ompR (an MDR gene) and vanR (encoding resistance against glycopeptides) (Fig. 2B, Table S3). Only the aminoglycoside resistance genes significantly varied according to the salinity gradient, decreasing at high salt concentrations (LM: p = 0.0386). No significant patterns were observed for the other resistance classes (LM: p ≥ 0.0606). Total ARG abundances were significantly lower for higher salinities (LM: p = 0.03205). Richness of the resistome did not significantly differ along the salinity gradient (NB-GLM: p = 0.113), even if a decreasing trend was observed (Fig. S3A). However, when looking at the beta-diversity, the salinity, as factor, explained 42.6% of the resistome composition (PERMANOVA: p = 0.0028). Again, in the compositional tree, samples clustered based on the salinity (Fig. S3B). The resistome composition resulted significantly and strongly correlated with the microbial community structure (MANTEL TEST: r = 0.88; p = 0.0028).

Fig. 2
figure 2

Relative abundances of antimicrobial resistance classes and genes per sample based on the reads analysis. The composition of samples according to A antimicrobial resistance classes and B top 25 most abundant ARGs. The top 25 most abundant ARG were filtered out, and based on their total sum, their relative abundance per depth was calculated

Full size image

Community stratification of ARG-carrying MAGs in the anchialine environment

In the Sarcophagus Cave’s anchialine system 418 MAGs (CheckM completeness ≥ 50%, contamination ≤ 10%) spanning 35 different phyla were identified (Fig. 3). We observed a shift in the community structure at the second sampling layer (SC2) in the cave, coinciding with the halocline and oxycline (Table S1). This resulted in strong community stratification, with anoxic layers dominated by MAGs attributed to Archaea, Patescibacteria, and Omnitrophota (Fig. S3).

Fig. 3
figure 3

Taxonomical composition on the phylum level of metagenome-assembled genomes (MAGs) which carry antibiotic-resistance genes (ARGs) at a specific depth (SC1-SC6)

Full size image

In this study, we aimed to explore the potential presence of antibiotic resistance genes (ARGs) in this unique ecosystem. Our initial screening for ARG presence across the 418 MAGs obtained from the depth profile was performed with 3 diverse tools, based on their distinct modes of search for ARG sequences. All strict CARD hits (290) were kept for further analysis, while from the loose hits only those confirmed with DeepARG and/or farGene were selected. This yielded 578 potential ARG sequences from 230 MAGs in the anchialine cave (Table S4). Compared to the overall MAG community, certain exceptions were observed, such as the phylum Aenigmatarchaeota, which lacked any representative carrying ARGs, and Iainiarchaoeta and Nanoarchaeota, with only one detected ARG each (Fig. S4). Interestingly, Campylobacterota was frequently annotated among the retrieved MAGs at all depths (except at SC1) and in MAGs harboring ARGs (Fig. S4). Furthermore, as depth increased, the proportion of ARG-containing MAGs with Pseudomonadota (synonym Proteobacteria), Bacteroidota and Actinomycetota (synonym Actinobacteriota) decreased (Fig. 3). On the other hand, the proportion of Patescibacteria and Omnitrophota increased at depths SC3–SC6 compared to SC1-SC2.

Among the ARG carriers, 26 phyla were identified, including three from the Archaea domain (Iainarchaeota, Nanoarchaeota, and Thermoplasmatota) (Fig. S5A). Patescibacteria (21.7%), Omnitrophota (20.4%), Pseudomonadota (19.6%), Bacteroidota (8.3%), Desulfobacterota (6.5%), Campylobacterota (3.9%), altogether accounted for 80% of total ARG-carrying MAGs. Particularly rich in ARGs were Desulfobacterota and Pseudomonadota (> 90% of their MAGs contained at least 1 ARG sequence). Most of the ARG-carrying MAGs (82%) were classified up to the family level, while 50% remained unclassified at the genus level (Table S4).

In terms of classes, ARG-carrying MAGs were assigned to a total of 61 classes. Bacterial classes such as Koll11 (Omnitrophota), Gammaproteobacteria (Pseudomonadota), ABY1 (Patescibacteria), Bacteroidia (Bacteroidota), Alphaproteobacteria (Pseudomonadota), and Paceibacteria (Patescibacteria) collectively carried 61% of all ARG-carrying MAGs in the anchialine environment (Fig. S5B).

Mechanisms and classes of antibiotic resistance in diverse phyla

Our analysis revealed a wide range of resistance classes carried by different phyla. Pseudomonadota (Alphaproteobacteria and Gammaproteobacteria) were found to carry genes encoding for the resistance against multidrug, fluoroquinolone, tetracycline, aminoglycoside, bacitracin, β-lactam, and macrolide-licosamide-streptogramin (MLS). In contrast, Patescibacteria and Omnitrophota had a nearly exclusive prevalence of the glycopeptide resistance class. Indeed, resistance against glycopeptide was the most dominant class (70% of ARG-carrying MAGs), followed by MDR (24% of ARG-carrying MAGs), aminoglycoside (21.7%), and fluoroquinolone/tetracycline (20.9% MAGs) resistance classes (Fig. 4).

Fig. 4
figure 4

The association between antibiotic resistance gene (ARG) resistance classes (according to CARD terminology) and taxonomic composition of the microbial community at the phylum level based on the ARG-carrying MAGs. Phyla of MAGs are displayed on the outer circos layer. Pseudomonadota are presented as classes (Gammaproteobacteria, Alphaproteobacteria and Magnetococcia). Archaea are presented as phyla Iainarchaeota, Nanoarchaeota, and Thermoplasmatota. MDR multidrug–resistance; confers resistance to at least three different drug classes)

Full size image

Multiple antibiotic resistance mechanisms were found with all detected ARGs falling into the following categories based on CARD terminology: antibiotic target alteration (246), antibiotic efflux (201), antibiotic inactivation (91), reduced permeability to antibiotic (18), antibiotic target protection (2), and a combination of these mechanisms (Table S5).

Diversity, distribution, and environmental drivers of ARGs

We obtained 578 potential ARGs which were manually checked for conserved domains (Table S6). Of the 578 ARGs identified, the majority (232 ARGs; 40%) were found in Pseudomonadota. These were split between Gammaproteobacteria (27%) and Alphaproteobacteria (13%). Patescibacteria accounted for 77 ARGs (13%), Omnitrophota for 49 ARGs (8.5%), and Desulfobacterota for 45 ARGs (7.8%). All other detected phyla represented less than 7% of the total number of ARG-carrying contigs (Table S4).

The top three most abundant ARGs present at all sampling depths were those from the vancomycin resistance gene cluster (vanTG and vanWI) and the efflux pump gene adeF (Table S4, Figs. 5, S6A). Overall, fourteen vancomycin resistance gene sequences were detected (they are part of the vancomycin resistance gene clusters).

Fig. 5
figure 5

Heatmap of antibiotic resistance genes (ARGs) for six sampling sites (SC1-SC6). The abundance of ARGs was quantified and normalized based on the abundance of a particular ARG-containing MAG at a specific sampling depth. Rows were clustered using the “complete” method with the pheatmap function (presented in Fig. S6A)

Full size image

The shallowest sampling depth (SC1) with surface freshwater, contains most ARGs and, in contrast to other depths, the most dominant resistance classes were aminoglycoside and MDR genes. At the second sampling depth, SC2, still above the halocline, resistances to peptide and fluoroquinolone/tetracycline were the dominant classes, while the other four sampling depths (under the halocline) were populated with glycopeptide, fluoroquinolone/tetracycline, and peptide resistance classes (Fig. 5). The diversity of ARGs also changes between different sampling depths, with SC1 having the most diverse ARG composition (with 19 unique ARGs) and SC6 the lowest (1 unique). Fifteen ARGs are shared across depths, representing a core anchialine resistome (Fig. S6B).

We used Pearson’s correlation analysis to investigate the relationship between salinity and ARGs across sampling depths. The analysis revealed a strong negative correlation between salinity and the relative abundance of ARGs in water layers, suggesting that the abundance of ARGs decreased as salinity levels increased (Pearson’s correlation; R =  − 0.95, p = 0.004, Fig. S7A). Moreover, there was a negative correlation between the relative abundance of ARGs and the richness of MAGs at the family level (Pearson’s correlation; R =  − 0.84, p = 0.037, Fig. S7B), while no correlation was observed at the genus level. Finally, we conducted VPA to quantify the relative contribution of salinity to the microbial community with ARGs (Fig. S8). VPA analysis revealed that the dissimilarity in the microbial community with ARGs along the depth was strongly related to salinity (43.2%).

The potential of MGEs in the dissemination of cave ARGs

In our initial analysis, we employed PlasFlow, Integron finder, and ISEScan to search for MGEs associated with ARG-carrying contigs. We found that a limited number of MGEs were detected, with 21 ARG-carrying contigs identified as of plasmid origin (3.9% of total contig count), 412 were of chromosomal origin and 105 were unclassified. The plasmid-associated ARGs were primarily attributed to Pseudomonadota and Bacillota (synonym Firmicutes) and primarily involved the antibiotic efflux resistance mechanism. We also identified a small number of CALIN elements (n = 4) lacking integrin—integrases (2–7 attc sites) associated with ARG-carrying sequences. Only in one case the detected CALIN element was located within the proposed 12 kb ARG flankings [40]. Complete integrons were not detected in ARG—carrying contigs. Insertion sequences (ISs) were found in 54 ARG—carrying contigs. In total, we classified 6 ARG (2 vanTG, 2 adeF, 1 vanWI, and 1 ugd) as IS-associated (11% of all detected ISs). To further explore the potential mobility of ARGs, we conducted an extended analysis using BLASTp against the mobileOG-db, a curated database of bacterial mobile genetic elements. Despite employing quite stringent e-value (1e-10), query coverage (80%), and identity (80%) thresholds, this approach revealed a much more comprehensive picture of MGEs within the ARG-carrying microbial community. We obtained 567 MGE-related protein sequences which were categorized into five main MGE categories: phage specific biological processes (175), integration and excision from one genetic locus to another (145), replication, recombination, or nucleic acid repair (98), element stability, transfer, or defense (77), and interorganism transfer (72) (Fig. 6A). There were 153 unique MGE-related proteins found in 37% of ARG-carrying MAGs (Fig. S9). The top 10 MGE-related proteins were hup, cas1, tnpA, groL, clpX, hsdM, hfq, clpA, thyA, with the unclassified (category integration/excision) being the most abundant (Table S7). They mostly belonged to Pseudomonadota (Gammaproteobacteria) and Campylobacterota (Fig. 6B). As in the case of ARGs, Pearson’s correlation analysis revealed a strong negative correlation between salinity and relative abundance of MGEs in water layers, suggesting that their abundance decreased as salinity levels increased (Pearson’s correlation; R =  − 0.89, p = 0.0016, Fig. S7C). Moreover, the absolute number of ARGs positively correlated with the presence or absence of a particular MGE type in a specific layer (Fig. S7D).

Fig. 6
figure 6

Distribution and taxonomic composition of mobile genetic elements (MGEs) in the anchialine cave. A Heatmap of MGE categories for six sampling sites (SC1–SC6). The abundance of MGE categories was quantified and normalized based on their abundance in a particular MAG at a specific sampling depth. Rows were clustered using the “complete” method with the pheatmap function. B The bar plots represent the relative abundance of bacterial phyla corresponding to the top 10 most abundant MGEs

Full size image

Potential novel β-lactamases and their relationship with host

Since β-lactamases are perceived as one of the most important antibiotic resistance genes considering human health risk, we analyzed the anchialine MAGs for their presence. We used farGene since it was shown to outcompete other similar programs for β-lactamase prediction. All farGene β-lactamase candidate sequences were confirmed in CARD hits. In total, 42 β-lactamases which belonged to 28 MAGs were found in the anchialine cave.

The β-lactamase–carrying MAGs belonged to 9 different phyla (Pseudomonadota, Bacteroidota, Acidobacteriota, Omnitrophota, Actinomycetota, Campylobacterota, Chlamydiota, Thermoplasmatota–Archaea, Verrucomicrobiota). Gammaproteobacteria (33%), Alphaproteobacteria (28,5%) and Bacteroidia (16,6%) were the most dominant classes. Approximately one-half (20/42) of β-lactamase-carrying MAGs were identified at the genus level and only one was annotated at the species level (A. xenomutans).

Among the detected β-lactamases 11, 21, 1, and 9 were classified into class A, B, C, and D, respectively (based on molecular Ambler’s classification system). A phylogenetic tree with β-lactamase representatives (reference sequences from the KEGG database) confirmed their position within the annotated classes (Fig. S10). Most β-lactamases belonged to the subclass B3 (18/42), with sequences similar to CAR-1, BJP-1, FEZ-1 and SPG-1 being the most frequent and present within Opitutaceae, Hyphomonadaceae_UBA7672, Vicinamibacterales, Caulobacter, Pseudohongiellaceae, Rhizomicrobium, and Alteromonas. Eleven β-lactamases were labeled as class A, among them one blaS1-like (64.6% sequence similarity with CARD ARO:3,005,111), a predominant β-lactamase in the species Mycobacterium smegmatis, was found in a Mycobaterium (Fig. 7). Class C β-lactamases were represented with only one sequence detected in Pararheinheimera (Enterobacterales), similar to ACT-20 from Enterobacter hormaechei (ARO:3,001,841). In total, three β-lactamases resistant to carbapenems (CAR-1, B3 class—30.1% identity with ARO:3,006,903; PER-6, class A class—44.4% identity with ARO: 3,002,368; and ACT-20, class C—54.5% identity with ARO: 3,001,841) were found in Enterobacterales, which are listed by WHO as critical priority for the development of new antimicrobials. Oxacillinases were the most dominant ARGs in class D β-lactamases and were found in Bacteroidales, Emticicia, Flavobacteriales_UBA10329, Jiella, Pseudohongiellaceae, Rhabdochlamydia, Sulfurimonadaceae, and Thiobacillaceae.

Fig. 7
figure 7

β-lactamases and MAGs phylogenetic trees. A β-lactamases phylogenetic tree based on maximum–likelihood (1000 ultrafast bootstrap approximation replicates) contains novel anchialine β-lactamase sequences and representative sequences from the Comprehensive Antibiotic Resistance Database (CARD). B Metagenome—assembled genomes from this study were aligned with the representative genomes and a phylogenetic tree was constructed using the Genome Taxonomy Database Toolkit (GTDB-Tk). β-lactamase sequences and MAGs detected in this study are coloured in blue

Full size image

Our comparison using CARD reference genes revealed that all identified β-lactamase sequences exhibited less than 70% pairwise identity with representative CARD β-lactamases (Table 1).

Table 1 Similarity ranges of detected β-lactamase sequences. CARD reference protein database was used for the blastp search
Full size table

When we searched for similarities using BLASTp and the NCBI database (Table S8), only two sequences from class A, namely Pararheinheimera and A. xenomutans, showed high pairwise identity (~ 99%) with two NCBI entries, MBN8446110.1 (Gammaproteobacteria bacterium) and WP_026948658.1 (Alcanivorax), respectively. Interestingly, almost half (19/42) of the β-lactamase sequences showed less than 70% of the identity to NCBI best hits. These diverse sequences were distributed across Pseudomonadota (12), Omnitrophota (2), Actinobacteriota (1), Bacteroidota (3), and Campylobacterota (1), as illustrated in Fig. 7.

Special attention should be given to β-lactamase sequences carried by MGEs [56]. The majority of the β-lactamase sequences in the anchialine cave were identified in MAGs containing MGE-related proteins. For example, the genera: Jiella, Mycobacterium, Alteromonas, PFJX01 (family Thiobacillaceae), Emticicia, and Pararheinheimera contained ten or more MGE-related proteins. Interestingly, a β-lactamase sequence similar to the OXA-198 (69% sequence identity with ARO:3,001,805) from the WHO-listed pathogen P. aeruginosa, was detected in a MAG belonging to Thiobacillaceae (g__PFJX01). In addition to the OXA-198, this particular MAG contained 13 different MGE-related proteins (belonging to categories: integration/excision, phage, replication/recombination/repair, and stability/transfer/defence) (Fig. 6A). Altogether, these results might suggest a potential for HGT of β-lactamase sequences within the cave environment.

Discussion

Antibiotic resistance genes (ARGs) have gained increased focus in recent years due to their potential implications for human health and economic stability [57, 58]. Antibiotic resistance is a natural and ancient phenomenon, and it is expected that the majority of resistance genes, including those still uncharacterized, are less likely to be found within pathogens or human commensals, rather they would be found in environmental bacteria [59, 60].

Anchialine cave environments are recognized for their distinctive and isolated ecological conditions, which stimulate unique biodiversity. Despite their scientific interest, these areas remain inadequately studied due to limited accessibility for researchers. Microbial communities in certain caves have been partially explored [13, 14, 61]. In the investigated cave, salinity played a crucial role in shaping the microbial community composition, with archaea being particularly abundant. Additionally, microbial diversity tended to increase with greater water depth and salinity, aligning partially with findings from previous studies [14].

The resistome of anchialine caves has yet to be comprehensively investigated. To our knowledge, this is the first study on the occurrence of ARGs, as well as the investigation and characterization (in terms of MGE content) of their host organisms at the metagenomic level. As such, it represents a significant contribution to our understanding of the prevalence and distribution of ARGs in these understudied environments. Read-based annotation revealed a clear pattern, with the bacA gene, which encodes bacitracin resistance, being significantly the most abundant. It was followed by MDR genes. These findings are consistent with results from other aquatic ecosystems, where, regardless of anthropogenic pressure, these genes are among the most prevalent ARGs [62, 63]. The third most abundant resistance genes were those encoding for peptide resistance, which were more prevalent in SC2, a transitional zone characterized by a combination of environmental stressors, where salinity begins to increase while dissolved oxygen decreases. This prevalence may be linked to the potential adaptive advantage provided by genes involved in target alteration and efflux resistance mechanisms. Notably, aminoglycoside resistance genes were the only class of resistance genes that differed significantly among the sampling sites, showing a clear inverse relationship between their abundance and salinity concentration, consistent with observations in other aquatic ecosystems [64].

As observed for the entire microbial community, salinity emerged as the primary driver of diversity in the ARG-carrying MAGs in the anchialine environment. Distinctions were noted among specific ARG-carrying MAGs within the water layers. Notably, Patescibacteria and Omnitrophota were primarily found beneath the halocline because of their adaptation to unique conditions such as higher salinity, low oxygen (SC3-SC6 are anoxic layers), darkness, and low nutrient content in the deeper layers of the cave [65, 66]. We found that these bacteria predominantly carry glycopeptide (vancomycin) resistance genes, strengthening the knowledge that Patescibacteria play an important role in spreading vancomycin resistance in groundwater environments [67]. There are contrasting findings regarding the influence of salinity on antibiotic resistance in different environments, which underlines the complexity of ARG dynamics [68,69,70]. The prevalence of the order Burkholderiales (Gammaproteobacteria), among the ARG-carrying contigs in our study aligns with their known abundance in fresh groundwaters [71]. Additionally, the abundance of ARGs (determined by read-based annotation) was influenced by the salinity levels showing lower abundance with increasing depth and salinity. This result was further confirmed by assembly analysis, showing a strong negative correlation between salinity and the abundance of ARGs in water layers. The freshwater and brackish niches at higher levels of the anchialine cave might be more influenced by groundwater, increasing the likelihood of introducing antibiotic-resistant strains from the soil. This can occur through the washing of upper layers of soil, known to be rich in antibiotic-resistant strains, and making it a potential source for the dissemination of such strains into aquatic ecosystems [66, 67]. Furthermore, this emphasizes the role of external environments in shaping the resistome of cave ecosystems.

Notably, within the resistome of the investigated cave, the discovery of β-lactamases in 9 diverse phyla (10 classes) inhabiting the anchialine cave habitat, together with low pairwise identity between the detected β-lactamases and previously identified sequences (< 70% pairwise identity), suggests that the cave may be a rich source of undiscovered and novel β-lactamase variants. Due to the substantial representation of unclassified genera in our data, we opted to utilize the Hidden Markov models farGene software [37] to broaden the possibility of finding novel lactamase sequences. This approach proved to be successful as it enabled the detection of β-lactamase sequences that would have otherwise been undiscovered if we had used only the more stringent (strict hits-based) similarity approach with the CARD-RGI tool. The detection of carbapenem-resistant β-lactamases (Alteromonas; CAR-1-like) in Enterobacterales [72] species within the anchialine cave is of particular interest. Carbapenem-resistant β-lactamases from Enterobacterales (CRE) can originate from various sources, including healthy individuals, soil, water, plants, dairy products, and raw meat [73], and are listed as critical priority pathogens [74]. Although the presence of these genes has been mostly studied in hospital environments and wastewater [75, 76], it is now well-established that CRE can be found in various aquatic environments such as rivers, lakes, marine and recreational water bodies [77, 78]. The significance of monitoring CRE in these environments to evaluate their spread in diverse environmental communities has already been brought to the attention [79]. Our finding of CRE in this isolated and unique ecosystem emphasizes the importance of continued research to unravel their role in natural ecosystems [80]. Additionally, the detection of another carbapenemase, OXA-198-like sequence, in non-pathogenic bacteria from the family Thiobacillaceae, which was similar (pairwise identity 69%) to the lactamase from the pathogenic bacteria Pseudomonas aeruginosa [81], suggests that non-pathogenic bacteria can be reservoirs for resistance genes. Although β-lactamase sequences identified in this study were found to be similar to other β-lactamases, their real resistance and potential environmental or human health risks cannot be determined solely from sequence similarity [82]. Despite this limitation, the discovery of novel β-lactamase variants in the anchialine cave presents an opportunity for natural products discovery and bioprospecting of new antibiotics, as it has been demonstrated that there is a high possibility of finding novel antimicrobial compounds in microbes that inhabit underground environments [83,84,85]. Our findings strengthen the argument for exploring unique habitats like anchialine caves for novel antimicrobial compounds.

It is recognized that mobile genetic elements (MGEs) are key facilitators of horizontal gene transfer (HGT) among bacteria, a process that plays a significant role in the widespread distribution of antibiotic resistance genes [5, 86]. Within our initial screening of MGEs, we found a minor proportion of ARGs (3.9% of ARG-carrying contigs) in the anchialine cave associated with plasmids, while the extended analysis using mobileOG-db, revealed a more comprehensive view of MGEs within the ARG-carrying MAGs. Our finding of 153 unique MGE-related proteins found in 37% of ARG-carrying MAGs suggests there is a significant potential for ARG mobility within the microbial community. Among the identified MGE categories, phage-specific biological processes and integration and excision mechanisms were predominant, in line with previous studies [87, 88]. Pseudomonadota and Campylobacterota were the primary carriers of MGE-related proteins, also comparable with recent studies [87, 89]. The presence of a wide variety of MGEs in the anchialine cave, and plasmids among them, suggests the potential horizontal transfer and spread of ARGs in this particular environment.

While the association between MGEs and ARGs is well established in environments impacted by human activity [90, 91], their presence and distribution patterns in natural environments have been less thoroughly explored. Some Antarctic genome sequencing studies were unable to detect any presence of MGE [92] and references therein], while others found a limited proportion even when using a hybrid approach for assembling MAGs [93]. These variations underscore the necessity for further investigation of the role of MGEs in the spread of antibiotic resistance within diverse natural environments. Our results offer a preliminary understanding of the mobility potential of ARGs in an anchialine environment. Additional techniques, such as long-read sequencing could help reveal in more detail the mobility potential of the anchialine ARGs and the mechanisms driving their spread within the microbial community.

Conclusions

This study substantially advances our understanding of the composition of the resistome associated with microbial communities living in the specific ecological settings of anchialine caves. The bacterial community, ARG-carrying MAGs, and the resistome showed distinct responses to the salinity gradient indicating adaptive strategies developed to thrive in these environmental niches. This finding underscores the anchialine cave’s potential to function as a significant reservoir for clinically relevant resistance traits, necessitating a need to evaluate broader ecological implications. The presence of a wide variety of MGEs, including plasmids, implies that HGT may facilitate the exchange and development of novel resistance traits. Furthermore, we showed that anchialine caves are potential reservoirs of novel β-lactamase variants, emphasizing the importance of these unique ecosystems as a source of previously undetected ARGs. Given the rising global threat of antibiotic resistance, our findings highlight the critical importance of research in such distinct environments to inform our understanding of the emergence and spread of ARGs. Future research should focus on employing advanced methodologies, such as long-read sequencing, and on more extensive sampling to elucidate the mechanisms driving the spread of ARGs within these distinctive ecosystems.

Availability of data and materials

The raw reads of 16S rRNA amplicon sequencing (SRX22885116, SRX22885117, SRX22885118, SRX22885119, SRX22885121 and SRX22885122) and the shotgun sequencing (SRX22885113, SRX22885126, SRX22885128, SRX22885130, SRX22885132 and SRX22885135) can be accessed under the Bioproject accession number PRJNA1052463. A simplified version of the R script is available online https://github.com/kkajan/ARG_anch_cave.

Abbreviations

ARG:

Antibiotic resistance gene

MGE:

Mobile genetic element

MAG:

Metagenome-assembled genome

HGT:

Horizontal gene transfer

References

  1. Murray CJ, Ikuta KS, Sharara F, Swetschinski L, Robles Aguilar G, Gray A, et al. Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis. Lancet. 2022;399:629–55.

    Article  CAS  Google Scholar 

  2. World Health Organization. Global antimicrobial resistance and use surveillance system (GLASS) Report 2022. Braz Dent J. 2022.

  3. Collignon P, Beggs JJ, Walsh TR, Gandra S, Laxminarayan R. Anthropological and socioeconomic factors contributing to global antimicrobial resistance: a univariate and multivariable analysis. Lancet Planet Heal. 2018;2:e398-405.

    Article  Google Scholar 

  4. Ellabaan MMH, Munck C, Porse A, Imamovic L, Sommer MOA. Forecasting the dissemination of antibiotic resistance genes across bacterial genomes. Nat Commun. 2021;12:2435.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Partridge SR, Kwong SM, Firth N, Jensen SO. Mobile genetic elements associated with antimicrobial resistance. Clin Microbiol Rev. 2018;31:e00088-e117.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. von Wintersdorff CJH, Penders J, Van Niekerk JM, Mills ND, Majumder S, Van Alphen LB, et al. Dissemination of antimicrobial resistance in microbial ecosystems through horizontal gene transfer. Front Microbiol. 2016;7:173.

    Google Scholar 

  7. Larsson DGJ, Flach CF. Antibiotic resistance in the environment. Nat Rev Microbiol. 2022;20:257–69.

    Article  CAS  PubMed  Google Scholar 

  8. Davies J, Davies D. Origins and evolution of antibiotic resistance. Microbiol Mol Biol Rev. 2010;12:9–16.

    Google Scholar 

  9. Makowska-Zawierucha N, Mokracka J, Małecka M, Balazy P, Chełchowski M, Ignatiuk D, et al. Quantification of class 1 integrons and characterization of the associated gene cassettes in the high Arctic–interplay of humans and glaciers in shaping the aquatic resistome. Ecol Indic. 2022;145:109633.

    Article  CAS  Google Scholar 

  10. Zhang S, Yang G, Hou S, Zhang T, Li Z, Liang F. Distribution of ARGs and MGEs among glacial soil, permafrost, and sediment using metagenomic analysis. Environ Pollut. 2018;234:339–46.

    Article  CAS  PubMed  Google Scholar 

  11. Singh AK, Kaur R, Verma S, Singh S. Antimicrobials and antibiotic resistance genes in water bodies: pollution, risk, and control. Front Environ Sci. 2022;10:830861.

    Article  Google Scholar 

  12. Bishop RE, Humphreys WF, Kršinic F, Sket B, Iliffe TM, Zic V, et al. “Anchialine” redefined as a subterranean estuary in a crevicular or cavernous geological setting. J Crustac Biol. 2015;35:511–4.

    Article  Google Scholar 

  13. Brankovits D, Pohlman JW, Niemann H, Leigh MB, Leewis MC, Becker KW, et al. Methane-and dissolved organic carbon-fueled microbial loop supports a tropical subterranean estuary ecosystem. Nat Commun. 2017;8:1835.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Elbourne LDH, Sutcliffe B, Humphreys W, Focardi A, Saccò M, Campbell MA, et al. Unravelling stratified microbial assemblages in australia’s only deep anchialine system, The Bundera Sinkhole. Front Mar Sci. 2022;9:1–15.

    Article  Google Scholar 

  15. Kajan K, Cukrov N, Cukrov N, Bishop-Pierce R, Orlić S. Microeukaryotic and Prokaryotic Diversity of Anchialine Caves from Eastern Adriatic Sea Islands. Microb Ecol. 2022;83:257–70.

    Article  CAS  PubMed  Google Scholar 

  16. Bhullar K, Waglechner N, Pawlowski A, Koteva K, Banks ED, Johnston MD, et al. Antibiotic resistance is prevalent in an isolated cave microbiome. PLoS ONE. 2012;7:e34953.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Goh KM, González-Siso MI, Sani RK. Genomics of extreme environments: unveiling the secrets of survival. Sci Rep. 2023;13:21441.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Calderón-Gutiérrez F, Sánchez-Ortiz CA, Huato-Soberanis L. Ecological patterns in anchialine caves. PLoS ONE. 2018;13:1–19.

    Article  Google Scholar 

  19. Bura-Nakić E, Helz GR, Ciglenečki I, Ćosović B. Reduced sulfur species in a stratified seawater lake (Rogoznica Lake, Croatia); seasonal variations and argument for organic carriers of reactive sulfur. Geochim Cosmochim Acta Pergamon. 2009;73:3738–51.

    Article  Google Scholar 

  20. Marguš M, Morales-Reyes I, Bura-Nakić E, Batina N, Ciglenečki I. The anoxic stress conditions explored at the nanoscale by atomic force microscopy in highly eutrophic and sulfidic marine lake. Cont Shelf Res Pergamon. 2015;109:24–34.

    Article  Google Scholar 

  21. Kajan K, Kirkegaard R, Pjevac P, Orlić S, Mehrshad M. Niche and spatial partitioning restrain ecological equivalence among microbes along aquatic redox gradient. BioRxiv. 2024. https://doi.org/10.1101/2024.08.09.607300.

    Article  Google Scholar 

  22. Parada AE, Needham DM, Fuhrman JA. Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ Microbiol. 2016;18:1403–14.

    Article  CAS  PubMed  Google Scholar 

  23. Apprill A, Mcnally S, Parsons R, Weber L. Minor revision to V4 region SSU rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquat Microb Ecol. 2015;75:129–37.

    Article  Google Scholar 

  24. Pjevac P, Hausmann B, Schwarz J, Kohl G, Herbold CW, Loy A, et al. An economical and flexible dual barcoding, two-step PCR approach for highly multiplexed amplicon sequencing. Front Microbiol. 2021;12:669776.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Callahan BJ, Sankaran K, Fukuyama JA, McMurdie PJ, Holmes SP. Bioconductor workflow for microbiome data analysis: from raw reads to community analyses. F1000Research. 2016;5:1492.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Arango-Argoty G, Garner E, Pruden A, Heath LS, Vikesland P, Zhang L. DeepARG: a deep learning approach for predicting antibiotic resistance genes from metagenomic data. Microbiome. 2018;6:23.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011;17:10–2.

    Article  Google Scholar 

  28. Li D, Liu CM, Luo R, Sadakane K, Lam TW. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics. 2015;31:1674–6.

    Article  CAS  PubMed  Google Scholar 

  29. Li H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics. 2018;34:3094–100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Kang DD, Li F, Kirton E, Thomas A, Egan R, An H, et al. MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ. 2019;7:e7359.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. Bioinformatics. 2013;29:1072–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25:1043–55.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Chaumeil PA, Mussig AJ, Hugenholtz P, Parks DH. GTDB-Tk: a toolkit to classify genomes with the genome taxonomy database. Bioinformatics. 2020;36:1925–7.

    Article  CAS  Google Scholar 

  35. Alcock BP, Raphenya AR, Lau TTY, Tsang KK, Bouchard M, Edalatmand A, et al. CARD 2020: antibiotic resistome surveillance with the comprehensive antibiotic resistance database. Nucleic Acids Res. 2020;48:D517–25.

    CAS  PubMed  Google Scholar 

  36. Hyatt D, Chen G-L, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Integrated nr database in protein annotation system and its localization. BMC Bioinf. 2010;11.

  37. Berglund F, Österlund T, Boulund F, Marathe NP, Larsson DGJ, Kristiansson E. Identification and reconstruction of novel antibiotic resistance genes from metagenomes. Microbiome Microbiome. 2019;7:1–14.

    Google Scholar 

  38. Krawczyk PS, Lipinski L, Dziembowski A. PlasFlow: predicting plasmid sequences in metagenomic data using genome signatures. Nucleic Acids Res. 2018;46:E35.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Néron B, Littner E, Haudiquet M, Perrin A, Cury J, Rocha EPC. IntegronFinder 2.0: identification and analysis of integrons across bacteria, with a focus on antibiotic resistance in Klebsiella. Microorganisms. 2022;10:700.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Nielsen TK, Browne PD, Hansen LH. Antibiotic resistance genes are differentially mobilized according to resistance mechanism. Gigascience. 2022;11:1–17.

    Article  CAS  Google Scholar 

  41. Xie Z, Tang H. ISEScan: automated identification of insertion sequence elements in prokaryotic genomes. Bioinformatics. 2017;33:3340–7.

    Article  CAS  PubMed  Google Scholar 

  42. Che Y, Yang Y, Xu X, Brinda K, Polz MF, Hanage WP, et al. Conjugative plasmids interact with insertion sequences to shape the horizontal transfer of antimicrobial resistance genes. Proc Natl Acad Sci U S A. 2021;118:e2008731118.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Katoh K, Misawa K, Kuma K, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. RStudio Team. RStudio: integrated development environment for R [Internet]. Boston, MA: PBC; 2020. Available from: http://www.rstudio.com/.

  45. R Core T. A Language and environment for statistical computing R. 2020; Available from: www.R-project.org.

  46. Oksanen r J, Simpson GL, Blanchet FG, Solymos P, Stevens MHH, Szoecs E, et al. Title community ecology package. Vegan Community Ecol Packag. 2022.

  47. Venables WN, Ripley BD. Statistics and computing: editorial. Mod Appl Stat with S [Internet]. 2002;12:7. Available from: https://www.stats.ox.ac.uk/pub/MASS4/.

  48. Wickham H, François R, Henry L, Müller K. A grammar of data manipulation [R package dplyr version 1.0.7]. Media. 2021.

  49. Wickham H. ggplot2 - Elegant graphics for data analysis | Hadley Wickham | Springer. Springer. 2017.

  50. Gao CH, Yu G, Cai P. ggVennDiagram: an intuitive, easy-to-use, and highly customizable R package to generate venn diagram. Front Genet. 2021;12:706907.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Kolde R. Package “pheatmap”: Pretty heatmaps. Version 1012. 2019. Available from: https://cran.r-project.org/web/packages/pheatmap/.

  52. Gu Z, Gu L, Eils R, Schlesner M, Brors B. Circlize implements and enhances circular visualization in R. Bioinformatics. 2014;30.

  53. Wilke C. cowplot: streamlined plot theme and plot annotations for “ggplot2”. R package version 1.1.3. 2024. Available from: https://wilkelab.org/cowplot/.

  54. Vries AD, Ripley BD. Package ‘ggdendro’ create dendrograms and tree diagrams using “ggplot2.” 2024; Available from: https://cran.r-project.org/web/packages/ggdendro/.

  55. Kassambara A. ggpubr: “ggplot2” based publication ready plots. R package version 0.2. https://CRANR-project.org/package=ggpubr. 2020.

  56. Bush K, Bradford PA. Epidemiology of β-lactamase-producing pathogens. Clin Microbiol Rev. 2020;33:1–37.

    Article  Google Scholar 

  57. Zhang Z, Zhang Q, Wang T, Xu N, Lu T, Hong W, et al. Assessment of global health risk of antibiotic resistance genes. Nat Commun. 2022;13:1–11.

    Google Scholar 

  58. Akram F, Imtiaz M, Haqul I. Emergent crisis of antibiotic resistance: a silent pandemic threat to 21st century. Microb Pathog. 2023;174:105923.

    Article  CAS  PubMed  Google Scholar 

  59. D’Costa VM, King CE, Kalan L, Morar M, Sung WWL, Schwarz C, et al. Antibiotic resistance is ancient. Nature. 2011;477:457–61.

    Article  PubMed  Google Scholar 

  60. Allen HK, Donato J, Wang HH, Cloud-Hansen KA, Davies J, Handelsman J. Call of the wild: antibiotic resistance genes in natural environments. Nat Rev Microbiol. 2010;8:251–9.

    Article  CAS  PubMed  Google Scholar 

  61. Gonzalez BC, Iliffe TM, Macalady JL, Schaperdoth I, Kakuk B. Microbial hotspots in anchialine blue holes: initial discoveries from the Bahamas. Hydrobiologia. 2011;677:149–56.

    Article  CAS  Google Scholar 

  62. Sivalingam P, Sabatino R, Sbaffi T, Corno G, Fontaneto D, Borgomaneiro G, et al. Anthropogenic pollution may enhance natural transformation in water, favouring the spread of antibiotic resistance genes. J Hazard Mater. 2024;475:134885.

    Article  CAS  PubMed  Google Scholar 

  63. Xie X, Chen B, Zhu S, Yang R, Yuan K, Yang Y, et al. Comparative analysis of characteristics of antibiotic resistomes between Arctic soils and representative contaminated samples using metagenomic approaches. J Hazard Mater. 2024;469:133943.

    Article  CAS  PubMed  Google Scholar 

  64. Liang H, Wang F, Mu R, Huang J, Zhao R, Li X, et al. Metagenomics analysis revealing the occurrence of antibiotic resistome in salt lakes. Sci Total Environ. 2021;790:148262.

    Article  CAS  PubMed  Google Scholar 

  65. Tian R, Ning D, He Z, Zhang P, Spencer SJ, Gao S, et al. Small and mighty: adaptation of superphylum Patescibacteria to groundwater environment drives their genome simplicity. Microbiome Microbiome. 2020;8:1–15.

    CAS  Google Scholar 

  66. Herrmann M, Wegner CE, Taubert M, Geesink P, Lehmann K, Yan L, et al. Predominance of Cand. Patescibacteria in groundwater is caused by their preferential mobilization from soils and flourishing under oligotrophic conditions. Front Microbiol. 2019;10:1–15.

    Article  Google Scholar 

  67. Maatouk M, Ibrahim A, Rolain J-M, Merhej V, Bittar F. Small and equipped: the rich repertoire of antibiotic resistance genes in candidate phyla radiation genomes. MSystems. 2021;6:e0089821.

    Article  PubMed  Google Scholar 

  68. Huang J, Zhu J, Liu S, Luo Y, Zhao R, Guo F, et al. Estuarine salinity gradient governs sedimentary bacterial community but not antibiotic resistance gene profile. Sci Total Environ. 2022;806:151390.

    Article  CAS  PubMed  Google Scholar 

  69. Bergeron S, Brown R, Homer J, Rehage S, Boopathy R. Presence of antibiotic resistance genes in different salinity gradients of freshwater to saltwater marshes in southeast Louisiana, USA. Int Biodeterior Biodegrad. 2016;113:80–7.

    Article  CAS  Google Scholar 

  70. Zhang YJ, Hu HW, Yan H, Wang JT, Lam SK, Chen QL, et al. Salinity as a predominant factor modulating the distribution patterns of antibiotic resistance genes in ocean and river beach soils. Sci Total Environ. 2019;668:193–203.

    Article  CAS  PubMed  Google Scholar 

  71. Ruiz-González C, Rodríguez-Pie L, Maister O, Rodellas V, Alorda-Keinglass A, Diego-Feliu M, et al. High spatial heterogeneity and low connectivity of bacterial communities along a Mediterranean subterranean estuary. Mol Ecol. 2022;31:5745–64.

    Article  PubMed  PubMed Central  Google Scholar 

  72. Caliskan-Aydogan O, Alocilja EC. A review of carbapenem resistance in Enterobacterales and its detection techniques. Microorganisms. 2023;11:1491.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Tilahun M, Kassa Y, Gedefie A, Ashagire M. Emerging carbapenem-resistant enterobacteriaceae infection, its epidemiology and novel treatment options: a review. Infect Drug Resist. 2021;14:4363–74.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. World Health Organization (WHO). Prioritization of pathogens to guide discovery, research and development of new antibiotics for drug-resistant bacterial infections, including tuberculosis. WHO Bull. 2017.

  75. Chia PY, Sengupta S, Kukreja A, Ponnampalavanar S, Ng OT, Marimuthu K. The role of hospital environment in transmissions of multidrug-resistant gram-negative organisms. Antimicrob Resist Infect Control. 2020;9:29.

    Article  PubMed  PubMed Central  Google Scholar 

  76. Islam MA, Islam M, Hasan R, Hossain MI, Nabi A, Rahman M, et al. Environmental spread of New Delhi metallo-β-lactamase-1-producing multidrug-resistant bacteria in Dhaka. Bangladesh. Appl Environ Microbiol. 2017;83:e00793.

    Article  PubMed  PubMed Central  Google Scholar 

  77. Cohen R, Paikin S, Rokney A, Rubin-Blum M, Astrahan P. Multidrug-resistant enterobacteriaceae in coastal water: an emerging threat. Antimicrob Resist Infect Control. 2020;9:169.

    Article  PubMed  PubMed Central  Google Scholar 

  78. Mahon BM, Brehony C, Cahill N, McGrath E, O’Connor L, Varley A, et al. Detection of OXA-48-like-producing Enterobacterales in Irish recreational water. Sci Total Environ. 2019;690:1–6.

    Article  CAS  PubMed  Google Scholar 

  79. Adegoke AA, Fatunla OK, Okoh AI. Critical threat associated with carbapenem-resistant gram-negative bacteria: prioritizing water matrices in addressing total antibiotic resistance. Ann Microbiol. 2020;70:1–13.

    Article  Google Scholar 

  80. Mills MC, Lee J. The threat of carbapenem-resistant bacteria in the environment: evidence of widespread contamination of reservoirs at a global scale. Environ Pollut. 2019;255:113143.

    Article  CAS  PubMed  Google Scholar 

  81. El Garch F, Bogaerts P, Bebrone C, Galleni M, Glupczynski Y. OXA-198, an acquired carbapenem-hydrolyzing class D β-lactamase from Pseudomonas aeruginosa. Antimicrob Agents Chemother. 2011;55:4828–33.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Mahfouz N, Ferreira I, Beisken S, von Haeseler A, Posch AE. Large-scale assessment of antimicrobial resistance marker databases for genetic phenotype prediction: a systematic review. J Antimicrob Chemother. 2020;75:3099–108.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. Avguštin JA, Etrič A, Ašić L. Screening the cultivable cave microbial mats for the production of antimicrobial compounds and antibiotic resistance. Int J Speleol. 2019;48:295–303.

    Article  Google Scholar 

  84. Zada S, Sajjad W, Rafiq M, Ali S, Hu Z, Wang H, et al. Cave microbes as a potential source of drugs development in the modern era. Microb Ecol. 2022;84:676–87.

    Article  CAS  PubMed  Google Scholar 

  85. Ghosh S, Kuisiene N, Cheeptham N. The cave microbiome as a source for drug discovery: reality or pipe dream? Biochem Pharmacol. 2017;134:18–34.

    Article  CAS  PubMed  Google Scholar 

  86. Benler S, Koonin EV. Recruitment of mobile genetic elements for diverse cellular functions in prokaryotes. Front Mol Biosci. 2022;9:1–15.

    Article  Google Scholar 

  87. Khedkar S, Smyshlyaev G, Letunic I, Maistrenko OM, Coelho LP, Orakov A, et al. Landscape of mobile genetic elements and their antibiotic resistance cargo in prokaryotic genomes. Nucleic Acids Res. 2022;50:3155–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Kim H, Kim M, Kim S, Lee YM, Shin SC. Characterization of antimicrobial resistance genes and virulence factor genes in an Arctic permafrost region revealed by metagenomics. Environ Pollut. 2022;294:118634.

    Article  CAS  PubMed  Google Scholar 

  89. Botelho J. Defense systems are pervasive across chromosomally integrated mobile genetic elements and are inversely correlated to virulence and antimicrobial resistance. Nucleic Acids Res. 2023;51:4385–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Dai D, Brown C, Bürgmann H, Larsson DGJ, Nambi I, Zhang T, et al. Long-read metagenomic sequencing reveals shifts in associations of antibiotic resistance genes with mobile genetic elements from sewage to activated sludge. Microbiome BioMed Central. 2022;10:1–16.

    Google Scholar 

  91. Han XM, Hu HW, Chen QL, Yang LY, Li HL, Zhu YG, et al. Antibiotic resistance genes and associated bacterial communities in agricultural soils amended with different sources of animal manures. Soil Biol Biochem Pergamon. 2018;126:91–102.

    Article  CAS  Google Scholar 

  92. Hwengwere K, Paramel Nair H, Hughes KA, Peck LS, Clark MS, Walker CA. Antimicrobial resistance in Antarctica: is it still a pristine environment? Microbiome BioMed Central. 2022;10:1–13.

    Google Scholar 

  93. Marcoleta AE, Arros P, Varas MA, Costa J, Rojas-Salgado J, Berríos-Pastén C, et al. The highly diverse Antarctic Peninsula soil microbiota as a source of novel resistance genes. Sci Total Environ. 2022;810:152003.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgments

T.V.Z and K.K. want to acknowledge the University Computing Centre SRCE, at University of Zagreb, for the computational resources offered by the Isabella computer cluster.

Funding

This work was partially supported under the project STIM–REI, Contract Number: KK.01.1.1.01.0003, a project funded by the European Union through the European Regional Development Fund—the Operational Programme Competitiveness and Cohesion 2014–2020 (KK.01.1.1.01). A.H. was supported by the 9th China-Croatia Science and Technology cooperation committee program (No. 9–21). R.S. was supported by the National Recovery and Resilience Plan (NRRP), Mission 4 Component 2 Investment 1.4—Call for tender No. 3138 of December 16, 2021, rectified by Decree n.3175 of December 18, 2021 of Italian Ministry of University and Research funded by the European Union—NextGenerationEU (code CN_00000033, Concession Decree No. 1034 of June 17, 2022 adopted by the Italian Ministry of University and Research, CUP B83C22002930006, Project title “National Biodiversity Future Center—NBFC”).

Author information

Author notes

  1. Tanja Vojvoda Zeljko and Katarina Kajan have equally contributed to this work.

Authors and Affiliations

  1. Division for Marine and Environmental Research, Ruđer Bošković Institute, Bijenička Cesta 54, 10000, Zagreb, Croatia

    Tanja Vojvoda Zeljko, Neven Cukrov, Marija Marguš & Nuša Cukrov

  2. Division of Materials Chemistry, Ruđer Bošković Institute, Bijenička Cesta 54, 10000, Zagreb, Croatia

    Katarina Kajan & Sandi Orlić

  3. Center of Excellence for Science and Technology-Integration of Mediterranean Region (STIM), Zagreb, Croatia

    Katarina Kajan & Sandi Orlić

  4. Croatian Biospeleological Society, 10000, Zagreb, Croatia

    Branko Jalžić

  5. CAS Key Laboratory of Urban Pollutant Conversion, Institute of Urban Environment, Chinese Academy of Sciences, Xiamen, 361021, China

    Anyi Hu

  6. Croatian Geological Survey, 10000, Zagreb, Croatia

    Tamara Marković

  7. Molecular Ecology Group (MEG), National Research Council of Italy (CNR), Water Research Institute (IRSA), Largo Tonolli 50, 28922, Verbania, Italy

    Raffaella Sabatino & Andrea Di Cesare

  8. National Biodiversity Future Center (NBFC), Piazza Marina 61, 90133, Palermo, Italy

    Raffaella Sabatino & Andrea Di Cesare

Authors
  1. Tanja Vojvoda Zeljko
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Katarina Kajan
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Branko Jalžić
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Anyi Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Neven Cukrov
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Marija Marguš
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Nuša Cukrov
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Tamara Marković
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Raffaella Sabatino
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Andrea Di Cesare
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Sandi Orlić
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceptualization: T.V.Z., K.K., S.O.; Field work and lab analysis: K.K., S.O., B.J., N.C., M.M., N.C., T.M.; Methodology: T.V.Z., K.K.; Data curation: T.V.Z., K.K., R.S., A.C.; Formal analysis: T.V.Z., K.K., R.S., A.C.; Writing-original draft: T.V.Z.; Writing—review & editing: T.V.Z., K.K., R.S., A.C., A.H., S.O., Funding acquisition: S.O., A.H.; Project administration: S.O. All the authors contributed to manuscript revision and approved the final manuscript.

Corresponding author

Correspondence to Sandi Orlić.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. 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://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Vojvoda Zeljko, T., Kajan, K., Jalžić, B. et al. Genome-centric metagenomes unveiling the hidden resistome in an anchialine cave. Environmental Microbiome 19, 67 (2024). https://doi.org/10.1186/s40793-024-00612-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40793-024-00612-2

Keywords

Environmental Microbiome

ISSN: 2524-6372