Kp MDR evolution is largely driven through acquisition of AMR genes on diverse mobilisable plasmids⁷ which are particularly prevalent among the subset of clones that have become globally disseminated and frequently cause hospital outbreaks⁸; e.g. clonal group (CG) 258 which is implicated in global spread of the Kp carbapenemases⁹. Kp pathogenicity is driven by a wide array of interacting factors^{10 –12} that are present in all strains, including the type III fimbriae (mrk) and the surface polysaccharides (capsule and lipopolysaccharide (LPS))^{12, 13} which exhibit antigenic variation between strains. The majority of hypervirulent Kp, distinguished clinically as causing invasive infections even outside the hospital setting¹⁴, are associated with just two^{4, 15} of the >130 predicted capsular serotypes¹⁶, K1 and K2, that are considered particularly antiphagocytic and serum resistant^{15, 17}. Hypervirulent Kp are also associated with high prevalence of several other key virulence factors; the rmpA/rmpA2 genes that upregulate capsule expression to generate hypermucoidy^{18, 19}; the colibactin genotoxin that induces eukaryotic cell death and promotes invasion to the blood from the intestines^{20, 21}; and the yersiniabactin, aerobactin and salmochelin siderophores that promote survival in the blood by enhancing iron sequestration^11,22–24.

Yersiniabactin synthesis is encoded by the ybt locus, which is usually mobilised by an integrative, conjugative element known as ICEKp. It is present in ∼40% of the general Kp population and seems to be frequently acquired and lost from MDR clones²⁵. Fourteen distinct ybt+ICEKp variants are recognised, one of which also carries the colibactin synthesis locus (clb)²⁵. In contrast, the salmochelin (iro), aerobactin (iuc) and rmpA/rmpA2 loci are usually co-located on a virulence plasmid^{26, 27}. These loci are much less common in the Kp population (<10% prevalence each) and until recently were rarely reported among MDR strains^33,5.

The reasons for the apparent separation of MDR and hypervirulence are unclear but there are growing reports of convergence from both directions, i.e. hypervirulent strains gaining MDR plasmids^28–33 and MDR strains gaining a virulence plasmid plus/minus an ICEKp^6,30,34. Most such reports are sporadic, but in 2017 Gu and colleagues described a fatal outbreak of MDR, carbapenem-resistant Kp belonging to CG258 that had acquired ICEKp in the chromosome plus iuc and rmpA2 on a virulence plasmid⁶. The report fuelled growing fears of an impending public health disaster in which highly virulent MDR strains may be able to spread in the community, causing dangerous infections that are extremely difficult to treat³⁵. However, there remain significant knowledge gaps about Kp evolution that limit our ability to understand the severity of this public health threat, and to evaluate the relative risks of convergence events.

Here we report a comparative analysis of genome evolutionary dynamics between MDR and hypervirulent clones, leveraging a curated collection of 2265 Kp genomes to identify 28 common clones with at least 10 genomes in each (total 1092 genomes in 28 clones, 10–266 genomes per clone, see Supplementary Tables 1,2 and Figure 1a,b). First we explored the distribution of key virulence loci and resistance determinants for 10 distinct antimicrobial classes across the total population (Methods). The presence of the virulence plasmid (with or without ICEKp) was negatively associated with the presence of acquired AMR genes: the majority (n=77/88) of genomes harbouring the virulence plasmid contained 0–1 acquired AMR genes, while the distributions were much broader for genomes without the virulence plasmid plus/minus ICEKp (median 1, interquartile range (IQR) 0-9, p < 1×10⁻¹⁰ for both pairwise Wilcoxon Rank Sum tests). Interestingly, among the genomes without the virulence plasmid the distribution of acquired AMR genes was slightly shifted towards higher numbers in genomes harbouring ICEKp (median 1 vs 1, IQR 0-8 vs 0-10, p < 1×10⁻⁸; see Supplementary Figure 1). Next we calculated the proportion of genomes harbouring the virulence and AMR loci within each of the 28 common clones (Figure 1c). Hierarchical clustering of the virulence locus data defined a group of 6 hypervirulent clones including the previously described CG23, CG86 and CG65³, each of which harboured the virulence plasmid-associated genes iuc, iro and/or rmpA/rmpA2 at high frequency (31– 100%). Hierarchical clustering of the AMR data defined a group of 8 MDR clones, including the outbreak associated CG258, CG15 and CG147⁸, each with a high frequency (≥56%) of genomes encoding acquired resistance determinants for ≥3 drug classes (in addition to ampicillin to which all Kp share intrinsic resistance via the chromosomally encoded SHV-1 beta-lactamase). As expected, AMR genes were rare among the hypervirulent clones, with the exception of CG25 in which 11 of 16 genomes harboured ≥4 acquired AMR genes. The iuc, iro and rmpA/rmpA2 loci were rare (<12% frequency) among the MDR clones and those not assigned to either group (‘unassigned’ clones); however, the ybt locus was frequently identified across the spectrum of clones as has been reported previously²⁵.

We used Gubbins³⁶ to identify putative chromosomal recombination imports within each clone and calculated r/m (the ratio of single nucleotide variants introduced by homologous recombination relative to those introduced by substitution mutations), which ranged from 0.02–25.50 (Figure 2a, Supplementary Table 1). With the exception of CG25, the hypervirulent clones generally exhibited lower r/m values (median 1.15), while the MDR clones trended towards higher values (median 5.47) although the differences were not statistically significant (Kruskall-Wallis test p = 0.07). Recombination events were not evenly distributed across chromosomes: in 19/28 clones ≥50% of the chromosome was not subject to any recombination events, while the maximum recombination load in each clone ranged from mean 1.1–47.6 events (Figure 2b, Supplementary Figure 2).

In many cases there was a major peak defining a recombination hot-spot at the capsule (K) and adjacent LPS antigen (O) biosynthesis loci (see e.g. CG258 in Figure 2b, and Supplementary Figure 2). Among the 17 clones with ≥1 detectable recombination hotspot (arbitrarily defined as mean recombination count of ≥5 per base calculated over non-overlapping 1000 bp windows), the galF K locus gene was ranked among the top 2% recombination counts in 16 clones (Figure 2c and Supplementary Figure 2). Consistent with these findings, 20 clones were associated with ≥3 distinct K loci and 11 clones were also associated with ≥3 O loci (Figure 2d and Supplementary Figures 3, 4). Together these data indicate that the capsule and LPS are subject to strong diversifying selection. However, this was not the case for the hypervirulent clones, which were associated with low K and O locus diversity: five out of the six had just one K and one O locus type (either KL1 or KL2, plus O1/O2v1 or O1/O2v2) and showed no evidence of recombination events affecting galF (Figure 2c,d Supplementary Figures 3, 4).

The key selective drivers for K/O locus diversity are not known, but the mammalian immune system is unlikely to play a major role since Kp live ubiquitously in the environment and are opportunistic rather than obligate human pathogens³⁷. Instead, phage and/or protist predation are likely candidates³⁸. Numerous capsule specific Kp phage have been reported^{39, 40} and ecological modelling supports a key role for phage-induced selective pressures in maintaining surface polysaccharide diversity in free-living bacteria⁴¹. The relative lack of diversity among the hypervirulent clones may suggest that they are not subject to the same selective pressures, perhaps indicating some sort of ecological segregation. This possibility is intriguing and could explain the separation of hypervirulence and MDR, by limiting opportunities for horizontal gene transfer between MDR and hypervirulent clones. Isolates representing both clone types have been identified among diverse host-associated niches^31–42,43 but it is not possible to determine any particular ecological preference due to the lack of systematic sampling efforts to-date. An alternative explanation is that the hypervirulent clones are subject to some sort of mechanistic limitation for chromosomal recombination, that in turn limits surface polysaccharide diversity and the acquisition of other chromosomally encoded accessory genes, as have recently shown to be frequently acquired by CG258 strains⁴⁴. If so, we may also expect a general trend towards lower gene content diversity in the hypervirulent clones.

To assess overall gene content diversity we conducted a pan-genome analysis using Roary⁴⁵. Jaccard gene content distances were generally lower for genome pairs within hypervirulent clones than the MDR or unassigned clones, suggesting the former have less diverse pan-genomes (p < 1×10⁻¹⁵ for each pairwise Wilcoxon Rank Sum test, Figure 3a). Supporting this trend, the hypervirulent clones were associated with comparatively shallow pan-genome accumulation curves (Figure 3b). In order to quantify the differences in these curves we fitted the pan-genome model proposed by Tettelin and colleagues⁴⁶, and derived an alpha value for each clone (Supplementary Table 1), whereby values <1 indicate an open pan-genome and >1 indicate a closed pan-genome. Consistent with previous data showing extensive gene content diversity within the Kp species⁵, all but two clones had alpha values below 1. The exceptions were hypervirulent CG380 (alpha = 1.06) and MDR CG231 (alpha = 1.16). There was a general trend towards higher alpha (i.e. less open) among the hypervirulent clones (median alpha = 0.77, IQR 0.70–0.89), although the difference was only statistically significant in comparison to the unassigned clones (median alpha = 0.52, IQR 0.51–0.54, p = 0.005) and not the MDR clones (median alpha = 0.62, IQR 0.59–0.69, p = 0.23, Figure 3b).

It is well known that large groups of accessory genes can be linked on the same mobile element (e.g. large conjugative MDR plasmids that are common in the MDR clones, or the virulence plasmids characteristic of the hypervirulent clones), so a single gain or loss event may have a large effect on gene-based measures such as pairwise Jaccard distances and accumulation curves. Hence we used a principal component analysis (PCA) to generate a metric that is less sensitive to the correlation structure in the gene content data (see Methods). The PCA transformed the accessory gene content matrix comprising 1092 genomes vs 39375 genes into coordinates in a 40-dimensional space. These 40 axes captured >60% of the variation in accessory gene content and were used to calculate the Euclidean distance of each genome to its clone centroid. The resulting distributions of distances provided further support that the MDR and unassigned clones display greater gene content variation than hypervirulent clones (Figure 3c; p < 1×10⁻¹⁵, 2 d.f., Kruskal-Wallis test; p < 1×10⁻¹⁵ for each pairwise Wilcoxon Rank Sum test) and suggest this is associated with a greater frequency of horizontal gene transfer events rather than a similar number of events introducing larger changes in gene content. In addition, the putative ancestry of accessory genes (see Methods) was more diverse among MDR and unassigned clones than the hypervirulent clones, supporting that the latter are subject to a more limited range of partners for horizontal gene transfer (Wilcoxon Rank Sum tests: hypervirulent vs MDR, p = 0.0027; hypervirulent vs unassigned, p = 1×10⁻⁴; MDR vs unassigned, p = 0.38; Figure 3d, Supplementary Table 1).

To further explore differences in common sources of accessory gene diversity, we assessed phage and plasmid diversity. For each genome we summed the length of genomic regions identified as phage by VirSorter⁴⁷ (range 0–221 kbp, Supplementary Table 2) and used a PCA of phage-associated gene content to calculate distance to clone centroids as for the total pan-genome (Figure 4a,b and Supplementary Figure 6). Hypervirulent clones showed similar phage load and diversity to the unassigned clones (Figure 4a,b; Wilcoxon Rank Sum tests, hypervirulent vs unassigned: load, p = 0.15; diversity, p = 0.15), whereas MDR clones were generally associated with higher load and diversity than both the hypervirulent and unassigned groups (p < 1×10⁻¹⁵ for each pairwise comparison). Although these analyses were dependent on the quality and breadth of the underlying viral sequence database which may be subject to species bias, we have no reason to expect that this would be skewed with respect to MDR over hypervirulent clones. Hence it is clear that Kp, and in particular the MDR clones, are subject to frequent attack by diverse phage.

Unfortunately it is not possible to reliably identify plasmid sequences from draft genome assemblies⁴⁸. Instead we used plasmid replicon and relaxase (mob) typing as indicators of plasmid load and diversity. Each genome contained 0–12 of 69 uniquely distributed replicon markers and 0–23 mob-positive assembly contigs (detected by screening against the PlasmidFinder⁴⁹ database and mob PSI-BLAST^{50, 51}, respectively; Supplementary Table 2). MDR and unassigned genomes harboured a greater number of replicon markers than hypervirulent genomes, largely driven by low replicon loads in CG23 that was overrepresented among the hypervirulent genomes (Figure 4c,d, Supplementary Figure 7; Wilcoxon Rank Sum tests: hypervirulent vs unassigned, p = 1.5×10⁻⁴; MDR vs unassigned, p < 1×10⁻¹⁵; MDR vs hypervirulent, p < 1×10⁻¹⁵). There were no significant differences between the hypervirulent and unassigned groups for counts of mob-positive contigs per genome, and comparatively small differences between the MDR and unassigned or hypervirulent groups (Wilcoxon Rank Sum test: MDR vs unassigned, p = 1×10⁻⁶; MDR vs hypervirulent, p <1×10⁻⁶, Supplementary Figure 7).

Comparison of effective Shannon’s diversity of replicon profiles indicated that the hypervirulent clones harbour less plasmid diversity than either of the unassigned or MDR clones (not driven solely by CG23, see Figure 4d, Supplementary Figure 7, Wilcoxon Rank Sum tests: hypervirulent vs MDR, p = 0.0013; hypervirulent vs unassigned, p = 0.0015; MDR vs unassigned, p = 0.44). Similar trends were seen for effective Shannon’s diversity of mob types but the differences were not statistically significant after Bonferroni correction for multiple testing (n=3 tests), a finding that is not surprising given that far fewer mob types have been defined and that only ∼48% completely sequenced Enterobacteriaceae plasmids deposited in GenBank could be mob typed (whereas ∼83% could be replicon typed)⁵⁰.

While these data are subject to the biases of the underlying databases within which clinically relevant (MDR and virulence) plasmids are overrepresented, they are also consistent with the findings above regarding overall pan-genome diversity. These data imply that MDR clones frequently acquire and lose plasmids, consistent with the high plasmid diversity reported previously for ST258^{52, 53} and several others⁷, and with data from recent investigations of Kp circulating in hospitals which showed that individual plasmids transferred frequently between clones^{54, 55}. In contrast, the hypervirulent clones were associated with comparatively low plasmid diversity, mirrored by generally narrow plasmid load distributions. Taken together these data imply that hypervirulent clones acquire novel plasmids infrequently but can stably maintain them. For example, we recently estimated that the virulence plasmid, which by definition is highly prevalent in these clones, has been maintained for >100 years in CG23³¹. In addition, laboratory passage experiments have shown that hypervirulent strains can maintain MDR plasmids introduced in vitro²⁹, and we showed that a horse-associated subclade of CG23 has maintained a single MDR plasmid for at least 20 years³¹.

The combination of infrequent plasmid acquisitions and limited chromosomal recombination suggests that hypervirulent clones may be subject to particular constraints on DNA uptake and/or integration. One possible explanation is that these Kp clones possess enhanced defences against incoming DNA such as CRISPR/Cas or restriction-modification (R-M) systems. However, our genome data reveals no significant differences in either system (see Supplementary Text, Supplementary Figures 8–11). Alternatively, the key virulence determinants themselves, or other proteins encoded on the virulence plasmid, may play a role. Two variants of the virulence plasmid predominate among hypervirulent clones and share limited homology aside from the iuc, iro and rmpA loci²⁷. It seems unlikely that a siderophore system would influence DNA uptake, however it is conceivable that upregulation of capsule expression by rmpA^{18, 56} may play a role by exacerbating the inhibitory effect of the capsule.

Capsule expression has been associated with a comparative reduction in Kp transformation frequency in vitro⁵⁷ and in a natural Streptococcus pneumoniae population⁵⁸. Additionally, the capsule is known to conceal the LPS⁵⁹ which, together with the OmpA porins, are considered key target sites for attachment of conjugative pili during the initial phases of mate-pair formation^{60, 61}. Hence we speculate that overexpression of the capsule in hypervirulent clones may result in a reduction of DNA uptake. Given that capsule types differ substantially in their thickness and polysaccharide composition^{18, 62}, it is also likely that their influence on DNA uptake is type dependent. The K2 capsule, which is associated with five of the 6 hypervirulent clones investigated here is considered among the thicker capsule types¹⁸ and thus may have a comparatively greater influence. We used our genome data to test this hypothesis by comparing the genomic diversity of KL2 and non-KL2 genomes within MDR CG15, the only clone with sufficient KL2 and non-KL2 genomes for comparison. CG15-KL2 genomes formed a deep branching monophyletic subclade consistent with long-term maintenance of KL2 for an estimated 34 years (see Supplementary Text and Supplementary Figures 12 and 13). This KL2 subclade showed a comparatively low rate of recombination (r/m 0.58 vs 6.75) and more limited gene content diversity than the rest of the clone (p = 0.0004, Supplementary Text and Supplementary Figure 12). However, gene-content diversity in the KL2 subclade of CG15 was higher than that of the hypervirulent clones (Supplementary Figure 12), perhaps due to the absence of the rmpA capsule upregulator. Thus the genome data support our hypothesis and should motivate future laboratory studies of this phenomenon; a task that will not be trivial given the low efficiency of in vitro transformation for wild-type Kp strains⁶³, the challenge of identifying suitable selective markers for distinguishing MDR strains, and the sensitivity of conjugation efficiencies to laboratory growth conditions⁶⁰. If confirmed, this would imply that hypervirulent clones are evolutionarily constrained by a key determinant of the hypervirulent phenotype, and as such are self-limited in their ability to adapt to antimicrobial pressure.

Regardless of the mechanisms, our data clearly show that hypervirulent Kp clones are less diverse than their MDR counterparts, and suggest that the rate of virulence plasmid acquisition by MDR clones will far exceed the rate of MDR plasmid acquisition by hypervirulent clones. This is particularly worrying from a hospital infection control perspective since many of the MDR clones investigated here appear well adapted to transmission and colonisation in the human population, and are frequent causes of hospital outbreaks^{7, 8}. Given the mounting evidence that MDR clones can carry multiple plasmids at limited fitness cost^64–66 and frequently exchange plasmids with other bacteria^{54, 55}, it seems these MDR clones may also be the perfect hosts for consolidation and onwards dissemination of MDR and virulence determinants. The greatest concern is that these determinants will be consolidated onto a single mobile genetic element; indeed mosaic Kp plasmids carrying AMR genes plus iuc and rmpA2 have already been reported in an MDR Kp clone³⁰, and Escherichia coli plasmids bearing iuc, rmpA and AMR genes have been detected in Kp²⁷. Whether these convergent strains and plasmids are fit and disseminating is not known. Recent experience with convergent carbapenem-resistant CG258 in China – which retrospective surveillance studies showed was already widely disseminated at the time of the outbreak report^{6, 67} – highlights the ease with which deadly strains can circulate unnoticed. As reports of convergent Kp strains continue to increase, the need for global genomic surveillance encompassing clone, AMR and virulence locus information⁶⁸ is clearly greater than ever.