A total of 788 Corynebacterium pseudotuberculosis genome sequences were analyzed in this study. Genomes were obtained from multiple sources, including publicly available sequences from GenBank (n = 233), a research isolate collection at UC Davis dating back to 1996 (Rhodes et al., 2015), and clinical isolates recovered at the UC Davis Veterinary Medical Teaching Hospital up to 2021. Additional isolates were collected at Texas A&M University between 2010 and 2017, bringing the combined total from UC Davis and Texas A&M to 380. The remaining genomes were obtained from collections in Brazil. Comprehensive metadata for all isolates is provided in Supplementary Table 1.

For UC Davis isolate collection, isolate frozen stocks (80°C) were used to inoculate vented culture tubes with 5 mL of tryptic soy broth, or brain heart infusion broth (BD Difco, Franklin Lakes, NJ, United States) supplemented with 0.2% (v/v) polysorbate 80. Cultures were incubated aerobically with shaking at 37°C for approximately 48 h. Bacteria were pelleted in microcentrifuge tubes by centrifugation at 16,000×g for 5 min, resuspended in phosphate buffered saline to wash, and pelleted again. The supernatant was removed, and pellets were stored at –80°C until DNA was extracted. DNA extraction was done with the Wizard Genomic DNA Purification kit (Promega, Madison, WI, United States). Cell pellets were lysed by suspending in Nuclei Lysis Solution from the kit, addition of 0.1 mm Mini-BeadBeater glass beads (BioSpec Products, Bartlesville, OK, United States), and bead beat for 15 s at 16,000 RPM on FastPrep-96 platform (MP Biomedicals, Santa Ana, CA, United States). Cell lysates were further processed following kit procedure. The DNA pellets were rehydrated in 100 μL of 10 mM Tris-HCl (pH 8.0). Genomic DNA quality was evaluated as described previously (Jeannotte et al., 2014; Kong et al., 2014). Briefly, purity for protein and organic contamination was assessed by Nanodrop One UV-Vis Spectrophotometer (ThermoScientific, Waltham, MA, United States), using the thresholds A260/280 ≥1.5 and A260/230 ≥1.5 for accepted samples. Genomic DNA integrity was evaluated by genomic DNA TapeStation (Agilent 4200, Santa Clara, CA, United States). The DNA was stored at –20°C until processed for Whole genome sequencing.

For the Brazilian isolate collection, genomic DNA was extracted using the Wizard^® Genomic DNA Purification Kit (Promega), according to the manufacturer’s instructions, as previously described (Sousa et al., 2025). DNA quantity and purity were assessed spectrophotometrically (260/280 ratio; NanoDrop 2000, Thermo Fisher Scientific) and verified by agarose gel electrophoresis. DNA quality was routinely assessed by spectrophotometric purity (260/280 ratio) and agarose gel electrophoresis. Concentrations were standardized across samples, and approximately 5 μg of high-quality DNA per isolate was prepared and submitted for sequencing.

Whole genome sequencing followed methods previously described for studies under the 100K Pathogen Genome Project (Weis et al., 2016; Chen et al., 2017; Weimer, 2017). Briefly, high-quality genomic DNA was used to construct sequencing libraries with 400–550 bp inserts by enzymatic sheering, followed by size selection to an average of 450 bp and sequenced to target 50x depth of coverage per genome. Sequencing was done using paired-end 150 short read method of the Illumina HiSeq X platform (San Diego, CA, United States). Raw sequence information is available on NCBI database under the 100K Pathogen Genome Bioproject (PRJNA203445).

Raw paired-end reads were trimmed to remove adapters and low-quality bases using Trimmomatic (v0.39) (Bolger et al., 2014) with the following parameters: ILLUMINACLIP:<adapters>:2:40:15, LEADING:2, TRAILING:2, SLIDINGWINDOW:4:15, and MINLEN:50. Read quality was assessed before and after trimming using FastQC (v0.11.9) (Andrews, 2010), and summary reports were compiled using MultiQC (v1.21) (Ewels et al., 2016). Potential PhiX or other sequencing contaminants were removed by aligning reads to a PhiX reference genome using Bowtie2 (v2.5.1) (Langmead and Salzberg, 2012), retaining only unmapped read pairs. De novo genome assemblies were generated using Shovill (v1.0.4) (Seemann, 2022) with default parameters and the SPAdes assembler. Assembly quality was assessed with CheckM2 (v1.0.1) (Chklovski et al., 2023), and sequencing depth was estimated using Mosdepth (v0.3.8) (Pedersen and Quinlan, 2018). Assemblies were considered high quality if they met the following criteria: sequencing coverage > 10×, estimated genome completeness above 95%, and contamination below 5%.

Genomic comparisons were conducted using Sourmash (v4.8.5) (Brown and Irber, 2016). For each genome, MinHash signatures were computed using scaled sketches (–scaled 10), and pairwise Jaccard distances (k = 31) were calculated and exported to a distance matrix. To visualize genome-wide similarity patterns, we used ComplexHeatmap (v2.21.1) (Gu et al., 2016) in R to generate a heatmap from the distance matrix, including sample annotations indicating the host and geographical locations of the genomes.

All genome assemblies were annotated using Prokka (v1.14.6) (Seemann, 2014) with default parameters. The annotated GFF files were then used as input for pangenome construction. Pangenome analysis was conducted using Roary (v3.13.0) (Page et al., 2015), which clusters orthologous genes across annotated genomes based on amino acid similarity (≥95% similarity). Separate analyses were performed for C. pseudotuberculosis biovars ovis and equi, followed by a combined analysis of all genomes. Roary was executed with the –mafft option to generate multiple sequence alignments of core genes.

The gene presence/absence matrix generated by Roary was used to investigate the gene discovery dynamics of the pangenome. We performed rarefaction analysis in R (v4.2.3) using the micropan (Snipen and Liland, 2015) package to evaluate how the number of unique genes increases as more genomes are sampled. For each of 100 bootstrap replicates, genomes were incrementally subsampled from one up to the total number of genomes (N). At each sampling step, the cumulative count of unique, non-redundant genes was recorded. The average number of genes across all replicates was then plotted against the number of genomes on a log–log scale to generate rarefaction curves according to Heaps’ law, expressed as:

where G(n) represents the number of unique genes observed in n genomes, k is a constant, and b is the gene discovery rate exponent. The model provided estimates of the exponent (b), its standard error (SE), and the coefficient of determination (R²), which were used to evaluate model fit and gene discovery dynamics. For comparative purposes, we applied the same rarefaction analysis and Heaps’ law modeling pipeline to estimate the exponent b for additional species, including Corynebacterium glutamicum, Mycobacterium tuberculosis, and Helicobacter pylori. Genomes were downloaded from GenBank and subjected to the same analysis. Results from other species (Lactococcus lactis, Pasteurella multocida, Mannheimia haemolytica, Escherichia coli, Salmonella spp., and Campylobacter spp.) were included from previous studies (Kaufman et al., 2020; Garzon et al., 2025), and are summarized in Table 1.

To investigate allelic variants of the biovar-specific petrobactin import system permease protein (fatD) from two distinct gene clusters in the Roary pangenome analysis—one exclusive to biovar equi and the other to ovis—we performed a sequence comparison. First, we extracted the corresponding locus tags for each genome from the gene_presence_absence.csv file generated by Roary (Page et al., 2015). Using the annotated gene coordinates from Prokka-generated GFF files (Seemann, 2014), we determined the contig, start, and end positions of each fatD homolog. These coordinates were formatted into a BED file and used to extract nucleotide sequences from the genome assemblies (.fna files) with SeqKit (Shen et al., 2016) via the subseq command: seqkit subseq –bed “$BED_FILE” “$FNA” > “$OUT_FASTA.” For each biovar, we created multiple sequence alignments of the extracted fatD sequences using MAFFT (Katoh and Standley, 2013), followed by consensus sequence generation with the cons command from the EMBOSS suite (Rice et al., 2000): cons -sequence ${BIOVAR}_fatD_aligned.fna -outseq ${BIOVAR}_consensus.fna.

Consensus sequences were then aligned using BLASTn to assess allelic variation. Although six nucleotide differences were identified between the ovis and equi fatD alleles, both translated into identical amino acid sequences, suggesting functional conservation across biovars.

To identify in silico virulence and antimicrobial resistance (AMR) genes, all genome assemblies were screened using ABRicate (v1.0.1) (Seemann, 2020). For virulence profiling, the VFDB was queried. For resistance screening, assemblies were compared against the CARD, NCBI AMRFinderPlus, ResFinder, ARG-ANNOT, and MEGARes databases. ABRicate was executed with default parameters for each database, and gene presence was summarized across all genomes using the –summary function.

Variant calling was performed using Snippy (v4.6.0) (Seemann, 2015). For biovar ovis isolates, sequencing reads were aligned to the C. pseudotuberculosis strain 1002B genome (Accession: GCF_001433475.1), whereas biovar equi isolates were aligned to strain 258 (Accession: GCF_000263755.5). Each sample was processed individually. Resulting outputs were merged with snippy-core to obtain a full alignment. Low-quality and gappy alignment regions were removed using snippy-clean_full_aln, and recombination events were detected and masked using Gubbins (v2.3.4) (Croucher et al., 2015). The recombination-filtered alignment was then reduced to core SNPs using SNP-sites (v2.5.1) (Page et al., 2016). A phylogenetic tree was constructed using FastTree (v2.1.11) (Price et al., 2009) under the GTR model with nucleotide input. The final tree was visualized and annotated using the Interactive Tree of Life (iTOL) (Letunic and Bork, 2021).

To identify single-nucleotide polymorphisms (SNPs) associated with host specificity in C. pseudotuberculosis, machine learning models were developed using the XGBoost algorithm (v1.7.8.1) (Chen and Guestrin, 2016). SNP presence/absence data were formatted as binary matrices, with each dataset annotated according to host classification. Models were implemented in R using the XGBoost, caret, and dplyr packages. Multi-class classification was performed using a softmax objective (multi:softprob), with model parameters set to eta = 0.3, max_depth = 6, and nrounds = 100, incorporating early stopping to prevent overfitting. Each dataset was partitioned into training (80%) and testing (20%) subsets using stratified sampling to preserve class proportions. Model performance was evaluated using confusion matrices and classification accuracy on the held-out test data. Feature importance scores were extracted from the trained models using built-in XGBoost methods. SNPs with high importance scores were prioritized as potential markers of host specificity. The resulting ranked SNPs were compared to findings from pangenome and core-genome SNP analyses to identify convergent genomic signatures linked to ecological adaptation.

To investigate the biological relevance of top-ranked SNPs associated with host species classification, functional annotation and pathway analysis were performed. Genes linked to high-importance SNPs identified by XGBoost were mapped to metabolic and regulatory pathways using BioCyc Pathway Tools (Karp et al., 2018). In parallel, protein-protein interaction networks were reconstructed using the genes of the top-ranked SNPs in the STRING database (v11.5) (Szklarczyk et al., 2021), enabling the identification of enriched functional modules.

A total of 788 good-quality genome assemblies of C. pseudotuberculosis were assessed in this study, comprising 408 isolates classified as biovar equi and 380 isolates classified as biovar ovis. K-mer–based genomic similarity profiling, using a sketch size of 10, enabled fine-scale resolution of whole-genome comparisons and revealed patterns consistent with the underlying population structure of the species. The genomes clustered into three distinct genomic groups: two corresponding to the classical ovis and equi biovars, and a third novel equi cluster composed predominantly of camelid-derived isolates (Figure 1). Sequence conservation within each genomic cluster was remarkably high, as reflected by pairwise Jaccard similarity values: 0.868–1 for the equi cluster, 0.888–1 for the ovis cluster, and 0.970–1 among camelid-associated strains. The narrow range of pairwise Jaccard similarity values indicates high genomic homogeneity within each cluster. The majority of biovar ovis isolates were associated with small ruminants, while most biovar equi isolates originated from equine samples, consistent with historical reports of host specificity for these biovars (Costa et al., 1998; Haas et al., 2017; do Nascimento Sousa et al., 2024). A total of 25 bovine-derived isolates were identified in the dataset. Of these, 22 were assigned to biovar equi and 3 to biovar ovis. 1 assigned to the biovar equi cluster with the camelid cluster. This supports previous observations that cattle can act as incidental hosts for both biovars, and have strains similar to all genome types (Almeida et al., 2017). While most genomes were consistent with expected biovar—host associations based on epidemiological reports, several exceptions indicated instances of cross-host infection. Specifically, seven ovis biovar genomes were isolated from equine hosts, and eight equi biovar genomes were obtained from small ruminants (four from caprine and four from ovine sources). Also notable in this comparison was the inclusion of six human-derived isolates: two equi strains from Romania and four ovis strains, including three from New Zealand and one from Oslo.

While there was distinct biovar-level genomic clustering and expected species-association, no meaningful population structure was observed related to geographic origin. Limited clustering was observed by country or US state, however, a high level of genomic similarity by Jaccard index (>0.87) was observed. The origins of the isolates included globally diverse locations such as Brazil, Germany, New Zealand, Romania, Switzerland and the United States. The lack of distinct subclusters related to geographic location coupled to the monolithic genome structure indicates that an unusually highly conserved genomic gene content exists for each biovar globally.

To quantify the relationship between genome sampling and gene diversity, we applied a rarefaction analysis based on Heaps’ law, as described in the methods. This model characterizes how the number of unique genes grows with increasing genome sampling, with the Heaps’ exponent (b) serving as a metric of pangenome openness. Higher values indicate more dynamic and diverse gene repertoires, while lower values suggest greater genomic conservation. As previously observed (Tettelin et al., 2005) the discovery of new genes follows a power law model (log-log scale), where an increase in the number of genomes sequenced leads to a proportional increase in the number of genes discovered. The overall gene discovery rate for C. pseudotuberculosis was remarkably low, with a fitted exponent of b = 0.1005. This implies that approximately 10 genomes are required to identify one additional unique gene (1/b ≅ 10), reflecting a highly conserved gene repertoire despite extensive sampling (n = 788). Subgroup analyses within C. pseudotuberculosis revealed subtle but informative differences. The equi biovar cluster (n = 385) showed a comparable exponent (b = 0.101), while the ovis biovar (n = 380) was slightly more conserved (b = 0.065), requiring over 15 genomes to discover one new gene. The most conserved subgroup was the camelid-associated cluster (n = 22), with a Heaps’ exponent of just b = 0.037, indicating that nearly 27 genomes would be needed per novel gene identified. These values collectively underscore the exceptional genomic stability of C. pseudotuberculosis, particularly in host-associated subgroups. This genomic stability points to an evolutionary path that is highly constrained, reflecting long-term adaptation to specific hosts and ecological niches. Rather than gaining or losing genes frequently, C. pseudotuberculosis appears to have refined existing genetic functions to optimize survival within its preferred environments. The closed nature of its pangenome supports the view of limited horizontal gene transfer, strong purifying selection, and a largely clonal population framework. Consequently, most diversification within the species arises from subtle sequence-level mutations and metabolic fine-tuning instead of major shifts in gene content. Together, these features suggest that the bacterium’s specialized pathogenic way of life depends on a stable and well-conserved genetic toolkit that ensures effective colonization, persistence, and transmission in its hosts.

In contrast, other pathogens exhibited substantially more open pangenomes. Close relatives such as Corynebacterium diphtheriae (b = 0.225) and C. glutamicum (b = 0.235) demonstrated higher gene discovery rates, with only ∼4 genomes required per new gene. Mycobacterium tuberculosis, despite its known clonal structure, exhibited a comparable exponent (b = 0.237), suggesting greater gene content variation. Among Gram-negative pathogens of veterinary importance, Pasteurella multocida (b = 0.321) and Mannheimia haemolytica (b = 0.326) showed even higher pangenome openness (Garzon et al., 2025), requiring ∼3 genomes per novel gene. Further comparisons with additional species such as Helicobacter pylori, Escherichia coli, Salmonella spp., Lactococcus lactis, and Campylobacter spp.— which ranged from b = 0.454 to 0.645— highlighting the contrast in pangenome dynamics. These organisms typically required only 1.5–2.2 genomes to yield a new gene, consistent with highly open pangenomes. To evaluate whether gene discovery rates differ between these Gram-positive and Gram-negative bacteria, we compared their mean discovery rates using the Kruskal—Wallis rank sum test. This non-parametric test yielded a chi-squared statistic of χ² = 3.33 (p = 0.068), indicating a nonsignificant trend toward higher gene discovery rates in Gram-negative species. Although this result is not statistically significant at the α = 0.05 level, the observed pattern aligns with previous reports that Gram-negative bacteria tend to have more open pangenomes, driven by larger accessory genomes and greater inter-strain gene content variability. Among Gram-positive taxa, C. pseudotuberculosis and Lactococcus lactis are outliers, exhibiting the lowest and highest gene discovery rates, respectively. Notably, C. pseudotuberculosis stands out as an extreme case of exceptionally limited gene discovery even in comparison to other Gram-positive taxa, reinforcing its characterization as a genomically stable species.

The contrast in gene discovery became even more pronounced when compared to highly diverse bacterial species such as Escherichia coli, Salmonella spp., and Campylobacter spp., which exhibited gene discovery rates approximately 4.5–6.5-fold higher than that of C. pseudotuberculosis. These findings highlight the species’ exceptionally clonal nature and limited genetic variability, consistent with previous reports describing its evolutionary stability.

In addition to metabolic capacity, virulence factors can also contribute to biovar-specific host adaptation. To assess the presence of virulence factors, genes that support a pathogen’s ability to infect a host and initiate disease, all 788 C. pseudotuberculosis genomes were analyzed using the Virulence Factor Database (VFDB). The phospholipase D gene (pld), which encodes a potent exotoxin, was found in all isolates, consistent with its known crucial role in C. pseudotuberculosis virulence. The diphtheria toxin gene (tox) was detected only in 11 isolates from buffalo hosts, in accordance with previous reports (Viana et al., 2017). No additional virulence-associated genes were identified by VFDB across the dataset, indicating that pld and tox represent the only detectable virulence factors under the applied screening criteria.

Subsequently, the genomes were screened for antimicrobial resistance (AMR) genes using multiple curated databases. While most genomes lacked significant matches to known AMR genes under standard thresholds, a subset of isolates carried resistance determinants. Specifically, the aminoglycoside resistance gene APH(3’)-IIa and the beta-lactamase gene TEM-116 were co-detected in seven isolates (phoP, SigmaE, sigB, sigH, sigM, sigD, and sigC), all derived from goat hosts in Brazil. In contrast, APH(3’)-IIIa was detected exclusively in a single, distinct isolate (Cp13), also from a goat in Brazil. All detected resistance genes exhibited high sequence identity (≥99.9%).

Building upon the observed high conservation of C. pseudotuberculosis genomes at both the genome content and SNP levels, machine learning was used to assess whether subtle patterns of nucleotide variation could distinguish isolates according to host species between biovars. The XGBoost algorithm was applied to the core genome SNP matrices with a total of nine host categories considered. Feature importance analysis revealed that classification was not driven by a single dominant variant but rather by the cumulative effect of multiple SNPs, each contributing modestly to model performance (Table 2). Feature importance scores reflect the relative contribution of each SNP to the model’s predictive accuracy during classification, rather than direct causal effects on host specificity.

The top-ranked SNP was a missense variant (c.442G>A; p.Val148Ile) within the oppB1 gene, encoding a subunit of an ABC-type dipeptide/oligopeptide/nickel transporter. Additional important variants included a missense mutation in rnhB (an RNase H enzyme), a synonymous variant in oppC4 (another oligopeptide transporter), and a missense variant in bioD1 (involved in biotin biosynthesis) (Table 2). These SNPs were prioritized based on their importance within the trained model and collectively capture nucleotide patterns that improve discrimination among host-associated isolate groups. Individually, none of these variants uniquely defines a host species; instead, host classification emerges from the combined signal of multiple SNPs distributed across conserved metabolic and transport genes.

Functional enrichment analysis of SNP-associated genes revealed a significant overrepresentation in the “Amino Acid Biosynthesis” ontology, which is a subclass of the broader “Biosynthesis” category (Figure 5A). Enrichment of amino acid processes suggests that a subset of core metabolic functions may be differentially modulated among isolates from different hosts.

Further analysis of protein-protein interactions for the gene products of top SNPs revealed two major functional modules (Figure 5B). The first comprised a densely interconnected cluster of oligopeptide permease (Opp) system components, including oppB1, oppC4, and oppA—F paralogs. The Opp system is a well-known ATP-binding cassette (ABC) transporter involved in peptide import and environmental sensing. This system has been associated to host adaptation in other bacterial pathogens, including Streptococcus pneumoniae and Listeria monocytogenes (Borezee et al., 2000). The second cluster was composed of biosynthetic enzymes, such as bioD1 (biotin biosynthesis), argD (arginine biosynthesis), and dapC (lysine pathway).

GO enrichment (Figure 5C) further supported these findings. Statistically significant terms included “Amino Acid Biosynthetic Process,” “Proteinogenic Amino Acid Biosynthesis,” and “Modification of Amino Acids Within Proteins,” all with FDR values < 1.8 × 10^–6. These pathways emphasize metabolic functions are associated with host-adaption in each biovar.