Molecular Characterization of Carbapenem Resistant Klebsiella pneumoniae and Klebsiella quasipneumoniae Isolated from Lebanon

Klebsiella pneumoniae is a Gram-negative organism and a major public health threat. In this study, we used whole-genome sequences to characterize 32 carbapenem-resistant K. pneumoniae (CRKP) and two carbapenem-resistant K. quasipneumoniae (CRKQ). Antimicrobial resistance was assessed using disk diffusion and E-test, while virulence was assessed in silico. The capsule type was determined by sequencing the wzi gene. The plasmid diversity was assessed by PCR-based replicon typing to detect the plasmid incompatibility (Inc) groups. The genetic relatedness was determined by multilocus sequence typing, pan-genome, and recombination analysis. All of the isolates were resistant to ertapenem together with imipenem and/or meropenem. Phenotypic resistance was due to blaOXA-48, blaNDM-1, blaNDM-7, or the coupling of ESBLs and outer membrane porin modifications. This is the first comprehensive study reporting on the WGS of CRKP and the first detection of CRKQ in the region. The presence and dissemination of CRKP and CRKQ, with some additionally having characteristics of hypervirulent clones such as the hypermucoviscous phenotype and the capsular type K2, are particularly concerning. Additionally, mining the completely sequenced K. pneumoniae genomes revealed the key roles of mobile genetic elements in the spread of antibiotic resistance and in understanding the epidemiology of these clinically significant pathogens.

The incompatibility (Inc) group FIIK is the most common plasmid replicon found in K. pneumoniae 10 . Other common plasmid replicons include the IncR and IncL/M 11 . Numerous antimicrobial resistance genes (ARGs) co-exist on the same plasmid and are transferred together, leading to treatment failure in cases of severe infections 11,12 .
In this study, we conducted an in-depth molecular characterization of 24 MDR, eight extensively-drug resistant (XDR) K. pneumoniae, and two MDR K. quasipneumoniae clinical isolates. These included all of the CRKP isolates subsampled from a larger dataset of K. pneumoniae collected from 2012 to 2015 during routine clinical microbiology screening in a 350-bed tertiary care hospital in Lebanon. We used whole-genome sequences (WGS) to investigate and compare their genetic relatedness, antimicrobial resistance, and the presence of virulence determinants. The genomes were used to build a genomic landscape and were placed into a comprehensive context through a comparison with isolates recovered from various clinical sources worldwide. The outcome of this study also provided new insights into the pathogenesis and resistance of the previously under-recognized human pathogen K. quasipneumoniae in Lebanon.

Results
Isolates characteristics. Isolates were subsampled from a pool of K. pneumoniae collected between 2012-2015 as part of routine clinical microbiology screening of patients admitted to the American University of Beirut Medical Center (AUBMC). The inclusion criterion was non-susceptibility, defined as resistant or intermediate phenotypes to at least one of the three clinically tested carbapenems, namely imipenem, ertapenem, or meropenem following the Clinical and Laboratory Standards Institute (CLSI) guidelines 13 . All of the isolates were non-susceptible to ertapenem, 85.3% (29/34) to imipenem, and 70.6% (24/34) Table S1a,b).

to meropenem (Supplementary
Resistance to meropenem was observed only when the isolate was simultaneously resistant to the remaining two carbapenems. In all cases where an isolate was resistant to only two carbapenems, these were ertapenem and imipenem. Regarding colistin, a last resort antibiotic used in the treatment of CRKP, 32.3% (11/34) of the isolates showed intermediate resistance according to the cut-off values proposed by Galani et al. (2008) 14 . 23% (8/34) of the isolates which were non-susceptible to at least one agent in all but two or fewer antimicrobial classes were classified as XDR (Supplementary Table S2) 15 .
PCR based ARG detection and ResFinder applied to the WGS revealed the presence of 54 distinct genes, each conferring resistance to one of the nine antibacterial classes: aminoglycosides, sulfonamides, quinolones, β-lactams, tetracyclines, fosfomycin, trimethoprim, chloramphenicol, and macrolides ( Fig. 1). In instances where a carbapenemase gene was absent, core chromosomal genes involved in carbapenem resistance, such as genes encoding for outer membrane porin proteins OmpK35 and OmpK36, were analyzed. Carbapenem resistance patterns detected through antimicrobial susceptibility testing were in agreement either with the detection of various carbapenemase genes or with modified chromosomal genes linked to carbapenem resistance.
In silico whole-genome and PCR-based analysis were used to assess the prevalence of ARGs. bla OXA-48 (61.2%; 21/34) was the most common carbapenemase detected, followed by bla NDM-1 (14.7%; 5/34) and bla NDM-7 (2.9%; 1/34). bla CTX-M-15 (41.2%; 14/34) was the most common ESBL variant detected. bla OKP-B-3 , a chromosomally-encoded β-lactamase, was detected in the K. quasipneumoniae isolates. Carbapenemase-encoding genes were absent in seven isolates (KP2, KP7-9, KP21, KP27, and KP32) showing phenotypic resistance to carbapenems associated with the presence of an ESBL coupled with porin loss or modification 6 . KP2, KP7-9 all carried multiple β-lactamases and mutated versions of the outer membrane porins OmpK35 and OmpK36 with at least seven and 13 mutations found in each, respectively. KP21 carried a bla SHV-11 ESBL and a truncated OmpK36 porin due to a premature stop codon in the ompK36 gene. This variation was also observed in KP27 and KP32. Since premature stop codons were only detected in ompK36, its non-synonymous changes within all isolates were assessed in an attempt to detect positive selection exerted on these mutations. Hence, SNAP (http://www.hiv.lanl.gov) was used to calculate synonymous and non-synonymous substitution rates based on Figure 1. Sequence types, antibiotic resistance genes, and Inc groups detected by PCR and in silico. ST, Sequence type. Classes of antibiotic resistance genes are marked as follows: A, aminoglycoside resistance genes; S, sulfonamide resistance genes; Q, quinolone resistance genes; B, β-lactam resistance genes; T, tetracycline resistance genes; F, fosfomycin resistance genes; R, trimethoprim resistance genes; C, chloramphenicol resistance genes; M, macrolide resistance genes; Inc Groups, plasmid incompatibility groups. the codon-aligned nucleotide sequences of an intact ompK36 sequence (accession # KY086538.1) against two other references (accession # KY086537.1 and MG577056.1) obtained from NCBI and the 34 ompK36 sequences retrieved from the isolates' WGS. Upon performing the analysis using SNAP, the number of observed and potential synonymous and non-synonymous substitutions were calculated and a plot was generated ( Supplementary  Fig. S1). These were then used to calculate the ratio of non-synonymous to synonymous substitutions (dN/dS) which was equal to 1.2589.
Plasmid studies. PCR-based replicon typing (PBRT) was performed and PlasmidFinder was applied to the WGS to successfully identify 13 Inc groups (Fig. 1). Using these two complementary methods, groups that were not covered by the PBRT kit, such as IncX4, were detected by PlasmidFinder, while many other missed plasmids through in silico analysis of draft genomes were confirmed by PBRT. At least three different Inc groups were identified in the majority of the tested isolates (91.2%; 31/34) (Fig. 1). The highest number of Inc groups (n = 7) was detected in KP1. The PLAsmid Constellation NETwork (PLACNETw) was used to reconstruct the plasmids from paired-end reads and it showed the presence of one big multi-replicon plasmid, carrying IncFIIK, IncFIA, and IncR, and three relatively smaller plasmids having IncA/C, IncL/M, and IncX1 replicons ( Supplementary Fig. S2). The IncP replicon was the only one missing from the short reads.
The IncX group of plasmids is relatively poorly studied compared to other Inc types 16 . The IncX4 plasmid, which carried a bla CTX-M-14b gene and detected for the first time in Lebanon, was only found in KP11. KP28 was positive for an IncX3 plasmid, and was the only isolate positive for bla NDM-7. To our knowledge, this is the first report of IncX3 carrying bla NDM-7 in Lebanon. Upon retrieving the complete plasmid sequence from the KP28 genome and BLASTing against GenBank, the IncX3 plasmid was found to be highly identical to previously published plasmids with different NDM subtypes. For instance, it was 99% similar to bla NDM-7 harboring pKpN01-NDM7 from K. pneumoniae (CP012990.1) isolated in Canada; 99% similar to bla NDM-7 encoding pM110_X3 from Escherichia coli (accession # AP018141.1) isolated in Myanmar; 97% similar to bla NDM-7 harboring pOM26-1 from E. coli (accession # KP776609.1) in the Arabian Peninsula, and 99% similar to pCRCB-101_1 (accession # NZ_CP024820.1) encoding bla NDM-5 found in Citrobacter freundii isolated in South Korea. Analysis of the genetic environment of IncX3 reflected the presence of only one copy of ∆ISAba125 in the KP28 plasmid (Fig. 2).
The wzi capsular gene was sequenced to assign a capsular type (K-type) to individual isolates. Based on that, thirteen different K-types were identified. 14.7% (5/34) of the isolates were typed as K15/K17/K50/K51/K52, while 26.5% (9/34) were not linked to any K-type. Due to the low resolution provided by sequencing a single wzi gene, the complete cps locus was targeted instead based on the WGS data using Kaptive. This resulted in the identification of 20 different K-types which were subsequently adopted throughout the manuscript 18 . No significant correlation was found upon comparing the STs, K-types, β-lactamases, incompatibility groups, and virulence features, highlighting the level of diversity among the isolates (Table 1).

Genome Statistics and Comparative Analysis.
On average, the 34 genomes consisted of 5,516,119 ±134,003 bp in 92 ± 16.98 contigs and had an average GC content of 57% ± 0.22% with a 30x reads coverage ( Supplementary Fig. S3). 594 ± 3.95 subsystems, 5,221 ± 178 coding sequences and 113 ± 7.22 RNAs were annotated on average by isolate (Supplementary Table S5). Comparative genome analysis using BRIG and KP1 (ST147) as the reference revealed deletion events involving membrane transport genes (e.g. ABC Transporter system, RND efflux system), fimbrial and type VI secretion system encoding genes. The type VI secretion system proteins were detected in six isolates. The circular genome based comparison may be overlooking key acquisition/ loss events in the remaining strains that are not present in the used reference genome (Fig. 5). Differences in the phage content of the isolates were observed upon using PHASTER annotation tool. In total, 43 different phages were identified. The number of detected phages ranged from one in KP14 and KP29 to 11 phages in KP18 with an average of 4.27 ± 2.23 phages per isolate. Salmonella phage RE-2010 (accession # NC_019488) was the most common phage, being identified in 38.2% (n = 13) of the isolates (Supplementary Table S6).
Pan-genome analysis revealed a total of 11,375 genes which were separated according to their prevalence within the isolates according to different categories including core (genes shared by 99-100% of the isolates), soft (genes shared by 95-99% of the isolates), and accessory genomes (genes shared by less than 95% of the isolates) ( Supplementary Fig. S4). 4,005 core and 7,214 unique accessory genes were identified. Isolates sharing the same accessory genes were clustered together (Fig. 6). All nodes had above 95% bootstrap support. A considerable correlation was not observed between sites of infection, age of patient, and overlap of resistance phenotypes. The phylogenetic analysis revealed close association between some of the ST-types (ST15, ST16, ST1799, and ST147) and their K-types, with the exception of one ST15 isolate (KP32) that clustered separately from the remaining ST15 isolates. This observed ST paraphyly is a direct implication of the robustness and value of whole-genome-based typing approaches over ST analysis, which is based only on seven housekeeping genes.
As an added layer of depth in the analysis, assessment of recombination events between the core genomes segregated the isolates based on predicted regional hotspots and unique single nucleotide polymorphisms (SNPs) denoted in Fig. 7 as red and blue vertical lines, respectively. 99,852 core genome SNPs, defined as base variations within the core genomes of the isolates and the eight reference sequences, were detected among the isolates and reference genomes and their alignment was used to generate a maximum-likelihood phylogenetic tree. This high resolution provided by WGS and SNP analysis has made such approaches the ultimate genotyping method currently available.
Using BLASTn search against the non-redundant nucleotide database we were able to associate the highly recombinant regions to known proteins such as: the P-type HAD ATPase IC family protein, 4-hydroxy-tetrahydrodipicolinate synthase, Glutamate/Aspartate Symporter, HipAB, and a protein of the Major Facilitator Superfamily.

Discussion
In this study, we performed an in-depth molecular characterization of 24 MDR, eight XDR K. pneumoniae, and two MDR K. quasipneumoniae clinical isolates. Using a whole-genome-based approach, we compared the genetic relatedness, antimicrobial resistance, and the virulome of the isolates under study along with other international clones. Significant variability was observed in the resistance and plasmid Inc groups among the same and different ST and K-types. None of the ARGs could be associated with a specific Inc group, except for bla OXA-48 , which was confined to IncL/M, and bla NDM-1, to the IncF replicons. Several virulence determinants were detected, with the yersiniabactin siderphore locus being notably carried either on ICEKp or IncFIB plasmid.
One of the key factors sustaining the successful global spread of carbapenemase encoding genes in Enterobacteriaceae is the transferability of plasmids 19 . The increase in carbapenem resistance is largely driven by self-conjugative plasmids 2 . Such plasmids carry multiple or single resistance determinants, such as bla NDM-1 and bla OXA-48 11 . OXA-48 (61.8%; 21/34) and NDM (17.6%; 6/34) were the only two detected carbapenemases. Many other studies from Lebanon have reported that NDM-1 and OXA-48 are the only circulating carbapenemases in the country, but these completely focused on the detection of individual genes without looking into the core and pan-genome, mobilome, virulome, and resistome [20][21][22][23] . bla OXA-48 was detected in isolates belonging to different clones (15 different STs and 15 different K-types).
The highest number of detected Inc groups in one isolate was seven (KP1), presented as one big multi-replicon plasmid (IncFIIK, IncFIA, and IncR) and three relatively smaller plasmids each having a single replicon (IncA/C or IncX1, or IncL/M) as determined by PLACNETw. The IncP replicon was missed from the WGS data. Multi-replicon plasmids have been increasingly detected with IncF plasmids being the most common type showing this feature, and plasmids with three distinct replicons, IncFIIK, IncR, and IncFII, were previously described by Coelho et al. 19 . Interestingly, the multi-replicon plasmid detected in this study carried a bla NDM-1 gene. Since plasmids might occasionally be replaced or lost by an incoming incompatible plasmid, a multi-replicon status allows the acquisition of plasmids having incompatible replicons by driving its replication by a compatible one 12 . This entails a broader host range replication and is certainly alarming as it can further facilitate widespread dissemination 12 .
The detection of a functionally intact bla CTX-M-14b gene in this study explains the resistance of KP11 to 2 nd , 3 rd , 4 th , and non-extended spectrum cephalosporins, with the exception of ceftazidime, a phenotype that was not seen in the remaining ST45 isolates.
To date 20 NDM variants were reported, with bla NDM-7 being recently described 24 . This variant differed from NDM-1 by two amino acid substitutions at positions Asp-130-Asn and Met-154-Leu 25 . bla NDM-7 shows a higher hydrolysis efficiency against carbapenems compared to its ancestor bla NDM-1 , which has a direct clinical and therapeutic implication, and often co-exists with other antibiotic resistance genes on the same IncX3 type conjugative plasmid 25 . IncX plasmids usually harbor either drug efflux system encoding genes or genes involved in resistance such as quinolone resistance, carbapenem resistance or extended-spectrum β-lactamases 8,26,27 . Unlike IncX1 and IncX4, the analysis of the complete genetic environment of bla NDM-7 on IncX3 in KP28 compared to a set of publicly available IncX3 plasmids revealed a typical backbone of IncX plasmids consisting of genes encoding replication (repB), partitioning (parB), maintenance (topB), and conjugal transfer (pil and tax) (Fig. 6) 8 . bla NDM-7 was carried by an IS5-∆ISAba125-bla NDM-7 -ble MBL -trpF-dsbC-IS26 genetic element. A divergence in the genetic environment mobilizing the bla NDM-7 in KP28 IncX3 was highlighted through the absence of a second ∆ISAba125 upstream of bla NDM-7 , as previously evidenced in pKpN01-NDM7 plasmid 28 . Earlier studies infer that IncX3 pOM26-1 plasmid might have evolved from pKpN01-NDM7 plasmid by the loss of one ∆ISAba125 element flanking IS5, which could also be the case for the identified plasmid in this study.
The genetic environment of the bla NDM-1 was also investigated, which in all isolates was carried on an IncF replicon in a chimeric structure with ble MBL encoding bleomycin resistance 29 . The bla NDM -ble MBL operon along with few neighboring genes was first formed in Acinetobacter spp. before being acquired by Enterobacteriaceae 30 . It is usually present on Tn125 flanked by two ∆ISAba125 29 . The downstream ∆ISAba125 was detected in all the isolates, whereas the upstream ∆ISAba125 was absent due to Tn3 insertion. This was also previously observed in K. pneumoniae and other Enterobacteriaceae 31 . Additionally, bla NDM-1 detected in KP18 was found on a multi-replicon plasmid possessing IncFIIK, IncFIB, and IncR, which were the only detected replicons. This  suggested the presence of a single big plasmid carrying a considerable number of ARGs. IncR replicons detected in K. pneumoniae are commonly associated with a variety of ARGs such as bla CTX-M-15 , bla TEM-1 , bla DHA-1 , bla OXA-1 , qnrS1, aac(6')-Ib, and qnrB4, all of which were found in KP18 2 . KP18 was typed as ST11, a SLV of ST258. ST258 is the predominant clone of KPC-producing K. pneumoniae and was responsible for its global dissemination 32 . ST11, sharing 80% of ST258 core genome, was also involved in the global spread of CTX-M-15 during the rise of ESBL-producing K. pneumoniae, and in recent years has acquired the bla NDM-1 33 . Accordingly, ST11 harboring bla NDM-1 is of a great concern since it has a superior epidemic potential and, once put under favorable conditions, can undergo quick clonal expansion and contribute to the successful global spread of bla NDM-1 7 . Loss or alterations of major porins OmpK35 and OmpK36 coupled with ESBL production is also associated with carbapenem resistance 34 . The underlying mechanisms of resistance involve either alterations through point mutations that could lead to protein modifications or premature truncation, or through insertion sequences that could lead to porin gene interruption. In seven of the isolates showing carbapenem resistance, both porin encoding genes were altered through point mutations, with three causing the premature truncation of the OmpK36. Premature stop codons were only detected in the ompK36 gene. Of note, all of these isolates were resistant to ertapenem, none showed resistance to meropenem and only one was resistant to imipenem. A similar finding was previously reported with carbapenem resistance involving porins and ESBLs being phenotypically observed against ertapenem only 35 . The presence of positive selection on these non-synonymous mutations was assessed using SNAP. Upon calculating the dN/dS ratio, it was shown that these mutations are positively selected for since the ratio was more than 1. Moreover, the overall number of non-synonymous mutations were shown to be much higher than that of synonymous mutations (Supplementary Fig. S1). This finding was expected since changes in porin proteins coupled with an ESBL or AmpC type of enzymes confer a survival advantage. Positive selection involving outer membrane porin proteins was also previously detected in E. coli, Shigella flexneri, Neisseria gonorrhea and N. meningitidis 36,37 . The yersiniabactin siderophore, mobilized by the integrative and conjugative element ICEKp, was one of the most common VFs found in almost half of the isolates in this study. The yersiniabactin locus was traced on diverse ICEKp types with some being localized on a plasmid highly similar to K. pneumoniae INF167_p0001 plasmid (Supplementary Table S4) 17 . Unlike recent reports, none of the ICEKp elements harbored a colibactin  or salmochelin synthesis loci 17 . The observed diversity of the detected ICEKp elements suggested its acquisition through horizontal gene transfer rather than stable vertical inheritance 17 . The presence of various iron-scavenging mechanisms can give rise to increased efficiency of iron uptake, which has been correlated with hypervirulent phenotypes 3 . Interestingly, KP25 carried a unique repertoire of VFs that included kfu, all, glx, ybb, and ylbE. The identification of kfu and all in KP25 with the capsular type K151 is a novel finding where it was previously only reported in strains having a K1 capsular type 38,39 . all is a relatively rare virulence gene found in organisms metabolizing allantoin to obtain carbon and nitrogen from the environment 4 . One out of the three NDM-1 producing XDR isolates, KP16, remained susceptible to only one drug class. The eight detected XDR isolates (ST15, ST16, ST147 and ST995) harbored on the average two VFs, with the exception of KP16 that carried none. The convergence of ARGs and VFs is a particularly alarming phenomenon for the possible emergence of untreatable invasive K. pneumoniae infections.
Maximum-likelihood clustering based on core genome SNPs and recombination hotspots was better at showing the evolutionary relatedness of the isolates; the branch length indicates the number of substitutions that an isolate has undergone relative to a common ancestor.
Analyzing the recombination hotspots of the isolates revealed several sequences involved in transport and resistance. One hotspot was found to be allocated within the P-type HAD ATPase IC family sequence. The P-type HAD ATPase IC is a transmembrane pump for ion and lipid transportation linked to phospholipids modification; recently gaining attention as unique drug targets in fungi 40 . Similarly, 4-hydroxy-tetrahydrodipicolinate synthase (DHDPS) was another hotspot identified, which previously was identified as a potential drug target that can be used against Mycobacterium tuberculosis 41 . Recombination events in the hipA/B gene and the genes encoding the proteins of the Major Facilitator Superfamily detected in this study are also of prime importance due to them being involved in bacterial resistance to fluoroquinolones, β-lactams, and quaternary ammonium compounds 42,43 .
To our knowledge, this is the first comprehensive whole-genome-based characterization of carbapenem resistant K. pneumoniae in Lebanon, and the first detection and description of K. quasipneumoniae in the region. Few studies have addressed in details the genomic content of K. quasipneumoniae isolated from a clinical setting. These studies tackled hypervirulence in community-acquired K. quasipneumoniae isolates [44][45][46] . Even fewer studies have described MDR K. quasipneumoniae 47,48 . We were the first to identify and report the presence of carbapenem resistant K. quasipneumoniae in the Middle East. The multi-replicon nature of the detected plasmids, the observed diversity of STs and K-types, and the coexistence of specific VFs with a very high number of ARGs all indicate the possible emergence of a population of untreatable invasive K. pneumoniae. This is a problem that is too important to be ignored. WGS allowed us to explore in-depth the genetic variation within the K. pneumoniae population and the genetic relatedness of the highly resistant isolates. Within clinical settings of Lebanese hospitals, traditional methods of identification and typing, such as API20E, disk diffusion, or PFGE, are still commonly used. When compared to WGS, these methods are known to suffer from low resolution and discriminatory power. In this study, a genome-based workflow was coupled with these traditional methods to overcome their limitations and culminated in an accurate population snapshot of carbapenem resistant K. pneumoniae isolated from hospitalized patients in Lebanon. The obtained data can be exploited to limit the cross-transmission of such pathogens in hospital settings, prevent outbreaks, and global dissemination.

Materials and Methods
Ethical Approval. Ethical approval was not required as clinical isolates were collected and stored as part of routine clinical care. Clinical isolates and patient records/information were anonymous and de-identified prior to analysis. Bacterial isolates. A total of 32 K. pneumoniae and two K. quasipneumoniae isolates collected from 2012 to 2015 were recovered at the Clinical Microbiology Laboratory at AUBMC and were designated as KP1-KP32 and KQ1-KQ2, respectively. AUBMC is a 350-bed tertiary care teaching hospital. The isolates were collected from various body sites such as urine, deep throat, skin wounds, including one from a C-section, among others (Supplementary Table S2). The mean patients' age was 58 ± 21 years old, with a range of 18 to 87 years with a 1:1 male to female ratio. Initial species identification was performed based on the API 20E system (bioMérieux, Marcy l'Etoile, France) and by the Matrix-Assisted Laser Desorption/Ionization Time of Flight (MALDI-TOF) system (Bruker Daltonik, GmbH, Bremen, Germany) following the manufacturer's instructions.
String test. The string test, a phenotypic screening test to assess the hypermucoviscous phenotype, was performed on fresh bacterial colonies, as previously described 3 . Colonies were considered positive if the viscous string was at least five mm in length as determined through stretching bacterial colonies on an agar plate using an inoculation loop.
Antimicrobial testing. Antimicrobial susceptibility testing was performed by the disk diffusion technique on Mueller-Hinton agar and it included a panel of 26 antibiotic disks belonging to 15 classes (Supplementary  Table S1). Additionally, carbapenems were tested by the E-test methodology (AB BIODISK, Solna, Sweden) to determine the minimal inhibitory concentrations (MICs) of ertapenem, imipenem, and meropenem. The obtained data was interpreted according to the CLSI guidelines (CLSI, 2017) and Galani et al. 14  Detection of resistance genes by PCR assays. The most common ESBLs and carbapenemases encoding genes, bla CTX-M , bla SHV, bla IMP , bla VIM , bla KPC , bla NDM-1 , and bla OXA-48-like , were amplified as previously described (Supplementary Table S7).
Plasmid replicon typing. Plasmid characterization was performed using the DIATHEVA PBRT kit (Diatheva, Fano, Italy) through a PCR-based replicon typing method consisting of eight multiplex PCR assays for the amplification of 25 replicons: A/C, B/O, FIA, FIB, FIB-M, FIC, FII, FIIK, FIIS, HI1, HI2, HIB-M, I1, I2, K,  L/M, N, P, R, T, U, W, X1, X2, and Y found in the family Enterobacteriaceae. Positive controls were included for all reactions. All PCR reactions were performed according to the manufacturer's instructions and visualized on a 2.5% agarose gel stained with ethidium bromide.
PCR amplification and sequencing of wzi gene. To determine the capsular type, wzi gene typing was performed using wzi-TR and wzi-TF primers as previously described (Supplementary Table S7). K-types were assigned using the Institute Pasteur database (http://bigsdb.pasteur.fr/klebsiella). The isolates that were associated with multiple or no K-types by wzi sequencing were re-analyzed using the publicly available Kaptive online tool (https://github.com/katholt/Kaptive) 18 . Kaptive makes use of WGS data as input and assigns capsular types by analyzing the entire cps locus.
Multilocus sequence typing (MLST). MLST was performed as described on the Institute Pasteur MLST database targeting seven housekeeping genes (rpoB, gapA, mdh, pgi, phoE, infB, and tonB) using primers with universal sequencing tails. Genes were sequenced using the universal oF and oR primer pair. STs were assigned using the Institute Pasteur database. Novel STs were submitted to the curator and assigned as new designations (www.pasteur.fr/mlst). For isolates showing identical pulsotypes or that were untypable by XbaI, PFGE was repeated using the secondary enzyme AvrI. PFGE profiles were analyzed with the BioNumerics software version 7.6.1 (Applied Maths, Sint-Martens-Latem, Belgium), with profiles assigned as different pulsotypes if three or more bands were different between the two of them. Pulsotypes were clustered based on the BioNumerics software analysis through dice correlation coefficients with an optimization of 1.5% and tolerance of 1.5%.

Pulsed-field gel electrophoresis (PFGE
Whole-genome sequencing. Genomic libraries were constructed using the Nextera XT DNA library preparation kit with dual indexing (Illumina). The libraries were sequenced on an Illumina MiSeq with 250 bp x 2 read length. Genome assembly was performed de novo using Spades Genome Assembler Version 3.6.0 49 . Quality control checks on the obtained raw sequence data was performed using FastQC version 0.11.5 50 .
Genome assembly, annotation and analysis. The assembled genomes were annotated using the RAST online server (http://rast.nmpdr.org) 51 . The Comprehensive Antibiotic Resistance Database (CARD) and ResFinder 3.0 available on the Center for Genomic Epidemiology website (www.genomicepidemiology.org) were used to determine the presence of resistance genes 52,53 . The presence of putative VFs was screened using the VF scheme available on http://bigsdb.pasteur.fr. PlasmidFinder 1.3 was used to determine the presence of plasmids in the genomic sequences 54 . PLACNETw was used to identify and analyze plasmids starting from raw reads by creating a network of contig interactions upon assembly of reads with Velvet 55 . The sequences encoding plasmids were then extracted and aligned to references obtained from NCBI. PHASTER was used to identify the phage content of the isolates 56 . SNAP (http://www.hiv.lanl.gov) was used to calculate synonymous and non-synonymous substitution rates of the ompK36 porin gene 57 .
Comparative genome analysis. Plasmid sequences were extracted and aligned with corresponding reference strains using BioNumerics software version 7.6.1. A circular genome comparison was performed using BRIG 58 .