Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Microbiol., 21 January 2026

Sec. Antimicrobials, Resistance and Chemotherapy

Volume 16 - 2025 | https://doi.org/10.3389/fmicb.2025.1728860

This article is part of the Research TopicA Molecular and Structural Approach to Deciphering and Combating Infectious PathogensView all 10 articles

Genetic identification of mutations and MLST types associated with decreased susceptibility to ceftriaxone in Neisseria gonorrhoeae


Xin ZhangXin Zhang1Hairui WangHairui Wang1Yixin GuYixin Gu1Xiaoli Chen,Xiaoli Chen1,2Guilan Zhou,Guilan Zhou1,3Chang LiuChang Liu1Liyin JiLiyin Ji1Rui XieRui Xie1Jianzhong ZhangJianzhong Zhang1Zhujun Shao
Zhujun Shao1*Maojun Zhang
Maojun Zhang1*
  • 1National Key Laboratory of Intelligent Tracking and Forecasting for Infectious Diseases, National Institute for Communicable Disease Control and Prevention, Chinese Center for Disease Control and Prevention, Beijing, China
  • 2Department of Clinical laboratory of Shengli Oilfield Central Hospital, Dongying Key Laboratory of Cell Biology, Dongying, China
  • 3Nanjing Center for Disease Control and Prevention, Nanjing, Jiangsu Province, China

Antimicrobial resistance (AMR) in Neisseria gonorrhoeae severely limits treatment options, with increasing resistance even to first-line and last-line ceftriaxone (CRO), posing a major global public health threat. In this study, we systematically identified 53 significantly different mutations between ceftriaxone-resistant and susceptible strains in multiple proteins through bioinformatics analysis. Among these, 33 mutations were identified for the first time, notably including the PorB Q143K via structural analysis. Minimum spanning tree (MST) analysis based on these mutations marked improved sensitivity and specificity for identifying ceftriaxone-resistant strains compared to traditional sequence typing of PenA, PonA, PorB, and MtrR (68.4% vs. 53.2%; 77.3% vs. 57.5%, respectively). Furthermore, analysis of PenA sequences from global 8,325 strains (470 MLST types) revealed that mutation frequencies at key PenA sites are highly associated with MLST types, with 34 high-frequency MLST types (STs) identified. The proportions of these 34 STs were 88.38% in 611 decreased susceptibility to ceftriaxone (CRO-DS) strains and 33.09% in 8,325 background strains, respectively, revealing an extremely significant association between 34 high-frequency STs and CRO-DS (P < 0.0001). In conclusion, this work provides further insights into the molecular mechanisms of CRO resistance while offering significant value for monitoring and predicting emerging CRO-DS-associated MLST types.

Introduction

Neisseria gonorrhoeae is the causative agent of gonorrhea, a globally prevalent sexually transmitted infection (STI), with approximately 82.4 million new infections in 2020 according to the World Health Organization (WHO) (WHO, 2024). The development of AMR in N. gonorrhoeae has narrowed treatment options, making ceftriaxone monotherapy or its combination with azithromycin the first-line and last-line treatments for uncomplicated gonorrhea (Workowski and Berman, 2010). However, in 2020, to address rising AMR to azithromycin, combination therapy was replaced by ceftriaxone monotherapy (500 mg, or 1 g for people ≥ 150 kg) in the United States (Cyr et al., 2020). Alarmingly, CRO-DS isolates have surged globally (Lin et al., 2022; Fifer et al., 2024; Xiu et al., 2024), and treatment failures due to high-level resistance have occurred (Unemo et al., 2012). A study of 4,113 N. gonorrhoeae isolates from Guangdong, China, found the prevalence of CRO-DS increased from 2.05% in 2016 to 16.18% in 2019 (Lin et al., 2022), highlighting the critical importance and urgency of CRO resistance surveillance and resistance mechanism research.

Ceftriaxone resistance in N. gonorrhoeae is multifactorial, with the genes penA (penicillin-binding protein 2), ponA (penicillin-binding protein 1A), porB (trimericporin PorB.IB), and mtrR (multidrug efflux system transcriptional repressor MtrR) being the most frequently analyzed and reported contributors (Lindberg et al., 2007; Lin et al., 2021; Liao et al., 2024). A few studies also report potential links between CRO-DS and genetic locus like the mtrCDE operator (efflux transporter MtrCDE operator) (Beggs et al., 2021), blaTEM (beta-lactamase TEM) ADDIN EN.CITE (Agbodzi et al., 2023; Qin et al., 2024), rpoB (DNA-directed RNA polymerase subunit beta), rpoD (RNA polymerase sigma factor RpoD) (Palace et al., 2020), and rpoH (RNA polymerase sigma factor RpoH) (Beggs et al., 2021). PenA is the paramount determinant of CRO resistance in N. gonorrhoeae, with two distinct types: non-mosaic and mosaic alleles acquired through horizontal gene transfer (Unitt et al., 2024). For non-mosaic penA, structural analysis shows multiple mutations at A501 increase resistance to CRO. Transformation with the A501T mutation increased CRO MIC twofold, statistically confirmed (Zhu et al., 2023). Similarly, transforming mosaic penA with A501V significantly increased CRO MIC (Tomberg et al., 2010). Additionally, G542S and P551L/S are also associated to CRO-DS (Whiley et al., 2010). For mosaic penA, A311V and T483S are reported to associate with CRO-DS (Liao et al., 2024). Furthermore, A311V serves as a key marker for detecting FC428-like mosaic strains (Wang et al., 2024).

Molecular typing is critical for studying genetic evolution and antimicrobial resistance surveillance in N. gonorrhoeae. MLST (multilocus sequence typing), NG-STAR (Neisseria gonorrhoeae Sequence Typing for Antimicrobial Resistance), and NG-MAST (Neisseria gonorrhoeae multi-antigen sequence typing) are the three primary typing schemes (Smolarczyk et al., 2021). MLST compares genetic variations in seven conserved housekeeping genes (abcZ, adk, fumC, aroE, pdhC, gdh, pgm) and exhibits high discriminatory power for global isolates (Bennett et al., 2007). The NG-STAR website (https://ngstar.canada.ca) offers a standardized method of classifying seven well-characterized genes (penA, mtrR, porB, ponA, gyrA, parC and 23S rRNA) associated antimicrobial resistance to three classes of antibiotics (cephalosporins, macrolides and fluoroquinolones) (Demczuk et al., 2017). NG-MAST is a highly discriminatory molecular typing method targeting a 490-bp porB fragment and a 390-bp tbpB (transferrin-binding protein subunit B) fragment (Bilek et al., 2007).

An increasing number of studies, especially in Asian regions, have associated CRO-DS strains with specific MLST types (STs). In Shenzhen, China (2014–2018), ST7365 (21.43%, 12/56) and ST7360 (18.18%, 10/55) significantly correlated with CRO-DS (Li et al., 2021). Across 20 cities in Guangdong, China (2021), CRO-DS was mainly associated with ST7363 (16%, 8/50), ST1903 (14%, 7/50), ST1901 (12%, 6/50), and ST7365 (10%, 5/50) (Liao et al., 2023). Among 463 isolates from Guangdong, Jiangxi, and Shanghai, China (2021), an 8.9% CRO-DS rate was detected, highly linked to ST7365 and ST1903 (Tang et al., 2023). In Hanoi, Vietnam (2023–2024), the CRO resistance rate was 27% (90/352 isolates), with ST1901 predominating (71.1%, 64/90) (Laumen et al., 2025). However, the reason why these specific STs are consistently linked to CRO-DS remains unclear. It appears not to be due to their high prevalence. For example, ST9363, one of the most globally prevalent STs, is rarely linked to CRO-DS, with reports mainly related to resistance to agents like azithromycin (Melendez et al., 2024). Therefore, this study systematically analyzed the genetic determinants of CRO resistance, MLST types, and CRO-DS in N. gonorrhoeae to investigate the molecular mechanisms of CRO resistance and identify STs associated with CRO-DS.

Materials and methods

Strain collection and screening

As of November 18, 2024, a total of 23,381 whole-genome sequences of N. gonorrhoeae isolates from 2005 to 2024 were downloaded from PubMLST database (https://pubmlst.org/) (Jolley et al., 2018). The European Committee on Antimicrobial Susceptibility Testing (EUCAST) defines an MIC of ≤ 0.125 mg/L as susceptible and an MIC of > 0.125 mg/L as resistant to ceftriaxone in N. gonorrhoeae (EUCAST, 2025). Among these, 94 isolates were CRO-R strains, 14,118 isolates were CRO-S strains, and 9,169 isolates had no available CRO resistance information.

Genomes used for the analysis of genetic determinants of CRO resistance were required to have ≤ 300 scaffolds and ≤ 1,000 N bases, resulting in the selection of 79 CRO-R strains. For CRO-S strains, 10 isolates were selected annually, preferably from different countries, to ensure geographic diversity. Due to the unavailability of CRO resistance data in 2023 and 2024, a total of 180 CRO-S strains were randomly selected. Including the susceptible reference strain FA1090, a total of 260 strains were selected (For detailed background information, see Supplementary Table 1). The genomes used for analyzing the correlation between PenA mutations and MLST types were selected from the 9,169 isolates not used in the previous analysis, all lacking cephalosporin resistance data. From these, we chose a total of 8,325 isolates with complete MLST information (for detailed background information, see Supplementary Table 2).

All the genomes of isolates were validated with the species identification module (taxonomy_wf) in Genome Database Taxonomy toolkit (GTDB-Tk) (Parks et al., 2022). All isolates were compared with a representative N. gonorrhoeae genome, and the average nucleotide identity (ANI) values were calculated for each isolate. A threshold of 0.95 ANI was applied to confirm that the genomes belonged to the N. gonorrhoeae species. The quality of the genomes was assessed using CheckM software v1.2.2 (Parks et al., 2015). To ensure the reliability of the genomes used for analyzing mutation sites, we applied stringent quality control criteria. Only the isolates with a completeness of ≥ 90% and contamination of ≤ 10% were considered suitable for further analysis.

Phylogenetic and pan-genome analysis

Core Single Nucleotide Polymorphism (SNP) calling was conducted using the Snippy 4.6.0 (https://github.com/tseemann/snippy), utilizing FA1090 (GenBank accession no. AE004969.1) as a reference strain. The Gubbins 2.4.1 was employed as the recombination-removal tool, enabling the extraction of pure SNP data by eliminating recombination events. The genomes underwent gene prediction and functional annotation using the Prokka pipeline v1.14.6 (Seemann, 2014). Core-pan genome analysis was conducted using the Roary pipeline 3.13.0 (Page et al., 2015), which utilizes annotated assemblies in GFF3 format obtained from Prokka results. The MAFFT 7.471 was then employed to perform multiple sequence alignment of these sequences. Subsequently, the phylogenomic tree was constructed using FastTree 2.1.10 (Price et al., 2009) with the maximum-likelihood (ML) algorithm. The phylogenomic tree was visualized using iTOL v7.1.1 (Letunic and Bork, 2024).

Analysis of genetic determinants of ceftriaxone resistance

Based on the reference sequence of each gene, a python script was used to extract the full-length amino acid sequences of the PenA, PonA, PorB, MtrR, blaTEM, RpoB, RpoD, RpoH, and a 27mer of mtrCDE operator, from 260 strains via BLAST alignment. Protein structure analysis was performed to identify key regions of each gene. AlphaFold 3 was used for 3D structure prediction with consistent seed values across all predictions to ensure comparability (Abramson et al., 2024). The crystal structure of the PenA-CRO complex was obtained from the PDB database (ID: 6XQV) (Fenton et al., 2021). Active pockets were predicted using POCASA 1.1 (Yu et al., 2010), and surface hydrophobicity was analyzed with the Python script color_h.py (Eisenberg et al., 1984). Structural alignment and visualization were carried out using PyMOL 2.5.2. WebLogo-3.7.9 was used to perform amino acid conservation analysis (Crooks et al., 2004). DeepPBS (Mitra et al., 2024) was used for geometric deep learning of protein-DNA binding specificity, with the parameter “Calculate heavy atom relative importance (RI) scores” enabled.

The 260 strains were divided into CRO-R (n = 79) and CRO-S (n = 181) groups. A python script was used to calculate the amino acid or nucleotide frequency at each site for each gene in both groups. For the penA, mosaicism was defined by alterations in amino acids 375–377, with mosaic and non-mosaic variants analyzed separately (Deng et al., 2019). Then, to identify mutation sites associated with CRO resistance or susceptibility, a two-tailed Fisher's exact test was performed for each differential site by R 4.4.1. A P value < 0.05 was considered statistically significant. Corrected P value (q value) was obtained using the Benjamini–Hochberg false discovery rate (FDR) approach. The Cleveland dot plot was generated using the R package tidyverse.

Molecular typing

MLST, NG-STAR and NG-MAST (v2.0) were defined by using the PubMLST database. The cgMLST scheme of 260 strains was prepared using chewBBACA v2.8.5 (Silva et al., 2018). A training file generated by Prodigal v2.6.3 using the genome FA1090. The RemoveGenes module was used to delete the detected 39 paralogous loci from the 2,741 loci. The TestGenomeQuality module was used to verify the reliability of the dataset quality for 260 genomes. The ExtractCgMLST module was run to determine the set of loci in the core genome for the loci presence thresholds of 99%. Core genome composed of 1361/2702 genes. All minimum spanning trees used for molecular typing were constructed using the MSTreeV2 algorithm and visualized using GrapeTree v2.1 (Windows version) (Zhou et al., 2018). Statistical analysis of mutation sites was performed using the Chi-square or Fisher's exact test, with significance assessed by the false discovery rate (FDR) and defined by q values: P q < 0.05 as statistically significant, q < 0.01 as highly statistically significant, and q < 0.001 as extremely statistically significant.

Analysis of the correlation between mutation frequencies at key PenA sites and MLST types

Eight thousand three hundred and twenty five genomes were annotated using Prokka, and sequences named “penicillin-binding protein 2” were exported from the faa files, resulting in 7678 non-mosaic penA alleles and 647 mosaic penA alleles. A cluster heatmap of PenA mutations sites and MLST types was generated using the R package pheatmap. A neighbor-joining tree was constructed using the FastME V2.0 algorithm in the GrapeTree software. The 647 mosaic penA alleles were assigned identification numbers through the pubMLST website. The nucleotide sequences of the 31 obtained types were then used to construct a maximum-likelihood tree using FastTree. A Sankey diagram of MLST, NG-STAR, and NG-MAST was created using the R package networkD3.

Results

Phylogenetic and pan-genome analysis of 260 strains employed for genetic determinants of ceftriaxone resistance analysis

Phylogenetic analysis of the 260 strains was conducted to trace evolutionary traits. Continent, country, gender, source, year, penA alleles, MLST, NG-STAR, NG-MAST, CRO-SIR (Susceptible, Intermediate, Resistant), CRO-MIC, CFM-SIR, CFM-MIC, and porB type are shown in Supplementary Figure 1. The 260 strains from 35 countries across six continents comprised 79 STs (ST1901, ST7363, ST9363, and ST1903 with ≥ 10 isolates each), 141 NG-STAR types, and 154 NG-MAST types.

Core-pan genome analysis was used to explore the genetic diversity of 260 strains, defining core genes as those present in at least 99% of the strains with ≥ 95% identity. The pan-genome contained a total of 5,823 genes, including 1,481 core genes. According to the topology of the phylogenomic tree based on core genes, these strains can be divided into six clades, named Clade 1 to Clade 6. Interestingly, most ceftriaxone-resistant strains (67/79, 84.81%) and penA mosaic strains (74/79, 93.67%) clustered in Clades 3, 4, and 5 (Figure 1A). In combination with the phylogenomic tree constructed from core homologous clusters, these six clades exhibited significantly different gene family structures from each other, especially in the gene clusters labeled in red as A and B (Figure 1B). Clade 3 tends to possess both clusters, Clade 4 tends to have A but not B, whereas Clade 6 tends to lack both.

Figure 1
Two-part graphical representation of bacterial clades and genetic data. Panel A shows a radial phylogenetic tree dividing strains into six clades, with markers indicating ceftriaxone resistance (R) and sensitivity (S), and penA allele types (mosaic, non-mosaic). Panel B presents a dendrogram of 260 strains aligned with a gene presence-absence matrix of 5823 clusters, highlighting specific gene clusters A and B. Clade colors correspond between panels.

Figure 1. Phylogenetic and pan-genome analysis. (A) The phylogenetic tree based on core genes. The first layer annotation represents ceftriaxone SIR (Susceptible, Intermediate, Resistant), and the second layer annotation represents penA mosaic types. (B) Pan-genome matrix. Blue indicates presence and white indicates absence. Clade color coding applies to both (A) and (B). Clades 3, 4, and 5 are represented in shades of orange, while Clades 1, 2, and 6 are shown in shades of green.

Further characterization identified a total of 73 genes (a = 47, b = 26) in the two clusters. The A cluster includes topA_2 (type I DNA topoisomerase), traC (type-IV secretion protein), parA (chromosome partitioning protein), dsbC (DsbC family protein), and yhaV (Ribonuclease toxin). The B cluster contains hin (DNA-invertase), traG (conjugal transfer protein), ssbB (single-stranded DNA-binding protein), trfA (plasmid replication initiator), trbM (Dutch type tetM conjugative plasmid), and traJ (conjugal transfer regulator). The remaining 62 genes are annotated as hypothetical proteins. Among the 11 identified genes, traC, parA, dsbC, traG, and ssbB are located on the horizontally transferred gonococcal genetic island (GGI, containing approximately 62 CDSs), and traC, parA, dsbC, and traG belong to the type IV secretion system genes encoded by GGI. In addition, the A cluster has a GC content of 45.1%, close to the reported GGI value of 44% (Hamilton et al., 2005), suggesting it contains more GGI-related genes. The B cluster has a GC content of 50.6%, close to the N. gonorrhoeae genome average of 51% (Dillard et al., 2011), indicating relatively fewer GGI-related genes.

Analysis of genetic determinants of ceftriaxone resistance

Based on literature review, we identified 8 proteins and one operator potentially related to CRO-R in N. gonorrhoeae. We conducted structural analysis to pinpoint the regions or sites most impactful for CRO resistance. The structure of PenA is divided into two parts, with Q241-L564 being the region directly interacting with CRO, thus requiring special focus, especially around the active site S310 (Supplementary Figure 2A). Further interaction analysis identified 11 amino acids forming hydrogen bonds with CRO: E307, A310 (S310), S362, N364, Y422, T483, T498, T500, R502, G546, and V547 (Supplementary Figure 2E). Notably, there are 33 amino acids within 6Å of CRO: E307, P308, G309, A310 (S310), A311, K313, R345, T347, H348, K361, S362, N364, T417, F420, G421, Y422, T483, K497, T498, G499, T500, A501, R502, K503, Y509, H514, T517, Y543, Y544, G545, G546, V547, V548. Additionally, PenA includes non-mosaic and mosaic types, differing by 40-50 amino acids in the Q241-L564 region. Their structures, active pocket locations (Supplementary Figure 2B), surface charges, and hydrophobicity (Supplementary Figures 2C, D) also differ significantly. Therefore, mutation sites in mosaic and non-mosaic PenA must be analyzed separately.

The standard NG-STAR typing for PorB only covers 10 amino acids (positions 117–126, Supplementary Figure 2F, purple). However, structural analysis reveals a core variable region of 49 amino acids in the PorB channel (Supplementary Figure 2F, yellow + purple; sequences and conservation in Supplementary Figure 2G), with significant structural differences between PorB1b and PorB1a, the latter associated with disseminated gonococcal infection (Cartee et al., 2022). We therefore recommend expanding PorB1b analysis to include these 49 amino acids. Meanwhile, most studies focus on MtrR and the A deletion in its promoter's−35 region, but the mtrCDE operator region (including the−35 region), which directly binds the inhibitory protein, also warrants attention. DeepPBS was used to analyze the binding specificity between MtrR and the mtrCDE operator, identifying 10 directly interacting amino acids and multiple key specificity-determining nucleotide positions (Supplementary Figure 2H).

The 260 strains were divided into ceftriaxone-resistant (CRO-R, n = 79) and ceftriaxone-susceptible (CRO-S, n = 181) groups. Amino acid frequencies in the “key regions” of each gene were statistically analyzed to assess the association between mutation sites and CRO resistance or susceptibility. Statistical analysis identified 53 distinct mutations with significant differences. These were considered potential genetic determinants of CRO resistance or susceptibility (Table 1). Cleveland dot plots were used to show the frequency differences of these mutations between resistant and susceptible groups, with newly identified mutations highlighted in red (Figure 2), including the following: mosaic PenA V279A, E285D, K288R, Q291R, A328P, S341P, T485I, A549T, P552V, K555Q, I556V, PorB A109V, A121S, N122K, F131Y, N134E, V135F, V135L, G140K, Q143K, Q143E, Q143G, Q143R, V151A, MtrR A39T, A40D, R44H, D79N, S183N, M197I, mtrCDE operator G15A, RpoB E1175D, RpoH A184T.

Table 1
www.frontiersin.org

Table 1. Potential genetic determinants of ceftriaxone resistance or susceptibility.

Figure 2
Scatter plot showing mutation sites and their frequency in resistant (red) and susceptible (blue) groups. The x-axis represents frequency percentages, while the y-axis lists mutation sites. Different gene categories are labeled: penA (non-mosaic and mosaic), ponA, porB, mtrR, mtrCDE operator, rpoB, rpoD, and rpoH.

Figure 2. The frequency difference analysis of mutation sites. Resistant group in red, susceptible group in blue. Line length corresponds to the magnitude of intergroup difference.

Molecular mechanism analysis of ceftriaxone resistance in Neisseria gonorrhoeae through molecular typing

To further investigate the molecular mechanisms of CRO resistance, we conducted a comparative analysis of various sequence typing methods. First, we compared MLST, NG-STAR, and cgMLST (core genome multilocus sequence) methods (Figures 3AC). As expected, NG-STAR (based on 7 resistance genes) and cgMLST (using 1,361 core genes) outperformed MLST in differentiating resistant and susceptible strains. However, the numerous genes and countless sites obscure which ones are critical for CRO resistance.

Figure 3
Six phylogenetic trees and one graph illustrate ceftriaxone susceptibility, denoted as S (green) and resistance, denoted as R (orange) across various groups A, B, and C. Each tree (A-F) visualizes relationships among samples with varying sizes of circular nodes representing the proportions of susceptibility and resistance. The graph (G) indicates mutation sites affecting susceptibility with markers for statistical significance, showing Group A (red), Group B (blue), and Group C (gray) along the x-axis.

Figure 3. Molecular typing and mutation site difference analysis associated with ceftriaxone resistance. Orange indicates ceftriaxone resistance, and green indicates susceptibility. (A) MLST of the 7 housekeeping genes. (B) NG-STAR of the 7 resistance genes. (C) cgMLST of the 1,361 core genes. (D) Sequence typing of 20 sites in PenA. (E) Sequence typing of 14 common sites in PenA, PonA, PorB, and MtrR. (F) Sequence typing of 49 sites in PenA, PonA, PorB, MtrR, mtrCDE operator, blaTEM, RpoB, RpoD, and RpoH. (G) Comparison of the percentage of different mutation sites among the three groups in Figure (F). For each mutation site, red, blue, and gray dots represent its percentage in Group A, B, and C, respectively. Lines between dots show comparisons between groups, with significance levels denoted as follows: ns (not significant, q ≥ 0.05), * (q < 0.05), ** (q < 0.01), and *** (q < 0.001).

Therefore, to identify which sites are crucial for CRO resistance, we compared the sequence typing of mutation sites across multiple genes. As shown in PenA sequence typing (Figure 3D), the sensitivity (percentage of resistant strains forming distinct branches) was only 36.7%, and the specificity (percentage of sensitive strains forming distinct branches) was only 12.7%. Therefore, PenA alone is insufficient for distinguishing CRO-R strains. In contrast, the combined sequence typing of PenA, PonA, PorB, and MtrR (Figure 3E) showed a sensitivity of 53.2% and a specificity of 57.5%, representing a clear improvement. Notably, sequence typing based on the 49 sites (Figure 3F, details in Supplementary Table 1) achieved a sensitivity of 68.4% and specificity of 77.3%, showing further improvement. Additionally, a large majority of non-mosaic resistant (Group A, red circle) and sensitive (Group C) strains are clustered by this method, with all mosaic resistant strains uniquely assigned to Group B (blue circle). Resistant strains in Groups A and B account for 89.9% of all resistant strains, demonstrating strong clustering efficacy.

The percentages of the previously identified resistance-associated mutations were compared among the three groups. The mutation sites showing significant differences, including PenA A501V/T, F504L, A510V, A516G, G542S, P551L/S; PonA L421P; PorB G120K, A121D/G, Q143K; MtrR H105Y; the−35A Del in the mtrR promoter; RpoD I229V; and RpoH A184T, are displayed in Figure 3G. These mutations effectively distinguished the clusters of resistant strains (Groups A and B) from sensitive strains (Group C), highlighting their critical role in differentiating CRO resistance and susceptibility. Unlike other mutations, PorB A121D/G and the first-identified PorB Q143K (red font) also show significant differences between Groups A and B, suggesting distinct mutation profiles in non-mosaic and mosaic resistant strains. Overall, these findings provide further insight into the molecular mechanisms of ceftriaxone resistance in N. gonorrhoeae.

Analysis of the correlation between mutation frequencies at key PenA sites and MLST types

Based on the above analysis, 20 distinct mutation sites in PenA were identified. Among them, H541N was associated with susceptibility, while F504L, A510V, and A516G were commonly mutated in both resistant and susceptible strains, and thus excluded. The final key PenA mutation sites were identified as A501V/T, G542S, and P551L/S in non-mosaic alleles, along with 13 sites in mosaic alleles. The mutation frequency at key PenA sites and the occurrence of mosaic penA alleles were statistically analyzed for each MLST type. The mutation frequency coefficient (k = 1) was defined as the overall mutation rate at each site, calculated by dividing the total number of mutations at each site by the total number of strains (counting mosaic and non-mosaic strains separately). The threshold of k = 2 was set for identifying high-frequency STs.

Among 8,325 strains (non-mosaic = 7,678; mosaic = 647), 129 STs exceeding the threshold were identified as high-frequency STs. A total of 34 MLST types with ≥ 5 isolates were defined as statistically significant high-frequency STs. 33 STs were identified from non-mosaic alleles, including ST7363, ST7822, ST1579, ST1583, ST1901, ST10314, ST7827, ST1893, ST11981, ST11706, ST7371, ST1600, ST1903, ST7823, ST13489, ST1902, ST7365, ST7367, ST8123, ST10313, ST11230, ST7360, ST14422, ST1927, ST13143, ST12982, ST1920, ST8153, ST16414, ST1590, ST6714, ST11231, and ST9902. The high-frequency STs corresponding to mosaic sites must also belong to the high-frequency mosaic types, resulting in the identification of only two types: ST7363 and ST13734.

Clustering analysis of mutation frequencies at key PenA sites was performed for all MLST types with ≥ 5 isolates (Figure 4). For ST1902 and ST9902, all three non-mosaic loci are high-frequency mutations. For ST1901, ST8123, and ST7360, two non-mosaic loci are high-frequency mutations, and their mosaic frequencies exceed the threshold. ST7363 is the only MLST type that belongs to both high-frequency non-mosaic and mosaic groups. Detailed mutation sites, frequencies, and the identifiers of the 129 STs are provided in Supplementary Table 3.

Figure 4
Clustered heatmap displaying mutation presence in different strains, labeled on the right. Columns represent various mutations such as A501V/VT, G542S, P551L/S, mosaic, and mosaic-mutation. Intensity ranges from light orange for lower values to dark orange for higher values, indicating mutation frequency, with a scale from 0 to 1.

Figure 4. Clustering heatmap of PenA Mutations and MLST types. A gradient color scale from white (0) to orange (1) indicates the mutation frequency values across the range of 0 to 1. Use red font for high-frequency MLST types, and display the exact frequency values that exceed the thresholds.

A neighbor-joining tree based on MLST was constructed to show the evolutionary relationships among STs (Figure 5). This analysis systematically illustrated the distribution patterns of mutation frequencies at key PenA sites across different MLST types. Strains with key mutation sites were classified as “non-mosaic-high” or “mosaic-high,” while the others were categorized as “non-mosaic-low” or “mosaic-low.” From this, we can observe that the mutation frequencies at key PenA sites are highly associated with MLST types. 34 high-frequency STs are labeled in red and highlighted in yellow, while some types show very low or even zero mutation frequency despite having large strain numbers (n > 100, labeled in black).

Figure 5
Phylogenetic tree diagram showing relationships between various strains, represented as circles. Circle size indicates prevalence; colors represent MLST categories: non-mosaic-low, non-mosaic-high, mosaic-low, and mosaic-high. A legend in the bottom right clarifies colors and category numbers. Lines connect related strains, highlighting genetic relationships.

Figure 5. Neighbor-joining tree of MLST for 8,325 isolates. Each circle represents an MLST type, with size indicating the number of isolates. High-frequency types are displayed in orange (non-mosaic) and pink (mosaic), while low-frequency types are shown in dark blue (non-mosaic) and light blue (mosaic).

Not all mosaic strains are CRO-R (Mlynarczyk-Bonikowska et al., 2022), but specific mosaic penA alleles, particularly in FC428-like strains, are strongly associated with CRO-R (Xiu et al., 2025). Therefore, for 647 penA mosaic strains, we must not only analyze the relationship between mutation sites and MLST types but also investigate the link between mutation sites and penA alleles. The phylogenetic tree based on 31 mosaic penA alleles from 647 isolates divided them into two groups (Figure 6). One group contained 1 to multiple resistance-associated mutations, with penA-60.001 harboring all 13 mutations and being significantly associated with CRO-R (Adamson et al., 2024). The other group lacked resistance-associated mutations, such as the most prevalent allele penA-34.001. This suggests that the 13 previously identified key mosaic mutation sites are highly correlated with mosaic penA alleles and CRO-R.

Figure 6
Phylogenetic tree and heatmap showing mutations in penicillin-binding protein 2 strains. The tree on the left displays evolutionary relationships. The heatmap on the right indicates mutation sites across strains, color-coded for specific mutations: V279A, E285D, K288R/Q, Q291R, A311V, A328P/L, S341P, T483S, T485I, A549T, P552V, K555Q, I566V. The number of strains is illustrated by different shades of blue to purple.

Figure 6. The phylogenetic analysis of 31 mosaic penA alleles from 647 isolates. Midpoint-rooted for visualization only. Mosaic penA alleles linked with resistance-associated mutations are shown in red font. The number of isolates is indicated by a gradient from light blue to purple. Pre-mutation amino acids are colored green, while post-mutation amino acids are colored orange or pink.

To evaluate the utility of the 34 high-frequency STs in ceftriaxone resistance surveillance, we validated them using 611 CRO-DS strains (MIC > 0.064 mg/L; see Supplementary Table 2 for background information) (Lindberg et al., 2007; Whiley et al., 2010; Lin et al., 2021) from all 23,381 isolates. As shown in Table 2, 88.38% (540/611; 95% CI: 85.85–90.91%) of CRO-DS strains belonged to 34 high-frequency STs, compared to 33.09% (2,755/8,325; 95% CI: 32.09–34.09%) of background strains. A chi-square test was used to compare the distribution of the 34 high-frequency STs between the two groups. The results showed that CRO-DS strains were significantly enriched in the 34 high-frequency STs (χ2 = 747.48, df = 1, P < 0.0001). Relative risk analysis demonstrated that CRO-DS strains were 2.67 times more likely to occur in high-frequency STs than background strains (RR = 2.67, 95% CI: 2.56–2.79). In summary, the 34 high-frequency STs demonstrated an extremely significant positive correlation with CRO-DS (P < 0.0001). 88.38% of CRO-DS strains were concentrated in the 34 high-frequency STs, which constitute only 7.23% (34/470) of the MLST types in the background strains. This indicates that the 34 high-frequency STs are strong predictors of CRO-DS risk. Moreover, strains belonging to the 34 high-frequency STs account for 33.09% of background strains, suggesting that approximately one-third of N. gonorrhoeae may carry a notably higher risk of CRO-DS. Meanwhile, the major NG-STAR and NG-MAST types corresponding to the 34 high-frequency STs are shown in Supplementary Figure 3, including NG-STAR ST90 (Demczuk et al., 2017) and NG-MAST ST1407 (Gianecini et al., 2019), both reported in association with CRO-DS. Therefore, continuous monitoring of 34 high-frequency STs and their corresponding NG-STAR and NG-MAST types is of great significance.

Table 2
www.frontiersin.org

Table 2. Distribution of 34 high-frequency MLST types among CRO-DS strains and background strains.

Discussion

Given the increasing global reports of CRO-DS isolates (Lin et al., 2022; Fifer et al., 2024; Xiu et al., 2024), the importance of researching first-line and last-line ceftriaxone resistance mechanisms and enhancing surveillance cannot be overstated. This study utilized genomic and resistance data from 23,381 isolates (as of November 18, 2024) in pubMLST, the largest global N. gonorrhoeae database. CRO-R and CRO-S strains were used to investigate the molecular mechanisms of CRO resistance. Strains lacking cephalosporin resistance data (forming a distinct dataset from the CRO-DS strains used for validation) were used to analyze relationships between mutation frequencies at key PenA sites, MLST types, and CRO-DS.

Pan-genome analysis identified two gene clusters, A and B, associated with GGI. The GGI secretes DNA into the environment via a type IV secretion system (T4SS) independent of direct cell-cell contact (Hamilton et al., 2005), while N. gonorrhoeae can acquire environmental DNA through natural transformation. This dual mechanism drives the rapid spread of antibiotic resistance. Notably, the GGI has been reported to significantly correlate with CRO-R (Harrison et al., 2016). These findings further suggest potential links between the A/B clusters and CRO-R, underscoring the importance of pan-genome analysis in studying gonococcal antibiotic resistance.

This study systematically identified 53 mutation sites with significant differences in genes associated with ceftriaxone resistance (CRO-R) and susceptibility (CRO-S) in N. gonorrhoeae through literature review, statistical analysis, and protein structure/interaction assessments. Among these, 33 mutations were newly reported, with the Q143K mutation in PorB (expanded from 10aa to 49aa by structural analysis) being highlighted for the first time as significantly linked to ceftriaxone resistance. Based on these mutation sites, we constructed a minimum spanning tree (Figure 3F) to investigate the molecular mechanisms of CRO resistance, which effectively segregated CRO-R and CRO-S strains. However, this typing method still has limitations. First, the limited number of mosaic strains (n = 79) and mosaic alleles (n = 12) may lead to incomplete identification of key mosaic mutation sites. Second, the largest branch of mosaic strains includes 19 susceptible strains (penA-34.001) and 16 resistant isolates (penA-34.001, n = 9; penA-42.001, n = 4; penA-34.009, n = 2; penA-34.01, n = 1), which could not be distinguished by these mutation site combinations, limiting further improvements in sensitivity and specificity. This suggests that there may be other, yet unidentified resistance mechanisms. Since the mechanisms of ceftriaxone resistance in N. gonorrhoeae involve a complex, multi-gene determined process, further rationally designed experiments are therefore needed for phenotypic validation. Nevertheless, these findings provide new insights for identifying CRO-R strains and understanding the molecular mechanisms of CRO resistance.

ST7363 and ST1901 are the most widely reported STs associated with CRO-DS globally (Thomas et al., 2021), likely due to their high global prevalence and the high-frequency mutations/accumulation at key PenA sites coupled with a remarkably high mosaic frequency. ST7363 exhibits 83% mutation frequency at non-mosaic sites A501V/T and P551L/S, with mosaic frequency and mosaic site mutation frequency at 34% and 45%, respectively. ST1901 shows 46% and 47% mutation frequencies at non-mosaic sites G542S and P551L/S, alongside a mosaic frequency of 62%. Although CRO resistance in N. gonorrhoeae is a complex multigenic process, analysis of mutation frequencies at key PenA sites still partially explains the molecular mechanisms underlying CRO resistance in ST7363 and ST1901.

Our study revealed for the first time that mutation frequencies at key PenA sites are highly correlated with MLST types, with substantial variation observed across different MLST types. Validation using 611 CRO-DS strains confirms that MLST types with high-frequency mutations at these key sites are more prone to develop CRO-DS strains. We know that different MLST types represent distinct genetic lineages. Therefore, we hypothesize that high-frequency and low-frequency mutation groups at key PenA sites may differ in their DNA replication and repair mechanisms, capacity to take up or transform exogenous DNA, or the local sequence context of the penA and other resistance-related genes. Studies have shown that local sequence context strongly impacts DNA replication accuracy in vivo (Schroeder et al., 2016). This may represent an important avenue for future research into the mechanisms underlying why specific MLST types are more prone to developing CRO-DS strains. However, considerable experimentation is still required to validate these hypotheses.

Our findings are strongly supported by previous reports. ST7363, ST1901, ST1903, ST7365, and ST7360 are widely reported to be associated with CRO-DS (Li et al., 2021; Liao et al., 2023; Tang et al., 2023; Laumen et al., 2025). Notably, in Shenzhen, China, ST11231, ST7365, and ST8123 exhibited notable CRO-R rates of 35.71%, 28.30%, and 13.04% during 2019–2020, while these STs exhibited 100% susceptibility to CRO in 2014–2018. The authors defined ST11231, with the highest CRO-R rate, as a key target for future surveillance (Wang et al., 2023). Separately, 8 of 10 CRO-R strains detected in the UK during 2021–2022 belonged to ST8123 (Day et al., 2022). All reported STs fall within the 34 high-frequency STs identified in our study. Additionally, among the 95 types exceeding the threshold but excluded from the 34 high-frequency STs, some, like ST16406 (Ouk et al., 2023; Fifer et al., 2024) and ST13871 (Trinh et al., 2022), are still widely reported to be associated with CRO-R. Therefore, it is necessary to increase the number of strains to verify whether these types belong to high-frequency STs. We recommend prioritizing these 34 STs for monitoring CRO resistance in N. gonorrhoeae and encourage researchers to adopt similar approaches using diverse isolates, along with more mutation sites, to identify additional high-frequency STs.

Data availability statement

The genomic data used in this study were all obtained from the PubMLST database (https://pubmlst.org/bigsdb?db=pubmlst_neisseria_isolates&page=query&genomes=1). Genomic background information and additional data supporting the findings of this study are available in the text and Supplementary material.

Author contributions

XZ: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. HW: Writing – review & editing. YG: Writing – review & editing. XC: Writing – review & editing. GZ: Writing – review & editing. CL: Writing – review & editing. LJ: Writing – review & editing. RX: Writing – review & editing. JZ: Writing – review & editing. ZS: Writing – review & editing. MZ: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – review & editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work was supported by Project for Gonorrhea Vaccine Development (33076), Study on Prevention and Diagnosis of Helicobacter pylori Infection in Tropical Area (2022KJCX62) and Study on Treatment Technique and Prevention Strategy of Helicobacter pylori Infection in Hainan Area (ZDYF2024SHFZ049).

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2025.1728860/full#supplementary-material

Supplementary Figure 1 | The phylogenetic tree of 260 isolates based on core SNPs.

Supplementary Figure 2 | Structural analysis of PenA, PorB, and MtrR.

Supplementary Figure 3 | The correspondence between MLST, NG-STAR, and NG-MAST.

Supplementary Table 1 | Background information and mutation sites of 260 strains.

Supplementary Table 2 | Background information of 8,325 background strains and 611 CRO-DS strains.

Supplementary Table 3 | Mutation sites and frequencies in PenA from 8,325 strains and identifiers of 129 MLST types.

References

Abramson, J., Adler, J., Dunger, J., Evans, R., Green, T., Pritzel, A., et al. (2024). Accurate structure prediction of biomolecular interactions with AlphaFold 3. Nature 630, 493–500. doi: 10.1038/s41586-024-07487-w

PubMed Abstract | Crossref Full Text | Google Scholar

Adamson, P. C., Hieu, V. N., Nhung, P. H., Whiley, D. M., and Chau, T. M. (2024). Ceftriaxone resistance in Neisseria gonorrhoeae associated with the penA-60.001 allele in Hanoi, Viet Nam. Lancet Infect. Dis. 24, e351–e352. doi: 10.1016/S1473-3099(24)00230-5

PubMed Abstract | Crossref Full Text | Google Scholar

Agbodzi, B., Duodu, S., Dela, H., Kumordjie, S., Yeboah, C., Behene, E., et al. (2023). Whole genome analysis and antimicrobial resistance of Neisseria gonorrhoeae isolates from Ghana. Front. Microbiol. 14:1163450. doi: 10.3389/fmicb.2023.1163450

PubMed Abstract | Crossref Full Text | Google Scholar

Beggs, G. A., Ayala, J. C., Kavanaugh, L. G., Read, T. D., Hooks, G. M., Schumacher, M. A., et al. (2021). Structures of Neisseria gonorrhoeae MtrR-operator complexes reveal molecular mechanisms of DNA recognition and antibiotic resistance-conferring clinical mutations. Nucleic Acids Res. 49, 4155–4170. doi: 10.1093/nar/gkab213

PubMed Abstract | Crossref Full Text | Google Scholar

Bennett, J. S., Jolley, K. A., Sparling, P. F., Saunders, N. J., Hart, C. A., Feavers, I. M., et al. (2007). Species status of Neisseria gonorrhoeae: evolutionary and epidemiological inferences from multilocus sequence typing. BMC Biol. 5:35. doi: 10.1186/1741-7007-5-35

PubMed Abstract | Crossref Full Text | Google Scholar

Bilek, N., Martin, I. M., Bell, G., Kinghorn, G. R., Ison, C. A., and Spratt, B. G. (2007). Concordance between Neisseria gonorrhoeae genotypes recovered from known sexual contacts. J. Clin. Microbiol. 45, 3564–3567. doi: 10.1128/JCM.01453-07

PubMed Abstract | Crossref Full Text | Google Scholar

Cartee, J. C., Joseph, S. J., Weston, E., Pham, C. D., Thomas, J. C., Schlanger, K., et al. (2022). Phylogenomic comparison of Neisseria gonorrhoeae causing disseminated gonococcal infections and uncomplicated gonorrhea in Georgia, United States. Open Forum Infect. Dis. 9:ofac247. doi: 10.1093/ofid/ofac247

PubMed Abstract | Crossref Full Text | Google Scholar

Crooks, G. E., Hon, G., Chandonia, J. M., and Brenner, S. E. (2004). WebLogo: a sequence logo generator. Genome Res. 14, 1188–1190. doi: 10.1101/gr.849004

PubMed Abstract | Crossref Full Text | Google Scholar

Cyr, S. S., Barbee, L., Workowski, K. A., Bachmann, L. H., Pham, C., Schlanger, K., et al. (2020). Update to CDC's Treatment Guidelines for Gonococcal Infection, 2020. Mortal. Wkly. Rep. 69, 1911–1916. doi: 10.15585/mmwr.mm6950a6

PubMed Abstract | Crossref Full Text | Google Scholar

Day, M., Pitt, R., Mody, N., Saunders, J., Rai, R., Nori, A., et al. (2022). Detection of 10 cases of ceftriaxone-resistant Neisseria gonorrhoeae in the United Kingdom, December 2021 to June 2022. Eurosurveillance 27:2200803. doi: 10.2807/1560-7917.ES.2022.27.46.2200803

PubMed Abstract | Crossref Full Text | Google Scholar

Demczuk, W., Sidhu, S., Unemo, M., Whiley, D. M., Allen, V. G., Dillon, J. R., et al. (2017). Neisseria gonorrhoeae sequence typing for antimicrobial resistance, a novel antimicrobial resistance multilocus typing scheme for tracking global dissemination of N. gonorrhoeae strains. J. Clin. Microbiol. 55, 1454–1468. doi: 10.1128/JCM.00100-17

Crossref Full Text | Google Scholar

Deng, X., Allan-Blitz, L. T., and Klausner, J. D. (2019). Using the genetic characteristics of Neisseria gonorrhoeae strains with decreased susceptibility to cefixime to develop a molecular assay to predict cefixime susceptibility. Sex. Health 16, 488–499. doi: 10.1071/SH18227

PubMed Abstract | Crossref Full Text | Google Scholar

Dillard, J. P., Woodhams, K. L., and Ramsey, M. E. (2011). The gonococcal genetic island and type IV secretion in the pathogenic Neisseria. Front. Microbiol. 2:61. doi: 10.3389/fmicb.2011.00061

PubMed Abstract | Crossref Full Text | Google Scholar

Eisenberg, D., Schwarz, E., Komaromy, M., and Wall, R. (1984). Analysis of membrane and surface protein sequences with the hydrophobic moment plot. J. Mol. Biol. 179, 125–142. doi: 10.1016/0022-2836(84)90309-7

PubMed Abstract | Crossref Full Text | Google Scholar

EUCAST (2025). The European Committee on Antimicrobial Susceptibility Testing. Breakpoint tables for interpretation of MICs and zone diameters. Version 15.0, 2025. European Committee on Antimicrobial Susceptibility Testing. Available online at: https://www.eucast.org/fileadmin/src/media/PDFs/EUCAST_files/Breakpoint_tables/v_15.0_Breakpoint_Tables.pdf (Accessed January 1, 2025).

Google Scholar

Fenton, B. A., Tomberg, J., Sciandra, C. A., Nicholas, R. A., Davies, C., and Zhou, P. (2021). Mutations in PBP2 from ceftriaxone-resistant Neisseria gonorrhoeae alter the dynamics of the β3–β4 loop to favor a low-affinity drug-binding state. J. Biol. Chem. 297:101188. doi: 10.1016/j.jbc.2021.101188

PubMed Abstract | Crossref Full Text | Google Scholar

Fifer, H., Doumith, M., Rubinstein, L., Mitchell, L., Wallis, M., Singh, S., et al. (2024). Ceftriaxone-resistant Neisseria gonorrhoeae detected in England, 2015–24: an observational analysis. J. Antimicrob. Chemother. 79, 3332–3339. doi: 10.1093/jac/dkae369

PubMed Abstract | Crossref Full Text | Google Scholar

Gianecini, R. A., Golparian, D., Zittermann, S., Litvik, A., Gonzalez, S., Oviedo, C., et al. (2019). Genome-based epidemiology and antimicrobial resistance determinants of Neisseria gonorrhoeae isolates with decreased susceptibility and resistance to extended-spectrum cephalosporins in Argentina in 2011–16. J. Antimicrob. Chemother. 74, 1551–1559. doi: 10.1093/jac/dkz054

PubMed Abstract | Crossref Full Text | Google Scholar

Hamilton, H. L., Domínguez, N. M., Schwartz, K. J., Hackett, K. T., and Dillard, J. P. (2005). Neisseria gonorrhoeae secretes chromosomal DNA via a novel type IV secretion system. Mol. Microbiol. 55, 1704–1721. doi: 10.1111/j.1365-2958.2005.04521.x

PubMed Abstract | Crossref Full Text | Google Scholar

Harrison, O. B., Clemence, M., Dillard, J. P., Tang, C. M., Trees, D., Grad, Y. H., et al. (2016). Genomic analyses of Neisseria gonorrhoeae reveal an association of the gonococcal genetic island with antimicrobial resistance. J. Infect. 73, 578–587. doi: 10.1016/j.jinf.2016.08.010

PubMed Abstract | Crossref Full Text | Google Scholar

Jolley, K. A., Bray, J. E., and Maiden, M. C. J. (2018). Open-access bacterial population genomics: BIGSdb software, the PubMLST.org website and their applications. Wellcome Open Res. 3:124. doi: 10.12688/wellcomeopenres.14826.1

PubMed Abstract | Crossref Full Text | Google Scholar

Laumen, J. G. E., Hieu, V. N., Nhung, P. H., Vandelannoote, K., Nguyen, T. T., Nguyen, T. T., et al. (2025). High prevalence of ceftriaxone-resistant Neisseria gonorrhoeae in Hanoi, Vietnam, 2023–2024. J. Infect. Dis. 232, e73–e77. doi: 10.1093/infdis/jiaf071

PubMed Abstract | Crossref Full Text | Google Scholar

Letunic, I., and Bork, P. (2024). Interactive Tree of Life (iTOL) v6: recent updates to the phylogenetic tree display and annotation tool. Nucleic Acids Res. 52, W78–W82. doi: 10.1093/nar/gkae268

PubMed Abstract | Crossref Full Text | Google Scholar

Li, Y., Li, Y., Xiu, L., Zeng, Y., Zhang, C., Sun, L., et al. (2021). Typing of Neisseria gonorrhoeae isolates in Shenzhen, China, from 2014 to 2018 reveals the shift of genotypes significantly associated with antimicrobial resistance. Antimicrob. Agents Chemother. 65:e02311–02320. doi: 10.1128/AAC.02311-20

Crossref Full Text | Google Scholar

Liao, Y., Xie, Q., Li, X., Yin, X., Wu, X., Liu, M., et al. (2023). Dissemination of Neisseria gonorrhoeae with decreased susceptibility to extended-spectrum cephalosporins in Southern China, 2021: a genome-wide surveillance from 20 cities. Ann. Clin. Microbiol. Antimicrob. 22:39. doi: 10.1186/s12941-023-00587-x

PubMed Abstract | Crossref Full Text | Google Scholar

Liao, Y., Xie, Q., Yin, X., Li, X., Xie, J., Wu, X., et al. (2024). penA profile of Neisseria gonorrhoeae in Guangdong, China: novel penA alleles are related to decreased susceptibility to ceftriaxone or cefixime. Int. J. Antimicrob. Agents 63:107101. doi: 10.1016/j.ijantimicag.2024.107101

PubMed Abstract | Crossref Full Text | Google Scholar

Lin, E. Y., Adamson, P. C., Deng, X., and Klausner, J. D. (2021). Establishing novel molecular algorithms to predict decreased susceptibility to ceftriaxone in Neisseria gonorrhoeae strains. J. Infect. Dis. 223, 1232–1240. doi: 10.1093/infdis/jiaa495

PubMed Abstract | Crossref Full Text | Google Scholar

Lin, X., Chen, W., Xie, Q., Yu, Y., Liao, Y., Feng, Z., et al. (2022). Dissemination and genome analysis of high-level ceftriaxone-resistant penA 60.001 Neisseria gonorrhoeae strains from the Guangdong Gonococcal antibiotics susceptibility Programme (GD-GASP), 2016–2019. Emerg. Microbes Infect. 11, 344–350. doi: 10.1080/22221751.2021.2011618

PubMed Abstract | Crossref Full Text | Google Scholar

Lindberg, R., Fredlund, H., Nicholas, R., and Unemo, M. (2007). Neisseria gonorrhoeae isolates with reduced susceptibility to cefixime and ceftriaxone: association with genetic polymorphisms in penA, mtrR, porB1b, and ponA. Antimicrob. Agents Chemother. 51, 2117–2122. doi: 10.1128/AAC.01604-06

PubMed Abstract | Crossref Full Text | Google Scholar

Melendez, J. H., Edwards, V. L., Tirado, A. M., Hardick, J., Mehta, A., Aluvathingal, J., et al. (2024). Local emergence and global evolution of Neisseria gonorrhoeae with high-level resistance to azithromycin. Antimicrob. Agents Chemother. 68, e00927–e00924. doi: 10.1128/aac.00927-24

PubMed Abstract | Crossref Full Text | Google Scholar

Mitra, R., Li, J., Sagendorf, J. M., Jiang, Y., Cohen, A. S., Chiu, T.-P., et al. (2024). Geometric deep learning of protein–DNA binding specificity. Nat. Methods 21, 1674–1683. doi: 10.1038/s41592-024-02372-w

PubMed Abstract | Crossref Full Text | Google Scholar

Mlynarczyk-Bonikowska, B., Kowalewski, C., Krolak-Ulinska, A., and Marusza, W. (2022). Molecular mechanisms of drug resistance and epidemiology of multidrug-resistant variants of Neisseria gonorrhoeae. Int. J. Mol. Sci. 23:10499. doi: 10.3390/ijms231810499

PubMed Abstract | Crossref Full Text | Google Scholar

Ouk, V., Pham, C. D., Wi, T., van Hal, S. J., and Lahra, M. M. (2023). The Enhanced Gonococcal Surveillance Programme, Cambodia. Lancet Infect. Dis. 23, e332–e333. doi: 10.1016/S1473-3099(23)00479-6

PubMed Abstract | Crossref Full Text | Google Scholar

Page, A. J., Cummins, C. A., Hunt, M., Wong, V. K., Reuter, S., Holden, M. T. G., et al. (2015). Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics 31, 3691–3693. doi: 10.1093/bioinformatics/btv421

PubMed Abstract | Crossref Full Text | Google Scholar

Palace, S. G., Wang, Y., Rubin, D. H. F., Welsh, M. A., Mortimer, T. D., Cole, K., et al. (2020). RNA polymerase mutations cause cephalosporin resistance in clinical Neisseria gonorrhoeae isolates. eLife 9:e51407. doi: 10.7554/eLife.51407

PubMed Abstract | Crossref Full Text | Google Scholar

Parks, D. H., Chuvochina, M., Rinke, C., Mussig, A. J., Chaumeil, P.-A., and Hugenholtz, P. (2022). GTDB: an ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy. Nucleic Acids Res. 50, D785–D794. doi: 10.1093/nar/gkab776

PubMed Abstract | Crossref Full Text | Google Scholar

Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P., and Tyson, G. W. (2015). CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055. doi: 10.1101/gr.186072.114

PubMed Abstract | Crossref Full Text | Google Scholar

Price, M. N., Dehal, P. S., and Arkin, A. P. (2009). FastTree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol. Biol. Evol. 26, 1641–1650. doi: 10.1093/molbev/msp077

PubMed Abstract | Crossref Full Text | Google Scholar

Qin, X.-L., Chen, Y., Wu, X.-Z., Chen, W.-T., Xue, Y.-H., Huang, J.-M., et al. (2024). Emerging epidemic of the Africa-type plasmid in penicillinase-producing Neisseria gonorrhoeae in Guangdong, China, 2013–2022. Emerg. Microbes Infect. 14. doi: 10.1080/22221751.2024.2440489

PubMed Abstract | Crossref Full Text | Google Scholar

Schroeder, J. W., Hirst, W. G., Szewczyk, G. A., and Simmons, L. A. (2016). The effect of local sequence context on mutational bias of genes encoded on the leading and lagging strands. Curr. Biol. 26, 692–697. doi: 10.1016/j.cub.2016.01.016

PubMed Abstract | Crossref Full Text | Google Scholar

Seemann, T. (2014). Prokka: rapid prokaryotic genome annotation. Bioinformatics 30, 2068–2069. doi: 10.1093/bioinformatics/btu153

PubMed Abstract | Crossref Full Text | Google Scholar

Silva, M., Machado, M. P., Silva, D. N., Rossi, M., Moran-Gilad, J., Santos, S., et al. (2018). chewBBACA: a complete suite for gene-by-gene schema creation and strain identification. Microb. Genom. 4:e000166. doi: 10.1099/mgen.0.000166

PubMed Abstract | Crossref Full Text | Google Scholar

Smolarczyk, K., Mlynarczyk-Bonikowska, B., Rudnicka, E., Szukiewicz, D., Meczekalski, B., Smolarczyk, R., et al. (2021). The impact of selected bacterial sexually transmitted diseases on pregnancy and female fertility. Int. J. Mol. Sci. 22:2170. doi: 10.3390/ijms22042170

PubMed Abstract | Crossref Full Text | Google Scholar

Tang, Y., Liu, X., Chen, W., Luo, X., Zhuang, P., Li, R., et al. (2023). Antimicrobial resistance profiling and genome analysis of the penA-60.001Neisseria gonorrhoeae clinical isolates in China in 2021. J. Infect. Dis. 228, 792–799. doi: 10.1093/infdis/jiad258

Crossref Full Text | Google Scholar

Thomas, J. C., Joseph, S. J., Cartee, J. C., Pham, C. D., Schmerer, M. W., Schlanger, K., et al. (2021). Phylogenomic analysis reveals persistence of gonococcal strains with reduced-susceptibility to extended-spectrum cephalosporins and mosaic penA-34. Nat. Commun. 12:3801. doi: 10.1038/s41467-021-24072-1

PubMed Abstract | Crossref Full Text | Google Scholar

Tomberg, J., Unemo, M., Davies, C., and Nicholas, R. A. (2010). Molecular and structural analysis of mosaic variants of penicillin-binding protein 2 conferring decreased susceptibility to expanded-spectrum cephalosporins in Neisseria gonorrhoeae: role of epistatic mutations. Biochemistry 49, 8062–8070. doi: 10.1021/bi101167x

PubMed Abstract | Crossref Full Text | Google Scholar

Trinh, T. M., Nguyen, T. T., Le, T. V., Nguyen, T. T., Ninh, D. T., Duong, B. H., et al. (2022). Neisseria gonorrhoeae FC428 subclone, Vietnam, 2019–2020. Emerg. Infect. Dis. 28, 432–435. doi: 10.3201/eid2802.211788

PubMed Abstract | Crossref Full Text | Google Scholar

Unemo, M., Golparian, D., Nicholas, R., Ohnishi, M., Gallay, A., and Sednaoui, P. (2012). High-level cefixime- and ceftriaxone-resistant Neisseria gonorrhoeae in France: novelpenAMosaic allele in a successful international clone causes treatment failure. Antimicrob. Agents Chemother. 56, 1273–1280. doi: 10.1128/AAC.05760-11

Crossref Full Text | Google Scholar

Unitt, A., Maiden, M., and Harrison, O. (2024). Characterizing the diversity and commensal origins of penA mosaicism in the genus Neisseria. Microb. Genom. 10:001209. doi: 10.1099/mgen.0.001209

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, D., Li, Y., Zhang, C., Zeng, Y., Peng, J., and Wang, F. (2023). Genomic epidemiology of Neisseria gonorrhoeae in Shenzhen, China, during 2019-2020: increased spread of ceftriaxone-resistant isolates brings insights for strengthening public health responses. Microbiol. Spectr. 11:e01728–e01723. doi: 10.1128/spectrum.01728-23

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, L., Li, Y., Xiu, L., Hu, L., Huang, J., Yong, G., et al. (2024). Multigeographic clinical assessment of a molecular diagnostic assay for detection of key codons to predict decreased susceptibility or resistance to cephalosporins in Neisseria gonorrhoeae. Antimicrob. Agents Chemother. 68:e0116524. doi: 10.1128/aac.01165-24

PubMed Abstract | Crossref Full Text | Google Scholar

Whiley, D. M., Goire, N., Lambert, S. B., Ray, S., Limnios, E. A., Nissen, M. D., et al. (2010). Reduced susceptibility to ceftriaxone in Neisseria gonorrhoeae is associated with mutations G542S, P551S and P551L in the gonococcal penicillin-binding protein 2. J. Antimicrob. Chemother. 65, 1615–1618. doi: 10.1093/jac/dkq187

PubMed Abstract | Crossref Full Text | Google Scholar

WHO (2024). Multi-Drug Resistant Gonorrhoea. World Health Organization. Available online at: https://www.who.int/news-room/fact-sheets/detail/multi-drug-resistant-gonorrhoea (Accessed July 4, 2024).

Google Scholar

Workowski, K. A., and Berman, S. (2010). Sexually transmitted diseases treatment guidelines, 2010. MMWR Recomm. Rep. 59, 1–110.

PubMed Abstract | Google Scholar

Xiu, L., Wang, L., Li, Y., Hu, L., Huang, J., Yong, G., et al. (2025). Molecular screening to track ceftriaxone-resistant FC428-like Neisseria gonorrhoeae strains' dissemination in four provinces of China, 2019 to 2021. Eurosurveillance 30:2400166. doi: 10.2807/1560-7917.ES.2025.30.6.2400166

PubMed Abstract | Crossref Full Text | Google Scholar

Xiu, L., Zhang, L., and Peng, J. (2024). Surge in ceftriaxone-resistant Neisseria gonorrhoeae FC428-like strains, Asia-Pacific Region, 2015–2022. Emerg. Infect. Dis. 30, 1683–1686. doi: 10.3201/eid3008.240139

PubMed Abstract | Crossref Full Text | Google Scholar

Yu, J., Zhou, Y., Tanaka, I., and Yao, M. (2010). Roll: a new algorithm for the detection of protein pockets and cavities with a rolling probe sphere. Bioinformatics 26, 46–52. doi: 10.1093/bioinformatics/btp599

PubMed Abstract | Crossref Full Text | Google Scholar

Zhou, Z., Alikhan, N.-F., Sergeant, M. J., Luhmann, N., Vaz, C., Francisco, A. P., et al. (2018). GrapeTree: visualization of core genomic relationships among 100,000 bacterial pathogens. Genome Res. 28, 1395–1404. doi: 10.1101/gr.232397.117

PubMed Abstract | Crossref Full Text | Google Scholar

Zhu, L., Liang, J., Zheng, Y., Chen, S., Xu, Q., Yin, S., et al. (2023). Combined mutations of the penA with ftsX genes contribute to ceftriaxone resistance in Neisseria gonorrhoeae and peptide nucleic acids targeting these genes reverse ceftriaxone resistance. J. Glob. Antimicrob. Resist. 35, 19–25. doi: 10.1016/j.jgar.2023.08.005

PubMed Abstract | Crossref Full Text | Google Scholar

Keywords: decreased susceptibility to ceftriaxone, multilocus sequence typing (MLST), mutation frequencies, mutations, Neisseria gonorrhoeae, penA (penicillin-binding protein 2)

Citation: Zhang X, Wang H, Gu Y, Chen X, Zhou G, Liu C, Ji L, Xie R, Zhang J, Shao Z and Zhang M (2026) Genetic identification of mutations and MLST types associated with decreased susceptibility to ceftriaxone in Neisseria gonorrhoeae. Front. Microbiol. 16:1728860. doi: 10.3389/fmicb.2025.1728860

Received: 21 October 2025; Revised: 25 December 2025;
Accepted: 29 December 2025; Published: 21 January 2026.

Edited by:

Eliana De Gregorio, University of Naples Federico II, Italy

Reviewed by:

Dmitry Gryadunov, Engelhardt Institute of Molecular Biology (RAS), Russia
Mona Mohamed Sobhy, Animal Production Research Institute (APRI), Egypt

Copyright © 2026 Zhang, Wang, Gu, Chen, Zhou, Liu, Ji, Xie, Zhang, Shao and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Maojun Zhang, emhhbmdtYW9qdW5AaWNkYy5jbg==; Zhujun Shao, c2hhb3podWp1bkBpY2RjLmNu

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.