The host-specific resistome in environmental feces of Eurasian otters (Lutra lutra) and leopard cats (Prionailurus bengalensis) revealed by metagenomic sequencing.

Investigation of feces of wildlife, which is considered as reservoirs, melting pots, vectors and secondary sources of antimicrobial resistance genes (ARGs), provides insights into the risks and ecology of ARGs in the environment. Here, we investigated microbiomes, virulence factor genes (VFGs) of bacterial pathogens, and resistomes in environmental feces of Eurasian otters (Lutra lutra) and leopard cats (Prionailurus bengalensis) using shotgun metagenome sequencing. As expected, the taxonomic compositions of bacteria were significantly different between the animals. Importantly, we found that the compositions of ARGs were also significantly different between the animals. We detected ARGs including iri, tetA(P), tetB(P), floR, sulII, strA, strB, tetW and tetY. Some of them were significantly more abundant in either of the host animals, such as strA, strB and tetY in Eurasian otters, and tetA(P), tetW and iri in leopard cats. We also found that some ARGs were selectively correlated to particular VFGs-related bacteria, such as tetA(P) and tetB(P) to Clostridium, and iri to Mycobacterium. We also found that there were positive correlations between Acinetobacter and ARGs of multiple antimicrobial classes. The host-specific resistomes and VFGs-related bacteria may be due to differences in the host's gut microbiome, diet and/or habitat, but further investigation is needed. Overall, this study provided important baseline information about the resistomes of the wildlife in Korea, which may help the conservation of these endangered species and assessment of human health risks posed by ARGs and bacterial pathogens from wildlife.

Introduction

The emergence of antimicrobial resistance genes (ARGs) is a global public health concern, and the gut microbiome is a prime reservoir of ARGs [1]. The resident bacterial community in the intestine does not only exchange ARGs among themselves but also transfer to transient bacteria (passing through the colon). Similarly, virulence factor genes (VFGs) that enable microorganisms to cause disease can also spread among bacteria, and can convert non-pathogenic bacteria into potential pathogens. Recently, wildlife has attracted attention as reservoirs, melting pots, vectors, and secondary sources of ARGs [2]. Additionally, wildlife has long been known as a source of zoonotic pathogens [3]. Several studies have confirmed that wildlife could act as bio-sentinels of ARGs and VFGs and may therefore be useful in identifying potential sources of ARGs and VFGs that can spread into the environment [[4], [5], [6], [7]].

The gastrointestinal tracts of animals harbor diverse microbial communities that play key roles in many vital processes of host animal species, such as digestion, development, energy metabolism, immune response, and protection from pathogens. Despite the importance of the gut microbiome in the health and disease of host animals, our understanding is still limited about the animal gut microbiome. There are some studies, however, which have been performed on captive animals [8,9]. The fecal microbiome in wildlife has received less focus due to the challenges associated with sample collection and metadata gathering.

The diversity, structure, and function of the animal gut microbiome are mainly shaped by the host's genetic (e.g., evolutionary history) and ecological factors (e.g., diet, geography, and behavior) [10,11]. A growing number of animal phylogenetic studies showed that closely related animal species harbor more similar gut microbiome than distant relatives [12]. Thus, the intestinal microbiomes of animals are known to be host-specific. However, the host specificity of the intestinal resistome between animals has not been well known. We hypothesize that the intestinal resistomes are also host-specific due to differences in their diet and behavioral patterns and the gut microbiomes.

Here, we characterized the environmental fecal microbiome, including VFGs, and resistome of the Eurasian otter (Lutra lutra) and leopard cat (Prionailurus bengalensis) using metagenomic sequencing. In the IUCN Red List, the Eurasian otter is listed as ‘Near Threatened’ [13] and leopard cat as ‘Least Concern’ [14]. The metagenomic sequencing provides information that could be useful in the conservation of these endangered animals. Moreover, it will provide insights into the human health risks posed by ARGs and VFGs from wildlife. Understanding the risks posed by wildlife is important as rapid urbanization in Korea increases the proximity of humans to wildlife.

Materials and methods

Fecal samples

We analyzed 18 fecal samples of Eurasian otters (N = 7) and leopard cats (N = 11). The samples in this study have been studied in our previous studies for analyses of diet [15] and zoonotic parasites [16]. The location of fecal sampling is shown in Fig. 1. The samples' metadata are also listed in Supplementary Table 1. The samples were stored at −80 °C until analysis.

Fig. 1

Sampling sites of feces of Eurasian otters and leopard cats.

DNA extraction and sequencing

DNA was extracted from fecal samples using the DNeasy PowerSoil Pro-Kit (Qiagen, Hilden, Germany) with modifications suggested earlier [9,15]. Extracted DNA was sent to Integrated Microbiome Resource (Halifax, Nova Scotia, Canada) for shotgun metagenome sequencing. Metagenomic sequencing libraries were constructed using the Nextera Flex kit (Illumina, San Diego, CA), and sequenced on Illumina NextSeq system with 150-bp paired-end chemistry. Raw sequence reads were uploaded to the NCBI Sequence Read Archive database under the BioProject accession number PRJNA772888.

Taxonomic and functional annotations

The low-quality reads and adapter sequences were removed using Trimmomatic v0.39 with default settings [17]. Then, paired-end reads were merged using FLASH v1.2.11 [18]. DIAMOND v2.0.4.142 was used to align merged reads against NCBI non-redundant protein database (downloaded in January 2020) using BLASTX with the default setting [19]. The taxonomic binning of BLASTX results were performed using the naïve lowest common ancestor algorithm of MEGAN v6.5.7 [20] with default parameters. For functional annotation, BLASTX results were mapped to the Kyoto Encyclopedia of Genes and Genomes mapping files (July 2020) in MEGAN.

Additionally, 16S rRNA gene fragments were extracted using SortMeRNA v2.1b [21] and aligned against the SILVA database (release 138) in the mothur [22]. Sequences with read length < 100 bp and base-call ambiguity were removed. The taxonomic annotation of high-quality 16S rRNA gene sequences was performed against the EzBiocloud database (May 2018) [23] using the Naïve Bayesian Classifier with a confidence threshold of 80% in mothur.

ARGs and VFGs analysis

The raw metagenomic reads were processed with BBmap [24]: residual adapters, phiX contaminants, and low-quality ends were trimmed and removed with bbduk.sh set to the following parameters: k = 25, ktrim = r, mink = 11, trimq = 10. To quantify the relative abundance of the ARGs, we used the Short Read Sequence Typing program [25] to map the quality-filtered sequence reads and cluster similar sequences against the antibiotic resistance gene database (ARG-ANNOT Nt V3) containing all known ARGs sequences [26]. For the VFGs analysis, the clean sequence reads were aligned to the Virulence Factors of Bacterial Pathogens database (November 2020) [27] using Bowtie2 v2.4.1 [28]. The relative abundances of ARGs and VFGs were calculated using the reads per kilobase of genome equivalents (RPKG) metric, using the MicrobeCensus package [29].

Statistical analysis

R version 4.0.0 was used. Non-metric multidimensional scaling (NMDS) plot was used to visualize the pairwise Bray-Curtis dissimilarities of taxonomic and functional profiles. A permutational multivariate analysis of variance (PERMANOVA, ‘adonis’ function in vegan R package) was carried out to test the effect of the host animals on Bray-Curtis dissimilarities. We used Wilcoxon rank-sum tests to assess the effect of the hosts on the relative abundance of bacterial taxa, ARGs, and VFGs associated with pathogens. Furthermore, we have computed Spearman rank correlations between VFGs grouped according to taxonomy and ARGs.

Results

Sequencing statistics

The shotgun metagenome sequencing yielded ~776 million raw reads from 18 fecal samples with an average of ~43 million reads per sample with coverage ranging from ~25 to 57 million reads per sample (Supplementary Table 2). A total of ~277 million merged sequences were obtained with an average of ~15 million sequences per sample and with coverage ranging from ~9 to 20 million sequences per sample.

Taxonomic and functional compositions

The most dominant phylum, based on the extracted 16S rRNA gene sequences, across all fecal samples were Firmicutes (&65%) (Fig. 2a). The most abundant genus was Clostridium (~11%) (Fig. 2b). Some genera were differentially abundant between the hosts, of these, Psychrobacter, Marinobacter, Jeotgalibaca, Epulopiscium, and Jeotgalicoccus were abundant in Eurasian otters, and Lactobacillus, Pseudomonas, Peptostreptococcus, Blautia, Collinsella, and Clostridioides in leopard cats. In addition, in Supplementary Table 3, the relative abundance of the 50 most dominant bacterial genera is shown, which include Acinetobacter and Mycobacterium. NMDS ordinations based on Bray–Curtis dissimilarity matrices showed that both taxonomic (PERMANOVA: Pseudo-F = 5.5, R2 = 0.26, P = 0.001) and functional (PERMANOVA: Pseudo-F = 4.6, R2 = 0.22, P = 0.005) profiles were clustered according to the host animals (Fig. 3). The results of functional profiles are summarized in Supplementary material and Supplementary Fig. 1.

Fig. 2

(a) Stacked barplot showing the relative abundance of bacterial phyla in fecal samples of Eurasian otters and leopard cats based on extracted 16S rRNA gene sequences. (b) Heat-map showing the relative abundance of 20 most dominant bacterial genera in fecal samples of Eurasian otter and leopard cat. Asterisks on taxon names denote significant differences between hosts (*P < 0.05, **P < 0.01 and ***P < 0.001; Wilcoxon rank-sum). The genera in red or blue letters represent those that were significantly more abundant in the leopard cat or Eurasian otter, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3

Non-metric multidimensional scaling (NMDS) plots based on Bray–Curtis dissimilarities of (a) taxonomic and (b) functional genes profiles. The results are based on total sequence reads obtained by metagenomic sequencing, not based on extracted 16S rRNA gene sequences.

Virulence factor genes (VFGs)

The relative abundance of sequences related to VFGs grouped based on the genus-level taxonomy varied greatly among the samples, with Clostridium being the most dominant (Fig. 4). The relative abundance of sequences related to the selected VFGs based on RPKG is listed in Supplementary Table 4. The abundant VFGs associated with Clostridium include CPE2281, pilD, pilB2, pilC2, and pilB (Supplementary Table 4). These VFGs are known to be associated with Clostridium perfringens [27]. The VFGs that are known to be associated with Acinetobacter baumannii, an opportunistic human pathogen, such as bfmR, ABK1_0097, ABZJ_00085, and ABZJ_00086 [27] were also detected (Supplementary Table 4).

Fig. 4

(a) Chord diagram showing the distribution of the relative abundance based on the reads per kilobase of genome equivalents (RPKG) of virulence factor genes (VFGs) grouped according to taxonomy (at genus level) in fecal samples of Eurasian otters and leopard cats.

Antimicrobial resistance genes (ARGs)

A total of 14 ARG classes were identified (Fig. 5a). The relative abundance of the selected ARGs is shown in Supplementary Table 5. The aminoglycosides resistance genes were most abundant across all samples. The abundant ARGs in each class include strA, strB and aph(3′)-Ic for the aminoglycoside (AGly) class, sulII for the sulfonamide (Sul) class, floR for the phenicols (Phe) class, iri for the rifampin (Rif) class, and tetA(P), tetB(P), tetW and tetY for the tetracyclines (Tet) class (Fig. 5b). The ARGs of the AGly class tend to be more abundant in Eurasian otters, with statistical significance observed in strA, strB, aph(6)-Id, sat-2A, aphA-2, catA8, aph-stph and vanRc (Fig. 5b). Conversely, the ARGs of the Tet class tend to be more abundant in leopard cats, with statistical significance observed in tetA(P), tetW, tet(40), and tet(44). However, tetY was significantly more abundant in Eurasian otters (Fig. 5b).

Fig. 5

(a) Chord diagram showing the distribution of the relative abundance based on the reads per kilobase of genome equivalents (RPKG) of antimicrobial resistance genes (ARGs) classes in fecal samples of Eurasian otters and leopard cats. (b) Heatmap showing the relative abundance based on RPKG of the selected 30 ARGs. The ARGs shown are selected from the 100 most abundant ARGs. Among them, 30 ARGs that are differentially abundant between the hosts and/or the most abundant are selected. Asterisk (*) represents statistically significantly different (P < 0.05; Wilcoxon rank-sum) between the hosts. The genes in red or blue letters represent those that were significantly more abundant in the leopard cat or Eurasian otter, respectively. Abbreviations: AGly, aminoglycosides; Bla, beta-lactamases; Gly, glycopeptides; Fcd, fusidic acid; MLS, macrolide-lincosamide-streptogramin; Fcyn, fosfomycin; Phe, phenicols; Flq, fluoroquinolones; Rif, rifampin; Sul, sulfonamides; Oxzin oxazolidinones; Tet, tetracyclines; and Tmt, trimethoprim. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Correlation between ARGs and VFGs

Fig. 6 shows the Spearman's rank correlation coefficients (ρ) of the relative abundance between VFGs grouped according to the genus-level taxonomy and ARGs. The correlations were found between ARGs and the VFGs-related bacteria in some combinations. For instance, ARGs positively correlated with Clostridium include tetA(P) (ρ = 0.74) and tetB(P) (ρ = 0.76) (Fig. 6 and Supplementary Table 6). The ARGs that were positively correlated with Mycobacterium include iri (ρ = 0.81). The ARGs that were positively correlated with Acinetobacter include strA (ρ = 0.75), strB (ρ = 0.83), aph(3′)-Ic (ρ = 0.79), sulII (ρ = 0.75), tetY (ρ = 0.63), msrE (ρ = 0.53), vanRc (ρ = 0.52), CARB-8_OU (ρ = 0.56), CARB-5 (ρ = 0.59) and CARB-10 (ρ = 0.59).

Fig. 6

Spearman's rank correlation coefficients of the relative abundance based on the reads per kilobase of genome equivalents (RPKG) between virulence factor genes (VFGs) grouped according to taxonomy (at genus level) and antimicrobial resistance genes (ARGs). The colour represents the correlation coefficient, and the area of each circle represents the absolute value of the correlation coefficient. The values of the correlation coefficient and statistical results are shown in Supplementary Table 6.

Discussion

Similar to a previous study that investigated ARGs in feces of güiña (Leopardus guigna), a feline species, in Chile [7], our study revealed that ARGs were present in the feces of the leopard cat (Prionailurus bengalensis) in Korea. Moreover, we found that ARGs were present in the feces of the Eurasian otter (Lutra lutra), a semi-aquatic musteline species. The leopard cats' scats were collected in inland areas, while the otters' spraints were collected in a coastal area (Fig. 1). These suggest that ARGs are widespread regardless of the type of wildlife habitat. However, habitat type may have some influence on differential distribution of certain ARGs discussed below.

The most abundant ARG class was aminoglycosides (AGly) (Fig. 5a). The AGly-class ARGs have been reported to be abundantly detected from feces of chickens, pigs, and humans in China [30]. The most common ARGs in this class were strA and strB (Supplementary Table 5). These two genes are linked and encode enzymes that inactivate aminoglycosides including streptomycin [31]. In the sulfonamides (Sul) class, sulII was most abundant (Supplementary Table 5). The co-occurrence of strA-strB and sulII is reasonable because they are often linked on the bacterial genomes [31]. Several studies reported the detection of these genes from domesticated animals [[32], [33], [34], [35]]. Antimicrobials such as sulfonamide and spectinomycin, an aminoglycoside antimicrobial, are frequently used in pig farms to prevent and treat bacterial infections [35]. It is speculated that the use of antimicrobials and the discharge of resistant bacteria and contaminated water from animal farms may be one cause of the spread of ARGs in wildlife. Additionally, natural factors may also contribute since some bacteria are known to naturally produce aminoglycosides [36].

The AGly-class ARGs tended to be more abundant in the Eurasian otter (Fig. 5a). Therefore, there may be some factors that contribute to the spread of AGly resistance genes in this semi-aquatic species. Interestingly, it is reported that a plasmid-encoded AGly resistance gene (aph(3′)-Ic) was identified in the fish pathogen Pasteurella piscicida [37]. Our study also detected a large amount of aph(3′)-Ic particularly from the feces of Eurasian otters (Fig. 5b and Supplementary Table 5). Therefore, it may be possible to hypothesize that the AGly resistance genes derived from fish-related bacteria (e.g., fish pathogens) may also contribute to the prevalence of AGly resistance genes particularly in the Eurasian otter that is known to feed mainly on fishes [15].

The strA, strB, and aph(3′)-Ic genes in the AGly class and the sulII gene in the Sul class were positively correlated with Acinetobacter (Fig. 6 and Supplementary Table 6). Moreover, several ARGs in other classes such as tetY (Tet), CARB-8_OU, CARB-5 and CARB-10 (Bla), msrE (MLS), and vanRc (Gly) were also positively correlated with Acinetobacter (Fig. 6 and Supplementary Table 6). The VFGs (e.g., bfmR) related to Acinetobacter baumannii [27], an opportunistic human pathogen [38], were also detected (Supplementary Table 4). The co-occurrence of ARGs in multiple antimicrobial classes and A. baumannii-related VFGs indicates the possible presence of multidrug-resistant (MDR) A. baumannii. MDR A. baumannii in particular poses a serious global public health threat [38]. MDR A. baumannii used to be prevalent only in hospitals, and rarely isolated from the community [39]. However, our result indicates that it may have now spread to wildlife populations, too.

The second most abundant ARG class was tetracyclines (Tet) (Fig. 5a). Tetracyclines were consumed most widely in Korea during 2003–2012 [40] as therapeutic drugs for humans and animals, and growth promoters for livestock, poultry, and aquaculture [35]. Of the Tet-class ARGs, tetA(P), and tetB(P) were the most abundant (Supplementary Table 5). The tetA(P) and tetB(P) genes are found in Clostridium spp. especially Clostridium perfringens [41], a zoonotic pathogen that causes infections in poultry, pigs, cows, sheep, and horses [42]. The detection of tetA(P) and tetB(P) indicates the presence of Clostridium spp. including C. perfringens. In fact, we found that the abundance of Clostridium was correlated with the abundance of VFGs associated with C. perfringens, with ρ = 0.74 and 0.76 for tetA(P) and tetB(P), respectively (Fig. 6 and Supplementary Table 6). The resistance of C. perfringens to tetracyclines is a serious public health problem. For instance, the majority of C. perfringens isolates from chicken and human patients in the U.S. were found to be resistant to tetracycline, and tetA(P) and tetB(P) were frequently detected from them [43]. In China, tetA(P) and tetB(P) are reported to have been frequently detected from C. perfringens isolated from chickens and pigs in farms [44].

Among the Tet class ARGs, tetW and tetY were also frequently detected (Supplementary Table 5). The tetW gene was abundant in the feces of leopard cats (Fig. 5b). The tetW gene, first identified in rumen bacteria isolated from bovine and sheep [45], has a wide host range and spreads widely in environments, e.g., animal farms [46]. One hypothesis why tetW was abundant in the feces of leopard cats may be due to dietary differences from Eurasian otters. The leopard cat feeds on mammals, which may include ruminants such as deer. In fact, we detected sequences assigned to Ruminantia, a group of ruminants, only in the feces of leopard cats (Supplementary Table 7 and Supplementary Data 1), indicating ruminants as their dietary item. Since tetW was first identified from ruminant bacteria, it may be possible to hypothesize that tetW is abundant in ruminants (e.g., deer) and can be more likely transmitted to predatory animals on them (e.g., leopard cat). Conversely, tetY was abundant in the feces of Eurasian otters (Fig. 5b). Interestingly, tetY is reported to have been detected from bacteria isolated from fishes in fish farms in which tetracycline is permitted to be added to their feed [47]. The authors of this study argued that farmed fish may contribute to the transmission of tetY to humans [47]. We speculate that a similar transmission mechanism may be possible for Eurasian otters due to their piscivorous character.

As discussed, some ARGs were detected differentially by the host animal, e.g., iri from leopard cats (Fig. 5b). The iri gene is known to confer resistance to rifampin [48], which is used to treat Mycobacterium infections such as tuberculosis. We also found that Mycobacterium was more abundant in the feces of leopard cats (Fig. 4) though it was not statistically significant. Therefore, we speculated that iri was associated with Mycobacterium, and in fact, a positive correlation was observed (ρ = 0.81) (Fig. 6 and Supplementary Table 6). It is well known that domestic cats can be infected with Mycobacterium, such as Mycobacterium bovis [49]. Our study revealed that Mycobacterium spp. are also present in the wild leopard cats (Fig. 4) and correlated well with the iri gene (Fig. 6). The source of the iri gene selectively found in leopard cats is unclear. One hypothesis is that it may be related to their scavenging behavior. The leopard cats scavenge carcasses [50]. Interestingly, a study indicates the behavior of scavenging carcasses died of bovine tuberculosis as a route of infection of carnivorous animals to M. bovis [51]. It is thought that leopard cats are more likely to be exposed to Mycobacterium and therefore to the co-existing iri gene since they are more likely to scavenge carcasses than Eurasian otters.

One of the limitations of our study is that we cannot exclude the possibility of environmental contamination of some fecal samples because of our sampling procedure in which the time from defecation to fecal sampling is unknown. It is inevitable since it is almost impossible to collect the feces of these elusive animals immediately after defecation. Therefore, the microbiome and resistome may have changed in some fecal samples due to postdefecation microbial succession [52]. However, the overall trend was likely conserved since the trends in bacterial composition observed were consistent with those observed in the fresh feces of captive animals by our previous study [9]. Specifically, we observed the enrichment of Jeotgalicoccus in the otter feces and Blautia and Collinsella in the leopard cat feces (Fig. 2b), which is in congruence with our previous finding in the rearing environment [9]. Furthermore, environmental contamination is unlikely to explain our results because the number of contaminating cells would be very small compared to mammalian feces, which contains 1000 times the number of microbial cells than in a typical environmental sample (e.g., soil) [53,54]. The results suggest that the small numbers of contaminating reads are not disproportionately contributing to the ARGs and VFGs patterns observed in our data. Moreover, from a public health perspective, it is not unreasonable to assess ARGs and VFGs in feces collected sometime after defecation. It is possible that human exposure to ARGs and pathogens sourced from environmental feces occurs long after defecation.

Conclusion

We found common features in most of the ARGs and VFGs profiles in environmental feces of Eurasian otters and leopard cats in Korea. However, certain ARGs and VFGs were differentially more abundant between either of the host animals possibly due to differences in the host's gut microbiome, diet, and/or habitat, but further investigation is needed. Wildlife is thought to serve as reservoirs, melting pots, vectors, and secondary sources of ARGs, and their importance as sentinels to monitor their risk has been also recognized [2]. This study provided several important baseline information about the resistomes of the wildlife, which could be useful in assessing the ecological and human health risk posed by ARGs. However, knowledge is still lacking about how ARGs are transmitted from wildlife to human and vice versa. Future research should include studies of wildlife and environmental reservoirs in order to identify sources of ARGs and elucidate their routes to (or from) wildlife from the perspective of One Health.

Funding

This work was supported by the Midcareer Researcher Program of the (2020R1A2C1004903) (NY).

Declaration of Competing Interest

None to declare.