Mobile Antimicrobial Resistance Genes in Probiotics

Even though people worldwide tend to consume probiotic products for their beneficial health effects on a daily basis, recently, concerns were outlined regarding the uptake and potential intestinal colonisation of the bacteria that they carry. These bacteria are capable of executing horizontal gene transfer (HGT) which facilitates the movement of various genes, including antimicrobial resistance genes (ARGs), among the donor and recipient bacterial populations. Within our study, 47 shotgun sequencing datasets deriving from various probiotic samples (isolated strains and metagenomes) were bioinformatically analysed. We detected more than 70 ARGs, out of which rpoB mutants conferring resistance to rifampicin, tet(W/N/W) and potentially extended-spectrum beta-lactamase (ESBL) coding TEM-116 were the most common. Numerous ARGs were associated with integrated mobile genetic elements, plasmids or phages promoting the HGT. Our findings raise clinical and public health concerns as the consumption of probiotic products may lead to the transfer of ARGs to human gut bacteria.


Introduction
Probiotics and probiotic products have gained a worldwide reputation and popularity in our everyday lives irrespective of cultural background, geographic location or social standards. Beneficial health effects assigned to probiotics have been reported in several studies [1]. What these studies have in common is that they state that microbes carried in probiotics must remain present in the intestinal tract for a shorter or longer period of time to exert the expected beneficial effects. Nevertheless, the success of colonisation depends on several factors, thus the certainty of its realisation varies from individual to individual [2]. Recently, however, the possibility of some unfavourable or sometimes even adverse effects of probiotic consumption have also been raised [3]. Several publications indicate that bacterial strains included in probiotic compounds, powders and capsules may contain antimicrobial resistance genes (ARGs) [4][5][6][7]. Recognising that ARGs may enter into the human body by food (e.g., probiotic products), studies on the genetic characteristics of microorganisms (including bacteria) used in the food chain have been recommended by European Food Safety Authority (EFSA) in recent years [8,9].
Genes, including ARGs of the probiotic bacteria, can be transmitted to bacteria within the intestinal tract of the consumers by horizontal gene transfer (HGT). If such ARGs are received by pathogenic bacteria, the efficacy of antibiotic therapy prescribed as medical intervention for the diseases they cause may lessen. HGT can take place by transformation, conjugation or transduction. All these processes have one important property in common, namely, a DNA fragment getting introduced to a recipient cell. Apart from transformation, by which any gene can be taken up by the bacterium from its environment, the routes of HGT require special active delivery processes. By conjugation, cell-to-cell contact provides the opportunity for a copy of a plasmid to translocate to a recipient bacterium [10]. In contrast, transduction negates the necessity for cell-to-cell contact, as in this case bacteriophages act as a conduit for shuttling genes among bacteria [11]. The genetic environment of the genes, possibly ARGs, involved in the transfer has a significant influence on the efficacy of the latter two HGT processes, i.e., on the genes' mobility. The transferability of ARGs is facilitated by the presence of mobility genes in their tight genetic environment. If the genes harbour on plasmids or prophages, the chance of their transfer increases. By probiotics with supposedly mobile ARGs, the likelihood of gene transmission to other bacteria in the intestinal tract increases. In our work, we aimed to gain insight into the mobility of ARGs in probiotics for human consumption using freely available samples sequenced by other research groups or by ourselves. Currently, the few accessible data on probiotic ARG mobility originate from studies with diverse methodologies [6,[12][13][14][15][16]. Therefore, we intended to analyse the next-generation sequencing datasets of different probiotics and probiotic isolated bacterial strains with a unified bioinformatics approach.

Materials and Methods
In this study, we followed the FAO/WHO definition of probiotics, i.e., live microorganisms, which confer a health benefit on the host when administered in adequate amounts [17].

Data
For the study, we selected freely available samples from the sequencing of probiotic products for human consumption or from bacterial strains isolated from such products. The details of analysed samples are listed in Table 1. One probiotic capsule was shotgun sequenced (PRJNA644361) by the authors. All further short read datasets were obtained from NCBI SRA repository.

DNA Extraction and Metagenomics Library Preparation for PRJNA644361
Total metagenome DNA of the probiotic capsule sample was extracted using the UltraClean Microbial DNA Isolation kit from MoBio Laboratories. The quality of the isolated total metagenome DNA was checked using an Agilent Tapestation 2200 instrument. The DNA sample was used for in vitro fragment library preparation. In vitro fragment library way prepared using the NEBNext Ultra II DNA Library Prep Kit for Illumina. Paired-end fragment reads were generated on an Illumina NextSeq sequencer using TG NextSeq 500/550 High Output Kit v2 (300 cycles). Primary data analysis (base-calling) was carried out with Bbcl2fastq software (v2.17.1.14, Illumina).

Bioinformatic Analysis
Quality based filtering and trimming of the raw short reads was performed by Trim-Galore (v.0.6.6, https://github.com/FelixKrueger/TrimGalore, accessed on 22 Match 2021), setting 20 as a quality threshold. Only reads longer than 50 bp were retained and taxonomically classified using Kraken2 (v2.1.1) [18] and a database created (24 March 2021) from the NCBI RefSeq complete archaeal, bacterial and viral genomes. For this taxon assignment, the -confidence 0.5 parameter was used to obtain more precise species level hits. The taxon classification data was managed in R [19] using functions of the packages phyloseq [20] and microbiome [21]. The preprocessed reads were assembled to contigs by MEGAHIT (v1.2.9) [22] using default settings. The contigs were also classified taxonomically by Kraken2 with the same database as above. From the contigs having more than 500 bp, all possible open reading frames (ORFs) were gathered by Prodigal (v2.6.3) [23].
The protein translated ORFs were aligned to the ARG sequences of the Comprehensive Antibiotic Resistance Database (CARD, v.3.1.1) [24,25] by Resistance Gene Identifier (RGI, v5.1.1) with Diamond [26] The ORFs classified as perfect or strict were further filtered with 90% identity and 90% coverage. All nudged hits were excluded. The integrative mobile genetic element (iMGE) content of the ARG harbouring contigs was analysed by MobileElementFinder (v1.0.3) [27]. Following the distance concept of Johansson et al. [27] for each bacterial species, those with a distance threshold defined within iMGEs and ARGs were considered associated. In the MobileElementFinder database (v1.0.2) for Escherichia coli, the longest composite transposon (cTn) was the Tn1681. In the case of this species, its length (24,488 bp) was taken as the cut-off value. For Lactococcus lactis, this threshold was the length of the Tn5721 transposon, 11,256 bp. For Enterococci, the database contained cTn, the Tn6246 (5147 bp) transposon, in E. faecium only. The same threshold was used for E. faecalis contigs. As the database neither contains species-level, nor genus-level cTn data for Bacillus, Bifidobacterium and Streptomyces species, a general cut-off value was chosen for the contigs of these species. This value was declared as the median of the longest cTns per species in the database (10,098 bp). The average nucleotide identity (ANI) was calculated for the region of iMGE and associated ARGs by FastANI (v1.32) [28]. The plasmid origin probability of the contigs was estimated by PlasFlow (v.1.1) [29]. The phage content of the assembled contigs was predicted by VirSorter2 (v2.2.1) [30]. The findings were filtered for dsDNAphages and ssDNAs. All data management procedures, analyses and plottings were performed in R environment (v4.0.4) [19].

Results
The analysis of the sequencing data from 20 isolates and 27 metagenomic (multimicroorganism) samples (Table 1) is summarised in three sections. Following the presentation of the bacteriome and the identified AGRs (resistome), predictions regarding the mobility potential of ARGs were also summarized based on genetic characteristics that may play a significant role in HGT. If integrated mobile genetic elements (iMGE) are identified in the sequence context of an ARG, its greater mobility can be assumed. The case is the same if the contig harbouring an ARG are derived as plasmid or phage originated. In the mobilome section, we summarise these results.

Resistome
The median length of the filtered contigs harbouring ARGs constructed by de novo assembly was 102,711 bp (IQR: 105,696). The number of ARGs found on the contigs ranged from 1 to 12. Besides 182 perfect ARG matches, a further 225 hits were classified strict (RGI) and met the criteria of having 90% coverage and 90% sequential identity.
ARGs were detected in all metagenomic samples and in few isolates ( Figure 2). The majority of isolates (s01, s02, s03, s04, s05, s06, s07, s08, s09, s10, s20) contained no ARG. The highest number of ARGs was found in samples s14-s19, obtained from sequencing six Escherichia coli strains isolated from the same probiotic product. It is important to highlight that we also found the H-NS gene in these samples which is not indicated in the figure, as its effect is anti-AMR. The most common ARGs were the rpoB mutants conferring resistance to rifampicin, TEM-116 and tet(W/N/W) genes, detected in 18, 15 and 13 samples, respectively.  Table 1.

Mobilome
The frequencies of iMGEs, phages and plasmids associated with ARGs by bacteria of origin are summarised in Figure 3.

Coexistence of ARGs and iMGEs
Based on the distance method proposed by Johansson et al. (2021) [27] iMGE associated ARGs were detected in three species (Bifidobacterium animalis, Enterococcus faecalis and Escherichia coli). In seven metagenomic samples (m01, m02, m03, m07, m16, m17, m24) we found tet(W/N/W) associated with ISBian1 insertion sequence on contigs classified as B. animalis originated. In two further samples (m02, m06) on E. faecalis originated contigs, tetM is linked to the transposon Tn6009. The ARG mdtG in the E. coli sample s14 and the ARG ugd in s15 are associated with IS3 and IS100, respectively. On two different contigs in the sample s17, multiple ARGs were detected with iMGE. One of them has the ISKpn24 associated with mdtE and mdtF. The other one has the IS102 linked to emrY, emrK, evgA and evgS genes. According to the average nucleotide identity (ANI) analysis most of the contig region of iMGE and associated ARGs had a high level of conservation (ANI > 97%). Nevertheless, both contigs classified as E. faecalis originated showed ANIs below 80%.

Phages
By phage prediction, only dsDNAphages were detected. One contig, classified as Bacillus subtilis from the m05 metagenomic sample, contained prophage harbouring gene aadK. One prophage in predicted Enterococcus faecalis originated contig was found in sample m04 having gene efrA. The same content was detected in sample m01 on contigs classified as E. faecalis. All three E. faecalis isolates (s11, s12, s13) contained contigs harbouring the gene efrA within a prophage. In sample m17, one E. coli classified contig had the gene TEM-116, while a Lactococcus lactis classified one carried the gene lmrD within a prophage. All the E. coli isolates contained contigs with prophages harbouring ARG. In the sample s17 and s19 the mdfA gene is presented within a prophage. The sample s15 contains contigs harbouring prophage with the gene marA, marR. The sample s16 harbours contigs with prophage having genes emrK, emrY, evgA. The gene ampC was found in sample s15, while the gene cpxA in samples s14 and s18 within prophages.

Discussion
The results presented demonstrate that the bacteria of probiotics may not only carry significant amounts of ARGs, but in numerous cases, those genes may also be mobile, thereby contributing to their spread to other bacteria and having possible consequences on the antibiotic treatment efficacy.
Bacterial genera identified in the metagenomic samples also appear in many probiotic related articles of the current international literature. Various species of Bacilli, Bifidobacteria, Enterococci, Lacticaseibacilli, Lactiplantibacilli, Lactobacilli, Lactococci, Ligilactobacilli, Limosilactobacilli and Streptococci are the core members of commercial probiotic bacterial communities [31][32][33][34][35][36][37][38][39][40]. Two identified bacterial genera (Sphingobacterium, Weizmannia) in the various samples are less frequent probiotic components. The possibility of exploiting Sphingobacteria in probiotic foods was previously mentioned based on the characterization of flour and batter samples of sorghum and pearl millet [41]. Members of the genus were detected by the high-throughput sequence analyses of fermented beverages [42]. Probiotic Weizmannia species (e.g., former Bacillus coagulans) have recently been reclassified [43], and have an unquestionable probiotic significance [44]. It is important to note that there may be notable differences in the gene pool between strains of particular species, so the results presented do not mean that all strains of a given species contain the genes identified here.
While at least one ARG was found in each metagenomic sample, less than half of the isolates contained any of them. No ARG was detected in Lacticaseibacillus rhamnosus, Lactiplantibacillus plantarum, Lactobacillus delbrueckii subsp. bulgaricus, Limosilactobacillus fermentum, Pseudomonas sp. RGM2144 or Streptococcus thermophilus. Contigs originating from Bacillus subtilis, Bifidobacterium animalis, Bifidobacterium bifidum, Bifidobacterium breve, Bifidobacterium longum, Enterococcus faecalis, Enterococcus faecium, Escherichia coli, Lactococcus lactis and Streptomyces albulus each contained at least one ARG.
The available literature was screened to evaluate our findings and gain reliable knowledge of the ARGs that could have been attached to bacteria at the species level. All ARGs found in Bacillus subtilis (aadK, B. subtilis mprF, B. subtilis pgsA with mutation conferring resistance to daptomycin, bmr, lmrB, mphK, vmlR, ykkC, ykkD) have previously been identified in B. subtilis and many of them were even reported to be specific for this species or the Bacillus genus [45][46][47][48][49][50][51]. In the Bifidobacterium genus, ARGs were associated with four species (B. animalis, B. bifidum, B. breve and B. longum). None of the B. animalis, B. bifidum, B. breve and B. longum related B. adolescentis rpoB mutants conferring resistance to rifampicin and tet(W/N/W) are specific for the identified species but both genes have previously been described in them [6,[52][53][54][55]. B. bifidum ileS conferring resistance to mupirocin reported in B. bifidum supposedly cannot be exclusively linked to this species of the genus, but it had been identified in it before [56]. Out of the Enterococcus faecalis deriving genes, dfrE was first identified in E. faecalis [57], but according to a recent study it is not exclusive to this species any more [58]. The genes efrA and efrB have been described in E. faecalis and E. faecium [59,60]. Gene emeA has only been identified in E. faecalis so far [59]. Apart from E. faecalis, lsaA has been attached to Streptococcus agalactiae, while tetM appears in a broad spectrum of bacterial species [61][62][63][64][65]. All three ARGs (AAC(6')-Ii, eatAv, msrC) associated with E. faecium have been previously published as appearing in this species, and the first two are even specific for it [66][67][68][69]. All ARGs originating from Escherichia coli (acrB, acrD, acrE, acrF, acrS, bacA, baeR, baeS, cpxA, CRP, emrA, emrB, emrK, emrR, emrY, eptA, E. coli acrA, E. coli acrR with mutation conferring multidrug antibiotic resistance, E. coli ampC beta-lactamase, E. coli ampC1 beta-lactamase, E. coli ampH beta-lactamase, E. coli emrE, E. coli GlpT with mutation conferring resistance to fosfomycin, E. coli marR mutant conferring antibiotic resistance, E. coli mdfA, E. coli soxR with mutation conferring antibiotic resistance, E. coli soxS with mutation conferring antibiotic resistance, evgA, evgS, gadW, gadX, kdpE, marA, mdtA, mdtB, mdtC, mdtE, mdtF, mdtG, mdtH, mdtM, mdtN, mdtO, mdtP, msbA, PmrF, TEM-116, TolC, ugd, YojI) have previously been described in this species and many of them are even specific to it, according to the Comprehensive Antibiotic Resistance Database (CARD) [24,25]. Gene lmrD, the only ARG deriving from Lactococcus lactis has been identified in this species along with some others [70,71]. Even though AAC(3)-IV has been identified in several studies [72,73], according to our knowledge this is the first time it has been detected in Streptomyces albulus.
Gene TEM-116, which is often referred to as a clinically significant extended-spectrum beta-lactamase (ESBLs), was the most frequently identified finding in our study. ESBLs are most commonly defined as the members of a ubiquitous enzyme family that is capable of conferring resistance to penicillins, first-, second-and third-generation cephalosporins and aztreonam, and of being impeded by beta-lactamase inhibitors such as clavulanic acid [74]. The 400 TEM variants that have been identified so far, can be disclosed in two clusters with one deriving from TEM-1 (the first TEM protein to be described) and one linked to TEM-116 as a progenitor [75]. In line with our findings, gene TEM-116 is reported to be present worldwide harbouring in the conjugative plasmids of a wide range of Gramnegative bacteria. Despite its wide geographical dissemination, establishment on multiple plasmids and centrality in the TEM family network indicating it is a naturally occurring enzyme with microbiologically proven ESBL characteristics [76,77], some concerns have arisen about its designation, after the gene was found in non-ESBL producing Klebsiella pneumoniae strains [78]. Moreover, commercial Taq polymerases used in PCRs may be contaminated with bla TEM−116 DNA which could lead to the erroneous identification of the gene in samples that do not actually contain it [79,80]. In our study, each sample in which this gene was detected originated from the same bioproject (PRJNA542229). As the samples come from different dietary supplements, one may interpret that this finding is an artefact or contamination as a consequence of some sample preparation steps. Nevertheless, as more detailed information on sample preparation is not available, this issue cannot be resolved.
As seen above, and as described in other publications [81], there is still a great deal of variation in details which need to be clarified by the interpretation of ARGs. Nevertheless, the suspicion that the identified ARGs may undermine the efficacy of several antibiotic classes, including acridine dye, aminocoumarins, aminoglycosides, benzalkonium chloride, carbapenems, cephalosporins, cephamycins, diaminopyrimidines, fluoroquinolones, fosfomycins, glycylcyclines, lincosamides, macrolides, monobactams, mupirocins, nitroimidazoles, nucleosides, oxazolidinones, penams, penems, peptides, phenicols, pleuromutilins, rhodamines, rifamycins, streptogramins, tetracyclines and triclosans raises some clinical considerations. According to the latest CDC report on antimicrobial use in the U.S., amoxicillin (penam), azithromycin (aminoglycoside), amoxicillin and clavunalic acid (penam, increased activity), cephalexin (cephalosporin) and doxycycline (tetracycline) are the most commonly administered compounds [82]. Moreover, based on the latest WHO report on global antimicrobial use, amoxicillin (penam), ciprofloxacin (fluoroquinolon), sulphametoxazole and trimethoprim are the most commonly prescribed oral drugs and ceftriaxone (cephalosporin), gentamicin (aminogylcoside) and benzylpenicillin (penam) are the most commonly used parenteral compounds in 4 surveyed countries of the African region. In six countries of the region of the Americas, amoxicillin (penam), cefalexin (cephalosporin) and doxycycline (tetracycline) are the antibiotics with the highest oral consumption rates and ceftriaxone (cephalosporin), oxacillin (penam) and gentamicin (aminogylcoside) are the ones with the highest parenteral use. In the European region, reports were made of 46 countries. Among orally administered antibiotics, amoxicillin (penam), amoxicillin and beta-lactamase inhibitors (penam, increased activity) and doxycycline (tetracycline) are the top 3 compounds, while ceftriaxone (cephalosporin), gentamicin (amynoglycoside), and cefazoline (cephalosporin) are the most common parenteral ones. Amoxicillin (penam), azithromycin (macrolide) and amoxicillin and beta-lactamase inhibitors (penam, increased activity) are the most commonly consumed oral antibiotics and ceftriaxone (cephalosporin), benzathine benzylpenicillin (penam) and procaine benzylpenicillin (penam) are the top 3 parenterally administered agents in the Eastern Mediterranean region. In the six surveyed countries of the Western Pacific region amoxicillin (penam), doxycycline (tetracycline) and amoxicillin and beta-lactamase inhibitors (penam, increased activity) are the most commonly prescribed oral antibiotics, while cefazolin (cephalosporin), ceftriaxone (cephalosporin) and cefuroxime (cephalosporin) are the most frequently used parenteral compounds [83]. Many of the most highly prioritized antibiotics could be affected by the presence of the detected ARGs. Meanwhile, out of the 15 antibiotic groups mentioned in the latest WHO report on critically important antimicrobials (CIA) for human medicine, nine (aminoglycosides, carbapenems and other penems, cephalosporins, glycylcyclines, macrolides, monobactams, oxazolidinones, penicillins of various cathegories, quinolones) could possibly be affected by the ARGs identified in the various samples [83].
It is important to underline that all the six E. coli isolates contained the gene H-NS, which plays a crucial role in the global gene regulation of various bacteria, including this species. The expression of a wide variety of genes is repressed by H-NS, and its deletion increases AMR and decreases drug accumulation. Even though this gene is stored in CARD [24,25], its functional effect is adverse to that produced by ARGs [84].
If ARGs are transmitted from probiotic bacteria to pathogenic bacteria within the consumer's body, they may reduce the effectiveness of antibiotic therapy on the diseases participating pathogenic bacteria cause. The execution of gene transfer processes is more likely among bacteria that are in close physical proximity to each other and if the ARGs are associated to a mobile genetic environment. According to our results a considerable number of ARGs, such as those which are iMGEs-linked or have resided in plasmids or prophages.
The co-occurence of tet(W/N/W) and ISBian1 is in line with the findings of Rozman et al. [6], according to which all genomes of B. animalis (subspecies lactis or animalis) (n = 42) available in 2019 contained this gene. Moreover, by the investigation of the mobility characteristics of tetW, out of the transposases belonging to the family of the insertion sequences, ISBian1 seemed to be subspecies dependent in B. animalis subsp. lactis and flanking tetW in the majority of the strains [6]. Our results of tetM linking to the transposon Tn6009 in E. faecalis is consistent with finding of Zangue et al. in South-African fecal samples [85].
In two samples, contigs harbouring tet(W/N/W) originating from Bifidobacterium longum and Bifidobacterium animalis were predicted to belong to plasmids. Several studies reported a wide prevalence of the tetW gene in Bifidobacteria [6,12,86,87]. While the co-occurrence of tetW and its flanking transposase is a common genetic feature of B. animalis, previous reports lack the identification of plasmids in B. animalis, even though the gene was associated with plasmids in other bacterial species [88]. Despite AAC(6')-Ii deriving from E. faecium being located in the chromosome in previous studies and it being defined as a chromosome-borne ARG on CARD [24,25,89], our research indicates it may take place in a plasmid. An E. faecium-associated contig contained gene msrC. According to the available literature, msrC is a chromosomal-encoded gene that is mentioned as an intrinsic property of E. faecium strains [24,25,90]. While the expected bacterial species of origin was confirmed, our finding raises the likelihood of the gene being connected to a plasmid as well. In 15 samples, E. coli-originated contigs harboured the gene TEM-116. Plasmid origin is a common feature of ESBL genes such as TEM-116 according to several publications and is often referred to as a feature to facilitate their quick spread [91][92][93]. In the E. coli isolate sample s15, one contig had the marA and marR genes. These widespread multiple antibiotic resistance genes had been identified on plasmids before [94]. The gene efrA harbouring in contigs with a prediction of phage origins were identified in all publicly available E. faecalis genome sequences by Panthee and colleges too, along with a large set of phages in the genomes [95].
As our results derive from in silico data analysis, it is only possible to describe the features that prove and facilitate presence and mobility of the genes. Whether or not the identified genes operate in the bacterial strains of a given probiotic cannot be determined. In order to clarify this, additional functional, e.g., gene expression studies, should be performed.
An important aspect to take into consideration by the interpretation of the ARG occurrence in probiotics is that constituent strains can often naturally be, or rendered multiresistant, so that they can be co-administered with oral antibiotics and reduce gastrointestinal side effects [96,97]. In our study we could not distinguish whether the examined samples contained the ARGs for this purpose. Moreover, as ARGs were found in the vast majority of the samples tested, not a negligible proportion of them, it is possible that the presence of ARGs in bacteria may also play a role in their probiotic effect. ARGs play a role in defence against antibiotics and may provide general fitness against specific toxic effects for bacteria [98,99]. One may make an analogy with earlier practice. In livestock farming, antibiotics have been widely used as feed supplements for yield enhancement on a purely empirical basis. By this practice, antibiotics have put pressure on the gut bacteria and selected for resistant strains. As a result, animal feed efficiency and production indicators have improved. When probiotics are consumed, the expectation is that the "good" microorganisms, bacteria will colonise the gut. In numerous animal husbandry areas (e.g., broiler chicken production), the producers try to achieve this by continuous probiotic feeding. If these probiotics also contain bacterial strains harbouring ARGs, they achieve very similar results as before with the selective effect of antibiotic utilisation. If it is true that certain ARGs are essential for the efficacy of probiotic bacteria, then the selection of strains should be carried out with consideration of the human health consequences. That is, bacterial strains that contain ARGs having no significant influence on human antimicrobial therapy efficiency should be used. However, based on our results, it can also be suggested that bacteria that do not contain ARGs at all can be used as probiotic components. To have a more detailed insight into this topic, several further studies would be needed. For instance, they could also focus on reducing the mobility of genes whose presence may be necessary for the probiotic nature of particular bacteria. Based on the results, we consider it essential to monitor the ARG content of probiotic preparations and their mobility characteristics in the fight against antimicrobial resistance.