Antimicrobial resistance (AMR) is a universal phenomenon the origins of which lay in natural ecological interactions such as competition within niches, within and between micro- to higher-order organisms. To study these phenomena, it is crucial to examine the origins of AMR in pristine environments, i.e., limited anthropogenic influences. In this context, epilithic biofilms residing in glacier-fed streams (GFSs) are an excellent model system to study diverse, intra- and inter-domain, ecological crosstalk. We assessed the resistomes of epilithic biofilms from GFSs across the Southern Alps (New Zealand) and the Caucasus (Russia) and observed that both bacteria and eukaryotes encoded twenty-nine distinct AMR categories. Of these, beta-lactam, aminoglycoside, and multidrug resistance were both abundant and taxonomically distributed in most of the bacterial and eukaryotic phyla. AMR-encoding phyla included Bacteroidota and Proteobacteria among the bacteria, alongside Ochrophyta (algae) among the eukaryotes. Additionally, biosynthetic gene clusters (BGCs) involved in the production of antibacterial compounds were identified across all phyla in the epilithic biofilms. Furthermore, we found that several bacterial genera (Flavobacterium, Polaromonas, Superphylum Patescibacteria) encode both atimicrobial resistance genes (ARGs) and BGCs within close proximity of each other, demonstrating their capacity to simultaneously influence and compete within the microbial community. Our findings help unravel how naturally occurring BGCs and AMR contribute to the epilithic biofilms mode of life in GFSs. Additionally, we report that eukaryotes may serve as AMR reservoirs owing to their potential for encoding ARGs. Importantly, these observations may be generalizable and potentially extended to other environments that may be more or less impacted by human activity.
IMPORTANCE Antimicrobial resistance is an omnipresent phenomenon in the anthropogenically influenced ecosystems. However, its role in shaping microbial community dynamics in pristine environments is relatively unknown. Using metagenomics, we report the presence of antimicrobial resistance genes and their associated pathways in epilithic biofilms within glacier-fed streams. Importantly, we observe biosynthetic gene clusters associated with antimicrobial resistance in both pro- and eukaryotes in these biofilms. Understanding the role of resistance in the context of this pristine environment and complex biodiversity may shed light on previously uncharacterized mechanisms of cross-domain interactions.


Today, antimicrobial resistance (AMR) has become a well-known threat to human health with an estimated number of 700,000 people per year dying of drug-resistant infections (1). The dramatic rise of antimicrobial resistance over the past decade has even led to the moniker “silent pandemic” (2). Therefore, AMR is often directly associated with human impacted environments with a global increase in resistant bacteria linked to the over- and misuse of antibiotics (3). However, contrary to public perception, AMR is a natural phenomenon, which has existed for billions of years (4). Long before the rather recent use of antibiotics in the clinical setting, microorganisms have used these, along with corresponding protective mechanisms, to establish competitive advantages over other microbes contending for the same environment and/or resources (5).
Microorganisms, in general, produce a range of secondary metabolites with diverse chemical structures, which in turn confer a variety of functions, including antibiotics (6). Such secondary metabolites, including metal transporters and quorum-sensing molecules (7, 8) are associated with enhanced fitness by driving growth, including by inhibiting the growth of other bacteria by depleting iron through increased capacity to compete for this essential resource. Consequently, many of these natural products have found their applications in industrial settings, but also as anti-infective drugs in human medicine (7, 9, 10). The biosynthetic pathways responsible for producing these specialized metabolites are encoded by locally clustered groups of genes known as “biosynthetic gene clusters” (BGCs). Typically, BGCs include genes for expression control, self-resistance, and metabolite export (11). They can, however, be further divided into various classes, including nonribosomal peptide synthetases (NRPSs), type I and type II polyketide synthases (PKSs), terpenes, and bacteriocins, alongside others (10). NRPSs and PKSs specifically have been of interest due to their known synthesis of putative antibiotics (12, 13). Furthermore, evidence suggests that within these BGCs at least one resistance gene conferring resistance can be found as a self-defense mechanism against the potentially harmful secondary metabolites encoded by the BGC (14). For instance, the tylosin-biosynthetic gene cluster of Streptomyces fradiae also encodes three resistance genes (tlrB, tlrC, and tlrD) (15), while in another report, Streptomyces toyacaensis, the vanHAX resistance cassette, is proximal to the vancomycin biosynthesis gene cluster, thereby encoding inherent resistance (16).
Remote and pristine microbial communities provide a potentially rich genetic resource to explore naturally occurring antibiotic resistance from the preantibiotic era. These phenomena may currently only be evident in a few environments with limited anthropogenic influence (e.g., permafrost, glaciers, deep sea, and polar regions). These ARGs and resistant bacteria evolving in pristine environments may therefore be considered the inherent antibiotic resistance present in the environment (5). Despite this, some recent reports suggest that even these remote microbial communities may be coming under the influence of anthropogenic antimicrobial resistance (17), thus necessitating the need to characterize these relatively pristine environments before their imminent decline due to climate change.
We have recently reported the genomic and metabolic adaptations of epilithic biofilms to environmental conditions in glacier-fed streams (GFSs) (18). For example, given the short flow season during glacial melt in summer, the incentive to reproduce quickly while conditions are favorable, is high. During these so-called windows of opportunity, the necessity for taxa to appropriate resources, especially when they are limiting as in GFSs, yields a competitive environment. Within these biofilms, we observed cross-domain interactions and internal recycling of nutrients between microorganisms to potentially mitigate the harsh nutrient and environmental conditions of the GFSs. Additionally, owing to complex biodiversity (19) and generally oligotrophic conditions (20) in GFSs, epilithic biofilms are ideal model systems for understanding BGCs and AMR in the context of cooccurring microorganisms. Our previous insights revealed that taxa such as Polaromonas, Acidobacteria, and Methylotenera have strong interactions with eukaryotes such as algae and fungi (21). The inherent diversity allows for understanding the influence of AMR in microbial interactions as evidenced by bacterial–fungal interactions leading to the discovery of penicillin by Alexander Fleming in 1928 (22). Along similar lines, Netzker et al. (23) reported that microbial interactions lead to the production of bioactive compounds, including antibiotics potentially shaping the microbial consortia within a community.
Here, to shed light on the role of AMR in shaping microbial communities within pristine environments, we used whole-genome shotgun sequencing (metagenomics) to investigate twenty-one epilithic biofilms from glacier-fed streams. These samples were collected from 8 GFSs spread across the Southern Alps in New Zealand and the Caucasus in Russia (Table S1). Herein, we found 29 categories of ARGs within the GFSs across both bacterial and eukaryotic domains. Importantly, most of the AMR was found in bacteria. We also identified antibacterial BGCs that were encoded both in bacterial and eukaryotic metagenome-assembled genomes (MAGs). Our findings demonstrate that microorganisms within biofilms from pristine environments not only encode ARGs, but that they may potentially influence several features of epilithic biofilms such as biofilm formation, community assembly, and/or maintenance, including conferring potential mechanisms for competitive advantages under extreme conditions.


Antimicrobial resistance in a pristine environment.

We characterized the resistomes of GFS epilithic biofilms and assessed the distribution of AMR in twenty-one epilithic biofilm samples, across 8 individual glaciers originating from the Southern Alps in New-Zealand (SA1, SA2, SA3, and SA4) and the Caucasus in Russia (CU1, CU2, CU3, CU4). In total, we identified a high number (n = 1,840) of ARGs within 29 categories of AMR, with similar AMR profiles observed across all GFSs (Fig. 1a), except for CU1, CU2, and SA3, where the differences were driven by elevated fluoroquinolone, peptide, glycopeptide, and phenicol resistance (Fig. S1). It is to be noted that while ARGs refer to the genes encoding specific resistance, AMR categories derived from metagenomic data in this context, typically reflect the functional potential associated with respect to the resistance encoded. Of the identified AMR categories, beta-lactam and multidrug resistance (i.e., resistance conferring protection against multiple antibiotic classes), followed by aminoglycoside resistance, were found to be highly abundant in all samples. We subsequently analyzed the diversity of ARGs within the various resistance categories and found beta-lactam resistance to represent the largest resistance category, contributing 930 unique ARGs to the resistome. This was followed by multidrug (179 ARGs) and aminoglycoside (176 ARGs) resistance (Table S2). In contrast, some resistance categories such as polymyxin and pleuromutilin resistance were only detected at very low levels within the epilithic biofilm resistomes.
FIG 1 Epilithic biofilms in GFSs harbor a diverse resistome. (a) Relative abundance of 29 AMR categories within 21 epilithic biofilms collected from four New Zealand Southern Alps (SA) and four Russian Caucasus (CU) GFSs. (b) Bar plots depicting the relative abundance of bacteria and eukaryotes encoding ARGs. (c) Phylum-level representation of the AMR abundances across bacteria and eukaryotes. Size of the closed circle indicates the normalized relative abundance (RNum_Gi; see Materials and Methods), and the color represents individual phyla.
We further investigated the contribution of microbial populations to the resistome and found contributions from both prokaryotes and eukaryotes (Fig. 1b). Prokaryotes within this study refer to bacteria alone, since archaea encoded an infinitesimal number of ARGs (<0.000001% RNum_GI; see Materials and Methods.), and therefore were excluded from further analyses. Among the eukaryotes, the phylum Ochrophyta (algae) was the dominant contributor and encoded most of the AMR categories (Fig. 1c; Fig. S2a). In bacteria, AMR was more evenly distributed, with most of the phyla encoding ARGs across all categories (Fig. 1c). However, members of the Alphaproteobacteria, Betaproteobacteria, and the Bacteroidetes/Chlorobi group encoded the highest overall ARG abundance (Fig. 1c; Fig. S2b). Additionally, AMR categories such as aminoglycoside, beta-lactam, glycopeptide, and rifamycin resistance (among others) were widely distributed in both bacteria as well as among the eukaryotes. On the other hand, categories such as aminocoumarin, bacitracin, and diaminopyrimidine resistance were found to be primarily encoded by bacteria.
Given the novelty of the identification of putative ARGs in eukaryotes, we applied additional checks to validate the ARG sequences. First, the eukaryotic ARG sequences were assessed for homology with known bacterial representatives from the CARD database. For example, both TETA (tetracycline) and ERM (MLS phenotype) sequences were the most representative for ARGs in ~95% of the eukaryotic MAGs and were therefore retrieved from the eukaryotic contigs. Additionally, these genes were chosen for assessment given the availability of different isoforms in CARD. Phylogenetic analyses subsequently revealed that several of the eukaryotic-encoded ARGs indeed shared considerable homology with their bacterial counterparts (Fig. 2), as identified by their grouping within clades. A similar phenomenon was observed when comparing the ARGs retrieved from the GFS bacteria with the CARD reference ARGs (Fig. S3). Interestingly, several eukaryotic ARGs were distantly related, raising the hypothesis that they may be naturally occurring enzymes that may be distantly related homologs of the reference ARGs. To account for this possibility, all the ARG sequences from the eukaryotes were compared against the UniProtKB database (Fig. S4a and Table S3). With both TETA and ERM, we found that a majority (>73%) of the prokaryotic ARGs had hits in the database with at least 50% identity covering 90% of the total sequence length. Only ~8% of the eukaryotic ARGs demonstrated 50% identity with the database sequences, albeit with high E values suggesting poor or partial homology with existing reference sequences. We manually verified the partially homologous sequences from eukaryotic TETA resistance genes and found weak homology (<58%) to the existing list of entries in UniProtKB (representative example: Fig. S4b). Among the amino acids with substantial homology in terms of their positions, we found that 141 to 159 aa from the TETA genes found in UniProtKB were partially homologous with those of the eukaryotes’ TETA genes. Nevertheless, it is at present unclear whether these amino acids form the active sites conferring tetracycline resistance. Due to the lack of functional information, we additionally used fARGene (24) to determine the likelihood that functional resistance genes were encoded by these genes and to subsequently classify any novel antibiotic resistance genes. Although fARGene classified and identified all bacterial ARGs, eukaryote-encoded ARGs were unclassified, likely due to their novelty or the lack of appropriate representatives in the database underlying fARGene.
FIG 2 Phylogenetic analyses of putative eukaryotic ARGs. (a) Phylogenetic tree depicting the relatedness and potential homology between ARG sequences from the reference database, i.e., CARD, with the TETA (tetracycline) resistance gene sequence retrieved from GFS eukaryotes. (b) Tree built using the ERM (MLS phenotype) resistance genes retrieved from the eukaryotes and CARD databases. Group in both figures indicates the origins of the sequence, i.e., CARD database or eukaryote.

Antibiotic biosynthesis pathways and biosynthetic gene clusters.

As described above, beta-lactam, multidrug, and aminoglycoside resistance were the most abundant resistance categories within GFS epilithic biofilms. This was not surprising as beta-lactams and aminoglycosides are natural and prevalent compounds (25, 26). Furthermore, multidrug resistance is typically conferred via efflux machineries which were also common in the GFS epilithic biofilms. These typically serve dual purposes in particular for protein export within most bacteria (27). Based on these results, it is therefore likely that pristine environments such as GFSs potentially reflect the spectrum of natural antibiotics and their resistance mechanisms, reinforcing their capacity to serve as natural baselines for assessing enrichments and spread of AMR.
To further understand if these encoded resistance genes reflected natural antibiotic pressure, we investigated pathways associated with antibiotic biosynthesis using the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (28). In total, we identified seven different pathways corresponding to the biosynthesis of macrolides (MLS), ansamycins, glycopeptides (vancomycin), beta-lactams (monobactam, penicillin, and cephalosporin), aminoglycosides (streptomycin), and tetracyclines, which were present in various abundances in all samples (Fig. S5a). Importantly, the identified antibiotic synthesis genes thereby corresponded to the resistance categories identified within the epilithic biofilms. Interestingly, in most of the GFSs, antibiotic biosynthesis was primarily encoded by bacteria spanning multiple phyla (Fig. S5b, c). Exceptions to these were GL11 and GL15 in which biosynthesis pathways were equally distributed among eukaryotes, specifically Ochrophyta, in addition to bacteria.
KEGG pathways described for the production of certain antibiotics tend to contain enzymes involved in far-upstream reactions. For example, the “streptomycin biosynthesis” pathway in KEGG ortholog (ko00521) is defined with 21 KOs, including general and widespread genes such as glucose metabolism. Therefore, to further validate our observations regarding the antibiotic biosynthesis pathways, we assessed the abundance of BGCs, in the epilithic bioficlms within the MAGs, that are known to encode genes for secondary metabolite synthesis, including antibiotics. We found six different structural classes of BGCs by annotating 537 medium-to-high-quality (>50% completion and <10% contamination) bacterial and 30 eukaryotic MAGs using antiSmash (29) and DeepBGC (30). Using this ensemble approach, we identified one or more BGCs in most bacterial (n = 490, ~91% of all bacterial MAGs) and eukaryotic (n = 28) MAGs. Of these BGCs, those annotated with an antibacterial function were dominant across the microbial populations, represented here by the MAGs, and were found across all phyla (Fig. 3a). Overall, a wider variety of BGCs associated with cytotoxic activity, inhibitory, and antifungal mechanisms were also identified in bacteria. Eukaryotes, on the other hand, encoded a high prevalence of antibacterial BGCs (~93% of all eukaryotic MAGs) (Fig. 3a). We further annotated those BGCs identified as antibacterial to determine their subtypes and found that most of them were “unknown” (Fig. 3b). However, other identified subtypes include ribosomally synthesized and posttranslationally modified peptides (RiPPs) such as bacteriocins, along with NRPs, PKs, and terpenes. In line with the abundance of aminoglycoside and beta-lactamase resistance genes, we found that ~19 of the MAGs encoded beta-lactamase gene clusters, while ~15% encoded aminoglycoside-associated BGCs, followed by BGCs of “unknown” classification.
FIG 3 Biosynthetic gene clusters indicate the resistome potential. (a) Heatmap depicting the overall abundance of BGCs identified across bacterial and eukaryotic MAGs. The respective phyla are listed on the left while the colored legend represents the taxonomic order. (b) In-depth characterization of the “antibacterial” BGCs found within all phyla and orders across medium-to-high quality MAGs. (c) Alluvial plots depicting the taxa where both BGCs and AMR were found adjacently on the same contig. Colours indicate the genera associated with the MAGs. (d) Genome plots indicating the location of the ARGs in relation to the biosynthetic gene clusters.
According to the resistance hypothesis (14), within or close to each BGC there is at least one gene conferring resistance to its encoded secondary metabolite. To test this, we assessed whether the MAGs encoding a BGC also encoded corresponding ARGs. In line with this hypothesis, we identified BGCs and their respective resistance genes in close proximity to each other through their localization on the same contig. Consequently, we identified various BGCs encoded together with ARGs in both the bacterial and eukaryotic MAGs. For example, we found that an antibacterial BGC was encoded by Flavobacterium spp. on the same contig as both MLS (i.e., macrolides, lincosamides, and streptogramin) and beta-lactam resistance genes (Fig. 3c). Incidentally, we also found that a candidate phyla radiation (CPR) bacterium (Aalborg-AAW-1; phylum Patescibacteria) also encoded both antibacterial BGC and MLS resistance on the same contig. Further closer inspection revealed that the ARGs were encoded within the BGCs in certain cases, while in other instances the ARGs directly flanked the BGC itself (Fig. 3d). Importantly, we found that the closely encoded ARGs and BGCs were localized within a distance of 100 bp to 2 Kb of each other.


Microbial reservoirs in pristine environments, with little to no impact from anthropogenic selection pressures, provide the opportunity to investigate the natural propensity and linked evolutionary origins of AMR. Here, by leveraging metagenomics on twenty-one epilithic biofilms, we assessed the resistomes of eight individual GFS epilithic biofilms.
To date, while many studies have looked for novel antibiotics and resistance genes in pristine environments such as the deep sea (31) or the polar regions (32), few have explored the full diversity of antibiotic resistance in such environments (33, 34). Van Goethem et al. (35) identified 117 naturally occurring ARGs associated with multidrug, aminoglycosidem and beta-lactam resistance in pristine Antarctic soils (17). Similarly, D’Costa et al. (4) identified a collection of ARGs encoding resistance to beta-lactams as well as tetracyclines and glycopeptides in 30,000-year-old Beringian permafrost sediments. In agreement with these previous studies, we identified 29 AMR categories, including the previously mentioned resistance categories, in the studied biofilm communities. Among these, the highest ARG abundance was associated with aminoglycoside and beta-lactam resistance. Our study further suggests that although the overall abundance differs, the epilithic resistome was highly similar in all GFSs, independent of origin (i.e., New Zealand or Russia). Furthermore, our results agree with the results obtained in other resistomes identified in pristine environments such as Antarctic soils and permafrost in terms of the identified ARGs. Unlike previous studies, where ARGs were primarily associated with bacteria, we report for the first time that AMR was associated with both bacteria and eukaryotes in various abundances in environmental samples, including GFSs. A previous study by Brown et al. (36) reported that the IRS-HR (isoleucyl-tRNA synthetase—high resistance) type gene conferring resistance against mupirocin was identified in Staphylococcus aureus. More importantly, they suggested that horizontal gene transfer led to the acquisition of IRS-HR genes by bacteria from eukaryotes (36). Despite these early reports, the contribution of eukaryotes to most resistomes, including from pristine environments, has largely been unexplored thus far. An exception to this was the report by Fairlamb et al. (37), who identified eukaryotic drug resistance, especially encoded by fungi (Candida and Aspergillus) and parasites (Plasmodium and Trypanosoma). However, most of these modes of resistance were highly specific toward particular drug treatments (37). Our results specifically revealed that taxa from the phylum Ochrophyta encoded resistance to 28 AMR categories, and this was also reflected in other (micro-)eukaryotes. Interestingly, phylogenetic analyses revealed shared homology between the ARGs from the GFSs compared to reference databases. There were also other ARGs that formed separate clades indicating that they may either be distantly related or previously unidentified resistance genes. While going beyond the scope of the present study, the distantly related ARGs would need to be experimentally validated as performed in previous studies on phage-borne ARGs (38) and novel ARGs (24). Importantly, dinoflagellates are known to produce biologically active natural compounds that have been implicated in ecological roles and thus have been targeted in bioprospecting drives (39). However, although the anti-inflammatory (40) and vasoconstrictive (41) activities of dinoflagellate-derived biomolecules have been well-defined, dinoflagellate-encoded BGCs have been largely unknown. It must be noted that our analyses provide insights on the functional potential of the identified ARGs, within both pro- and eukaryotes. Substantial additional wet-lab experiments are required to validate the function of the identified ARGs, especially within the identified eukaryotic genomes.
Apart from encoded resistance mechanisms, microalgae such as Ochrophyta have been of interest as a source of (new) antimicrobial compounds (42, 43). Martins et al. reported that extracts from different microalgae may potentially serve not only as antimicrobial agents, but also as anticancer therapeutics. However, our present results suggest that these taxa may also serve as environmental reservoirs for AMR itself. It is however presently unclear whether this phenomenon confers advantages with respect to niche occupation and protection against bacterial infection as well as whether the eukaryotes are sensitive to the antibiotics produced by them. Our attempts at trying to elucidate putative or potential functions using tools such as fARGene and through databases such as UniProtKB were unsuccessful. Although we found partial homology of the eukaryote-encoded TETA genes to those found in UniProtKB (141 to 159 aa), it is unclear whether the homologous regions represent an active site and/or the expression of the ARG. This circumstance further highlights the lack of sufficient knowledge highlighting the need for additional wet-lab validations of our findings.
Studies delving into the origins of AMR have reported that fecal pollution may explain ARG abundances in anthropogenically impacted environments (44). This phenomenon was also observed by Antelo et al. (45) and others (46) who detected ARGs in soils in Antarctica, especially in proximity to scientific bases. A similar trend manifested in our analyses, where quinolone (synthetic antibiotic) resistance was found in the glacier streams from the Southern Alps. Interestingly, we observed that the quinolone resistance was higher in a glacier that was accessible to hikers along the Cascade Saddle, a famous trekking hot spot. This observation is along the lines of a recent report by Hwengwere et al. suggesting that the concentrations of anthropogenically derived ARGs are greater in regions closer to the research stations in Antarctica (17). Although it is plausible that some of the GFSs sampled in our study may indeed be under anthropogenic influence, albeit minimal, the majority of them are considered to be pristine environments. In this context, AMR is likely derived from natural antibiotics produced by microorganisms as a competitive advantage given the oligotrophic conditions. Microorganisms acquire resistance either as a protective measure against other microorganisms (47, 48) or as a self-defense mechanism to prevent inadvertent suicide by damaging metabolites (14). Accordingly, we found both antibiotic biosynthesis pathways and BGCs within the epilithic resistomes. We identified pathways for the biosynthesis of glycopeptides, beta-lactams, and aminoglycosides, among others, concurrent with the high abundance of ARGs against said antibiotics. Additionally, we identified BGCs with a predicted antibacterial function in both eukaryotes and bacteria. While a limited number of studies such, as Waschulin et al. (49) and Liao et al. (50), have shown BGCs in pristine environments, none of these studies have contextualized the cooccurrence of BGCs with AMR. Hence, we not only found that most of our MAGs contain BGCs, of which many have an antibacterial function, but also found all MAGs to encode multiple resistance genes. Additionally, we found several BGCs closely localized to ARGs on the same contig, thereby indicating an immediate self-defense mechanism against the encoded secondary metabolites. This agrees with the resistance hypothesis highlighted by Tran et al. stating that a gene conferring resistance to potentially harmful metabolites produced by the organism is to be found within the BGC-encoding operons (14). We also observed that the recently identified CPR bacteria (51) (in our case, phylum Patescibacteria) not only encoded AMR but also harbored genes associated with the production of molecules with antibacterial effects. Although Patescibacteria have been identified in oligotrophic environments (52, 53) with carbon and/or nutrient limitations similar to those observed for GFSs, it is plausible that their ability to survive with minimal biosynthetic and metabolic pathways may indeed depend on the expression of BGCs and AMR. At the time of this writing, a preprint by Maatouk et al. (54), described the presence of ARGs across publicly available CPR bacterial genomes. In addition, we report the identification of AMR within GFS-derived CPR genomes, likely as a means of competitive inhibition against other taxa. Alternatively, biofilms may also allow for collective resistance, tolerance, and exposure protection to antibacterial compounds (55). The AMR and BGCs encoded by most phyla may therefore affect cooperation and/or interactions associated with nutrient exchange, leading to the privatization of public goods (55). Such a phenomenon may be achieved due to the competition within taxa, both at the intra- and interspecies levels, via secretion of toxins (48) and occupying spatial niches (56, 57) thereafter. Furthermore, Stubbendieck and Straight previously highlighted the multifaceted effects of bacterial competition which include the potential taxation and subsequent increase in bacterial fitness (58). Thus, the in-situ competition within multispecies biofilms may allow for cross-phyla and cross-domain interactions while simultaneously increasing the overall fitness of the endogenous epilithic microbial community. Alternatively, these interactions or lack thereof may shape the overall community, including spatial organization (59), especially in energy-limited systems such as the GFSs.


Epilithic biofilms are an integral and key mode of survival in the extreme environment of glacier-fed streams. Herein, we report that these biofilms provide critical insights into the naturally occurring resistome. Our findings shed light on the survival mechanisms via the resistome, including potential ecological dimensions of ARGs and BGCs within the GFS microbial communities. Furthermore, we reveal the congruence of genes encoding both BGCs and AMR, in both bacteria and eukaryotes. More importantly, we highlight for the first time the comprehensive AMR profile of CPR bacteria and of (micro-)eukaryotes. Collectively, our results highlight underlying resistance mechanisms, including BGCs, employed in “biological warfare” in oligotrophic and challenging glacier-fed stream ecosystems.


Sampling and biomolecular extractions.

Eight GFSs were sampled in early to mid-2019 from the New Zealand Southern Alps and the Russian Caucasus, respectively, for a total of 21 epilithic biofilms (Table S1). The biofilm samples were collected from each stream reach, depending on the stream geomorphology, i.e., the presence of boulders, as described previously by Busi et al. (18). One to three biofilm samples were collected per reach (Table S1), taken using sterilized metal spatulas to scrape rocks, followed by their immediate transfer to cryovials. Samples were immediately flash-frozen in liquid nitrogen and stored at −80°C until DNA was extracted. DNA from the epilithic biofilms was extracted using a previously established protocol (60) adapted to a smaller scale due to relatively high DNA concentrations. DNA quantification was performed for all samples with the Qubit dsDNA HS kit (Invitrogen).

Sequencing and data processing for metagenomics.

Random shotgun sequencing was performed on all epilithic biofilm DNA samples after library preparation using the NEBNext Ultra II FS library kit. Fifty nanograms of DNA was enzymatically fragmented for 12.5 min, and libraries were prepared with six PCR amplification cycles. An average insert of 450 bp was maintained for all libraries. Qubit was used to quantify the libraries followed by sequencing at the Functional Genomics Centre Zurich on a NovaSeq (Illumina) using a S4 flowcell. The metagenomic data were processed using the Integrated Meta-omic Pipeline (IMP v3.0; commit# 9672c874 available at https://git-r3lab.uni.lu/IMP/imp3) (61). IMP’s workflow includes preprocessing, contig assembly, genome reconstruction (metagenome-assembled genomes, i.e., MAGs), and taxonomic and additional functional analysis of genes based on custom databases in a reproducible manner (61). Eukaryotic binning was performed in addition to the processing with the IMP workflow as described in Busi et al. (18). Briefly, eukaryotic MAGs were binned using CONCOCT (62), using contigs that were identified as eukaryotic, including those not mapping to bacteria or viruses in the maxikraken2 database available at https://lomanlab.github.io/mockcommunity/mc_databases.html. The eukaryotic taxonomies were further assigned using the PhyloDB and MMETSP databases associated with EUKulele (commit# fb8726a; available at https://github.com/AlexanderLabWHOI/EUKulele). Subsequently, contigs binned in eukaryotic MAGs were independently validated using WHOkaryote (63), which distinguishes eukaryotic and prokaryotic contigs in metagenomes based on gene structure. Consensus taxonomy per contig was then used for downstream analyses and association with ARGs.

Identification of antimicrobial resistance genes, antibiotic biosynthesis pathways, and BGCs.

For the prediction of ARGs, the IMP-generated contigs were used as input for PathoFact (64). Identified ARGs were further collapsed into their respective AMR categories in accordance with the Comprehensive Antibiotic Resistance Database (CARD) (65). PathoFact utilizes two distinct tools: DeepARG-LS, a deep-learning model, and RGI, based on homology and SNP models, to identify ARGs, therefore possibly also detecting resistance genes within eukaryotic genomic fragments. Subsequently, the raw read counts per open reading frame (ORF), obtained from PathoFact, were determined using FeatureCounts (66).
To identify pathways for the biosynthesis of antibiotics, we assigned KEGG orthology (KOs) identifiers to the ORFs using a hidden Markov model (HMM) (67) approach using hmmsearch from HMMER 3.1 (68) with a minimum bit score of 40. Additionally, we linked the identified KOs to their corresponding KEGG orthology pathways and extracted the pathways annotated as antibiotic biosynthesis pathways by KEGG. Both the identified ARGs and KEGG pathways were then further linked to associated bacterial taxonomies using the contig classification. We further identified BGCs within the MAGs using antiSMASH (ANTIbiotics & Secondary Metabolite Analysis SHell) (29) and annotated these using deepBGC (30). Owing to the possibility that a single tool such as deepBGC may not be ideal for identifying BGCs, we additionally used GECCO to (69) identify and annotate BGCs. Since GECCO was trained on the MiBiG database (70), we expected to improve our detection sensitivity with the additional approach. BGCs that were identified by both tools, including those unique to each, were used for downstream analyses. To link BGCs and ARGs, we linked the resistance genes to the assembled contigs, followed by identifying the corresponding bins (MAGs) to which said contigs belonged. We subsequently extracted the genomic coordinates from the individual GenBank (.gbk) files for each of the BGCs, based on the deepBGC and GECCO analyses. The genomic coordinates were used to establish the location of the BGCs with respect to the ARGs identified via PathoFact. For the analyses described above, we only considered ARGs and BGCs within a maximum distance of 2 genes from each other, which corresponded to up to 2.5 Kb. This was based on a previous report by Conway et al. (71), who reported that, on average, there are 1.98 genes in a bacterial operon.

Validation and phylogenetic analyses of putative ARGs.

To ensure the accuracy of the ARGs identified within both pro- and eukaryotes and assess their homology to sequences within existing databases, the respective ARG sequences were retrieved and the following analyses performed. The retrieved sequences were cross-referenced using DIAMOND (72) against the UniProtKB database (73) to determine if the putative ARGs found within pro- and eukaryotes encoded other functions. To further determine how the putative ARGs from GFSs relate to their homologs, the existing ARG database (i.e., CARD), a phylogenetic analysis was performed. For the phylogenetics analyses, all ARG nucleotide sequences were collected from the CARD database and grouped by classification. Pro- and eukaryotic ARGs from the GFS samples with similar classifications were collated respectively, into multifasta files along with the reference sequences, and a multiple sequence alignment was performed using MAFFT (74) with default settings. The aligned sequences were subsequently trimmed with trimAl (75) using “-gt 0.9” as a threshold for gap filtering. Lastly, IQ-TREE 2 (76) was used to construct a phylogenetic tree using trimmed multiple sequence alignment with at least 10 FASTA sequences. The built trees were visualized using the ggtree package (77) using version 3.6 of the R statistical software package (78).

Data analysis.

The relative abundance of the ORFs was calculated based on the RNum_Gi method described by Hu et al. (79). Figures for the study, including visualizations derived from the taxonomic and functional analyses, were created using the R statistical software and the tidyverse package (80). Alluvial plots were generated using the ggalluvial package (81), while heatmaps were generated using the ComplexHeatmap package (82) developed for R. The corresponding visualization and analysis code is available at https://github.com/susheelbhanu/rock_biofilm_amr.

Data availability.

The Biosample accession IDs listed under Table S4 can be found on NCBI under the BioProject accession number PRJNA781406. The analysis code for IMP and downstream analys is detailed at https://git-r3lab.uni.lu/susheel.busi/nomis_pipeline. Binning and manual refinement of eukaryotic MAGs was done as described here: https://git-r3lab.uni.lu/susheel.busi/nomis_pipeline/-/blob/master/workflow/notes/MiscEUKMAGs.md. All visualization and analysis code is available at https://git-r3lab.uni.lu/susheel.busi/rock_amr and https://github.com/susheelbhanu/rock_biofilm_amr.


We gratefully acknowledge the laboratory support from Emmy Marie Oppliger at EPFL and Lea Grandmougin, Janine Habier, and Laura Lebrun at the University of Luxembourg. We also acknowledge the key input from Rashi Halder at the LCSB Sequencing Platform regarding library preparation. We thank Patrick May and Anna Heintz-Buschart for their role in expanding the IMP workflow and also Cedric Christian Laczny for the crucial insights into metagenomic processing. The computational analyses were performed at the HPC facilities at the University of Luxembourg (https://hpc.uni.lu) (83).
We declare that we have no competing interests.
This research has been supported by The NOMIS Foundation to T.B. and the Swiss National Science Foundation (CRSII5_180241) supporting S.B.B. L.d.N. and P.W. are supported by the Luxembourg National Research Fund (FNR; PRIDE17/11823097) awarded to P.W.
Authors’ contributions: S.B.B., L.d.N., P.W., and T.J.B. conceived the project. P.P. extracted DNA, and S.B.B. and P.P. prepared the metagenomic libraries for sequencing. S.B.B. and L.d.N. conceptualized and performed the data analyses. S.B.B. and L.d.N. wrote the manuscript with P.W. and T.B., with significant input and editing from all coauthors.

Supplemental Material

File (spectrum.04069-22-s0001.pdf)
File (spectrum.04069-22-s0002.xlsx)
File (spectrum.04069-22-s0003.xlsx)
File (spectrum.04069-22-s0004.xlsx)
File (spectrum.04069-22-s0005.xlsx)
File (spectrum.04069-22-s0006.xlsx)
ASM does not own the copyrights to Supplemental Material that may be linked to, or accessed through, an article. The authors have granted ASM a non-exclusive, world-wide license to publish the Supplemental Material files. Please contact the corresponding author directly for reuse.


Stanton IC, Bethel A, Leonard AFC, Gaze WH, Garside R. 2020. What is the research evidence for antibiotic resistance exposure and transmission to humans from the environment? A systematic map protocol. Environ Evid 9:12.
Balasegaram M. 2021. Learning from COVID-19 to tackle antibiotic resistance. ACS Infect Dis 7:693–694.
Wright GD. 2007. The antibiotic resistome: the nexus of chemical and genetic diversity. Nat Rev Microbiol 5:175–186.
D'Costa VM, King CE, Kalan L, Morar M, Sung WWL, Schwarz C, Froese D, Zazula G, Calmels F, Debruyne R, Golding GB, Poinar HN, Wright GD. 2011. Antibiotic resistance is ancient. Nature 477:457–461.
Scott LC, Lee N, Aw TG. 2020. Antibiotic resistance in minimally human-impacted environments. Int J Environ Res Public Health 17:3939.
Tyc O, Song C, Dickschat JS, Vos M, Garbeva P. 2017. The ecological role of volatile and soluble secondary metabolites produced by soil bacteria. Trends Microbiol 25:280–292.
Chen R, Wong HL, Kindler GS, MacLeod FI, Benaud N, Ferrari BC, Burns BP. 2020. Discovery of an abundance of biosynthetic gene clusters in shark bay microbial mats. Front Microbiol 11:1950.
Demain AL, Fang A. 2000. The natural functions of secondary metabolites, p 1–39. In Fiechter A (ed), History of modern biotechnology I. Springer Berlin Heidelberg, Berlin, Germany.
Newman DJ, Cragg GM. 2016. Natural products as sources of new drugs from 1981 to 2014. J Nat Prod 79:629–661.
Medema MH, Kottmann R, Yilmaz P, Cummings M, Biggins JB, Blin K, de Bruijn I, Chooi YH, Claesen J, Coates RC, Cruz-Morales P, Duddela S, Düsterhus S, Edwards DJ, Fewer DP, Garg N, Geiger C, Gomez-Escribano JP, Greule A, Hadjithomas M, Haines AS, Helfrich EJN, Hillwig ML, Ishida K, Jones AC, Jones CS, Jungmann K, Kegler C, Kim HU, Kötter P, Krug D, Masschelein J, Melnik AV, Mantovani SM, Monroe EA, Moore M, Moss N, Nützmann H-W, Pan G, Pati A, Petras D, Reen FJ, Rosconi F, Rui Z, Tian Z, Tobias NJ, Tsunematsu Y, Wiemann P, Wyckoff E, Yan X, et al. 2015. Minimum information about a biosynthetic gene cluster. Nat Chem Biol 11:625–631.
Martinet L, Naômé A, Deflandre B, Maciejewska M, Tellatin D, Tenconi E, Smargiasso N, de Pauw E, van Wezel GP, Rigali S. 2019. A single biosynthetic gene cluster is responsible for the production of bagremycin antibiotics and ferroverdin iron chelators. mBio 10:e01230-19.
Martínez-Núñez MA, López y López VE. 2016. Nonribosomal peptides synthetases and their applications in industry. Sustainable Chemical Processes 4:13.
Ridley CP, Lee HY, Khosla C. 2008. Evolution of polyketide synthases in bacteria. Proc Natl Acad Sci USA 105:4595–4600.
Tran PN, Yen M-R, Chiang C-Y, Lin H-C, Chen P-Y. 2019. Detecting and prioritizing biosynthetic gene clusters for bioactive compounds in bacteria and fungi. Appl Microbiol Biotechnol 103:3277–3287.
Cundliffe E, Bate N, Butler A, Fish S, Gandecha A, Merson-Davies L. 2001. The tylosin-biosynthetic genes of Streptomyces fradiae. Antonie Van Leeuwenhoek 79:229–234.
Kwun MJ, Hong H-J. 2014. Genome sequence of Streptomyces toyocaensis NRRL 15009, producer of the glycopeptide antibiotic A47934. Genome Announc 2:e00749-14.
Hwengwere K, Paramel Nair H, Hughes KA, Peck LS, Clark MS, Walker CA. 2022. Antimicrobial resistance in Antarctica: is it still a pristine environment? Microbiome 10:71.
Busi SB, Bourquin M, Fodelianakis S, Michoud G, Kohler TJ, Peter H, Pramateftaki P, Styllas M, Tolosano M, De Staercke V, Schön M, de Nies L, Marasco R, Daffonchio D, Ezzat L, Wilmes P, Battin TJ. 2022. Genomic and metabolic adaptations of biofilms to ecological windows of opportunity in glacier-fed streams. Nat Commun 13:2168.
Battin TJ, Besemer K, Bengtsson MM, Romani AM, Packmann AI. 2016. The ecology and biogeochemistry of stream biofilms. Nat Rev Microbiol 14:251–263.
Battin TJ, Wille A, Sattler B, Psenner R. 2001. Phylogenetic and functional heterogeneity of sediment biofilms along environmental gradients in a glacial stream. Appl Environ Microbiol 67:799–807.
Busi SB, Bourquin M, Fodelianakis S, Michoud G, Kohler TJ, Peter H, Pramateftaki P, Styllas M, Tolosano M, De Staercke V, Schön M, de Nies L, Marasco R, Daffonchio D, Ezzat L, Wilmes P, Battin TJ. 2021. Genomic and metabolic adaptations of biofilms to ecological windows of opportunities in glacier-fed streams. bioRxiv.
Gaynes R. 2017. The discovery of penicillin—new insights after more than 75 years of clinical use. Emerg Infect Dis 23:849.
Netzker T, Flak M, Krespach MK, Stroe MC, Weber J, Schroeckh V, Brakhage AA. 2018. Microbial interactions trigger the production of antibiotics. Curr Opin Microbiol 45:117–123.
Berglund F, Österlund T, Boulund F, Marathe NP, Larsson DGJ, Kristiansson E. 2019. Identification and reconstruction of novel antibiotic resistance genes from metagenomes. Microbiome 7:52.
Krause KM, Serio AW, Kane TR, Connolly LE. 2016. Aminoglycosides: an overview. Cold Spring Harb Perspect Med 6:a027029.
Tahlan K, Jensen SE. 2013. Origins of the β-lactam rings in natural products. J Antibiot (Tokyo) 66:401–410.
Borges-Walmsley MI, McKeegan KS, Walmsley AR. 2003. Structure and function of efflux pumps that confer resistance to drugs. Biochem J 376:313–338.
Kanehisa M, Goto S. 2000. KEGG: Kyoto Encyclopedia of Genes and Genomes. Nucleic Acids Res 28:27–30.
Blin K, Shaw S, Kloosterman AM, Charlop-Powers Z, van Wezel GP, Medema MH, Weber T. 2021. antiSMASH 6.0: improving cluster detection and comparison capabilities. Nucleic Acids Res 49:W29–W35.
Hannigan GD, Prihoda D, Palicka A, Soukup J, Klempir O, Rampula L, Durcak J, Wurst M, Kotowski J, Chang D, Wang R, Piizzi G, Temesi G, Hazuda DJ, Woelk CH, Bitton DA. 2019. A deep learning genome-mining strategy for biosynthetic gene cluster prediction. Nucleic Acids Res 47:e110.
Tortorella E, Tedesco P, Palma EF, January GG, Fani R, Jaspars M, et al. 2018. Antibiotics from deep-sea microorganisms: current discoveries and perspectives. Mar Drugs 16:355.
McCann CM, Christgen B, Roberts JA, Su J-Q, Arnold KE, Gray ND, Zhu Y-G, Graham DW. 2019. Understanding drivers of antibiotic resistance genes in High Arctic soil ecosystems. Environ Int 125:497–504.
Yuan K, Yu K, Yang R, Zhang Q, Yang Y, Chen E, Lin L, Luan T, Chen W, Chen B. 2019. Metagenomic characterization of antibiotic resistance genes in Antarctic soils. Ecotoxicol Environ Saf 176:300–308.
Centurion VB, Delforno TP, Lacerda-Júnior GV, Duarte AWF, Silva LJ, Bellini GB, Rosa LH, Oliveira VM. 2019. Unveiling resistome profiles in the sediments of an Antarctic volcanic island. Environ Pollut 255:113240.
Van Goethem MW, Pierneef R, Bezuidt OKI, Van De Peer Y, Cowan DA, Makhalanyane TP. 2018. A reservoir of “historical” antibiotic resistance genes in remote pristine Antarctic soils. Microbiome 6:40.
Brown JR, Zhang J, Hodgson JE. 1998. A bacterial antibiotic resistance gene with eukaryotic origins. Curr Biol 8:S1–S12.
Fairlamb AH, Gow NAR, Matthews KR, Waters AP. 2016. Drug resistance in eukaryotic microorganisms. Nat Microbiol 1:16092.
Moon K, Jeon JH, Kang I, Park KS, Lee K, Cha C-J, Lee SH, Cho J-C. 2020. Freshwater viral metagenome reveals novel and functional phage-borne antibiotic resistance genes. Microbiome 8:75.
Sánchez-Suárez J, Garnica-Agudelo M, Villamil L, Díaz L, Coy-Barrera E. 2021. Bioactivity and biotechnological overview of naturally occurring compounds from the dinoflagellate family symbiodiniaceae: a systematic review. ScientificWorldJournal 2021:1983589.
Kita M, Ohishi N, Konishi K, Kondo M, Koyama T, Kitamura M, Yamada K, Uemura D. 2007. Symbiodinolide, a novel polyol macrolide that activates N-type Ca2+ channel, from the symbiotic marine dinoflagellate Symbiodinium sp. Tetrahedron 63:6241–6251.
Nakamura H, Asari T, Ohizumi Y, Kobayashi J, Yamasu T, Murai A. 1993. Isolation of zooxanthellatoxins, novel vasoconstrictive substances from the zooxanthella Symbiodinium sp. Toxicon 31:371–376.
Silva A, Silva SA, Carpena M, Garcia-Oliveira P, Gullón P, Barroso MF, Prieto MA, Simal-Gandara J. 2020. Macroalgae as a source of valuable antimicrobial compounds: extraction and applications. Antibiotics (Basel) 9:642.
Martins RM, Nedel F, Guimarães VBS, da Silva AF, Colepicolo P, de Pereira CMP, Lund RG. 2018. Macroalgae extracts from antarctica have antimicrobial and anticancer potential. Front Microbiol 9:412.
Karkman A, Pärnänen K, Larsson DGJ. 2019. Fecal pollution can explain antibiotic resistance gene abundances in anthropogenically impacted environments. Nat Commun 10:80.
Antelo V, Giménez M, Azziz G, Valdespino-Castillo P, Falcón LI, Ruberto LAM, Mac Cormack WP, Mazel D, Batista S. 2021. Metagenomic strategies identify diverse integron-integrase and antibiotic resistance genes in the Antarctic environment. Microbiologyopen 10:e1219.
Hernández F, Calısto-Ulloa N, Gómez-Fuentes C, Gómez M, Ferrer J, González-Rocha G, Bello-Toledo H, Botero-Coy AM, Boıx C, Ibáñez M, Montory M. 2019. Occurrence of antibiotics and bacterial resistance in wastewater and sea water from the Antarctic. J Hazard Mater 363:447–456.
Reygaert WC. 2018. An overview of the antimicrobial resistance mechanisms of bacteria. AIMS Microbiol 4:482–501.
Granato ET, Meiller-Legrand TA, Foster KR. 2019. The evolution and ecology of bacterial warfare. Curr Biol 29:R521–R537.
Waschulin V, Borsetto C, James R, Newsham KK, Donadio S, Corre C, Wellington E. 2021. Biosynthetic potential of uncultured Antarctic soil bacteria revealed through long-read metagenomic sequencing. ISME J.
Liao L, Su S, Zhao B, Fan C, Zhang J, Li H, et al. 2019. Biosynthetic potential of a novel Antarctic Actinobacterium Marisediminicola antarctica ZS314T revealed by genomic data mining and pigment characterization. Mar Drugs 17:388.
Hug LA, Baker BJ, Anantharaman K, Brown CT, Probst AJ, Castelle CJ, Butterfield CN, Hernsdorf AW, Amano Y, Ise K, Suzuki Y, Dudek N, Relman DA, Finstad KM, Amundson R, Thomas BC, Banfield JF. 2016. A new view of the tree of life. Nat Microbiol 1:16048.
Tian R, Ning D, He Z, Zhang P, Spencer SJ, Gao S, Shi W, Wu L, Zhang Y, Yang Y, Adams BG, Rocha AM, Detienne BL, Lowe KA, Joyner DC, Klingeman DM, Arkin AP, Fields MW, Hazen TC, Stahl DA, Alm EJ, Zhou J. 2020. Small and mighty: adaptation of superphylum Patescibacteria to groundwater environment drives their genome simplicity. Microbiome 8:51.
Vigneron A, Cruaud P, Langlois V, Lovejoy C, Culley AI, Vincent WF. 2020. Ultra-small and abundant: candidate phyla radiation bacteria are potential catalysts of carbon transformation in a thermokarst lake ecosystem. Limnol Oceanogr Lett 5:212–220.
Maatouk M, Ibrahim A, Rolain J-M, Merhej V, Bittar F. 2021. Small and equipped: the rich repertoire of antibiotic resistance genes in Candidate Phyla Radiation genomes. bioRxiv.
Bottery MJ, Pitchford JW, Friman V-P. 2021. Ecology and evolution of antimicrobial resistance in bacterial communities. ISME J 15:939–948.
Bottery MJ, Passaris I, Dytham C, Wood AJ, van der Woude MW. 2019. Spatial organization of expanding bacterial colonies is affected by contact-dependent growth inhibition. Curr Biol 29:3622–3634.e5.
Schluter J, Nadell CD, Bassler BL, Foster KR. 2015. Adhesion as a weapon in microbial competition. ISME J 9:139–149.
Stubbendieck RM, Straight PD. 2016. Multifaceted interfaces of bacterial competition. J Bacteriol 198:2145–2155.
Estrela S, Brown SP. 2018. Community interactions and spatial structure shape selection on antibiotic resistant lineages. PLoS Comput Biol 14:e1006179.
Busi SB, Pramateftaki P, Brandani J, Fodelianakis S, Peter H, Halder R, Wilmes P, Battin TJ. 2020. Optimised biomolecular extraction for metagenomic analysis of microbial biofilms from high-mountain streams. PeerJ 8:e9973.
Narayanasamy S, Jarosz Y, Muller EEL, Heintz-Buschart A, Herold M, Kaysen A, Laczny CC, Pinel N, May P, Wilmes P. 2016. IMP: a pipeline for reproducible reference-independent integrated metagenomic and metatranscriptomic analyses. Genome Biol 17:260.
Alneberg J, Bjarnason BS, de Bruijn I, Schirmer M, Quick J, Ijaz UZ. 2013. CONCOCT: Clustering cONtigs on COverage and ComposiTion. arXiv. http://arxiv.org/abs/1312.4038.
Pronk LJU, Medema MH. 2022. Whokaryote: distinguishing eukaryotic and prokaryotic contigs in metagenomes based on gene structure. Microb Genom 8:mgen000823.
de Nies L, Lopes S, Busi SB, Galata V, Heintz-Buschart A, Laczny CC, May P, Wilmes P. 2021. PathoFact: a pipeline for the prediction of virulence factors and antimicrobial resistance genes in metagenomic data. Microbiome 9:49.
Alcock BP, Raphenya AR, Lau TTY, Tsang KK, Bouchard M, Edalatmand A, Huynh W, Nguyen A-LV, Cheng AA, Liu S, Min SY, Miroshnichenko A, Tran H-K, Werfalli RE, Nasir JA, Oloni M, Speicher DJ, Florescu A, Singh B, Faltyn M, Hernandez-Koutoucheva A, Sharma AN, Bordeleau E, Pawlowski AC, Zubyk HL, Dooley D, Griffiths E, Maguire F, Winsor GL, Beiko RG, Brinkman FSL, Hsiao WWL, Domselaar GV, McArthur AG. 2020. CARD 2020: antibiotic resistome surveillance with the comprehensive antibiotic resistance database. Nucleic Acids Res 48:D517–D525.
Liao Y, Smyth GK, Shi W. 2013. featureCounts: An efficient general-purpose program for assigning sequence reads to genomic features. arXiv. http://arxiv.org/abs/1305.3347.
Yoon B-J. 2009. Hidden Markov models and their applications in biological sequence analysis. Curr Genomics 10:402–415.
Eddy SR. 2011. Accelerated profile HMM Searches. PLoS Comput Biol 7:e1002195.
Carroll LM, Larralde M, Fleck JS, Ponnudurai R, Milanese A, Cappio E, Zeller G. 2021. Accurate de novo identification of biosynthetic gene clusters with GECCO. bioRxiv.
Kautsar SA, Blin K, Shaw S, Navarro-Muñoz JC, Terlouw BR, van der Hooft JJJ, van Santen JA, Tracanna V, Suarez Duran HG, Pascal Andreu V, Selem-Mojica N, Alanjary M, Robinson SL, Lund G, Epstein SC, Sisto AC, Charkoudian LK, Collemare J, Linington RG, Weber T, Medema MH. 2020. MIBiG 2.0: a repository for biosynthetic gene clusters of known function. Nucleic Acids Res 48:D454–D458.
Conway T, Creecy JP, Maddox SM, Grissom JE, Conkle TL, Shadid TM, Teramoto J, San Miguel P, Shimada T, Ishihama A, Mori H, Wanner BL. 2014. Unprecedented high-resolution view of bacterial operon architecture revealed by RNA sequencing. mBio 5:e01442-14–e01414.
Buchfink B, Xie C, Huson DH. 2015. Fast and sensitive protein alignment using DIAMOND. Nat Methods 12:59–60.
UniProt Consortium. 2021. UniProt: the universal protein knowledgebase in 2021. Nucleic Acids Res 49:D480–D489.
Katoh K, Standley DM. 2013. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol 30:772–780.
Capella-Gutiérrez S, Silla-Martínez JM, Gabaldón T. 2009. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25:1972–1973.
Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, Lanfear R. 2020. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol 37:1530–1534.
Yu G, Smith DK, Zhu H, Guan Y, Lam TT-Y. 2017. Ggtree: an r package for visualization and annotation of phylogenetic trees with their covariates and other associated data. Methods Ecol Evol 8:28–36.
The R Development Core Team. 2013. R: a language and environment for statistical computing. R Core Team, Vienna, Auatria. https://www.yumpu.com/en/document/view/6853895/r-a-language-and-environment-for-statistical-computing.
Hu Y, Yang X, Qin J, Lu N, Cheng G, Wu N, Pan Y, Li J, Zhu L, Wang X, Meng Z, Zhao F, Liu D, Ma J, Qin N, Xiang C, Xiao Y, Li L, Yang H, Wang J, Yang R, Gao GF, Wang J, Zhu B. 2013. Metagenome-wide analysis of antibiotic resistance genes in a large cohort of human gut microbiota. Nat Commun 4:2151.
Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H. 2019. Welcome to the tidyverse. J Open Source Softw 4:1686.
Brunson J. 2020. ggalluvial: layered grammar for alluvial plots. J Open Source Softw 5:2017.
Gu Z, Eils R, Schlesner M. 2016. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics 32:2847–2849.
Varrette S, Bouvry P, Cartiaux H, Georgatos F. 2014. Management of an academic HPC cluster: the UL experience. International Conference on High Performance Computing Simulation (HPCS): 959–967. https://www.semanticscholar.org/paper/Management-of-an-academic-HPC-cluster%3A-The-UL-Varrette-Bouvry/333f6aa940113abc9ac7f779bef6bcef7eee60e2.

Information & Contributors


Published In

cover image Microbiology Spectrum
Microbiology Spectrum
Volume 11Number 114 February 2023
eLocator: e04069-22
Editor: Jeffrey A. Gralnick, University of Minnesota Twin Cities
PubMed: 36688698


Received: 8 October 2022
Accepted: 22 December 2022
Published online: 23 January 2023


  1. glacier-fed streams
  2. metagenomics
  3. antimicrobial resistance
  4. biosynthetic gene clusters
  5. cross-domain interactions
  6. freshwater ecosystems



Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
Laura de Nies
Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
Paraskevi Pramateftaki
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Massimo Bourquin
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Tyler J. Kohler
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Stilianos Fodelianakis
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Grégoire Michoud
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Michail Styllas
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Matteo Tolosano
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Vincent De Staercke
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Martina Schön
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland
Valentina Galata
Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
Tom Battin
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland


Jeffrey A. Gralnick
University of Minnesota Twin Cities


Susheel Bhanu Busi and Laura de Nies contributed equally to this work. Author order was determined by a game of “Sequence.”
Paul Wilmes and Susheel Bhanu Busi are the senior authors.
The authors declare no conflict of interest.

Metrics & Citations



  • For recently published articles, the TOTAL download count will appear as zero until a new month starts.
  • There is a 3- to 4-day delay in article usage, so article usage will not appear immediately after publication.
  • Citation counts come from the Crossref Cited by service.


If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. For an editable text file, please select Medlars format which will download as a .txt file. Simply select your manager software from the list below and click Download.

View Options

Figures and Media






Share the article link

Share with email

Email a colleague

Share on social media

American Society for Microbiology ("ASM") is committed to maintaining your confidence and trust with respect to the information we collect from you on websites owned and operated by ASM ("ASM Web Sites") and other sources. This Privacy Policy sets forth the information we collect about you, how we use this information and the choices you have about how we use such information.
FIND OUT MORE about the privacy policy