The complete genome sequence of C. difficile was retrieved from the GenBank database (GenBank) (https://www.ncbi.nlm.nih.gov/genbank/). Both plasmids and chromosomes were considered for separate analyses. Supplementary Table 1 provides detailed information, including accession numbers, BioProject numbers, and metadata for the sequences.

Both chromosomal and plasmid sequences were annotated using Prokka (version 1.14.6), a rapid annotation tool for prokaryotic genomes (Prokka GitHub) (Seemann, 2014). Prokka uses comprehensive pipelines and automated annotation to identify coding sequences, RNA genes, and other genomic features. Default settings were used, and annotations were verified against known databases to ensure accuracy and consistency.

MLST was conducted using the command-line tool MLST (version 2.19.0) (MLST GitHub) (Jolley and Maiden, 2010). This tool characterizes isolates based on internal fragment sequences of multiple housekeeping genes. The resulting allelic profiles were used to classify each C. difficile strain into specific STs to improve our understanding of genetic diversity and epidemiological patterns.

Pangenome analysis and clustering were performed using Roary (version 3.13.0), a pipeline for analyzing prokaryotic pangenome (Roary GitHub) (Page et al., 2015). Roary determines the presence or absence of genes across different strains and helps define core and accessory genomes. Default parameters were applied, and the results were visualized to illustrate the genomic diversity of the C. difficile strains.

The identification of antibiotic resistance genes and virulence factors, including specific C. difficile-associated toxins such as tcdA, tcdB, tcdC, tcdE, tcdR, cdtA, cdtB, and cdtR, was achieved through an integrated approach using both established and custom databases. The Comprehensive Antimicrobial Resistance Database (CARD) (https://card.mcmaster.ca/home) (version 3.2.5) (https://card.mcmaster.ca/home) (Jia et al., 2017) was used to detect resistance genes. In contrast, the Virulence Factor Database (VFDB, version 2023) (http://www.mgc.ac.cn/VFs/) (Chen et al., 2016) is a primary source of virulence genes. In addition, custom datasets were meticulously crafted to enhance the precision of detecting the critical C. difficile toxins. The ABRicate software tool was instrumental in this analysis, enabling thorough screening of genomic sequences and ensuring accurate identification of resistance genes and virulence factors, including specific toxins of interest (Chen et al., 2016).

The text with hyperlinks added to each package or software application is described as follows. To analyze the genetic data of C. difficile and examine the relationship between sequence length and year of isolation, we used the R package (version 4.3.1). Essential packages facilitate data handling, visualization, and sequence analysis. Pearson’s correlation coefficients were calculated using the Hmisc package (version 5.0-0) to assess trends between chromosomal sequence lengths and collection dates, and ggplot2 (version 3.4.0) was used to create scatter plots and trend lines. The repel package (version 0.9.1) improves plot readability by preventing label overlap. Data manipulation was streamlined using dpyr (version 1.1.2) and read (version 2.1.4), whereas the tidyverse (version 2.0.0) ensured a cohesive analysis workflow. Biostrings (version 2.70.1) provided robust tools for biological sequence operations, and Entrez (version 1.2.3) enabled efficient genomic data retrieval from NCBI. Complex datasets were visualized using Pheatmap (version 1.0.12) for heat maps and PHYLOViZ Version 2 to analyze and plot MLST data with associated resistance genes.

We analyzed 19,711 whole-genome sequences of C. difficile, including 207 isolates at the complete genome level, 39 at the chromosome level, 1,918 at the scaffold level, and 17,547 at the contig level. On average, each isolate contained 3,915 annotated genes, including 3,778 protein-coding genes and 65 pseudogenes, with an average genome length of 4.2 Mbp. The N50 values were 272 Kbp for contigs and 657 Kbp for scaffolds. These isolates, spanning 56 countries from 2000 to 2024, were predominantly human-derived (16,303 isolates), with the remainder from non-human hosts (3,407 isolates). Detailed genomic data can be found in Supplementary Data Sheet 1 .

The analysis identified 366 distinct sequence types (STs) among 19,654 C. difficile isolates, reflecting significant genetic diversity. The most common STs were ST1 (1,326 isolates), ST2 (1,141 isolates), ST11 (893 isolates), and ST42 (763 isolates). North America, particularly the United States, exhibited the highest ST diversity with 195 unique STs, Europe with 125 STs, Asia with 105 STs, Oceania with 55 STs, and South America with 28 STs. Detailed geographical distributions are provided in Supplementary Data Sheet 1 , with visual representations in Figure 1 .

The temporal analysis over 23 years based on 11,062 isolates revealed a consistent prevalence of ST1 and ST11, with an increasing prominence of ST11 post-2010. The frequencies of other STs, such as ST37, ST42, and ST54, increased during the early 2010s, indicating shifts in the genetic landscape likely influenced by clinical practices, antibiotic use, and environmental factors. The data suggest a trend towards increasing ST diversity from 2010 onwards, highlighting a broad genetic pool within C. difficile populations ( Figure 2 ).

The minimum-spanning tree (MST) visualization of sequence types (STs) and their geographic distribution, generated using PHYLOViZ analysis, provides a detailed overview of the genetic relationships and global spread of C. difficile isolates. As shown in Figure 3 , each node represents a distinct ST, with the size of the node proportional to the number of isolates. For instance, ST1 and ST2, the largest nodes, indicate a higher number of isolates, while smaller nodes represent less prevalent STs. The connections between nodes demonstrate the genetic relatedness of the STs, with closely related STs clustered together. Notably, ST1 is found to be widely distributed across Europe, North America, and Asia, while ST3 is predominantly present in Europe. The pie charts within each node provide additional insights into the geographic distribution, with countries like the United States, France, and the UK showing significant representation of specific STs. This MST analysis highlights the widespread presence of certain STs and their clustering based on genetic similarity, offering valuable information on the global geographic spread and diversity of C. difficile strains.

Further analysis of sequence metrics across different STs showed a significant inverse relationship between chromosomal sequence length and collection year, with longer sequences, such as ST947 and ST54, possibly conferring adaptive advantages ( Figure 4A ). The total gene count ( Figure 4B ), protein-coding gene count ( Figure 4C ), and pseudogene count also exhibited negative correlations with collection year ( Figure 4D ), supporting the hypothesis of genome contraction and optimization in response to environmental pressures. Supplementary Figure 1 provides additional details on sequence type length trends.

Pan-genome analysis of 239 C. difficile isolates, specifically those with assembly levels classified as complete genomes or chromosomes, revealed substantial genetic diversity and adaptability within the species. The pan-genome comprised 16,329 genes, which were categorized into four distinct groups: the core genome, containing 1,670 genes found in 99-100% of strains; the software genome, comprising 517 genes present in 95-99% of strains; the shell genome, comprising 2,278 genes occurring in 15-95% of strains; and the cloud genome, containing 11,864 genes found in 15% of strains. The core genome is responsible for essential cellular functions, which are conserved across all strains. In contrast, shell and cloud genomes highlight the remarkable capacity of these species to adapt to diverse environments and ecological niches. Supplementary Figure 2 illustrates this distribution, emphasizing the extensive variability in non-core genes, likely contributing to the bacterium’s evolutionary flexibility and environmental success.

A comprehensive analysis of 19,711 C. difficile isolates revealed widespread and emerging antibiotic resistance patterns. Resistance genes were categorized into three frequency classes: high-frequency (>10% of isolates), medium-frequency (1-10% of isolates), and low-frequency (<1% of isolates) ( Figure 5A ). High-frequency genes, such as cdeA (99.46% of isolates), cplR (99.63%), and nimB (99.67%), were nearly ubiquitous, suggesting they confer significant survival advantages. Other notable high-frequency genes included bla _CDD-1 (69.27%) and bla _CDD-2 (30.32%). Genes associated with tetracycline and macrolide resistance, such as tet(M) (18.62%) and ermB (19.09%), indicated well-established resistance mechanisms.

Medium-frequency genes, like ermQ (0.2%) and aph(2’’)-IIa (0.18%), suggested emerging resistance, with tet(O) (0.39%) and cfr(B) (0.73%) indicating the potential for increased prevalence under selective pressure. While rare, low-frequency genes, such as blaCTX-M-2, cfxA3, cfxA6, and blaZ (each found in 0.005% of isolates), could pose future risks under certain conditions.

The antibiotic resistance (AMR) gene distribution between human and non-human hosts revealed shared and host-specific patterns ( Figure 5B ). Common AMR genes, such as ccd-1, cdeA, cplR, and nimB, were prevalent in both host groups, indicating their critical role in C. difficile survival across diverse environments. Host-specific genes included aac(6’)-Ii, dfrE, efmA, and rpoB2, found exclusively in non-human hosts at low prevalence (<0.1%). Conversely, tet(T), vanR, and vanS from the vanA cluster were detected only in human isolates, potentially linked to human-specific antibiotic exposure. Genes associated with reduced vancomycin susceptibility, particularly those from the vanG cluster, were found in both host groups at low frequencies (<0.1%), with a slight predominance in human hosts (24 isolates) compared to non-human hosts (five isolates), raising concerns about the potential for the spread of reduced vancomycin susceptibility. Additionally, tet(M) was more prevalent in non-human hosts (6.0%) than in human hosts (4.7%), suggesting higher selection pressure in non-human environments.

Gene presence patterns (GPPs) provided insights into the distribution of genes associated with resistance to critical antibiotics like metronidazole, vancomycin, and tigecycline ( Supplementary Figure 3 ). Metronidazole resistance, mediated by nim genes, was observed in GPP 471 (associated with ST11) and GPP 462 (linked to ST48), with the presence of nimB indicating potential emerging resistance. Vancomycin resistance, mediated by van gene clusters, was identified in GPP 40 (associated with STs 15, 1, and 651) and GPP 173 (associated with ST3). Despite their low frequency, multiple van genes suggest robust resistance mechanisms in these patterns. The presence of van genes in isolates from 2007 to 2019 in the United States and the United Kingdom and their association with STs, such as ST1, ST3, ST21, ST63, and ST43, highlight the potential for the spread of reduced vancomycin susceptibility.

A broader analysis highlighted key antibiotic resistance patterns across the dataset. GPP 1, present in 50.2% of the isolates and associated with ST2, ST42, and ST8, was positive for cdeA, cplR, nimB, and bla _CDD-1. GPP 2, found in 17.47% of the isolates and linked to STs1, ST55, and ST37, contained tet(M), ermB, bla _CDD-2, and aph(3’)-IIIa. GPP 3, observed in 3.48% of the isolates and associated with ST1, ST15, and ST54, included the genes catA1, aad(6), sat-4, tet(O), and ant(6)-Ia, demonstrating the prevalence and combinations of specific resistance genes within the C. difficile population and highlighting the complex landscape of antimicrobial resistance.

Temporal analysis ( Figure 5C ) revealed significant shifts in the prevalence of antibiotic-resistance genes over time. Genes like cdeA and nimB showed increased prevalence, indicating a growing resistance trend. However, regression analysis ( Figure 5D ) suggested that the overall increase in AMR genes per isolate was not statistically significant (R² = 0, p = 0.36), reflecting the dynamic nature of resistance evolution. Further analysis ( Supplementary Figure 4 ) revealed that 97% of antibiotic resistance in C. difficile was chromosomal, with only 3% linked to plasmid-encoded genes. Chromosomal resistance was dominated by antibiotic inactivation (38.9%) and efflux pump mechanisms (38.1%), while on plasmids, efflux mechanisms (43.9%) and inactivation (24.6%) were significant, illustrating the bacteria’s capacity for rapid adaptation and resilience against antibiotics.

An analysis of C. difficile isolates revealed significant variability in the prevalence of specific virulence factor genes ( Figure 6A ). Highly prevalent genes, including zmp1 (99.62%), groEL (99.60%), and rpoB/rpoB2 (99.60%), were nearly universally present across the isolates. Other critical genes, such as CD0873 (99.57%), CD2831 (99.52%), and cwp84 (99.51%), also showed high prevalence, underscoring their role in the pathogenicity and survival of C. difficile. Moderately prevalent genes, like cwp66 (54.61%) and slpA (79.02%), exhibited a less consistent presence. Rare virulence factors, such as cwpV, iapC, and motA, were found in only 0.01% of the isolates, suggesting strain-specific or condition-dependent occurrence, highlighting the need for ongoing monitoring to understand their roles in pathogenicity.

The distribution of virulence genes between human and non-human hosts revealed similarities and differences ( Figure 6B ). Common virulence genes like CD0873, CD2831, cwp84, and zmp1 were prevalent in 11% of isolates from both host types, indicating their essential role in C. difficile virulence across environments. The gene cbpA had a slightly higher prevalence in human isolates (9.1%) than in non-human isolates (8.5%), suggesting a broad distribution across host species. Some genes, like cwp66, were more common in non-human isolates (6.7%) than humans (6.1%). In contrast, others, such as gap and bp, were more prevalent in non-human hosts (4.1% vs 2.8% in humans), indicating possible environmental pressures or host-pathogen interactions. Certain genes, including acm, ahpC, pilH, pilX, and sgrA, were exclusive to non-human hosts but at very low prevalence (<0.1%). At the same time, bvrR and motA were found only in human isolates, suggesting host-specific virulence factors critical for C. difficile adaptation.

Trend analysis of virulence factor gene prevalence ( Figure 6C ) indicated a slight but consistent decline in the average number of virulence factor genes per C. difficile isolate from 2000 to 2024. The regression analysis ( Figure 6D ) suggested a negative trend (y = 119.2 - 0.06x), with a moderate correlation (R² = 0.32) but insignificant results (p = 0.11). This decline may reflect shifting evolutionary pressures on C. difficile, possibly due to changes in clinical practices, antibiotic use, and environmental factors influencing virulence trait selection.

The distribution of virulence factors genes presence patterns (GPPs) across C. difficile isolates ( Supplementary Figure 5 ) revealed associations with specific sequence types (STs). The most prevalent pattern, GPP 1, included 37.43% of isolates, with eight genes commonly associated with ST2, ST42, and ST8, suggesting that these STs may dominate certain environments or host populations, contributing to widespread virulence. GPP 2, representing 23.96% of isolates, included nine positive genes linked to ST3, ST54, and ST35, highlighting their significant role in C. difficile virulence.

Other GPPs, such as GPP 3 (11.88%), GPP 4 (7.79%), and GPP 5 (4.29%), demonstrated the genetic diversity and complexity within C. difficile populations, suggesting multiple virulence strategies employed by different strains. Rare patterns like GPP 111 and GPP 112, representing only 0.01% of isolates, involved many positive genes, indicating specialized or emerging strains with unique virulence profiles, potentially marking new evolutionary paths or niche adaptations.

Further analysis identified various virulence genes, including tcdE, tcdB, tcdR, and tcdA, in isolates from countries collected between 2001 and 2023, associated with multiple STs. Among the isolates, 95.97% were toxigenic, possessing one or both toxin genes (tcdA and tcdB), while 4.03% were non-toxigenic. Specific STs, including ST15, ST747, ST3, ST83, and ST26, consistently tested positive for all genes except tcdC, suggesting these STs are likely highly toxigenic or pathogenic. These findings underscore the critical role of certain STs and gene presence patterns in shaping the virulence potential of C. difficile populations.

An analysis of toxin gene prevalence in C. difficile isolates revealed high occurrences of several key virulence factors ( Figure 7A ). The tcdE gene was the most prevalent, found in 91.26% of isolates (n = 17,989), followed by tcdC (89.67%, n = 17,675) and tcdB (89.06%, n = 17,555). These genes are widespread across the C. difficile population, underscoring their critical roles in pathogenicity. The regulatory gene tcdR was present in 86.20% of isolates (n = 16,990), while the major toxin gene tcdA was present in 64.19% of isolates (n = 12,653), indicating a somewhat lower but still significant prevalence.

In addition to the tcd genes, binary toxin genes were also assessed. The cdtR gene was found in 82.59% of isolates (n = 16,280), whereas cdtB and cdtA were less prevalent, found in 26.26% (n = 5,176) and 26.19% (n = 5,163) of isolates, respectively. These findings, depicted in Figure 7B , highlight the varying prevalence of different toxin genes within C. difficile populations, reflecting the bacterium’s diverse virulence strategies.

The host-specific distribution of toxin genes revealed distinct patterns between human and non-human hosts ( Figure 7C ). Binary toxin genes cdtA and cdtB were significantly more prevalent in non-human hosts (41.83% and 41.97%, respectively) than in human hosts (25.77% and 25.83%), suggesting these genes are more commonly retained in strains infecting non-human hosts, potentially due to differing environmental or selective pressures. Primary toxin genes, including tcdA, tcdB, tcdC, tcdE, and tcdR, were more prevalent in human isolates. Specifically, tcdA was present in 73.68% of human isolates and 48.07% of non-human isolates. tcdB was detected in 97.23% of human and 92.54% of non-human isolates, while tcdC and tcdE were more prevalent in human hosts (97.74% and 99.19%, respectively) compared to non-human hosts (93.93% and 97.12%). The regulatory gene cdtR was slightly more prevalent in human-associated strains (89.82%) compared to non-human isolates (87.59%).

Toxin gene presence patterns (GPPs) in C. difficile isolates, as shown in Supplementary Figure 6 , revealed significant diversity in toxin gene profiles across different sequence types (STs). The most prevalent pattern, GPP1, observed in 42.78% of isolates, was associated with ST2, ST42, and ST8, and included six positive genes (tcdR, tcdB, tcdE, tcdA, cdtR, and cdtB), indicating that these STs harbor a comprehensive set of virulence factors contributing significantly to C. difficile pathogenicity. Other notable patterns, such as GPP 2 (16.38% of isolates) and GPP 3 (11.17%), involved different STs but also displayed a robust array of toxin genes. Conversely, rare patterns like GPP 10 (0.68%) and GPP 11 (0.44%), as well as GPP 70, GPP 69, and GPP 68, each representing only 0.01% of isolates, highlighted the genetic variability and potential niche specialization within the C. difficile population. This variation in gene expression across different GPPs underscores the complex evolutionary dynamics, with some patterns becoming dominant while others remain rare, potentially indicating emerging or declining strains.

Trends in toxin gene prevalence over time were also examined ( Figure 7C ). The tcdE, tcdC, tcdB, and tcdR genes exhibited relatively stable prevalence rates, consistently appearing in a high percentage of isolates throughout the study period. In contrast, tcdA showed a slight increase in prevalence, suggesting its rising presence in recent C. difficile strains. Interestingly, the binary toxin genes cdtB and cdtA demonstrated a marked decline in prevalence over time, indicating that these virulence factors may be losing prominence within the C. difficile population. This decline could reflect shifts in selective pressure, possibly due to changes in environmental conditions, host immune responses, or antibiotic use patterns favoring other virulence mechanisms.

Overall, the analysis of toxin gene trends provides critical insights into the evolving virulence landscape of C. difficile. A regression analysis of the average toxin gene counts ( Figure 7D ) indicated a slight decline in toxin gene prevalence from 2000 to 2024 (y = 99.28 - 0.05x). However, the R² value of 0.01 and a p-value of 0.16 suggest this trend is not statistically significant. This implies that while there was a downward trend, the change in toxin gene prevalence over time was minimal and may not significantly impact the overall virulence profile of C. difficile populations.