We performed a retrospective study by collecting non-repeated K. pneumoniae isolates from patients with invasive infections at a hospital specializing in liver disease (Beijing, China) from May 2018 to August 2021. K. pneumoniae isolates were obtained from cultures of blood, ascites, abscess drainage, biliary tract fluid, pleural fluid, and bronchoalveolar lavage aspirate, and identified by an automated Vitek II system (bioMerieux, Balmes-Les-Grottes, France). Afterwards, whole genome sequencing-based species identification with assembled genomes was performed using Kleborate v2.2.0 (Lam et al., 2021). Clinical information was retrieved from electronic medical records, including basic demographic characteristics, underlying diseases, antimicrobial agent exposures, the site of infection, and the use of invasive devices. Hospital-acquired infections were defined as incident infections with an onset after 48 hours of admission.

To determine the HMV phenotype, a string test was performed by touching a colony grown overnight on a blood agar plate at 37°C with a loop and pulling it up. Strains exhibiting a mucoid string with a length of ≥5 mm were defined as HmKp. Otherwise, they were classified as non-HmKp. Each string test was repeated three to five times.

Antimicrobial susceptibility testing was performed by the VITEK 2 (bioMérieux) with Vitek2 Gram Negative ID cards (AST-N335; bioMérieux) (ticarcillin, piperacillin, ceftazidime, cefoperazone, cefepime, aztreonam, imipenem, meropenem, amikacin, tobramycin, ciprofloxacin, levofloxacin, doxycycline, minocycline, trimethoprim) and broth dilution susceptibility testing (ertapenem). K. pneumoniae ATCC700603 was used as a quality control strain. Results were interpreted according to the recommendations of the Clinical and Laboratory Standards Institute 2020 (CLSI, 2020).

Genomic DNA was extracted using a QIAamp DNA Mini Kit (#51306 Qiagen) following the manufacturer’s recommendations. The Illumina nova seq6000 and MGISEQ-2000 platforms were used for sequencing. De novo assembly was performed using Unicycler v0.5.0 (Wick et al., 2017). Centrifuge v1.0.4 was used to classify contigs and remove contaminating sequences (Kim et al., 2016).

Gene prediction and the annotation of assembled sequences were performed via Prokka v1.12 (Seemann, 2014). STs were determined via MLST v2.16.1 and reconfirmed using stringMLST v0.6.3 based on reads. For any strain that did not match an existing ST number, we uploaded the sequencing data to the Klebsiella MLST database (https://bigsdb.pasteur.fr/cgibin/bigsdb/bigsdb.pl?db=pubmlst klebsiellas eqdef&page=login) to obtain a new ST number. Capsular K serotypes were identified by whole-genome data using Kleborate v2.2.0 (Lam et al., 2021).

The collection of Prokka-annotated sequences was analyzed with Roary v3.13.0 (Page et al., 2015) to infer core and pangenomes. Maximum-likelihood phylogenetic tree analyses were generated using RAxML-NG v1.0.1 (Kozlov et al., 2019) with 1000 bootstrap replicates and the specific model selected by ModelTest-NG v0.1.7 (Darriba et al., 2020). The tree was visualized using the Figtree v1.4.4 graphical viewer (http://tree.bio.ed.ac.uk/software/figtree/) and iTOL (https://itol.embl.de.).

The presence of resistance and virulence genes was predicted using the Comprehensive Antibiotic Resistance Database (Jia et al., 2017) and Virulence Factor Database (Chen et al., 2016), respectively, in the ABRicate v0.8.13 program. A heatmap was generated using Tbtools software (Chen et al., 2020). Platon v1.6 (Schwengers et al., 2020) was used to classify plasmid contigs and confirmed plasmids via plasmidID v1.6.5(https://github.com/BU-ISCIII/plasmidID.wiki. git).

Pangenome analysis was performed using Roary v3.13.0 (Page et al., 2015). Rarefaction and accumulation curves were created using R version 4.2.1. Functional annotation of genes was done on RAST using the SEED subsystems approach. For RAST annotation, nucleotide files were uploaded to RAST by default features (RAST annotation scheme: RASTtk, automatically fix errors, fix frameshifts, build metabolic model, backfill gaps, turn on debug: yes, verbose level: 0).

All data were analyzed using R version 4.2.1. The chi-square test or Fisher’s exact test were used to analyze categorical variables. Continuous variables were presented as means ± standard deviation (SD) and were compared using Student’s t-test. Normally distributed continuous variables were expressed as means ± SD and compared using Student’s t-test. Non-normally distributed continuous variables were presented as medians with interquartile ranges (IQR) and compared using the Mann–Whitney U test. Univariate and multivariable logistic regressions were performed to explore the risk factors of HmKp infection. All variables with a value of p<0.1 within univariate analysis were included in the following multiple logistic regression model. Stepwise regression was used to identify statistically significant predictors. p<0.05 was considered as statistically significant.

A total of 203 K. pneumoniae isolates were collected; 90 (44.3%) exhibited a positive string test and were designated as HmKp. The proportions of HmKp isolates among invasive K. pneumoniae isolates in 2018, 2019, 2020, and 2021 were 35.3%(12/34), 56.0% (28/50), 40.2%(27/67), and 44.2%(23/52), respectively. The χ2 test for trend indicated that the annual proportions of HmKp isolates were similar over the four-year timespan ( Figure 1A ). The four departments with the highest prevalence of HmKp were Adolescent Hepatology (5/5,100.0%), Emergency (3/3,100.0%), Uninfected Hepatology (3/4, 75.0%), and Transplantation Medicine (6/8,75.0%) ( Figure 1B ). HmKp was isolated most often from the bloodstream (43/90, 48%), followed by ascites (19/90, 21%) and liver abscesses (15/90,17%). Non-HmKp isolates were obtained most often from the bloodstream (42/113,37%), ascites (31/113,27%), and bile (25/113, 22%) ( Figure 1C ).

Demographic and clinical features of the 203 infected patients are listed in Table 1 . The prevalence of diabetes mellitus (58.9% vs 36.3%; p=0.001) and liver abscess (26.7% vs 8.8%; p=0.001) were significantly higher in the HmKp than in the non-HmKp group. Solid malignancy (22.2% vs 41.5%; p=0.004) and digestive diseases (10.0% vs 27.4%; p=0.002) were more prevalent in the non-HmKp group. Additionally, there was a lower proportion of patients with HmKp infections undergoing drainage (32.2% vs 47.8%; p=0.025); p=0.035) compared to patients with non-HmKp infections. Fever was more common in the HmKp group than in non-HmKp group (74.4% vs 60.2%; p=0.032). Fewer patients were treated with penicillin in the HmKp group than in the non-HmKp group (7.8% vs 19.5%; p=0.018).

Subsequent to the above analysis, all factors from univariate analysis with a p<0.1 were included in the multivariate model. In the multivariate analysis, diabetes mellitus (odds ratio [OR]=2.20, 95% confidence interval [CI]=1.20-4.05; p=0.010) and liver abscess (OR=2.93, 95%CI=1.29-7.03; p=0.012) were associated with increased odds of HmKp infection. Solid malignancy (OR=0.45, 95%CI=0.23-0.85; p=0.016) and digestive diseases (OR=0.38,95%CI=0.15-0.86; p=0.025) were associated with markedly decreased odds of HmKp infection ( Table 2 ).

K. pneumoniae isolates were highly diverse, comprising 87 STs and 54 serotypes. Eight of the STs (ST6105, ST6106, ST6107, ST6108, ST6109, ST6110, ST6111, and ST6112) were first identified in this study. The most frequent STs in the HmKp group were ST23(25/90, 27.8%) and ST65 (15/90, 16.7%). In the non-HmKp group, the most frequent ST was ST11 (12/113, 10.6%) ( Figure S1 ). The most frequent serotypes in the HmKp group were KL2 (31/90, 34.4%) and KL1(27/90, 30.0%), which accounted for >50% of isolates ( Figure 2 ). Interestingly, most KL1 strains belonged to ST23 (31/35, 88.6%), and 25 of KL1-ST23 strains were HmKp.

A total of 126 virulence genes were recorded, with KL1 isolates having the highest repertoire (120/126, 95.2%). Thirteen virulence genes located on the capsular polysaccharide synthesis (cps) region were identified only in the KL1 serotype. The prevalence of 79 genes was statistically higher in HmKp; these included genes encoding colibactin; allantoin utilization; aerobactin, yersiniabactin, salmochelin, and lipopolysaccharide (LPS) synthesis; the type VI secretion system (tssD, tssF, tssJ, tssK, tssL, tssM, fyuA); and the gene cluster for capsule production (rmpA, rmpA2, ugd, gnd, manB, manC, gmd, wcaG, wcaH, wcaI, wcaJ, wcsS, wcsT, wclY, wza, wzb, wzc, wzi, wzx, magA) (p<0.05) ( Figure 2 ). In addition, 112 HvKp isolates were defined by various combinations of peg-344, iroB, iucA, _prmpA, or _prmpA2, and included 88 HmKp. In the HmKp group, only two isolates were not HvKp.

HmKp isolates were generally more susceptible than non-HmKp isolates to antimicrobial agents that included meropenem, imipenem, aztreonam, cefepime, ceftazidime, ciprofloxacin, levofloxacin, piperacillin, tobramycin, doxycycline, minocycline, and ticarcillin (p<0.05) ( Table S2 ). Only one of 14 carbapenem-resistant isolates was an HmKp (ST11-KL64). Notably, the prevalence rates of SHV-type ESBL genes (SHV-11, SHV-190, SHV-207) were higher in the HmKp group than in the non-HmKp group (p<0.05) ( Figure S2 ).

The core genome-based phylogenetic tree is shown in Figure 3 . Klebsiella variicola subsp. Variicola, kv291 collected from the same hospital during the study period was used as an outgroup to root the tree. A cluster was defined as a large branch with the same serotype and sequence type. Most isolates had clearly distinct phenotypes in the cluster, such as HMV and carbapenem resistance. According to these criteria, we identified 9 clusters. Cluster 5 was the largest, and was comprised entirely of KL1 strains. KL2 isolates clustered in different sub-branches with different STs (Cluster 2: KL2 & ST86; Cluster 3: KL2 & ST380; Cluster 7: KL2). HmKp strains were scattered in different clusters whereas carbapenem-resistant (CPKP) strains were grouped closely in Cluster 9 (ST11).

To evaluate the phylogenetic relationships of the 203 K. pneumoniae isolates in public genomes, we selected 116 K. pneumoniae strains from PATRIC to construct a phylogenetic tree together with the 203 isolates in our study using the method described above. We downloaded complete gene sequences of 116 K. pneumoniae isolates causing invasive infections worldwide from 2018-2021. Our 203 isolates did not cluster together but were scattered within the clusters of public genomes. HmKp isolates were dispersed throughout the phylogenetic tree, consistent with our phylogenetic analysis ( Figure S3 ).

The pangenome of 203 K. pneumoniae isolates consisted of 22435 genes, of which 3652 (16.3%) belonged to the core genome, and 18783 (83.7%) were accessory ( Figure 4 ). Rarefaction and accumulation curves were indicative of an open pangenome with a well-defined core, and suggested that additional genomes were unlikely to impact the size of the core genome much further. We conducted additional analysis to distinguish genes between the HmKp and non-HmKp groups. We applied the accuracy rate as an indicator, and found that 17 of 22435 genes were highly associated with HmKp ( Table 3 ). When a particular gene was used as a biomarker for HmKp, accuracy rate was calculated by dividing the sum of the number of isolates carrying the gene in the HmKp group and the number of isolates without the gene in non-HmKp group divided by the total number of isolates (Ye et al., 2016).

The 17 predicted HmKp-associated genes were located on a 21889 bp DNA fragment carried by the pK2044 plasmid and also detected on a 191041 bp pM186-like plasmid from kpn191. Among the 17 genes, rmpA is a known HmKp associated gene, and 7 genes (iroB, iroC, iroD, iroN, fecA, fecI, fecR) are related to iron-acquisition. The liver abscess-causing K. pneumoniae (LAKP)-associated gene pagO was also encoded in this fragment.

Classification of genes into discrete functional units provides valuable insight into how resources are allocated to each function. We functionally annotated the K. pneumoniae genome using RAST, resulting in the classification of genes into 25 subsystems. To assess HmKp related functions, median numbers of genes of each group were compared at the RAST subsystem level. As shown in Figure 5 , sixteen categories showed significant difference between HmKp and non-HmKp isolates.

Most notably, the number of genes in the “Iron acquisition and metabolism” category was significantly higher in HmKp than in non-HmKp isolates. The numbers of genes in the “Virulence, Disease and Defense”, “Regulation and Cell signaling,” “Stress Response,” “Cofactors, Vitamins, Prosthetic Groups, Pigments,” “Fatty Acids, Lipids, and Isoprenoids,” “Nitrogen Metabolism,” “ Amino Acids and Derivatives,” and “Phosphorus Metabolism” categories were significantly greater in HmKp compared to non-HmKp isolates. In addition, non-HmKp isolates had significantly more members of genes in 6 other categories (“Cell Wall and Capsule,” “Phages, Prophages, Transposable elements, Plasmids,” “RNA Metabolism,” “Respiration,” “Miscellaneous” “Metabolism of Aromatic Compounds,” and “Carbohydrates”) than HmKp strains ( Figure 5 ).