The whole genome sequence sets from Paraclostridium spp. (i.e., P. bifermentans and P. benzoelyticum) and P. sordellii were downloaded from the NCBI genome Refseq database. The isolation sources of each strain were collected from the metadata table downloaded from the Refseq database or by searching the literature. The genomic sequences were evaluated for completeness and contamination using CheckM v1.0.12 (Parks et al., 2015). The standard lineage workflow “lineage_wf” from CheckM was performed to summarize general sequence features such as completeness, contaminations, genomic size, GC contents, coding density, and predicted rRNAs. Genomic sequences with completeness less than 70% or contamination greater than 10% were eliminated from the sequence sets. A total of 27 Paraclostridium spp. genomic sequences comprising one P. benzoelyticum strain JC272 and 26 P. bifermentans strains, and 60 P. sordellii genomic sequences were finally collected for further analysis (Supplementary Table S1).

Each genome sequence in the genome sets was annotated using Prokka v1.14.5 with default parameters. Aligned core and accessory genome sequences of the genome sets were extracted separately from the resultant gff files using the Roary v3.11.2 pangenome pipeline (Page et al., 2015). The phylogenomic tree of each alignment (i.e., core and accessory) was constructed using FastTree v2.1.11 with the generalized time-reversible model. The average nucleotide identity (ANI) between all the collected sequences was computed using fastANI v1.33 (Jain et al., 2018), and the genome distance (GD) was calculated using the phylonium v1.2 (Klötzl and Haubold, 2020). The tanglegram between the core-genome tree and accessory-genome tree was visualized using the dendextend v1.15.2 R package. The ANI matrix and GD matrix were visualized using the pheatmap v1.0.12 R package. Finally, the core-genome tree with general features was visualized using iTol v6 (Letunic and Bork, 2021).

Plasmid sequences in the genome sets were predicted using (1) PlasForest v1.2,¹ a machine learning-based tool that uses sequence homology, and (2) BLAST-based method with an e-value lower than 1e⁻⁵ and a sequence similarity higher than 90% and query coverage higher than 50% by aligning to the documented plasmids from Paraclostridium spp., P. sordellii and C. perfringens, and (3) circular contigs with no dnaA gene and no rRNA coding sequences although not predicted as plasmids based on the sequence homology methods. Pangenome features of either P. bifermentans or P. sordellii were extracted using the Roary v3.11.2 pangenome pipeline from the gff files produced by Prokka v1.14.5. Subsequently, the extracted pangenome features were analyzed and visualized using the R package Pagoo v0.3.12 (Ferrés and Iraola, 2021). The core genome was defined as the set of CDSs that are present in over 95% of genomes. All other CDSs were classified as the accessory-genome, which was further subdivided into the cloud genome (i.e., those present in less than 5% of genomes) and the shell genome (i.e., those present in 5–95% of genomes). The functional genome annotation was conducted using eggNOG-mapper v2.1.5, which comprehensively annotated the Carbohydrate-Active enZYmes (CAZy), Cluster of Orthologous Groups (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Cantalapiedra et al., 2021). Virulence factor coding genes, including toxin genes and antimicrobial resistant genes, in the pangenome were predicted using PathoFact v1.0 with the “complete” parameter (de Nies et al., 2021).

Prior to characterizing the genomic features, we carefully checked the isolation source of the strains. Interestingly, despite the ubiquitous presence of P. sordellii strains in the environment, none of the strains that had genomes sequenced are isolated from the environment, while 15 originated from animal sources and 45 from human sources. In contrast, strain P. benzoelyticum JC272 was isolated from the environment, and 6, 6, and 13 out of the 26 P. bifermentans strains were isolated from environmental, animal, and human sources, respectively (Figure 1A). Our analysis of the whole genome tree for both P. bifermentans (P.b) and P. sordellii (P.s) did not reveal any correlation between isolation source and within-species phylogeny, indicating that isolation source did not impact the phylogenetic relationships within these species (Figure 1A). Then, the collected genome assemblies were compared on assembly level, genome size, GC content, completeness, contamination, and predicted genes (Figure 1A). All of the genomes showed high completeness (> 99%), 4 out of the 87 genomes had medium contamination (5% < contamination <10%), and only 1 was assembled into more than 500 scaffolds. These features suggest that the genome set is of high quality. However, only 4 genomes were assembled at the completed level (2 from P.b and 2 from P.s), and 3 at the chromosomal level, and the others at the scaffold or contigs levels (Figure 1A). Of note, the rRNA sequence numbers showed high variability between strains within and between species. This is likely due to the limitations of short-read sequencing and assembly methods in accurately assembling the multicopy rRNA genes in bacteria. Analysis of completed genomes revealed that P.b and P.s possess a high number of rRNA genes, ranging from 48 to 51, which suggests a high protein synthesis capacity for these two species (Figure 1A). Additionally, these two species showed similar coding density, ranging from 0.84 ~ 0.87 (Figure 1B). The genome size and GC content of P.b are significantly higher than those of P.s (Figure 1B). For the following reasons, five genome sequences were excluded from further functional analyses: (1) they were identified as outliers in genome size and were assembled from metagenomes (i.e., MAG072, L3_131_244G1_dasL3_131_244G1_concoct_40_sub); (2) they were identified as belonging to a different species (i.e., strain JC272 of P. benzoelyticum); (3) they were identified as outliers in both genome size and GC content (i.e., BIOML-A1, and BIOML-A2) (Figure 1B).

As shown by the core-genome tree constructed from P.b and P.s, P. bifermentans is significantly divergent from P. sordellii (Supplementary Figure S1). The tangled dendrogram of the trees is displayed in Figure 2A. In detail, P.b and P.s are clearly separated by the topological structure of both the core-genomic and accessory-genomic tree, from which the tanglegram displayed highly consistence (entanglement = 0.11), underscoring a distinct divergence in phylogeny between P.b and P.s across both the core-genome and accessory-genome (Figure 2A). Furthermore, the heatmap and hierarchical clustering from ANI and GD further reveal the phylogenomic divergence of P.b and P.s at the species level (Figures 2B,C), highlighting a robust separation between the two species and substantiating the evolutionary divergence at both the core and accessory genomic levels.

To compare the coding capacity of P.b and P.s, we performed pangenome analyses using Roary. The pangenome of 24 P. bifermentans genomes and 58 P. sordellii genomes consisted of 10,030 and 11,076 orthologous genes, respectively. The pangenome curves for both species formed a “C” shape, and the alpha index of the heaps law models being lower than 1 (Figure 3A) (Tettelin et al., 2008), which revealed that both P.b and P.s carried open pangenomes (Figure 3A). The pangenome curves of P.b and P.s were best fitted for the positive-power model, with the formulas y P . b = 3284.7 χ 0.3451 R 2 = 0.9971 and y P . s = 2948.1 χ 0.3153 R 2 = 0.9915 , respectively, indicating that the pangenome of P.b is more plastic than that of P.s ( y P . b > y P . s ) (Figure 3A). Moreover, the core-genome size of P.b (2384) is comparatively smaller than that of the P.s (2526) (Figures 3B,C), which is further corroborated by the core-genome curves (Figure 3A). As the average genomic size of P.b is larger than that of P.s (Figure 2), the comparison in core-genome size between P.b and P.s further supported that P.b had a more plastic pangenome than that of P.s by carrying a greater number of accessory CDSs. The core-genome curve showed that 1949 orthologous genes are shared by all P.b strains, and 1912 orthologous genes are shared by all P.s strains, with the strain-specific genes in P.b and P.s are 3,673 and 4,095, respectively (Figures 3A–D). The core-genomes of P.b and P.s accounted for only 23.77% (2,384/10030) and 22.81% (2,526/11076) of the respective pangenome (Figures 3B,C). The core-genome-coded COG functions between P.b and P.s were similar in COG categories and counts, while the accessory-genome-coded COG functions showed high differences by COG counts (Figure 3D). P.b’s accessory-genome, for example, encoded more COG functions in nutrient metabolism, protein synthesis and energy conversion (e.g., “C Energy production and conversion,” “G Carbohydrate transport and metabolism,” “I Lipid transport and metabolism,” “J Translation, ribosomal structure and biogenesis,” “K Transcription,” “O Posttranslational modification, protein turnover, chaperones,” “P Inorganic ion transport and metabolism,” “T Signal transduction mechanisms” and “U Intracellular trafficking, secretion, and vesicular transport”), while that of P.s encoded more COG functions in cell cycle and host adaptation (e.g., “D Cell cycle control, cell division, chromosome partitioning,” “L Replication, recombination and repair,” “M Cell wall/membrane/envelope biogenesis,” “N Cell motility,” “S Function unknown” and “V Defense mechanisms”) (Figure 3D).

We further explored the virulence factor coding capacity of P.b and P.s to compare their pathogenic potentials. The virulence gene profiles between P.b and P.s were divergent, enabling effective discrimination between the two species (Supplementary Figure S2a). Specifically, each P.s strain carried an average of 71 ± 3 virulence genes, while each P.b strain carried an average of 70 ± 4 virulence genes. A total of 176 and 167 virulence genes were carried by the pangenomes of P.s and P.b, respectively. Among them, 29.34% (49/167) and 31.25% (55/176) belonged to the core-genomes of P.b and P.s, respectively, accounting for approximately 70% of the virulence genes in each strain (Supplementary Figure S2B, Figure 4A). Furthermore, P.s carried 39 virulence genes encoding secreted toxins in the pangenome, and 14 of them were in the core-genome, whereas that of P.b was 28, and only 9 of them were in the core-genome. Since we had confirmed in this study that P.b carried a much opener pangenome than that of P.s (Figure 3A), we can deduct from the results here that the core-genome of P.b may encode less secreted-toxin genes than that of P.s, with the available P.b genomes increased. In addition, 44.91% (75/167) of the virulence genes in P.b were from MGEs, 19 of which are core genes, and over half of the MGE-related virulence genes (40/75) are phage-originated (Figure 4A). Similar findings were observed in P.s (Supplementary Figure S2B). The MGEs, especially the phage-originated MGEs, play vital roles in the expansion of the virulence gene pool in the pangenomes of P.b and P.s (Figure 4A, Supplementary Figure S2B). Regarding virulence gene families, most of them were shared by the pangenomes of P.b and P.s, but the gene counts varied. For example, the P.b pangenome coded more Zinc-dependent phospholipase C, UDP-glucose 4-epimerase C-term subunit, PLD-like domain, and Nitroreductase family virulence factors, while that of P.s coded more toxin A/B, SpaB C-terminal domain, and enterotoxin D virulence factors (Figure 4B). Additionally, there are several species-specific toxin gene families, such as enterotoxin C, pretoxin HINT domain, and Clostridium neurotoxin translocation domain in P.b, and Zinc carboxypeptidase, and Clostridium enterotoxin in P.s (Figure 4C). A detailed virulence gene co-occurrence network between P.b and P.s is shown in Supplementary Figure S3.

We further compared the toxin gene counts derived from different sequence types to estimate the virulence coding capacity of P.b and P.s. The results showed no significant difference in the pangenome-coded toxin gene counts between P.b and P.s. However, the coregenome-coded toxin genes in P.s were higher than those in P.b (Figures 5A–C). Moreover, both chromosomal- and phage-originated toxin gene counts in the pangenome or coregenome of P.s were higher than that in P.b (Figures 5A,B), which displayed an opposite pattern in the accessory-genome (Figure 5C). It is noteworthy that the plasmid-originated toxin genes were higher in both the pangenome and coregenome of P.b than in P.s (Figures 5A,B). As a result, the plasmid-originated sequences played a significant role in shaping the toxin gene profile of P.b’s coregenome, while the phage-originated sequences played a significant role in that of the accessory-genome (Figures 5B,C).

Furthermore, we compared the presence of the common reported clostridial toxin genes in P.b and P.s. The results showed that the pangenomes carried a total of eight toxin gene homologs, out of which five, namely colA, pfo, plc, toxB/tcsL, and iap, were shared by both species. The former three genes were core genes in both species that coded for proteins involved in host cell surface attachment and cell membrane damage, indicating that P.b and P.s may share similar mechanisms for initiating infections (Table 1). The toxB/tcsL and iap that encode host cell intracellular toxicity, which may lead to lethal and dermonecrotic effects, are presented more frequently in P.s than in P.b, indicating that the P.s is more toxic than the P.b (Table 1). The remaining three genes, namely toxA/tcsH, plcB and cpe, were specifically carried by P.s, coding functions of hemorrhagic toxin, pathogen’s vacuole escape, and cell pore-forming. Furthermore, the plcB gene was a core gene in P.s but absent in P.b (Table 1).

To elucidate the role of plasmids in the virulence of P. b, the plasmid sequences in the genomes of P.b were predicted and extracted using the methods described in the methods section. Plasmid pPbm14_8 from strain Cbm was depleted in further analyses because of that the pPbm14_8 was identified as a chromosome sequence by aligning to the NCBI-nt database (Supplementary Table S2). By aligning the reported circular plasmid sequences in P.b and P.s strains, we then found that all plasmids between species are divergent, while some of those are highly conserved within-species (Supplementary Figure S4a). Further analyses on the plasmid sequence similarity between C. perfringens, P.b and P.s also revealed high within-species conservation (Supplementary Figure S4b). Finally, eight conserved plasmid groups in P.b were identified, ranging from 25 kb to198 kb in length (Figures 6A–F, Supplementary Figure S5). Among them, group_a, group_b and group_c are the most conserved groups with a presence in 25–33.33% P.b. strains, which were carried by 8, 8 and 6 P.b strains, respectively (Figures 6A–C). Plasmids belonging to group_a are the longest and two of them (i.e., pHD0315_2–1 and DSM14991 unnamed1) were assembled at the completed level (Figure 6A). The longest plasmid in group_a is from strain BSD_D6, the NODE_4 in the draft genome, which is 198,761 bp, 40 kb longer than other plasmids in this group (Figure 6A). This group of plasmids encode about 150–180 CDSs, of which one third can be annotated, and is the best annotated plasmid group in P.b (Figure 6A). The reference plasmid in this group carries eight potential virulence genes, of which five (i.e., Cupin_2, nitroreductase family, puuR, corC_2, and tlyC) are core genes in P.b (Figures 4A, 6A). It is noteworthy that the plasmids in this group are essential for regulating bacterial growth and adaptation through the encoding of various functions related to organic sulfate reduction (asrABC, nirC), 5,6-dimethylbenzimidazole synthesis (nox_2), styrene degradation (gatA), aminoacyl-tRNA biosynthesis (aspC, gatA), etc. (Figure 6A). Additionally, this group of plasmids encode five virulence genes, two of which, namely cat and drrA_2, confer resistance to chloramphenicol and doxorubicin. Group_f represents a rare but typical virulence plasmid group in P.b that is exclusively carried by two strains, namely parabai and Cbm, isolated from anopheline-endemic areas in Malaysia and Brazil, respectively. This plasmid group encodes two pesticidal toxins (Cry16Aa and Cry17Aa) and one novel neurotoxin, which are flanked by Tn3 and IS1182 family transposases, respectively (Figure 6F). Notably, the Tn3 transposase is also encoded by other plasmid groups (e.g., group_b and group_c) (Figures 6B,C). While group_b and group_e plasmids are less annotated, some of the genes they carry, such as pqqD, srtB and dtpT, are related to pyrroloquinoline quinone synthesis (Ludueña et al., 2017), surface anchoring (Weiss et al., 2004), and salt stress protection (Wouters et al., 2005). Conversely, group_d, group_g and group_h plasmids are poorly annotated (Figure 6D, Supplementary Figures S5g,h). Interestingly, strains carrying plasmids clustered in the same conservation group were isolated from distant geographic locations (Figure 6G), which may be clustered by the isolation-source and have co-evolved with their habitats, as we found that the group_a and group_b plasmids are predominantly carried by strains isolated from feces (Figure 6G). Furthermore, strains with completed genome sequences (i.e., HD0315_2 and DSM14991) carried more plasmids than those with draft genome sequences. Unfortunately, over half of the strains with genome sequences are not isolated from clinical settings (Figure 6G).