Chlamydia trachomatis is an obligate intracellular bacterium and causes the most common bacterial sexually transmitted infection (STI) worldwide, with 129 million new cases estimated annually, according to the World Health Organization (WHO) ¹. The classification of C. trachomatis strains has historically relied on ompA-genotyping, a method based on sequencing the gene encoding the major outer membrane protein (MOMP), which is highly variable and immunodominant ². Fifteen major ompA-genotypes are defined: A, B/Ba, C, D/Da, E, F, G, H, I/Ia, J/Ja, K, L1, L2, L2b and L3 ³, which correlate with three main clinical disease groups. Specifically, genotypes A-C are associated with ocular infections and trachoma, constituting the leading cause of preventable blindness in the world ^4,5. Genotypes D-K are associated with urogenital infections, the most frequent manifestation of C. trachomatis, which are often asymptomatic, particularly in women, contributing to underdiagnosis and increased risk of pelvic inflammatory disease (PID), ectopic pregnancy, and/or infertility ^5,6. Finally, ompA-genotypes L1-L3 are associated with lymphogranuloma venereum (LGV), characterised by its invasive nature and severe clinical manifestations such as inguinal lymphadenopathy and proctitis, often observed in high-risk populations such as men who have sex with men (MSM) ^5,7–9.

In 2003, an outbreak of L2b C. trachomatis strains, first detected in Europe, marked the re-emergence of LGV as a significant public health concern, prompting renewed attention to its surveillance and diagnosis strategies ^13–15. Consequently, rapid PCR assays targeting potential LGV-specific molecular markers, such as the 36 bp deletion in the pmpH gene, were introduced and widely adopted for diagnostic and epidemiological purposes ^16–19. While these assays improved routine detection of LGV in clinical settings, their inability to capture circulating diversity, for example, to identify novel, hybrid, or recombinant strains, has raised concerns about potential misclassifications and limited epidemiological insight ^10–12. Other complementary methods, such as multi-locus sequence typing (MLST), have been proposed to improve strain discrimination beyond rapid LGV detection and ompA-genotyping ^9,20,21. However, their uptake in routine diagnostic settings remains reduced ²². In fact, most laboratories do not even perform ompA-genotyping, and practical constraints associated with MLST have limited its widespread adoption. Indeed, MLST requires samples with a relatively high bacterial load and provides only modest resolution, limiting its effectiveness for epidemiological surveillance, outbreak detection, and transmission tracking ^20,23.

In this context, whole-genome sequencing (WGS) has emerged as the decisive tool for studying C. trachomatis evolution and transmission. With its superior resolution, WGS can identify and track new clones, reveal recombination events, and uncover microevolutionary changes with significant clinical and epidemiological relevance, all at a finer scale than traditional typing methods. Foundational WGS studies ^7,24,25 revealed that C. trachomatis genome is characterized by a high protein-coding density, few pseudogenes and few non-essential genes (in contrast to other obligate intracellular pathogens, such as Rickettsia prowazekii) ^26,27. Also, it presents a nearly overlapping core-genome [i.e., genes shared by (virtually) all strains of a given species] and pan-genome [i.e., all genes described in a given species], with a single non-chromosomal element, a ~ 7.5 Kbp plasmid that co-evolved with the chromosome ^5,28,29 and is known to play a key functional role in regulating chromosomal genes ^30–38. Subsequent landmark WGS studies by Harris et al. (2012) ⁵ and Hadfield et al. (2017) ²⁹ unveiled the presence of four deep-branching, evolutionarily distinct lineages, which broadly correlate with the three historical C. trachomatis disease groups (LGV, ocular and urogenital), namely: the LGV lineage, the “T1 urogenital” lineage (essentially composed of prevalent genotypes E and F), the “T2 urogenital” lineage (mainly enrolling non-E/F genital genotypes) and the “T2 ocular” lineage. Following this ancient diversification, C. trachomatis genomes have been extensively reshaped by frequent homologous recombination and marked by recent clonal expansions ^5,29,39. Comparative genomic analyses have further demonstrated that recombination occurs both intragenically and across dispersed genomic regions, involving strains with similar or distinct tissue tropisms, and with no absolute barriers to genetic exchange ^5,29,39–42. Opportunities for recombination appear higher among strains infecting the same tissue niche, with higher gene flow observed within urogenital lineages compared to LGV or ocular lineages ⁴⁰. Not unexpectedly, WGS has enabled the identification of epidemiologically and clinically relevant clones that challenged conventional molecular typing strategies, such as the L2b epidemic strain associated with a proctitis outbreak among MSM, including the recombinant L2b/D-Da strains ^{29,39,43–45}. More recently, genomic analyses have described the emergence of a novel LGV strain, designated L4 due to its ompA signature distinct from the classical L1-L3 genotypes (45). First identified in South America ^46,47, L4 has shown signs of potential increasing international circulation ⁴⁸. These findings underscore the added value of implementing high-resolution WGS-based surveillance to adequately and timely capture the evolving landscape of C. trachomatis, including the emergence of sublineages with potential clinical and epidemiological relevance. Still, routine implementation of WGS in C. trachomatis surveillance remains very limited, primarily due to technical and operational constraints associated with its obligate intracellular nature. In order to overcome these limitations, advances including targeted enrichment approaches such as RNA- or DNA-based hybrid capture have enabled successful WGS directly from clinical samples with sufficient bacterial load allowing for culture-independent, high-resolution genomic surveillance ^29,49–51. However, these protocols are often costly, labor-intensive, and not widely standardized across laboratories. Moreover, the lack of harmonized analytical pipelines and curated genomic typing schemas for C. trachomatis, like those applied for other human bacterial pathogens, such as Listeria monocytogenes or Salmonella enterica, limits the practical integration of WGS into public health systems ^52,53. For instance, core-genome MLST (cgMLST) schemas are well-established for outbreak investigation and surveillance of many bacterial pathogens, but very little explored for C. trachomatis. Although a cgMLST schema has been developed for C. trachomatis in 2018 ⁵⁴, being publicly available in the PubMLST platform ⁵⁵, to our knowledge, its application has been limited to the original study. Its design for online usage, together with the absence of benchmark studies with state-of-the-art allele callers for local deployment, like the widely used chewBBACA ⁵⁶, may justify its lower adoption.

Recognizing the value of cgMLST for advancing C. trachomatis molecular surveillance, this study introduces a novel cgMLST schema tailored for local deployment, providing a standardized and scalable framework for comparing genome data across laboratories and studies. By analysing over 1200 genome sequences, this work provides new insights into the global phylogenetic structure and lineage diversity of C. trachomatis, representing one of the most extensive genomics efforts to date. Beyond public health applications, the generated cgMLST resources, like the in-depth allelic mapping of current C. trachomatis global diversity, might be of great potential to accelerate research, thereby bridging surveillance and scientific investigation.

C. trachomatis samples collection and selection for WGS

The National Reference Laboratory (NRL) for STIs, at the National Institute of Health Doutor Ricardo Jorge (INSA), Lisbon, Portugal, maintains a vast collection of C. trachomatis-positive clinical specimens collected through routine diagnostics or voluntarily sent by collaborating laboratories across the country since 1990. By the end of 2021, a total of 5824 samples had been successfully ompA-genotyped, as described in Lodhia et al., 2025 ⁵⁷. Among these, 1188 samples were further subtyped for LGV-specific ompA subvariants ⁴⁸, and 502 samples were subjected to multiplex sequencing to survey genetic markers of antimicrobial resistance (AMR) and infer strain-specific genomic backbones ⁵⁸. To select samples for WGS, a dataset of 850 DNA C. trachomatis-positive samples was assembled (Supplementary file S1). These samples had been processed using two distinct approaches: 15 DNA samples derived from clinical specimens that were cultured in McCoy cells (standard conditions, as previously described ⁵⁹, and 835 DNA extracts obtained directly from C. trachomatis-positive clinical specimens that were further sub-selected for targeted enrichment WGS protocol. The genomic DNA of the 15 cultured samples was obtained using a scale‑down version of the urografin‑based purification protocol and subjected to standard WGS, as described by Pereira et al. ⁶⁰. For the 835 C.trachomatis-positive clinical DNA samples, bacterial load was quantified using an absolute quantification real-time PCR (qPCR) assay targeting the number of C. trachomatis copies, as previously described by Gomes et al., 2006 ⁶¹. This method enabled absolute quantification of bacterial genome copies per µL of DNA. Based on our previous experiences with SureSelect targeted enrichment approach ^29,50, indicating a positive correlation between higher bacterial load and successful genome recovery, samples with ≥ 10⁴ genome copies in the input volume (7 µL) were prioritized for sequencing. A total of 246 samples were evaluated for dsDNA quantification (using Qubit, ThermoFisher Scientific). A final subset of 125 samples, passing both the bacterial load (≥ 10⁴ genome copies in 7 µL) and DNA quantity (> 10ng in 7 µL) minimum thresholds defined, was selected for targeted WGS. Additionally, 23 samples not meeting the ≥ 10⁴ genome copies threshold were also included, given their epidemiological/clinical interest and acceptable DNA quantity to proceed with the protocol. Full details on bacterial load and DNA quantification are provided in Supplementary file S1.

Targeted C. trachomatis whole-genome capture and sequencing directly from clinical samples

Whole-genome sequencing was performed using a targeted enrichment approach based on the SureSelect XT HS technology (G9702-90000, version F0, September 2022, Agilent Technologies, Santa Clara, CA, United States), optimized for C. trachomatis, employing custom-designed RNA oligonucleotide baits targeting the C. trachomatis genome, as previously described by Borges et al., 2019 ⁴³, using a 1:5 bait dilution. Library preparation began with the fragmentation of high-quality DNA using the Agilent’s SureSelectXT HS Low input Enzymatic Fragmentation kit and was further carried out following the manufacturer’s instructions. Libraries fragments size, concentration and molarity were determined on a Fragment Analyzer system (Agilent, Santa Clara, CA, United States). Libraries were then sequenced on Illumina MiSeq or NextSeq platforms (Illumina, San Diego, CA, United States), applying paired-end multiplexed sequencing. Descriptive and laboratory data, including sequencing metrics, are detailed in Supplementary file S1.

Quality control and de novo genome assembly

The analysis of the WGS data of C. trachomatis relied on an assembly-based approach targeting C. trachomatis sequencing reads. First, Quality Control (QC) and trimming of sequencing reads was performed with INNUca pipeline v4.2.2 (https://github.com/B-UMMI/INNUca), using FastQC v0.11.5 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and Trimmomatic v0.38 ⁶². Kraken v2.1.3 ⁶³ was used to classify trimmed reads using Kraken2 standard 16Gb database downloaded on September 4th, 2024. Trimmed reads classified by Kraken2 as Chlamydiales (taxon ID: 204428) were then extracted using the extract_kraken_reads.py script of Kraken Tools v1.2 ⁶⁴ and using the “--include-children” argument. Samples without sequencing reads classified as belonging to this taxon ID or respective “children” were discarded from the analysis. For the remaining, the filtered sequencing reads classified as belonging to this taxon ID or respective “children” were assembled with INNUca pipeline v4.2.2 (https://github.com/B-UMMI/INNUca), using SPAdes v3.14 ⁶⁵ and Pilon v1.23 ⁶⁶. Only assemblies with at least 90% and no more than 110% of the expected genome size (compared with the reference genome NC_000117.1: 1 042 519 bp) were used for further analysis. At a later stage, the number of loci called in the novel cgMLST schema (see subsection “Creation of a novel cgMLST schema for C. trachomatis”) was also used as a proxy for assembly quality, and all isolates with less than 95% loci called were excluded.

Genetic clustering using the novel cgMLST schema was obtained for all possible threshold levels with ReporTree v2.5.4 ⁸², providing the desired allele matrix and metadata as input, and using the single-linkage or the GrapeTree ⁸³ clustering method, depending on the analysis. Samples with less than 95% of the cgMLST loci called were excluded with the ‘--loci-called’ argument. The arguments ‘--columns_summary_report’ and ‘--metadata2report’ were used to characterize the different genetic clusters or metadata groups, respectively, with the available metadata (Supplementary file S2). ReporTree was also used to obtain count and frequency tables with the distribution of a given metadata variable of interest over another one (selected compilation in Supplementary file S4, and full compilation in Zenodo repository https://doi.org/10.5281/zenodo.17177579).

In order to establish a hierarchical nomenclature that allows monitoring the main circulating clusters, a first threshold grouping the whole dataset into four clusters, corresponding to the four expected lineages ^5,29, was determined based on the results obtained for the hierarchical clustering of the whole dataset (Supplementary file S5). Then, to define thresholds with higher discriminatory power, cluster stability regions, i.e. threshold ranges providing similar clustering, were determined with comparing_partitions_v2.py (https://github.com/insapathogenomics/ComparingPartitions) ⁸², following a previously described approach ^84–87. Stability regions were defined as at least 5 subsequent thresholds yielding an Adjusted Wallace Coefficient (AWC) ≥ 0.98 (Supplementary file S6). Two stability phases at high resolution were selected to cluster isolates into high- and low-level WGS-based genogroups. An additional ReporTree run was then performed to classify each isolate into a specific genogroup by requesting an hierarchical nomenclature code attributed based on the hierarchical clustering results at the thresholds corresponding to lineage, high- and low-level genogroup classification, using a threshold that fell between the edges of the respective stability region. The final code has the structure LIN-HG-LG, where “LIN” is a lineage code (LGV for LGV lineage, T2O for T2 ocular, T2G for T2 urogenital and T1G for T1 urogenital), and “HG” and “LG” correspond to the high- and low-level genogroup classification, respectively. The genogroup codes correspond to the cluster or singleton nomenclature attributed by ReporTree ⁸²: “Cn” when the isolate belongs to a cluster and “Sn” when the isolate is a singleton, with “n” corresponding to the respective cluster/singleton number.

The assessment of cluster congruence between two different typing approaches (e.g., cgMLST-based and SNP-based clustering analysis) was performed using a previously developed methodology ^84,85,87 with EvalTree v1.0.0 (https://github.com/insapathogenomics/ACDTree/tree/main/EvalTree), providing the ReporTree output folder of each method (obtained as described in the “Cluster detection and characterization” section) as input. Briefly, EvalTree compares the genetic clusters at each possible threshold levels in the two methods (available in ReporTree “partitions” table) and calculates a congruence score (CS) between them based on the Adjusted Rand and Adjusted Wallace coefficients ^84,85,87. Inter-pipeline corresponding points, i.e. the allele/SNP threshold in one pipeline that provided the most similar clustering results in the other, were assessed only considering comparisons rendering CS ≥ 2.85.

Comparison with C. trachomatis PubMLST cgMLST schema

To compare the resolution and clustering congruence between the novel cgMLST and the cgMLST schema available in PubMLST (https://pubmlst.org/) ^54,55, which was designed for online usage within this platform, each of the 817 loci of the latter was downloaded on August 6th, 2025. For its local deployment with chewBBACA allele caller ⁵⁶, the schema was prepared with chewBBACA v3.10 PrepExternalSchema module using default settings and providing the same Prodigal training file as used for the novel schema creation (see section “Creation of a novel cgMLST schema for C. trachomatis”). Afterwards, allele calling with chewBBACA v3.10 AlleleCall module was iteratively performed using the curated dataset (n = 1230) (Supplementary file S2) until no new alleles were detected (5 rounds). The final allele matrix was used for cluster analysis, as described in the “Cluster detection and characterization” section.

To compare the resolution and clustering congruence of the cgMLST approach with traditionally applied SNP-based methods, three comparative exercises were conducted. The first analysis involved comparison with the 53 complete C. trachomatis chromosome sequences used for the initial construction of the cgMLST schema. For this, the full chromosome sequences were aligned using the progressive alignment algorithm implemented in Mauve v2.3.1 ⁸⁸. A core alignment was then extracted by retaining only regions where all chromosome sequences aligned over at least 500 bp. The resulting core genome alignment was manually inspected and corrected, yielding 1 007 702 nucleotide positions (~ 97% of the chromosome), including 16 532 single-nucleotide variant (SNV) sites.

The second analysis included the entire dataset comprising 1230 genome assemblies that passed QC criteria (see Supplementary file S2). These assemblies were aligned using Parsnp v2.0.5 ⁸⁹, with default parameters, except for the –C parameter, which was set to 2000 to maximize reference coverage. The resulting core SNV alignment used in this evaluation contained 20 941 variant sites.

Finally, targeted analyses were conducted using a reference-based read mapping and SNP calling approach for two LGV clusters of epidemiological interest: the epidemic L2b strain ^5,7,29,69 and the emergent L4 strain ⁴⁶. For this, QC-filtered reads from samples identified by cgMLST within each cluster (L2b or L4) (Supplementary file S2) were mapped using Snippy v4.6.0 (https://github.com/tseemann/snippy; --mapqual 30 --mincov 5 --minfrac 0.9 --basequal 20) against representative genome sequences for each strain: the complete genome of L2b/UCH-1/proctitis (AM884177.2/AM886279.1) ⁷, and the assembly of an L4 isolate from 2018 (the earliest year with an L4 genome analyzed in the present study) selected based on best contig metrics (number and N50) among those with available reads (ENA accession: ERR13385367), respectively. Core genome SNV alignments were then extracted from the Snippy-derived whole alignment using the “alignment processing module” of ReporTree v2.5.4 ⁸², with “--site-inclusion” 0.9 (i.e., only keeping variant positions with < 10% missing data). For the L2b cluster, all recombination-driven SNVs within the previously identified ompA containing recombinant region (NC_010280; approximate positions 55221–59461) of the recombinant L2b/D-Da strain ^43,45 were additionally excluded. The final core SNV alignments, comprising 156 SNVs for L2b and 16 SNVs for L4, were analyzed using single-linkage hierarchical clustering in ReporTree v2.5.4 (65) and maximum-likelihood phylogenetic reconstruction with IQ-TREE v2.1.4 ⁹⁰, employing the automatically detected best-fit model (K3P + ASC) and 1000 bootstrap replicates.

Compilation of a WGS dataset with the global C. trachomatis genomic diversity

With the aim of creating a novel cgMLST schema for C. trachomatis and explore the global species diversity, a large and diverse C. trachomatis genome dataset with 1173 assemblies (publicly available or newly generated) was compiled and curated (Supplementary file S2, details in Material and Methods). Given the broad lack of investment in WGS of this prevalent human pathogen (e.g., only 83 of 1173 genomes had available collection year of 2015 onwards), this study further sought to advance national and international genomic surveillance by sequencing more clinical isolates from the Portuguese NRL collection. For this, 850 C. trachomatis-positive DNA samples were evaluated for eligibility for targeted-enrichment WGS by estimating bacterial genome copies via qPCR and quantifying dsDNA in the required input volume (Supplementary file S1). In total, 125 samples met both selection criteria (≥ 10⁴ genome copies and > 10 ng DNA in 7 µL), being considered as eligible, while a small subset (n = 23) not fulfilling the copy number threshold was still selected for targeted-enrichment WGS (Supplementary figure S1A). At the end, 45 (30.41%) of the samples sequenced with a “targeted-enrichement” approach yielded complete genome assemblies (Supplemental table S1), with samples with higher bacterial loads generally producing a greater proportion of “on-target” reads (Supplementary figure S1B). Indeed, although the genome assembly success was 34.4% (43/125) for samples fulfilling the adopted threshold of ≥ 10⁴ copies in WGS input, this value increases, for example, to 64.7% (22/34) when considering samples with ≥ 10⁵ copies in WGS input. This pattern was reflected in the proportion of “on-target” reads, which had a median percentage of 48.84% (IQR: 39.60–70.17) among samples yielding complete genomes (n = 45), compared with 5.16% (IQR: 1.43–12.43) observed in samples that failed to meet the criteria for inclusion in the final genome dataset (n = 103) (Supplementary figure S1B). Besides the 45 genomes recovered directly from DNA samples, complete genome assemblies were also obtained for 12 cultured clinical isolates (Supplementary file S1), which were included in the final curated dataset (Supplementary file S2). Overall, this effort led to the inclusion of 57 new C. trachomatis genomes from Portugal, increasing by 17% the proportion of European isolates in the dataset and by 56% the amount of anorectal samples. The final dataset comprises 1230 C. trachomatis genome assemblies, spanning isolates collected between 1957 and 2023, from different clinical origins, and representing at least 26 countries (Supplementary file S2).

Development of a novel cgMLST schema for large-scale genomic analysis and monitoring of C. trachomatis diversity and evolution

To leverage a high resolution WGS-based typing method that can be easily deployed in local settings for routine genomic monitoring of C. trachomatis diversity and evolution, in this study, a novel cgMLST schema was developed for this species. Schema construction was based on 53 complete chromosomal sequences representing the species diversity ^5,67 and 42 complete plasmid sequences (listed in Supplementary file S2). The final version was populated iteratively using the curated dataset (n = 1230) (Supplementary file S2) until no new alleles were detected (3 rounds). The final cgMLST schema is composed of 846 loci (838 CDSs from the chromosome and 8 from the plasmid), covering ~ 85–88% of expected genome size and representing a total of 26,528 unique alleles (median: 27; min: 2; max: 182 per locus, Supplementary figure S2). When compared with the local deployment of the previously available PubMLST schema (designed primarily for online usage) ⁵⁴, the novel schema showed highly congruent clustering but a slight increase in resolution power (Supplementary figure S2). For instance, clusters defined at 10 allelic differences (ADs) in the novel schema corresponded to 9 ADs (CS = 2.999) in PubMLST schema (Supplementary file S5). Most importantly, the novel schema achieved a higher proportion of loci called per sample (median 99.2%; range 95.1–99.9%) than PubMLST schema (median 96.7%; range 92.4–97.8%). These results demonstrate the enhanced resolution and typeability of the novel schema in local laboratory settings using the widely used chewBBACA allele caller ⁵⁶. In addition, the performance of the novel schema was tested by comparing two popular clustering algorithms, single-linkage hierarchical clustering and GrapeTree MSTreeV2 ⁸³. As expected ⁸⁷, they presented a high cluster congruence at all resolution levels, (Supplementary file S5), demonstrating the flexibility for implementation of this cgMLST schema in the commonly used workflows for bacterial genomic surveillance ⁹¹.

To further assess its discriminatory power and ability to capture C. trachomatis genomic population structure, the novel C. trachomatis cgMLST schema was also benchmarked with SNP-based approaches. First, it was compared with a curated core genome alignment of the set of 53 complete C. trachomatis genomes (Supplementary file S2), covering ~ 97% genome and enrolling 16 532 SNVs. This exercise demonstrated that, despite inferring lower genetic distances between the samples (as expected), the novel cgMLST approach is able to capture the expected species population structure, presenting absolute clustering congruence with the core genome alignment at several thresholds of resolution (Supplementary figure S2 and Supplementary file S5). To further support these observations, a second exercise was performed with the larger and more diverse dataset of C. trachomatis (n = 1230 assemblies) (Supplementary file S2), comparing the novel cgMLST with an assembly-based SNP approach with Parsnp ⁸⁹, which yielded a core SNP alignment of 20 941 SNVs. The results not only corroborated the previous observation but also demonstrated that the novel cgMLST provides similar, or even better, resolution than the assembly-based SNP approach (Supplementary figure S2). Indeed, both methods discriminate > 900 profiles for the 1230 isolates, showcasing the high discriminatory power of this novel schema (Supplementary file S5).

The implementation of cgMLST-based hierarchical nomenclature systems, increasingly adopted in genomic surveillance of multiple bacterial pathogens ^92,93, can add substantial value by enabling stable and comparable tracking of main circulating genogroups across studies and laboratories. Therefore, in order to leverage a surveillance-oriented nomenclature for the novel cgMLST schema, C. trachomatis clustering patterns were explored towards the identification of hierarchical thresholds that could be adequate for nomenclature design. First, the minimum threshold required to segregate the global hierarchical clustering tree (Fig. 1) into the expected C. trachomatis four main lineages (“T1 urogenital”, “T2 ocular”, “T2 urogenital” and LGV) ^5,29 was identified, corresponding to 449 ADs (Supplementary files S5). This threshold fell within a large clustering stability region (Supplementary file S6), which supported the rationale of having a first nomenclature level with direct correlation with lineage classification. This analysis further highlighted a large stability region between 63 and 86 ADs, which could be used to define genogroups for longitudinal surveillance, and another stability region between 23 and 30 ADs, which opened the possibility of extending the genogroup definition to a more discriminatory level. Based on these results, an intuitive cgMLST-based nomenclature system for C. trachomatis was designed, consisting on a hierarchical code that reflects the isolate clustering at lineage level (475 ADs), high-level genogroup (75 ADs) and low-level genogroup (25 ADs) (see Material and Methods), thus capturing the global population structure at different resolution levels. In total, 473 different nomenclature codes were attributed to the 1230 C. trachomatis isolates, covering 179 and 473 high- and low-level genogroups, respectively (Figs. 1 and 2, and Supplementary file S2). Of note, taking advantage of the clustering congruence analysis described above, the GrapeTree thresholds providing the most congruent results with these hierarchical code levels selected for nomenclature were identified and confirmed to also fall within stability regions, specifically 525 ADs for lineage classification, 86 ADs for high-level genogroup and 27 ADs for low-level genogroup (Supplementary file S5). These data provide the basis to directly implement this novel nomenclature system into workflows relying on any of these two clustering methods (Supplementary file S5).

The novel cgMLST schema and associated resources for immediate applicability into standardized, high-resolution WGS-based routine genomic surveillance are publicly available at Zenodo repository (https://doi.org/10.5281/zenodo.17177579). An adapted version of the schema is also available at chewie-NS (https://chewbbaca.online/) ⁹⁴.

Capturing the global population structure of C. trachomatis and major genotype-lineage incongruences by cgMLST

The cgMLST clustering results of the studied 1230 C. trachomatis genomes were used to explore the genomic diversity of this species. As mentioned before, this analysis confirmed the expected clustering into the four previously defined lineages (Figs. 1 and 2). As expected, the lineages strongly - but not exclusively - correlate with the three historical “disease groups” (ocular, genital and LGV) inferred from traditional ompA classification (Figs. 1–3 and Supplementary file S2). For instance, T1 (n = 376) and T2 urogenital (n = 281) lineages are composed mostly of ompA genital strains (97.3% and 93.2%, respectively), although both also included ompA ocular strains (2.7% and 6.8%, respectively). While T1 lineage predominantly includes the prevalent genotypes E (61.2%) and F (20.7%), the T2 urogenital lineage shows a broader ompA-genotype composition, comprising genotypes G (31.7%), D (20.3%), K (19.2%), J (7.5%), B (6.8%), Ia (6.4%), H (6.0%), F (1.1%), I (0.7%), and Da (0.4%). By other side, the T2 ocular lineage (n = 371) consisted exclusively of ompA ocular strains (100.0%), distributed as follows: A (67.7%), Ba (31.3%), B (0.5%) and C (0.5%). Finally, the LGV lineage (n = 202) included predominantly ompA LGV strains (94.6%), namely L1 (8.9%), L1/L2 hybrid (3.0%), L2 (41.6%), L2b (34.2%), and L3 (1.5%), and the emerging L4 (5.4%), but also atypical strains carrying non-LGV ompA-genotypes (5.4%) (Fig. 3). These included 10 isolates of the recently described recombinant L2b/D-Da strain ^43,45, as well as a genotype E strain (accession ERR351523) displaying a highly atypical recombinant profile (detailed below).

In total, 30 samples (besides the recombinant L2-L2b/D-Da strain characterized elsewhere ^43,45) displayed major incongruence between their ompA-genotype and phylogenomic lineage (Figs. 1–3 and Supplementary file S2). As described above, the most striking profile corresponded to a sequence (accession ERR351523; genogroup LGV-S3-S5) carrying an ompA-genotype E and yet clustering within the LGV lineage (Figs. 1–3). Comparative allelic profiling suggested that this strain likely resulted from the transfer of at least 229 cgMLST loci from a T1 strain into an LGV genomic background, enrolling the exchange of at least seven distinct genomic regions (each spanning a minimum of 3 genes) via homologous recombination (Fig. 2, Supplementary files S7 and S8). Noteworthy, within the LGV lineage, a distinct group of 22 strains (composing the high-level genogroups LGV-C4, LGV-C7 and LGV-C10) that, while retaining an LGV “L2” ompA-genotype, exhibit substantial mosaicism LGV-T1 was also identified (Fig. 2, Supplementary files S7 and S8). Several of the exchanged regions, particularly the large regions between loci CT_374-CT_402, CT_440-CT446 and CT_542-CT650, partially overlap those identified in ERR351523, suggesting they might represent potential recombination hotspots between T1 and LGV strains (Fig. 2).

Other notable inter-lineage ompA exchanges include 19 strains carrying ocular ompA-genotype B, but clustering within the T2 urogenital lineage. The cgMLST analysis shows that they are not confined to a single monophyletic (clonal) branch but are instead dispersed across different high-level genogroups, consistent with multiple and independent ocular-urogenital ompA exchange events. These include three genogroups composed exclusively of recombinant B isolates (T2G-C60, T2G-C84 and T2G-S77), as well as two genogroups that also harbor isolates with “genital” ompA (T2G-C57, containing ompA-genotypes B, D/Da, G and H, and T2G-C61, containing ompA-genotypes B and D) (Figs. 1 and 3B). Noteworthy, only the high-level genogroup T2G-C84 includes isolates from the past decade, among them two recombinant B anorectal isolates (from 2019 and 2023) from Portugal, possibly representing a recent sublineage that is currently circulating. Comparative allelic profiling (Fig. 2, Supplementary files S7 and S8) suggests that these recombinant T2G-C84 isolates possess a genomic backbone largely consistent with the T2 lineage, with only a limited set of loci (n = 126) carrying strong evidence of inter-lineage recombination (i.e., having alleles exclusively shared between ocular and urogenital T2 lineages).

Regarding the detection of ompA-genotypes Ba (n = 2) and C (n = 8) within the T1 lineage (Fig. 3A), these isolates belong to two divergent and highly clonal genogroups (T1G-C25-C55 and T1G-C17-C20, respectively) exclusively composed of ocular ompA-genotype isolates from the 1980s and 1990s. Moreover, among other infrequent genotype-lineage combinations (Fig. 3A), an anorectal G strain (SRR25447256) within T1 (genogroup T1G-C16-C17) is noted, since this genotype is typically associated with the T2 urogenital lineage (Figs. 1 and 3). Conversely, three ompA-genotype F strains (genogroup T2G-C78-C134) appeared as tightly clustered within the T2 urogenital lineage, although this genotype is normally predominant in T1 (Fig. 3) ^5,29,80.

Finally, a closer examination of the ompA allelic profile distribution by lineage (Supplementary figure S3 and Supplementary file S4) revealed additional notable trends. For instance, although “D-like” genotypes are generally reported in both T1 and T2 urogenital lineages (23), with the exception of one sequence in each group, D ompA alleles occurred exclusively in T2, while Da alleles were restricted to T1. Similarly, all ompA-genotype J strains (ompA alleles 14 and 26) were found in the T2 urogenital lineage, whereas genotype Ja strains (ompA alleles 56 and 110) appeared exclusively in T1. Noteworthy, cases where all strains of a given ompA-genotype belonged to the same lineage and formed a well-defined, genetically homogeneous cluster at cgMLST level were rare (Fig. 1). These cases were mainly observed in genotypes with very low representation in the dataset, such as genotype I (n = 2, genogroup T2G-C83-C149) and L3 (n = 3, genogroup LGV-C9-C9), or in cases of clonal expansion of recently emerged, epidemiologically relevant LGV strains, namely the the epidemic “L2b” strain (genogroup LGV-C1-C1) ^5,7,29,69, which includes the recombinant strain with hybrid L2-L2b/D-Da ompA (alleles 63 and 77) (n = 10) ^43,45, and the emergent L4 strain (genogroup LGV-C5-C4; ompA alleles 49 and 89) (n = 11) ⁴⁶. A more in-depth analysis within the LGV lineage is presented in the next section.

Exploring C. trachomatis epidemiologically-relevant genogroups and novel insights for genomic surveillance

The cgMLST analysis of 202 C. trachomatis genomes belonging to the LGV lineage (Supplementary file S2) revealed considerable intra-lineage diversity, allowing clear and straightforward identification of both previously described and emerging subvariants (Fig. 1) ^{5,7,29,70,95–99}. Of greatest current epidemiological and clinical relevance, two well-defined LGV clusters stand out: the epidemic “L2b” strain ^5,7,29,69 and the emergent L4 strain ⁴⁶.

Regarding the epidemic “L2b” strain, the cgMLST analysis effectively captured its monophyly ^5,29, clustering all isolates (n = 112) at a single-linkage threshold of 23 ADs, thus placing them all within the same low-level genogroup (LGV-C1-C1) (Fig. 1). Within this genogroup, there was a clear separation of a subcluster containing the earliest L2b isolates (n = 2) identified in the USA in the 1980s ^7,97–99, and another comprising all genomes (n = 10) of the recombinant L2b/D-Da strain, which grouped within a 10 AD threshold, consistently with its ongoing global expansion ^43,45. Notably, a typical core SNP maximum likelihood analysis (mapping against the L2b/UCH-1 reference genome) also segregated the same subclusters of interest (Supplementary figure S4), showing that both genome-scale typing approaches present comparable discriminatory power in capturing microevolution-driven divergence during L2b clonal expansion. In the studied dataset, within the “L2b” epidemic clone, three main ompA-genotypes were observed, namely L2b (n = 69), L2 (n = 33), and the hybrid L2-L2b/D-Da (n = 10), comprising seven, four and two distinct ompA alleles, respectively. Noteworthy, although most isolates genotyped as “L2” clustered within a separate sub-branch, L2 sequences were also found in other sub-branches. This pattern, which was observed in both SNP (Supplementary figure S4) and cgMLST (Fig. 1) analyses, supports the previously proposed hypothesis that the reported reversion of the ompA “L2b” signature” (mutation A485G) ^100,101 had occurred independently on multiple occasions throughout the epidemic. When examining the L2b phylogeny at an in-depth sub-cluster level, particularly the placement of ompA-genotype LGV subvariants (i.e., harboring distinct ompA alleles), which were recently compiled ⁴⁸, all “L2bV5” isolates (n = 17) were found to be clustered together in the core SNP phylogeny (Supplementary figure S4), a pattern less evident in the cgMLST clustering analysis (Fig. 1). L2bV5 is the second most frequent L2b subvariant in the genome dataset (after the prototype “L2b”) and was also frequent in a large LGV ompA-genotyping study in Portugal (18.5%, 2007–2023) ⁴⁸, supporting its epidemiological relevance. By suggesting a L2bV5 monophyletic nature, these results indicate that its characteristic ompA allele (allele 48 in cgMLST, with G271A and C493A mutations leading to Ala91Thr and His165Asn protein changes ⁴⁸) can still be a reliable marker for tracking its international spread. Of note, a GyrA mutation potentially associated with fluoroquinolone resistance (Ser83Ile; allele 36 of the locus "CT_189_gyrA") was previously identified in LGV clinical samples, including recombinant L2b/D-Da (42, 44). In the present genome dataset it was only observed in one genome of that recombinant strain, corroborating its potential rarity in vivo.

Regarding the emergent L4 strain, cgMLST robustly clustered all L4 genomes (n = 11) within a deeply separated and highly homogeneous genetic cluster within LGV lineage, consistent with its recent description (66). Indeed, the L4 genogroup (LGV-C5-C4) could be segregated from other LGV genomes at a cgMLST threshold as high as 229 ADs. As with L2b, the intra-subcluster diversity detected by cgMLST was very low, and of similar magnitude to that obtained with core SNP analysis, with all L4 genomes presenting a pairwise median distance of 8 ADs and 8 SNPs.

In addition to the two LGV clusters of highest contemporary relevance (“L2b” and “L4”), cgMLST analysis within the LGV lineage supported the expected high diversity of ompA-genotyped L1 or L1-like genomes, as previously observed ^5,29. Indeed, despite they predominantly grouped into a single monophyletic clade (genogroup LGV-C3-C3; n = 15), other two largely divergent genogroups were observed: i) LGV-C8-C8 with three L1 sequences, including the historical prototype strain L1/440; and ii) LGV-C6-C6 with 6 genomes presenting a hybrid L1/L2 ompA-genotype from South Africa ^5,102 (Fig. 1).

Finally, regarding the analysis of potentially epidemiologically relevant genogroups within non-LGV lineages, this study was limited not only by scarce metadata (e.g., missing sampling dates for 497 isolates), which hinders temporal resolution, but also by the paucity of confidently dated contemporary genomes. Indeed, the dataset only included 73 non-LGV genomes dated from the last decade (since 2015), comprising one from China, sequences from a recent study in the Fiji Islands ^51,103, and the new genomes generated in Portugal. However, it was still possible to extract some insights that shed light on strains that have recently circulated or are currently in circulation. For instance,the recent non-LGV Portuguese isolates were split into two singletons and four clusters at the low-level genogroup assignment. Besides the genogroup described above that harbored two recombinant “B” anorectal isolates from Portugal (T2G-C84-C152), the remaining three genogroups included strains from multiple countries, namely genogroups T1G-C20-C28 (Portugal and unknown origin; ompA-genotype Da), T2G-C66-C122 (Sweden, Portugal, and Argentina; ompA-genotype G), and T2G-C69-C125 (Sweden, Portugal, and the UK; ompA-genotype J), possibly reflecting strains with international dissemination.

An additional added-value of cgMLST is its utility to rapidly explore associations between specific alleles and clinically or epidemiologically relevant traits, whether phylogenomic (e.g., lineage) or phenotypic (e.g., tissue tropism or disease). Therefore, in this study, the distribution of alleles across the 1230 C. trachomatis samples and the 846 cgMLST loci was straightforwardly assessed and evaluated by following a novel approach that maps the allele predominance and “lineage exclusive/shared” profile across the single-linkage tree generating an “allelic heatmap” (“allelic_heatmap.py”, see Material and Methods) (Fig. 2). This allowed the detection of potential homologous recombination events and enabled the streamlined investigation of genetic mosaicism of isolates of interest, such as those showing major genotype–lineage incongruences (addressed above). An interactive HTML version of this heatmap is provided as Supplementary file S7, and a complete list of loci and their allelic classification, including “lineage exclusive/shared” profile and category association is compiled in Supplementary file S8. Additionally, by developing a complementary approach (“allelic_distribution.py”, see Material and Methods) that assesses the allele distribution across any metadata variable, such as lineage (Supplementary file S9), this study also streamlined the identification of loci with category-specific alleles, like those frequently absent or displaying pseudogene-like profiles in particular lineages (explored in the next section). In practice, these innovative tools and complementary output tables enabled several levels of insight with relevance for both surveillance and research activities.

Indeed, the systematic analysis and coloring scheme of the allelic heatmap (Fig. 2, Supplementary file S7) allowed the straightforward distinction between highly conserved and variable loci, as well as the identification of all cgMLST loci with alleles classified as “exclusive” of a given lineage or “shared” between isolates of two or more lineages (Table 1 and Supplementary files S8), thereby supporting the detection of candidate targets for diagnostics or functional studies. For instance, 214 loci presented a single and exclusive allele in LGV, 49 of which being called in all samples (Table 1 and Supplementary file S8), thus representing potential robust markers for LGV identification. Moreover, 124 genes with alleles exclusive to the emergent L4 strain were detected, 120 of which fully conserved (i.e., unique allele across all L4), thus constituting potential targets of future wet- or dry-lab assays for its direct tracking.

Although distinguishing between signals of common ancestry and those of recombination being always a challenge, the comprehensive allelic heatmap (Fig. 2 and Supplementary file S7) concordantly reflected the known C. trachomatis evolutionary history, marked by ancient diversification into four lineages, followed by contemporary mixing, as discussed elsewhere ^{5,7,29,40,41,72–74,104}. Indeed, on the one hand, it clearly highlighted the presence of shared alleles between lineages that diverged from a more recent common ancestor (T2 ocular and T2 urogenital), reflecting alleles retained from that common ancestry rather than ongoing exchange. On the other hand, it supported the absence of strict barriers to genetic exchange within the species, demonstrating that homologous recombination has occurred repeatedly within and between lineages, with a higher frequency being observed among lineages infecting similar tissue niches (e.g., urogenital and anorectal) compared with the more specialized LGV and ocular lineages. As explored previously ^5,29,71–73, a distinctive genetic feature of the LGV lineage is the markedly high allele sharing with the prevalent urogenital T1 lineage compared with other lineages. In this study, the precise identification of all these shared profiles (loci and alleles) may serve as a starting point for selecting candidate genes to explore previously raised hypotheses about the ecological success of T1 strains, in which LGV-T1 shared antigenic and genetic traits could have contributed to the remarkable expansion and persistence of the T1 lineage ⁴⁰. Ultimately, it is noteworthy that recombination signals were found to be extended across nearly the entire genome (Fig. 2, Supplementary files S7 and S8), with only 20 cgMLST loci showing exclusively lineage-specific alleles (no evidence of inter-lineage exchange) (Supplementary file S8). Moreover, only 69 loci displayed inter-lineage allelic exchanges confined to T1 and T2 urogenital lineages, with no evidence of exchange across the three main disease-associated groups (LGV, ocular, genital), suggesting fully concordant tropism-specific allelic signatures.

In practice, the vast allelic database, and associated tools for rapid link of cgMLST loci/alleles to specific traits, may serve as a useful resource to identify molecular signatures for epidemiological surveillance or outbreak tracing (e.g., for the development of multiplex sequencing approaches or rapid computational classification of genomes), as well as for research activities (e.g., functional studies to identify determinants of tissue tropism, virulence, etc.).

To advance WGS-based surveillance of C. trachomatis and improve the routine analysis of its population structure and diversity, this study assembled the largest curated C. trachomatis genome dataset to date and developed a novel cgMLST schema for the species. While WGS offers transformative resolution for studying genomic diversity, evolutionary dynamics, and epidemiology, far surpassing traditional molecular methods such as ompA-genotyping or MLST ^{5,9,21,29,42,43,49,50}, its systematic use in C. trachomatis surveillance remains limited. The main challenges are wet-lab constraints associated with the obligate intracellular nature of the pathogen, uneven sampling, and the absence of standardized genome-scale typing frameworks ^21,45,49,50.

On the wet-lab side, obtaining complete C. trachomatis genomes from culture or directly from clinical material is still technically demanding in routine laboratory settings, due to the overwhelming presence of host nucleic acids in extracted DNA samples ^{49–51,111,112}. In line with these challenges, despite our efforts to pre-screen the bacterial load of hundreds of C. trachomatis-positive clinical samples by quantitative PCR, only a small fraction met the criteria for targeted enrichment WGS, and fewer than one-third of the sequenced samples yielded complete genome assemblies. Enrichment success was variable (Supplementary figure S1), and some samples with high bacterial loads above the selection threshold (≥ 10⁴ genome copies in the WGS input) still failed to generate complete genomes. This suggests that successful genome recovery is multifactorial and cannot be explained by bacterial load alone ^49–51. Technical and biological variables, such as DNA integrity, sample storage conditions, presence of PCR inhibitors, and the ratio of bacterial to host DNA all likely influence enrichment outcomes ^49–51. In addition, technical aspects of library preparation (including fragmentation size, probe design and coverage, and hybridization conditions) can directly affect capture efficiency and the proportion of on-target reads ^43,49,50. While this reinforces the need for continued community efforts to optimize and improve culture-independent WGS protocols ^51,113, our quantitative criteria for sample selection may provide useful guidance for prioritizing samples in future studies. Opportunities for optimization include refining capture probe panels to maximize coverage across diverse lineages, adjusting hybridization protocols to improve specificity, or adopting complementary pathogen enrichment strategies. Nonetheless, the new genomic data contributed by this study (57 genomes from Portugal collected between 1995 and 2023) represent an important addition, expanding the representation of European isolates in the global dataset by 17% and anorectal samples by 56%.

On the analytical side, this study advances C. trachomatis molecular surveillance through the development of a novel cgMLST schema, designed for local deployment in both research and surveillance settings, providing a standardized and scalable approach to compare genome data across laboratories and studies. Indeed, as vastly consolidated for other important bacterial pathogens ^52,53, the availability of robust, user-friendly, and portable cgMLST schemes can be pivotal for integrating WGS into routine surveillance, and our results underscore its potential to accelerate this transition by delivering an efficient and harmonized framework for genomic monitoring and research. Taking advantage of the most comprehensive C. trachomatis genome collection to date, covering the currently known species diversity, we constructed a novel cgMLST schema of 846 loci spanning chromosomal and plasmid genes. The schema was tested with commonly used clustering algorithms (single-linkage hierarchical clustering and GrapeTree ^82,83, which showed strong overall cluster congruence, demonstrating its flexibility for implementation in commonly used bacterial genomic surveillance workflows. Using the 1230 genomes analysed here, cgMLST recapitulated the expected four deep evolutionary lineages of C. trachomatis (Figs. 1 and 2) ^5,29,41, enabling a straightforward identification of multiple instances where the ompA-genotype was incongruent with the ancestral phylogenomic lineage (Fig. 3). The most striking and previously undocumented example was that of an “E” strain clustering within the LGV lineage, with the in-depth allelic profiling (Fig. 2) revealing multiple large regions imported from a T1 lineage strain. At the same time, cgMLST achieved very high resolution to resolve intra-lineage diversity and detect mosaic profiles, thereby extending and refining beyond previous large-scale SNP-based evolutionary analysis ^51,54. Importantly, the novel cgMLST design includes a proposal of hierarchical genogroup nomenclature, which might streamline a future harmonized tracking of major circulating genogroups across laboratories and regions. For instance, within the LGV lineage, cgMLST robustly grouped the epidemic L2b clone into a genetically homogeneous cluster (genogroup LGV-C1-C1), consistent with its clonal expansion across MSM networks in Europe since the early 2000s ^5,43,100. At the same time, cgMLST clearly segregated the recently described L4 strain into a well-supported cluster (genogroup LGV-C5-C4), in agreement with recent genomic reports ^46,113. This demonstrates the cgMLST ability to discriminate and capture microevolution-driven divergence within epidemic clones and demarcate emergent lineages, illustrating its immediate translational value for genomic surveillance. Of note, these observations relied on comparisons with typical core SNP analysis (mapping against a reference genome, followed by maximum likelihood phylogenetics), which confirmed that both genome-scale typing approaches (novel cgMLST and core SNP) achieve comparable resolution. Despite large-scale C. trachomatis genomic analyses having, thus far, been almost exclusively performed with SNP-based approaches ^5,29,50,54, their routine application to large and diverse datasets is computationally demanding (especially if reference-based mapping is used) and typically does not deliver a comprehensive or readily exploitable view of the recombination landscape. In contrast, as demonstrated in the present study, cgMLST not only ensures high discriminatory power, but also inherently buffers the noise introduced by recombination, while still being able to more effectively capture its signal, and thereby providing a clearer and more consistent framework for routine phylogenomic analysis of C. trachomatis. Notwithstanding, it should be noted that cgMLST can face limitations in resolving the direction of evolution (or transmission route) among highly similar genomes. Unlike SNP-based phylogenies, which can infer ancestry and evolutionary relationships, the rooted-ancestry–independent cgMLST framework is less suited for such investigatory fine-scale resolution purposes.

Beyond phylogenetic resolution, cgMLST provided a versatile framework to link allelic variation with traits of epidemiological and research relevance. By implementing an innovative allelic heatmap approach across the 846 loci and 1230 samples (Fig. 2, and Supplementary files S7 and S8), we were able to comprehensively explore allele distributions. Although distinguishing between signals of common ancestry and those of recombination is always challenging, our simple yet comprehensive strategy allowed to identify lineage-specific and shared alleles, thereby pinpointing potential recombination events, confirming the existence of loci disrupted or pseudogenized in specific lineages, and highlighting candidate molecular signatures for lineage or disease-group. Also, all genomes were screened for the presence of a 74 bp insertion located upstream of CT_105, confirming its 100% LGV-specificity (Fig. 1) and potential suitability for the development of rapid LGV/non-LGV discriminatory tests, in line with previous observations by multiplex NGS of hundreds of clinical samples ⁵⁸. This in-depth exploration of multiple C. trachomatis genomic fingerprints, linking loci and alleles to traits of epidemiological and research relevance, illustrates how cgMLST can identify robust molecular signatures for both practical surveillance applications (e.g., development of diagnostic panels, rapid genome classification, outbreak tracing) and research directions (e.g., functional studies on determinants of tissue tropism or virulence). For instance, we exemplified that this innovative and flexible approach can be extended at different resolution levels, like at genogroup level, as demonstrated by the comprehensive identification of 124 loci with alleles exclusive to the emergent L4 strain (genogroup LGV-C5-C4), which represent promising candidates for the development of L4-specific wet- or dry-lab assays. This need was further underscored by the observation that a 9 bp insertion in pmpH that has been used to specifically detect the epidemic L2b strain ¹⁰ would most likely yield positive results also for the emergent L4 strain, potentially jeopardizing their separate tracking in countries relying solely on rapid tests for LGV molecular surveillance. This, together with a currently low (and potentially declining) adoption of ompA genotyping and MLST, further corroborate that public health surveillance laboratories should invest in setting routine WGS as a critical asset for genomic epidemiology of this important pathogen.

In summary, this work underscores the added value of the developed cgMLST for C. trachomatis genomics, providing a portable and standardized schema that preserves deep lineage structure, while offering enhanced resolution to capture sublineages, recombination, and clinically relevant clonal expansions. By integrating seamlessly with tools, such as chewBBACA ⁵⁶, ReporTree ⁸² and the newly developed scripts, the schema enables automated clustering for prospective tracking, while promoting inter-laboratory comparability. More broadly, the allelic database and associated analytical outputs, including heatmaps and searchable tables, create a rich resource for rapid link of genomic loci to epidemiological traits, bridging research and public health applications. Ultimately, this study delivers both methodological innovation and actionable genomic resources, paving the way for more harmonized, high-resolution WGS-based surveillance of C. trachomatis worldwide.