We retrieved all public KP genome assemblies from BV-BRC database as of September 4, 2024 (Olson et al., 2023). The search criteria used were: “GENOMES = Klebsiella pneumoniae”, “isolation country = China”, “genome status = WGS”, “genome quality = good” and “host group = human”. The corresponding metadata for these genomes was also obtained from the BV-BRC database and cross-verified manually using the NCBI GenBank database. Quality assessments for downloaded genomes were calculated with QUAST v5.2.0 and fastANI v1.33, respectively (Gurevich et al., 2013; Musiał et al., 2024). Genomes with over 95% average nucleotide identity and at least 80% coverage relative to the reference genome HS11286 (NCBI RefSeq: NC_016845) were selected for further analysis.

Kleborate v3.0.8 identified MLST profiles, virulence loci such as ybt, clb, iro, iuc, and rmpA/A2, predicted capsule types (serotype), and detected acquired antimicrobial resistance (AMR) genes, including blaESBLs and blaCarbs (Lam et al., 2021). We defined isolates carrying at least one of blaCarbs as CRKP. Isolates possessing both the iucA and rmpA/A2 genes were categorized as hvKP (Hu et al., 2024). Isolates that met both criteria were defined as CRhvKP (Chen et al., 2024).

The core genome alignment and SNP calling (cgSNP) were conducted using Parsnp v1.2 from the HarvestTools suite (Treangen et al., 2014). A maximum likelihood phylogenetic tree based on the concatenated alignment was inferred using FastTree v2.1 with the GTR+GAMMA model (Price et al., 2010). The resulting phylogenetic tree with metadata was visualized using the interactive Tree of Life (iTOL) web application (Letunic and Bork, 2021).

For cgMLST allele calling, the chewBBACA v3.3.9 was employed with a public 2,358-gene typing scheme derived from 14,254 genomes of K. pneumoniae, K. variicola, and K. quasipneumoniae, as provided by the RIDOM Nomenclature Server (https://www.cgmlst.org) (Silva et al., 2018). Pairwise cgMLST distances between strains were calculated using cgmlst-dists v0.4.0 (https://github.com/tseemann/cgmlst-dists), based on the core genes present in more than 95% of the collected genomes. A minimum spanning tree (MST) was constructed using GrapeTree v1.5.0 with the MSTv2 algorithm (Zhou et al., 2018). The cgMLST-based genomic clustering network was visualized and analyzed with Cytoscape v3.10.2 (Shannon et al., 2003).

Plasmid sequences were identified and reconstructed from the genome assemblies with MOB-Suite v3.19 (Robertson and Nash, 2018). This process involved using MOB-typer for relaxase and replicon typing, generating MOB-cluster codes, and determining host range. Plasmid mobility was predicted based on relaxase (MOB), mating pair formation (MPF) complex, and oriT gene annotations. In summary, plasmids were classified as conjugative if they had both relaxase and MPF, mobilizable if they contained either relaxase or oriT but lacked MPF, and non-mobilizable if they had neither relaxase nor oriT (Smillie et al., 2010).

Pairwise distances between plasmids were calculated using Mash v2.2, which estimates sequence similarity by converting sequences into fixed-length MinHash sketches (Ondov et al., 2016). The Mash distance (ranging from 0 for nearly identical sequences to 1 for completely dissimilar sequences) was used to assess similarity (1-Mash distance), with a k-mer length of 13 and a sketch size of 5,000.

For plasmid community detection, the Louvain algorithm was employed to identify clusters by optimizing modularity through iterative expectation-maximization, implemented via the Python module python-Louvain v0.16 (Blondel et al., 2008). We calculated the Louvain partition for the network and selected nodes from communities containing at least 10 members.

Plasmids were annotated using Prokka v1.13.4 (Seemann, 2014). Insertion sequences (ISs) were detected by digIS v1.2 (Puterová and Martínek, 2021). Pangenome analysis was performed by Panaroo v1.5.0 (Tonkin-Hill et al., 2020). Abricate v1.0.1 (https://github.com/tseemann/abricate) was used to identify acquired AMR genes from the CARD database and virulence factors (VFs) from the VFDB database, using a threshold of 90% sequence identity and 90% coverage (Liu et al., 2022; Alcock et al., 2023). Circular alignment of the plasmids was performed and visualized with BRIG v0.95 (Alikhan et al., 2011).

Annotations for integron integrases, attC/attI sites, and integron promoters (P_c/P_int) were generated using IntegronFinder v2.0.5 with default parameters, except for the options “–local-max” (Néron et al., 2022). IntegronFinder classifies integron elements into three types: (1) Complete integrons, which include an integrase and at least one attC site; (2) In0 elements, consisting of an integrase without any nearby attC site; and (3) CALIN elements, clusters of attC sites that lack associated integron integrases, comprising at least two attC sites. The AMR genes located on the integron-associated gene cassettes were also detected by Abricate v1.0.1 based on the CARD database.

We collected 5,036 KP genomes with metadata from China after screening a database on September 4, 2024. Genotyping with Kleborate identified 3,449 isolates as CRKP and 2,212 as hvKP, with 1,538 classified as CRhvKP (44.6% of CRKP and 69.5% of hvKP; Supplementary Table S1 ). These CRhvKP isolates were reported from 26 province-level divisions, primarily concentrated in the central and eastern regions ( Figure 1A ). The top five regions by isolate count were Zhejiang (n=629), Shanghai (n=288), Beijing (n=127), Jiangsu (n=77), and Henan (n=63). The earliest sequenced CRhvKP in China was isolated from Zhejiang in 2005. Over the past 20 years, CRhvKP isolates have significantly increased, peaking in 2021 before declining.

According to the virulence detection of Kleborate, 1,439 (93.6%) isolates had a virulence score of 4, indicating the presence of both aerobactin and yersiniabactin ( Figure 1B ). Additionally, 63 strains (4.0%) had the highest score of 5, carrying aerobactin, colibactin, and yersiniabactin, while only 36 strains (2.3%) scored 3, indicating the presence of aerobactin alone. Notably, since 2019, CRhvKP strains carrying only aerobactin have become increasingly rare.

bla _KPC were detected in 1,224 (92.1%) of collected CRhvKP isolates and were the most prominent blaCarb. bla _OXA-48-like was the second most frequent (n=273, 17.8%), followed by bla _NDM (n=53, 3.4%), bla _IMP (n=26, 1.7%) and bla _VIM (n=1, 0.1%). 39 strains carried multiple blaCarbs: 34 encoded bla _KPC and bla _NDM, 2 had bla _KPC and bla _OXA-48-like, 2 had bla _NDM and bla _OXA-48-like, and 1 had bla _KPC and bla_IMP. There has been a significant increase in the proportion of strains encoding bla _OXA-48-like since 2021 ( Supplementary Figure S1 ). Furthermore, 96.9% (n=1,491) of collected CRhvKP isolates contained resistance determinants for 3 or more classes of antibiotics, and the majority (n = 1,212, 78.8%) contained blaESBLs. Over the past five years, the proportion of CRhvKP strains co-harboring both blaCarbs and blaESBLs has risen annually, from 70.1% in 2017 to 92.3% in 2023 ( Figure 1C ).

To gain molecular insights, we reconstructed a phylogenetic tree for all 1,538 CRhvKP isolates, mapping it with ST, KL type, blaCarb class, blaESBL presence, collection year, and location ( Figure 2A , Supplementary Table S1 ). The analysis showed that strains clustered by STs and KL types, identifying 36 STs, with the top five being ST11 (68.7%), ST15 (21.1%), ST23 (1.9%), ST268 (1.5%), and ST65 (1.2%). Although ST11 was the most prevalent, it was genetically distant from ST15, ST23, and ST268. The proportion of ST15 increased over the past five years, reaching 50.8% in 2021 ( Figure 2B ). We identified 134 capsular types (KLs), with the top five being KL64 (54.4%), KL112 (15.6%), KL47 (13.3%), KL24 (5.1%), and KL1 (2.1%) ( Figure 2C ).

Notably, ST11 and ST15 had a virulence score of 4 (>99%), while ST23, ST268 and ST65 had a score of 5 in 96.6%, 87% and 50% of cases, respectively ( Figure 2D ). Besides, 100% of ST65, 93.1% of ST23 and 56.5% of ST268 strains encoded the iroBCDN virulence gene cluster. In contrast, only 6.2% of ST11 and 0.3% of ST15 strains carried these genes ( Supplementary Figure S2 ). Moreover, the number of acquired AMR classes varied significantly among the major STs ( Figure 2E ). The proportion of strains carrying more than six antibiotic resistance classes was highest in ST11 at 83.6%, followed by ST65 (61.1%), ST268 (60.9%), and ST15 (55.2%), while the proportion in ST23 was only 17.2%.

Using a threshold of 15 cgMLST allele differences (Lan et al., 2020), 1,369 isolates (89.0%) were grouped into 131 clonal clusters, ranging in size from 2 to 180 members, and ordered from largest to smallest, from C1 to C131 ( Figure 3A , Supplementary Table S2 ). Compared to the clustering rates of ST11 (91.2%) and ST15 (95.7%), ST23 had a lower clustering rate of 41.4%. Clustered samples of over 10 strains were detected in 9 province-level divisions, ranked by clustering rate: Guangdong (94.1%, 48/51), Zhejiang (92.2%, 599/650), Henan (90.5%, 57/63), Beijing (88.2%, 112/127), Anhui (88.1%, 37/42), Shanghai (87.8%, 253/288), Jiangsu (87.0%, 67/77), Sichuan (86.1%, 31/36) and Shandong (83.7%, 36/43), suggesting the severe CRhvKP epidemic in these areas.

Strains from 42 clusters were sampled from 20 province-level divisions, with two major clusters, C1 and C2, each containing over 100 members. C1 included 180 strains collected from eastern (Shanghai, Zhejiang, Jiangsu), northern (Beijing), central (Anhui, Hunan), and western (Sichuan, Yunnan) regions over 18 years (2005–2023). C2 comprised 164 strains from eastern (Shanghai, Zhejiang, Jiangsu, Fujian), northern (Beijing), central (Henan), and western (Sichuan) regions, collected over 7 years (2015–2022). A bipartite network analysis indicated that eastern China, especially Shanghai, Zhejiang, and Jiangsu, served as a transmission hub for CRhvKP, linking strains from 13 regions, excluding Hebei, Inner Mongolia, Chongqing and Jilin ( Figure 3B ). Notable connections also emerged from Beijing (linked to Hebei and Inner Mongolia) and Sichuan (linked to Chongqing).

A total of 8,645 plasmids were detected based on MOB-Suite, with the number of plasmids per strain ranging from 1 to 11. The number of plasmids in ST15 strains was significantly higher than in the other four main STs (Kruskal–Wallis test P < 0.01, Supplementary Figure S2 ). Plasmids ranged in length from 1.2 to 365.0 kb. Twelve types of replicon families were detected from these plasmids, with the IncF being the most prevalent replicon (41.3%, 3569/8645).

A total of 1,496 hv plasmids carrying hv genes were identified in 1,494 strains, and 1,546 CR plasmids encoding blaCarbs were found in 1,519 strains ( Supplementary Table S3 ). Additionally, 42 strains carried both hv and blaCarb genes on a single plasmid. Using the Louvain algorithm with a Mash distance threshold of 0.95 (Matlock et al., 2021), we identified 10 plasmid communities (PC1 to PC10), with cluster sizes ranging from 15 to 1,260 ( Figure 4 ). PC1 (n=1,254), PC4 (n=227) and PC5 (n=126) were broadly distributed among ST11, ST15, ST23 and other STs, whereas PC2 was predominantly found in ST11 (97.7%, 979/1002), and PC3 was primarily associated with ST15 (96.4%, 240/249). Plasmid mobility analysis showed that most PCs were conjugative or mobilizable, except for PC1 and PC9, which were largely non-mobilizable.

CR and hv plasmids were assigned to different PCs, with CR plasmids showing greater diversity across 9 PCs, where 64.8% (1002/1546) clustered in PC2. Conversely, 83.8% (1254/1496) of hv plasmids were predominantly in PC1, with an additional 15.1% (226/1496) in PC4. Ten different replicons were identified in both CR and hv plasmids, with IncF (74.5%, 2236/3002), IncH (43.5%, 1306/3002), IncR (31.0%, 930/3002), and Col-like (14.9%, 447/3002) being the most common. Specifically, PC1 and PC4 predominantly featured IncF/IncH, PC2 was mainly composed of IncF/IncR, and PC5 included a mix of Col-like/IncF/IncR. PC10 was characterized by IncF/IncU, while PC3, PC6, PC7, PC8 and PC9 were primarily represented by Col-like, IncL/M, IncN, IncU and IncU, respectively ( Supplementary Figure S3 ).

A total of 40 CR-hv convergent plasmids were identified, carrying both blaCarb (38 bla _KPC-2, 2 bla _NDM-1) and hv genes: 32 were classified as PC1, six as PC4, one as PC2, and one as other ( Supplementary Table S3 ). PC1 was found in ST11, ST23, ST268, ST660 and ST1660, comprising hybrid plasmids from seven replicon types (IncF, IncH, IncR, IncA/C, IncN, IncU and IncX). In contrast, PC2 and PC4 were only detected in ST11 and were solely IncF/IncR hybrids. Since their first identification in Zhejiang isolates in 2013, CR-hv convergent plasmids have been consistently detected through 2022, primarily in eastern China (85.0%, n=34). Strains carrying these plasmids belong to 13 clonal clusters, with C2 being the largest at 37.5% (n=15), indicating a crucial role of these plasmids in spreading both CR and hv traits.

To investigate the genetic features of CR-hv convergent plasmids, we selected two complete circular IncF/IncR plasmid sequences from PC1 and PC2 as references: CP119178 (226.6 kb) and CP095248 (206.0 kb). Sequence alignments revealed that both plasmids resulted from fusion events in small homologous regions. CP119178, a non-mobilizable plasmid, formed by fusing a 133.3 kb sequence from the IncF/IncH plasmid CP154274 (PC1) with a 93.6 kb sequence from the IncF/IncR plasmid CP119178 (PC2), linked by a 328 bp region encoding a hypothetical protein ( Figure 5A ). In contrast, CP095248, a conjugative plasmid, merged a 118.3 kb sequence from CP107298 (PC1) with an 88.0 kb sequence from CP099415 (PC2), connected by a 1,345 bp IS5075 (IS110-like) transposon and a 255 bp region encoding a hypothetical protein ( Figure 5B ). Additionally, the replicons (IncF and IncR) in both CP154274 and CP119178 were located on segments derived from the fused PC2 plasmid.

We identified hv and blaCarb genes on the chromosomes of several strains in addition to plasmids. A total of 28 strains carried hv genes (iucABCD, iutA and rmpA2) on their chromosomes, with 41 having them exclusively located there. Among these, 21 were ST11, 6 were ST2237, and one belonged to ST1544. Notably, a ST11 strain named XHKPN391, isolated from Shanghai in 2018, carried both bla _KPC-2 and hv genes on its chromosome (CP066915), with hv genes found only there. In the clonal cluster (C10) containing XHKPN391, none of the other 30 members had bla _KPC-2, although 15 strains carried hv genes.

In ST11 and ST2237, the hv genes are situated within a 40.8 kb region, bordered by ISEc36 (IS3-like) and IS1663 (IS110-like) ( Figure 6A ). In contrast, ST1544 features hv genes within a 35.6 kb segment, flanked by two IS903 transposons. BLAST results showed that the two chromosomal hv regions on ST11 and ST2237 closely matched the PC1 plasmid (99.9% identity, 99% coverage), while the chromosomal hv region on ST1544 was highly homologous to the PC4 plasmid (99.9% identity, 95% coverage), featuring IS903 upstream and IS1663 downstream. Additionally, comparative analysis revealed that these hv regions shared 30 core genes, which make up 52-83% of each region.

Two ST11 strains, XHKPN391 and SM117, carried blaCarbs on both their chromosomes and plasmids. XHKPN391 harbored two bla _KPC-2 genes on its chromosome (CP066915) and one on a plasmid, while SM117 contained one bla _KPC-2 on its chromosome (CP130662) and both a bla _KPC-2 and bla _NDM-5 on two plasmids. The bla _KPC-2 in both strains were located within genetic regions flanked by two IS26 (IS6-like) transposons, measuring 13 kb and 21 kb, respectively ( Figure 6B ). A BLAST analysis showed that the IS26-flanked bla _KPC-2 regions in XHKPN391 and SM117 shared high synteny and identity (>99%) with plasmids from ST11 isolates IR12079_1 (CP097691) and Kp337 (CP154492), respectively. The core genetic structure for bla _KPC-2 in these replicons was characterized by IS26, ISKpn6, bla _KPC-2, ISKpn27, tnpR, and IS26.

A total of 1,336 class 1 integrons were identified across 1,300 (84.5%) CRhvKP isolates, including 725 complete integrons, 448 In0, and 163 CALIN elements. Integrons were found in more than 75% of ST11, ST15, ST268, and ST65 isolates, while the detection rate in ST23 strains was significantly lower, at 44.8% ( Figure 7A ). The predominant integron type in ST11, ST15 and ST65 strains was the complete integron, while ST23 and ST268 were primarily characterized by In0 and CALIN elements, respectively ( Figure 7B ). We identified 14 AMR gene families from the mobile gene cassettes within these integrons, conferring resistance to various antibiotics and antiseptics ( Figure 7C ). Notably, four integrons in ST65 carried bla _IMP, and one in ST23 contained bla _VIM, underscoring their role in blaCarb dissemination. Furthermore, the integron-related AMR gene profiles varied significantly across different STs. For example, ST11 mainly carried ANT(3’’), dfr and qacEΔ1, while ST15 primarily harbored AAC(6’) and arr.