A Whole-genome Sequencing Analysis of Neisseria gonorrhoeae Isolates in China: An Observational Study.

BACKGROUND: Tracking the spread of the Neisseria gonorrhoeae strains with decreased susceptibility or resistance to cephalosporins is a major priority for global surveillance programmes. Whole-genome sequencing (WGS) has been widely used by increasing countries in North America, Europe, and Pacific to determine the decreased susceptible or resistance determinants of Neisseria gonorrhoeae, track the spread of these determinants throughout the gonococcal population at national or regional level. However, no studies to date have examined the genomic epidemiology of gonorrhea in Asia where the antimicrobial resistant strains of Neisseria gonorrhoeae appears to have emerged before disseminating the strains globally.
METHODS: We obtained clinical isolates and data from the China Gonococcal Resistance Surveillance Programme (China-GRSP) from 2012 to 2013. We sequenced the genomes of 435 clinical isolates of Neisseria gonorrhoeae, including 112 (25.6%) isolates with decreased susceptibility to ceftriaxone (Cfx-DS). We assessed the association between antimicrobial resistance genotype and phenotype. We also compared our data with the whole genome data of the isolates from the USA and the UK in the GenBank.
FINDINGS: The most prevalent MLST STs in our gonococcal population were MLST ST7827 (n = 74), followed by ST7365 (n = 58), ST1600 (n = 38), ST7367 (n = 35), and ST7363 (n = 29). MLST ST1901 which was reported as the predominant ST in the US was not found in our population. A total of 2512 strains, including additional 2077 published NG strains, were further included for phylogenetic analysis. It generated two distinct lineages - lineage 1 and lineage 2. Analysis of MLST ST1901 in the database indicate that most of MLST ST1901 isolates in the lineage2.6 were Cfx-DS isolates while all isolates in the lineage 2.1 were sensitive to ceftriaxone (77/110 vs. 0/13; p < 0.001). ST1901/lineage 2.6 is a ceftriaxone resistant clone which cannot distinguished by MLST genotyping. In the isolates from our study, the MICs of ceftriaxone for ST7363/lineage 2.6 isolates ranged from 0.008-0.125 mg/L (mean ± SD; 0.054 ± 0.043 mg/L) while those for ST7363/lineage 2.8 isolates ranged from 0.032-0.250 mg/L (0.134 ± 0.085 mg/L). All ST7363/lineage 2.8 isolates contained penA mosaic alleles.
INTERPRETATION: To our knowledge, current study is the first WGS-based analysis of gonococcal population at national level in Asia. China harbors the different predominant clones associated with decreased susceptibility to ceftriaxone from those clones circulated in other regions. The findings from the study can be not only used as baseline data for future studies in China but also contributable to our understanding on spread of N. gonorrhoeae and its resistant strains at regional and global levels. FUNDING: The Chinese Academy Medical Sciences (CAMS) Initiative for Innovative Medicine.

Research in context

Evidence before this study

We searched PubMed for publications up until 15 June 2018 with the terms (“Neisseria gonorrhoeae” OR “gonorrhea” OR “gonorrhea”) AND (“whole genome sequencing” OR “WGS”), and references and subsequent citations (identified using Google Scholar) were also reviewed. Genomic epidemiological studies using whole genome sequencing (WGS) of Neisseria gonorrhoeae to investigate the spread of drug resistant strains have been conducted at a national level in countries within North America (the USA, and Canada), Europe (the UK, and Ireland), and the Pacific (New Zealand and Australia). To our knowledge, there have been no published national studies that used WGS of N. gonorrhoeae in any countries in Asia. Previous studies have used WGS to determine the resistance determinants of N. gonorrhoeae, track the spread of these determinants throughout the N. gonorrhoeae population at national or regional level, and investigate the local outbreaks of gonorrhea-resistant strains, highlighting the advantages and potential of WGS to be applied in strengthening the surveillance of gonococcal AMR at local and international scales. No study to date has systematically applied WGS database to compare the predominant clones across the countries.

Added value of this study

To our knowledge, we report the first WGS-based epidemiological study on AMR of N. gonorrhoeae at national level in Asia. We describe the genomic distribution of the Chinese gonococcal population which can be used as baseline data for a longitudinal genomic surveillance in China. We demonstrate a high number of different sequence types, including newly types, in our study population and a different distribution of NG-MAST STs and lineages from that in the US and the UK, and the specific ST predominated isolates (ST7363/lineage 2.8), harboring the penA mosaic gene, to be associated with high MIC levels to ceftriaxone.

Implications of all the available evidence

WGS has been widely used by increasing countries into their N. gonorrhoeae surveillance programmes to enhance our understanding of the distribution of gonococcal AMR clones in different regions in the world. While the WGS-based surveillance data have revealed the dissemination of the gonococcal strains with decreased susceptibility or resistance to ceftriaxone globally, the genomic heterogeneity among these gonococcal strains exists across the countries. Comparisons of findings from genomic epidemiological studies in a geographic context could provide evidence aimed at the early identification of transmission and spread of N. gonorrhoeae resistant strains at local, regional and global levels.

Alt-text: Unlabelled Box

Introduction

Gonorrhea, caused by the Gram-negative bacterium Neisseria gonorrhoeae (NG), remains a major public health concern globally. The World Health Organization (WHO) estimated in 2012 that there were 78 million new gonococcal infections among people 15–49 years old annually, 45% of which occurred in the WHO Western Pacific Region (WPR) [1]. As the most populous country in the region, China contributes substantially to the number of people with gonorrhea in the WPR. The reported incidence of gonococcal infections in China increased by 20.7% from 2016 to 2017 (unpublished data from the National Center for Sexually Transmitted Disease [STD] Control in China). Case management through correct diagnosis and effective treatment has been listed as priority component for immediate action to effectively control sexually transmitted infections (STIs) including gonorrhea. Unfortunately, rapid development and spread of gonococcal antimicrobial resistance (AMR) have been substantially compromising the effectiveness of NG treatment with the antibiotics and consequently threatening the NG control strategies [2]. Of particular concern is the fact that decreased susceptibility of or resistance to extended-spectrum cephalosporins (ESCs), usually the last line of available monotherapy, has been reported in all WHO regions [3].Global efforts to monitor trends in gonococcal AMR, to improve the quality, comparability, and timeliness of gonococcal AMR data across countries, and to assess resistance patterns in key populations at highest risk for gonococcal AMR have been established through the Enhanced Gonococcal Antimicrobial Surveillance Program [4]. Molecular methods have been introduced as powerful tools to enhance the current gonococcal AMR surveillance program in which NG multi-antigen sequence typing (NG-MAST) and multilocus sequence typing (MLST) are widely used to identify the gonococcal AMR strains. Recently, the use of the large-scale next-generation whole-genome sequencing (WGS) technology has been shown to provide a better resolution of strains for comparative studies. The WGS-based epidemiological studies on gonococcal AMR have been conducted in the North American [5], [6], [7], European [8], [9] and Pacific countries [10], [11] in recent years. To date, however, no WGS-based studies have been published from any countries in Asia where gonococcal AMR appears to have originally emerged before disseminating globally [12], [13]. The current study was aimed to describe the genomic diversity of gonococcal population in China and to compare our results with the genomic data from the USA and the UK in the GenBank.

Materials and Methods

Clinical Isolates

A total of 438 NG isolates were collected from patients (one isolate per patient) attending the clinic sites in eleven provinces within the China Gonococcal Resistance Surveillance Programme (China-GRSP) during 2012 and 2013. These isolates were a subset of our previous study [14], [15]. The strains were randomly selected at the ration of 1 (resistant of decreased susceptible):3 (susceptible). Geographic distribution of the participating provinces and their corresponding numbers of NG isolates available for the study are shown in Fig. 1.

Fig. 1

Geographic locations of the provinces where the N. gonorrhoeae isolates were collected for the study.

The number of isolates from each province is given in parentheses.

Antimicrobial Susceptibility Testing

Antimicrobial susceptibility to ceftriaxone was determined using agar dilution method according to the recommendations from the WHO [16]. The minimum inhibitory concentration (MIC) was defined as the lowest concentration of an antibiotic that inhibited the growth of gonococci. We categorized MICs as ceftriaxone susceptible (MIC ≤ 0.06 mg/L, Cfx-S) or decreased susceptibility (MIC ≥ 0.125 mg/L, Cfx-DS) in accordance with the interpretive criteria of the Clinical and Laboratory Standards Institute [16], [17]. For quality assurance, all the central STD laboratories participated in the external quality assurance program of the WHO Western Pacific GASP through the National STD Reference Laboratory in Nanjing [15].

DNA Extraction

DNA was extracted from bacterial suspensions using the QIAxtractor DX Kits (Qiagen, Hilden, Germany) on a QIAxtractor automated genomic DNA extraction instrument, according to the manufacturer's instructions (Qiagen, Hilden, Germany).

Sequencing, Alignment and SNP Calling

Genomic DNA was sequenced on the Illumina HiSeq 2500 platform, generating single end reads of 101 bp in length. For each sample duplicate reads were removed by custom Perl scripts (available from http://www.mgc.ac.cn/Resources/NGscripts.tgz). Further quality control was conducted using the NGSQC Toolkit with a cutoff of Q20 [18]. Valid reads were then aligned to the reference genome sequence of FA1090 (GenBank accession No. NC_002946) using the Burrows-Wheeler algorithm as implemented in BWA [19]. For each strain ~ 3.9 (2.6–7.1) millions of valid reads were produced to yield ~ 99.2% (98.8%–99.8%) coverage of the reference genome with an average depth of 172X (79X-306X). SNPs were identified with a minimum depth of 10 × and a consensus quality score of 50 using SAM tools [20]. SNPs located within repetitive regions and prophages of the FA1090 genome identified by RepeatMasker (http://www.repeatmasker.org/) or PHAST [21], were excluded. Mixed base calls were considered valid only if the numbers of the most abundant (n1) and the second most abundant (n2) nucleotides at each SNP in each strain satisfied the criteria n1/n2 ≧ 5. All sequences generated in this study have been submitted to GenBank under the accession number SRP133345.

Phylogenetic Analysis

Additional 2077 published NG strains, including 18 complete genomes (Table S1) and 2059 raw sequencing data from two previous studies [22], [23] were further included for phylogenetic analysis (Table S2). A concatenate superset of SNPs relative to FA1090 was generated across all newly sequenced strains and 2077 published NG genomes. SNP sites with missing data in over 5% of the strains within the dataset were removed. To avoid the potential effects of homoplasy of drug resistance-associated mutations in phylogeny, SNPs located in the known genes/regions, including penA, porB, mtrR and mtrR promoter were further excluded from the dataset for phylogenetic tree construction. Four representative strains of N. meningitidis were included as outgroups: Z2491, MC58, FAM18 and 053442 (GenBank accessions Nos.: NC_003116, NC_003112, NC_008767 and NC_010120). The refined SNP set was used to construct the maximum-likelihood phylogeny using RAxML under the GTR gamma substitution model as previously described [24], [25]. It is well established that maximum-likelihood methods generally outperform distance and parsimony methods over a broad range of realistic conditions [26], [27]. The reliability of each node was tested via a bootstrap analysis on 100 resampled datasets. The online iTOL platform was used for further phylogenetic tree visualization and annotation [28].

Molecular Typing

To determine the epidemiological relatedness, NG MLST and NG-MAST were performed for all isolates according to the protocols outlined on the respective database web site (MLST, http://pubmlst.org/neisseria/, NG-MAST, http://www.ng-mast.net). Furthermore, penA genotypes were also determined according to previous report [29]. Sanger sequencing was used to genotype the isolates. For MLST data, eBURST analysis was conducted to discern the evolutionary patterns and explore the founding genotypes. The input STs were subdivided into groups under the most stringent group definition, which means STs within the same group must share at least six identical alleles with at least one other ST in the group. Default setting was used to define subgroup founders, which means subgroup founders must have three direct links to other STs. Reliability of the founding genotypes was assessed by bootstrap resampling procedures. The eBURST analyses were performed as per the instructions at http://eburst.mlst.net.

Statistical Analysis and Ethics Approval

Student's t-test, and chi-square or Fisher's exact test were used as appropriate. A p-value less than 0.05 was considered statistically significant. Corrected p-value (q-value) was obtained using the Benjamini–Hochberg false discovery rate (FDR) approach. Statistical analysis was performed using SPSS for Windows (version 16.0; SPSS, Chicago, Illinois). The study protocol was reviewed and approved by the Medical Ethics Committee of the Chinese Academy of Medical Sciences (CAMS) Institute of Dermatology and the National Center for STD Control at Nanjing, China (approval number 2011-LS-003).

Role of the Funding Source

The funder of the study had no role in study design, data collection, data analysis, data interpretation, or writing of the report. The corresponding authors had full access to all the data in the study and had final responsibility for the decision to submit for publication.

Results

The demographic characteristics of 438 patients (one isolate per patient) were summarized in Table 1. Briefly, 90.6% were males and 8.7% were females. The mean age was 35.81 ± 11.68 (standard deviation, SD) years (ranging from 15 to 86 years). Of the 438 patients, 84.7% reported heterosexual orientation. Majority of the patients (91.3%) recognized their infections with NG from the local sexual contacts. These demographic and behavioral characteristics are not different between groups with Cfx-S and Cfx-DS isolates.

Table 1

Characteristics of patients included in the study.

Characteristic	Value	Susceptibility to ceftriaxone		P value*
		Decreased susceptible(n = 112)	Susceptible (n = 326)
Age at diagnosis, years
Range	15–86	20–78	15–86	0.53
Mean (± SD)	35.81 ± 11.68	35.21 ± 10.97	36.02 ± 11.88
Sex
Male	397(90.6%)	104	293	0.294
Female	38(8.7%)	7	31
Unknown	3(0.7%)	1	2
Nationality
Han	418(95.4%)	111	307	0.085
Other	18(4.1%)	1	17
Unknown	2(0.5%)	0	2
Sexual orientation
Heterosexual	371(84.7%)	103	268	1.000
Homosexual	6(1.4%)	1	5
Unknown	61(13.9%)	8	53

NG-MAST and MLST Genotyping

Three isolates failed in whole genome sequencing and were excluded from further analysis (Table S3). Four hundred and thirty-five isolates (99.3% of 438 isolates) were assigned to 60 different MLST sequence types (STs), of which 14 were newly reported in the Pubmlst database. The eBURST analysis result showed that all 60 STs could be subdivided into one group and ST7365 was predicted as the founding genotype (Fig. S1). These new MLST STs represented 3.9% (17/435) of the isolates, most of which (14/17) were found in the isolates collected from the coastal areas where more cases of gonorrhea were reported than in other areas in the country [30].The most prevalent MLST STs (≥ 20 isolates) were MLST ST7827 (n = 74), ST7365 (n = 58), ST1600 (n = 38), ST7367 (n = 35), and ST7363 (n = 29). For MLST data of 4287 N. gonorrhoeae isolates in PubMLST database, eBURST analysis was conducted. 566 identified STs could be subdivided into six groups and most of which (517/566) were found in group 1. In group 1, ST1901, the predominant ST in the US, was predicted as the founding genotype, but was not found in our population (Fig. S2).There were 284 NG-MAST STs (including 4 new NG-MAST STs) identified in our study and the most prevalent STs were NG-MAST ST2318 (n = 17), followed by ST4846 (n = 9), ST2083 (n = 8) and ST1866 (n = 7).

Genomic Epidemiology

By comparing the phylogenetic information of N. gonorrhoeae isolates in our study with the whole genome data of the isolates from the USA and the UK in the GenBank (Fig. 2), a total of 2512 NG strains (Dataset 1 in Table S2) generated two distinct lineages - lineage 1 and lineage 2. The lineage 1, which represents about 19% of the strains, evolves earlier than the lineage 2 and represents most of the UK strains. The lineage 2 represents the strains predominantly circulated worldwide. In addition, the diversity in the lineage 2 is significantly higher than that in the lineage 1. The majority of the lineage 2 strains can be further grouped into 9 sub-lineages. We listed lineage determinant SNPs for lineage 1 and lineage 2 and sub-lineage specific SNPs of lineage 2 in Tables S4 and S5. The most abundant MLST STs are mapped into the phylogenetic tree and well correlate with the demarcation of lineages and sub-lineages (Fig. 2). Compared with the isolates with AMR information from the US (Dataset 2 in Table S2), all Cfx-DS strains from our study distribute in all 9 sub-lineages of the lineage 2 (Fig. S3) while those from the USA (n = 95) distribute only in the lineage 1 and the lineage2.1 (Fig. 3). Analysis of MLST ST1901 indicates that 110 and 13 ST1901 strains were assigned to the lineage 2.1 and the lineage 2.6, respectively. Most of MLST ST1901 isolates in the lineage 2.6 were Cfx-DS isolates while all isolates in the lineage 2.1 were sensitive to ceftriaxone (77/110 vs. 0/13; p < 0.001). ST1901/lineage 2.6 is a ceftriaxone resistant clone which cannot be distinguished by MLST genotyping. In the isolates from our study, the MICs of ceftriaxone for ST7363/lineage 2.6 isolates ranged from 0.008–0.125 mg/L (mean ± SD; 0.054 ± 0.043 mg/L) while those for ST7363/lineage 2.8 isolates ranged from 0.032–0.250 mg/L (0.134 ± 0.085 mg/L). The differences in MICs of ceftriaxone for our isolates in lineage 2.6 vs 2.8 are statistically significant (p < 0.01) (Table S6). All ST7363/lineage2.8 isolates contained penA mosaic alleles.

Fig. 2

Phylogenetic analysis. Four representative strains of N. meningitidis were included as outgroups.

Maximum-likelihood tree of 2512 N. gonorrhoeae isolates based on 69,374 genome-wide SNP sites. Four representative strains of N. meningitidis (Z2491, MC58, FAM18 and 053442)were used as outgroups. Lineage 1 (L1) and 9 sub-lineages of lineage 2 (L2) are marked with background shadows. The proportions of predominate MLST types of each lineage/sub-lineage are highlighted by internal pie charts with size proportional to the total number of isolates. The external color strips range from 1 (inner) to 3 (outer). Strip 1, isolates with predominate MLST types of the current dataset. Strip 2, origin of each strain with color code for countries: deep pink, China; sky blue, USA; light green, UK. Strip 3, distribution of N. gonorrhoeae strains isolated from southeast (orange) and northwest (blue) of China. The scale is in the units of mutations per site.

Fig. 3

Phylogenetic analysis of 635 NG isolates.

A degraded maximum-likelihood tree of 635 N. gonorrhoeae strains with available cephalosporin antibiotics susceptibility information. Four reference completes, FA1090, FA19, MS11 and NCCP11945 are also included. Colored dots at the leaves denote cephalosporin antibiotic resistance (red) or sensitive (blue) isolates. The external color strip indicates the sources of each isolates: deep pink, the current study; sky blue, USA cohort; gold, representative strains from WHO; gray, reference genomes. The scale is in the units of mutations per site.

Ceftriaxone Resistant Determinants

In addition to detecting the mutations in four ceftriaxone resistant determinants (penA, porB, mtrR and pilQ) (Table S7), 18 amino acids mutations in these four genes were newly found to be statistically associated with Cfx-DS in our study (Table 2). Characteristics of ceftriaxone resistant isolates were described in Table S8. A total of 1480 SNPs which were not reported in previous studies were statistically associated with Cfx-DS in our study (Q < 0.01).

Table 2

Suspected ceftriaxone genetic resistance determinants.

Suspected resistant determinants	Cfx-DS isolates	Cfx-S isolates	P-value	Q-value
PBP2
D100E	105/207	23/428	1.14E-40	0.00135
V159A	105/207	15/428	4.68E-46	0.00134
N172S	105/207	15/428	4.68E-46	0.00134
Q213E	102/207	15/428	3.37E-44	0.00133
A278V	101/207	15/428	1.38E-43	0.00132
D284E	97/207	14/428	7.27E-42	0.00132
R287K	90/207	14/428	1.07E-37	0.00131
R290Q	35/207	8/428	1.55E-12	0.00130
H541N	102/207	80/428	1.36E-15	0.00128
PORB
T40R	115/207	125/428	1.37E-10	0.00315
A54S	131/207	135/428	2.97E-14	0.00185
V75I	90/207	81/428	6.24E-11	0.00296
T89S	102/207	61/428	2.78E-21	0.00186
S217 N	97/207	107/428	3.21E-08	0.00490
A272V	109/207	110/428	2.1E-11	0.00268
S295D	94/207	97/428	4.66E-09	0.00415
MTRR
H105Y	131/207	210/428	0.00076	0.01500
PILQ
S333N	201/207	382/428	0.00072	0.01492

Discussion

As a key technology for providing a far greater degree of resolution and ruling out potential transmission links that may be inferred by the traditional genotyping methods (such as MLST and NG-MAST) [31], WGS has been widely used to determine the resistance determinants of N. gonorrhoeae, track the spread of these determinants throughout the N. gonorrhoeae population at national or regional level, and investigate the local outbreaks of gonorrhea-resistant strains [32], [33].Eyre et al. used WGS approach to predict antibiotic MICs for N. gonorrhoeae among specimens from England, the USA and Canada, and demonstrate this approach that allows reliable MIC prediction for gonorrhea antimicrobials including cefixime, and azithromycin [34]. Grad et al. performed a genomic epidemiological study of N. gonorrhoeae with reduced susceptibility to cefixime in the USA and found that WGS approach could help slow the transmission of antibiotic-resistant gonorrhea [35]. De Silva et al. carried out a WGS-based survey to track AMR of N. gonorrhoeae and gonorrhea transmission in the UK by including the WGS data from the USA for comparative analysis and demonstrate transmission of the infection locally, nationally and internationally [36]. Lee et al. conducted a genomic epidemiological study in New Zealand and identified several clusters of isolates with raised MICs to both ceftriaxone and cefixime, suggestive of de novo acquisition of the reduced susceptibility to cephalosporins, with subsequent transmission within clusters [10]. However, all of these published studies were conducted in countries within Europe, North America or Pacific. To our knowledge, the current survey in China is the first WGS-based epidemiological study on AMR of N. gonorrhoeae at national level in Asia. In this study, we have shown that our population covers a high number of different sequence types, including newly types, which were not found in China before, mostly detected from population in the coastal areas which are usually the destinations of national and international migrations. We have also shown the distribution of NG-MAST STs and lineages in our study population, which are different from that in the US and the UK. These findings could be used as baseline data for a longitudinal genomic surveillance of gonococcal infections and resistance in China.

Gonococcal AMR appears to have emerged in Asia before disseminating globally from the 1960s onward [37]. During the past decade, several ceftriaxone-resistant N. gonorrhoeae strains had been reported to be referred to F89 in France in 2010 and Spain in 2011, A8806 in Australia in 2013, GU140106 in Japan in 2014, and FC428 and FC460 in Japan in 2015 [38] since the high-level resistant strain to ceftriaxone, which was referred to as “superbugs” (H041), was firstly identified in Japan in 2009 [39]. Fortunately, these strains were considered not to have sustained transmission nationally or internationally before 2016. However, the first case with treatment failure of dual therapy of ceftriaxone and azithromycin reported in the UK in June 2016 was infected with the N. gonorrhea strain (ST6800) identical to a strain spreading in Japan that has shown reduced susceptibility to ceftriaxone and azithromycin [40], [41].Recently, the gonococcal isolates that had substantive similarity to the previously described FC428 strain in Japan have been reported from Australia, Canada and Denmark [38], [42], [43].Thus, it is likely that this strain appears to have been circulating regionally and globally. It has been recognized that a key step in management of gonococcal AMR globally should be to strengthen global collaboration in identifying gonococcal WGS from Asian strains, some of which are likely to spread to the rest of the world, which has an important implication for enriching the international genomic database by providing N. gonorrhoeae genome information of gonococcal isolates from Asia.

MLST ST1901 was the most prevalent ST among the N. gonorrhoeae isolates in many countries including in Japan and responsible for most of the decreased susceptibility to ESCs in Europe in previous studies [44], [45], [46]. However, this specific strain was not found in our study in which the most prevalent ST was ST7827 although it is not clear if the ST1901 strain was prevalent in China before this study. Although a high prevalence (around 10%) of Cfx-DS has been reported in China [15], treatment failure of ceftriaxone has not been documented in China. It is not quite clear whether no reporting of the treatment failure in China is primarily due to prevalent clone of ST7827 rather than ST1901 or the current surveillance system was not sensitive to capture the treatment failure. It is noted that close monitoring of the predominant type of ST7827 in China and other countries for determining if it will follow the pattern of other prevalent STs to be prevalent within the country or spread out of the country.

Earlier studies showed that decreased susceptibility to ESCs is primarily clonal, and the majority of ESC-DS isolates belong to only a few major lineages [47], [48].However, some of previous studies have shown heterogeneity among their gonococcal isolates with decreased antibiotic susceptibility [49], [50]. Although differences in the numbers of isolates included into the studies and the discrimination methods used for differentiate the strains across the studies may contribute to the perplexing inconsistency, further studies in a national and international collaborative manners may be important. In our study, there were specific ST predominated isolates (ST7363/lineage 2.8), harboring the penA mosaic gene, with high MIC levels to ceftriaxone (0.134 ± 0.085 mg/L). It is unclear whether this sub-lineage is a newly emerging clone to be related to decreased susceptibility or resistance to ceftriaxone. In addition, some amino acid mutations and consequently SNPs were detected in our study population to be statistically associated with Cfx-DS. However, such associations and their mechanisms should be further explored in addition to continuously monitoring the emergence and spread of strains with these mutations.

There are some potential limitations to this study. First, the incidence of gonorrhea is most likely highly underestimated in China due to the suboptimal diagnostics, screening, reporting and epidemiological surveillance performed in many settings. Although the current study was a nationwide survey, only 11 of the 31 provinces in the country participated in the study. The representativeness of these participating clinics in the 11 provinces in China is a concern. Second, as the majority of our study population were men and the majority of clinical isolates came from urethral specimens, it is unclear to what extent the specimens investigated in the current study are representative of those from women or from other anatomic sites (rectal, and pharyngeal sites) of infection in men. Third, as all the isolates were collected from symptomatic patients at clinics, it is unclear to what extent the study population is representative of those asymptomatic cases in the community. Despite these limitations, this study outlined information regarding the molecular epidemiology of gonococcal infections and resistance in China. Our work highlights the advantages and potential of WGS to be applied at scale in the current surveillance of gonococcal AMR in China. Enhanced surveillance, including genomic epidemiological studies, has been included as one of the priorities in the ROADMAP plan to address research needs for gonococcal AMR in China [51].

In conclusion, to our knowledge, we report the first nationwide study on WGS-based analysis of gonococcal infections and resistance in China as well as in Asia. The findings from this study could not only be used as the baseline data for future studies in China but also be contributable to our understanding on spread of N. gonorrhoeae and its resistant strains at regional and global levels.

Contributors

YPY, JPP, XSC, and QJ designed, initiated, and coordinated the study. YPY, and XSC were responsible for coordination of the China-GRSP and management of data gathered through programme. XQD, HPZ, WMG, BYZ, GY, NZ, LHH, WLC, ZJZ, FW and QZ coordinated the collection of information and clinical isolates in the local clinics. JPP, and SCC did the laboratory analyses. JPP, JY, BL developed the WGSA, CZ, LSX, JD, LLS and YFZ conducted the laboratory work. JPP, YPY and SCC analyzed and interpreted the data. YPY, JPP, and SCC wrote a first draft of the paper. XSC made critical revision on the first draft. All authors read, commented on, and approved the final manuscript.

Declaration of Interests

We declare no competing interests.