Comprehensive genome data analysis establishes a triple whammy of carbapenemases, ICEs and multiple clinically relevant bacteria.

Carbapenemases inactivate most β-lactam antibiotics, including carbapenems, and have frequently been reported among Enterobacteriaceae, Acinetobacter spp. and Pseudomonas spp. Traditionally, the horizontal gene transfer of carbapenemase-encoding genes (CEGs) has been linked to plasmids. However, given that integrative and conjugative elements (ICEs) are possibly the most abundant conjugative elements among prokaryotes, we conducted an in silico analysis to ascertain the likely role of ICEs in the spread of CEGs among all bacterial genomes (n=182 663). We detected 17 520 CEGs, of which 66 were located within putative ICEs among several bacterial species (including clinically relevant bacteria, such as Pseudomonas aeruginosa, Klebsiella pneumoniae and Escherichia coli). Most CEGs detected within ICEs belong to the IMP, NDM and SPM metallo-beta-lactamase families, and the serine beta-lactamase KPC and GES families. Different mechanisms were likely responsible for acquisition of these genes. The majority of CEG-bearing ICEs belong to the MPFG, MPFT and MPFF classes and often encode resistance to other antibiotics (e.g. aminoglycosides and fluoroquinolones). This study provides a snapshot of the different CEGs associated with ICEs among available bacterial genomes and sheds light on the underappreciated contribution of ICEs to the spread of carbapenem resistance globally.

Data Summary

All the bacterial genomes scanned in this study have been deposited previously in the National Center for Biotechnology Information (NCBI) genome database and are listed in the supplementary tables. The 66 extracted ICEs (in fasta format) and the outputs for the profile HMMs scanned on the 386 putative MGEs identified in this study have been deposited on figshare at https://figshare.com/projects/_Comprehensive_genome_data_analysis_establishes_a_triple_whammy_of_carbapenemases_ICEs_and_multiple_clinically-relevant_bacteria/78369.

Carbapenems are commonly used to treat severe infections in humans. Resistance is often mediated by carbapenemases. These enzymes degrade carbapenems and are frequently present in plasmids. Here, we demonstrate that common carbapenemase-encoding genes (CEGs) found in clinical isolates (e.g. bla KPC, bla GES, bla IMP, bla NDM, bla VIM) can also be located within integrative and conjugative elements (ICEs). CEG-bearing ICEs belong to three mating pair formation families. These mobile elements may be particularly important in bacteria where plasmids do not seem to play a significant role in the spread of antibiotic resistance genes, such as spp. This study considerably expands our knowledge of the repertoire of CEGs-bearing ICEs among clinically relevant bacterial pathogens, such as , and .

Introduction

Due to the importance of carbapenems for the treatment of severe infections in humans, the World Health Organization (WHO) stated that these antibiotics should be reserved for infections caused by multidrug-resistant Gram-negative bacteria in humans [1]. Recently, the same agency presented a list of bacterial pathogens for which new antibiotics research and development are urgently required, and the top priority pathogens were the carbapenem-resistant strains of , and [2].

The evolution of carbapenem resistance in bacteria is often driven by the horizontal gene transfer (HGT) of carbapenemase-encoding genes (CEGs) [3, 4]. Carbapenemases are beta-lactamases that are able to hydrolyze carbapenems as well as most other beta-lactam antibiotics. These enzymes are members of serine beta-lactamases classes A and D, and the class B metallo-beta-lactamases [5]. The CEGs are often located on integrons or transposons that themselves target mobile genetic elements (MGEs) such as plasmids [3, 4], which makes the dissemination of these genes unpredictably complex within bacterial communities. Recently, it was demonstrated that another type of MGE, the integrative and conjugative elements (ICEs), are likely to play a significant role as vehicles for the dissemination of CEGs among [6]. Besides genes conferring antibiotic resistance, ICEs may harbour additional cargo genes that provide an adaptive advantage over other elements. Some of these examples include the presence of the siderophore yersiniabactin encoded within the ICEKp in hypervirulent clonal group CG23 from [7]; the Tn5252-related ICEs carrying bacteriocin clusters in [8]; the type I-C CRISPR-Cas systems identified within pKLC102-like ICEs in [9]; and the type III restriction–modification systems from SXT/R391-related ICEs in spp. [10].

ICEs are self-transmissible MGEs that can integrate into and excise from the genome (like transposons and phages) and can exist as circular, sometimes replicable, extrachromosomal elements and be transferred by conjugation (like some plasmids) [11-14]. ICEclc from [15], SXT from [16], pKLC102 from [17] and Tn4371 from [18] are among the most well studied ICEs. ICEs appear to have a bipartite lifestyle that shifts between vertical and horizontal transmission [12, 19, 20]. HGT by conjugation requires three main components: a relaxase (MOB), a mating pair formation (MPF) system and a type IV coupling protein, with the last two forming a spanning-membrane multi-protein complex named the type IV secretion system (T4SS) [21]. To date, eight MPF classes have been proposed (B, C, F, FA, FATA, G, I and T), based on the phylogeny of VirB4, the only ubiquitous protein among the T4SS. The MPFT is widely distributed in both conjugative plasmids and ICEs, while MPFF is more prevalent in plasmids and MPFG on ICEs [11].

Given that ICEs have been identified in most bacterial clades and have been proposed to be more prevalent than conjugative plasmids [11], we conducted an in silico analysis to explore the distribution of CEG-bearing ICEs among all sequenced bacterial genomes available in the National Center for Biotechnology Information (NCBI). Our results demonstrate that CEG-bearing ICEs belong to three MPF families and are primarily located in several clinically relevant bacterial pathogens. Our analysis highlights the importance of investigating these elements thoroughly as important vehicles for the spread of antibiotic resistance (AR), particularly with respect to carbapenems.

Methods

Bacterial genome and carbapenemase search

In Fig. 1, we present the workflow used in this study, from the acquisition of bacterial genomes to the identification and characterization of putative ICEs. We retrieved all bacterial genomes available in the NCBI Reference Sequence Database (RefSeq, accessed on 21 March 2020), including complete and draft genome sequences, using ncbi-genome-download v0.2.12 (https://github.com/kblin/ncbi-genome-download). We downloaded over 6000 curated AR protein sequences from the AMRfinder database (https://ftp.ncbi.nlm.nih.gov/pathogen/Antimicrobial_resistance/AMRFinderPlus/database/3.6/2020-01-22.1/) [22] and built an in-house database only including the proteins that code for a carbapenemase (n=1014, Table S1, available in the online version of this article). We then blasted the genomes against the extracted carbapenemases using diamond v0.9.29.130 (http://www.diamondsearch.org/index.php) [23], using minimum 100 % identity and subject cover and with the sensitive mode enabled.

Fig. 1.

Overview of the workflow followed in this study. All assemblies available in NCBI RefSeq were downloaded and blasted against an in-house database of carbapenemases using diamond blastx (step 1). NCBI annotated proteins from CGE-bearing genomes were then extracted (step 2) and used for the identification of relaxase and serine or tyrosine recombinase (step 3). Search of directed repeats and delimitation of putative ICEs was also performed. CONJscan was used to identify the MPF family of each element. We then looked for AR genes, restriction–modification systems, CRISPR arrays and their associated (Cas) proteins, as well as secondary metabolites within extracted ICEs. We also characterized the functional annotations of their proteomes and the MLST of the genomes carrying a CEG-bearing ICE. Abbreviations: AR, antibiotic resistance; CEG, carbapenemase-encoding gene; HMM, hidden Markov model; ICE, integrative and conjugative element; MPF, mating-pair formation; RM, restriction–modification.

Tracing ICEs among the bacterial genomes

The RefSeq protein files from the CEG-bearing genomes identified by diamond were extracted. We used the hmmsearch function of the HMMER3 software package v3.3 (http://hmmer.org/) [24] to search the proteomes against the standalone version of MOBfamDB, a curated hidden Markov models (HMM) relaxase database (https://castillo.dicom.unican.es/mobscan_about/) [25]. We also used this function to search the pfam v33.0 database for tyrosine or serine recombinase accessions numbers (Pfam IDs PF00589 and PF07508). The hmmsearch command was used with default parameters and an E-value threshold of 0.01. The CEG-bearing genomes with relaxase and integrase hits were further analysed. We used the Find Repeats tool from Geneious Prime 2020.0.4 (https://www.geneious.com) to inspect the hits for direct repeats. To delimit CEG-harbouring ICEs, we manually scanned candidate terminal regions with direct repetitions of the 3′ end from tRNA genes located next to the integrase-encoding gene. When no tRNA gene was identified next to this gene, we scanned the presence of direct repeats next to the integrase-encoding gene and next to candidate terminal regions. To assist in identifying putative terminal regions, we looked for blocks of DNA with variation in GC content. To predict the MPF families, the translated coding sequences of delimited ICEs were analysed on the standalone CONJscan module of MacSyFinder v1.0.5 (https://github.com/gem-pasteur/macsyfinder) [26, 27]. To identify the multi-locus sequence type of the genomes containing CEG-bearing ICEs, we used mlst v2.16.1 (https://github.com/tseemann/mlst), which scans the genomes against PubMLST typing schemes (https://pubmlst.org/) [28].

Characterization of the CEG-bearing ICEs

Screening of AR genes among ICEs was performed using amrfinder v3.6.10 (https://github.com/ncbi/amr) [22]. The genetic platforms involved in the acquisition of CEGs by ICEs were annotated using Galileo AMR (https://galileoamr.arcbio.com/mara/) (Arc Bio, Cambridge, MA, USA) [29]. We ran our extracted ICEs against REBASE (http://rebase.neb.com/rebase/rebase.html) to look for restriction–modification systems [30]. We used CRISPRCasFinder (https://crisprcas.i2bc.paris-saclay.fr/CrisprCasFinder/Index) to look for CRISPR (clustered regularly interspaced short palindromic repeats) arrays and their associated (Cas) proteins within ICE sequences [31]. Secondary metabolite biosynthetic gene clusters were traced using antismash v5.1.2 (https://antismash.secondarymetabolites.org/) [32]. We used eggNOG-mapper v2 (http://eggnog-mapper.embl.de/) for functional annotation based on orthology assignments of the ICE proteomes [33].

Results

Carbapenemase-encoding genes are mainly found in proteobacteria

We retrieved a total of 182 663 bacterial genomes from NCBI (16 798 complete genomes and 165 865 genomes assembled at the chromosome, scaffold or contig level). We identified a total of 17 520 CEGs, with 1422 CEGs on 1236 complete genomes (including 512 chromosomes and 724 plasmids) and 16 098 CEGs on 16 038 draft genomes (Table S2). We identified a total of 377 carbapenemase variants among the 17 520 hits. Our results show that CEGs are mostly located on and dominated by clinically relevant pathogens such as , , and (Table S2). These genomes encode a wide diversity of carbapenemases, including OXA-23 (15.6%, n=2 739/17 520), KPC-2 (13.2%), KPC-3 (8.9%) and NDM-1 (6.7%). As we are tracing MGEs integrated in the chromosome, the 724 plasmid hits in complete genomes were excluded from the analysis. Among the hits on the 16 038 draft genomes, we filtered out sequences with the word ‘plasmid’ present on the fasta header (n=131, Table S3). To maximize the chances of detecting entire ICEs, we filtered out sequences shorter than 40 kb (n=10 050). All the excluded ones are available for analysis in Table S4. The remaining sequences from draft genomes (n=5 857) were inspected for the presence of CEG-harbouring ICEs.

A large proportion of CEG-bearing ICEs belong to three families and target clinically relevant Gram-negative bacteria

We identified a total of 66 putative ICEs, including 42 newly characterized elements associated with 17 different CEGs (Table 1 and Fig. 2). We could predict the boundaries from 55 of these elements (Table S5). The terminal region of the remaining 11 putative ICEs could not be determined due to a fragmented contig or assembly gaps within the sequence. Nearly half of the putative ICEs (48.5%, n=32/66) were integrated at the 3′ end of a tRNAGly gene. Integration next to random genes was also observed (Fig. 3 and Table S5). The bacterial hosts housing these elements belong to 23 sequence types (STs) (Table S5).

Table 1.

Diversity and characterization of carbapenemase-encoding genes in integrative and conjugative element-associated genomes

MPF family	CEG	Integrase type	Relaxase type	Bacterial species	MGEs flanking the CEGs	ICE length (kb)*	References
MPFG (n=43)	AFM-1 (n=1)	INT_P4_C	MOBH	Bordetella trematum	Flanked by IS91 family ISs	130	This study
	DIM-1 (n=1)	INT_P4_C	MOBH	Pseudomonas aeruginosa	Class I In flanked by IS6100	89	This study, [6]
	GES-5 (n=5)	INT_P4_C	MOBH	Pseudomonas aeruginosa	Class I In within Tn3-like Tn	93–116	This study, [6]
	GES-24 (n=1)	INT_P4_C	MOBH	Pseudomonas aeruginosa	Class I In within Tn3-like Tn	63	This study
	IMP-1 (n=4)	INT_P4_C, INT_Rci_Hp1_C	MOBH	Pseudomonas aeruginosa	Class I In within Tn3-like Tn	76–109	This study, [6]
	IMP-13 (n=9)	INT_P4_C, INT_Rci_Hp1_C	MOBH	Pseudomonas aeruginosa	Class I In within Tn3-like Tn	65–89	This study, [6]
	IMP-14 (n=2)	INT_P4_C	MOBH	Achromobacter xylosoxidans	Class I In within Tn3-like Tn	106–122	This study
	IMP-16 (n=1)	INT_P4_C	MOBH	Pseudomonas monteilii	Class I In within Tn3-like Tn	86	This study
	IMP-54 (n=1)	INT_P4_C	MOBH	Pseudomonas aeruginosa	Class I In within Tn3-like Tn	91	This study
	KPC-2 (n=1)	INT_P4_C	MOBH	Pseudomonas aeruginosa	Tn4401	115	This study, [43]
	NDM-1 (n=12)	INT_P4_C	MOBH	Pseudomonas aeruginosa, Pseudomonas asiatica, Morganella morganii	Next to (or flanked by) IS91 family ISs	97–167	This study, [44, 45]
	VIM-2 (n=1)	INT_P4_C	MOBH	Pseudomonas aeruginosa	Class I In within Tn3-like Tn	65	This study
	VIM-4 (n=4)	INT_P4_C	MOBH	Pseudomonas aeruginosa, Klebsiella pneumoniae, Alcaligenes faecalis	Class I In within Tn3-like Tn	88–102	This study
MPFT (n=16)	AFM-1 (n=1)	INT_P4_C	MOBP	Stenotrophomonas maltophilia	Flanked by IS91 family ISs	77	This study
	KPC-2 (n=1)	INT_Rci_Hp1_C	MOBP	Pseudomonas sp.	Next to ISKpn6 and IS26	62	This study, [6]
	NDM-1 (n=3)	INT_P4_C, INT_Rci_Hp1_C	MOBP	Pseudomonas aeruginosa, Pseudomonas asiatica	Flanked by IS91 family ISs	73–74	This study, [35]
	SPM-1 (n=11)	INT_Rci_Hp1_C	MOBP	Pseudomonas aeruginosa	Flanked by ISCR3-like elements	44–58	This study, [36]
MPFF (n=3)	IMP-6 (n=1)	DNA_BRE_C	MOBF	Escherichia coli	Class I In within Tn3-like Tn	118	This study
	KPC-2 (n=2)	DNA_BRE_C	MOBF	Klebsiella pneumoniae	Tn4401	75	This study
Incomplete MPFG (n=2)	IMP-8 (n=1)	DNA_BRE_C	MOBH	Enterobacter cloacae	Class I In within Tn3-like Tn	124	This study, [46]
	NDM-1 (n=1)	INT_P4_C	MOBH	Pseudomonas aeruginosa	Flanked by ISCR3-like elements	98	This study
Incomplete MPFT (n=1)	KPC-2 (n=1)	INT_Rci_Hp1_C	na	Achromobacter sp.	Next to ISKpn6 and IS26	43	This study
Incomplete MPFF (n=1)	NDM-5 (n=1)	DNA_BRE_C	na	Escherichia coli	Flanked by ISCR1 and IS26	54	This study, [47]

*Some sequences are not complete. For more details, please refer to Table S5.

CEG, carbapenemase-encoding gene; DNA_BRE_C, DNA breaking–rejoining enzymes, C-terminal catalytic domain; ICE, integrative and conjugative elements; In, integron; INT_P4_C, P4 integrase, C-terminal catalytic domain; INT_Rci_Hp1_C, shufflon-specific DNA recombinase Rci and bacteriophage Hp1_like integrase, C-terminal catalytic domain family; MPF, mating pair formation; na, no hit; ST, sequence type.

Fig. 2.

Sankey diagram showing the contribution of different MPF families to the spread of CEGs among several bacterial genomes. The left, centre and right axes represent the association between the identified carbapenemases, MPF type and the bacterial genus, respectively, while the width of each connection is proportional to the number of positive hits. Abbreviations: Ax, Achromobacter xylosoxidans; Af, Alcaligenes faecalis; Bt, ; Ec, Escherichia coli; Kp, Klebsiella pneumoniae; Mm, ; Pas, ; Pm, ; Psp, sp.; Sm, .

Fig. 3.

Schematic representation of the different insertion site/integrase/CEG/MPF/relaxase profiles identified in this study. The figure is not to scale and the relative position of the different modules (integration, conjugation, accessory CEGs) is for illustrative purposes to show the various relationships observed. Abbreviations: CDS, Coding Sequence; CEGs, carbapenemase-encoding genes; DNA_BRE_C, DNA breaking–rejoining enzymes, C-terminal catalytic domain; INT_P4_C, P4 integrase, C-terminal catalytic domain; INT_Rci_Hp1_C, shufflon-specific DNA recombinase Rci and bacteriophage Hp1_like integrase, C-terminal catalytic domain family; MPF, mating-pair formation; tRNA, transfer RNA gene.

sp.

sp.

Using the CONJscan module of MacSyFinder we identified the MPF family for 62 hits (incomplete MPF classes were predicted for the remaining 4 hits) and we noted that these hits belong to 3 families: MPFG (69%, n=43/62), MPFT (26%) and MPFF (5%) (Fig. 2, Tables 1, S5 and S6). In our results, the MPFG class was only associated with MOBH, the MPFT class with MOBP and the MPFF class with MOBF (Fig. 3). All ICEs identified here carried a tyrosine recombinase, with the majority of them (56%, n=37/66) belonging to the P4 integrase, C-terminal catalytic domain family (INT_P4_C). The shufflon-specific DNA recombinase Rci and bacteriophage Hp1_like integrase, C-terminal catalytic domain family (INT_Rci_Hp1_C) and the DNA breaking–rejoining enzymes, C-terminal catalytic domain family (DNA_BRE_C) integrases were also identified in our collection (36 and 8 %, respectively) (Tables 1 and S5). INT_P4_C and INT_Rci_Hp1_C integrases were associated with ICEs belonging to the MPFG and MPFT classes, while DNA_BRE_C was found on MPFF ICEs (Fig. 3). ICEs from the MPFG and MPFT classes were particularly promiscuous, being responsible for the spread of several CEGs of the metallo-beta-lactamase family such as bla NDM-1, bla SPM-1 and bla IMP variants among clinically relevant pathogens such as , and . ICEs of the MPFF class carrying blaIMP-6 or bla KPC-2 were restricted to and . The bla SPM-1 gene was exclusively identified in and in ICEs of the MPFT class (Table 1).

We analysed the types of integrase, relaxase and MPF classes present among four model ICEs: ICEclc, pKLC102, SXT and Tn4371. The MPFG-INT_P4_C ICEs identified here are related to ICEclc, since this ICE belongs to the same class and carries a MOBH relaxase and also an INT_P4_C integrase. The MPFT-INT_Rci_Hp1_C ICEs belong to the Tn4371 family, which also carries the MOBP relaxase. No conserved domain family could be attributed for SXT; however, the MPFF ICEs reported here should be related to this model ICE, which also uses a MOBH relaxase.

We also identified 386 hits encoding an integrase and a relaxase in the vicinity of CEGs (Table S7). For these hits, however, we could not predict if the CEG is located on a plasmid or an ICE, since the contig is fragmented and tracing the boundaries of the element is not possible, or the sequences have assembly gaps that make this prediction challenging. Some plasmids may also encode a tyrosine or serine recombinase, and some ICEs may encode replicases and partition systems that are typical of plasmids [34], which can hinder the accurate prediction of genetic platforms when the sequence has poor quality or is highly fragmented due to short-read sequencing approaches.

A variable repertoire of CEG-bearing integrons and transposons target ICEs

We identified 17 CEG variants among the 66 putative ICEs, dominated by bla NDM-1 and bla SPM-1 (Table 1 and Fig. 2). Insertion sequences (ISs; e.g. ISCR3-like elements) were frequently linked to the acquisition of bla SPM-1 and bla NDM-1 [35, 36], while bla IMP, bla VIM, and bla GES were found on class I integrons frequently integrated into Tn3 family transposons [6]. The bla KPC-2 gene was typically found within Tn4401-like transposons, which are capable of conferring a high frequency of transposition [37]. The recently identified AFM-1 metallo-beta-lactamase (GenBank accession number MK143105.1) was identified here in two ICEs inserted in and genomes (Table 1). We also found bla NDM-1 genes in ICEs integrated into the genomes of a recently proposed species, which is spreading in hospital settings in Myanmar [38]. Besides CEGs, the ICEs identified in this study also harbour genes conferring resistance to other antibiotics, such as aminoglycosides, fluoroquinolones, macrolides and tetracyclines (Table S8), widening the spectrum of transmissible AR genes selectable by carbapenems due to linkage.

Acquisition of additional traits by ICEs, including competitive weapons such as bacteriocins and siderophores

Besides genes conferring AR, the CEG-bearing ICEs identified here harbour other cargo genes that may confer a selective advantage to the ICE host. We found DUF692 domains typical of bacteriocin producing genes among the six MPFG ICEs from strains and the N15-01092, 1334/14 and ST773 strains (Table S5). All these ICEs carry a bla NDM-1 gene and a similar copy of the bacteriocin-encoding gene. Curiously, the bacteriocin producing gene was not identified in the MPFT ICE from strain MY569. Additionally, we found the siderophore aerobactin operon within an ICE in strain E302 (Table S5). This operon is usually found in enterobacterial plasmids, and was also identified in a pathogenicity island in uropathogenic strain CFT073 [39]. We identified no CRISPR-Cas systems among the ICEs here identified. Nearly half of them (47.0%, n=31/66) carried complete or incomplete restriction–modification systems belonging to types II, III and IV (Table S9).

The majority of the proteins encoded within the 66 ICEs refer to replication, recombination, transcription and intracellular trafficking functions (Fig. 4). Several proteins, however, encoded for unknown functions (34.0%, n=1 615/4 754, Table S10), highlighting the lack of knowledge concerning the ICE proteome.

Fig. 4.

Grouped bar chart representing the incidence of each eggNOG function broken down by category (represented by different colours) in the CEG-bearing ICE sequences.

Discussion

We have set out to comprehensively identify the CEGs among all bacterial genomes deposited in the NCBI database and the CEG-bearing ICE sequences. Our study considerably expands our knowledge of the repertoire of CEGs-bearing ICEs. We uncovered 66 putative ICEs that may be involved in HGT of CEGs amongst bacterial genomes. To expand our predictions, we also used the CONJscan module of Macsyfinder to trace the MPF families likely to be involved in HGT. Our analysis on the co-occurrence of relaxases with MPF families (Tables 1 and S5) is in agreement with the combinations observed by Guglielmini and colleagues [11]; all ICEs belonging to the MPFG class carried a MOBH relaxase, and the MOBP and MOBF relaxases were linked to MPFT and MPFF, respectively. All of the MPFG class ICEs described here present a MPF class/relaxase/integrase profile that resembles that from ICEclc. Even though pKLC102 is also a representative of the MPFG class and carries a MOBH relaxase, it uses a DNA_BRE_C integrase instead of a INT_P4_C. The absence of CEG-carrying ICEs from the pKLC102 family has already been reported [6].

The scenario observed for the acquisition of the most important CEGs by ICEs (ISs for bla NDM-1 and bla SPM-1, class I integrons for bla IMP and bla VIM and Tn4401 for bla KPC-2) resembles that of plasmids [3] and provides additional support for the notion that the line separating these elements is blurred [13, 14]. We now show that besides plasmids, this promiscuous repertoire of integrons and transposons frequently targets ICEs of different MPF families. Even though CEGs might spread rapidly worldwide, local selection is likely required for them to reach fixation, as can be seen for the clonal expansion of ST277 harbouring bla SPM-1 [36]. Surprisingly, we noted that this gene was not detected beyond the same clonal lineage. Indeed, all hits were identified in ST277 strains from Brazil and within a MPFT ICE family, indicating that these ICEs are transferring vertically and/or horizontally within this STs. Understanding the limitations on HGT of this family of ICE may be translatable to other, more transferable, ICEs and could underpin a control strategy to prevent the spread of these elements in the future.

Although bla KPC and bla OXA were the most frequently identified CEGs within the analysed genomes (Table S2), we only found five ICEs carrying bla KPC-2 and two ICEs with bla OXA genes (Table S8). The bla KPC genes are mostly located in Tn4401 transposons that target plasmids, while the bla OXA genes are frequently associated with ISs that tend to target spp. chromosomes and plasmids [3]. These genes may take advantage of the copy number and the higher genetic plasticity of plasmids [34]. This plasticity may increase the rate at which novel mutations appear and the high copy number may amplify the effect due to the increased gene dosage [40].

In addition to AR, ICEs can be involved in other adaptive traits such as carbon source utilization, symbiosis, restriction–modification, and siderophore and bacteriocin synthesis [14]. Since bacteria commonly inhabit highly competitive environments, the production of specific secondary metabolites (such as bacteriocins and siderophores) may confer a selective advantage to the host [41]. We speculate that the presence of these metabolites within the ICEs here characterized may promote their stability by preferentially selecting for cells harbouring the ICE. CRISPR-Cas systems are rarely found on MGEs, and the type I-C systems carried within pKLC102-related ICEs are one of the few examples [9]. Since none of the ICEs described here belong to the pKLC102 family, the absence of CRISPR-Cas systems within our dataset was expected.

One caveat of our studies is that these results do not yet expose the complete set of CEG-bearing ICEs present in all bacteria. There is an inherent bias in the number of times we detect a particular CEG in certain bacterial genomes, as some are over-represented in the database compared to others. It is possible that certain observations will flatten out as more genomes are analysed. Fewer than 10 % of the bacterial genomes currently present in the NCBI database are complete. This is a major drawback, since the putative ICEs present in draft genomes tend to be fragmented due to the presence of repetitive regions that are not resolved using short-read sequencing. Further, the relaxase and the T4SS encoded by ICEs resemble those of plasmids [11, 13, 42]. Plus, it is possible that ICEs and plasmids have swapped conjugation modules throughout their evolutionary history [11]. We believe that a more thorough exploration of this issue, especially regarding the precise delimitation of ICEs, will be an important further step toward an improved understanding of the contribution of these elements to bacterial adaptation and the evolution of AR.

While we have chosen to focus on CEG-bearing genomes, our computational approach can be applied to trace other relevant AR genes and other cargo genes that may confer a selective advantage to the ICE host. Leveraging knowledge linking the accurate prediction of ICE sequences to the carriage of AR genes will not only improve our understanding of HGT, but may also uncover potential approaches to tackle the spread of AR.

Click here for additional data file.