
==== Front
Cell Rep
Cell Rep
Cell Reports
2211-1247
Cell Press

S2211-1247(24)00941-0
10.1016/j.celrep.2024.114602
114602
Article
Population genomics uncovers global distribution, antimicrobial resistance, and virulence genes of the opportunistic pathogen Klebsiella aerogenes
Feng Yu 12
Yang Yongqiang 12
Hu Ya 12
Xiao Yuling 3
Xie Yi 3
Wei Li 4
Wen Hongxia 12
Zhang Linwan 5
McNally Alan 6
Zong Zhiyong zongzhiy@scu.edu.cn
127∗
1 Center of Infectious Diseases, West China Hospital, Sichuan University, Chengdu, China
2 Center for Pathogen Research, West China Hospital, Sichuan University, Chengdu, China
3 Laboratory of Clinical Microbiology, Department of Laboratory Medicine, West China Hospital, Sichuan University, Chengdu, China
4 Department of Infection Control, West China Hospital, Sichuan University, Chengdu, China
5 Department of Clinical Research Management, West China Hospital, Sichuan University, Chengdu, China
6 Institute of Microbiology and Infection, College of Medical and Dental Science, University of Birmingham, Birmingham, UK
∗ Corresponding author zongzhiy@scu.edu.cn
7 Lead contact

12 8 2024
27 8 2024
12 8 2024
43 8 11460213 2 2024
13 6 2024
23 7 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an open access article under the CC BY-NC license (http://creativecommons.org/licenses/by-nc/4.0/).
Summary

Klebsiella aerogenes is an understudied and clinically important pathogen. We therefore investigate its population structure by genome analysis aligned with metadata. We sequence 130 non-duplicated K. aerogenes clinical isolates and identify two inter-patient transmission events. We then retrieve all publicly available K. aerogenes genomes (n = 1,026, accessed by January 1, 2023) and analyze them with our 130 genomes. We develop a core-genome multi-locus sequence-typing scheme. We find that K. aerogenes is a species complex comprising four phylogroups undergoing evolutionary divergence, likely forming three species. We delineate remarkable clonal diversity and identify three worldwide-distributed carbapenemase-encoding clonal clusters, representing high-risk lineages. We uncover that K. aerogenes has an open genome equipped by a large arsenal of antimicrobial resistance genes. We identify two genetic regions specific for K. aerogenes, encoding a type VI secretion system and flagella/chemotaxis for motility, respectively, both contributing to the virulence. These results provide much-needed insights into the population structure and pan-genomes of K. aerogenes.

Graphical abstract

Highlights

• Klebsiella aerogenes is a species complex comprising multiple taxa

• A core-genome multi-locus sequence-typing scheme is established

• Three carbapenem-resistant high-risk lineages are distributed worldwide

• A T6SS and a motility-encoding region are specific and contribute to the virulence

Feng et al. provide the genome sequence of 130 Klebsiella aerogenes clinical isolates with analysis of 1,026 additional publicly available genomes. The analysis uncovers K. aerogenes as a species complex with remarkable clonal diversity, identifies three worldwide-distributed carbapenemase-encoding clonal groups, and detects two specific genetic regions contributing to the virulence.

Keywords

Klebsiella
Klebsiella aerogenes
genome analysis
antimicrobial resistance
cps
cgMLST
Published: August 12, 2024
==== Body
pmcIntroduction

Klebsiella aerogenes, sometimes (incorrectly) called Klebsiella mobilis in the literature,1 is a bacterium of clinical significance.2 K. aerogenes has been previously known as Enterobacter aerogenes, but genome-based analysis has shown that K. aerogenes is actually more closely related to Klebsiella pneumoniae than Enterobacter species.3 K. aerogenes is able to cause a variety of infections such as bloodstream infections, pneumonia, and urinary tract infections, the majority of which are healthcare associated.4 Outbreaks due to K. aerogenes have been reported, particularly in intensive care units (ICUs).5,6 Infections due to K. aerogenes could lead to poor outcome for affected patients.7 For instance, previous studies have uncovered a 13.5%–28% mortality rate in patients with K. aerogenes infections.7,8,9,10

As a Klebsiella species, K. aerogenes is a member of the ESKAPE multi-drug-resistant pathogens.11 Compared with other Klebsiella species, K. aerogenes has two hallmark features. First, K. aerogenes has an intrinsic, chromosomally located ampC gene, which encodes AmpC, a class C β-lactamase. The carriage of intrinsic ampC is different from that in all other Klebsiella species, which have class A β-lactamase (e.g., SHV and OXY) genes. Class C β-lactamases are found in the mnemonic SPICE group of pathogens (Serratia, Pseudomonas, indole-positive Proteus, Citrobacter, and Enterobacter).12 Hyperproduction of AmpC β-lactamases can rapidly confer resistance to most β-lactams during therapy with third-generation cephalosporins.13 Notably, acquired genes encoding various plasmid-borne AmpC β-lactamases,14 extended-spectrum β-lactamases,15,16,17 and carbapenemases4,18,19,20,21,22 have also been well reported in K. aerogenes. Second, the motility of K. aerogenes is a distinguishing feature within the typically non-motile Klebsiella genus.23 Despite the clinical significance, unlike other ESKAPE pathogens, K. aerogenes is understudied. A previous analysis of 97 publicly available K. aerogenes genomes attempted to determine the population structure and pan-genomes of this species.24 However, such a limited number of genomes is not adequate for analyzing pan-genomes and uncovering the population structure of K. aerogenes. The dissemination of strains across different geographic locations also needs to be examined to inform infection control.

To address the knowledge gap, we performed a genome-based study aligned with clinical data or metadata. We sequenced 130 non-duplicated K. aerogenes clinical isolates and retrieved all publicly available K. aerogenes genomes (n = 1,026), making up a dataset of 1,156 genomes for pan-genome analysis and comparative genomics. These analyses allow us to unfold the complex taxonomy, diverse clonal structure, various antimicrobial resistance determinants, and virulence-related species-specific genetic clusters of K. aerogenes. These findings highlight the clinical relevance and the remarkable genomic plasticity of this opportunistic pathogen. We further established two schemes for high-resolution, convenient strain typing based on core-genome multi-locus sequence typing (cgMLST) and the nomenclature of its intrinsic AmpC-encoding gene, respectively. These two schemes will be invaluable for facilitating future studies of K. aerogenes.

Results

Genome sequencing of 130 K. aerogenes clinical isolates uncovered a wide spectrum of diseases and identified clonal transmission

We sequenced 130 non-duplicated clinical isolates of K. aerogenes from our hospital in Chengdu, China (Data S1). The 130 isolates were recovered from various types of clinical samples, among which blood (n = 64, 49.2%) was the most common, followed by sputum or respiratory aspirations (n = 14, 10.8%). Among the 64 cases of bloodstream infection due to K. aerogenes, 37 (57.8%) were primary without known sources, and 27 (42.2%) were secondary to urinary tract infection (n = 9), pneumonia (n = 8), bile duct infection (n = 4), intra-abdominal infection (n = 4), and wound or abscess (n = 2). This illustrates that K. aerogenes is associated with a wide range of infections including bloodstream infection, meningitis, peritonitis, pneumonia, pleurisy, urinary tract infection, and wound infection.

The 130 isolates could be assigned to 82 sequence types (STs) including 38 newly identified in this study, exhibiting a great diversity of clonal background. In our collection, most STs (n = 67, 81.7%) comprised a single isolate only and therefore did not exhibit intra-hospital transmission. Among all STs, ST4 was the most common type (n = 15, 11.5%), followed by ST14 (n = 11), ST93 (n = 6), and ST208 (n = 6). We determined core-genome single-nucleotide polymorphisms (SNPs; of note, all SNPs hereafter refer to core-genome SNPs) between isolates of the same STs (Data S2). The 15 ST4 isolates had a pairwise SNP distance ranging from 86 to 43,192, suggesting a non-clonal background. Using SNP typing, we uncovered three pairs of isolates belonging to the same clone. Two ST249 isolates with no SNPs were recovered from two samples that were processed consecutively in the microbiological lab on the same day, as evidenced by the sample IDs (the last four digits, 1135 and 1136). The two samples were collected from different patients who had stayed in different wards without any overlaps. K. aerogenes grew in large numbers from the first sample, while it grew in small quantity from the second one. The above findings pointed to a laboratory contamination from the first sample to the second sample rather than inter-patient transmission of ST249 isolates. In contrast, two ST14 isolates differing by four SNPs were recovered from samples collected on different days from two patients who been had hospitalized in the same patient room of a gastrointestinal surgery ward for 11 days. Two ST93 isolates differing by two SNPs were recovered from samples collected on different days, and the corresponding two patients were cared for by the same ICU medical team for 4 days despite being in different patient rooms. These findings suggest nosocomial transmission of a common strain, ST14 or ST93, in both cases.

K. aerogenes is a species complex comprising three species

We retrieved all publicly accessible genomes labeled K. aerogenes from GenBank, NCBI (accessed on January 1, 2023). We identified a total of 1,864 unique BioSample numbers of K. aerogenes. After quality control, we discarded 858 genomes from further analysis due to the absence of genomic data (n = 393), failure to pass quality control (n = 82), derivatives of metagenomic sequencing (n = 55), or species misidentification (not belonging to K. aerogenes, n = 308) (Figure 1). Therefore, we included 1,026 genomes of K. aerogenes, making a total of 1,156 (130 from this study; Data S3) for analysis. Among the 1,156 isolates, 1,121 had a specified geographic location and were recovered from 47 countries across six continents, namely North America (n = 454), Asia (n = 331), Europe (n = 226), Oceania (n = 54), South America (n = 34), and Africa (n = 22).Figure 1 Data flow chart of this study

KAC, K. aerogenes complex; KOC, K. oxytoca complex; KPC, K. pneumoniae complex; QC, quality control. Numbers of genomes are shown. We sequenced 130 KAC genomes and included all available KAC genomes from NCBI after quality control (accessed by January 1, 2023; n = 1,026), making a total of 1,156 KAC genomes for analyzing the population structure and developing a cgMLST scheme. For comparative genomics, we included 746 KAC genomes and 5,030 of other Klebsiella species as representatives after dereplication.

As illustrated by the SNP-based phylogenomic tree (Figure 2), K. aerogenes comprises four phylogroups. The vast majority (n = 989, 85.6%) of K. aerogenes isolates including the type strain KCTC 2190T (GenBank: CP002824) belong to phylogroup 1. Notably, one isolate (SRA: ERR7448948) is an outlier of phylogroup 1 and had a 97.0% average nucleotide identity (ANI) and a 75.6% in silico DNA-DNA-hybridization (isDDH) value with the type strain KCTC 2190T (Table S1). Phylogroup 2 comprised 88 isolates, which had a 98.6%–100% ANI and an 89.9%–100% isDDH between each other but a 96.8%–97.2% ANI and a 73.6%–76.0% isDDH with the type strain KCTC 2190T. This indicates that isolates of phylogroup 2 indeed belong to the species K. aerogenes but may represent a subspecies distinct from phylogroup 1. In contrast, phylogroup 3 comprised 54 isolates with an intragroup ANI of 98.3%–100% and an intragroup isDDH of 87.4%–100% but a 95.7%–96.0% ANI and a 65.6%–66.7% isDDH compared to the type strain KCTC 2190T. The ANI values fall into the 95%–96% inconclusive zone for species demarcation,25,26 but the isDDH values are well below the 70% cutoff to define species.27 It is therefore likely that phylogroup 3 represents a novel species. Similarly, the 24 isolates of phylogroup 4 had an intragroup ANI of 98.5%–100% ANI and an intragroup isDDH of 89.7%–100% but a 96.0%–96.2% ANI and a 68.2%–68.9% isDDH with the type strain KCTC 2190T. The ANI values are slightly higher than the 96% cutoff for defining species, while the isDDH values are still below the 70% species-defined cutoff, suggesting that phylogroup 4 is also likely to represent a new species. Notably, the ANI values between isolates of phylogroup 4 and those of phylogroup 3 are 95.8%–96.3% and the isDDH values were 66.8%–69.0%, similar to those between phylogroup 4 and phylogroup 1. This suggests that phylogroup 4 may be an intermediate species between K. aerogenes (sensu stricto, represented by phylogroup 1) and the new species represented by phylogroup 3. The above findings suggest that K. aerogenes may comprise three species (K. aerogenes sensu stricto and two new species corresponding to phylogroups 3 and 4, respectively) and two subspecies (corresponding to phylogroups 1 and 2). As such, we refer to K. aerogenes as K. aerogenes complex (KAC) in the following text.Figure 2 The radial phylogenomic tree of KAC

All available genomes of KAC from NCBI (accessed by January 1, 2023) after quality control and the 130 genomes in this study were included. The tree was inferred from concatenated SNPs found in core genes shared by at least 99% genomes using IQ-TREE under GTR+GAMMA model with ascertainment bias corrected and 1,000 bootstraps. This tree uncovers four phylogroups (phylogroups 1–4) with a notable outlier (SRA: ERR7448948) close to phylogroup 1.

To further explore the taxonomic positions, we also performed phenotypic tests for 18 KAC isolates comprising four or five per phylogroup. Consistent with previous findings,28,29 all tested KAC isolates regardless of phylogroups were positive for β-galactosidase, lysine decarboxylase, ornithine decarboxylase, citrate utilization, and Voges-Proskauer reaction and negative for arginine dihydrolase, H2S production, urea hydrolysis, deaminase, indole production, and gelatinase (Table S2). In addition, all tested KAC isolates were able to utilize adonitol, amygdalin, arbutin, D-arabinose, D-arabitol, D-cellobiose, D-fructose, D-galactose, D-glucose, D-lactose (bovine origin), D-maltose, D-mannitol, D-mannose, D-melibiose, D-sorbitol, D-raffinose, D-ribose, D-trehalose, D-xylose, esculin ferric citrate, gentiobiose, glycerol, inositol, L-arabinose, L-fucose, L-rhamnose, salicin, and sucrose, but not amidon (starch), glycogen, L-arabitol, potassium gluconate, potassium 2-ketogluconate, potassium 5-ketogluconate, or xylitol (Table S2). Notably, isolates of phylogroup 4 can be clearly differentiated from isolates of the other three phylogroups by their ability to utilize D-tagatose or the inability of fermenting methyl-α-D-glucopyranoside (Table S2). In contrast, isolates of phylogroup 3 exhibited highly similar patterns with those of phylogroups 1 and 2 but can be differentiated from isolates of the other three phylogroups by their weak ability to utilize N-acetylglucosamine (Table S2). This supports the suggestion that phylogroups 3 and 4 represent two novel species within the complex.

KAC exhibits a remarkable clonal diversity comprising 366 STs and 214 clonal clusters

The 1,156 KAC isolates could be assigned to 366 STs including 176 newly identified herein, which comprise the aforementioned 38 new STs from our 130 isolates. Phylogroups 1, 2, 3, and 4 comprised 292, 36, 23, and 14 STs, respectively, while the outlier strain (SRA: ERR7448948) belonged to ST623, a new ST comprising only this isolate at present. ST93 (phylogroup 1) was the most common type comprising 213 isolates, followed by ST4 (phylogroup 1, comprising 89 isolates), both of which had a worldwide distribution (Figure 3). In contrast, ST14 (phylogroup 1), which was common in our collection, comprised two additional isolates in the NCBI, one from a non-specified location of China and the other with unknown geographic information. This suggests that ST14 is likely a lineage largely restricting to our local settings. Notably, 216 STs comprised a single isolate only and there were only 15 STs containing ten or more isolates.Figure 3 The circular phylogenomic tree of KAC with isolate information

All available genomes of KAC from NCBI (accessed by January 1, 2023) after quality control and the 130 genomes in this study were included. The year, continent, and source of recovery of each isolate are indicated. Sequence types comprising 10 or more isolates are shown. The tree was inferred from concatenated SNPs found in core genes shared by at least 99% genomes using IQ-TREE under GTR+GAMMA model with ascertainment bias corrected and 1,000 bootstraps.

To understand the clonal structure, we also assigned the 1,156 KAC isolates into clonal clusters (CCs) based on whole-genome-based phylogeny with recombination sites accounted for and cluster partitioning. A total of 214 CCs were assigned with designations based on the descending order of the number of isolates identified within each group (Figure S1). Notably, isolates of the same ST may be assigned to multiple CCs, while isolates of different STs could be clustered in the same clonal group (see Figure 4 for an overview, Figure S2 for the top three STs and CCs, and Data S4 for each ST and CC). For instance, ST93 isolates were assigned to two CCs (CC1 [n = 156] and CC3 [n = 57]), while CC1 also comprises isolates of ST475, ST513, ST515, ST567, ST571, ST575 ST585, and ST597, and CC3 also consists of those of ST388, ST453, ST511, ST596, and ST598 (Figure S2 and Data S4). Similarly, ST4 was partitioned into CC2 (n = 60), CC13 (n = 14), CC22 (n = 10), and CC35 (n = 5), while CC2 also comprises ST550 and ST594 (Figure S2 and Data S4).Figure 4 STs and CCs

This chord diagram illustrates the reassignment of traditional sequence types (left) to new CC (right). Only the top 15 STs and CCs are individually labeled, while all others are collectively categorized into “others.” For the two predominant STs, ST4 and ST93, their respective paths to the designated CCs are distinctively highlighted with dashed lines.

The overlaps of STs and CCs, the diffuse distribution of isolates of the same ST, and the wide range of SNPs among isolates of the same ST (e.g., 86 to 43,192 SNPs between ST4 isolates) indicate that the current MLST scheme cannot accurately reflect the clonal structure of KAC and need to be modified due to the discordance arising from recombination events (see below for details). CCs may be a better option to address the clonal structure and study strain transmission.

The discordance between STs and CCs is largely driven by recombination

To investigate the discordance between STs and CCs, we performed an ancestral recombination analysis using fastGEAR30 on two representative STs (ST4 and ST93) and their associated CCs. The analysis revealed that both ST4 and ST93 had numerous ancestral recombination blocks. Specifically, CC13 and CC35 of ST4 appear to have originated from CC2, with increased genetic distance likely due to the exchange of genetic regions between strains other than ST4, without affecting the seven MLST loci. Conversely, the most divergent lineage, CC22 of ST4, involves an originally non-ST4 strain that received several genetic regions from strains of CC2 and CC13, including three MLST loci (dnaA, gyrB, and leuS), resulting in convergence into the same ST4 but with a distinct genomic background (Figure S3). Instances such as ST458 and ST474, which share similar genetic backgrounds with phylogenetic neighbors, appear to be assigned different STs due to small-scale variations affecting MLST loci.

For strains associated with ST93, apart from several scattered STs, the main population of ST93 is divided into CC1 and CC3. Unlike the observations in ST4, significantly larger ancestral recombination blocks seem to be the predominant source responsible for converging two distantly related populations into the same ST (Figure S3). A possible explanation could be that these two CCs share the same locus leuS, with a large ancestral recombination block from CC1 harboring the other six loci, replacing the corresponding genomic region of CC3, resulting in ST93. However, since fastGEAR assumes the larger population as the donor when reporting, the actual direction of recombination remains uncertain. It is also possible that both CC1 and CC3 were initially the same population, but the latter received large genomic regions from an unidentified origin, which coincidentally shared the same locus leuS without affecting the other six loci. While this led to the same ST, it significantly elevated the genetic distance between the non-recombined and recombined populations CC1 and CC3, respectively.

Many KAC CCs have international or intercontinental distribution and comprise multiple country-specific clones

We categorized CCs as international or intercontinental based on their presence in at least three countries or at least three continents, respectively. As a result, 58 CCs were identified as international and 50 as intercontinental. Notably, CC1, CC2, and CC3 demonstrated worldwide distribution, whereby CC1 was found in 17 countries spanning five continents, CC2 in 11 countries across four continents, and CC3 in nine countries over five continents (Data S4). In addition, for the two most prominent sequence types, ST4 and ST93, all ST93 isolates and most (67.42%, 60/89) ST4 isolates belonged to CC1, CC2, or CC3 (Data S4). We therefore focused on these three CCs and performed further investigations on the clonal relatedness based on SNPs.

Within each of these CCs, we aligned isolates against their corresponding reference complete chromosomes, namely G7 (GenBank: CP011539) for CC1, Y1 (GenBank: CP045870) for CC2, and FDAARGOS_513 (GenBank: CP033817) for CC3, as shown in Figure S4. Using the phylogenomic assessment of CC1, CC2, and CC3, we uncovered a notable pattern in strain distribution, elucidated by a heatmap (Figure S4) that quantifies the SNP differences among isolate pairs. We classified isolates exhibiting a maximum of 22 SNPs as clonal, indicating a recent common ancestry or an active transmission. This clonal cutoff was established based on the following procedures and considerations. First, we determined the average nucleotide substitution rate for isolates of CC1, CC2, or CC3, respectively. However, the root-to-tip analysis revealed a weak correlation between genetic distance and time for all three clusters (Figure S5), posing challenges for convergence within a sensible length of Markov chain Monte Carlo (MCMC). The weak temporal signal could be explained by numerous recombination events (as shown in Figure S6) as well as sampling bias due to missing data and uneven sampling over time. Nevertheless, we obtained the substitution rate of the entire chromosome for CC2 isolates, which is 6.68 per genome per year, with a 95% credible interval ranging from 3.08 to 10.6. This translates to a maximum of 21.2 SNPs per year, which we round to 22 SNPs, representing the maximum genetic distance observed between any two CC2 strains that diverged from a common ancestor within 1 year. Second, as true transmission events may be masked by an inadequate sampling rate, we considered isolates that share a common ancestor within a year as clonal. Notably, the cutoff of 22 SNPs to define clonal is similar to the 21-SNP cutoff that is commonly used for K. pneumoniae.31 We identified a total of 27 clones (16 of CC1, six of CC2, and five of CC3) with an average SNP distance of 9.43, ranging from 0 to 26 SNPs (Data S5). These clones were predominantly associated with human hosts, with exceptions of a single environmental clone and two clones with unknown source. At the fine-scale level of clones, we did not detect isolates with international spread; instead, we observed domestic clones circulating within individual countries (Data S5). Notably, there are 15 clones of the three CCs circulating in the United States (US), which may be due to the sampling bias toward sequencing more American isolates but could suggest the US as a hotspot for the dispersal of KAC isolates.

A robust cgMLST scheme for KAC typing was developed

To enhance the resolution of strain typing, we attempted to establish a robust cgMLST scheme. We identified a total of 3,430 genetic loci shared by 99% of KAC genomes. We selected loci for cgMLST under the guidance of a comprehensive analysis utilizing three distinct indices, with the inter-intra clusters variation index proving to be the most crucial in determining loci inclusion. We applied an optimal threshold of 0.9 for the inter-intra clusters variation index in conjunction with a cutoff value of 432 (see STAR Methods for details) and then developed a cgMLST scheme consisting of 2,067 loci. The robustness of this scheme was evidenced by its statistical performance, yielding a Fowlkes-Mallows (FM) index of 97.85% and a balanced accuracy (BA) of 97.87%, with 95.75% sensitivity and 100% specificity. The FM index measures the similarity between the inferred clustering (based on cgMLST) and the true clustering (based on phylogenetic analysis), with a value close to 100% indicative of a high level of agreement between the two clustering methods.32 BA considers both sensitivity and specificity, giving equal importance to both, with a value close to 100% indicating that the scheme (our cgMLST) is highly accurate in correctly classifying strains.33 The sensitivity and specificity measure the proportion of true positives (correctly identified strain pairs) out of all actual positives (true strain pairs) and that of true negatives (correctly identified single strains) out of all actual negatives, respectively.34 These metrics supported the high degree of precision and reliability in strain categorization for KAC.

KAC has multiple capsular types that are strongly associated with STs, CCs, and phylogroups

The capsule is an important feature of Klebsiella including KAC, and a well-defined typing scheme would greatly facilitate linking virulence and phage therapy to specific KAC isolates once their associated phenotypes have been determined. We identified a total of 24 primary capsular polysaccharide synthesis (cps) STs, denoted KL1 to KL24, ranging from 21,315 bp to 36,017 bp, with the highest pairwise coverage of 65.81% as determined using the Kaptive algorithm (Figure S7 and Data S7). Additionally, we reconstructed 26 variants harboring truncated or inserted versions of the corresponding primary KLs, ranging from 16,388 bp to 37,966 bp, based on assembly graphs. Notably, KAC cps KL17 is nearly identical to K. pneumoniae cps KL62, as visualized by its greatest distance from other KLs (Figure S7).

Upon determining the reference sequence for each KL, we typed all 1,156 KAC strains in the dataset (Data S6). Among these, five strains harbored KAC KL17; four strains were untypable due to gaps or assembly dead ends in the cps region; and eight strains were only typable based on assembly graphs. The results showed that the top three primary KLs, KL1, KL2, and KL3, accounted for almost 60% of the cases, with 441, 142, and 99 strains, respectively. KL3 and KL8 exhibited the widest distribution across all phylogroups, with KL3 identified with 11 variant forms resulting from various additions and deletions of genetic regions by insertion sequences. The presence of KAC KL17, which shares the same cps loci with KL62 of K. pneumoniae, demonstrates the high mobility of the cps region, capable of switching between different lineages or even species. Surprisingly, no unique KL was identified for phylogroups 2, 3, and 4, suggesting that these phylogroups may have originated from phylogroup 1 and diverged from the K. aerogenes sensu stricto population due to accumulated genetic distances likely introduced by recombination. To enhance accessibility, all KL loci have been adapted to the Kaptive algorithm and made publicly available.

To determine whether cps locus typing aligns with predefined groupings, we performed association tests between primary KLs and STs, CCs, and phylogroup classifications. The results revealed a strong association between KL types and these grouping methods (p < 0.01), with Cramér’s V indicating the strongest correlation with MLST (Figure S7). This is expected, given that the number of categories generated by MLST is significantly higher than those by CCs or phylogroups. Although discordance between STs and CCs was observed, a single nucleotide difference in one of the seven MLST loci can create a new ST category, resulting in more distinct classifications. This higher granularity in MLST provides an advantage in Cramér’s V estimation, leading to stronger observed associations.

KAC carries a diversity of intrinsic ampC genes, likely formed by recombination

As already mentioned, KAC typically has an intrinsic ampC gene. Among the 1,156 genomes, 1,150 contained a complete ampC gene while the remaining six had an interrupted or incomplete ampC gene. Given the absence of a nomenclature system for AmpC of KAC, we proposed AER (from “aer” in K. aerogenes) as the name of these AmpC β-lactamases. The 1,150 complete blaAER ampC genes encode a total of 196 AER proteic variants, assigned AER-1 to AER-196 (Data S8) sharing ≥97.64% amino acid identity (Data S9). We inferred a phylogenetic tree based on blaAER gene sequences and another tree based on AER amino acid sequences (Figure 5). Both trees illustrated that blaAER alleles variants are largely clustered within phylogroups but with notable exceptions. A few blaAER genes of isolates belonging to a phylogroup were clustered together with isolates of another phylogroup (Figure 5). This illustrates inter-strain and inter-phylogroup transfers of blaAER.Figure 5 Phylogenetic trees of KAC ampC and AER

(A) Phylogenetic trees of KAC ampC genes.

(B) Phylogenetic trees of KAC AER enzymes.

The discovery of certain blaAER alleles clustering with those from different phylogroups, in the absence of mobilization evidenced by insertion sequences or mobile genetic elements, led to the hypothesis that homologous recombination may account for this unexpected phylogenetic pattern. To investigate this, we conducted an analysis on the blaAER gene and its adjacent genetic regions. We examined the chromosome of KAC and found that a 14,500-bp region between glaR (encoding a DNA-binding transcriptional repressor) and blaAER is conserved across KAC isolates. We used the recombination predictive model, implemented via Gubbins, and found a significant incidence of recombination events impacting both the blaAER and the adjacent helix-turn-helix-domain-containing gene (Figure S8). We then performed phylogenetic analyses by inferring phylogenetic trees both with and without accounting for recombination sites to assess the extent to which recombination had influenced phylogeny. By this side-by-side comparison of the phylogenetic trees, we revealed pronounced shifts in the placement of isolates among phylogroups due to recombination events. In particular, we noted that isolates within phylogroup 1 appear to acquire regions from other phylogroups, markedly affecting their positions in the trees (Figure S9). These recombination events, combined with other mechanisms of sequence variation, resulted in a total of 185 unique synonymous SNPs and 133 unique missense SNPs, along with three in-frame deletions and one in-frame insertion, across the entire blaAER gene. Taken together, these variations generated 196 AER variants within only 1,150 KAC strains, demonstrating its significant diversity (Data S10).

KAC possesses an arsenal of important acquired antimicrobial resistance genes and carries multiple virulence factors

Although there is a possible selection bias in favor of carbapenem-resistant strains for genome sequencing with deposit in databases, we detected carbapenemase-encoding genes in complete form in 217 (18.8%) genomes, including eight from our collection (Data S3; see Figure S10 for a heatmap of antimicrobial resistance genes alongside a KAC phylogenomic tree). KPC was the most common carbapenemase, seen in 105 genomes (KPC-2 in 79 and KPC-3 in 26), four of which also had NDM, a metallo-β-lactamase (MBL). NDM was the second most common carbapenemase, seen in 55 genomes (NDM-1 in 35 and NDM-4, -5, -6, -7, and -9 in the remainder). OXA-48-like carbapenemases were seen in 42 genomes (OXA-48 in 37, OXA-181 in four, and OXA-484 in one). IMP (IMP-1, -4, and -26), VIM (VIM-1, and -2), and NMC-A were seen in 13, five, and one genomes, respectively. We also aligned the distribution of carbapenemase-encoding genes with clones defined above. We therefore identified two major concerning clones, each comprising more than three isolates, namely clone 1 (comprising four isolates from Sweden) of CC1 encoding OXA-48 and clone 23 (comprising five isolates from Brazil) of CC3 encoding KPC-2 (Data S5). We also identified three other clones carrying carbapenemase-encoding genes, which comprise only two isolates for each, namely clone 14 (seen in the US) of CC1 encoding KPC-2, clone 15 (in Switzerland) of CC1 encoding NDM-5, and clone 21 (in the US) of CC2 encoding KPC-3 (Data S5).

Extended-spectrum β-lactamase (ESBL)-encoding genes were present in 119 genomes (10.3%) including eight carrying genes encoding two ESBLs (two CTX-M variants or a CTX-M variant plus SHV-12). Twelve CTX-M variants (CTX-M-1, -2, -3, -9, -12, -14, -15, -55, -59, -65, -68, and -199) were detected in 83 genomes with CTX-M-15 as the most frequent variant (seen in 41 genomes). Five SHV (SHV-5, -7, -12, -30, and -31) ESBL variants were detected in 37 genomes, with SHV-12 present in 31. Other ESBLs were GES-1 (n = 1), SFO-1 (n = 2), and TEM-24 (n = 2). Among the 217 genomes carrying carbapenemase-encoding genes, 67 (30.9%) encode ESBLs. Plasmid-borne colistin resistance mcr genes were detected in 24 genomes (2.1%) including mcr-1 in 12, mcr-9 in 11, and mcr-10 in one. Unlike carbapenemase-encoding genes, we did not detect clones comprising at least three isolates carrying ESBL or mcr genes, implying less importance of their carriage in dissemination.

A total of 52 virulence genes were detected (Table S3 and Data S3; see Figure S10 for a heatmap alongside KAC phylogenomic tree). These genes encode biofilm formation (mrkH), allantoin utilization (allS, arcC, fdrA, glxK, and ylbE), microcin (mceABGH), and iron acquisition (clbA-R, fyuA, iroBCDN, irp1, irp2, iucABCD-iutA, kfuABC, and ybtAEPQTUX). Among these genes, iutA was seen in all KAC genomes and kfuABC (1,152, 99.65%), mrkH (1,147, 99.22%), and iroBCDN (n = 1,065, 92.13%) were present in the vast majority of the 1,156 genomes (Table S3). Of note, rmpA (a regulator of the hypermucoviscous phenotype) was not found in KAC, while rmpA2 (another regulator) was only present in one isolate (SRA: SRR19409809), which belongs to ST432 of phylogroup 2 and was recovered from a human urine sample in the US. This isolate also contains the complete operon of aerobactin-encoding genes iucABCD and iutA, therefore meeting the criteria for defining hypervirulence for K. pneumoniae.35 Alarmingly, similarity searches predict that these virulence determinants reside on a plasmid of over 380 kb in size, which additionally carries an ESBL gene (blaCTX-M-15) and a carbapenemase gene (blaNDM-5) (Data S11) and has been described as both hypervirulent and multi-drug-resistant.36,37 Notably, this large plasmid has only been reported in K. pneumoniae strains, making this the first piece of evidence to demonstrate a KAC strain acquiring such a hypervirulent and multi-drug-resistant plasmid. However, the exact structure of this plasmid requires complete plasmid sequencing for confirmation.

KAC has an open genome revealed by pan-genome analysis

We identified a significant correlation between the core-genome branch length and the frequency of gene-exchange events including both the acquisition and loss of genes (predominantly, sharing accessory genes) from analyzing the pan-genome of KAC (p = 1.45e−89), with a core branch coefficient of 1.000033 (Figure S11). This coefficient is significantly different from zero (p = 2.2e−16), suggesting that the KAC pan-genome is actively evolving and should be regarded as an open genome. We also performed another, traditional, analysis using a saturation curve of unique gene clusters extending infinitely under Heaps’ law model and identified an α value of 0.338. As an α value of <1 signifies that new genes are continuously discovered with the sequencing of additional genomes, reflecting a high level of genetic diversity within the species,38 the α value of 0.338 indicates an open pan-genome (Figure S12). Notably, we found a significantly higher frequency of gene gain/loss events per unit of branch length for terminal branches compared to that of the internal ones (p = 4.18e−33) in the phylogenomic tree. This implies that gene gain and loss in KAC are primarily driven by mobile genetic elements that emerged relatively recently and have not yet widely disseminated in the population. This conclusion is also supported by the observation that finer or smaller-scale groupings, such as CCs and country-specific groups, showed significantly stronger associations with the functional traits linked to virulence, antimicrobial resistance, stress tolerance, and plasmid replicon, compared to broader groupings such as phylogroups. Specifically, we found 3,010 significant associations (p < 0.05) with CCs and 1,071 significant associations with country-specific groups, whereas only 205 significant associations were found when using phylogroups as the grouping method. These findings highlight the crucial role of mobile genetic elements in shaping subpopulations of KAC, as these functional determinants are common passengers of mobile genetic elements (Figure S13 and Data S12).

Further investigation into horizontal gene transfer among the four KAC phylogroups has generated notable insights. Phylogroup 1 harbors the largest accessory genome, but this is likely due to sampling biases as over 85% of the sampled isolates belong to this group, leading to a more diverse genomic representation. Importantly, the rate of horizontal gene transfer relative to the core-genome phylogeny in phylogroup 1 is not significantly higher than that of other phylogroups, supporting the idea that its larger accessory genome is a result of sampling bias rather than increased horizontal gene transfer. In contrast, by comparing gene gain/loss events between different pan-genome datasets, phylogroup 2 exhibits the highest rate of horizontal gene transfer, indicated by the highest core estimate, which is significantly greater than those observed in phylogroups 1 (p = 5.09e−9), 3 (p = 6.16e−3), and 4 (p = 0.0105). This indicates that phylogroup 2 is undergoing the most rapid evolution among all. Furthermore, there is no significant difference between the core estimates of horizontal gene transfer rates among phylogroups 1, 3, and 4. Additionally, there is no significant difference in the terminal branch rates and dispersion parameters across all four phylogroups. Given that the number of genes involved in each horizontal gene transfer event is indicated by the dispersion parameter, the lack of significant differences in this parameter suggests that the size of each gene-exchange event is consistent across all phylogroups. Terminal estimates, which represent the frequency of recent gene-exchange events, also showed no significant differences, indicating that the recent rate of gene-exchange events is stable across all phylogroups. Consequently, since both the size and frequency of gene-exchange events are consistent across all four phylogroups and considering that mobile genetic elements are the primary drivers of subpopulation diversification, we conclude that the evolutionary pattern of KAC, as evidenced by the emergence of subpopulations, remains uniform among all four phylogroups.

KAC contains two specific regions encoding flagella/chemotaxis and a type i3 T6SS

We performed comparative genomics using a set of tools to identify genes and pathways specific to the complex to enhance the appreciation of KAC. In the comparative genomics study, we included 746 KAC genomes and 5,030 of other Klebsiella species as representatives after dereplication. We identified 87 genes encoding products that met the stringent criteria of a minimum amino acid difference of 50% (i.e., <50% amino acid identity) between KAC genomes and other Klebsiella ones (Data S13). Within this subset, 60 genes were assigned with a Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology number, and 48 were found to be associated with at least one KAC pathway (Data S13). Notably, two pathways, flagellar assembly (involving 36 genes) and bacterial chemotaxis (involving eight genes), exhibit significant enrichment within the KAC-specific gene set (q value <0.01). These pathways share five genes and are adjacent on the chromosome, spanning a 58,771-bp region encoding approximately 66 open reading frames. This region corresponds to nucleotide positions 3,352,788 to 3,411,558 on the chromosome of K. aerogenes KCTC 2190T (GenBank: CP002824). We refer to this region as KAC-specific flagella region hereafter. Additionally, we identified a distinct 26,632-bp region, mapped between nucleotide positions 4,141,861 and 4,168,492 on the same chromosome. Although this region did not show significant enrichment in the over-representation analysis, 12 of the 24 genes in this region are members of the 87 KAC-specific genes (Data S13). The 12 genes encode for T6SS regulation or structural components rather than effectors. This region was characterized as a type i3 T6SS gene cluster39 and, notably, T6SS of KAC has not been included in the KEGG pathway database. Interestingly, type i3 T6SS was not present in non-KAC Klebsiella species, which contain type i2 T6SS instead (Data S14), suggesting horizontal gene transfer of T6SS genes into the KAC. The KAC-specific flagella and T6SS regions are highly prevalent in KAC, with only one strain (GenBank: GCA_011604725.1) lacking the flagella loci and another strain (GenBank: GCA_003075515.1) missing the T6SS region, in addition to two other strains (GenBank: GCA_001571545.2 and SRA: SRR6656655) having a truncated T6SS region. These two regions are distributed across the complex without significant difference of nucleotide identity in the main components of the flagella and T6SS systems except for several hypothetical proteins (Figure S14).

We systematically screened the presence of similar T6SS and flagella regions across 325,005 genomic assemblies from the NCBI RefSeq database, with results at per species or per genus level, as provided by the Genome Taxonomy Database (GTDB), shown in Data S15 (for the T6SS) and Data S16 (for the flagella region). We also performed homology analysis in an attempt to identify the potential origin of KAC-specific T6SS and flagella region.

KAC-specific T6SS has a common origin with that of Pluralibacter

We found a 2.8% overall detection rate for the KAC-specific T6SS across 325,005 genomic assemblies. Apart from K. aerogenes, 8,644 assemblies, spanning 159 species across 29 genera, exhibited at least 50% coverage for the KAC-specific T6SS region, with identities of the entire region ranging from 66.91% (Rouxiella) to 93.20% (Escherichia) (Data S14). However, only 0.03% (12 out of 38,803) of Escherichia contains such a T6SS region, suggesting acquisition from a foreign origin in this genus. Similarly, the second highest identity (80.65%) was seen in non-KAC Klebsiella, but only one of 21,523 genomes contain such a T6SS, indicative of a foreign origin. The next most similar matches (an average identity of 77.36%–80.00% and an average coverage of 80.70%–84.41%) were present in five genera, in all of which >50% isolates contained such a T6SS region, namely Pluralibacter (26 out of 26, 100%), Franconibacter (9 out of 12, 75%), Pseudescherichia (21 out of 32, 65.63%), Enterobacter (5,071 out of 5,270, 96.22%), and Cronobacter (714 out of 715, 99.86%) (Data S14). We identified 84 isolates from 60 GTDB-classified species containing the similar T6SS operon in homology analysis. Notably, in the analysis of T6SS operon, representatives from Escherichia, non-KAC Klebsiella, Scandinavium, and Pluralibacter closely clustered with the reference K. aerogenes, sharing a most recent common ancestor (MRCA) distinct from other isolates (Figure 6). Like those of Escherichia and non-KAC Klebsiella as mentioned above, such a T6SS region of Scandinavium has a foreign origin, as it was only present in only one of the 15 genomes of the genus. Conversely, all Pluralibacter genomes contain such a T6SS, suggesting its intrinsic nature in this genus. It is therefore likely that the T6SS regions of Pluralibacter and KAC have a common, yet-to-be-identified, origin.Figure 6 Homology assessment of KAC-specific T6SS region

This phylogenetic tree is based on concatenated amino acid sequences from 24 genes within the T6SS region. It includes 84 representative isolates, each with an average coverage of at least 80% relative to the T6SS region in the reference strain, K. aerogenes KCTC 2190T (highlighted in red on the tree). Bootstrap support values are illustrated by the colored branches on a gradient from red (70%) to green (100%). A heatmap that displays the similarity of each locus, with a color scale bar on the right indicating the similarity index, is plotted alongside the tree. The order of loci in the heatmap reflects the actual genetic arrangement of the T6SS operon. Representatives from Escherichia, non-KAC Klebsiella, Scandinavium, and Pluralibacter closely clustered with the reference K. aerogenes, sharing a most recent common ancestor distinct from other isolates.

KAC-specific flagella region originated from Serratia

We detected the flagella region similar (≥50% coverage and 60% identity) to the KAC-specific one in 0.69% of the 325,005 genomic assemblies encompassing 1,791 assemblies across 40 species of ten genera, with average identity of the entire region ranging from 73.61% to 84.28% (Serratia) (Data S16). The majority of these isolates were found in Serratia (n = 1,422) and its sister Serratia genera (Serratia_A, Serratia_B, Serratia_F, Serratia_I, and Serratia_J) defined by the GTDB (n = 316), accounting for 97.04% of the positive cases. Notably, in the GTDB there is a genus named Serratia in conjunction with multiple sister genera named Serratia_A to Serratia_J. Also of note, genomes of Serratia_C, Serratia_D, Serratia_E, and Serratia_H were not included in the NCBI RefSeq database due to failure to pass quality control. Within Serratia and its sister genera, the prevalence of KAC-specific flagella was over 90%, except for the putative genera Serratia_I (present in only 1/3) and Serratia_G (absent from the only one available genome). We identified 46 isolates from 32 GTDB-classified species containing the flagella operon in homology analysis. In the phylogenetic tree, KAC-specific flagella operon is clustered within those of Serratia (Figure 7), suggesting that the KAC-specific one is likely to originate from Serratia. We observed notable genetic diversities between flgL (encoding a flagellar hook-filament junction protein) and fliR (encoding one of the six flagellar export apparatuses) at species level, even within KAC (Figure S15). Additionally, we identified another hotspot of genetic diversity at strain level between fliE (encoding a structural protein of the flagellum hook-basal body complex) and fliT (encoding a putative filament) (Figure S15).Figure 7 Homology assessment of KAC-specific flagella region

This phylogenetic tree represents an analysis based on the concatenated amino acid sequences of 67 genes in the flagella region. A total of 46 representative isolates, each compared to the flagella region of the reference strain K. aerogenes KCTC 2190T, marked in red, with a minimum average coverage of 60%, are included. The tree’s bootstrap supports are denoted by a color gradient, ranging from red (70%) to green (100%). A heatmap is displayed alongside the phylogeny, illustrating the similarity of each locus, ordered by their actual genetic location in the flagella operon. KAC-specific flagella operon is clustered within those of Serratia, suggesting that the KAC-specific one is likely to originate from Serratia.

Discussion

In this study, we determined the population structure of K. aerogenes and demonstrate that K. aerogenes is a species complex comprising three taxa, exhibiting a great diversity of clonal background, causing a wide spectrum of infectious diseases, disseminating between patients and between humans and animals, carrying various determinants mediating resistance to clinically important antimicrobial agents such as carbapenems and colistin, and harboring two specific genetic regions related to virulence compared with other Klebsiella species. We believe that our study sheds a new light on this understudied clinically important pathogen, which may help to raise awareness toward this clinically significant pathogen and could provide useful information for infection control and clinical management.

We uncovered that K. aerogenes species is undergoing significant evolutionary divergence and forms four well-separated phylogroups. Phylogroups 1 and 2 encompass most genome sequenced isolates available publicly and belong to a single species, K. aerogenes (sensu stricto). In contrast, phylogroup 3 can be defined as an independent species according to criteria based on genome comparison. Further, the remaining phylogroup (phylogroup 4) is slightly closer to phylogroup 1 than phylogroup 3 and should also be regarded as a separate species. The above findings illustrate that K. aerogenes is actually a complex of three taxa.

We detected that the current MLST may not be reliable to address KAC population structure and uncovered that the discordance between STs and CCs is largely driven by recombination. These findings highlight the complexity and dynamics of genetic recombination in KAC, emphasizing the need for comprehensive genomic analyses to understand the evolutionary mechanisms driving the discordance between STs and CCs. Such observation in our study provides evidence that the discriminatory power of the current MLST scheme is significantly compromised by recombination events in the KAC population, particularly due to the uneven distribution of the seven loci across the genome where the proximity of loci such as dnaA and gyrB, and fusA and rplB, facilitates their co-transfer into the target genome, contributing to the observed discordance between STs and CCs.

As for KL types, due to the lack of phenotype data, the actual observable differences between different types remain unclear, preventing the establishment of a link between virulence and cps at this stage. This warrants further investigation in the future. Additionally, KL typing from genomic data can be greatly affected by poor sequencing and assembly quality. The combination of insertion sequences is a major force shaping the cps region, and those without a single intact cps region, although assigned a KL in this study, may actually represent a novel conformation of cps region, indicating a new primary locus or at least a new variant form of the primary KL. Consequently, the number of variant KLs should be interpreted with caution, as fewer variants could result from technical issues in whole-genome sequencing. Nevertheless, all reference sequences of cps loci in the KAC genomes are compatible with Kaptive typing and are available at https://github.com/fengyuchengdu/KAC_cps_reference.

We identified that KAC has two specific genetic regions encoding T6SS and flagella, respectively. The distinct divergence of over 50% in amino acid sequences for these genes between KAC and other Klebsiella species points to unique origins for the two genetic regions. Given that species of genus Klebsiella do cluster together at genome level, the introduction of KAC-specific T6SS and flagella may have occurred after speciation of the genus but prior to the diversification of KAC. A notable finding is the prevalence of type i3 T6SS in KAC species, a widely spread T6SS type in Citrobacter, Enterobacter, Kosakonia, Pluralibacter, and Phytobacter (Data S14), contrasting with the type i2 T6SS commonly found in other Klebsiella species. This characteristic is a hallmark of KAC species within Klebsiella genus and suggests horizontal gene transfer of T6SS genes. This consideration of horizontal gene transfer is further reinforced by the presence of type i3 T6SS in a tiny subset of Escherichia isolates (12 out of 38,803, 0.03%), which displays an exceptionally high identity of 93.20% with the KAC-specific T6SS. As there were only 0.03% with such a T6SS and the 12 isolates were recovered relatively recently (1990–2021), it is unlikely that Escherichia is an ancient donor of the T6SS. As mentioned in results, all Pluralibacter genomes contain an intrinsic type i3 T6SS sharing an MRCA with the KAC-specific T6SS. This indicates that the T6SS regions of Pluralibacter and KAC have a common origin, but the exact origin remains to be identified. The motility of K. aerogenes, another hallmark within the typically non-motile Klebsiella genus,40 may be attributed to its unique flagella operon. This motility, often associated with chemotaxis-driven swimming in several pathogenic bacteria, plays crucial roles in colonization and infection.41 The similarity of KAC flagella genes to those in Serratia and its sister genera suggests a possible horizontal gene transfer of the flagella operon from Serratia. Notably, this is firstly identified in a previous study of a KAC clinical isolate.42 In our study, we extend these results and significantly improve robustness by comparing all available KAC genomes against 325,562 RefSeq genomes in the public database. The flagellar region identified in our study is bounded by flhD and fliA, the same as have been previously reported,42 but the exact length of this region varies, due to notable genetic diversities between flgL and fliR, and fliE and fliT. However, identifying a specific donor species within Serratia is challenging, due to the close phylogenomic distances among its members. Intriguingly, the phylogenomic positioning of KAC flagella amid Serratia species implies that the horizontal gene transfer event might have occurred during the early stages of Serratia’s speciation.

We provide new genomic evidence that KAC is an important pathogen able to spread in and between healthcare institutions in view of the identification of inter-patient transmission and the presence of globally disseminated CC1 (ST93), CC2 (ST4), and CC3 (ST93) lineages. The global dissemination of ST4 and ST93 has also been highlighted in a recently published analysis of 565 KAC genomes.43 Nevertheless, among our 130 cases with KAC, we only identified two events of inter-patient transmission, involving two patients for each event. The relatively less recovery from clinical samples and the limited in-hospital transmission suggest that unlike its close relative K. pneumoniae, KAC has not yet evolved to become a major problem in healthcare settings at present. Further, by analyzing publicly available genomes, we identified clonal transmission only in individual countries without international transmission of a common clone at present. As such, it is not too late to pay close attention to this important pathogen by rigorous surveillance, detailed comparison of its genomic features with other Klebsiella species, and enhanced studies on its pathogenicity and antimicrobial resistance.

We developed a cgMLST scheme for facilitating future studies of KAC as a more accessible alternative of the SNP-based phylogenomic analysis, which requires advanced expertise in bioinformatics and bacterial biology and is computationally demanding and labor intensive.44 Notably, the genome’s diverse evolutionary pressures lead to an uneven distribution of allele numbers across various loci,45,46 influenced by both negative and positive selection. Certain loci that demonstrate excessive conservativeness, manifesting minimal variation among clusters, along with those exhibiting high diversity within a single cluster due to positive selection, can potentially compromise the accuracy of cgMLST. To address this, we applied a two-tiered process consisting of initially filtering loci that could undermine the accuracy, then optimizing the cutoff for precise CC assignment. This methodical approach ensures equilibrium between the inclusion of significant genetic markers and the exclusion of loci that might lead to inaccuracies via exaggeration of the true genetic distance between strains, thereby enhancing the reliability and precision of cgMLST. Using such an approach, we established a cgMLST scheme for KAC that appears to be robust based on a large set of genomes and rigorous quality control for genomes and genes. This cgMLST scheme exhibits a high degree of precision and reliability in strain categorization and therefore provides a useful tool for studies of evolution, clonal structure, and strain transmission of the complex.

As already mentioned, KAC has an intrinsic AmpC-encoding gene, and its high-level expression mediates resistance to most β-lactams13 but spares cefepime, carbapenems, cefiderocol, and agents containing newer β-lactamase inhibitors such as ceftazidime-avibactam.47 Notably, a proportion of KAC isolates also carries genes encoding ESBL, expanding the resistance spectrum to cefepime. Conversely, the production of AmpC leads to the failure of phenotypical detection of ESBL using synergic tests of a cephalosporin (commonly cefotaxime) and a classical β-lactamase inhibitor (e.g., clavulanic acid).48 Nevertheless, routine ESBL testing is not recommended by the Clinical and Laboratory Standards Institute.49 Carbapenems are stable to AmpC and ESBLs; however, some KAC isolates have carbapenemase-encoding genes. For these isolates, cefiderocol or agents containing newer β-lactamase inhibitors such as ceftazidime-avibactam (in combination with aztreonam for MBL-producing isolates) typically remain active.50

Limitations of the study

We are aware of limitations of this study. First, we studied a collection of KAC isolates recovered from our hospital with incomplete epidemiological investigation data such as the lack of environmental sampling to uncover whether transmission of a common strain between patients was due to direct inter-patient interaction or indirectly via a contaminated source. Second, we did not perform phenotypic testing for ESBLs, AmpC, and carbapenemases and therefore were unable to investigate their phenotype-genotype link. Third, we analyzed publicly available genomes, while the sampling and subsequent genome sequencing are largely biased and may not provide comprehensive coverage to study the population structure of KAC. Some regions may be under-represented, potentially affecting the observed distribution patterns. Fourth, we used the mutation rate calculated from CC2 of phylogroup 1, but such a rate may vary among different lineages. However, the numerous recombination events in the KAC genome disrupted the temporal signal, not allowing for accurate dating of the phylogeny for the entire complex. CC2 was the largest cluster that could be dated with a sensible length of MCMC, providing the best available approximation of the annual SNP accumulation rate for KAC genomes. Fifth, we did not perform experiments to verify the exact function of the two KAC-specific genetic regions, which is beyond the scope of this study. Despite the limitations, our analyses and findings provide much-needed novel insights into the population structure, distribution, transmission, antimicrobial resistance, and virulence of K. aerogenes and highlight that this species complex has remarkable genomic plasticity with specific features contributing to its virulence.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Bacterial and viral strains	
	
Klebsiella aerogenes clinical isolates	This manuscript	GenBank: PRJNA922913	
K. aerogenes KCTC 2190T	NCBI	GenBank: CP002824.1	
	
Critical commercial assays	
	
VITEK II	bioMérieux	N/A	
	
Deposited data	
	
130 K aerogenes genome assemblies	This manuscript	GenBank: PRJNA922913	
cgMLST scheme of K. aerogenes complex	This manuscript	https://github.com/fengyuchengdu/KAC_cgMLST_scheme (https://doi.org/10.5281/zenodo.12749494)	
cps references of K. aerogenes complex	This manuscript	https://github.com/fengyuchengdu/KAC_cps_reference (https://doi.org/10.5281/zenodo.12749512)	
264 K aerogenes genome assemblies	NCBI	NCBI RefSeq database	
1,829 K aerogenes raw sequencing reads	NCBI	NCBI Sequence Read Archive (SRA)	
17,681 non-K. aerogenes genome assemblies	NCBI	NCBI RefSeq database	
325,562 available bacterial genome assemblies	NCBI	NCBI RefSeq database	
	
Software and algorithms	
	
SPAdes (v3.15.3)	GitHub	https://github.com/ablab/spades	
MLST (v2.23.0)	GitHub	https://github.com/tseemann/mlst	
AMRFinderPlus (v3.10.23)	GitHub	https://github.com/ncbi/amr	
Kaptive (v3.0.0)	GitHub	https://github.com/klebgenomics/Kaptive	
Kleborate (v2.4.1)	GitHub	https://github.com/klebgenomics/Kleborate	
Snippy (v4.6.0)	GitHub	https://github.com/tseemann/snippy	
SeqKit (v2.8.0)	GitHub	https://github.com/shenwei356/seqkit	
FastANI (v1.33)	GitHub	https://github.com/ParBLiSS/FastANI	
CheckM (v1.2.2)	GitHub	https://github.com/Ecogenomics/CheckM	
BLAST	NCBI	https://blast.ncbi.nlm.nih.gov	
Gubbins (v3.3.5)	GitHub	https://github.com/nickjcroucher/gubbins	
fastbaps (v1.0.8)	GitHub	https://github.com/gtonkinhill/fastbaps	
BactDating (v1.1.2)	GitHub	https://github.com/xavierdidelot/BactDating	
chewBBACA (v3.3.2)	GitHub	https://github.com/B-UMMI/chewBBACA	
PIRATE (v1.0.5)	GitHub	https://github.com/SionBayliss/PIRATE	
Panstripe (v0.3.0)	GitHub	https://github.com/gtonkinhill/panstripe	
HAllA (v0.8.20)	GitHub	https://github.com/biobakery/halla	
Assembly Dereplicator (v0.3.2)	GitHub	https://github.com/rrwick/Assembly-Dereplicator	
Prokka (v1.14.5)	GitHub	https://github.com/tseemann/prokka	
dplyr (v.1.14)	GitHub	https://github.com/tidyverse/dplyr	
clusterProfiler (v4.10)	GitHub	https://github.com/YuLab-SMU/clusterProfiler	
GTDB-Tk (v2.3.2)	GitHub	https://github.com/Ecogenomics/GTDBTk	
SecReT6 (v3)	SJTU	http://db-mml.sjtu.edu.cn/SecReT6/	
MAFFT (v7.5207)	GitHub	https://github.com/Schaudge/MAFFT	
IQ-TREE (v2.2.68)	GitHub	https://github.com/Cibiv/IQ-TREE	
R package	CRAN	https://www.R-project.org	
ggplot2 (v3.4.4)	CRAN	https://cran.r-project.org/web/packages/ggplot2/index.html	
ggtree (v3.10)	GitHub	https://github.com/YuLab-SMU/ggtree	
gggenomes (v0.9.12.9000)	GitHub	https://github.com/thackl/gggenomes	

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Zhiyong Zong (zongzhiy@scu.edu.cn).

Materials availability

Strains used in this study will be shared by the lead contact upon request.

Data and code availability

• The genome assemblies included in this analysis have been deposited at the National Center of Biotechnology Information (NCBI) and are available as of date of publication. BioProject accession numbers are listed in the key resources table. The cgMLST scheme of K. aerogenes complex generated in this study has been deposited to https://github.com/fengyuchengdu/KAC_cgMLST_scheme (https://doi.org/10.5281/zenodo.12749494). The cps references of K. aerogenes complex defined in this strudy has been deposited to https://github.com/fengyuchengdu/KAC_cps_reference (https://doi.org/10.5281/zenodo.12749512).

• All original code used to produce the assembly, analyses, results, and visualization have been deposited at GitHub and is publicly available as of the date of publication.

• Any additional information required to reanalyze the data reported in this work paper is available from the lead contact upon request.

Experimental model and study participant details

Bacterial strains

All K. aerogenes clinical isolates collected between 2018 and 2022 as part of routine care in West China Hospital were included. Preliminary species identification and in vitro antimicrobial susceptibility testing were performed using Vitek II (bioMérieux, Marcy-l'Étoile, France). Primary and secondary bloodstream infection was defined according to the Bloodstream Infection Event protocol of National Healthcare Safety Network (NHSN, https://www.cdc.gov/nhsn/psc/bsi/index.html).

Ethics approval

This study was approved by the Ethical Committee of West China Hospital with informed consent being waived. Patients were consented with the care provided by the hospital and its staff, which included sample collections for managing their infection or suspected infection. All patient data was anonymized.

Method details

Whole genome sequencing

All of the 130 isolates were subjected to whole genome sequencing using HiSeq X10 platform (Illumina, San Diego, CA, USA) with reads assembled in contigs using SPAdes v3.15.351 under the isolate mode. Isolates were assigned to sequence types (ST) using MLST 2.23.052 to query pubmlst (https://pubmlst.org/organisms/klebsiella-aerogenes), were screened for β-lactamase-encoding genes using AMRFinderPlus v3.10.23,53 and were subjected to identification of virulence factors using Kleborate v2.4.1.54 Single nucleotide polymorphisms (SNPs) between isolates of the same STs were determined using Snippy v4.6.0 (https://github.com/tseemann/snippy) with the reference randomly selected from each ST.

Phylogenomic analysis

All included genomes (n = 2,093), each with a unique BioSample number, were subjected to a series of quality control and genome profiling measures. All bioinformatic settings and thresholds were set to the default unless otherwise specified. Among the 2,093 genomes, 264 were assembled and 1,829 had only short reads without an assembly available. Therefore, short reads of the 1,829 genomes were retrieved from SRA of NCBI and were trimmed to remove adapters and reads shorter than 30 bp using Trimmomatic v.0.39.55 For each of the 1,829 genome, maximum coverage of 150 × sequencing depth was subsampled using SeqKit v2.8.056 and assembled into a draft genome using SPAdes v3.15.3.51

The precise species identification for each genome was performed using FastANI v1.3357 by comparing genomes with the type strain of each species from the genus Klebsiella. The completeness and heterogeneity of genomes were assessed using CheckM v1.2.2.58 STs were determined using MLST 2.23.052 to query pubmlst. Antimicrobial resistance genes were identified using AMRFinderPlus v3.10.23.53 Considering that the virulence factors of KAC have not been well defined, we screened the 1,156 genomes for known virulence factors of K. pneumoniae (n = 80, https://bigsdb.pasteur.fr/cgi-bin/bigsdb/bigsdb.pl?db=pubmlst_klebsiella_seqdef&page=downloadAlleles&scheme_id=4&render=1)59 using BLASTn with a minimum coverage of 70%.

Generating the true phylogeny and clustering isolates into clonal clusters

The study analyzed 1,156 isolates, aligning their assembled genomes to a reference chromosome for pseudo-multiple alignment using Snippy, which generates pseudo-reads from assemblies and maps them onto the reference. This process facilitated the development of a whole-genome-based phylogeny, representing the true homology of the isolates. The analysis accounted for recombination sites by excluding them during the phylogenetic inference, which was performed using the maximum-likelihood method with the GTR-GAMMA model in Gubbins v3.3.5.60 Cluster partitioning across three levels was executed using fastbaps v1.0.8,34 leveraging the custom-inferred tree and the prior "optimized-baps". To strike a balance between enhanced resolution and avoiding over-partitioning, clonal clusters (CCs) were delineated based on the second-level results, as illustrated in Figure S1

Determination of SNP cutoffs to define clones

The study estimated the average nucleotide substitution rate in a target group utilizing the Bayesian approach, as implemented in BactDating v1.1.2,61 with the “mixedcarc” model and 107 iterations of the Markov Chain Monte Carlo (MCMC). Briefly, a phylogenomic tree was constructed, integrating the recombination model through Gubbins v3.3.5.60 This tree was then used for dating nodes based on the reported sampling time of each isolate. Results showing an effective sample size (ESS) of at least 200 were considered reliable. We defined the maximum allowed distance between two isolates, allowing them to be considered gemonically linked, as the number of differences they could accrue within a defined period. This study applied varying time intervals to modulate the stringency of the analysis: for isolates collected locally, a maximum duration was set to one season, reflecting the robust sampling rate; conversely, for analysis incorporating isolates from the public database, we extended this period to a year, accommodating potential underestimations of true phylogenetic relationships due to lower sampling rates. This gave us maximum differences of 6 and 22 nucleotide substitutions as the thresholds to determine epidemiological relatedness for analyzing locally collected isolates and for isolates in the clonal clusters, respectively. Isolates were subsequently clustered based on the abovementioned thresholds. The inclusion of an isolate in a specific group did not necessarily require satisfying the pairwise distance criterion. Instead, an isolate was considered part of a group as long as it was epidemiologically linked to at least one other isolate within that same group.

Development of a robust cgMLST scheme for KAC typing

The cgMLST scheme for KAC was developed using the chewBBACA v3.3.262 protocol (https://chewbbaca.readthedocs.io/en/latest/user/tutorials/chewie_step_by_step.html), with a meticulous approach to ensure accuracy. Initially, only isolates with complete chromosomes (n = 37) were included to mitigate the risk of pseudo-genes arising from mis-assemblies and to exclude elements associated with plasmids. Additionally, paralogous genes were identified and removed. A conservative threshold of 95% was applied to exclude non-core loci at this stage. Following this, alleles for the remaining isolates (n = 1,119) were called. With the addition of new isolates to the study, the threshold for locus exclusion was raised to 99%, culminating in a final scheme comprising 3,430 loci. This progressive refinement ensures the cgMLST scheme remains robust and representative of KAC.

Next, we applied a two-tiered process comprising filtering of loci that could undermine the accuracy and subsequent optimizing the cutoff for precise clonal clusters (CC) assignment. To evaluate which loci should be excluded from the cgMLST scheme, we conducted a comprehensive evaluation encompassing three key aspects. These include Cluster-Wide Conservativeness Index, a metric defined as the ratio of the number of clusters sharing the same allele to the total number of clusters aiming to capturing loci that are too conservative across different clusters, Phylogroup-Wide Conservativeness Index, a metric measuring the proportion of phylogroups sharing the same allele relative to the total number of phylogroups, providing insight into the conservativeness of loci on a broader phylogenetic scale, and Inter-Intra Clusters Variation Index, a metric calculated as the adjusted p-values (Benjamini-Hochberg) when the mean difference rate in allele variation between clusters exceeds that within clusters by a specified threshold, facilitating in locating loci that show significant variation within clusters.

For each locus, these indices were calculated, and loci were filtered based on a combination of serial thresholds applied to these indices. These thresholds varied from 0 to 0.3 for the cluster-wide conservativeness index, from 0 to 0.5 for the phylogroup-wide conservativeness index, and from 0.1 to 0.9 for the inter-intra clusters variation index. Applying these criteria, alongside controls for each index, resulted in the generation of 350 distinct cgMLST profiles, each comprising a varying number of loci for subsequent analysis and evaluation. This methodical approach ensured a rigorous selection of loci, aiming to optimize the accuracy and reliability of the cgMLST scheme.

The refinement of schemes and cutoffs involved a meticulous comparative analysis between cgMLST-derived clustering and phylogeny-based clustering for isolate pairs (hereafter designated as the inferred and the true clustering, respectively). The evaluation leveraged two principal statistical metrics, the Fowlkes–Mallows Index (FM),32 which measures the congruence between the clustering methods and converges on 1 as the inferred clustering approaches the true clustering, and the Balanced Accuracy (BA),33 which quantifies accuracy by assigning equal importance to sensitivity and specificity. These metrics were derived from four foundational parameters, which were adapted from those defined previsouly34: True Positives (TP), denoting the count of isolate pairs that co-clustered in both inferred and true clustering; True Negatives (TN), pertaining to the count of isolates forming singular clusters in both methods; False Positives (FP), representing the isolate pairs that co-clustered in the inferred clustering but not in the true clustering; and False Negatives (FN), indicating the isolate pairs that co-clustered in the true clustering but were not matched in the inferred ones. This analytical framework facilitates a robust assessment of the performance of each cgMLST scheme, ensuring the selection of an optimal cutoff for accurate strain classification.

The iterative refinement of cutoff values commenced with pairwise cgMLST distances for each scheme. As indicated above, we compared 350 schemes with different combinations of loci, generated by applying serial thresholds to the three indices of the original cgMLST scheme. We then iteratively tested the cutoff between 1/10 to 1/2 of the total loci for each profile. The criterion for determining the optimal scheme-cutoff combination incorporated the highest FM, BA, and the inclusivity of loci. For each profile and cutoff combination, we evaluated the FM and BA, looking for the best pair of profile scheme and cutoff that maximized both metrics, with a preference for the one with more alleles if these two metrics were the same. The chosen combination, which exhibited superior concordance between the inferred and true clustering, was subsequently selected as the definitive scheme. This stringent methodology underpins the selection of a cutoff (the number of allelic differences for defining two strains belonging to the same CC) that most accurately reflects the genuine phylogenetic relationships among the isolates, thereby increasing the reliability of the cgMLST approach in discerning strain relatedness.

The selection of loci for the cgMLST profiles was guided by a comprehensive analysis utilizing three distinct indices, with the inter-intra clusters variation index proving to be the most crucial among the indices (FM and BA as the other two) in determining loci inclusion. An optimal threshold of 0.9 for this index facilitated the development of a scheme consisting of 2,067 loci. We proposed that this scheme in conjunction with a cutoff value of 432, could robustly classify KAC isolates into established clonal groups. The cgMLST scheme has been deposited at https://github.com/fengyuchengdu/KAC_cgMLST_scheme with the required dataset for predefined CCs available as Data S17.

Analysis for discordance between STs and CCs

Pseudo whole-genome alignments of ST4 and ST93 strains were generated as described above, using the complete chromosomes of strain Y1 (GenBank: CP045870.1) and strain G7 (GenBank: CP011539.1) as references. Ancestral recombination, defined as that affecting all isolates in a lineage, were identified using fastGEAR30 with the predefined partitioning of the population into clonal CCs. The identified recombination blocks were plotted alongside a phylogenetic tree using ggplot2 v3.4.463 in R package.

Capsular typing for KAC and association tests with other typing methods

The sequences of capsular polysaccharide synthesis (cps) loci of KAC species were extracted through the following steps.64 Genome annotation of all strains was performed using Prokka v1.14.5.65 The genomic region between the two hallmark genes, galF and ugd, which flank the cps region, was extracted if both genes were present on the same DNA segment. Coding proteins from the identified cps regions were clustered using PIRATE v1.0.5,66 and cps regions were grouped based on gene patterns. For each group, the sequence with NCBI RefSeq annotations representing the majority of the genetic landscape was selected as the reference. These sequences, representing distinct cps types in KAC, were iteratively evaluated by being as references for typing with Kaptive v3.0.0.67 Sequences with no missing of expected genes, a length difference <100 bp, an identity of ≥98%, and a coverage of ≥95% were assigned the same KL types. Others were treated as variants of the primary KL unless typing was hindered by the low quality of the query genomes, resulting in untypeable or lower confidence. The relative distances of primary KLs were calculated using Mashtree, and associations between various groupings and KLs were assessed using the chi-square test with Cramér’s V to determine the strength of the associations using R packages. The genomic synteny map of primary KLs were visualized using clinker v0.0.28.68 All reference sequences of KLs have been deposited at https://github.com/fengyuchengdu/KAC_cps_reference.

Analysis for open or close genomes

A gene presence-absence matrix was constructed using PIRATE v1.0.5 with default settings. The evolutionary dynamics of gene exchange through horizontal gene transfer across 1,156 KAC genomes were examined via Panstripe v0.3.0,69 under the assumption of a Tweedie distribution and based on the phylogeny inferred by Gubbins. This analysis was also performed to each phylogroup, using the group-specific gene presence-absence matrix and the subtree tailored to each group, derived from tree partitioning. Additionally, the resulting gene presence-absence matrix was analyzed using Heaps’ law model in R to determine the key parameter α, which was used to infer the openness of the genome through the saturation curve in a traditional manner. Hierarchical association tests were performed using HAllA v0.8.20 (https://github.com/biobakery/halla) between pre-defined groups (phylogroups, clonal clusters, sampling locations and continents) and common characteristic determinants (antimicrobial resistance, virulence, stress tolerance and plasmid replicons).

Comparative genomics

From an initial collection of 1,156 genomes, a total of 746 assemblies were selected to represent the pan-genome of KAC, using Assembly Dereplicator v0.3.2 (https://github.com/rrwick/Assembly-Dereplicator) with a distance cutoff of 0.001. Additionally, a separate dataset comprising 17,681 non-KAC assemblies was retrieved from NCBI RefSeq database. This dataset underwent a stringent quality assessment, identical to the one applied to KAC assemblies and survived with 11,076 genomes. Dereplication was then performed and yielded 5,030 genomes representing the pan-genome of non-KAC.

The comparative genomics analysis incorporated a total of 5,776 assemblies, comprising 746 KAC and 5,030 non-KAC genomes, maintaining an approximate ratio of 1:6. The analysis commenced with the annotation of all genomes using Prokka v1.14.5.65 This was followed by the clustering of coding sequences with PIRATE v1.0.5. A preliminary definition for a KAC-specific gene cluster was established as alleles that were present in 99% of KAC genomes and completely absent in non-KAC genomes. To validate the uniqueness of these initially identified KAC-specific alleles, a comprehensive pairwise comparison of their amino acid sequences across all 5,776 genomes was conducted using BLAST+ v.2.15.0.70 The results were summarized and filtered using R package dplyr v.1.14 (https://github.com/tidyverse/dplyr). Ultimately, a gene was classified as KAC-specific only if the minimum identity within KAC genomes was at least 50% higher than the maximum identity observed between KAC and non-KAC genomes.

Pathway enrichment

Pathway enrichment analysis was conducted focusing on K. aerogenes KCTC 2190T (GenBank: CP002824.1), which serves as the reference model in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. For this purpose, the identified KAC-specific genes of this isolate were subjected to KEGG pathway annotation, and statistical tests of over-representation were performed using clusterProfiler v4.10.71

Similarity searching in NCBI RefSeq database

The similarity searching process commenced by retrieving all available assemblies from the NCBI RefSeq database (n = 325,562 as of 01-12-2023). Taxonomy assignments for these genomes were performed using GTDB-Tk v2.3.272 with database 08-RS214 released on 28-04-2023. Notably, Raoultella was included in Klebsiella in GTDB. A total of 557 assemblies, deemed low-quality according to criteria outlined on the GTDB website (https://gtdb.ecogenomic.org/faq), were discarded, resulting in a total of 325,005 queries for the downstream BLAST analysis. The individual overall coverage and identity of target regions were summarized using custom python scripts, where the former was calculated by first concatenating the segments of the subject with a minimum of 100-bp alignment and overlaps accounted for, followed by dividing the length of this concatenated sequence by the total length of the subject, and the later was calculated by a weighted average based on alignment length. Ultimately, a query genome was considered to possess the subject sequence if its overall coverage in the query genome exceeded 50%.

T6SS prediction

Type VI secretion system (T6SS) gene clusters were predicted and classified based on BLASTP algorithm against experimentally validated records in SecReT6 v3.73 We also extracted all records of T6SS typing from complete bacterial genomes as curated by SecReT6 v3.73 Our analysis focused on evaluating the prevalence of different T6SS types within the genus Klebsiella and the highly related genera such as Enterobacter and Pluralibacter, which were identified based on their similarity to KAC-specific T6SS. To ensure taxonomic accuracy, we applied the ANI method, using a 95% cutoff for species delineation.

Homology analysis of full-length flagella and T6SS regions

The homology analysis focused on species classified by GTDB that exhibited an average coverage of 80% or higher for the KAC-specific T6SS region, and 60% or higher for the KAC-specific flagella region. Assemblies were initially clustered by species and subsequently dereplicated using Assembly Dereplicator with a distance cutoff of 0.03. This step effectively minimized genome redundancy while preserving species diversity. The KAC-specific flagella and T6SS of K. aerogenes KCTC 2190T served as the references, and their start and endpoints on the chromosome were manually adjusted, if necessary, to encompass the entire operon. These corrected sequences were used to identify the equivalent regions in each representative assembly, using SeqKit v2.8.056 with coordinates derived from BLASTn.

After annotation, the coding sequences were clustered into 67 gene clusters for flagella and 24 for T6SS. These clusters were then aligned by their amino acid sequences using MAFFT v7.5207.74 A phylogenetic tree was constructed from these concatenated aligned peptides, employing the LG model with rate heterogeneity and 1,000 ultrafast bootstraps, using IQ-TREE v2.2.68.75 Tree rooting was conducted using a non-reversible model for amino acids and the similarity index was defined as the ratio of the BLASTP score of the query protein to the maximum score in an identical match. Visual representations, including trees, heatmaps, and gene context synteny, were generated using R packages ggplot2 v3.4.4,63 ggtree v3.1076 and gggenomes v0.9.12.9000.77

Biochemical assays for isolates of each phylogroup

We tested 18 KAC isolates comprising all of phylogroups 3 (n = 4) and 4 (n = 5) and randomly selected ones of phylogroups 1 (n = 4) and 2 (n = 5) in our collection for enzyme production and metabolic capabilities on multiple carbon sources. The tests were performed using API 20E and API 50CHB kits (bioMérieux), following the manufacturer’s instructions and E. coli ATCC 25922 was used as the reference.

Quantification and statistical analysis

In the root-to-tip analysis, temporal signals were assessed using permutation tests on CC1, CC2, and CC3, which included 164, 63, and 63 strains, respectively. Each permutation test involved 10,000 permutations to compute the p-value, and R2 was calculated by fitting a linear regression model. These analyses were performed using the R package BactDating. For the development of cgMLST, Wilcoxon tests with the alternative hypothesis "greater" were used to compare two groups, determining if one group was significantly greater than the other by a specified value, with p-values adjusted using continuity correction in R. Groups here were defined by the combinations of grouping methods and alleles (e.g., clusters-alleles, phylogroups-alleles). To examine the association between 24 KL types and three grouping methods (365 STs, 213 CCs, and five phylogroups), Chi-squared tests were performed, with p-values calculated via Monte Carlo simulation. The strength of association was also measured using Cramér’s V, calculated with the R package vcd. In pan-genome analysis, p-values and coefficients for number of gene gain/loss events versus branch distance were calculated using the R package Panstripe, where coefficients represented the slope of gene gain/loss per branch unit under a Tweedie model. For saturation curve analysis, 1,156 strains with a total of 33,332 loci were fitted to Heaps’ law model using the R package minpack.lm. To test associations between 829 genetic traits (e.g., antimicrobial resistance, virulence factors, and plasmid replicons) and 310 groupings (e.g., CCs, STs, and geographical locations), hierarchical all-against-all association testing was conducted using the program HAIIA, employing the NMI metric and Benjamini-Hochberg correction for p-values. For gene enrichment analysis, 48 genes involved in KAC pathways were tested for pathway enrichment using Fisher’s exact test with the R package clusterProfiler, and p-values were adjusted using Benjamini-Hochberg correction. All settings not specifically mentioned were kept as default. All p-values were corrected for multiple testing where applicable, and a p-value of <0.05 was considered significant for all tests used in this study.

Supplemental information

Document S1. Figures S1–S15 and Tables S1–S3

Data S1. The 130 KAC isolates in this study

Data S2. Pairwise SNPs between isolates of the same ST within the 130 KAC isolates in this study

Data S3. The 1,156 KAC genomes

Data S4. ST and CC assignations for the 1,156 KAC genomes

This dataset comprises five spreadsheets to demonstrate: (1) the number of isolates for each ST and CC; (2) the ST and CC assignation for each isolate; (3) the matches between ST and CC; (4) the country distribution of each CC; (5) the continent distribution of each CC.

Data S5. Clone assignation within CC1, CC2, and CC3 and pairwise SNPs

This dataset comprises five spreadsheets to show: (1) the number of isolates and the pairwise SNP range for each clone; (2) the country distribution and the carriage of genes encoding carbapenemases or ESBLs for each clone; (3–5) numbers of pairwise SNPs within CC1, CC2, and CC3, respectively.

Data S6. Capsular typing for all 1,156 KAC isolates

Data S7. The 24 capsular (KL) types identified in KAC and the pairwise identity and coverage

This dataset comprises two spreadsheets to demonstrate: (1) the summary of pairwise comparison of identity and coverage between each KL types; (2) the detailed comparison of identity and coverage between each KL types.

Data S8. The 196 AER AmpC β-lactamases

Data S9. Pairwise amino acid identities of AER AmpC β-lactamases

Data S10. Mutations of blaAERampC genes and amino acid substitutions of AER

This dataset comprises two spreadsheets to demonstrate: (1) the non-synonymous mutations of each blaAER gene; (2) all mutations of blaAER in all KAC genomes.

Data S11. Virulence genes on a plasmid in a KAC isolate (SRA: SRR19409809)

Data S12. Associations of KAC subpopulations with the functional traits predicted by the determinants of virulence, antimicrobial resistance, stress tolerance and plasmid replicon

This dataset comprises four spreadsheets to demonstrate: (1) association of all subpopulations as clonal clusters, country-specific groups, and phylogroups with virulence, antimicrobial resistance, stress tolerance, and plasmid replicon; (2) association of clonal clusters with virulence, antimicrobial resistance, stress tolerance, and plasmid replicon; (3) association of country-specific groups with virulence, antimicrobial resistance, stress tolerance, and plasmid replicon; (4) association of phylogroups with virulence, antimicrobial resistance, stress tolerance, and plasmid replicon.

Data S13. KAC-specific genes and pathways

This dataset comprises three spreadsheets to demonstrate: (1) 87 genes encoding products with >50% amino acid difference between KAC genomes and other Klebsiella ones; (2) 60 genes able to be assigned with a KEGG Orthology (KO) number; (3) 48 genes associated with at least one KAC pathway.

Data S14. The distribution of T6SS type and number in KAC and selected genera or complex of the Enterobacteriaceae

Data S15. The distribution of KAC-specific T6SS and its similar (coverage, >80%) ones in bacterial genera and species defined by GTDB

This dataset comprises two spreadsheets to demonstrate: (1) the distribution in bacterial genera; (2) the distribution in bacterial species.

Data S16. The distribution of KAC-specific flagella region and its similar (coverage, >60%) regions in bacterial genera and species defined by GTDB

This dataset comprises two spreadsheets to demonstrate: (1) the distribution in bacterial genera; (2) the distribution in bacterial species.

Data S17. The predefined KAC CCs used for cgMLST

Document S2. Article plus supplemental information

Acknowledgments

We are grateful to Ms. Yanling He for performing the biochemical tests for isolates of the four phylogroups. The work was supported by grants from the 10.13039/501100001809 National Natural Science Foundation of China (project nos. 81861138055 and 82372290 ); the 10.13039/501100012166 National Key Research and Development Program of China (grant nos. 2022YFC2303900 and 2023YFC2308800 ); the 10.13039/501100000265 Medical Research Council -UKRI (MR/S013660/1 ); 10.13039/501100018542 Natural Science Foundation of Sichuan Province (2023NSFSC1467 ); the 10.13039/501100013365 West China Hospital, Sichuan University (1.3.5 project for disciplines of excellence, grant nos. ZYYC08006 and ZYGD22001 ); and 10.13039/501100000272 National Institute for Health and Care Research (NIHR) 10.13039/501100018952 Birmingham Biomedical Research Centre (NIHR203326 ).

Author contributions

Z.Z. designed the study. Y. Xiao, Y. Xie, L.W., and H.W. collected the isolates and subjected them to genome sequencing. Y.F., Y.Y., Y.H., L.Z., A.M., and Z.Z. analyzed the data and genome sequences. Y.F., A.M., and Z.Z. wrote the paper.

Declaration of interests

The authors declare no competing interests.

Supplemental information can be found online at https://doi.org/10.1016/j.celrep.2024.114602.
==== Refs
References

1 Tindall B.J. Sutton G. Garrity G.M. Enterobacter aerogenes Hormaeche and Edwards 1960 (Approved Lists 1980) and Klebsiella mobilis Bascomb et al. 1971 (Approved Lists 1980) share the same nomenclatural type (ATCC 13048) on the Approved Lists and are homotypic synonyms, with consequences for the name Klebsiella mobilis Bascomb et al. 1971 (Approved Lists 1980) Int. J. Syst. Evol. Microbiol. 67 2017 502 504 27902205
2 Davin-Regli A. Lavigne J.P. Pages J.M. Enterobacter spp.: update on taxonomy, clinical aspects, and emerging antimicrobial resistance Clin. Microbiol. Rev. 32 2019 00002 00019
3 Chavda K.D. Chen L. Fouts D.E. Sutton G. Brinkac L. Jenkins S.G. Bonomo R.A. Adams M.D. Kreiswirth B.N. Comprehensive genome analysis of carbapenemase-producing Enterobacter spp.: new insights into phylogeny, population structure, and resistance mechanisms mBio 7 2016 020933
4 Davin-Regli A. Pagès J.M. Enterobacter aerogenes and Enterobacter cloacae; versatile bacterial pathogens confronting antibiotic treatment Front. Microbiol. 6 2015 392 26042091
5 Georghiou P.R. Hamill R.J. Wright C.E. Versalovic J. Koeuth T. Watson D.A. Lupski J.R. Molecular epidemiology of infections due to Enterobacter aerogenes: identification of hospital outbreak-associated strains by molecular techniques Clin. Infect. Dis. 20 1995 84 94 7727676
6 De Gheldre Y. Maes N. Rost F. De Ryck R. Clevenbergh P. Vincent J.L. Struelens M.J. Molecular epidemiology of an outbreak of multidrug-resistant Enterobacter aerogenes infections and in vivo emergence of imipenem resistance J. Clin. Microbiol. 35 1997 152 160 8968898
7 Wesevich A. Sutton G. Ruffin F. Park L.P. Fouts D.E. Fowler V.G. Jr. Thaden J.T. Newly Named Klebsiella aerogenes (formerly Enterobacter aerogenes) Is Associated with Poor Clinical Outcomes Relative to Other Enterobacter Species in Patients with Bloodstream Infection J. Clin. Microbiol. 58 2020 005822
8 Laupland K.B. Edwards F. Harris P.N.A. Paterson D.L. Significant clinical differences but not outcomes between Klebsiella aerogenes and Enterobacter cloacae bloodstream infections: a comparative cohort study Infection 51 2023 1445 1451 36881325
9 Song E.H. Park K.H. Jang E.Y. Lee E.J. Chong Y.P. Cho O.H. Kim S.H. Lee S.O. Sung H. Kim M.N. Comparison of the clinical and microbiologic characteristics of patients with Enterobacter cloacae and Enterobacter aerogenes bacteremia: a prospective observation study Diagn. Microbiol. Infect. Dis. 66 2010 436 440 20071128
10 Alvarez-Marin R. Lepe J.A. Gasch-Blasi O. Rodríguez-Martínez J.M. Calvo-Montes J. Clinical characteristics and outcome of bacteraemia caused by Enterobacter cloacae and Klebsiella aerogenes: more similarities than differences J Glob Antimicrob Resist 25 2021 351 358 33964492
11 Boucher H.W. Talbot G.H. Bradley J.S. Edwards J.E. Gilbert D. Rice L.B. Scheld M. Spellberg B. Bartlett J. Bad bugs, no drugs: no ESKAPE! An update from the Infectious Diseases Society of America Clin. Infect. Dis. 48 2009 1 12 19035777
12 Carroll K.C. Munson E. Butler-Wu S.M. Patrick S. Point-Counterpoint: What's in a Name? Clinical Microbiology Laboratories Should Use Nomenclature Based on Current Taxonomy J. Clin. Microbiol. 61 2023 e0173222
13 Preston K.E. Radomski C.C. Venezia R.A. Nucleotide sequence of the chromosomal ampC gene of Enterobacter aerogenes Antimicrob. Agents Chemother. 44 2000 3158 3162 11036041
14 Lee S.H. Jeong S.H. Park Y.M. Characterization of blaCMY-10 a novel, plasmid-encoded AmpC-type β-lactamase gene in a clinical isolate of Enterobacter aerogenes J. Appl. Microbiol. 95 2003 744 752 12969288
15 Lavigne J.P. Bouziges N. Chanal C. Mahamat A. Michaux-Charachon S. Sotto A. Molecular epidemiology of Enterobacteriaceae isolates producing extended-spectrum β-lactamases in a French hospital J. Clin. Microbiol. 42 2004 3805 3808 15297534
16 Dumarche P. De Champs C. Sirot D. Chanal C. Bonnet R. Sirot J. TEM derivative-producing Enterobacter aerogenes strains: dissemination of a prevalent clone Antimicrob. Agents Chemother. 46 2002 1128 1131 11897606
17 Biendo M. Canarelli B. Thomas D. Rousseau F. Hamdad F. Adjide C. Laurans G. Eb F. Successive emergence of extended-spectrum β-lactamase-producing and carbapenemase-producing Enterobacter aerogenes isolates in a university hospital J. Clin. Microbiol. 46 2008 1037 1044 18234876
18 Franolic I. Bedenić B. Beader N. Lukić-Grlić A. Mihaljević S. Bielen L. Zarfel G. Meštrović T. NDM-1-producing Enterobacter aerogenes isolated from a patient with a JJ ureteric stent in situ CEN Case Rep 8 2019 38 41 30141138
19 Li Y. Wang Q. Xiao X. Li R. Wang Z. Emergence of blaNDM-9-bearing tigecycline-resistant Klebsiella aerogenes of chicken origin J Glob Antimicrob Res 26 2021 66 68
20 da Silva K.E. de Almeida de Souza G.H. Moura Q. Rossato L. Limiere L.C. Vasconcelos N.G. Simionatto S. Genetic Diversity of Virulent Polymyxin-Resistant Klebsiella aerogenes Isolated from Intensive Care Units Antibiotics (Basel) 11 2022 1127 36009996
21 Kayama S. Yu L. Kawakami S. Yahara K. Hisatsune J. Yamamoto M. Yamamoto K. Shimono N. Kibe Y. Kiyosuke M. Sugai M. Emergence of blaNDM-5-Carrying Klebsiella aerogenes in Japan Microbiol. Spectr. 10 2022 e0222221
22 Vargas J.M. Moreno Mochi M.P. Nuñez J.M. Mochi S. Cáceres M. Del Campo R. Jure M.A. Emergence and clonal spread of KPC-2-producing clinical Klebsiella aerogenes isolates in a hospital from northwest Argentina J. Med. Microbiol. 72 2023 001635
23 Podschun R. Ullmann U. Klebsiella spp. as nosocomial pathogens: epidemiology, taxonomy, typing methods, and pathogenicity factors Clin. Microbiol. Rev. 11 1998 589 603 9767057
24 Passarelli-Araujo H. Palmeiro J.K. Moharana K.C. Pedrosa-Silva F. Dalla-Costa L.M. Venancio T.M. Genomic analysis unveils important aspects of population structure, virulence, and antimicrobial resistance in Klebsiella aerogenes FEBS J. 286 2019 3797 3810 31319017
25 Richter M. Rosselló-Móra R. Shifting the genomic gold standard for the prokaryotic species definition Proc. Natl. Acad. Sci. USA 106 2009 19126 19131 19855009
26 Rossello-Mora R. Amann R. Past and future species definitions for Bacteria and Archaea Syst. Appl. Microbiol. 38 2015 209 216 25747618
27 Ciufo S. Kannan S. Sharma S. Badretdin A. Clark K. Turner S. Brover S. Schoch C.L. Kimchi A. DiCuccio M. Using average nucleotide identity to improve taxonomic assignments in prokaryotic genomes at the NCBI Int. J. Syst. Evol. Microbiol. 68 2018 2386 2392 29792589
28 Farmer J.J. 3rd Davis B.R. Hickman-Brenner F.W. McWhorter A. Huntley-Carter G.P. Asbury M.A. Riddle C. Wathen-Grady H.G. Elias C. Fanning G.R. Biochemical identification of new species and biogroups of Enterobacteriaceae isolated from clinical specimens J. Clin. Microbiol. 21 1985 46 76 3881471
29 Davin-Regli A. Lavigne J.P. Pages J.M. Enterobacter spp.: Update on Taxonomy, Clinical Aspects, and Emerging Antimicrobial Resistance Clin. Microbiol. Rev. 32 2019 0002 0019
30 Mostowy R. Croucher N.J. Andam C.P. Corander J. Hanage W.P. Marttinen P. Efficient Inference of Recent and Ancestral Recombination within Bacterial Populations Mol. Biol. Evol. 34 2017 1167 1182 28199698
31 David S. Reuter S. Harris S.R. Glasner C. Feltwell T. Argimon S. Abudahab K. Goater R. Giani T. Errico G. Epidemic of carbapenem-resistant Klebsiella pneumoniae in Europe is driven by nosocomial spread Nat. Microbiol. 4 2019 1919 1929 31358985
32 Fowlkes E.B. Mallows C.L. A Method for Comparing Two Hierarchical Clusterings J. Am. Stat. Assoc. 78 1983 553 569
33 Owusu-Adjei M. Ben Hayfron-Acquah J. Frimpong T. Abdul-Salaam G. Imbalanced class distribution and performance evaluation metrics: A systematic review of prediction accuracy for determining model performance in healthcare systems PLOS Digit. Health 2 2023 e0000290
34 Tonkin-Hill G. Lees J.A. Bentley S.D. Frost S.D.W. Corander J. Fast hierarchical Bayesian analysis of population structure Nucleic Acids Res. 47 2019 5539 5549 31076776
35 Russo T.A. Marr C.M. Hypervirulent Klebsiella pneumoniae Clin. Microbiol. Rev. 32 2019 e00001-00019
36 Starkova P. Lazareva I. Avdeeva A. Sulian O. Likholetova D. Ageevets V. Lebedeva M. Gostev V. Sopova J. Sidorenko S. Emergence of Hybrid Resistance and Virulence Plasmids Harboring New Delhi Metallo-β-Lactamase in Klebsiella pneumoniae in Russia Antibiotics (Basel) 10 2021 691 34207702
37 Falcone M. Tiseo G. Arcari G. Leonildi A. Giordano C. Tempini S. Bibbolino G. Mozzo R. Barnini S. Carattoli A. Menichetti F. Spread of hypervirulent multidrug-resistant ST147 Klebsiella pneumoniae in patients with severe COVID-19: an observational study from Italy, 2020-21 J. Antimicrob. Chemother. 77 2022 1140 1145 35040981
38 Tettelin H. Riley D. Cattuto C. Medini D. Comparative genomics: the bacterial pan-genome Curr. Opin. Microbiol. 11 2008 472 477 19086349
39 Boyer F. Fichant G. Berthod J. Vandenbrouck Y. Attree I. Dissecting the bacterial type VI secretion system by a genome wide in silico analysis: what can be learned from available microbial genomic resources? BMC Genom. 10 2009 104
40 Carabarin-Lima A. León-Izurieta L. Rocha-Gracia R.D.C. Castañeda-Lucio M. Torres C. Gutiérrez-Cazarez Z. González-Posos S. Martínez de la Peña C.F. Martinez-Laguna Y. Lozano-Zarain P. First evidence of polar flagella in Klebsiella pneumoniae isolated from a patient with neonatal sepsis J. Med. Microbiol. 65 2016 729 737 27283194
41 Zhou B. Szymanski C.M. Baylink A. Bacterial chemotaxis in human diseases Trends Microbiol. 31 2023 453 467 36411201
42 Diene S.M. Merhej V. Henry M. El Filali A. Roux V. Robert C. Azza S. Gavory F. Barbe V. La Scola B. The rhizome of the multidrug-resistant Enterobacter aerogenes genome reveals how new "killer bugs" are created because of a sympatric lifestyle Mol. Biol. Evol. 30 2013 369 383 23071100
43 Morgado S. Fonseca É. Freitas F. Caldart R. Vicente A.C. In-depth analysis of Klebsiella aerogenes resistome, virulome and plasmidome worldwide Sci. Rep. 14 2024 6538 38503805
44 Wang Z. Gu C. Sun L. Zhao F. Fu Y. Di L. Zhang J. Zhuang H. Jiang S. Wang H. Development of a novel core genome MLST scheme for tracing multidrug resistant Staphylococcus capitis Nat. Commun. 13 2022 4254 35869070
45 Maiden M.C.J. Jansen van Rensburg M.J. Bray J.E. Earle S.G. Ford S.A. Jolley K.A. McCarthy N.D. MLST revisited: the gene-by-gene approach to bacterial genomics Nat. Rev. Microbiol. 11 2013 728 736 23979428
46 Didelot X. Maiden M.C.J. Impact of recombination on bacterial evolution Trends Microbiol. 18 2010 315 322 20452218
47 Philippon A. Arlet G. Labia R. Iorga B.I. Class C β-Lactamases: Molecular Characteristics Clin. Microbiol. Rev. 35 2022 e0015021
48 Drieux L. Brossier F. Sougakoff W. Jarlier V. Phenotypic detection of extended-spectrum β-lactamase production in Enterobacteriaceae: review and bench guide Clin. Microbiol. Infect. 14 2008 90 103 18154532
49 CLSI Performance Standards for Antimicrobial Susceptibility Testing; 34th Informational Supplement. M100-S34 2024 Clinical and Laboratory Standards Institute
50 Tamma P.D. Aitken S.L. Bonomo R.A. Mathers A.J. van Duin D. Clancy C.J. Infectious Diseases Society of America 2023 Guidance on the Treatment of Antimicrobial Resistant Gram-Negative Infections Clin. Infect. Dis. 2023 ciad428 10.1093/cid/ciad1428 37463564
51 Bankevich A. Nurk S. Antipov D. Gurevich A.A. Dvorkin M. Kulikov A.S. Lesin V.M. Nikolenko S.I. Pham S. Prjibelski A.D. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing J. Comput. Biol. 19 2012 455 477 22506599
52 Larsen M.V. Cosentino S. Rasmussen S. Friis C. Hasman H. Marvig R.L. Jelsbak L. Sicheritz-Pontén T. Ussery D.W. Aarestrup F.M. Lund O. Multilocus sequence typing of total-genome-sequenced bacteria J. Clin. Microbiol. 50 2012 1355 1361 22238442
53 Feldgarden M. Brover V. Gonzalez-Escalona N. Frye J.G. Haendiges J. Haft D.H. Hoffmann M. Pettengill J.B. Prasad A.B. Tillman G.E. AMRFinderPlus and the Reference Gene Catalog facilitate examination of the genomic links among antimicrobial resistance, stress response, and virulence Sci. Rep. 11 2021 12728
54 Lam M.M.C. Wick R.R. Watts S.C. Cerdeira L.T. Wyres K.L. Holt K.E. A genomic surveillance framework and genotyping tool for Klebsiella pneumoniae and its related species complex Nat. Commun. 12 2021 4188 34234121
55 Bolger A.M. Lohse M. Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data Bioinformatics 30 2014 2114 2120 24695404
56 Shen W. Le S. Li Y. Hu F. SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation PLoS One 11 2016 e0163962
57 Jain C. Rodriguez-R L.M. Phillippy A.M. Konstantinidis K.T. Aluru S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries Nat. Commun. 9 2018 5114 30504855
58 Parks D.H. Imelfort M. Skennerton C.T. Hugenholtz P. Tyson G.W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes Genome Res. 25 2015 1043 1055 25977477
59 Bialek-Davenet S. Criscuolo A. Ailloud F. Passet V. Jones L. Delannoy-Vieillard A.S. Garin B. Le Hello S. Arlet G. Nicolas-Chanoine M.H. Genomic definition of hypervirulent and multidrug-resistant Klebsiella pneumoniae clonal groups Emerg. Infect. Dis. 20 2014 1812 1820 25341126
60 Croucher N.J. Page A.J. Connor T.R. Delaney A.J. Keane J.A. Bentley S.D. Parkhill J. Harris S.R. Rapid phylogenetic analysis of large samples of recombinant bacterial whole genome sequences using Gubbins Nucleic Acids Res. 43 2015 e15 25414349
61 Didelot X. Croucher N.J. Bentley S.D. Harris S.R. Wilson D.J. Bayesian inference of ancestral dates on bacterial phylogenetic trees Nucleic Acids Res. 46 2018 e134
62 Silva M. Machado M.P. Silva D.N. Rossi M. Moran-Gilad J. Santos S. Ramirez M. Carriço J.A. chewBBACA: A complete suite for gene-by-gene schema creation and strain identification Microb. Genom. 4 2018 e000166
63 Wickham H. ggplot2: Elegant Graphics for Data Analysis 2016 Springer-Verlag New York
64 Holt K.E. Lassalle F. Wyres K.L. Wick R. Mostowy R.J. Diversity and evolution of surface polysaccharide synthesis loci in Enterobacteriales Isme j 14 2020 1713 1730 32249276
65 Seemann T. Prokka: rapid prokaryotic genome annotation Bioinformatics 30 2014 2068 2069 24642063
66 Bayliss S.C. Thorpe H.A. Coyle N.M. Sheppard S.K. Feil E.J. PIRATE: A fast and scalable pangenomics toolbox for clustering diverged orthologues in bacteria GigaScience 8 2019 giz119
67 Wick R.R. Heinz E. Holt K.E. Wyres K.L. Diekema D.J. Kaptive web: user-friendly capsule and lipopolysaccharide serotype prediction for Klebsiella genomes J. Clin. Microbiol. 56 2018 001977-18
68 Gilchrist C.L.M. Chooi Y.-H. clinker & clustermap.js: automatic generation of gene cluster comparison figures Bioinformatics 37 2021 2473 2475 33459763
69 Tonkin-Hill G. Gladstone R.A. Pöntinen A.K. Arredondo-Alonso S. Bentley S.D. Corander J. Robust analysis of prokaryotic pangenome gene gain and loss rates with Panstripe Genome Res. 33 2023 129 140 36669850
70 Camacho C. Coulouris G. Avagyan V. Ma N. Papadopoulos J. Bealer K. Madden T.L. BLAST+: architecture and applications BMC Bioinf. 10 2009 421
71 Wu T. Hu E. Xu S. Chen M. Guo P. Dai Z. Feng T. Zhou L. Tang W. Zhan L. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data Innovation 2 2021 100141
72 Chaumeil P.A. Mussig A.J. Hugenholtz P. Parks D.H. GTDB-Tk v2: memory friendly classification with the genome taxonomy database Bioinformatics 38 2022 5315 5316 36218463
73 Zhang J. Guan J. Wang M. Li G. Djordjevic M. Tai C. Wang H. Deng Z. Chen Z. Ou H.Y. SecReT6 update: a comprehensive resource of bacterial Type VI Secretion Systems Sci. China Life Sci. 66 2023 626 634 36346548
74 Katoh K. Misawa K. Kuma K.i. Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform Nucleic Acids Res. 30 2002 3059 3066 12136088
75 Minh B.Q. Schmidt H.A. Chernomor O. Schrempf D. Woodhams M.D. von Haeseler A. Lanfear R. IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era Mol. Biol. Evol. 37 2020 1530 1534 32011700
76 Xu S. Li L. Luo X. Chen M. Tang W. Zhan L. Dai Z. Lam T.T. Guan Y. Yu G. Ggtree: A serialized data object for visualization of a phylogenetic tree and annotation data iMeta 1 2022 e56 38867905
77 Hackl T. Ankenbrand M.J. van Adrichem B. Gggenomes: A Grammar of Graphics for Comparative Genomics 2023 https://github.com/thackl/gggenomes
