Genomic surveillance of Neisseria gonorrhoeae in the Philippines, 2013–2014

Antimicrobial-resistant Neisseria gonorrhoeae is a major threat to public health and is of particular concern in the Western Pacific Region, where the incidence of gonorrhoea is high. The Antimicrobial Resistance Surveillance Program (ARSP) has been capturing information on resistant gonorrhoea since 1996, but genomic epidemiology studies on this pathogen are lacking in the Philippines. We sequenced the whole genomes of 21 N. gonorrhoeae isolates collected in 2013–2014 by ARSP. The multilocus sequence type, multiantigen sequence type, presence of determinants of antimicrobial resistance and relatedness among the isolates were all derived from the sequence data. The concordance between phenotypic and genotypic resistance was also determined. Ten of 21 isolates were resistant to penicillin, ciprofloxacin and tetracycline, due mainly to the presence of the blaTEM gene, the S91F mutation in the gyrA gene and the tetM gene, respectively. None of the isolates was resistant to ceftriaxone or cefixime. The concordance between phenotypic and genotypic resistance was 92.38% overall for five antibiotics in four classes. Despite the small number of isolates studied, they were genetically diverse, as shown by the sequence types, the N. gonorrhoeae multiantigen sequence typing types and the tree. Comparison with global genomes placed the Philippine genomes within global lineage A and led to the identification of an international transmission route. This first genomic survey of N. gonorrhoeae isolates collected by ARSP will be used to contextualize prospective surveillance. It highlights the importance of genomic surveillance in the Western Pacific and other endemic regions for understanding the spread of drug-resistant gonorrhoea worldwide.

The increase in N. gonorrhoeae infections resistant to front-line antibiotics triggered a global action plan from WHO to control the spread and impact of gonococcal resistance and a call for international collaborative action, especially in the Western Pacific Region. 1 The WHO Antimicrobial-resistant Neisseria gonorrhoeae is a major threat to public health and is of particular concern in the Western Pacific Region, where the incidence of gonorrhoea is high. The Antimicrobial Resistance Surveillance Program (ARSP) has been capturing information on resistant gonorrhoea since 1996, but genomic epidemiology studies on this pathogen are lacking in the Philippines.
We sequenced the whole genomes of 21 N. gonorrhoeae isolates collected in 2013-2014 by ARSP. The multilocus sequence type, multiantigen sequence type, presence of determinants of antimicrobial resistance and relatedness among the isolates were all derived from the sequence data. The concordance between phenotypic and genotypic resistance was also determined.
Ten of 21 isolates were resistant to penicillin, ciprofloxacin and tetracycline, due mainly to the presence of the bla TEM gene, the S91F mutation in the gyrA gene and the tetM gene, respectively. None of the isolates was resistant to ceftriaxone or cefixime. The concordance between phenotypic and genotypic resistance was 92.38% overall for five antibiotics in four classes. Despite the small number of isolates studied, they were genetically diverse, as shown by the sequence types, the N. gonorrhoeae multiantigen sequence typing types and the tree. Comparison with global genomes placed the Philippine genomes within global lineage A and led to the identification of an international transmission route.
This first genomic survey of N. gonorrhoeae isolates collected by ARSP will be used to contextualize prospective surveillance. It highlights the importance of genomic surveillance in the Western Pacific and other endemic regions for understanding the spread of drug-resistant gonorrhoea worldwide.
Antimicrobial-resistant Neisseria gonorrhoeae is a major threat to public health and is of particular concern in the Western Pacific Region, where the incidence of gonorrhoea is high. The Antimicrobial Resistance Surveillance Program (ARSP) has been capturing information on resistant gonorrhoea since 1996, but genomic epidemiology studies on this pathogen are lacking in the Philippines.
We sequenced the whole genomes of 21 N. gonorrhoeae isolates collected in 2013-2014 by ARSP. The multilocus sequence type, multiantigen sequence type, presence of determinants of antimicrobial resistance and relatedness among the isolates were all derived from the sequence data. The concordance between phenotypic and genotypic resistance was also determined.
Ten of 21 isolates were resistant to penicillin, ciprofloxacin and tetracycline, due mainly to the presence of the bla TEM gene, the S91F mutation in the gyrA gene and the tetM gene, respectively. None of the isolates was resistant to ceftriaxone or cefixime. The concordance between phenotypic and genotypic resistance was 92.38% overall for five antibiotics in four classes. Despite the small number of isolates studied, they were genetically diverse, as shown by the sequence types, the N. gonorrhoeae multiantigen sequence typing types and the tree. Comparison with global genomes placed the Philippine genomes within global lineage A and led to the identification of an international transmission route.
This first genomic survey of N. gonorrhoeae isolates collected by ARSP will be used to contextualize prospective surveillance. It highlights the importance of genomic surveillance in the Western Pacific and other endemic regions for understanding the spread of drug-resistant gonorrhoea worldwide.
Villamin et al N. gonorrhoeae surveillance in the Philippines antimicrobials representing four different classes, namely penicillin (PEN), ciprofloxacin (CIP), tetracycline (TCY), ceftriaxone (CRO) and cefixime (CFX), by disc diffusion and gradient diffusion (Etest, BioMerieux). The minimum inhibitory concentrations were interpreted as resistant, intermediate or susceptible according to the interpretative criteria in the Performance Standards for Antimicrobial Susceptibility Testing (26th edition) of the Clinical and Laboratory Standards Institute (CLSI). 12 All isolates were screened for β-lactamase production on cefinase paper discs (BD BBL). Azithromycin was not included in the panel of antibiotics, because, at the time isolates were collected, it was not in the treatment guidelines of the Philippines and no CLSI breakpoint was available.

DNA extraction and whole-genome sequencing
DNA was extracted from a single colony of each isolate with a Wizard® Genomic DNA Purification Kit (Promega), and the quantity and quality were determined with a Quantus Fluorometer (Promega) with picogreen and a NanoDrop 2000c spectrophotometer (Thermo Fisher Scientific). The DNA extracts were then shipped to Wellcome Trust Sanger Institute for sequencing on the Illumina HiSeq platform with 100-bp paired-end reads. Raw sequence data were deposited in the European Nucleotide Archive under the study accession number PRJEB17615. Run accessions are provided on the Microreact projects.

Bioinformatics analysis
Genome quality was evaluated according to metrics generated from assemblies, annotation files and the alignment of the isolates to the reference genome of N. gonorrhoeae strain TCDC-NG08107 (accession CP002441.1), as previously described. 13 Annotated assemblies were produced with the pipeline previously described. 14 We included 21 high-quality N. gonorrhoeae genomes in this study.
The MLST sequence types (STs) and NG-MAST types, as well as the presence of AMR determinants (known genes or mutations) and clustering of the isolates according to genetic similarity (tree), were predicted in silico from genome assemblies with Pathogenwatch (https://www.sanger.ac.uk/tool/pathogenwatch/). 15 In parallel, the evolutionary relations among the isolates were inferred from single nucleotide polymorphisms Gonococcal Antimicrobial Surveillance Programme has operated in the Western Pacific and South-East Asian regions since 1992, but surveillance of gonococcal antimicrobial resistance (AMR) remains limited in the Asia-Pacific region. 6 In a recent report, 18 of 21 countries in the Asia and the Pacific reported isolates with decreased susceptibility to ceftriaxone and/or isolates resistant to azithromycin between 2011 and 2016. 6 The Antimicrobial Resistance Surveillance Program (ARSP) of the Philippines Department of Health has been contributing AMR surveillance data to the Western Pacific Gonococcal Antimicrobial Surveillance Programme since 1996 and did not confirm isolates with decreased susceptibility or resistance to these antibiotics during 2011-2016, 6 while high gonococcal resistance rates against other first-line antibiotics have long been reported ( Fig. 1). Continuous surveillance is thus key to detecting potential emergence or introduction of resistance to current treatment options.
Molecular methods for defining the epidemiology of gonococci include both N. gonorrhoeae multiantigen sequence typing (NG-MAST) 7 and multilocus sequence typing (MLST), 8 although NG-MAST is more widely used to investigate specific gonococcal AMR phenotypes. 8 Whole-genome sequencing (WGS) was recently shown to provide better resolution and accuracy than NG-MAST or MLST. 9 Good understanding of the population structure and the mechanisms of resistance of N. gonorrhoeae in the Philippines would allow detection of high-risk clones associated with high-risk groups and contribute to the clinical management of gonococcal-related diseases and the creation of policies to prevent the spread of drug resistance. 10, 11 Here, we describe the results of the first genomic survey of gonococcal isolates in the Philippines.

Bacterial isolates
A total of 51 N. gonorrhoeae isolates were collected at ARSP sentinel sites in 2013 and 2014 ( Table 1). Of the 36 isolates referred to the ARSP reference laboratory for confirmation, 22 isolates from seven sentinel sites were resuscitated and submitted for WGS.

Antimicrobial susceptibility testing
All N. gonorrhoeae isolates in this study were tested at the ARS reference laboratory for susceptibility to five N. gonorrhoeae surveillance in the Philippines Villamin et al (SNPs) by mapping the paired-end reads to the reference genome of N. gonorrhoeae strain FA 1090 (accession AE004969.1), as described in detail previously. 13 Mobile genetic elements described in the FA 1090 genome 16 were masked in the alignment of pseudogenomes with a script available at https://github.com/sanger-pathogens/ remove_blocks_from_aln. SNPs were extracted with snp-sites, 17 and a maximum likelihood phylogenetic tree was generated from 7518 SNP positions with RAxML, 18 the generalized time-reversible model and the GAMMA method of correction for among-site rate variation, with 500 bootstrap replicates.
To complement the Pathogenwatch AMR results, known AMR determinants were identified from raw sequence reads with ARIBA 19 and a curated database of known resistance genes and mutations available at https://github.com/martinghunt/ariba-publication/tree/ master/N_gonorrhoeae/Ref. The combined genotypic predictions of AMR (test) were compared with the phenotypic results (reference), and the concordance between the two methods was computed for each of five antibiotics (105 total comparisons). Isolates found to be resistant or to have reduced susceptibility (intermediate) were pooled as non-susceptible for comparison purposes. 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  Villamin et al N. gonorrhoeae surveillance in the Philippines mtrR promoter and the non-synonymous substitutions in ponA (L421P) and porB (A121D), which may contribute to high-level penicillin resistance. 22 The partially assembled bla TEM gene was also detected in this genome but was not considered in the concordance analysis. We also identified the non-synonymous mutation in the folP gene (R228S), which confers resistance to sulfonamides in 20 of 21 genomes; however, isolates are not routinely tested for resistance to this antibiotic class.
The concordance between antimicrobial susceptibility testing results and genotypic predictions ( Table 3) was >95% for all antibiotics except tetracycline (66.67%), resulting in an overall concordance of 92.38% ( Table 3). The discrepancies were attributed to seven false-positive results (isolates with a susceptible phenotype but with known AMR determinants in their genomes), all of which contained the mutation V57M in the rpsJ gene; two isolates also carried the tetM gene. 23 over the total number of isolates assessed (expressed as a percentage).
The maximum likelihood tree, genotyping results and AMR predictions and the metadata collected from the sentinel sites were visualized with Microreact. 20 To contextualize the genomes from this study with publicly available global genomes, we combined them with two surveys available on Pathogenwatch, a European survey of 1054 genomes 10 and a global survey of 395 genomes. 21

Ethics statement
Ethical approval is not applicable, as we used archived bacterial samples processed by ARSP. No identifiable data were used in this study.

Demographic and clinical characteristics of the N. gonorrhoeae isolates
The 21 genomes included in this study represented seven sentinel sites, with Vicente Sotto Memorial Medical Center (VSM) contributing the most isolates (n = 11). The highest incidence was in the age group 15-24 years (47.6%, n = 10), followed by the age groups 5-14 years (28.6%, n = 6), 25-34 years (14.3%, n = 3) and 1-4 years (9.5%, n = 2, Table 2). The numbers of isolates from females and males were almost equal (n = 11 and n = 10, respectively). The most frequent specimen source was the vagina (n = 8) for female patients and penile discharge (n = 5) for males. All the patients were outpatients (n = 21).

Concordance between phenotypic and genotypic AMR
Isolates were tested for susceptibility to five antibiotics representing four classes. In line with the resistance trends shown in Fig. 1, the most prevalent resistance profile was PEN CIP TCY, identified in 10 isolates from five sentinel sites and linked mainly to the presence of the bla TEM gene, the S91F mutation in GyrA and the tetM gene, respectively (Table 3). One penicillin-resistant isolate (13ARS_DMC0024) harboured three mutations, the -57delA mutation in the  contained isolates resistant to ciprofloxacin harbouring the gyrA_S91F mutation, clade II was characterized by the presence of gyrA_D95G alone or in combination with the one mutation in parC, while clade III was characterized by the presence of gyrA_D95A with one or two mutations in parC in all but one isolates. Clades II and III showed different geographical distributions, although both were present at the VSM and ZPH sentinel sites, which submitted the most isolates. The genome from sentinel site Southern Philippines Medical Center (DMC) with a deletion in the mtrR promoter (-57delA) was found on a separate branch (I) in the tree and also carried a different complement of AMR determinants, indicating that it is genetically distinct from the others (Fig. 2).

N. gonorrhoeae from the Philippines in the global context
The Philippine genomes were contextualized with two recently published collections 10,21 available in Pathogenwatch. A recent global collection of 395 genomes from 58 countries (including the Philippines) showed two major lineages with different evolutionary strategies. Most of the genomes in this study were found within a subclade of lineage A, a multidrug-resistant lineage associated with infection in high-risk sexual networks (Fig. 3A), and mixed with genomes from Europe, Pakistan and South-East Asia, including genomes previously isolated in the Philippines

Genotypic findings
In silico genotyping A total of 15 different STs were identified, with only four (9364, 10316, 1582 and 8133) represented by more than one isolate. Nine genomes were assigned to eight known NG-MAST types and 12 genomes to nine novel types. Only three sentinel sites, namely Corazon Locsin Montelibano Memorial Regional Hospital (MMH), VSM and Zamboanga Del Norte Medical Center (ZPH), submitted more than one isolate, and all were represented by almost as many STs as isolates submitted. Details of the numbers and the most common STs and NG-MAST types found at each sentinel site are shown in Table 4.

Population structure of N. gonorrhoeae in the Philippines
The diverse gonococcal population was represented by a tree with three deep-branching clades and no clear geographical signal (Fig. 2). Clades II and III were characterized by a different repertoire of AMR genes and mutations. Clade II contained mostly isolates susceptible to or with reduced susceptibility to tetracycline and with the V57M mutation in the rpsJ gene alone, while clade III was composed of isolates resistant to tetracycline and also containing the tetM gene. Similarly, while both clades Villamin et al N. gonorrhoeae surveillance in the Philippines Table 4. Distribution of isolates, STs, NG-MAST types, resistance profiles and resistance genes and mutations at the seven sentinel sites. The genetic resistance mechanisms for all the isolates from each site are listed.    two genomes is consistent with isolates from the same location, 10 suggesting an epidemiological link and a route of international transmission. The gonococcal reference laboratory at the Norwegian Institute of Public Health confirmed that the isolate was also resistant to penicillin and ciprofloxacin and susceptible to extended-spectrum cephalosporins (tetracycline not tested) and that the Norwegian male had visited the Philippines and claimed to have contracted gonorrhoea during his stay (personal communication).

DISCUSSION
WGS showed that the N. gonorrhoeae genomes from the Philippines are genetically diverse and carry a variety of AMR determinants, such as chromosomal mutations and acquired genes. The concordance between phenotypic (1998 and 2008 21 ). The genome from DMC was found within a separate subclade with genomes from Europe, India and Pakistan, which also carried the mtrR_-57delA promoter deletion.
We further contextualized our isolates with 1054 genomes from 20 countries collected in 2013 in a European survey. 10 Notably, the genome of strain 14ARS_VSM0347 isolated from a female on 25/04/2014 in the Philippines was highly similar to that of ECDC_GC_088 isolated from a male on 31/12/2013 in Norway (Fig. 3B), with no SNP differences found in the Pathogenwatch core genome (1542 genes). One SNP was identified in the referencebased alignment of the two pseudogenomes, confirming that these two genomes are highly similar. The two strains also shared the same complement of AMR genes and mutations (Fig. 3B) and genotypic resistance was high (>95%) for most antibiotics but only 66% for tetracycline. Susceptible isolates carrying only the rpsJ_V57M mutation, which confers low-level tetracycline resistance, have been reported previously. 24 The two susceptible isolates with a full-length tetM gene reported by Pathogenwatch were re-tested by disc diffusion and their susceptibility confirmed. The discrepancy, which could be explained for example by the lack of expression of the gene, errors in susceptibility testing or low-level DNA contamination, will be further investigated. Although N. gonorrhoeae remains largely susceptible to extended-spectrum cephalosporins and azithromycin in the Philippines, we identified one isolate (13ARS_DMC0024) with a -57delA mutation in the mtrR promoter, which results in overexpression of the MtrCDE efflux pump and, in combination with other mutations, can increase the minimum inhibitory concentrations of azithromycin (as well as penicillin and tetracycline). 4, 22 Prospective surveillance with WGS can detect in real-time the acquisition of additional mutations that could result in decreased susceptibility.
A recent report suggested rapid recent intercontinental transmission of gonorrhoea, with common introductions from Asia to the rest of the world. 21 In support of this finding, the Philippine genomes, most of which were within a subclade of global lineage A, were interspersed with those from other countries (Fig. 3A). In addition, we found evidence of an introduction event from the Philippines to Norway associated with travel of a Norwegian male to the Philippines (Fig. 3B).
Global lineage A is associated with infection in highrisk sexual networks. 21 The combination of WGS with epidemiological information can reveal transmission routes and risk factors, which can be used to design better control measures. 25 The small size of our retrospective data set (n = 21) and the linked epidemiological data did not permit any inferences about sexual networks or risk factors, which is a limitation of this study. The number of reported isolates by ARSP has, however, since increased to more than 100 per year, and data on risk factors are also collected, which will allow a more comprehensive analysis of the population diversity and of risk factors in future reports.
Our results represent the first genomic survey of N. gonorrhoeae isolates collected by ARSP and will constitute the background for contextualizing continuous prospective surveillance. In addition, it highlights the im-