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.

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).

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.

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.

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.

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.

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.

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.

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 3A–C). 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.

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.

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.

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).

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.

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 (χ² = 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.