Whole-genome sequencing as an improved means of investigating Neisseria gonorrhoeae treatment failures

Additional keywords: antimicrobial resistance, public health, surveillance.

Antimicrobial-resistant (AMR) Neisseria gonorrhoeae has become a major global public health concern. Ceftriaxone is considered the last remaining empirical monotherapy for gonorrhoea, but isolates exhibiting decreased susceptibility to ceftriaxone are widespread¹ and there are recent reports of transmissible ceftriaxone-resistant strains from several countries.^2–4 Although relatively rare, there have also been a handful of reports regarding treatment failure using ceftriaxone monotherapy.^5–10 A dual therapy approach, combining azithromycin with ceftriaxone, as recommended by the World Health Organization (WHO), has recently been adopted in numerous countries as a means of delaying the emergence and spread of resistant gonorrhoea.¹¹ Although there is some suggestion that dual therapy has stabilised or even decreased cephalosporin resistance,¹² there have been concomitant widespread increases in azithromycin-resistant N. gonorrhoeae, as well as reports of strains now exhibiting resistance to both ceftriaxone and azithromycin.^2,3 There have also been at least two confirmed treatment failures using the dual therapy approach.^3,13

Among the strategies proposed in the WHO global action plan to minimise the impact of AMR gonorrhoea is the need to develop protocols to standardise and verify gonorrhoea treatment failures.¹⁴ At the laboratory level, this typically involves using phenotypic or genotypic tools including multilocus or N. gonorrhoeae multi-antigen sequence typing (MLST and NG-MAST respectively) to evaluate pre- and post-treatment isolates that are taken approximately 1 week apart. Isolates that are indistinguishable based on these traditional methods are considered treatment failures, whereas distinguishable isolates are considered consistent with reinfection. The major limitation of traditional gonococcal typing approaches is that they can have limited discriminatory power and may fail to distinguish different strains. Conversely, many target highly variable regions, like the gonococcal outer membrane protein which is encoded by the porB gene that may potentially mutate during therapy.

The NG-MAST method is a widely used gonococcal genotyping tool that involves DNA sequencing of the gonococcal transferrin-binding protein (which is encoded by tbpB) and porB gene, corresponding with 390 and 490 bp fragments, respectively.¹⁵ The method has been used extensively to examine treatment failures.¹⁶ In 2011, when ceftriaxone monotherapy was the recommended treatment for uncomplicated gonorrhoea, NG-MAST was used to investigate two suspected treatment failures in Australia, both involving pharyngeal infection.⁸ Upon having their initial samples taken, the individuals were asked to abstain from any further sexual intercourse and return for a ‘test of cure’ approximately 1 week later to determine whether the infection had resolved or not. Although isolates from Case 1 were indistinguishable by NG-MAST and were consistent with clinical information indicating treatment failure, the isolates from Case 2 were distinguished by NG-MAST based on an 18-bp deletion in the porB sequence of the post-treatment isolate and a difference in the ceftriaxone minimum inhibitory concentration (MIC; albeit within one doubling dilution). The patient from Case 2 denied sexual contact in the follow-up period, raising questions over reliability of the NG-MAST results. Ultimately, we categorised Case 2 as a suspected treatment failure.

The issues associated with Case 2 highlight the potential limitations of established laboratory methods for investigating gonococcal treatment failures; notably such issues have been observed in earlier suspected treatment failure cases.⁷ As identified by the WHO, addressing these issues is a priority, particularly in an environment where ceftriaxone-resistant strains are now emerging and rapid identification of treatment failures would be needed to facilitate timely public health responses. Whole-genome sequencing (WGS) is now being widely used to enhance N. gonorrhoeae epidemiological investigations, including the ability to determine direct or indirect transmission links between patients¹⁷ and assessing minimal nucleotide substitutions between multiple isolates from the same patient.^17,18 A recent study by Harris et al.¹⁹ comprehensively showed the potential limitations of NG-MAST for examining N. gonorrhoeae molecular epidemiology, including the inability of NG-MAST to accurately associate isolates with clinically relevant phylogenetic clades attained via WGS analyses. By providing single-nucleotide resolution, WGS enables better discrimination between strains considered to be linked via NG-MAST.²⁰ In this study, we used Cases 1 and 2 (described above) as examples to examine the utility of WGS to assess treatment failure. Furthermore, we sought to address whether short-read data provided sufficient resolution to aid in interpretation of these cases.

This investigation assessed two pharyngeal isolates from each individual. The first isolate was collected before the administration of antibiotics (pre-treatment isolate) and the second after antibiotic therapy (post-treatment isolate). The pre-treatment isolates were collected upon patients’ initial visit in March and July of 2011 (Case 1 and 2 respectively). The post-treatment isolates were acquired after patients’ return visits 1 week later, at which time the pharynx was still N. gonorrhoeae positive in both cases.

All four pharyngeal isolates underwent short-read Illumina (San Diego, CA, USA) sequencing, and the two post-treatment isolates underwent additional long-read sequencing using Pacific Biosciences (PacBio; Menlo Park, CA, USA) to generate high-quality genomes for each case. Only one isolate from each case was selected for long-read sequencing because a closely related reference genome is sufficient for identifying variants among isolates in cases of treatment failure. It should be noted that the Illumina sequencing for the pre- and post-treatment isolates was conducted at different times during the course of this study, and that the isolates were subsequently sequenced at different institutions.

Briefly, for post-treatment isolates, original stock cultures were recovered from –70°C storage and plated onto LB agar, which was then incubated at 37°C with 5% CO₂ in air for 24 h. A single colony was selected and subcultured onto LB agar and incubated under the conditions described above. This culture was used for both PacBio and Illumina sequencing. Genomic (g) DNA was extracted using the Ultraclean Microbial DNA Isolate Kit (GeneWorks, Adelaide, SA, Australia) according to the manufacturer’s instructions. These post-treatment DNA extracts underwent sequencing using a PacBio RSII sequencer with P6-C4 chemistry by an external sequencing service provider (Doherty Institute, The University of Melbourne, Vic., Australia). For Illumina sequencing of the post-treatment isolates, libraries were prepared using the Nextera XT DNA kit (Illumina) and sequenced on the NextSeq 500 using the High Output v2 kit (Illumina) at the Australian Centre for Ecogenomics, University of Queensland.

The pre-treatment isolates were cultured as described above, and DNA was extracted using the DSP DNA Mini Kit on the QIAsymphony SP (Qiagen, Valencia, CA, USA). The pre-treatment isolates were sequenced at Queensland Health Forensic and Scientific Services (Coopers Plains, Qld, Australia), using the same Illumina method as above, with the exception that the Mid Output v2 kit (Illumina) was used.

The PacBio raw sequence data were de novo assembled using the hierarchical genome assembly process (HGAP 3.0,²¹ using SMRT Analysis v2.3.0 by PacBio; https://www.pacb.com/products-and-services/analytical-software/smrt-analysis/; accessed 25 August 2019) and polished using Quiver (https://github.com/PacificBiosciences/GenomicConsensus; accessed 25 August 2019),²¹ both using default settings. The assembled chromosome for each isolate was circularised by screening the respective 5′ and 3′ ends to identify overlapping sequences using Contiguity (https://mjsull.github.io/Contiguity/; accessed 25 August 2019), which were then manually trimmed. The circularised genomes underwent additional reiterative polishing with Pilon v1.22 (https://github.com/broadinstitute/pilon; accessed 25 August 2019)²² by using the Illumina reads to resolve single nucleotide insertions and deletions. Genome assemblies were further corrected by manual inspection of the read pileups. Complete chromosomes were then annotated using Prokka v1.12 (https://github.com/tseemann/prokka; accessed 25 August 2019).²³

The Illumina raw sequence reads were assessed for quality using FastQC v0.11.4 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/; accessed 26 August 2019) and were subsampled to 1 100 000 reads if necessary using seqtk (https://github.com/lh3/seqtk; accessed 26 August 2019) to achieve comparable genome coverage for all isolates. A minimum read length of 100 bases was kept, and bases with a quality score <10 were removed using Trimmomatic v0.36 (https://github.com/timflutre/trimmomatic; accessed 26 August 2019).²⁴ De novo assembly was achieved using SPAdes v3.10.1 (https://github.com/ablab/spades; accessed 26 August 2019)²⁵ with the ‘–careful’ flag and kmers 21, 33, 55, 77 and 99. Contigs displaying less than 10-fold coverage were removed. Typing and identification of resistance markers for each isolate were assessed in silico (see Methods section in Supplementary Material).

SHRiMP (http://compbio.cs.toronto.edu/shrimp/; accessed 26 August 2019) as implemented in Nesoni v0.132 (https://github.com/Victorian-Bioinformatics-Consortium/nesoni; accessed 26 August 2019) was used with default settings to determine single nucleotide variants (SNVs) between the pre- and post-treatment isolates for each case, using their respective complete genome as a reference. Core SNVs were defined using Nesoni’s nway function with the ‘–require-all’ flag. SNVs associated with erroneous mapping of reads to repetitive regions and those found in regions of low coverage (<10-fold coverage) were subsequently removed. The Artemis Comparison Tool (http://sanger-pathogens.github.io/Artemis/ACT/; accessed 26 August 2019)²⁶ was used to manually inspect read pileups and to detect structural variants.

For each case, short reads were also aligned using the same approach to both close and distantly related, publicly available N. gonorrhoeae reference genomes determined from a pairwise distance matrix (see Methods section in Supplementary Material). Genome coverage plots using reference genomes of varying genetic relatedness can be seen in Figure S1, available as Supplementary Material to this paper. A phylogeny comprising isolates from both cases and publicly available genomes was generated as described previously,¹⁸ and further details are given in the Results section of the Supplementary Material.

The number of core genes was determined using Roary v3.12 (https://sanger-pathogens.github.io/Roary/; accessed 26 August 2019)²⁷ for de novo assembled isolates from Case 2 and included either the close or distant publicly available reference genome (see the Results section in the Supplementary Material).

For Case 1, genome data have been deposited under BioProject PRJNA451380, which comprises the pre-treatment isolate TFG-B1 and post-treatment isolate TFG-B2 (GenBank accession CP032429). For Case 2, genome data have been deposited under BioProject PRJNA451379, which comprises pre-treatment isolate TFG-A1 and post-treatment isolate TFG-A2 (GenBank accession CP032398). A summary of methylation profiles for both TFG-B2 and TFG-A2 is given in Table S1, available as Supplementary Material to this paper.

This study was approved by the South Eastern Sydney Local Health District Human Research Ethics Committee.

The overall TFG-B2 PacBio post-treatment genome size was 2 224 596 bp, which comprised a 2 220 291 bp chromosome and a 4305 bp cryptic plasmid. Both pre- and post-treatment isolates were MLST 1901 and NG-MAST 225. They both comprised a non-mosaic penicillin-binding protein 2 Type XII sequence, which is encoded by the penA gene, a single nucleotide (A) deletion in the promoter of the multiple transferable resistance Regulator, which is encoded by the mtrR gene, wild-type 23S rRNA alleles and two previously characterised non-synonymous nucleotide substitutions in porB (Table 1). No SNVs were detected between the pre- and post-treatment isolates after excluding base calls at both repetitive and low coverage regions in the TFG-B2 reference genome (Table S2). However, a comparison of the pre- and post-treatment draft genomes revealed the deletion of a 2.3 kb region involving the gonococcal pilin subunit encoded by pilE in the pre-treatment isolate. Loss of pilE in the pre-treatment isolate was confirmed using read mapping. Notably, one region excluded from the SNV analysis (positions 2110232–2110249 in TFG-B2) was located upstream of pilE, and is known to comprise portions of a partial pilin gene copy including a conserved region known as cys2 (Figure S2).²⁸ The TFG-B1 SNVs identified at this region were well supported by read coverage (Table S2), suggesting this could be a genuine difference between isolates. However, highly conserved cys2 regions are also present immediately preceding the hypervariable tail encoded by the 3′ end of pilE and multiple silent copies encoded by pilS loci distributed in other parts of the genome, suggesting that mismapping accounts for these differences (Figure S3).

Table 1.  Summary of treatment failures cases
AZI, azithromycin; CEF, ceftriaxone; MLST, multilocus sequence type; mtrR, multiple transferable resistance Regulator; NG-MAST, Neisseria gonorrhoeae multiantigen sequence type; PBP2, penicillin-binding protein 2; porB, major outer membrane porin; WT, wild type

The overall size of the TFG-A2 PacBio post-treatment genome was 2 225 966 bp, which comprised a 2 221 661 bp chromosome and a 4305 bp cryptic plasmid. Both isolates were MLST 1901, but had different NG-MAST profiles. The post-treatment isolate was identified as NG-MAST 225, but the pre-treatment isolate exhibited a novel NG-MAST profile due to the 18-bp deletion in porB. The isolates shared the same resistance markers as reported above for Case 1 (summarised in Table 1). Owing to the pilE observation for Case 1, pilE was also specifically interrogated. A 2.7 kb region that included pilE was found to be missing from the pre-treatment isolate, indicating a deletion of this gene. No other SNVs were detected between the pre- and post-treatment isolates after excluding calls at repetitive and low-coverage regions in the TFG-A2 reference genome (Table S3). Notably, one of the excluded regions (positions 76865–76969) corresponded to a phase variable repeat tract associated with the opacity-associated protein encoded by opa.

For Case 2, we also developed real-time polymerase chain reaction (PCR) assays to confirm the 18-bp deletion was present in the pre-treatment isolate (Table S4) because this deletion was previously reported in the post-treatment isolate.⁸ Briefly, we used both high-resolution melting and allele-specific PCR assays (Table S4) to detect the deletion. These PCR assays were applied to the DNA prepared for sequencing along with the original suspensions used in the initial investigation of these cases. This PCR testing confirmed a mistake had been made in the original article, and that the 18-bp deletion was in the pre-treatment isolate for Case 2. As a result, we are preparing a correction to Sexual Health.

Because not all laboratories have the capacity to rapidly generate complete reference genomes on a case-by-case basis, we investigated whether short-read sequence data alone were sufficient for assessing treatment failures. Fifteen high-quality N. gonorrhoeae complete reference genomes were obtained from GenBank with two (32867 and FA1090; GenBank accessions CP016015.1 and AE004969.1 respectively) selected for further comparison based on their genetic distance to isolates from Case 1 and Case 2 (Figures S4, S5). Reference 32867 was found to be the most closely related strain to Case 1 and Case 2, differing by 312 and 315 single nucleotide polymorphisms (SNPs) respectively, and shares the same MLST and NG-MAST profile. Reference FA1090 (MLST 1899 and NG-MAST 773) was one of the most distantly related strains, differing from either case by over 4000 SNPs, and is from a different ancestral lineage.

Read mapping to reference 32867 identified minor variation between the pre- and post-treatment isolates from Case 1 and Case 2, differing by 1 and 21 SNVs respectively. When using the distantly related reference FA1090, the pre- and post-treatment isolates for Case 1 and Case 2 differed by 73 and 82 SNVs respectively. However, following the removal of SNVs associated with repetitive or low-coverage regions, pre- and post-treatment isolates from either case were found to be indistinguishable. It should be noted that variants called between pre- and post-treatment isolates from Case 2 include the 18-bp deletion in porB. The location of these SNVs for each case can be found in Tables S2 and S3 for Cases 1 and 2 respectively.

In the context of 15 complete N. gonorrhoeae reference genomes, the phylogeny revealed that the isolates in Case 1 were very closely related to isolates in Case 2, with only 28 nucleotide substitutions separating them (Figure S4). To further explore the relatedness between these cases, an additional phylogeny was constructed using a global collection of 66 publicly available N. gonorrhoeae clinical isolates (Table S5)^{17,19,29–31} with the same MLST (1901) and NG-MAST (225) profiles (Figure S6). The complete genome of N. gonorrhoeae 32867 was also included as a reference. The treatment failure cases cluster with 25 isolates (from several countries) in a distinct, highly supported subclade (86% bootstrap support) separated from the other isolates by a single defining nucleotide substitution. This was a non-synonymous substitution (K94N) in the EamA-like family transporter protein AS012_00330 (as annotated in reference genome 32867, base position 62729 of the complete sequence).

In this study we compared pre- and post-treatment isolates from two previous cases (one confirmed and one suspected) of N. gonorrhoeae ceftriaxone treatment failure as examples to assess the feasibility and utility of WGS. Given the limited number of reports of N. gonorrhoeae ceftriaxone treatment failure for which pre- and post-treatment isolates have been available (only 13 cases documented to date),^3,5–10,13 the isolates from these two cases provided a rather rare opportunity to assess the value of using WGS compared with NG-MAST. We also emphasise that having appropriate laboratory methods to help investigate potential cases of ceftriaxone treatment failure are of increasing importance, particularly as ceftriaxone resistance is continuing to emerge.⁴

For Case 1, for which both the previous clinical and NG-MAST data were consistent with treatment failure, no differences were detected between the two isolates using WGS, aside from a pilE deletion observed only in the pre-treatment isolate. This is consistent with previous findings elsewhere showing that multiple isolates taken from the same patient,^17,18 or known sexual contacts,³² have very limited core differences. For Case 2, WGS confirmed the 18-bp deletion in porB (albeit in the pre- rather than post-treatment isolate; a correction is being prepared for the relevant journal regarding this mistake), but otherwise no other genomic differences were noted aside from a pilE deletion similar to that observed in Case 1. In the specific context of trying to confirm or exclude treatment failures, Case 2 highlights the benefits of using the enhanced resolution of WGS over NG-MAST by being able to compare the whole genome of two isolates, rather than two highly variable regions that undergo frequent antigenic variation. Of note is that there have only been five likely ceftriaxone treatment failure cases reported from Australia,^5,7,8 all in the past 10 years, and in two of these cases (Case 2 herein and the case reported by Chen et al.⁷), changes in porB hindered NG-MAST interpretation. In both cases, additional DNA sequencing and phenotypic AMR results were compared to help interpret the NG-MAST data. We contend that such ad hoc investigations are far from ideal and that WGS offers a new opportunity to standardise laboratory protocols used for verifying gonorrhoea treatment failure.

Although WGS provides enhanced resolution, it is important to map read data to an appropriate reference genome to provide the most accurate comparison.³³ Coverage plots can help assess regions where there is a significant drop in coverage, which are indicative of unreliable areas to call SNVs and often include repetitive regions. Selecting a distantly related genome will likely result in more false-positive SNVs and can reduce the likelihood of detecting true differences as the core genome becomes smaller (see Results section of the Supplementary Material). To generate an accurate SNV profile it is essential to filter out these false-positives,³³ which are often the result of mismapping within repetitive regions. This can be achieved by masking these repetitive regions within the reference genome. Furthermore, there may be near-identical repetitive regions among the genomes of interest that are present in fewer copies in the reference genome, which will also cause false-positive SNVs. This was observed when using the distantly related strain (FA1090), indicating that the use of a closely related reference may reduce the likelihood of this issue arising. If further investigations of repeat regions are needed, either long-read sequencing or targeted resequencing can help explore these regions. However, even with the advantage of a closely related reference in this study, we still needed to exclude certain regions from analysis. For example, for Case 2 we excluded a region encompassing the phase variable repeat tract of the opa gene, for which there are several highly similar copies in the N. gonorrhoeae genome that are known to independently phase vary on or off depending on repeat tract length.³⁴ Similarly, for Case 1, regions residing near pilE were excluded (positions 2110232–2110249). PilE contributes to pathogenesis via adhesion and is known to undergo high rates of antigenic variation due to recombination of multiple pilS loci into pilE, which comprises both hypervariable and highly conserved regions.³⁵

Independent deletions of pilE in the pre-treatment isolates for both cases is also intriguing, because may it suggest a link between pilE and treatment failure. However, there is limited evidence of clinical isolates lacking pilE.³⁶ Although high rates of antigenic variation are often associated with immune evasion, the loss of pilE could be a favourable adaptation because it may not be necessary to use pili once inside the host. Further studies are needed to explore these questions, which are outside the scope of the present study. Although PilE is involved in initial attachment to host cells, the capacity of pilE deletion mutants to cause urethral infection has been demonstrated previously.³⁷ Swanson et al.³⁸ have previously shown that both piliated and non-piliated colonies can be cultured from a single isolate. This could explain the pilE discrepancy for the cases of treatment failure, whereby non-piliated colonies were selected for when initially culturing these isolates. However, it is not possible to discern whether this loss occurred during infection or in vitro.

Overall, short-read sequencing was appropriate to help in the interpretation of these treatment failure cases. A critical advantage for sequencing treatment failures is that it allows us to rapidly assess relatedness to other cases. Although Cases 1 and 2 were identified within 5 months of each other in 2011 and shared a similar genotype and phenotypic profile, the original article treated them as two distinct events.⁸ However, the present study shows that the overall diversity between isolates in Cases 1 and 2 was minimal, clustering together in a distinct subclade with an additional 25 globally diverse MLST 1901 lineage isolates, within the broader MLST 1901 phylogeny. Surprisingly, this subclade was defined by a single non-synonymous substitution in the EamA-like family transporter protein AS012_00330. AS012_00330 is highly conserved within N. gonorrhoeae but its function is unknown. In Escherichia coli K12, a homologue of EamA is encoded by ydgE, which encodes for an efflux pump associated with the export of different metabolites of the cysteine pathway.³⁹ AS012_00330 bears little resemblance to EamA from K12 (BLASTP 38% amino acid sequence identity, 14% sequence coverage). However, AS012_00330 does possess three transmembrane domains, suggesting it may be involved with the transport of specific substances across the membrane. Although it is tempting to speculate that the observed non-synonymous substitution in AS012_00330 may contribute to ceftriaxone treatment failures, there is no indication in the literature to suggest that the other 25 isolates that shared this mutation were, in fact, associated with treatment failure. Further work is warranted to determine the function of AS012_00330 and its contribution, if any, to N. gonorrhoeae ceftriaxone treatment failure.

NG-MAST 225 is prevalent here in Australia⁴⁰ and elsewhere,^41,42 with no further treatment failure cases having been reported with this sequence type. However, it is important to remain aware of potential cases, particularly those now failing dual therapy associated with azithromycin and ceftriaxone.^3,13 However, this has become increasingly challenging with nucleic acid amplification tests now commonplace for diagnosing gonorrhoea.⁴³ Therefore, only a minority of gonorrhoea infections have an appropriate culture (and associated phenotypic testing), suggesting there could be an underestimation of treatment failures, particularly in areas where culture data are either limited or non-existent.

Having now used WGS to confirm that the Case 2 isolates were indistinguishable (except for the porB and pilE deletions) raises further questions regarding Case 2 and the significance of porB. We are intrigued by the fact that the pre-treatment isolate harboured the 18-bp porB deletion and that this isolate was of a novel and presumably rare NG-MAST type. It is therefore somewhat counterintuitive that the pre-treatment isolate could then mutate into NG-MAST 225, which is a widely reported NG-MAST profile.^40,44 Although the patient could have been reinfected before follow-up, they denied sexual intercourse during that period of time. While we still rely on individuals being truthful about when they engage in sexual intercourse, WGS provides enhanced resolution to aid in the investigation of treatment failures by being able to assess numerous genomic characteristics, including SNVs, structural variants and recombinant regions. This resolution is not achieved with traditional typing tools. We hypothesise that this patient may, in fact, have harboured a mixture of the above isolates (with and without the 18-bp porB deletion, and possibly the pilE deletion) before treatment, and that following treatment NG-MAST 225 was selected for. Unfortunately, we do not have the original swab samples available to test this hypothesis. Furthermore, although we only selected a single colony to culture, De Silva et al.¹⁷ have shown there is very minimal within-host diversity when independently sequencing a single anatomical site.

The 18-bp deletion is also intriguing in the context that the PorB protein is an essential outer membrane porin for maintaining N. gonorrhoeae viability that is subject to immune pressure and plays a key role in antimicrobial resistance;⁴⁵ it is these factors that influence its high variability, particularly in the extracellular loops of the protein.⁴⁶ Interestingly, the observed 18-bp deletion sits within extracellular loop 5, which likely (although not experimentally confirmed for this precise deletion) affects the binding of complement regulatory proteins, thereby conferring resistance to complement-mediated killing.⁴⁷ Future experimental investigation is warranted to determine whether this precise deletion has any effect on ceftriaxone MIC. The changes in porB loop 3, involving amino acids G120 and A121, are well documented, and only appear to increase resistance in the presence of mtrR mutations.^45,48 Although known resistance determinants affecting ceftriaxone were reported (porB and mtrR) for both cases, the only site to remain positive was the pharynx. This may suggest the primary reason for these ceftriaxone treatment failures is associated with suboptimal drug penetration into the oropharyngeal tissue.⁴⁹

The WHO has indicated the need to develop standardised protocols for verifying N. gonorrhoeae treatment failure. Here, we propose that assessing the whole genome compared with two highly variable genes will provide both a more confident predictor for determining treatment failures and evidence of the immediate effect that WGS can have in an important and direct clinical application for N. gonorrhoeae. Although short-read sequencing is suitable for assessing these cases, it is important to use a closely related reference genome to reduce the number of false-positive SNVs that need to be excluded and provides a more accurate comparison. Moving forward, sequencing these treatment failures can provide a wealth of genomic information, including key resistance genes, and facilitate rapid comparison with similar cases elsewhere. Where available, we recommend using WGS for all suspected N. gonorrhoeae treatment failures.

The authors declare no conflicts of interest.