From 83 KPC isolates carrying the bla_KPC gene collected from nosocomial infections in Minas Gerais state (Southeast Brazil) in 2011 we selected the isolate named KPC-142, which displayed the most divergent DNA pattern by RAPD cluster analysis (Figure S1). This isolate was initially assigned as K. pneumoniae KPC-142 and susceptibility tests were performed according to CLSI guidelines 2017, using Vitek 2 system and AST-N239 cards (bioMérieux, Inc., Durham, NC) according to the manufacturer's instructions and using Gram-negative identification (GN ID) cards (Bobenchik et al., 2017).

The bacterial isolate was grown overnight on nutrient agar and the DNA was extracted using the method originally described by Coimbra et al. (1999). The isolate KPC-142 is maintained in the certified strains collection of Ezequiel Dias Foundation (FUNED), Belo Horizonte, Brazil.

For the comparative genomic analyses and determination of average nucleotide identity (ANI), we used a bacterial strain panel representative of each Klebsiella phylogroup (Holt et al., 2015) that included: (i) nine K. quasipneumoniae subsp. similipneumoniae (phylogroup II-B) (HKUOPA4, HKUOPJ4, HKUOPL4, KPC-142, ATCC 700603, 07A044 T, MGH 44, 193_KOXY, and 385_ECLO), (ii) nine K. quasipneumoniae subsp. quasipneumoniae (phylogroup II-A) (UCICRE14, 01A030 T, FI HV 2014, MGH96, 18A069, ARLG-2711, PO1285, AK_SD_007, and 21_GR_13), nine K. pneumoniae (phylogroup I) (DSM 30104 T, Kp13, 1084, CAV1193, HS11286, KCTC 2242, KP617, MGH 78578, and NTUH-K2044), and nine K. variicola (phylogroup III) (BZ19, GJ1, UCICRE10, MGH 40, MGH 20, UCI 18, BIDMC 61, BIDMC88, and MGH 80), totaling 36 bacterial genomes. An additional 31 isolates were included for phylogenetic studies (see section MLST classification and phylogeny).

Genomic DNA was sequenced on an Illumina HiSeq 2500 sequencer using Nextera XT paired-end run with a 500-bp insert library at the High-Throughput Sequencing Platform of the Oswaldo Cruz Foundation (Fiocruz, Rio de Janeiro, Brazil). The assembly of reads into scaffolds, based on 3,804,017 reads, was accomplished using a combination of Newbler v 2.6 (Roche Inc.) and SPAdes 3.10.0 (Bankevich et al., 2012) programs. 99.99% of the bases on the assembled genome had an average Phred quality of >40, calculated by Newbler. Gaps intra- and inter-scaffolds were resolved using a three-phase approach: (1) by extracting consensus sequences that did not form gaps in one of the assemblies; (2) by performing local assemblies of reads falling in gap termini; and (3) using the GapFiller tool (Boetzer and Pirovano, 2012). Further details of each step are presented in the Figure S2.

Genome annotation was performed using the System for Automated Bacterial Integrated Annotation (SABIA) (Almeida et al., 2004). The steps performed by this platform consist of an initial gene prediction process using Glimmer (Delcher et al., 1999) and GeneMark (Besemer et al., 2001) including a start codon correction routine based in the multiple alignment of similar proteins identified by BLASTP against NCBI proteins database (www.ncbi.nlm.nih.gov/protein). CDSs identified after this process were annotated applying an automated annotation pipeline, where each open reading frame (ORF) is submitted to similarity searches using both nucleotide and amino acid sequences by Basic Local Alignment Search Tool (BLAST) (Altschul et al., 1990) against KEGG (www.genome.jp/kegg), NCBI-nr (www.ncbi.nlm.nih.gov/protein) and UniProtKB/Swiss-Prot (www.uniprot.org) databases, in addition to the prediction of protein domains and important catalytic sites using InterPro database (Finn et al., 2017), and the results are made available on the screen for the assessment of expert users. In this process, only CDSs with predicted lengths over 50 aa and at least one BLASTP hit against any of the four databases mentioned above were considered. The transfer RNAs (tRNAs) were detected by tRNAscan-SE (Lowe and Eddy, 1997) and the annotation of ribosomal RNA genes was carried out by RNAmmer (Lagesen et al., 2007).

Further bioinformatics analyses included BLAST searches against AtlasT4SS database (Souza et al., 2012)/Resfinder database (Zankari et al., 2012)/ARDB-Antibiotic Resistance Genes Database (Liu and Pop, 2009)/Comprehensive Antibiotic Resistance Database (Jia et al., 2017) and Virulence Factors database (VFDB) (Chen et al., 2005), which were performed for identification of type 4 secretion systems (T4SS), acquired antimicrobial resistance genes and bacterial virulence factors, respectively. Typing of plasmids was carried out by in silico detection using PlasmidFinder database (Carattoli et al., 2014).

Average nucleotide identity (ANI) was calculated using BLASTn (ANIb) in JSpeciesWS (Richter et al., 2016). Reference type strains of K. quasipneumoniae belonging to subspecies K. quasipneumoniae subsp. quasipneumoniae 01A030 (GenBank accession no. GCA_000751755) and K. quasipneumoniae subsp. similipneumoniae 07A044 (GenBank accession no. GCA_000613225) and type strain K. pneumoniae subsp. pneumoniae DSM 30104 (GenBank accession no. GCA_000281755), in addition to strain K. variicola GJ1 (GenBank accession no. GCA_001989495.1), were used to perform all-versus-all comparisons with a bacterial strain panel including KPC-142. An ANI threshold of ≥96% or greater was considered to delineate species boundaries as it correlates well to DNA-DNA hybridization studies (Goris et al., 2007; Federhen et al., 2016).

Comparative genomic and gene cluster analyses were performed using the Gview tools (Petkau et al., 2010) and Proteinortho (Lechner et al., 2011), respectively. The comparisons involved the bacterial strain panel. In Proteinortho, co-orthologous clustering was calculated using BLASTp, with an identity/coverage cutoff of >90% and E-values < 10⁻⁵. BLASTn Atlas was generated with Gview, applying as parameters identity >70% and E-values < 10⁻⁵ (Petkau et al., 2010).

The bla_KPC−2-harboring 10,951 bp plasmid of KPC-142 (pKQPS142b) was compared with plasmids of K. pneumoniae strain Kp13 (pKp13d), K. pneumoniae strain A60136 (p60136), and E. coli (IncQ pRSF1010) by BLASTn. The genomic context of the highly similar genes to pKQPS142b were obtained using the Gview program (Petkau et al., 2010), with manual editing.

Kaptive tool (Wyres et al., 2016) was used for classification of capsular loci (cps cluster) and wzi and wzc allele typing was performed by querying the predicted nucleotides sequences of both genes against the Institut Pasteur MLST database [bigsdb.pasteur.fr/klebsiella/].

The allele sequences for seven housekeeping genes (gapA, infB, mdh, pgi, phoE, rpoB, and tonB) and sequence types (STs) were assigned by using the MLST (http://bigsdb.pasteur.fr/klebsiella/klebsiella.html) and NCBI (www.ncbi.nlm.nih.gov) databases.

In order to investigate the phylogenetic placement of KPC-142, the seven housekeeping genes were concatenated for each of the 36 genomes used in a bacterial strain panel, in addition to K. quasipneumoniae subsp. quasipneumoniae, K. quasipneumoniae subsp. similipneumoniae, and K. pneumoniae strains, which sequence typing were closer to that identified to KPC-142. K. variicola was used as outgroup. The sequences were aligned in MAFFT using the accurate LINSI strategy (Katoh and Standley, 2013). Poorly aligned regions were removed in the final alignment by Gblocks program (Castresana, 2000), allowing positions with a gap in less than 50% of the sequences. Twenty-four percent of the 12,219 original base positions were retained in the final alignment.

Maximum likelihood tree for the matrix consisting of the seven concatenated loci was constructed with PhyML (Guindon et al., 2009), comprising 67 nucleotide sequences. jModelTest v. 2.1.10 (Darriba et al., 2012) was used to select the best-fit model of nucleotide substitution according to the corrected Akaike information criterion measure. The evolutionary history was inferred using GTR+I+G as nucleotide substitution model, which was chosen as best-scoring in the earlier step, with 1,000 replicates of a nonparametric bootstrap as clade support. Tree editing and annotation were performed using interactive Tree of Life (iTOL) (Letunic and Bork, 2016). A phylogenetic network using the same aligned sequences as input was inferred using the NeighborNet algorithm and untransformed distances (uncorrected-p) in SplitsTree4 v. 4.14.6 (Huson and Bryant, 2006). The pairwise homoplasy index (phi) (Bruen et al., 2006) was calculated within SplitsTree4 to evaluate if the recombination events were statistically significant.

The Bioproject accession number for KPC-142 is PRJNA383559 and the complete sequences of the chromosome and plasmids pKQPS142a and pKQPS142b have been deposited in GenBank/NCBI under accession numbers CP023478, CP023479, and CP023480, respectively.

We have sequenced to closure the genome of K. quasipneumoniae subsp. similipneumoniae KPC-142. This isolate was chosen for sequencing following a RAPD analysis of a panel of multiple bla_KPC−2-producing Klebsiella clinical isolates (Figure S1), in which KPC-142 presented a very divergent DNA pattern. The KPC-142 genome consists of one 5,217,996 bp chromosome with 57.84% G+C content (4,856 protein-coding genes) and two plasmids, namely pKQPS142a with 189,707 bp and 50.71% G+C content (199 protein-coding genes) and pKQPS142b with 10,951 bp and 57.32% G+C content (12 protein-coding genes). Also, KPC-142 contains eight copies of 16S rRNAs and 23S rRNA and nine copies of 5S rRNA.

A panel of 35 representative Klebsiella isolates belonging to K. quasipneumoniae, K. pneumoniae and K. variicola was chosen to comparatively investigate the genome of isolate Kqps142. The genomic characteristics of all 36 isolates (including Kqps142) are shown in Table 1, and the selected strains include the type species of the three Klebsiella species compared (and the two K. quasipneumoniae subspecies), and cover a wide range of sequence types (STs; including those closer to Kqps142), and variable composition of resistance determinants.

We calculated pairwise average nucleotide identity (ANI) of KPC-142 chromosome sequence including type strains K. pneumoniae subsp. pneumoniae DSM 30104, K. quasipneumoniae subsp. quasipneumoniae 01A030, K. quasipneumoniae subsp. similipneumoniae 07A044, K. variicola GJ1 and the remaining strains of the bacterial panel described in the material and methods (Table S1).

The results pointed that, although KPC-142 strain was originally reported as K. pneumoniae before sequencing, it did not fit genomic species boundaries threshold set at 96%, and displayed ANI with the K. pneumoniae DSM 30104 type strain of 93.7% (Table S1). Rather, comparison to the type strains of K. quasipneumoniae subsp. quasipneumoniae (ANI of 96.52%) and K. quasipneumoniae subsp. similipneumoniae (ANI of 99.04%) revealed their relatedness (Table S1). These results are in line with those previously reported by Brisse et al. (2014), further supporting the classification of KPC-142 as subspecies Klebsiella quasipneumoniae subsp. similipneumoniae (hereafter referred as Kqps142). Misidentification of Klebsiella species by conventional clinical microbiology laboratory techniques was evidenced recently by Long et al. (2017), who demonstrated that almost 30% of tentatively identified K. pneumoniae strains in their study were in fact K. quasipneumoniae. Furthermore, through this ANI analysis, it can be inferred that the K. pneumoniae strains are slightly more similar to the K. variicola strains (ANI of 94.38%, median among the compared Kv strains) than to the K. quasipneumoniae (Kqq + Kqs) strains (ANI of 93.60%, median among the compared Kqq and Kqs strains) (Table S1). This observation correlates with the findings of Holt et al. (2015), where the K. pneumoniae and K. variicola phylogroups are closer to each other than K. pneumoniae and K. quasipneumoniae phylogroups (Holt et al., 2015).

Regarding the plasmid pKQPS142a with ~190 Kbp, searches using the PlasmidFinder database classified it as a conjugative IncFIB(K) plasmid, more similar to pKPN3 plasmid (50% coverage query, 98% maximum nucleotide identity) that has very narrow host range, limited to a predominantly human-associated sub-clade of Klebsiella (Kaplan et al., 2015). While pKPN3 plasmid is conjugative, pKQPS142a appears to have lost the genes for T4SS, and has an average of 50.71% G+C content, which is lower than the K. quasipneumoniae median of 57% (reported for 103 assemblies currently available in the NCBI/Genome database (www.ncbi.nlm.nih.gov/genome). This might be due to the considerable rearrangements with several mobile genetic elements throughout the plasmid sequence. Interestingly, pKQPS142a has multiple regions with high identity (70% query coverage, 98% maximum nucleotide identity) to plasmid pKPN-332 of K. pneumoniae KPNIH39, which was recovered from nosocomial infection during long-term patient colonization (Conlan et al., 2016) (Figure 1). Several key features are in common between pKQPS142a and pKPN-332, such as the RepFIB-like origin of replication (repA/repC), the ter loci (terZABCDEF and terX) that confer resistance to tellurium as well as resistance to bacteriophage and colicins (Taylor et al., 2002), a gene encoding for creatinase (creatine amidinohydrolase), an important medical enzyme that has been used for clinical diagnosis of renal function (Liu et al., 2015), the copper-resistance operon (pcoABCDRS) that encode an efflux system to both remove the toxic metal in excess and cytoplasmic copper management (Brown et al., 1995), fibrinolysin (plasminogen activator, pla gene) and the silver resistance genes (sil) (Figure 1). The predicted amino acid sequence of fibrinolysin (pla gene KPC142_05211) was highly similar (BLASTp >80%) to the plasminogen activator protease of Yersinia pestis CO92, which has a proteolytic activity responsible for the invasive character of plague and is considered an important virulence factor in this species (Lähteenmäki et al., 2001). In Kqps142 the genetic context of pla gene is characterized by two upstream transposases (KPC142_1834 and KPC142_1835) homologous to IS3 family of Yersinia enterocolitica (NC_008800), which could be involved in recombination events leading to spread of this Yersinia virulence factor among Klebsiella species.

The sil genes represent a major concern, since silver is commonly used in several types of medical devices, e.g., wound dressings, implants, catheter, and endotracheal tubes. While Klebsiella spp. are not primary wound pathogens, they have been linked with ventilator-associated pneumonia (Torres et al., 2009). In a recent Swedish tertiary hospital study with restricted consumption of silver-based products, the presence of sil operon genes was observed at high frequency, exclusively for members of the Enterobacteriaceae family and most common among Enterobacter cloacae and K. oxytoca isolates (Sütterlin et al., 2017). Furthermore, in the case of K. pneumoniae the silver resistance could be selected in a single step (i.e., by a single point mutation) from organisms harboring a repressed sil operon (Randall et al., 2015). Collectively these studies and the results showed here for a nosocomial isolate of K. quasipneumoniae subsp. similipneumoniae emphasize the advisement for the controlled use of silver-based products as well as the monitoring of silver resistance in hospital environments.

On the other hand, pKQPS142a could be considered as a multidrug plasmid, since it also carries the arsenical resistance operon (ars) that encode a specific detoxification pathway for arsenic extrusion (Achour-Rokbani et al., 2010) and a sugE gene encoding for quaternary ammonium compound efflux pump that at high-level expression leads to resistance to a subset of toxic quaternary ammonium compounds (Chung and Saier, 2002) (Figure 1). Also, we highlight another virulence factor encoded on pKQPS142a, i.e., the citrate-dependent iron(III) transport system FecABCDE, FecR, and FecI, where citrate is usually considered an external siderophore (Hussein et al., 1981).

Of special interest in Kqps142 is the plasmid pKQPS142b carrying bla_KPC−2. We identified that this small plasmid has homology with the broad host-range IncQ1 plasmid pRSF1010 of E. coli (Scholz et al., 1989) in the regions involved in plasmid replication and mobilization, i.e., mob and rep genes (Figure 2). Also, pKQPS142b shows an extensive region of high identity (99% identity, from mobA, mobB, repression protein F, repA, repB to aph(3')-VIa) with the IncQ1-like p60136 plasmid of K. pneumoniae strain A60136. The p60136 carries a novel β-lactamase, namely Brazilian Klebsiella carbapenemase-1 (BKC-1) and was isolated in São Paulo, Brazil in 2009 (Nicoletti et al., 2015). While p60136 carries BKC-1, pKPC142b harbors bla_KPC−2 (KPC142_06049), flanked by an upstream ISKpn6 transposase from IS1182 family and a downstream transposon resolvase Tn3. This plasmid region has around 3 Kbp and showed high similarity to the homologous region of pKp13d of K. pneumoniae Kp13 isolated in Southern Brazil in 2009, whose complete genome has been previously reported by our group (Ramos et al., 2014) (Figure 2).

Nicoletti et al. demonstrated that p60136 is a non-transferable plasmid by conjugation from E.coli transformant cells, but it can be mobilized at a high frequency by helper conjugative plasmids (Nicoletti et al., 2015). However, we found no evidence of T4SS genes in the sequence of the largest plasmid pKQPS142a that would support plasmid-mediated conjugation functions provided in trans.

On the other hand, we found an integrative conjugative element ICEPm1 that is conserved in Proteus mirabilis, Providencia stuartii, and Morganella morganii (Flannery et al., 2009) and it was exclusive when compared to a panel of 35 strains representative of each Klebsiella phylogroup (Table 1 and see R7 in Figure 4). In Kqps142, the ICEPm1 is a 92,220 bp region that extends from a phage integrase encoding gene (KPC142_03753) to a chromosome partitioning-related encoding gene (KPC142_05483), which contains a core segment showing homology to a T4SS for DNA transfer (KPC142_03739 to KPC142_03716). Interestingly, and in a fashion similar to plasmids, ICEs have the ability to transfer other mobile elements in trans (Meyer, 2009). Another ICE, namely R391/SXT family of ICEs of Vibrio cholerae O139 has been demonstrated to mobilize the IncQ plasmid RSF1010 in an oriT-independent manner. This and further studies suggested that other mechanisms must account for this unexpected oriT-independent mobilization of RSF1010, which are still under investigation (Hochhut et al., 2000; Poulin-Laprade et al., 2015). Also, ICEPm1 and STX are compatible, since both ICEs were found in a single bacterium of P. mirabilis HI4320 (Flannery et al., 2011). These findings raise the question whether the ICEPm1 is capable of transferring the pKQPS142b (IncQ plasmid) in trans, further aiding in the dissemination of bla_KPC−2.

Additionally, in pKQPS142b the inverted right repeat (IRR) of transposon Tn3 was identified between ISKpn6 and bla_KPC−2, but as previously reported in pKp13d of K. pneumoniae Kp13 no inverted left repeat (IRL) was detected. This can be regarded as evidence of horizontal spread of bla_KPC−2 gene between K. pneumoniae and K. quasipneumoniae by recombination events involving ISKpn6 and the flanking Tn3-family sequences and facilitated by the broad host-range mobilizable IncQ1 plasmid. In addition, the bla_KPC−2-carrying IncQ pKQPS142b plasmid found here is highly relevant in view of few reports of plasmids associated to the spread of such relevant resistance gene (Carattoli, 2011; Nicoletti et al., 2015; Mollenkopf et al., 2017).

Multilocus sequence typing (MLST) analysis revealed a novel combination of known K. pneumoniae MLST alleles (gapA−18, infB−22, mdh−56, pgi−61, phoE−74, rpoB−38, tonB−99), with ST2736 being the closest sequence type, differing only in pgi (allele 16) and phoE (allele 11) loci, with ST1032, ST1361, ST1413, ST1703, ST2119, and ST2137 differing in three alleles.

Maximum likelihood phylogenetic analysis performed using the concatenated MLST genes showed the existence of three main groups. The first, Clade I, comprising all K. variicola sequences; Clade II grouped only K. pneumoniae strains, while Clade III was composed by strains from both K. quasipneumoniae subspecies, as well as some K. pneumoniae isolates (Figure 3A). The phylogenetic placement of Kqps142 was closer to Clade III strains, which included all K. quasipneumoniae strains, reinforcing its classification as K. quasipneumoniae. Kqps142 grouped with K. quasipneumoniae subsp. quasipneumoniae AJ055, an ST2119 isolate differing by three MLST alleles. Also in Clade III were strains with STs differing by few alleles compared to Kqps142, such as K. pneumoniae ERKP033 (ST1032), K. pneumoniae BK44389 (ST1703), and K. quasipneumoniae subsp. quasipneumoniae K38An (ST2137). In order to better elucidate the complex relationships between these strains, a split decomposition analysis was performed, which yielded a network-like structure suggestive of high levels of recombination in the studied population (Figure 3B and Figure S3). The pairwise homoplasy index (phi) was significant (p-value = 0.0), indicating strong statistical support for the recombination events in the split network. In this network, the three main clades identified in the maximum likelihood tree were recapitulated, however, the relationships within Clade III could be better delineated, as this branch was separated into Clade III-A (grouping both subspecies of K. quasipneumoniae and some K. pneumoniae strains) and Clade III-B (containing only K. quasipneumoniae subsp. quasipneumoniae strains) (Figure 3B). These latter strains could be regarded as closer to the type strain of this subspecies (strain 01A030), while those in Clade III-A, including Kqps142, probably underwent more recent recombination events leading to their clustering in the split network.

In silico capsule typing was performed using a standard nomenclature to refer to the capsular locus type (K-locus are designated as KL) (Wyres et al., 2016) and comparisons against the K. pneumoniae BIGSdb database (http://bigsdb.web.pasteur.fr). This identified that Kqps142 belongs to KL16 type, containing 22 genes with high similarity to the best match reference locus (97.2% identity and 100% coverage) with manB, manC, and wzc allele type 17. The capsular region also carries a novel wzi, a variant that differ to the most closely related allele 60 in two bases: ⁴⁴⁴T⁴⁴⁴G (base depth of coverage 63 without variants) ⁴⁴⁷C⁴⁴⁷T (base depth of coverage 63 without variants), raising the possibility that if this wzi allele is indeed novel, it could associate to the K16 capsular locus.

The antimicrobial susceptibilities of Kqps142 are summarized in Table 2. The isolate presented resistance to all tested beta-lactam antibiotics (ampicillin, ampicillin-sulbactam, cefepime, cefoxitin, ceftazidime, ceftriaxone, cefuroxime, ertapenem, imipenem, meropenem, and piperacilin-tazobactam). Susceptibilities were identified for amikacin (aminoglycoside), gentamicin (aminoglycoside), ciprofloxacin (fluoroquinolone), tigecycline, and colistin (polymyxin E). Correlation between phenotypic antibiotic susceptibility and the resistome, particularly for beta-lactam resistance, showed that besides KPC-2 (gene locus KPC142_06049, found in pKQPS142b), Kqps142 also harbors the chromosomal gene encoding for OKP-B-6 (KPC142_02249). While pKQPS142b carries aph(3')-Via (synonymous to aphA6), the MICs of the aminoglycoside antibiotics tested in this strain were low (amikacin, MIC 8 mg/L and gentamicin, MIC ≤ 1 mg/L). Similarly, Yoon et al. also observed an unexpected amikacin susceptibility of the aphA6-carrying Acinetobacter guillouiae strains, and the authors correlated this phenotype with low level of gene expression (Yoon et al., 2014), which remains to be tested in the isolate KPC-142.

We compared the chromosomal gene content between Kqps142 and strains belonging to K. quasipneumoniae subsp. similipneumoniae (phylogroup II-B), K. quasipneumoniae subsp. quasipneumoniae (phylogroup II-A), K. pneumoniae (phylogroup I), and K. variicola (phylogroup III), supporting a panel of 36 strains representative of each Klebsiella phylogroup (Holt et al., 2015).

Proteinortho co-orthologous clustering analysis revealed a total of 1,049 core genes conserved in all 36 genomes (Table S2), whereas we found 1,597 genes that were shared by ≥ 97% of the genomes that composed our bacterial panel of strains belonging to phylogroups KpI, KpII, and KpIII. We identified a pangenome of 18,330 unique protein-coding sequences among the 36 analyzed Klebsiella genomes. These results are commensurable with those shown by Holt et al. (2015), who identified a pangenome of 29,886 unique protein-coding sequences among 328 K. pneumoniae genomes and revealed an open pangenome (Holt et al., 2015). Concerning the “exclusive” genes, a total of 9,795 genes had no orthologos among the 36 genomes (yielding a median of 272 “exclusive” genes per genome). Particularly, this comparative analysis revealed 181 “exclusive” genes in the genome of isolate Kqps142, most of which are genes of unknown function, or genes related to ICEPm1 and capsule biogenesis (wzb, wzc, wzx, and wzy) (Table S2).

Major chromosomal features, common or exclusive, among compared Klebsiella strains are depicted in Figure 4. Genome Atlas shows that Kqps142 harbors seven main regions of genomic plasticity in comparison to other 35 Klebsiella genomes, almost all of which with variations in G+C content (except R5). R3 presents the capsular polysaccharide biosynthesis genes, while R1, R2, R4, and R6 are composed by ORFs involved in phage integration. R5 includes the allantoin operon and R7 contains phage sequences and the integrative conjugative element ICEPm1.

Virulence genes identified in the 36 genomes from each Klebsiella phylogroup, including the isolate Kqps142, are shown in Figure 5. These included genes associated with the biosynthesis of siderophores, fimbriae, lipid A, capsule, or genes encoding for microcin, allantoinase, and the ferric uptake operon kfuABC, which have been reported as virulence factors of K. pneumoniae on the basis of murine models of infection (Paczosa and Mecsas, 2016). Most of these genes were detected in all Klebsiella phylogroups, since the majority of compared Klebsiella strains were isolated from human infections (Table 1). However, a few genes were identified in just one or two Klebsiella phylogroups, namely those genes encoding for the siderophore yersiniabactin (found in five K. pneumoniae and one K. variicola), capsular gmd gene (found in only one K. pneumoniae and three K. quasipneumoniae subsp. quasipneumoniae), rfbB gene (found only in six K. pneumoniae), and the N-acetyl-neuraminic acid gene cluster (found exclusively in all K. quasipneumoniae subsp. similipneumoniae) (see also Table 3).

Table 3 shows the gene content related to the most relevant metabolic pathways among 36 strains representative of each Klebsiella phylogroup, including the isolate Kqps142. Regarding the protein-encoding genes shared only by K. quasipneumoniae subsp. similipneumoniae strains are of special interest the N-acetyl-neuraminic acid gene cluster encompassing nanTEAR (KPC142_01465 to KPC142_01468) and nanK (KPC142_01464), which are shared by the nine analyzed strains of K. quasipneumoniae subsp. similipneumoniae (Table 3, Figure 4 and Table S2). Blin et al. (2017) showed that the ability to use N-acetyl-neuraminic acid (Neu5Ac) could be a phenotypic marker for K. quasipneumoniae subsp. similipneumoniae (Blin et al., 2017). Neu5Ac is especially abundant in the epithelial mucus of eukaryotic host and bacteria can scavenge this canonical sialic acid, for either catabolism and/or sialylation of their cell surface, which serves as a key determinant of pathogenesis (Haines-Menges et al., 2015) (see also Figure S4).

The metabolism of allantoin is used by bacteria as a source of carbon and nitrogen and is considered an important determinant of virulence in K. pneumoniae (Paczosa and Mecsas, 2016), since the operon genes from allS to allD are upregulated in hypervirulent (HV) K. pneumoniae strains compared to classical strains (Liu et al., 1986; Chou et al., 2004). The presence of allantoin operon genes among the 36 analyzed Klebsiella strains was very rare, limited only to Kqps142, K. quasipneumoniae subsp. quasipneumoniae strains FIHV2014 and PO1285 and K. pneumoniae strains 1084 and NTUH K2044 (Table 3 and Figure 4). Indeed, strain NTUH K2044 causes liver abscess and meningitis (Wu et al., 2009), strain FIHV2014 was described as a hypermucoviscosity-positive clinical isolate (Arena et al., 2015) and although 1084 was described as a hypermucoviscosity-negative clinical isolate it was capable of inducing liver abscesses in a murine model (Lin et al., 2011, 2012). On the other hand, while strains PO1285 (NCBI BioProject PRJNA351846) and Kqps142 are clinical isolates from human infections, their virulence in vivo remains to be investigated.

The genes involved in the metabolism of D-tagatose were not found among nine K. quasipneumoniae subsp. quasipneumoniae genomes and were rare among Kqq, Kqs, and Kpn strains (Table 3 and Figure 4). These genes encode for tagatose-bisphosphate aldolase subunit GatY (EC number 4.1.2.40) (KPC142_01383) and the subunit GatZ (KPC142_01387) (Table S2) that functions as a chaperone-like for the proper and stable folding of GatY (Brinkkötter et al., 2002). Fermentation and substrate assimilation tests were previously performed showing that K. quasipneumoniae subsp. similipneumoniae strains produce acid from D-tagatose and can grow on this sugar (Brisse et al., 2014; Blin et al., 2017).

The tcuRABC genes that encode functions needed for the utilization of tricarballylate as a carbon and energy source were absent only in the nine analyzed K. pneumoniae strains (Table 3 and Figure 4). These genes were described in the Salmonella enterica serovar Typhimurium LT2 (Lewis et al., 2004) and the utilization of tricarballylic acid was experimentally demonstrated as very rare in K. pneumoniae (Blin et al., 2017). A putative operon involved in the ability to use D-psicose was more frequent in the K. pneumoniae strains and it was also identified in two K. quasipneumoniae subp. quasipneumoniae genomes (UCICRE14 and FIHV2014, Table 3), which are in agreement with reported phenotypic findings (Blin et al., 2017).

An operon (phnUTWX) involved in the metabolism of 2-aminoethylphosphonate (2-AEP) as phosphorus source was only absent in the nine analyzed K. variicola strains (Table 3 and Figure 4). In Enterobacter aerogenes and Salmonella enterica subsp. enterica serovar Typhimurium, 2-AEP catabolism might be controlled as part of the Pho regulon and these bacteria can use 2-AEP as an alternative source of phosphate only under conditions of P limitation (Jiang et al., 1995), and in Pseudomonas putida str. NG2 and Pseudomonas aeruginosa str. PAO1 2-AEP can be mineralized as the sole C, N and P source (Ternan and Quinn, 1998; McGrath et al., 2013).

In bacteria, the tripartite tricarboxylate transporter (TTT) TctCAB is involved in the citrate uptake, which then can be utilized as a carbon and energy source by bacterial cell. The system was well-studied in Bordetella pertussis (Antoine et al., 2003), Salmonella typhimurium (Widenhorn et al., 1989), and Corynebacterium glutamicum, a nonpathogenic species that is of interest due to its biotechnological importance as a producer of L-glutamate and L-lysine (Brocker et al., 2009). Here, we found the tctCAB operon genes were restricted to strains belonging to K. quasipneumoniae species (Table 3 and Figure 4).

Also, protein-encoding genes involved in two pathways, i.e., 1,2- propanediol (1,2-PD) and ethanolamine (Eut) utilization were found in almost all compared strains (except for K. pneumoniae DSM 30104, Table S2 and Figure 4). Both the coenzyme B12-dependent catabolism of 1,2-PD (PduL EC: 2.3.1.222, KPC142_03505) and B12-dependent degradation of ethanolamine (eut operon) use the bacterial microcompartment (BMC) structures encoded by pduA (KPC142_03506) and eutN/ccmL (KPC142_03510), respectively. The BMC function is to optimize biochemical pathways by confining toxic by-products or volatile metabolic intermediates (Bobik et al., 2015). The ethanolamine and 1,2-PD degradations seem to be important to enteric pathogenesis (Thiennimitr et al., 2011) and their relevance in the pathogenesis or fitness of Klebsiella species still remains to be experimentally examined.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.