Whole-Genome Sequencing and Comparative Genomics Analysis of a Newly Emerged Multidrug-Resistant Klebsiella pneumoniae Isolate of ST967

ABSTRACT Klebsiella pneumoniae is a common cause of hospital- and community-acquired infections globally, yet its population structure remains unknown for many regions, particularly in low- and middle-income countries (LMICs). Here, we report for the first-time whole-genome sequencing (WGS) of a multidrug-resistant K. pneumoniae isolate, ARM01, recovered from a patient in Armenia. Antibiotic susceptibility testing revealed that ARM01 was resistant to ampicillin, amoxicillin-clavulanic acid, ceftazidime, cefepime, norfloxacin, levofloxacin, and chloramphenicol. Genome sequencing analysis revealed that ARM01 belonged to sequence type 967 (ST967), capsule type K18, and antigen type O1. ARM01 carried 13 antimicrobial resistance (AMR) genes, including blaSHV-27, dfrA12, tet(A), sul1, sul2, catII.2, mphA, qnrS1, aadA2, aph3-Ia, strA, and strB and the extended-spectrum β-lactamase (ESBL) gene blaCTX-M-15, but only one known virulence factor gene, yagZ/ecpA, and one plasmid replicon, IncFIB(K)(pCAV1099-114), were detected. The plasmid profile, AMR genes, virulence factors, accessory gene profile, and evolutionary analyses of ARM01 showed high similarity to isolates recovered from Qatar (SRR11267909 and SRR11267906). The date of the most recent common ancestor (MRCA) of ARM01 was estimated to be around 2017 (95% confidence interval [CI], 2017 to 2018). Although in this study, we report the comparative genomics analysis of only one isolate, it emphasizes the importance of genomic surveillance for emerging pathogens, urging the need for implementation of more effective infection prevention and control practices. IMPORTANCE Whole-genome sequencing and population genetics analysis of K. pneumoniae are scarce from LMICs, and none has been reported for Armenia. Multilevel comparative analysis revealed that ARM01 (an isolate belonging to a newly emerged K. pneumoniae ST967 lineage) was genetically similar to two isolates recovered from Qatar. ARM01 was resistant to a wide range of antibiotics, reflecting the unregulated usage of antibiotics (in most LMICs, antibiotic use is typically unregulated.) Understanding the genetic makeup of these newly emerging lineages will aid in optimizing antibiotic use for patient treatment and contribute to the worldwide efforts of pathogen and AMR surveillance and implementation of more effective infection prevention and control strategies.

C ommonly associated with the environment (1), plants (2), animals (3,4), and humans (5), Klebsiella pneumoniae has emerged as a major global public health threat due to becoming a leading cause of nosocomial and community infections over the past decade. This opportunistic pathogen, responsible for a wide range of infections (6)(7)(8)(9)(10), has also been associated with multidrug resistance (MDR). A recent survey conducted in 41 countries by the Central Asian and the European Surveillance of Antimicrobial Resistance network (11) revealed that .50% of the isolates recovered from patients (predominantly in Southern and Eastern Europe) were resistant to thirdgeneration cephalosporins and .25% were resistant to carbapenems. Moreover, mortality rates caused by MDR K. pneumoniae have been reported to be as high as 50% or higher in clinical settings (12,13). This is further exacerbated by new emerging clones, such as sequence types 101 (ST101), ST307, and ST147, which have acquired virulence factors and antimicrobial resistance (AMR) genes, contributing to their success as persistent nosocomial infections (14,15). One of the key antibiotic resistance gene families found in extended-spectrum b-lactamase (ESBL)-producing K. pneumoniae isolates that have been dominating nosocomial infections globally is the bla CTX-M-type (16,17). In addition, bla CTX-M-15 has been identified as the most predominant type and has rapidly disseminated in K. pneumoniae isolates recovered from humans (18,19) and animals (3,20).
The driving force for the development of multidrug resistance is selective pressure, which is due to continuous exposure to multiple antibiotics in settings where they are overused and/or misused, such as health care and agriculture (21). In many low-and middle-income countries (LMICs), public exposure to MDR bacteria is increasingly high due to uncontrolled use of "antibiotic growth promoters" in animals and the ability to buy antibiotics "over-the-counter" (22). In addition, there are few clinical guidelines for the use of antibiotics for human treatment in animals. Moreover, severe financial constraints, a lack of diagnostic facilities, and inadequate pathogen surveillance systems make infection prevention and control (IPC) extremely challenging, which encourages the development of multidrug-resistant strains (23).
The spread of MDR K. pneumoniae is a global public health issue. Understanding the genetic makeup of new emerging lineages of such pathogens will enable us to optimize antibiotic use for patient treatment and contribute to the worldwide efforts of tackling antibiotic resistance and design appropriate IPC measures (24,25). Wholegenome sequencing (WGS) has been widely used in studying the global transmission and molecular mechanisms of K. pneumoniae (5,26,27). However, reports on wholegenome sequencing and population genetics of K. pneumoniae from LMICs are scarce (28). A systematic review we conducted for Armenia, a low-and middle-income country, showed that only one study has analyzed secondary data on AMR for the period of 2016 to 2019 (29). We recently conducted a small-scale pilot study to identify the genetic diversity of methicillin-resistant Staphylococcus aureus (MRSA) in a teaching hospital in Armenia, which provided insights into a previously unrecognized diversity of MRSA clones in Armenia, including pandemic and sporadic lineages seen internationally (30). Further WGS studies identified a novel CC8 MRSA clone circulating in a hospital setting in Yerevan, the capital of Armenia. Cross transmission/contamination due to the failure of infection control and prevention strategies was evident as MRSA isolates belonging to the same lineages were recovered from patients and hospital environments (31).
Here, we report for the first time the whole-genome sequencing analysis of a K. pneumoniae isolate, ARM01, recovered from a patient in Armenia and provide insights into its genomic features. Further comparative genomics analysis revealed its temporal and spatial phylogenetic dynamics. hospitals in Armenia between January 2019 and August 2019 (see Table S1 in the supplemental material). Whole-genome sequencing revealed that four out of eight isolates belonged to sequence type 307 (ST307), with other four isolates each belonging to ST37, ST147, ST807, and ST967. All isolates were resistant to ampicillin, amoxicillinclavulanic acid, ceftazidime, and cefepime, but were sensitive to meropenem (Table 1). In addition to the above antibiotics, the ARM01 isolate belonging to ST967 described further in this study was also resistant to norfloxacin, levofloxacin, and chloramphenicol (Table 1).
The core SNP phylogenetic tree showed that ARM01 was phylogenetically closely related to two human isolates recovered from Qatar (SRR11267906 and SRR11267909), forming a single monophyletic cluster ( Fig. 1). Moreover, assessing the 3,474 accessory genes identified in the pangenome, we found that ARM01 shared many of the same accessory genes found in two isolates recovered from Qatar (r = 0.77 for SRR11267909 and r = 0.72 for SRR11267906) (see Fig. S1A in the supplemental material). The close genetic similarities of ARM01 and the two Qatar isolates (SRR11267906 and SRR11267909) were further confirmed based on genome-wide average nucleotide identity (ANI) analysis. The Hadamard values of ARM01 with SRR11267909 and SRR12267906 were 0.988 and 0.984, respectively ( Fig. S1B and Table S3). In addition, hierarchical clustering of plasmid replicon patterns revealed that ARM01 shared the same replicon, IncFIB(K)(pCAV1099-114), with the Qatar isolates ( Fig. S2 and Table S4). The above findings indicate that the Armenian isolate (ARM01) was genetically closely related to two Qatar isolates (SRR11267906 and SRR11267909).
Evolutionary origins of K. pneumoniae ARM01. As to date no whole-genome sequencing studies have been conducted to characterize K. pneumoniae recovered from patients or hospital settings in Armenia, we wanted to estimate the date of the most recent common ancestor (MRCA) of ARM01 and infer its evolutionary origins with those phylogenetically closely related K. pneumoniae isolates. Therefore, we reconstructed a maximum clade credibility (MCC) tree using BEAST (Fig. 2). One isolate obtained from the Pathogenwatch database was removed from the analysis due to the absence of the   (Table S5). Finally, a spatial phylogeographic reconstruction was performed (Fig. S3), confirming that the Armenian isolate was genetically related to the two isolates recovered from Qatar. The Bayes factor for the transmission from Qatar to Armenia was 0.578, with a posterior probability of 0.122, which may be due to the small sample size. Comparative analysis of antimicrobial resistance and virulence genes. The AMR gene profile of ARM01 was compared to those of the 22 K. pneumoniae ST967 isolates obtained from the Pathogenwatch database and previously published studies (32). In total, 48 AMR genes were identified in the resistome of these isolates. A total of 22 out of 23 K. pneumoniae ST967 isolates possessed genes conferring resistance to three or more classes of antibiotic ( Fig. 1; Table S6). Only one of the isolates (ERR5530481), which was recovered from a broiler, harbored one AMR gene (bla SHV-27 gene). A number of AMR genes, including strA (15/23), strB (15/23), qnrS1 (15/23), mphA (12/23), sul1 (15/23), sul2 (17/23), tet(A) (16/23), bla TEM-1D (13/23), and bla SHV-27 (19/23), were detected in most of the K. pneumoniae ST967 isolates, including ARM01, suggesting that multidrug resistance was preserved in the ST967 lineage recovered from a human host. Moreover, we detected 15 genes associated with b-lactam antibiotics in all isolates, including 6 ESBL-associated genes conferring resistance to third-generation cephalosporins (bla CTX-M-3 , bla CTX-M-14 , bla CTX-M-15 , bla SHV-12 , bla SHV-12 , and bla SHV-27 1178H) and 3 genes conferring resistance to carbapenems (bla NDM-1 , bla KPC-2 , and bla CTX-M-33 ). Among the identified ESBL-associated genes, bla CTX-M- 15 and bla CTX-M-3 were the most common. bla SHV27 was present in most of the K. pneumoniae ST967 isolates (19/23 [82.6%]). On average, each K. pneumoniae ST967 isolate possessed 12 AMR genes resistant to various classes of antibiotics, which was consistent with the number of AMR genes we found in ARM01 (n = 13) (Fig. 1). We also found that ARM01 shared many of the same resistance genes with two Qatar isolates (SRR11267906 and SRR11267909) for b-lactamases (bla SHV-27 ), ESBL (bla CTX-M-15 ), trimethoprim (dfrA12), sulfonamides (sul1 and sul2), phenicol (catII.2), macrolides (mphA), fluoroquinolones (qnrS1), and aminoglycosides [aadA2 and aph3 (3)-IIa]. ARM01 also shared aminoglycoside resistance genes (strA and strB) with one of the Qatar isolates, SSR11267909 ( Fig. 1; Fig. S4).
We further compared virulence genes identified in K. pneumoniae ARM01 with those detected in the other 22 K. pneumoniae ST967 isolates (Table S7). Overall, 3 virulence genes were identified. One known virulence factor, yagZ/ecpA, was detected in all K. pneumoniae ST967 isolates. ERR7672084 (recovered from Switzerland) and ERR4782243 (recovered from Philippines) possessed additional aerobactin-associated genes: iucC and iucB.

DISCUSSION
Whole-genome sequencing (WGS) surveillance of K. pneumoniae is imperative due to its increased resistance to antibiotics and high mortality rate (12). Its rapid dissemination as a nosocomial and community-associated pathogen globally has led to public health crises (1, 10). As a tool, WGS can help advance our knowledge and understanding of pathogen global dissemination issues (33). However, WGS surveillance of K. pneumoniae is still scarce in LMICs, leaving essential knowledge gaps in genomic and epidemiological monitoring of this species. To date, only a few studies of K. pneumoniae ST967 have been reported (34,35), with none from Armenia. To the best of our knowledge, this is the first study reporting on whole-genome sequencing of K. pneumoniae ST967 recovered from a patient in Armenia.
In this study, we report the whole-genome sequencing and comparative analysis of a multidrug-resistant (resistant to 7 out of 11 antibiotics tested) K. pneumoniae ST967 isolate, designated ARM01. WGS analysis revealed that ARM01 carried 13 AMR genes conferring resistance to b-lactamase (bla CTX-M-15 and bla SHV-27 ), trimethoprim (dfrA12), tetracycline [tet(A)], sulfonamides (sul1 and sul2), phenicols (catII.2), macrolides (mphA), fluoroquinolones (qnrS1), and aminoglycosides (aadA2, aph3-Ia, strA, and strB) ( Fig. 1; see Table S6 in the supplemental material). Notably, there was a strong correlation between the results of phenotypic antibiotic susceptibility and the antibiotic resistance predicted by the AMR genes. Multidrug resistance was preserved throughout the K. pneumoniae ST967 lineage, with 22 out of 23 isolates harboring multiple AMR genes predicted to be resistant to three or more classes of antibiotics.
Hospital outbreaks of MDR K. pneumoniae, particularly in neonatal units, have often been associated with CTX-M types of ESBL producers (10,36). Due to the risk of easy dissemination of AMR genes to other bacterial species through mobile genetic elements (MGEs), such as transposons and integrons, ESBLs have become a serious issue (37). In this study, 73.9% (17/23) of K. pneumoniae ST967 isolates were positive for ESBL production, with bla CTX-M-3 (7/23) and bla CTX-M-15 (7/23) being the most common ESBL genes detected. A recent global surveillance program that was carried out in 48 countries from 2015 to 2019 showed that the Middle East/Africa had the highest prevalence of ESBL non-carbapenem-resistant Enterobacteriaceae (non-CRE) K. pneumoniae isolates, ranging from 32.1% to 39.0%, and that, in Eastern Europe, the rates ranged from 22.0% to 33.1%, with Lithuania having the highest prevalence rate at 55.2% (38).
As a CTX-M-1 group member, bla CTX-M-3 was originally disseminated in Poland and later (in the 2000s) in Europe (39), and it has only been reported on a few occasions (40,41). Interestingly, all bla CTX-M-3 isolates in this study were collected in Asia (Japan, China, and Philippines). Both belong to the CTX-M-1 group: bla CTX-M-15 differed from bla CTX-M-3 by a single amino acid change (an asparagine-to-glycine substitution in position ABL238) (42). bla CTX-M-15 has widely spread around the whole world and is known to hydrolyze many of the cephalosporins. bla CTX-M-15 was reported to be associated with epidemic MDR K. pneumoniae isolates (ST15, ST147, and ST336) (43). In addition, bla CTX-M-15 is prevalent in 93.7% (89/95 isolates) of K. pneumoniae ST307 isolates globally (44). bla SHV-27 was shown to be frequently associated with CTX-M-type ESBLs (45). bla SHV-27 was commonly found in K. pneumoniae, Escherichia coli, and Enterobacter cloacae isolates recovered from human clinical samples and companion animals (46)(47)(48). In this study, bla SHV-27 was detected both in ARM01 as well as in K. pneumoniae isolates belonging to the ST967 lineage (19/23 isolates). Moreover, ARM01 displayed nearly identical AMR and plasmid profiles to one of the Qatar isolates (SRR11267909), suggesting their close genetic relatedness.
In contrast to the hypervirulent K. pneumoniae isolates, ARM01 did not carry genes for any of the classical virulence factors, such as yersiniabactin (ybt), colibactin (clb), salmochelin (iro), aerobactin (iuc), or hypermucoidy (rmpA and rmpA2) (49,50). These virulence-associated genes have been associated with invasive infections and were frequently found in hypervirulent K. pneumoniae isolates (51). The absence of some typical virulence determinants (such as rmpA, rmpA2, iuc, and iro) has previously been shown to result in nonhypervirulent phenotypes (44,52), indicating the relatively low pathogenicity of ARM01. None of the K. pneumoniae ST967 isolates, with the exception of ERR4782243 and ERR7672084, contained the canonical hypervirulence markers; however, all K. pneumoniae ST967 isolates acquired the yagZ/ecpA gene. As a gene coding for an E. coli common pilus adherence factor, yagZ/ecpA has been reported to be involved in cell adhesion and biofilm formation (53). It was commonly found in K. pneumoniae isolates belonging to various sequence types previously recovered from patients (54). Similar to this, another study showed that 96% (66/69) of the hospitalacquired K. pneumoniae isolates harbored yagZ/ecpA and that there was a 94% phenotypic correlation of ECP production, which was essential for the formation of the adhesive structure (55).
Time-calibrated phylogenetic analysis estimated the MRCA of K. pneumoniae ST967 lineage to be 2004, implying that the K. pneumoniae ST967 subgroup was a newly emerged clone that had spread across countries over the last 2 decades. Comparative analysis revealed that ARM01 was genetically similar to two Qatar isolates (SRR11267906 and SRR11267909). SRR11267906 and SRR11267909 represent case sequencing data recovered from the same hospital (35) with 64 SNP differences, whereas the numbers of SNP differences between ARM01 and SRR11267909 and SRR11267906 were 27 and 45, respectively. It has previously been reported that two bacterial isolates may be considered closely related clones if they possess a relatively low SNP distance (fewer than 50 SNPs) due to the short divergence time (56,57). S. David et al. (58) showed that the most frequent minimum SNP distance for carbapenem-resistant K. pneumoniae ST258/512 (found in 32 different European countries) was 45 SNPs and was considered an international spread. Using spatial phylogenetic reconstruction of K. pneumoniae ST967, we reconstructed its transmission dynamics between countries, which indeed shed light on the high genetic relatedness of ARM01 and the Qatar isolates (SRR11267906 and SRR11267909). Undoubtably, there is an urgent need for additional data collection and analysis of this newly emerged K. pneumoniae lineage in Armenia, but also further molecular epidemiology investigations are warranted to reveal its evolutionary origins.
Conclusions. This study revealed the genomic characteristics of ARM01, an isolate belonging to the newly emerged K. pneumoniae ST967 lineage. Multilevel analysis found it to be genetically highly similar to two Qatar isolates, which shared and had descended from a common ancestor in 2017. ARM01 was resistant to a wide range of antibiotics despite harboring just one known virulence-associated gene. Additional surveillance studies are warranted to further understand the prevalence and molecular epidemiology of K. pneumoniae and other Gram-negative pathogens circulating in health care and community settings in Armenia and to identify possible risk factors, antimicrobial resistance associated with these isolates, genetic diversity, and disease burden in this country to inform National Infection Prevention and Control Policy.
Although our data demonstrate compelling evidence of characteristics of ARM01 and phylogeny of the K. pneumoniae ST967 lineage, there are several limitations to our study. The main limitation of this study is the small sample size both received for this study and those publicly available data. The fact that only one out of the eight K. pneumoniae isolates received for this study was ST967 is also a major limitation as we could not fully decipher the population structure of K. pneumoniae ST967 in Armenia. Furthermore, only 22 K. pneumoniae ST967 data were publicly available, and therefore, comparative genomics analysis of this recently emerged subtype remains to be explored. In addition, epidemiological data were lacking; we do not know if there were any associated risk factors, i.e., travel abroad.

MATERIALS AND METHODS
Identification and antibiotic susceptibility testing. Eight K. pneumoniae isolates were received from the Medical Microbiology laboratories of three hospitals in Armenia between January 2019 and August 2019. All isolates were recovered from various clinical specimens (urine, sputum, throat, blood, and stool) of hospitalized patients. The isolates were identified using a matrix-assisted laser desorption ionization-time of flight mass spectrometry (MALDI-TOF-MS) as described previously (59).
Whole-genome sequencing and assembly. Genomic DNA was extracted using the TIANamp bacterial DNA kit (Tiangen, China), and paired-end sequencing libraries were constructed using Nextera XT DNA sample preparation kits or the TruSeq DNA HT sample prep kit (Illumina, USA) following the manufacturer's instructions. Whole-genome sequencing was performed using the Illumina HiSeq platform, with a minimum coverage of 100-fold and generating 151-bp sequence reads. The quality of raw reads was assessed using FASTQC v.0.11.9 (https://github.com/s-andrews/FastQC) and trimmed using Trimmomatic v.0.38 to eliminate adaptors and low-quality sequences (61). Trimmed reads were de novo assembled using Velvet v. 1.2.10 (62), with contigs less than 500 bp being removed. Prokka v.1.14.5 was used for draft genome annotation (63).
Temporal and spatial phylogenetic reconstruction. To reconstruct the time-calibrated phylogeny of all isolates, we conducted a Bayesian analysis using BEAST v.1.8.4 package (76). jModelTest v.2.1.10 (77) was used to determine the best-fit model and parameters of nucleotide substitution. A coalescent Bayesian skyline model for population growth was performed, with a relaxed clock and a general time-Newly Emerged Multidrug-Resistant K. pneumoniae ST967 Microbiology Spectrum reversible (GTR) nucleotide substitution model. BEAST was run across four independent iterations, each set to a Markov chain Monte Carlo (MCMC) of 10 million generations with samples taken at every 1,000 generations. LogCombiner was used to combine output files from four independent runs. A combined log file was analyzed with Tracer v.1.7.2. The merged tree file was used to build the maximum clade credibility (MCC) tree, which was then displayed by FigTree v.1.4.4. The Bayesian stochastic search variable selection (BSSVS) method was used to establish a Bayesian factor (BF) test to infer the possibility of spatial geographic transmission based on discrete geographic information. The MCC tree with phylogeographic reconstruction was generated with 10% burn-in using the Tree Annotator. SpreaD v1.0.7 (78) was used to visualize the location-annotated MCC tree.
Pan-genome analysis. The annotated genome assemblies were used to construct the pan-genome analysis using Roary v.3.13.0 (79). Accessory genes were classified as genes found in less than 99% of all isolates. To compare the accessory gene profiles among all isolates, a hierarchically clustered heat map was constructed using the Pearson product-moment correlation coefficient in R v.3.6.2 (R Core Team; https://www.R-project.org) according to the presence (1) and absence (0) of all accessory genes.
Data availability. The short-read sequencing data generated in this study were deposited in the ENA database under the accession no. ERR9882335. Individual accession numbers for the genome sequencing data are included in Table S2.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only. SUPPLEMENTAL FILE 1, PDF file, 0.6 MB. SUPPLEMENTAL FILE 2, XLSX file, 0.04 MB.

ACKNOWLEDGMENTS
This work was supported by the Royal Society QR GCRF funding. H.V.M. contributed to conceptualization and design of the study, data analysis, writing, and critical review of the manuscript. J.S. contributed to WGS data analysis (bioinformatics), writing, and critical review of the manuscript. R.C. contributed to reviewing data analysis and the manuscript. M.M.T.-S. contributed to sample collection, laboratory work, and reviewing the manuscript. N.K. contributed to sample collection and laboratory work. J.C. contributed to whole-genome sequencing. L.Z. contributed to spatial phylogenetic analysis. T.J. contributed to review of the manuscript. All authors read and approved the final manuscript.
We declare no conflict of interest.