The antimicrobial resistance landscape of Neisseria gonorrhoeae in New Zealand from November 2018 to March 2019 and the role of sexual orientation in transmission.

Data Summary

Sequencing reads have been deposited in the National Centre for Biotechnology Information (NCBI) database under BioProject number PRJNA738299. Supporting data, code and protocols have been provided within the Methods or as Supplementary Material. This study is based on a client report prepared by the same authors for the Ministry of Health New Zealand: https://surv.esr.cri.nz/PDF_surveillance/Antimicrobial/Gono/NgonoSurvey2019_FINAL.pdf (Ministry of Health Report Number FW21002). Supplementary Material can be found on Figshare: https://doi.org/10.6084/m9.figshare.16458153.v1.

Antimicrobial resistance rates (AMR) among pathogens of human concern are constantly on the rise and a major threat to public health. is of particular concern with current treatment guidelines threatened by rising rates of azithromycin resistance. Our study comprises 314 . isolates collected from across New Zealand over 4 months in 2018/2019. Combining phenotypic and genotypic approaches, we identify the current status of antimicrobial resistance rates in New Zealand, highlighting that azithromycin resistance rates remain low compared to other countries. We also use SNP analysis and sexual orientation data to identify clusters of transmission, confirming other studies in finding that clusters of mixed sexual orientation are important factors to consider in driving disease transmission. An analysis of pairwise genetic distance was further used to show genomic coherency within sequence types (STs) for each of the three typing schemes [multi-locus sequencing typing (MLST), multi-antigen sequence typing (NG-MAST), and sequence typing for antimicrobial resistance (NG-STAR)]. Our results confirm the superiority of the NG-STAR scheme in terms of determining evolutionary history. These analyses together provide a comprehensive snapshot of the population structure, antimicrobial resistance rates and diversity of circulating in New Zealand during that time period, allowing careful consideration of the appropriateness of the public health response.

Introduction

Rising antimicrobial resistance rates are a global problem for many pathogens of public health concern threatening effective treatment, management and control of bacterial infections [1]. , causative agent of gonorrhoea, is the second most common sexually transmissible disease. infection rates are on the rise world-wide, in the USA alone rates of reported gonorrhoea cases have increased 82.6% since 2009, with an estimated 86.9 million world-wide cases in 2016 [2-4]. In New Zealand, the national rate has increased from 75 cases per 100,000 population in 2012 to 123 per 100,000 population in the 12 months ending in June 2019, an increase of 64.0% (https://www.esr.cri.nz/our-services/consultancy/public-health/sti/, accessed February 23, 2021).

Asymptomatic infections are important sources of transmission, which most often occur as vaginal infections in women, but also for rectal, pharyngeal, and some urethral infections in men [5]. These unresolved infections can lead to impaired control efforts, as well as severe sequelae such as infertility, eye infections of neonates or enhanced transmission of HIV [6-8]. This makes men, who have sex with men (MSM), and sex workers particularly high-risk groups, as asymptomatic and untreated infections can lead to increased incidence rates [9, 10]. A first link between population risk groups was described recently, where men who have sex with men and women (MSMW) have been identified as a possible bridging population between heterosexual and MSM groups using hierarchical single-linkage clustering of pairwise SNP differences between isolates [11]. Such insights are important for appropriate implementation of prevention and control methods for the rapid spread of gonorrhoea.

has evolved resistance to all discovered and recommended first-line antibiotics used to treat it, from sulfonamides in the 1940s, to high-level penicillin resistance mediated by β-lactamase plasmids and chromosomal mutations in the 1980s and treatment failures with currently recommended antibiotic classes macrolides and cephalosporins [7, 12–14]. For the past decade, the recommended guidelines suggest a dual antimicrobial therapy of azithromycin and extended-spectrum cephalosporins for uncomplicated gonorrhoea cases (ESC, i.e. cefixime or ceftriaxone) [15-17]. While treatment recommendations in New Zealand have remained the same, the USA and UK have adapted their guidelines to monotherapy using cephalosporins due to the increasing rates of decreased susceptibility to cephalosporins and resistance to azithromycin [18, 19].

Continued surveillance is key to track the emergence and spread of AMR rates in gonococci. However, diagnostic laboratories world-wide have adopted more cost-effective and faster culture independent testing approaches, such as nucleic acid amplification testing (NAAT) targeting the gonococcal porA pseudogene or the multicopy opa genes being now the standard of care [20]. The limited availability of cultures to do enhanced surveillance has led to a gap in antimicrobial susceptibility testing and cultures available to perform WGS. Alternative methods, such as multivariate regression modelling to predict minimum inhibitory concentrations (MICs) from genetic determinants of antibiotic resistance [21] or the development of primer-panels targeting typing and antimicrobial resistance genes for are underway (Straub et al., unpublished), but until then we rely on regular culture-based surveillance studies.

Our aim was to [1] determine AMR rates of in New Zealand using WGS and antibiotic susceptibility testing [2], gain insights into the population structure [3], identify clusters of transmission with sexual orientation data and [4] identify temporal variation of AMR rates in New Zealand by comparing findings with a national survey performed in 2014/2015 [22, 23].

Methods

Isolate collection and patient information

Between 5 November 2018 and 10 March 2019, diagnostic laboratories were requested to refer cultured isolates to the Institute of Environmental Science and Research (ESR) for a national survey. A total of 315 isolates were submitted and epidemiological data such as the national health index (NHI) number, sex, date of birth and the site of the specimen were supplied with each isolate. District Health Boards (DHBs) are the regional entities providing public health services in New Zealand, with 20 representative DHBs. The patient’s domicile, and therefore DHB was inferred from the location of the referring laboratory unless the referring laboratory specified otherwise. Sexual risk group data (MSM, heterosexual males or females) was provided for 192 isolates.

Phenotypic susceptibility testing

Susceptibility to azithromycin, cefixime, ceftriaxone, ciprofloxacin, ertapenem, gentamicin, penicillin, spectinomycin and tetracycline were determined by agar dilution according to the methods of the Clinical and Laboratory Standards Institute (CLSI, 2012) [24]. Ertapenem was purchased from TokuE (Bellingham, WA, USA). All other antibiotic pure substances were purchased from Sigma-Aldrich (Saint Louis, MO, USA). β-lactamase production was determined using the chromogenic cephalosporin compound nitrocefin. MIC breakpoints were interpreted following European Committee on Antimicrobial Susceptibility Testing (EUCAST) guidelines [25] (Table S1). Decreased susceptibility to ceftriaxone (CTRX) and cefixime (CFIX) was defined as MIC ≥0.06 and ≥0.12 mg l−1, respectively, with the definition of the ceftriaxone MIC in line with the Australian Gonococcal Surveillance Programme (AGSP) [26] and the cefixime MIC based on definitions by Deng et al. [27].

DNA extraction and whole-genome sequencing

Bacterial cultures were grown on chocolate agar (Fort Richard Laboratories, New Zealand) for 18 h at 35 °C with 5% CO2. DNA extractions from single colony picks were performed using the High Pure PCR Template Preparation Kit (Roche Molecular Biochemicals, IN, USA) following the manufacturer’s instructions. Isolates were sequenced using 2×150 bp paired-end reads on the Illumina NextSeq platform at ESR using Illumina libraries (Nextera-XT DNA library preparation) and protocols.

Bioinformatics

Quality control and genome assembly

Raw sequencing reads were put through the nullarbor2 pipeline (https://github.com/tseemann/nullarbor). Species confirmation and identification of potential contaminants was performed using Centrifuge 1.0.4 [28] and the refseq microbial database (built October 2018). One sample was excluded due to the presence of contamination, hence a total of 314 isolates were retained for downstream analysis. SPAdes 3.13.1 [29] was chosen for de novo assembly and annotation was performed using Prokka 1.14.5 [30]. The international reference strain FA1090 (accession number NC_002946.2) was used as a reference for the analysis.

In silico molecular typing

The WGS draft genome assemblies were used for in silico determination of sequence types (STs) for seven housekeeping genes (abcZ, adk, aroE, fumC, gdh, pdhC, pgm) using mlst v2.17.6 (https://github.com/tseemann/mlst). Further, NG-MAST typing of por and tbpB was performed using ngmaster v0.5.5 (http://www.ng-mast.net accessed 4 June 2020) [31]. For alleles not present in the NG-MAST database, new allele numbers and/or STs were assigned, commencing with ‘21,000’. For antimicrobial resistance determinant sequence typing, the Sequence Typing for Antimicrobial Resistance (NG-STAR) STs [32] were identified using pyngSTar (https://github.com/leosanbu/pyngSTar, alleles and profiles from https://ngstar.canada.ca/ accessed 7 October 2020). PHYLOViZ 2.0 [33, 34] was used to compute a minimum spanning tree (using the goeBURST algorithm) to calculate relatedness of allelic profiles for the three typing schemes. Clonal complexes (CC) were defined by single locus variations in the allelic profile.

Characterization of AMR determinants

Assemblies were further screened for genetic determinants of antimicrobial resistance (plasmids, presence/absence of resistance genes and chromosomal point mutations) as listed in Table S2.

Plasmids

A screen for the presence of bla TEM and tetM genes was performed using abricate v.0.8.7 (https://github.com/tseemann/abricate) using the National Centre for Biotechnology Information (NCBI) AMR database (built February 2020). For determination of plasmid types, an in silico PCR (Ipcress, exonerate v2.2.0), was run using the primers from Turner et al. [35] to distinguish between the Dutch and American tetM plasmid type. Three types are known for the bla TEM plasmids (Asian, African, Toronto/Rio type) and the primers from Palmer et al. [36] were used to distinguish bla TEM types.

Mutations in 23S rRNA gene

The number of copies and presence of point mutations in the 23S rRNA gene was determined by mapping the raw reads to a reference with three masked copies of 23S, forcing mapping to a single copy of genome NCCP11945 using snippy v4.3.6 (https://github.com/tseemann/snippy). Copy numbers were determined using a custom python script (https://github.com/kwongj/ng23S-mutations).

Chromosomal mutations

Screening for chromosomal point mutations was performed using Ariba v. 2.14.4 [37] with a custom database created for (https://github.com/chrstraub/Ngono_Straubetal2021). An in silico PCR was performed to test for the presence of ermF, ermB, ermC, ereA, ereB and mefAB using existing primers [38, 39].

Variant calling and phylogenetic tree inference

Raw sequencing reads were trimmed to remove adaptor sequences and low-quality bases using Trimmomatic v.0.39 [40]. Trimmed reads were aligned to the reference FA1090 using bwa mem v0.7.17 [41], duplicate sequences were marked using Picard v2.23.6 [42]. SNPs were called using Freebayes v1.3.2 [43] under a haploid model with a minimum mapping quality of 30, minimum coverage of 10× and a minimum alternate allele fraction of 0.7. Further SNP filtering was performed using vcftools v0.1.16 [44]. Repeat regions were identified using MUMmer v3.23 [45] and subsequently masked. Prophages present in the reference FA1090 were also masked in the alignment [46]. This filtering approach resulted in a total of 17,437 core SNPs present in 100% of isolates.

Recombinant parts of the genome were identified with the core multi-sequence alignment provided as input and recombinant regions were removed using gubbins v2.3.4 [47], reaching tree convergence after ten iterations. The final core SNP alignment (.filtered_polymorphic_sites.fasta) was tidied up using snp-sites [48] to remove any other sites than ACGT (including Ns). This alignment was used as input for maximum-likelihood tree building using IQ-TREE v1.6.9 [49], after having identified the best fitting nucleotide substitution model using ModelFinder Plus [50]. IQ-TREE was run using 1,000 ultrafast bootstrap iterations [51].

A genomic dataset consisting of the 314 isolates from this study and the publicly available dataset of 401 isolates used in the New Zealand survey from 2014/2015 (BioProject number PRJNA394216) [23] was analysed using the same pipeline as described above. A core SNP alignment of 19,506 SNPs for the 716 isolates (including reference strain FA1090) were identified after filtering and prior to removal of recombination events identified with gubbins after 40 iterations.

Population structure

Genetic clusters were determined using the R package rhierbaps v1.1.3 implementing the hierBAPS algorithm [52]. Clustering was performed using two levels in the hierarchy and K=65 as the upper bound for the number of clusters.

Transmission analysis

Pairwise SNP-distances between isolates were calculated based on the recombinant-free core-genome alignment using the dist.gene function (R package ‘ape’ [53]). Hierarchical single-linkage clustering was performed using the hclust function and then filtered based on a SNP cut-off of 10 using function cutree (R package ‘stats’). Sexual risk group data was available for 192 isolates (61.1%) distinguishing between men who have sex with men (MSM), men who have sex with females (heterosexual) and females who have sex with males. Isolates for which sexual risk group data was not available were included in the clustering analysis and clusters are classified as being associated with a single sexual risk group even if isolates with unknown risk behaviour are grouped together.

Defining threshold for transmission clusters

No information about sex partners was collected for the study, hence no threshold determined by transmission events within sexual relationships could be used. A SNP threshold of 13 SNPs was shown for with-in host diversity of [54], as well as a more conservative SNP-threshold of ten SNPs (based on SNP differences between epidemiologically linked partners) was applied to determine transmission events for city and state-wide transmission networks in Australia and the USA [11, 55]. This conservative SNP-cut-off of ten was used in this study to determine possible transmission clusters.

Statistics

All statistical analyses were performed in R v 4.0.3 [56]. Testing for correlation of phylogenetic structure with binomial traits was performed using the package ‘caper’ [57], which implements Fritz and Purvis D [58].

Results

Epidemiological characteristics

During the survey collection period from November 2018 to March 2019, there were 1,554 laboratory notifications of infection [59]. A total of 315 cultures were received from participating diagnostic laboratories across New Zealand for sequencing and susceptibility testing. One isolate was removed from the study due to contamination. The majority (47.8%) of samples came from the Auckland region, the most densely populated city in New Zealand, where the highest number of gonorrhoea cases are recorded. The remaining isolates were obtained from outside Auckland (Table S3). Nearly three quarters (n=230, 73.3%) of isolates were from male patients and the urethra was the most commonly sampled body site (n=157, 50.0%), followed by anorectal (9.9%) or penile (9.2%) samples. For the 80 isolates from female patients, 74 were isolated from vaginal or cervix samples.

High-quality genome sequences were obtained for 314 isolates and obtained draft genome assemblies had a mean number of 94 contigs (±6.71 sd) and a median N50 of 49,810 bp (Table S4). Genome coverage to the FA1090 reference was 95.7% (±0.2% sd) with median depth of coverage of 135 (IQR 77).

A total of 98.4% of isolates were susceptible to azithromycin, with five having MICs above the EUCAST ECOFF (MIC >1 mg l−1). For the antibiotics cefixime and ceftriaxone, 99.0 and 99.4% of the samples were susceptible, with three and two isolates showing decreased susceptibility: MIC ≥0.12 mg l−1 and ≥0.06 mg l−1, respectively. Higher levels of resistance were observed for tetracycline (35.0%), ciprofloxacin (25.2%) and penicillin (8.3%) (Table 1). β-lactamase activity was found for 19/26 (61.5%) penicillin-resistant isolates. Multi-drug resistance was observed for 18 isolates in this study, with one isolate having tested resistant to four antibiotics (azithromycin, ciprofloxacin, penicillin and tetracycline) and 17 isolates were resistant to three antibiotics (ciprofloxacin, penicillin and tetracycline), of which one isolate additionally had decreased susceptibility to cefixime and ceftriaxone (Fig. 1).

Table 1.

Results of phenotypic susceptibility testing of 314 isolates

Antibiotic	MIC range (mg l−1)	MIC50 (mg l−1)	MIC90 (mg l−1)	Percent (no.)
				Susceptible	Intermediate/ decreased susceptibility*	Resistant
Penicillin	0.06–16	0.25	1	1.3 (4)	90.4 (284)	8.3 (26)
Cefixime	0.004–0.12	0.008	0.03	99.0 (311)	1.0 (3)	0.0 (0)
Ceftriaxone	0.002–0.06	0.004	4	99.4 (312)	0.6 (2)	0.0 (0)
Ertapenem†	0.004–0.12	0.016	0.03	–	–	–
Ciprofloxacin	0.002–8	0.004	4	74.8 (235)	0.0 (0)	25.2 (79)
Tetracycline	0.06–16	0.75	32	50.0 (157)	15.0 (47)	35.0 (110)
Spectinomycin	8–32	32	32	100.0 (314)	0.0 (0)	0.0 (0)
Azithromycin	0.03–8	0.12	0.5	98.4 (309)	–	1.6 (5)
Gentamicin†	2–16	8	8	–	–	–

* The ‘decreased susceptibility’ category applies to cefixime and ceftriaxone only.

† No interpretive standards available for ertapenem and gentamicin.

Fig. 1.

Population structure of circulating in New Zealand in 2018–2019. The maximum-likelihood phylogeny is based on 314 . isolates and the reference strain FA1090. The ML tree was built using iqtree and is based on a recombinant-free core genome alignment of 5,059 sites. The inner ring represents the sex of the patient, the outer ring illustrates the first level of clustering for Bayesian analysis of population structure (BAPS).

–

–

–

–

–

–

–

In silico genotyping

In silico typing was performed by determining allelic profiles and STs for three typing schemes: MLST, NG-MAST and NG-STAR.

MLST based on seven housekeeping genes (abcZ, adk, aroE, fumC, gdh, pdhC, pgm) revealed a total of 44 STs, with 38 previously described STs encompassing 305 isolates (97.1%) and an additional six novel STs identified for nine isolates (Table S5). The most common STs were ST8156 (14.0%) and ST13914 (11.2%). Grouping isolates based on single-locus variations in allelic profiles revealed four clonal complexes (CCs) comprising 230 isolates and 10 singleton MLST STs (84 isolates).

A total of 89 NG-MAST types were found, of which 44 STs were novel (Table S5). The most common sequences types were ST21017 (9.2%), ST13489 (8.3%), ST21010 (6.7%) and ST2400 (6.1%). A total of 55 NG-MAST types were unique and a further 11 STs were found in two isolates each. Grouping isolates based on single-locus variation in allelic profiles resulted in 14 clonal complexes (291 isolates) and 15 singleton NG-MAST STs (23 isolates).

Diverse NG-MAST types were found in isolates with decreased susceptibility to extended cephalosporins or azithromycin resistance. Three distinct NG-MAST types were identified among isolates with decreased susceptibility to extended cephalosporins, which belonged to distinct clonal complexes. ST1513 was found for two isolates, one of which had decreased susceptibility to both ceftriaxone and cefixime and one isolate with decreased susceptibility to only cefixime. ST10386 was identified for an isolate with decreased susceptibility to ceftriaxone and ST18982 was found for another isolate with decreased susceptibility to cefixime. Two isolates with decreased cefixime susceptibility and same NG-MAST type were referred from the Auckland region and Waikato. The five phenotypically azithromycin-resistant isolates had NG-MAST types distinct from each other and distinct from NG-MAST types found in isolates with decreased susceptibility to ceftriaxone and cefixime.

Typing for antimicrobial resistance using the NG-STAR scheme found 62 NG-STAR types with 20 new STs identified (Table S5). The most common STs were ST442 (13.7%), ST4006 (9.6%) and ST231, ST158 with 7.0% each. Grouping isolates based on single-locus variation in allelic profiles resulted in nine clonal complexes (194 isolates) and 31 singletons (120 isolates).

Genetic basis of antimicrobial resistance

All isolates underwent genome screening for presence or absence of plasmids and of genes known to be involved in reduced susceptibility or resistance to antibiotics, as well as chromosomal point mutations (Table S2). The genetic basis of antimicrobial resistance identified are summarized in Table 2 along with correlation to the resistance phenotype.

Table 2.

Number of isolates carrying plasmids or point mutations associated with antibiotic resistance and further correlation with resistant phenotypes

Antimicrobial	Gene (protein)	Mutation	No. (%)	Phenotype resistant no. (%)	Phenotype intermediate resistance*/decreased susceptibility†
Penicillin total N=314	bla TEM-1	–	18 (5.7 %)	18 (100 %)	–
	bla TEM-135	–	1 (0.32 %)	1 (100 %)	–
	mtrR promoter	–	122 (39.9 %)	4 (3.3 %)	118 (96.7 %)
	mtrR (MtrR)	A39T	113 (36.0 %)	19 (16.8 %)	94 (83.2 %)
		G45D	25 (8.0 %)	–	25 (100 %)
	ponA (PBP1)	L421P	99 (31.5 %)	10 (10.1 %)	89 (89.9 %)
	penA (PBP2)	A502V	4 (1.3 %)	3 (75 %)	1 (25 %)
		G543S	8 (2.5 %)	5 (62.5 %)	3 (37.5 %)
		G545S	2 (0.6 %)	1 (50 %)	1 (50 %)
		I312M	3 (1.0 %)	1 (66.6 %)	2 (33.4 %)
		P552S	2 (0.6 %)	1 (50 %)	1 (50 %)
		V316T	3 (1.0 %)	1 (66.6 %)	2 (33.4 %)
	penB (PorB1b)	A121D/N	45 (14.3 %)	7 (15.6 %)	38 (84.4 %)
		G120K/N	49 (15.6 %)	10 (20.4 %)	39 (79.6 %)
	pilQ (PilQ)	E666K	–	–	–
Ciprofloxacin total N=314	gyrA (GyrA)	D95G/A/Y	78 (25.2 %)	78 (100 %)	na
		S91F	79 (25.2 %)	79 (100 %)	na
	parC (ParC)	D86N	19 (6.1 %)	19 (100 %)	na
		E91K	3 (1.0 %)	3 (100 %)	na
		S87I/N/R	25 (8.0 %)	25 (100 %)	na
	parE (ParE)	–	–	–	na
	norM promoter	–	1 (0.3 %)	1 (100 %)	na
Tetracycline total N=314	tetM	–	80 (25.5 %)	80 (100 %)	–
	rpsJ	V57M	193 (61.5 %)	110 (57 %)	47 (24.4 %)
	mtrR promoter	–	122 (39.9 %)	65 (53.3 %)	10 (8.2 %)
	mtrR (MtrR)	A39T	113 (36.0 %)	43 (22.3 %)	22 (11.4 %)
		G45D	25(8.0 %)	1(4 %)	–
	penB (PorB1b)	A121D/N	45 (14.3 %)	31 (68.9 %)	12 (26.7 %)
		G120K/N	49 (15.6 %)	35 (71.4 %)	12 (24.5 %)
	pilQ (PilQ)	E666K	–	–	–
Azithromycin total N=314	23S	C2611T	4 (1.3 %)	4 (100 %)	na
	mtrR promoter	–	122 (39.9 %)	1 (0.8 %)	na
	mtrR (MtrR)	A39T	113 (36.0 %)	0 (0 %)	na
		G45D	25 (8.0 %)	1 (4 %)	na
	ermB	–	–	–	na
	ermC	–	–	–	na
	ermF	–	–	–	na
	macAB promoter	–	–	–	na
	mef	–	–	–	na
Ceftriaxone/	mtrR (MtrR)	A39T	113 (36.0 %)	–	1 (0.9 %)
Cefixime total N=314		G45D	25 (8.0 %)	–	1 (4 %)
	ponA (PBP1)	L421P	99 (31.5 %)	–	3 (3 %)
	penA (PBP2)	A502V	4 (1.3 %)	–	1 (25 %)
		G543S	8 (2.5 %)	–	–
		G545S	2 (0.6 %)	–	2 (100 %)
		I312M	3 (1.0 %)	–	2 (66.7 %)
		P552S	2 (0.6 %)	–	1 (50 %)
		V316T	3 (1.0 %)	–	2 (66.7 %)
	penB (PorB1b)	A121D/N	45 (14.3 %)	–	3 (6.7 %)
		G120K/N	49 (15.6 %)	–	3 (6.1 %)
Spectinomycin total N=314	16S	C1192T	–	–	na
	rpsE	–	–	–	na

*lntermediate resistance applies to penicillin and tetracycline.

†Decreased susceptibility applies to ceftriaxone and cefixime.

na, not applicable.

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

–

Penicillin

The majority of isolates (n=284, 90.4%) showed intermediate susceptibility (MIC 0.12–1.0 mg l−1), four were completely susceptible, and 26 were resistant to penicillin. The bla TEM-1 gene was identified in 18 isolates, including all isolates with high-level resistance to penicillin (MIC >4 mg l−1) and one isolate with a MIC of 16 mg l−1 carried the bla TEM-135 gene. Further investigation of plasmid types showed all isolates with bla TEM-1 likely were of the African type, with the exception of two isolates, which failed to produce a PCR product likely due to contig breaks in genome assemblies. The isolate with the bla TEM-135 allele was associated with the Rio/Toronto plasmid type.

The point mutation L421P in PBP1 (penicillin-binding protein 1 a) was found for 99 isolates (33.5%). No isolates carried any mutations in the pilQ gene encoding the pore-forming secretin pilQ of the type IV pili. Mutations in the mtr repressor result in overexpression of the mtrCDE efflux pump and 43.0% of isolates possessed a mutation in mtrR (A39T, n=113 or G45D, n=25). A single base-pair deletion in the mtrR promoter (-A) was present in 122 isolates. Three isolates had no mapping reads to the mtrR promoter, likely due to sequencing errors. A total of 48 isolates had the G120K mutation in porB1, an outer membrane porin, of which 44 isolates carried additionally either A121D (n=42) or A121N (n=2). One isolate carried point mutations A121N and G120N in porB1b. For 4.1% of isolates (n=13), no sequencing reads mapped to porB1b and for 5.7% (n=18) porB1b assembled with low coverage (<40x) with no variants identified.

Tetracycline

Phenotypic tetracycline resistance was observed in 35.0% of isolates (n=110), with 15.0% of isolates (n=47) having an intermediate phenotype. The presence of the tetM gene and plasmid was confirmed for all high-level resistant isolates (MIC >8 mg l−1, n=80), distinguishing between the American (n=6) and Dutch type (n=74). A total of 193 isolates carried the V57M mutation in rpsJ (30S ribosomal protein S10) (61.6%) that reduces affinity of tetracycline to the 30S ribosome. No known mutations conferring resistance were observed for pilQ in any isolate. A total of 44.0% of isolates have mutations in mtrR (n=113 A39T mutation, n=25 G45D).

Spectinomycin

All isolates were susceptible to spectinomycin and none of the isolates carried either the C1192T mutation in the 16S rRNA gene or mutations in the rpsE gene conferring resistance to spectinomycin.

Azithromycin

Five isolates showed resistance to azithromycin with an azithromycin MIC greater than the 1 mg l−1 ECOFF. A total of four isolates carried mutations in the 23S rRNA gene with C2611T present in four copies for three isolates and in three copies for one isolate. No acquired azithromycin resistance genes (ermF, ermB, ermC, mefA/mefE, ereA, ereB) were found in any isolate. No isolates had mutations in rplD, tandem duplications for rplV or promoter mutations in macAB. However, one isolate with an azithromycin MIC greater than 1 mg l−1 did not carry any mutation associated with resistance. No sequencing reads for this isolate mapped to the mtrR promoter, so it remains possible that the isolate has the -A deletion in the mtrR promoter.

Ciprofloxacin

Ciprofloxacin resistance was identified for 78 isolates (25.0%) with a MIC >0.06 mg l−1. Mutations in gyrA reduce quinolone binding to DNA gyrase with 79 isolates carrying S91F mutations alongside either D95G (n=37), D95A (n=41) or D95Y (n=1). No mutations in parE were detected. A total of 44 isolates carried mutations in parC, which reduces quinolone binding to topoisomerase IV subunit A, with mutations in D86N (n=19), E91K (n=3), S81I (n=1), S87R (n=16) and S87N (n=17). Surprisingly, one isolate had a mutation in norM (MIC=8 mg l−1), alongside additional mutations in parC (S87N), gyrA (S91F and D96G).

Cefixime and ceftriaxone

Three isolates had decreased susceptibility to cefixime (≥0.12 mg l−1) and two isolates had decreased susceptibility to ceftriaxone (≥0.06 mg l−1). One isolate with decreased susceptibility to both cefixime and ceftriaxone had the same NG-MAST, MLST and NG-STAR ST as an isolate with reduced susceptibility to cefixime. Both carried the mosaic XXXIV penA allele with mutations in G545S, I312M and V316T and also had mutations in porB1b (A121N, G120K) and ponA (L421P). Of the remaining two isolates, one isolate with reduced susceptibility to ceftriaxone carried mutations in porB1b (G120K and A121D), ponA (L421P) and penA (A502V, P552S), whereas the other isolate with reduced susceptibility to cefixime had only the mtrR A39T mutation that is associated with beta-lactam resistance.

Others

The folP mutation (R228S) resulting in resistance to sulphonamides was carried by 272 isolates (86.6%). Resistance to rifampicin (H553N in rpoB) was detected in 21.0% of isolates (n=66). The highly conserved porA pseudogene is a commonly used target for molecular detection of gonorrhoeae, however there have been reports of test failures due to the presence of meningococcal porA in isolates [60]. All isolates in this study possessed the gonococcal porA.

Population structure

A maximum-likelihood (ML) phylogeny (Fig. 1) of the 314 isolates and the reference genome FA1090 was inferred based on the core-genome alignment of 5,059 sites (including 3,871 informative sites) after removal of recombinant regions. The phylogeny was significantly structured by sexual behaviour (P<0.001, Fritz and Purvis D). Isolates with resistance to penicillin, ciprofloxacin and tetracycline are clustering together (Fig. S1), while azithromycin-resistant isolates were in distant parts of the tree, with the exception of two isolates (Fig. S2). A hierarchical-clustering analysis identified 15 clusters largely in accordance with phylogenetic grouping, with the exception of cluster 11 and 8, which were split up in different parts of the tree (Fig. 1). While two clusters were only associated with males, the rest was mixed clusters of both sexes, however there was significant association of sex and sexual orientation with cluster IDs (P<0.001, Fisher’s exact test). Levels of diversity of circulating are similar (Fig. 2), when comparing isolates with data from a previous survey undertaken in New Zealand in 2014/2015 [23].

Fig. 2.

Maximum-likelihood phylogeny for the combined dataset from both national surveys. In total, 314 isolates from this dataset were combined with 401 isolates from the 2014/2015 study [23]. ML phylogeny was built based on a core-genome alignment of 4813 sites with 3269 informative sites. All samples from the 2014/2015 survey are highlighted in yellow and the isolates from 2018/2019 are illustrated in blue (inner ring). The outer ring illustrates the seven mixed transmission clusters, with isolates from 2014/2015 and 2018/2019.

Transmission analysis

A total of 49 transmission clusters were identified that comprised 85.7% (269/314) of isolates, with a median size of three isolates (range 2–29) (Table S6). Ten clusters were associated with MSM only with the largest being cluster 1 with 16 isolates (9=MSM, 7=NA, where no sexual orientation data was available). The other seven largest clusters (>nine isolates) were a mixture of sexual risk groups with MSM dominating in clusters 30 and 4 (35.7 and 45.0%, respectively) (Fig. 3), suggesting bridging across sexual networks. Clusters 52, 23 and 13 were only associated with heterosexual orientation (both male and female). The age range was more restricted in heterosexual clusters (median 26, range 17–38), compared to MSM, heterosexual and female mixed clusters (median 29, range 0–63), or the MSM only cluster (median 26, range 18–48). Clustering analysis including the isolates from the 2014/2015 national survey identified 96 clusters (SNP threshold=10, excluding singleton clusters). Seven clusters comprised isolates from both time periods (Fig. S3), indicating that transmission of strains might be ongoing, despite the time difference in sampling.

Fig. 3.

Transmission cluster analysis based on pairwise SNP-distances and a cut-off threshold of ten SNPs. (a) Highlighting clusters with>=9 isolates. (b) Cluster composition of patients sexual orientation. Cluster composition of patients sexual orientation. (c)Age range of patients whose isolates were assigned to cluster and sexual orientation highlighted. Age range of patients whose isolates were assigned to cluster and sexual orientation highlighted.

Diversity of isolates in sequence types

Finally, the pairwise genetic distance between isolates of the same sequence type were plotted for each of the typing schemes used in this study (MLST, NG-MAST, NG-STAR), highlighting the fact that some sequence-type schemes are more heterogenous compared to others (Fig. 4), which was especially observed for the MLST-typing scheme with 0–805 SNP difference between isolates of the same ST. Lower levels of heterogeneity were found for NG-MAST- and NG-STAR-typing schemes with pairwise SNP differences between 0 and 105 for NG-MAST, after exclusion of one outlier isolate in ST11461 with 836 and 841 SNPs' pairwise difference to the other isolates in that ST, and 0 and 153 for NG-STAR STs.

Fig. 4.

Pairwise SNP distance identified within MLST, NGSTAR and NGMAST sequence types. From each calculated pairwise SNP distance, both isolates had to be of the same type (MLST, NGSTAR or NGMAST) to be grouped. Only sequence types with n>1 are displayed.

Discussion

We report the findings of a national survey of circulating in New Zealand in 2018/2019. A detailed analysis of variation in genes coding for antimicrobial resistance was undertaken and correlated with phenotypic susceptibility testing. Patterns of transmission were determined combining WGS and epidemiological data to improve our understanding of transmission networks and antimicrobial resistance. Additionally, this study presented the unique opportunity to investigate contemporary temporal changes of diversity and resistance rates, by comparing our study with a national survey undertaken in New Zealand in 2014/2015 [22, 23].

Overall, we found a strong correlation between the antibiotic resistant phenotype and the genetic determinants associated with resistance, with the exception of one isolate that was resistant to azithromycin, but had no known associated mutation. Penicillinase-producing had high-level resistance to penicillin (MIC >4 mg l−1), conferred by a TEM-1 type β-lactamase plasmid. Importantly, one isolate with a penicillin MIC of 16 mg l−1 carried a bla TEM-135 type plasmid, which carries a single nucleotide substitution at position 135 (M182T) [61], and was found to be of the Rio/Toronto type. While bla TEM-135 is common in Australia [62], this type had not been observed in New Zealand before, suggesting that the patient might have had a travel history with Australia or had contact with someone who had recently travelled outside New Zealand. This is of particular concern, because resistance to extended-spectrum cephalosporin can evolve from a single nucleotide mutation [63]. One isolate had the promoter mutation in the norM gene encoding the sodium gradient-dependent MATE family efflux pump NorM [64], which is associated with increased resistance to fluoroquinolones. The isolate had a MIC of 8 mg l−1 to ciprofloxacin, though mutations in parC (S87N) and gyrA (S91F and D96G) were also present.

Comparing the observed AMR resistance rates to the earlier national survey completed in 2014/2015 [23] and applying the same MIC thresholds as used in our current study shows that azithromycin resistance rates have increased slightly to 1.6% in 2018/2019 (Table S7). This is still lower than the observed azithromycin resistance rates in Australia of 9.5% [65]. No isolates with high-level, acquired azithromycin resistance have been found in New Zealand. Reduced susceptibility rates to cefixime (1.0%) and ceftriaxone (0.6%) decreased since 2014/2015 and the ceftriaxone rate is comparable to Australia (1.1%) [63]. We observed an increase in resistance to tetracycline (35.0%), but a decrease in resistance to penicillin (8.3%) and ciprofloxacin (25.2%). The trend of decreasing rates of ciprofloxacin resistance has also been described in Australia [65]. Overall, the New Zealand rates of azithromycin resistance and decreased susceptibility to ceftriaxone remain lower than that in most countries that are part of the WHO Global Gonococcal Antimicrobial Surveillance Programme (GASP) [66], which suggests that our current treatment guidelines of a dual therapy of azithromycin and ceftriaxone remain effective.

We found a significant association of phylogeny with sexual behaviour groups (P<0.001, Fritz and Purvis D), and clustering analysis (NG-MAST and BAPS) also showed significant associations with sexual behaviour and transmission clusters (P<0.001, Fisher’s exact test). No isolates with NG-MAST ST1407, a multi-drug-resistant sequence type, which has caused particular concern in Europe [67], were detected in New Zealand. However, NG-MAST ST1513 defined for two isolates with decreased susceptibility to ESC, has been shown to belong to the same NG-MAST group as ST1407, when grouped according to relative distances on a phylogenomic tree [68]. Temporal variation in prevalent STs was prevalent in New Zealand when comparing NG-MAST types of 2018/2019 (ST21017, ST13489, ST21010, ST2400) to the 2014/2015 survey (ST4186, ST2400), which has been described in other countries [67]. Comparing the most common NG-MAST types circulating in NZ with the six most common types in Australia [11] found no overlap, but variations in prevalence of NG-MAST types between countries is also common [69].

Whilst nine out of 49 identified transmission clusters were solely associated with MSM, the other clusters included heterosexual patients, as well as MSM, suggesting bridging of sexual transmission networks between heterosexual, bisexuals and homosexuals. This is in accordance to earlier findings in New Zealand in 2014/2015, where MSM were not the sole driving force of gonorrhoea transmission [23], despite incidence rates being highest among MSM networks [70]. This cross-over in sexual networks is increasingly being described overseas with mixed clusters of both MSM and heterosexuals [11, 55] and an important driver for the spread of AMR to consider. NG-MAST types shared a close evolutionary history, with a median nucleotide divergence within a single ST from 1.0 to 836.0 SNPs (overall median 5.0), followed by NG-STAR types with a median 0.0 to 52.5 (overall median 5.0), whereas for MLST types, there is a lot more variation (median range 2.0 to 384.5 with an overall median of 13.0). This underlines previous findings of the NG-STAR typing scheme being superior to the two schemes in terms of genomic inherency [71].

Two isolates with azithromycin resistance were assigned to the same transmission cluster, which coincides with the typing results, as both had the same MLST and NG-STAR STs, but different NG-MAST STs. The remaining three azithromycin-resistant isolates did not cluster with any other isolates included in the study, suggesting the independent evolution of resistance to azithromycin or introduction from overseas.

Several limitations need to be considered when interpreting the results of this study: with the trend towards a nucleic acid-based diagnosis of gonorrhoea in New Zealand, no cultures were available for 77.9% of gonorrhoea cases that were notified during the survey period [59]. The geographic spread and number of included isolates are not necessarily consistent with notification rates and overall there was an over-representation of isolates from the Auckland region (47.77%), so the results might not reflect the national diversity. Males were also over-represented compared to laboratory notifications received, thought to have resulted from MSM being more likely to attend a sexual health clinic and have enhanced data collected [59]. Information about sexual behaviour was available for only 61.0% of the samples included in the study. Travel history data was not available, but may represent a risk factor for resistant , as shown for the first discovery of bla TEM-135 in New Zealand, while this plasmid type is common in Australia.

Molecular detection methods, such as nucleic acid amplification tests (NAAT), which are widely used for diagnostic purposes in New Zealand, typically include the gonococcal porA pseudogene as a target, because it is absent in commensal species and divergent enough from . Although false negative tests have been reported overseas due to the presence of the meningococcal porA gene [60], all surveyed isolates were confirmed to carry a gonococcal porA.

In conclusion, despite decreasing resistance rates for the antibiotics ceftriaxone, penicillin and ciprofloxacin and confirmation that despite an increase from 0.5–1.6% in azithromycin resistance compared to 2014/2015, the overall resistance rates remain low in New Zealand. However, the situation needs to be monitored and continued surveillance of gonococcal AMR and molecular epidemiology using culturing is required in New Zealand, particularly in regards to the concerning rates of resistance to azithromycin and ESC overseas and in Australia, given the extensive travel occurring between the two countries. This is particularly relevant for adapting treatment guidelines, if resistance rates reach the 5.0% of prevalence in the circulating population [15].

Click here for additional data file.

Click here for additional data file.