Genomic surveillance of Pseudomonas aeruginosa in the Philippines, 2013–2014

Pseudomonas aeruginosa is an opportunistic pathogen that often causes nosocomial infections resistant to treatment. Rates of antimicrobial resistance (AMR) are increasing, as are rates of multidrug-resistant (MDR) and possible extensively drug-resistant (XDR) infections. Our objective was to characterize the molecular epidemiology and AMR mechanisms of this pathogen. We sequenced the whole genome for each of 176 P. aeruginosa isolates collected in the Philippines in 2013–2014; derived the multilocus sequence type (MLST), presence of AMR determinants and relatedness between isolates; and determined concordance between phenotypic and genotypic resistance. Carbapenem resistance was associated with loss of function of the OprD porin and acquisition of the metallo-β-lactamase (MBL) gene blaVIM. Concordance between phenotypic and genotypic resistance was 93.27% overall for six antibiotics in three classes, but varied among aminoglycosides. The population of P. aeruginosa was diverse, with clonal expansions of XDR genomes belonging to MLSTs ST235, ST244, ST309 and ST773. We found evidence of persistence or reintroduction of the predominant clone ST235 in one hospital, and of transfer between hospitals. Most of the ST235 genomes formed a distinct lineage from global genomes, thus raising the possibility that they may be unique to the Philippines. In addition, long-read sequencing of one representative XDR ST235 isolate identified an integron carrying multiple resistance genes (including blaVIM-2), with differences in gene composition and synteny from the P. aeruginosa class 1 integrons described previously. The survey bridges the gap in genomic data from the Western Pacific Region and will be useful for ongoing surveillance; it also highlights the importance of curtailing the spread of ST235 within the Philippines.

P. aeruginosa surveillance in the Philippines Chilam et al ment options. These reports coincide with multi-locus sequence type (MLST) ST235, 9-11 the predominant global epidemic clone. The metallo-β-lactamase (MBL) genes bla VIM and bla IMP -usually associated with integrons carrying multiple resistance determinants -have been identified in ST235 P. aeruginosa isolates from Asian countries. [12][13][14] While the resistance rates and profiles of P. aeruginosa in the Philippines have been well characterized, 15, 16 the molecular epidemiology and AMR mechanisms of this pathogen remain largely unknown. Whole-genome sequencing (WGS) can identify transmission patterns, AMR mechanisms and the source of HA infections. 17 In this study, we characterized the clonal relatedness and resistance determinants of P. aeruginosa isolates from the ARSP using WGS.

Bacterial isolates
A total of 7877 P. aeruginosa isolates were collected and tested for resistance by the ARSP from January 2013 to December 2014. Of the 443 and 283 isolates referred to the Antimicrobial Resistance Surveillance Reference Laboratory (ARSRL) for confirmation in 2013 and 2014, respectively, 179 isolates from 17 sentinel sites were selected for WGS, as previously described. 18 Briefly, 113 isolates of carbapenemase-producing P. aeruginosa were selected; also included were 66 available isolates that were susceptible to all antibiotics tested. We used a proxy definition for "infection origin", whereby initial infection isolates collected in the community or on either of the first 2 days of hospitalization were categorized as communityacquired (CA), and isolates collected on hospital day 3 or later were categorized as hospital-acquired (HA). 19

Bioinformatics analysis
Genome quality was evaluated based on metrics generated from assemblies, annotation files and the alignment of the isolates to the reference genome of P. aeruginosa strain LESB58 (accession FM209186), as previously described. 18 Assemblies were produced from short-read Illumina data 18 and from long-read PacBio data with the HGAP v4 pipeline (Pacific Biosciences). A total of 176 isolates yielded high-quality P. aeruginosa genomes and were included in this study.
We derived the MLST of the isolates from the whole genome sequences. The sequence types (ST) were determined from assemblies with Pathogenwatch (https:// pathogen.watch/) and with MLSTcheck v1.007001, and from sequence reads with ARIBA 22 and the P. aeruginosa database hosted at PubMLST. 23 The MLST calls were curated, as previously described. 18 Integrons were detected in the genome assemblies with IntegronFinder. 24 Evolutionary relationships between the 176 isolates were inferred from core single-nucleotide polymorphism (SNP). A core gene alignment was performed with Roary v3.11.3, using the mafft aligner option and minimum percentage identity for blastp of 90%. Evolutionary relationships between 169 isolates from groups 1 and 2 were inferred from SNPs by mapping the paired-end reads to the reference genomes of P. aeruginosa strains LESB58 (ST146, FM209186) or NCGM2_S1 (ST235, AP012280.1). 18 Mobile genetic elements (MGEs) were masked in the alignment of pseudogenomes with a script available at https://github.com/sanger-pathogens/ remove_blocks_from_aln. For the phylogenetic analysis of ST235 genomes, recombination regions detected with Gubbins 25 were also removed. Alignments of SNPs were inferred with snp-sites v2.4.1, 26 and were used to compute pairwise SNP differences between isolates from different patients (minimum n = 3) belonging to the same or to different hospitals, using a script from https:// github.com/simonrharris/pairwise_difference_count. Maximum likelihood phylogenetic trees were generated with RAxML, 27 based on the generalized time reversible (GTR) model with GAMMA method of correction for among-site rate variation and 100 bootstrap replications.
To contextualize the Philippine genomes, we downloaded, assembled and quality controlled global P. aeruginosa genomes with linked geographical and temporal information, collected mainly between 2007 and 2017, for which raw Illumina paired-end sequence data were available at the ENA. A tree of 904 genomes was inferred with FastTree 28 from an alignment of 549 126 SNP positions, obtained after mapping the reads to the complete genome of strain LESB58 and masking regions with MGEs. A tree of 96 global ST235 genomes was inferred with RAxML from an alignment of 1993 SNP sites obtained after mapping the genomes to the complete genome of strain NCGM2-S1, and masking MGEs and recombination regions.
Known AMR determinants were identified with ARIBA 22 and a curated database of known resistance genes and mutations, 29 the Comprehensive Antibiotic Resistance Database, 30 and a custom database of mutations in the quinolone resistance-determining region of the gyrA/B and parC/E genes described for P. aeruginosa. 4 The output for the porin gene oprD was inspected to detect loss-of-function mutations. The oprD sequences were extracted from the whole-genome draft assemblies with blastn, using the oprD sequence from strain PAO1 (accession NC_002516.2, genome positions 1043982-1045314) as a query, then translated in silico to inspect the integrity of the coding frames. A 444 or Chilam et al P. aeruginosa surveillance in the Philippines 442 amino-acid protein that included a START and a STOP codon was considered functional.
The genomic predictions of AMR derived from the presence of known AMR genes and mutations (test) were compared with the phenotypic results (reference), and concordance was computed for each of six antibiotics (1056 total comparisons). Isolates with either a resistant or an intermediate phenotype were considered non-susceptible. An isolate with the same outcome for both the test and reference (i.e. both susceptible or both non-susceptible) was counted as a concordant isolate. Concordance was the number of concordant isolates as a percentage of the total number of isolates assessed.
All project data, including inferred phylogenies, AMR predictions and metadata were made available through Microreact.

Ethics statement
Ethical approval is not applicable. This study uses archived bacterial samples processed by the ARSP. No identifiable data were used in this study.

Demographic and clinical characteristics of the P. aeruginosa isolates
Of the 179 P. aeruginosa isolates submitted for WGS, 176 passed quality control and were confirmed in silico as P. aeruginosa ( Table 2). Patients were aged from under 1 to 96 years, with 27% (n = 47) of the isolates from patients aged 65 years or older. Fifty-eight per cent (n = 102) of the isolates were from HA infections. In terms of specimen type, 53% (n = 94) of isolates were from respiratory samples (tracheal aspirates and sputum).

Concordance between phenotypic and genotypic AMR
Isolates were tested for susceptibility to nine antibiotics representing five classes ( Fig. 1A-C, Table 3). Most isolates were non-susceptible to carbapenems (n = 100), 10 isolates were susceptible to carbapenems but resistant to other antibiotics, and 66 isolates were susceptible to all nine antibiotics ( Table 1). CA infections a Invasive isolates were considered as those obtained from specimen types marked with an asterisk (*). P. aeruginosa surveillance in the Philippines Chilam et al were more frequently associated with susceptible isolates and HA infections with resistant isolates (Fig. 1D, twotailed Fisher's exact test P = 0.000002).
Of the 18 isolates resistant to imipenem and meropenem but not to other β-lactam antibiotics, 17 carried both loss-of-function disruptions in the OprD porin, and disruptions or known non-synonymous mutations in the NalC (A186T, G71E, S209R) and/or NalD (S32N) regula- tors of the MexAB-OprM multidrug efflux pump, suggesting that their resistance is due to a combination of reduced influx and increased efflux of the carbapenem antibiotics ( Fig. 1E). Among the 81 carbapenem-resistant isolates that were also resistant to third-generation cephalosporin ceftazidime and/or fourth-generation cephalosporin cefepime, 67 isolates carried acquired MBL genes bla VIM-2 (n = 61 genomes), bla VIM-6 (n = 1), bla IMP-26 (n = 4) or bla NDM-1 (n = 1), while five carried disrupted oprD genes plus acquired extended-spectrum β-lactamase (ESBL) genes bla PER-1 (n = 3), bla CTX-M-15 (n = 1) or AmpC-like gene bla DHA-1 (n = 1). The remaining eight isolates harboured other β-lactamase genes, but their carbapenemresistance mechanisms remain uncharacterized. Of the 76 isolates susceptible to carbapenems, 75 carried either a full-length OprD porin (444 amino acids) without any known mutations, or a 442 amino acid-long OprD protein with an intact reading frame, while one isolate was missing the STOP codon in the oprD gene.
The overall phenotypic and genotypic concordance was 93.27% for the six antibiotics analysed (

Population structure of P. aeruginosa in the Philippines
The phylogenetic tree of 176 genomes from the Philippines comprises three major groups, 31 group 1 (n = 64) including PA14, group 2 (n = 105) including PAO1 and the more distantly related group 3 (n = 7) including PA7 ( Fig. 2A). All three groups included carbapenemresistant isolates and susceptible isolates, though most isolates in group 2 were susceptible (n = 39, 60.9%) and most in group 1 were resistant (n = 75, 71.4%, Fig. 2B).
The population of P. aeruginosa comprises a limited number of widespread clones selected from a diverse pool of rare, unrelated genotypes that recombine at high frequency. 32 A phylogenetic tree of 169 genomes from groups 1 and 2 showed that the clonal expansions were mostly within the major group 1 -represented by ST235, ST309, ST773 and ST313 (Fig. 2B) -found across multiple hospitals and resistant to multiple antibiotics. Most of the XDR isolates (n = 61, 87%) were found in ST235, ST244, ST309 and ST773, and most (n = 44, 62.8%) carried bla VIM (an MBL that can degrade all anti-pseudomonal β-lactamases except for aztreonam), 1 AAC(6')-Ib (an aminoglycoside acetyltransferase conferring resistance to tobramycin and amikacin), and the non-synonymous mutation T83I on GyrA associated with resistance to fluoroquinolones.
The higher prevalence of ST235 prompted us to look further at this clone. The phylogenetic tree of 49 ST235 isolates comprised two distinct clades with different geographic distribution (Fig. 2C). Clade I (n = 10) was represented in five hospitals in the Luzon (north) and Visayas (central) island groups, while clade II (n = 39) was represented in 10 hospitals from north to south of the country. The phylogeographic structure of the tree and the relatedness between genomes showed evidence of dissemination of ST235 between hospitals. Within clade Ib (Fig. 2C), one genome from hospital NKI differed from two genomes from hospital BRH by seven and eight SNPs, respectively. Within clade IIb (Fig. 2C), the genetic differences between isolates from the same hospital (mean pairwise SNP differences 36.41 ± 20.84, range 0-64) were not significantly different to those between isolates from different hospitals (mean 45.36 ± 8.12, range 29-61, Mann-Whitney U test z-score = -1.49145, P = 0.13622). The close relationships and the common repertoire of resistance genes between isolates from different hospitals support inter-hospital transmission.
The genomes from the hospital VSM (n = 24) formed at least three clusters within clade IIb, two of which exhibited discrete temporal distribution (VSM-2 and VSM-3, Fig. 2C), suggesting that they could represent hospital outbreaks. In agreement with this, the genomes from different patients within clade VSM-3 differed by an average of 11.55 pairwise SNPs (range 0-24). We also identified isolates within VSM-3 that were collected nine or more months apart (Fig. 2C), suggesting that ST235 can either persist in or be reintroduced to the hospital environment.
The distribution of acquired resistance genes and mutations showed that resistance determinants differed between clades I and II, with patterns that were consistent with the acquisition of multiple genes simultaneously by mobile genetic elements. Long-read sequencing of isolate 14ARS-VSM0870, representative of the XDR resistant profile CAZ FEP IPM MEM TZP GEN TOB AMK CIP (marked with an asterisk on Fig. 2C), revealed the acquisition of bla VIM-2 , bla OXA-10 , catB3, aadA1 (ANT(3")-Ia) and acc(6')-Ib within a class 1 integron integrated in the chromosome at position 977 774 (Fig. 2D). The ciprofloxacin resistance gene P. aeruginosa surveillance in the Philippines Chilam et al

P. aeruginosa from the Philippines in the global context
We placed the genomes from our retrospective collection in the global context of 904 contemporary P. aeruginosa public genomes. This collection of public genomes represented 17 countries and 178 STs, with more than 60% of the genomes being from Europe (n = 373) and the United States of America (USA) (n = 205). The Philippine genomes were found throughout the tree, indicating that the P. aeruginosa population captured in our survey largely represents the global diversity of this pathogen. Notably absent were the epidemic clones ST111 and ST175 (Fig. 3A), which, together with ST235, are responsible for MDR and XDR nosocomial infections worldwide.
A more detailed tree of 96 ST235 genomes of global distribution showed three major clades: clade 1 was represented by isolates from Japan, the Philippines, Poland and Thailand (n = 2); clade 2 showed the broadest geographic distribution across four continents and also included isolates from this study (n = 3); clade 3 comprised exclusively isolates from the Philippines (n = 44, Fig. 3B), which raises the possibility that this lineage of ST235 is characteristic to the Philippines; however, introductions from the other globally dispersed lineages may also occur, as shown in clades 1 and 2.

DISCUSSION
In this study, we combined WGS and laboratory-based surveillance to characterize susceptible and resistant P. aeruginosa circulating in the Philippines in 2013 and 2014, with a particular emphasis on resistance to carbapenems, which increased in the years preceding this survey. Drug-resistant P. aeruginosa infections are difficult to treat, resulting in poor patient outcomes. In a tertiary hospital in Manila, severity of illness and mortality rates were significantly higher among patients infected with drug-resistant P. aeruginosa than among those infected with susceptible isolates, while median duration of hospital stay was significantly longer. 33 P. aeruginosa strains exhibit a complex interplay between resistance mechanisms, both intrinsic and acquired. 34 The current gaps in understanding of some of these mechanisms were reflected in the variable concordance between phenotypic and genotypic resistance for the different antibiotics, even for those antibiotics belonging to the same class (aminoglycosides). Our characterization of the carbapenem resistance showed a combination of diverse known mechanisms, from inhibition of antibiotic influx into the cell, to upregulation of antibiotic efflux out of the cell and carbapenemase enzymes. The concordance between phenotypic and genotypic predictions of AMR was high for the carbapenems, but it required a degree of curation of results that is not practical within public health settings.
There are clear limitations in the genomic predictions of AMR for P. aeruginosa. First, publicly available, curated databases are not comprehensive of all the known mechanisms. We found no mutations leading to upregulation of the chromosomal cephalosporinase AmpC (bla PAO ), but an exhaustive search would require additional analyses. Second, the regulatory pathways of some mechanisms are not fully understood, such as those that regulate AmpC. 34, 35 Third, extensive manual curation of some of the predictions is needed to ensure accuracy, for example of the loss-of-function mutations in the oprD gene.
The most prevalent clone in our data set was ST235 (27.8% of the isolates, n = 49), found throughout the Philippines. ST235 is a well-characterized international epidemic clone causing drug-resistant nosocomial outbreaks. 32 Isolates carrying bla VIM-2 and belonging to ST235 were reported from Malaysia, the Republic of Korea and Thailand. 13 Using WGS, we showed evidence of potential localized hospital outbreaks of ST235, as well as of persistence or reintroduction of this clone within one hospital. The number of SNP differences between genomes of isolates from different patients (0-24) were consistent with those reported for a persistent outbreak of P. aeruginosa in a hospital in the United Kingdom of Great Britain and Northern Ireland. 36 We also showed evidence of transfer of ST235 between hospitals, with isolates from different hospitals separated by as few as seven SNPs. Patient transfer between hospitals is not common in the Philippines, but the sampling for this   36 We also showed evidence of transfer of ST235 between hospitals (Fig. 2C), with isolates from different hospitals separated by as few as seven SNPs. Patient transfer between hospitals is not common in the Philippines, but the sampling for this study only allows us to hypothesize about a possible role of the community, animals or the environment in the spread of this clone.
It was previously proposed that ST235 emerged in Europe in about 1984, coinciding with the introduction of fluoroquinolones, and then disseminated to other regions via two independent lineages, acquiring resistance determinants to aminoglycosides and β-lactams locally. 14 Simultaneous acquisition of resistance to multiple antibiotics via integrons, transposons and integrative conjugative elements is well described in P. aeruginosa, ( 36 ) and is apparent in the distribution of resistance genes in our genomes (Fig. 2C). We have shown an example of a class I integron carrying six resistance genes in the genetic background of ST235 (Fig. 2C-D). While this integron shared some features with others previously described in P. aeruginosa, 13,32 such as the 5' and 3' conserved curated databases are not comprehensive of all the known mechanisms. We found no mutations leading to upregulation of the chromosomal cephalosporinase AmpC (bla PAO ), but an exhaustive search would require additional analyses. Second, the regulatory pathways of some mechanisms are not fully understood, such as those that regulate AmpC. 34, 35 Third, extensive manual curation of some of the predictions is needed to ensure accuracy, for example of the loss-of-function mutations in the oprD gene.
The most prevalent clone in our data set was ST235 (27.8% of the isolates, n = 49), found throughout the Philippines. ST235 is a well characterized international epidemic clone causing drug-resistant nosocomial outbreaks. 32 Isolates carrying blaVIM-2 and belonging to ST235 were reported from South Korea, Malaysia and Thailand. 13 Using WGS, we showed evidence of potential localized hospital outbreaks of ST235, as well as of persistence or reintroduction of this clone within one hospital. The number of SNP differences between genomes of isolates from different patients (0-24) were consistent with those reported for a persistent outbreak  ( study only allows us to hypothesize about a possible role of the community, animals or the environment in the spread of this clone.
It was previously proposed that ST235 emerged in Europe around 1984, coinciding with the introduction of fluoroquinolones, and then disseminated to other regions via two independent lineages, acquiring resistance determinants to aminoglycosides and β-lactams locally. 14 Simultaneous acquisition of resistance to multiple antibiotics via integrons, transposons and integrative conjugative elements is well described in P. aeruginosa, 36 and is apparent in the distribution of resistance genes in our genomes. We have shown an example of a class 1 integron carrying six resistance genes in the genetic background of ST235. While this integron shared some features with others previously described in P. aeruginosa, 13,32 such as the 5' and 3' conserved segments, 37 the gene composition and synteny was different, supporting the hypothesis of local acquisition of resistance.
Country-specific ST235 lineages have been reported previously, 11,14 confirming that country-wide clonal expansions may occur in the context of the global circulation of this clone. A previous longitudinal study showed VIM-2-positive ST235 spreading throughout Belarus, Kazakhstan and the Russian Federation, albeit without the resolution of whole genome data. 38 The contextualization of our genomes with international ST235 genomes showed a distinct cluster of Philippine genomes with limited genetic variability, suggesting the clonal expansion and geographic dissemination of this lineage across the Philippines. Alternatively, this could be explained by the limited representation of the Western Pacific Region in the collection of global genomes, highlighting the need for public genome data with more even geographical coverage. Our retrospective survey contributed to bridging this gap by making raw sequence data available on public archives.
In conclusion, our detailed description of the epidemiology and resistance mechanisms of ST235 in the Philippines suggests that the burden of XDR P. aeruginosa infections in the Philippines may be largely driven by a local lineage of the international epidemic clone ST235. A recent study in a hospital in Jakarta, Indonesia analysed the population composition of P. aeruginosa before and after a multifaceted infection control intervention, with the relative abundance of ST235    textualization of our genomes with international ST235 genomes showed a distinct cluster of Philippine genomes ( Fig. 3B) with limited genetic variability (Fig. 2C), suggesting the clonal expansion and geographic dissemination of this lineage across the Philippines. Alternatively, this could be explained by the limited representation of the Western Pacific region in the collection of global genomes, highlighting the need for public genome data with more even geographical coverage. Our retrospective survey contributed to bridging this gap by making raw sequence data available on public archives. segments, 37 the gene composition and synteny was different, supporting the hypothesis of local acquisition of resistance.
Country-specific ST235 lineages have been reported previously, 11,14 confirming that country-wide clonal expansions may occur in the context of the global circulation of this clone. A previous longitudinal study showed VIM-2 positive ST235 spreading throughout Belarus, Kazakhstan and the Russian Federation, albeit without the resolution of whole genome data. 38 The con-  (1) NalC/D LOF (14) None (1)  almost halved in the 10 months post-intervention. 39 This highlights the importance of hospital infection control and of preventive measures to contain the spread of this high-risk clone.