P. aeruginosa isolates (n = 635) were collected between March-2018 and April-2019 as part of an antimicrobial resistance surveillance program that was initiated in 2018 by the Infectious Diseases Research Department (IDRD) at the National Guard Health Affairs (NGHA) to monitor the prevalence and trends of resistance in a variety of clinically important pathogens. The program involves the monthly collection of the first 10–30 non-duplicate consecutive isolates of each surveyed bacteria identified in the laboratories of NGHA medical cities located in Riyadh, Jeddah, Al Madinah, Dammam and Al Ahsa. The collection included in this study comprised 162 isolates referred from King Abdulaziz Medical City—Riyadh (Centre—Riyadh province, 1,500 bed facility), 275 from King Abdulaziz Medical City—Jeddah (West—Makkah province, 750 bed), 67 from Prince Mohammed Bin Abdul Aziz Hospital—Al Madinah (West—Al Madinah province, 215 bed), 86 from King Abdulaziz Hospital—Al Ahsa (East—Eastern province, 300 bed) and 45 from Imam Abdulrahman Al Faisal Hospital—Dammam (East—Eastern Province, 100 bed). Isolates were recovered from urine (208/635, 32.8%), respiratory (204/635, 32.1%), blood (101/635, 15.9%), wound (67/635, 10.6%) and other specimens (55/635, 8.7%); they were referred at an overall average of 45 (range 19–80) isolates per month (Table 1). Species identification and antimicrobial susceptibility testing were determined with the VITEK II system. Minimum inhibitory concentrations (MICs) of colistin were confirmed using the micro-broth dilution method. MICs were interpreted according to CLSI breakpoints.

Species identity of the isolates was confirmed with a PCR targeting the species-specific oprL gene. Presence of genes encoding IMP, GES, KPC, NDM, OXA-48-like, VIM, BEL, PER, and VEB β-lactamases were screened by PCR using the primers described in Supplementary Table 1.

Genomic DNA from all isolates (n = 45) was extracted with the MagnaPure compact system (Roche, Switzerland) and prepared for sequencing with the Nextera XT DNA library preparation kit (Illumina, United Kingdom) according to the manufacturer’s instructions. Sequencing was performed on the Miseq instrument using the 2 × 300 paired-end protocol. Of these, nine isolates were further sequenced on the Oxford Nanopore MinION using the ligation sequencing kit according to the manufacturer’s instructions (Oxford Nanopore Technologies, United Kingdom). Genome assemblies using the Illumina reads alone or in combination with the Nanopore long-reads were generated using Unicycler 0.4.8 (Wick et al., 2017). Multilocus sequence type (MLST) was determined in silico using the mlst-v2.18.1 software. Genes, mutations associated with antimicrobial resistance, and virulence factors were detected with Abricate 0.9.8¹ or Genefinder.² In order to put the analysis into an international context, all P. aeruginosa paired-end Illumina sequenced genomes (n = 16,337) deposited before December 2020 in the NCBI sequence read archive database³ were retrieved and the quality of corresponding reads was assessed with the FastQC software. Relatedness of recovered genomes with read coverage above 20x and those generated in this study was inferred using a single nucleotide polymorphisms (SNPs)-based approach by mapping reads against the publically available sequences of strain PAO1 (NC_002516) or the fully closed genome of strain ST235-MPA32 generated in this study. SNPs were first identified with Snippy 4.4.5⁴ to create a full core genome alignment that was latter checked for recombination events using Gubbins 2.4.1 (Croucher et al., 2015). Filtered alignments were then used to construct the phylogenetic trees using RaXML with the default option of Gubbins. Phylogenetic trees were annotated using iTOL v6 and Microreact tools (Argimón et al., 2016; Letunic and Bork, 2021). SNP locations were determined using the annotated genome of the reference strain PAO1 (AE004091.2). Genome assemblies were annotated with prokka 1.14.6 and their gene contents were compared with Roary (Page et al., 2015). The integrative and conjugative elements (ICE) and direct environments of carbapenemases and extended-spectrum β-lactamases in fully closed genomes was determined manually. Presence of these elements in each sequenced genome was later determined by checking the depth of coverage of reads mapped across the full sequence of each mobile element.

P. aeruginosa clinical isolates were referred over a 14-month period from five hospitals located in the eastern, western and central regions of Saudi Arabia. Susceptibility testing showed that 143/635 (22.52%) isolates were resistant to both imipenem and meropenem (MIC ≥ 8 mg/L) (Table 1). Otherwise, the isolates were variably resistant to ciprofloxacin (16.22%, range 4.44–20.90%), ceftazidime (14.17%, range 6.67–19.14%), cefepime (8.53%, 4.44–11.90%), piperacillin/tazobactam (16.33%, range 4.65–33.75%), amikacin (6.93%, range 2.33–11.94%), gentamicin (8.06%, range 2.35–13.43%) and tobramycin (7.61%, range 0–15.15%) but remained highly sensitive to colistin (96.17%, range 90.91–100%) (Table 1). Resistance to carbapenems was highest in respiratory isolates (40.69%) followed by blood (27.72%), urine (12.02%) then wound (8.96%) swab isolates. Extensively drug-resistant isolates, remaining only sensitive to colistin, accounted for 3.78% (24/635) of all isolates and were predominantly obtained from respiratory specimens (14/24, 58.33%). Regional and individual hospital variations in susceptibility patterns and frequencies were observed for most tested antibiotics (Table 1). Characteristically, hospitals hosting critical patients (i.e., Jeddah and Riyadh) experienced the highest levels of resistance across all tested antibiotics.

PCR screening detected genes encoding GES (n = 18), VIM (n = 6), NDM (n = 1), and OXA-48-like (n = 1) β-lactamases, explaining resistance to imipenem and meropenem in 18.18% (26/143) of carbapenem-resistant isolates. Genes encoding VEB (n = 6) and PER (n = 2) extended-spectrum β-lactamases (ESBLs) were detected in only eight isolates, including three of the six VIM-positive isolates. Most common GES- and VIM-positive isolates were referred from all regions at different sampling periods (i.e., 11 out of 14 sampling months) of the survey.

To further investigate the molecular mechanisms of resistance to carbapenems, all carbapenemase and ESBL producers identified by PCR (n = 31) and a set of randomly selected isolates exhibiting resistance to imipenem and meropenem but with no acquired β-lactamases (n = 14) were whole genome sequenced on the Illumina Miseq system.

Sequence analyses showed that all genomes carrying genes encoding GES type β-lactamases, including GESs -5 (n = 16), −15 (n = 1) and −1 (n = 1), belonged to ST235. The later accounted for 5.03% (822/16,337) of all genomes retrieved from the public databases, and of which, the majority (71.29%, 586/822) carried genes encoding acquired carbapenemases, including IMP alone (n = 264) or in combination with NDM (n = 1) or OXA-48-like (n = 1), VIM alone (n = 133) or in combination with OXA-48-like (n = 3), GES (n = 174), KPC (n = 9), and OXA-48-like (n = 1) (Table 2). Full genome SNP-based phylogeny against the P. aeruginosa strain PAO1 reference genome grouped nearly all published GES-positive ST235 genomes and those generated in this study in a distinct cluster with the exception of three GES-20 producers which grouped apart (Figure 1). In contrast, the two other common VIM and IMP carbapenemases were distributed across multiple clusters suggesting that they were acquired through multiple events. The phylogenetic tree constructed based on SNP calls relative to the fully closed ST235 MPA32-genome confirmed the clustering of the GES-positive genomes in a well distinct clade (Figure 2). Clustered GES-positive genomes originated from at least six different countries, including Australia, China, Germany, Japan, Indonesia and Pakistan (Figure 1). Phylogeny showed that these producers were related to each other, with at most 120 SNPs to distinguish between them. More specifically, the GES-positive genomes (n = 18) from this study clustered tightly with others from Germany (n = 26) in a distinct subgroup with at most 45 SNPs to distinguish them from each other. GES-positive genomes were separated from the remaining ST235 genomes by at least 29 SNPs, of which, 14 were non-synonymous in genes classified as transcriptional regulators, metabolic genes and hypothetical proteins whereas four were located in the promoter or potential regulatory regions of genes encoding the porin OprO, global transcriptional regulator IscR or hypothetical proteins (Supplementary Table 2). Of these, only IscR has been shown to regulate genes involved in iron homeostasis, resistance to oxidants and pathogenicity (Romsang et al., 2014; Saninjuk et al., 2019). Comparison of the gene contents of ST235 genomes identified 4,373 core and 1,086 soft-core genes that were shared among 98.94–100% and 94.87–98.94% of genomes studied, respectively. Other 25,228 genes were found to be present within one to a maximum of 536 genomes (i.e., < 94.87%) thus showing a relatively wide genetic variability in the accessory genome of this ST. Overall there were a very limited number of genes (n = 27) which presence or absence distinguished the GES-positive cluster from the remaining ST235 but these mainly encoded phage or hypothetical proteins. Only one gene showing homology to the transcriptional repressor NrdR lacked in all GES producers but was present in the majority of remaining ST235 genomes. The latter, shown to regulate the ribonucleotide synthesis may grant the adaptability to thrive in different environments (Crespo et al., 2015). Otherwise, the presence of genomic features encoding the virulence factors gathered in the virulence factor database (VFDB) were comparable across all ST235 genomes. The majority of the genes and gene clusters previously shown to be involved in the species virulence were identified in nearly all (i.e., 98–100%) ST235 genomes with the exception of wzz and wzy genes that are involved in the B-band lipopolysaccharide O antigen synthesis and the pyoverdine outer membrane receptor fpvA. In accordance with previous studies, all ST235 genomes including those carrying bla_GES carried the ExoU toxin-encoding gene (Sato et al., 2003; Treepong et al., 2018).

Sequence examination of the fully reconstructed genomes of two GES-5 producers co-located the β-lactamase gene with acc(6′)-Ib, aph(3′)-XV, aadA6, florR, sul1, tetG, and catB genes in a type I integron that was embedded in a 95 kb integrative conjugative elements (ICE) inserted in the chromosome downstream the tRNA-Gly gene (Figure 3). Mapping of the short reads to the full sequence of this element confirmed its presence in nearly all publically available (166/174, 95.4%) and newly sequenced (18/18, 100%) GES-positive genomes. Of the remaining, five (5/174, 2.87%) genomes carrying bla_GES–1 (n = 1) and bla_GES–5 (n = 4) lacked the region located between position ∼65.1 and ∼72.1 kb and comprising the floR and tetG genes whereas the three genome carrying bla_GES–20 (3/174, 1.73%) and clustering apart had the entire element missing (Figure 3). One GES-negative isolate (i.e., RPA66) sequenced in this study belonged to ST235 clustered also away, thus confirming the relatedness of the GES-producing ST235 isolates (Supplementary Tables 3, 4).

Sequence analysis showed that the majority of VIM-positive genomes (5/6, 83.3%) carried bla_VIM–2 and belonged to ST233 (n = 3) or ST357 (n = 2). The remaining isolate (1/6, 16.7%) harbored bla_VIM–28 and belonged to ST111. Isolates carrying bla_VIM–2 were from three hospitals located in the three different regions included in the study (Supplementary Table 3). Here also, nearly all published ST233 genomes (90/94, 95.74%), which originated from at least four different countries, including Germany, Japan, Spain and United States, harbored bla_VIM–2 and were highly related with at most 47 SNPs to distinguish them from each other. Isolates from this study clustered in two sub groups according to their geographical origins, nevertheless with 29 SNPs to distinguish them from each other (Figure 4). Sequence analysis located the bla_VIM–2 of all newly sequenced ST233 isolates in a 4.4 kb class I integron comprising aac(6′)-I1, dfrB-5, and aac(3′)-Id resistance genes (100% identity, accession number AY943084.1), and which, in turn, was embedded in a transposon similar to those found in published P. aeruginosa genomes (e.g., CP056774.1 from nucleotide positions 5,258,523–5,270,622) (Figure 5). Screening of publically available bla_VIM–2-positive ST233 genomes showed that nearly all (89/90, 98.89%) had reads covering the entire 4.4 kb integron cassette and of which the majority (72/90, 80%) had also reads covering more than 95% of the entire 12 kb transposon. On the other hand, only a handful of publically available ST357 genomes harbored the bla_VIM–2 (8/229, 3.49%) and had mainly IMP carbapenemases (46/229, 20.09%) and in particular the IMP-7 variant (42/46, 95.45%). The bla_VIM–2 of newly sequenced ST357 genomes was detected in an integron type I mobile element embedded in a complex transposon containing all the resistance genes identified in this isolate (Figure 5). In contrast, the bla_VIM–28 in the ST111 isolate was located on a 350 kb plasmid in a integron type I in association with aac(6′)-I.

Of the remaining carbapenemase-positive isolates, one carried bla_OXA–232 and belonged to ST244 while one had bla_NDM–1 and belonged to ST773. Genome assemblies located the bla_OXA–232 gene on a 6.1 kb ColKP3-type non-conjugative plasmid similar to previously published plasmid pColKP3 (accession number CP036323). Otherwise, the bla_NDM–1 was co-located with floR2, rmtB, sul1, and tet(G) in an ICE element of approximately ∼117 kb that was highly similar to the clc-like ICE recently identified in P. aeruginosa (accession number MK497171) (Figure 3).