Antimicrobial Resistance and Genomic Characterization of Six New Sequence Types in Multidrug-Resistant Pseudomonas aeruginosa Clinical Isolates from Pakistan

Pseudomonas aeruginosa (P. aeruginosa) is a major bacterial pathogen associated with a variety of infections with high mortality rates. Most of the clinical P. aeruginosa isolates belong to a limited number of genetic subgroups characterized by multiple housekeeping genes’ sequences (usually 5–7) through the Multi-Locus Sequence Typing (MLST) scheme. The emergence and dissemination of novel multidrug-resistant (MDR) sequence types (ST) in P. aeruginosa pose serious clinical concerns. We performed whole-genome sequencing on a cohort (n = 160) of MDR P. aeruginosa isolates collected from a tertiary care hospital lab in Pakistan and found six isolates belonging to six unique MLST allelic profiles. The genomes were submitted to the PubMLST database and new ST numbers (ST3493, ST3494, ST3472, ST3489, ST3491, and ST3492) were assigned to the respective allele combinations. MLST and core-genome-based phylogenetic analysis confirmed the divergence of these isolates and positioned them in separate branches. Analysis of the resistome of the new STs isolates revealed the presence of genes blaOXA-50, blaPAO, blaPDC, blaVIM-2, aph(3′)-IIb, aac(6′)-II, aac(3)-Id, fosA, catB7, dfrB2, crpP, merP and a number of missense and frame-shift mutations in chromosomal genes conferring resistance to various antipseudomonal antibiotics. The exoS, exoT, pvdE, rhlI, rhlR, lasA, lasB, lasI, and lasR genes were the most prevalent virulence-related genes among the new ST isolates. The different genotypic features revealed the adaptation of these new clones to a variety of infections by various mutations in genes affecting antimicrobial resistance, quorum sensing and biofilm formation. Close monitoring of these antibiotic-resistant pathogens and surveillance mechanisms needs to be adopted to reduce their spread to the healthcare facilities of Pakistan. We believe that these strains can be used as reference strains for future comparative analysis of isolates belonging to the same STs.


Introduction
The spread of antibiotic resistance is widely recognized but the data regarding their presence, source, and significance is still limited. The clinical care of patients have been complicated due to the emergence and spread of multidrug-resistant (MDR) and extensively drug-resistant (XDR) bacterial pathogens [1]. The emergence of new clones of hypervirulent antibiotic-resistant Pseudomonas aeruginosa (P. aeruginosa) strains pose serious clinical concerns. The development of drug resistance in persistent P. aeruginosa infections is primarily due to accumulation of pathoadaptive mutations or acquisition of antibiotic resistance genes (ARGs) through horizontal gene transfer [2]. The spread of these antibioticresistant P. aeruginosa strains is increasing within communities and hospital environments, leading to severe infection cases and eventual therapeutic impasse [3]. Therefore, the surveillance and genomic investigation of this pathogen is essential to better understand the epidemiology and evolution of antimicrobial resistance [4].
The advent of high-throughput sequencing technologies along with bioinformatics analysis has provided deep insights into bacterial pathogenesis including the polymorphism associated with adaptation to specific niches and increased antimicrobial resistance [5]. Multi-locus sequence typing (MLST) is a strain-typing scheme that focuses exclusively on conserved DNA sequences of internal fragments of multiple (usually 5-7) housekeeping genes [6] and is one of the widely accepted method with growing interest in re-application of this scheme in post-genomic era.
Several recent reports have called immediate attention to the emergence and spread of new sequence types (STs) and clones of antibiotic-resistant P. aeruginosa, with high epidemic risk in hospitals worldwide [7]. Several successful multi-locus sequence types of this pathogen are more prevalent worldwide despite its non-clonal epidemic population structure. For instance, clones with sequence types ST175, ST235, ST253, ST111, and ST274, in addition to the Liverpool epidemic strains (LES-1, LESB58), first appeared at only one location and later spread globally, causing high mortality rates [8]. Moreover, discrimination of bacterial isolates is more ideal on the basis of MLST for the comparative analysis of strains, regardless of their source or region. The development and standardization of the MLST database have also enabled the comparative analysis of allelic profiles and the identification of new/unique sequence types of bacterial pathogens [9]. MLST data has been widely used for epidemiological investigations to infer population structure of bacterial isolates [10].
In this study, we sequenced 160 non-duplicate P. aeruginosa isolates from a tertiary care hospital in Pakistan, and found six isolates with new allelic combinations of seven housekeeping genes while performing in-silico MLST analysis. We aimed to investigate the key genetic determinants behind high antibiotic resistance and virulence in these isolates. The phylogeny and clonal relationship of the new ST study isolates to the globally disseminated P. aeruginosa population was further investigated. MLST-based phylogenetic investigation also confirmed the uniqueness of the MLST profile of these isolates and positioned them in separated branches. The reported strains can be used as a reference strains for isolates belonging to the same STs for comparative genomic investigations.

Antibiotic Susceptibility
Antimicrobial susceptibility testing was carried out following Clinical and Laboratory Standards Institute (CLSI) M100-S27 guidelines [11]. According to the Multi-resistance classification guidelines, PA_64, PA_65, PA_107, and PA_152 strains were identified as MDR, whereas PA_88 and PA_141 were identified as XDR strains and grouped into the poor outcome strains (Table 1). All strains were sensitive to Colistin. The strain PA_88 was confirmed as a carbapenemase and metallo-beta-lactamase producer by the carbapenem inactivation assay (CIA).

Sequencing Statistics and Pattern of Gene Distribution
The genome size of the study strains ranges from 6.2 to 6.5 Mbp. The median N50 value was 308,371 bp (interquartile range (IQR), 204,338 to 411,114 bp), while the median number of contigs (length > 500 bp) were 68.6 (IQR, 48 to 102). The FastANI tool further confirmed species-level identification of all isolates with percent identity values of 98.5% for all assemblies to the reference strain PAO1. Detailed sequence statistics and genome features are summarized in Table 2. The quality criteria were met for all strain sequences. A total of 7517 orthologs were detected in the six genomes of this study and the reference strains (PAO1 and UCBPP-PA14) (Supplementary Table S1). Among these orthologs, 5065 genes belonged to the core genes, 1113 genes to the shell genes (genes present in two or more strains) and 1339 genes to the cloud genes (genes only found in a single strain).
Several pseudomonas-derived cephalosporinase (PDC) variants (formerly known as AmpC) were also detected, i.e., PDC-266 in PA_64, PDC-127 in PA_65, PDC-109 in PA_107, PDC-3 in PA_152 and PDC-66 in PA_88 and PA_141 ( Figure 1). Indeed, none of the isolates harbor the wild-type PAO1 AmpC sequence (PDC-1). The CATB family gene catB7, which gives resistance to chloramphenicol and fosA gene causing resistance to fosfomycin, was present in all the studied strains. Interestingly, all isolates harbored the mutated pmrA gene of the pmrA-pmrB two-component system, except PA_152. In addition, we observed several strain-specific antibiotic resistance genes, i.e., dfrB5, aac(6 )-II, aac(3)-Id, blaVIM-2, and adeF in PA_88; merP in PA_65; and crpP in the strain PA_152. In addition, several types of non-synonymous mutations were observed in the core genome of all new ST isolates compared to the reference genome PAO1 (Table 4). These mutations include single nucleotide polymorphisms (SNPs), multi-nucleotide polymorphisms (MNPs), complex variations, deletions, and insertions. The total variations in all the isolates ranged from 75,309 in isolate PA_88 to 26,583 in isolate PA_107. There was a median of 47,390 variations in the genomes of new ST isolates. Based on core genome phylogeny grouping, both XDR isolates within group 3 (PA_88 and PA_141) had the most variations (75,309-75,245) and SNPs (66,629). In terms of the number of SNPs and the specimen type, cystic fibrosis (CF) strain PA_88 had exceptionally high numbers of SNPs (66,754) and variant complexes (7325). Relatively more variants were observed in isolates from CF and bronchitis specimens. Non-synonymous mutations were also assessed in genes conferring antibiotic resistance to P. aeruginosa isolates (Supplementary Table S2). There were exceptionally high non-synonymous variations in antibiotic resistance conferring genes in both XDR isolates (PA_88, PA_141). The highest number of mutations in efflux-related genes were observed in mexM, Irf A, mexD, opmE, and mexK. The ampC from the antibiotic inactivation group of shortlisted genes had 12 non-synonymous mutations in both poor outcome strains PA_88 and PA_141. While several mutations were observed in antibiotic target alteration genes in all isolates, remarkably high variations were found in both XDR strains. This includes high mutations in the rosC gene (n= 44 non-synonymous mutations in PA_141 and n = 33 non-synonymous mutations in PA_88). We also found abundant OprD mutations ≥8 in all isolates of this study. All other mutations in antibiotic resistance-associated genes were random without any specific association to sequence type, phylogeny, or susceptibility to antibiotics.
Additionally, several virulence-related genes associated with the type III secretion system, type IV pili, pyoverdin and phenazine biosynthesis, fimbrial biogenesis, alginate production and regulation, flagellar biosynthesis, twitching motility and type IV secretion system were detected in all six strains ( Figure 2). Several variations were observed in virulence-associated genes with more significant variations in a set of effector proteins (toxins) of the type III secretion system (exoU, exoS, exoT, and exoY). Among the flagellar genes, i.e., fliC, fliD, fleI/flag, and fliT, more variations were detected in XDR isolates. The presence of another significant virulence gene, namely pldA (phospholipase D gene), was observed with no clear association to the source or antibiotic resistance profile of the strains. Notable sequence variations in the pvdE, lepA, and toxA genes were also observed in the new ST P. aeruginosa isolates compared to the PAO1 reference strain. All strains, except PA_64, were harboring multiple insertion sequences. These include ISPa6, ISPa1, ISPa32, and ISPa1328 in strain PA_65; ISPa6 and ICE(Tn4371)6041 in carbapenem resistant strain PA_88; ISPa5 in PA_107; ISPa86 in PA_141; ISPa6, ISPa5, ISPa2, ISPa32, ISPa11 and composite transposon cn_2386_ISPa11 in strain PA_152. Several questionable, incomplete, and intact prophages were also detected in the sequenced strains (Supplementary  Table S3). Among these, Pseudo_F10, Pseudo_Phi2, and Pseudo_YMC11/02/R656 were detected in more than three strains. Intact Phage sequence of Pseudo_Pf1 (14.5 kb) was detected in the PA_64 MDR strain. Many other questionable (score 70-90) and incomplete (score < 70) phage replicons were also detected in all the six strains. A few ARGs were also found in phage DNA fragments of MDR isolates, i.e., Pseudo_phi297, Pseudo-PhiCTX, and Klebsi_ST147_VIM1phi7, suggesting the putative role of these phages in the horizontal gene transfer of antibiotic resistance genes. Furthermore, strain-specific prophages, Pseudo_Bacill_G and Pseudo_Dobby were detected in the PA_64, PA_152, and PA_141 strains. Interestingly, more intact prophage fragments were detected in the genomes of MDR isolates compared to the XDR isolates.

Global Phylogenetic Analysis
The core genome SNP-based phylogenetic analysis clustered three strains, namely PA_64, PA_107, and PA_152, into three sub-lineages within group 1, which contains the most commonly studied strain PAO1 [12] and some cystic fibrosis strains such as LESB58 and DK2 [13,14] (Figure 3). The other MDR strain, PA_65, clustered into group 2 with the well-known virulent strain UCBPP-PA14 [15] and MDR Indian ocular isolate VRFPA04 [16]. Both XDR isolates PA_88 and PA_141 clustered into a separate clade with the Carbapenem resistant MDR P. aeruginosa reference strain AR_0446 from the USA. Additionally, the phylogenetic analysis of the concatenated sequences of conserved internal fragments of seven housekeeping genes revealed a tree consisting of two main phylogroups. The isolates PA_107 and PA_64 were clustered into group 1 with the PAO1 reference strain, while PA_65, PA_152, PA_141, and PA_88 clustered into group 2 with virulent MDR strain UCBPP-PA14 (Figure 4).

Figure 4.
A maximum likelihood tree based on variations in housekeeping genes of 174 global Pseudomonas aeruginosa strains (available at NCBI at the time of analysis, 20 January 2020) and six new sequence type study strains. The inner-colored ring represents the region of isolation of P. aeruginosa strains, while the outermost colored ring represents different sequence types (STs). The tree was refined and annotated using iTOL (https://itol.embl.de/).

Discussion
P. aeruginosa infections in hospitalized patients remains an important issue as it is associated with high morbidity and mortality rates in immunocompromised patients [17]. The financial burden to patients and healthcare providers for drug-resistant microorganisms is challenging. The problem is worsening as the bacteria are developing resistance much faster than the development of new drugs [18]. Recent improvements in microbial wholegenome sequencing (WGS) have rendered this technique applicable in clinical microbiology laboratories for infectious disease control and tracking the epidemiology of pathogens [19]. We performed WGS on a collection of antibiotic-resistant P. aeruginosa isolates collected over one year period (Sep-2016-Sep-2017) from a tertiary care hospital in Pakistan (data not shown). During in-silico MLST analysis, we found six P. aeruginosa strains (PA_64, PA_65, PA_88, PA_107, PA_141, and PA_152) having unique allelic profiles of seven housekeeping genes (acsA, aroE, guaA, mutL, nuoD, ppsA, and trpE) used for multi-locus sequence typing of P. aeruginosa strains. Following verification, the new sequence types were assigned to the respective allelic combinations; ST3493 (PA_64), ST3494 (PA_65), ST3472 (PA_88), ST3489 (PA_107), ST3491 (PA_141), and ST3492 (PA_152). A comparative genomic analysis of these new sequence type isolates was performed to identify evolutionary relationships, genomic divergence, profiles of acquired antibiotic resistance, and virulence factors to better understand the new sequence type strains.
The strains were observed to have diverse antibiotic susceptibility profiles based on which four isolates (PA_64, PA_65, PA_107, and PA_152) were classified as MDR and two (PA_88 and PA_141) as XDR isolates. However, all strains were found as phenotypically sensitive to colistin. Both XDR strains were involved in chronic respiratory system infections as P. aeruginosa is the main pathogen with a role in chronic respiratory infections and pulmonary complications [20].
The cumulative genetic information within a set of bacterial genomes represents the pan-genome and its size increases with the increasing number of strains used for pangenome estimation [21]. The pan-genome analysis revealed a total of 7517 pan-genes out of these, 5065 genes were common in ≥99% strains and represent the core genome for the strains of the current study. Previous studies have also reported the core genome sizes of 5316 [22], 5021 [23], and 5233 [24] in different P. aeruginosa strains. The other studies have also used smaller sets of genomes (five to seventeen genomes) and the results are comparable to our study.
The genotypic resistance profile of all strains was in line with the phenotypic resistance profile. In common with the other reference strains of P. aeruginosa, two beta-lactams (blaOXA-50 and blaPAO) and one each for the fosfomycin (fosA), aminoglycosides (aph(3 )-IIb), amphenicols (catB7), and polymyxin B (basS-pmrB) resistance genes were present in all strains of this study. The beta-lactam resistance gene blaOXA-494, may relate to carbapenem resistance in PA_107 in the absence of the carbapenemase encoding gene, while blaVIM-2 was likely responsible for the carbapenem resistance in strain PA_88 [25,26]. CIA results of meropenem-resistant phenotypes also had concordance with the genotypic profile of carbapenem resistance in PA_88 and PA_141 strains.
As expected, we observed the presence of more resistance genes (n = 16) in XDR than in MDR (n = 11) strains of this study. Wide genetic diversity, adaptability, and high antibiotic resistance attributed solely to the acquired antibiotic resistance genes were observed in the PA_88 strain. The PA_88 strain possess five unique resistance genes, namely PDC-66, aac(6 )-II, aac(3)-Id, dfrB5, and pmrA/pmrB, that confer resistance to betalactams, aminoglycosides, trimethoprim, and colistin, respectively. The aminoglycoside modifying enzyme encoding gene aac(6 )-II catalyzes the acetylation of gentamicin and is a significant determinant of gentamicin and tobramycin resistance but not of amikacin in P. aeruginosa [27]. Fortunately, all strains were found to be sensitive to colisitin; however, we found mutated pmrA(L71R) gene in all studied strains except in PA_152. High-level colistin resistance in clinical P. aeruginosa isolates have been reported due to mutations in the pmrA-pmrB two-component system [28]. Several chromosomally mediated AmpC-type variants were identified in all isolates. The presence of AmpC variants suggests a possible role in the reduced susceptibility of strains towards ceftazidime antibiotics. The role of PDC variants have been previously reported for broadening the enzymes' hydrolytic spectrum and facilitating the degradation of ceftazidime antibiotic and reduced susceptibility to cefepime and imipenem [29,30]. Interestingly, none of the variants have been found to be involved in ceftolozane-tazobactam resistance. Only the PA_88 strain was found resistant to the ceftolozane-tazobactam antibiotic, which is most likely attributed to the presence of acquired beta-lactamase blaVIM-2 in the PA_88 strain. These findings are consistent with the previous reports [31]. Moreover, mutations in the nalC (PA3721) gene were also observed in all isolates except in strain PA_152. The mutations in the nalC gene have been considered significant for the overexpression of the MexAB-OprM efflux pump in environmental P. aeruginosa isolates, leading to aztreonam-resistant phenotypes [32]. In addition, we also observed two unique genes, namely merP and crpP, in MDR isolates PA_65 and PA_152, respectively. The merP gene is involved in mercury resistance while crpPgene encodes a novel protein capable of conferring high-level resistance to ciprofloxacin [33]. Studies have also reported the spread of crpP-like genes in a large group of P. aeruginosa isolates from Europe [34]. The acquisition of these crpP genes has been likely mediated by the acquisition of PAGI-like elements [34]. Several virulence factors including exoU, pldA, pvdE, and toxA in ciprofloxacin-resistant (crpP gene-harboring) isolate PA_152 were also present. VFDB search showed that all isolates carried an arsenal of virulence genes involved in invasion, colonization, and extensive tissue damage [35]. All instances of a gene absence were manually confirmed by examining the orthologs from widely studied strains recommended by PGDB [36]. Out of 147 shortlisted virulence-associated genes, notable variations were found among 16 genes. These include T3SS toxins such as exoenzyme T, exoenzyme S, exoenzyme U, and exoenzyme Y, destabilizing the host defense and signaling system [37]. Among the T3SS toxin, PA14 exoU genes were taken as a reference because they are not present in PAO1. The exoU gene was identified only in PA_65 and PA_152 strains. P. aeruginosa strains carrying exoU and multiple antibiotic resistance genes are more virulent and may lead to poor clinical outcomes [38]. The exotoxin encoding exoS (identity ≥ 70%) and exoT (100% identity in PA_64, PA_65, PA_88, PA_152, and 70% identity in PA_107 and PA_141) were the most prevalent toxins in all the strains [39]. We identified exoS in all MDR and XDR strains, whereas exoU was identified only in MDR strains (PA_65, PA_152). Previous studies also reported that exoU is carried by genomic islands and strains possessing the exoU gene showed larger accessory genomes [40]. Relatively larger accessory genome size was observed for both exoU positive isolates, i.e., 1879 and 1816 accessory genes in PA_65 and PA_152, respectively, were observed compared to the accessory genomes of exoU negative isolates (Supplementary Table S4). Seven flagellar genes, namely fliC, fliD, fliS, fleI/flag, flgM, flgN, and fliT, were detected in MDR strains (PA_64, PA_65, PA_107, PA_152) with more sequence similarity (90-99%) to PAO1 compared to XDR strains (PA_88, PA_141). Both XDR isolates had altered fliC, fliD, fliS, fliT, and flag genes that may affect the flagellar function. The flagellar genes are essential for swimming, swarming, and twitching motility of bacteria [41]. Although non-flagellated strains are defective in acute functions, studies have reported that mutations in flagellin-encoding genes may lead to the loss of flagella an important antiphagocytic strategy in chronic infection isolates [42]. Recent studies on CF isolates have shown the down regulation or complete absence of the fliC gene [43]. High sequence similarity (up to 99%) for type VI secretion system gene pldA (phospholipase D gene) was observed between PAO1 and PA_152, and PA_65 and PA_88. The pldA gene is believed to promote chronic infections. Another notable variation was also observed in the pvdE gene, a pyoverdin synthesis precursor and an essential gene in the virulence of P. aeruginosa [44]. Three strains (PA_88, PA_65, and PA_107), irrespective of their source of isolation and MDR/XDR strain profile, had the PAO1 homolog of the pvdE gene, which is involved in enhancing P. aeruginosa invasion by increasing the expression of exoS [45]. Further studies will suggest the possible role of pvdE gene variants in enhancing P. aeruginosa pathogenesis.
Core genome SNP-based phylogenetic analysis clustered all study strains into three main phylogroups. It is evident that the P. aeruginosa population mainly clustered into three phylogenetically distinct groups [46]. Three MDR strains (PA_64, PA_107, and PA_152) were clustered into group 1 and the MDR strain PA_65 was clustered into group 2 (Figure 3). Similarly, both XDR isolates PA_88 and PA_141 were found in a separate cluster with the P. aeruginosa strain AR_0446. This strain could be a taxonomic outlier of P. aeruginosa group 3, alongside the PA7 strain [47]. The PA7 genome has considerable divergence from other P. aeruginosa strains [48].
The MLST-based phylogenetic investigation and distribution of the study isolates with publicly available P. aeruginosa global representative genomes also confirmed the uniqueness of the MLST profile of new sequence type isolates and positioned them in separate branches (Figure 4). Rooted and un-rooted trees generated for all genes' loci and concatenated gene fragments revealed that all six strains belonged to new sequence types and roughly reflect the geographical location of isolates. Our results show that these strains belong to diverse STs and are not similar to the previously described clinical epidemic isolates. The geoBURST analysis represents the population snapshot of the studied isolates with all STs (n= 3162) reported in the P. aeruginosa MLST database, grouped into different clonal complexes (CCs). A clonal complex is composed of at least two or more STs with single-locus variants (SLVs), in which six out of seven loci are shared with at least one other ST member of the complex [49]. ST3493 was found as a SLV of ST3083 and is part of CC1494. ST3494, ST3492, and ST3489 act as singletons. A singleton is defined as an ST that is not grouped into any previously described CC. ST3472 (PA_88) was found as a SLV of ST516 and ST3037 and was identified in CC90. ST3491 (PA_141) was found in CC179 with ST2042 and ST2446. ST2446 has recently been reported as a new ST of antimicrobial-resistant P. aeruginosa in Chinese minks infected by hemorrhagic pneumonia [50]. Interestingly, more intact prophage fragments were detected in MDR compared to XDR isolates.

Genomic DNA Isolation and Illumina Whole-Genome Sequencing
A suspension of~10 colonies from overnight incubated blood agar plates were used for genomic DNA isolation using the QIAamp BiOstic Bacteremia DNA Isolation Kit (Qiagen, Hilden, Germany). The extracted DNA quality was checked by gel electrophoreses using 1% agarose gel and quantified using the Qubit fluorometer dsDNA BR Assay (Invitrogen). Subsequently, 0.5 ng of DNA was used for the preparation of Nextera Illumina sequencing libraries (Illumina, San Diego, CA, USA) [54]. Whole-genome sequencing was performed on the Illumina Nextseq high-output sequencer to obtain~2 million 2 × 150 bp pairedend reads.

Serotyping and Multi-Locus Sequence Typing
Serotypes were determined using the Pseudomonas aeruginosa serotyper (PAst) program to BLAST search O-specific antigen (OSA) biosynthesis gene clusters in query genomes [61].
In-silico multi-locus sequence typing was performed using the MLST 2.0 online server (https://cge.cbs.dtu.dk/services/MLST/). The relationship of new ST isolates with the existing STs in the MLST database was assessed using geoBURST [62]. The geoBURST analysis also provided information about the clonal complexes of the new STs with the other globally prevalent STs.

Identification of Antibiotic Resistance Genes, Virulence Factors and Mobile Genetic Elements (MGEs)
In-silico identification of the acquired antibiotic resistance genes (ARGs) was performed by submitting assembled genomes to the Resfinder and CARD web servers [63,64]. Antibiotic resistance genes are often associated with mobile genetic elements that help their mobility and spread in a bacterial community. The online tool MobileElementFinder [65] was used to analyze the mobilome integrated into the isolate's genomes. Additionally, the sequence of 145 known virulence genes of PAO1 and exoU virulence gene of UCBPP-PA14 (NC_008463.1) were retrieved from the Virulence Factor Database (VFDB) (Supplementary Table S6) to determine their presence in the genomes of the new ST isolates [66]. The BLAST Ring Image Generator (BRIG) [67] tool was used to generate an image that shows the presence/absence and percent identity of virulence genes in the genomes of isolates under study. The PHASTER online server was used for the identification of intact prophage sequences [68].

Variant Calling
Several types of non-synonymous variations were called using snippy 4.6.0 [69]. On the basis of the literature review, a set of 73 genes related to antibiotic resistance in P. aeruginosa were shortlisted to manually examine the presence of non-synonymous SNPs to predict genotypic variations in the resistome (Supplementary Table S2). Only high-quality, non-synonymous SNPs were used for interpretation

MLST and Core-Genome SNPs-Based Phylogenetic Analysis
In order to obtain phylogroup information, 174 publicly available complete genomes of P. aeruginosa were downloaded from the NCBI and Pseudomonas genome database (PGDB) (Supplementary Table S7) [36]. All complete reference genomes were re-annotated using prokka and pangenome analysis of complete genomes (n = 174) along with the study genomes (n = 6) was performed using Roary version 3.8.0 [70]. The core-genome alignment file (ALN file) from the roary output was converted into a maximum likelihood tree with FastTree v2.1.10 and the resulting file was uploaded and annotated using iTOL [71]. Additionally, the MLST-based phylogenetic tree was constructed that strictly focuses on the variation in conserved gene fragments of seven housekeeping genes (acsA, aroE, guaA, mutL, nuoD, ppsA, and trpE). The concatenated sequences of seven housekeeping genes were aligned using MUSCLE and the maximum likelihood (ML) tree was constructed in MEGAX with 1000 bootstrap iteration [72]. The tree was visualized and annotated on the iTOL server [71].

Conclusions
This study reports the emergence of novel sequence types (ST3493, ST3494, ST3472, ST3489, ST3491, and ST3492) in MDR and XDR P. aeruginosa clinical strains isolated from Pakistan. All the new ST strains contain two beta-lactamases, namely blaOXA-50 and blaPAO, as well as fosfomycin resistance gene fosA, aminoglycoside resistance gene aph(3 )-IIb, amphenicols resistance gene catB7, and mutated basS (pmrB) conferring resistance to polymyxin B. A number of non-synonymous mutations in chromosomal genes conferring resistance to major classes of antipseudomonal antibiotics were also observed. The strains also harbor several important virulence factors including exoU, pldA, pvdE, and toxA, with more notable sequence variations in fliC, fliD, fleI/flag, and fliT genes among XDR isolates. These findings provide important information regarding emerging P. aeruginosa clones in Pakistan; however, further studies are needed to better understand the evolution of these new STs strains. We believe that these strains will be used as reference strains for future comparative genomic analyses of other P. aeruginosa isolates belonging to these STs.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/antibiotics10111386/s1. Table S1: Annotation of the 7517 orthologs found in the genomes of six P. aeruginosa isolates and reference strain PAO1. Table S2: Non-Synonymous SNPs detected in 73 genes related to antibiotic resistance in six new ST isolates using PAO1 as reference genome. Table S3: List of phages identified in P. aeruginosa New sequence types genomes. Table S4: Core, Accessory and unique genes characteristics of the new sequence type isolates compared to the P. aeruginosa reference strain. Table S5: Susceptibility Profile antibiotic panel with Cut-off points for Pseudomonas species following CLSI guidelines. Table S6: List of virulence genes of reference strain PAO1 used in the study. Table S7: List of all reference strains used in the study along with their geographic region of isolation.