Resequencing and characterization of the first Corynebacterium pseudotuberculosis genome isolated from camel

Samples, quality assessment, and taxonomy

The genome sequences of 142 C. pseudotuberculosis strains were obtained from the NCBI RefSeq Database (https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=1719). From those, we removed 10 mutant strains (SigH, SigmaE, SigM, sigC, T1, Cp13, sigB, SigD, phoP, SigK) and strains 1002 and DSM 20689 due to those being the same strains as 1002B and ATCC19410, respectively. A total of 130 strains remained for downstream analysis (Data S1). The genome assemblies in fasta, gbff and gff format were retrieved using NCBI Datasets v15.6.1 (https://github.com/ncbi/datasets) with an input file containing the list of RefSeq assembly accession numbers. The current assembly of Cp162 is on version 3. The first and second assemblies of Cp162 (GCF_000265545.1 and GCF_000265545.2) were added for completeness of the dataset, but they were not included in the species analysis due to potential sequencing errors and misassemblies. CheckM2 v1.0.2 (https://github.com/chklovski/CheckM2) was used to evaluate the completeness and contamination of the genome sequences, while GUNC v1.0.5 (Orakov et al., 2021) was used to identify chimeric contigs, both with standard parameters. Taxonomic classification was performed using GTDB-Tk v2.3.0 with database R214 (Chaumeil et al., 2022) with the “Classify workflow” (classify_wf) and skipping ANI screen (–skip_ani_screen).

Analysis of Cp162 genome assemblies

We compared the three versions of Cp162 assemblies for completeness and contamination, size, gene content, and synteny. The number of genes was collected from the GBFF file of each genome. Differences in gene content were identified using Panaroo v1.3.3 (Tonkin-Hill et al., 2020) for gene clustering, with the parameters “–remove-invalid-genes” and “–clean-mode strict”. An in-house script was used for identifying exclusive genes (Data S2). Synteny was verified using Mauve v20150226 (Darling et al., 2004) with progressiveMauve algorithm.

Genome characterization was performed on the latest genome assembly (GCF_000265545.3). For mobile elements, prophages were predicted using the online tool PHASTER (https://phaster.ca/) (Arndt et al., 2016), while Genomic Islands (GEIs) were predicted using GIPSy v1.1.3 (Soares et al., 2016) and C. glutamicum (NZ_CP025533.1) as a reference genome. Virulence genes were predicted using PanViTa v1.1.3 (Rodrigues et al., 2023) with the VFDB database (-vfdb) (Liu et al., 2022). As PanViTa uses the core dataset of VFDB, we modified it to use the full dataset (Data S3) Antimicrobial resistance genes were predicted using PanViTa with the CARD database (-card) (Alcock et al., 2023). CRISPR-Cas systems were predicted using the online tool CRISPRCasFinder v1.1.2 (https://crisprcas.i2bc.paris-saclay.fr/CrisprCasFinder/Index) (Couvin et al., 2018) with standard parameters. A circular map comparing the three assemblies of Cp162 was built using BRIG v0.95 (Alikhan et al., 2011). The GBFF file was used as input for PHASTER, GIPSy and PanViTa, while the nucleotide fasta file was used as input for Mauve and BRIG.

Species-level analysis

To build a phylogenomic tree, we used Panaroo to identify the shared protein-coding genes across the 130 C. pseudotuberculosis isolates, and the outgroup C. ulcerans NCTC 7910 (GCF_900187135.1) and perform a multiple sequence alignment (MSA) using MAFFT (Katoh et al., 2005). The parameters used were “–remove-invalid-genes–clean-mode strict -a core–core_threshold 0.95–aligner mafft”. The phylogenetic inference was built from the MSA using IQ-TREE2 v2.0.7 (Minh et al., 2020) with standard parameters, which includes the maximum-likelihood (ML) method in which the best-fit model of nucleotide substitution selected automatically by ModelFinder (Kalyaanamoorthy et al., 2017). Support values were calculated using ultrafast bootstrap approximation with 1,000 replicates (-B 1000) (Hoang et al., 2018). The tree was visualized and annotated using the online tool Interactive Tree of Life (iTOL) v6.8 (https://itol.embl.de/).

A pangenome is a set of non-redundant genes that composes the repertoire of all genomes of a species (Tettelin et al., 2005). The pangenome and distribution of genes across all the 130 strains were identified using Panaroo due to its feature of “refinding” genes that were not annotated due to annotation artifacts. In this software, the genes are classified by frequency in the categories core genes (99% ≤ strains ≤ 100%), softcore genes (95% ≤ strains < 99%), shell genes (15% ≤ strains < 95%), cloud genes (0% ≤ strains < 15%) and total genes (0% ≤ strains ≤ 100%) (Tonkin-Hill et al., 2020). The pangenome development was calculated using Heap’s Law formula, implemented in the R package Micropan v2.1 (https://github.com/larssnip/micropan) to estimate whether it is an open or closed pangenome using 10,000 permutations (nper = 10,000). The gene clusters were annotated using eggNOG-mapper v2.1.9 with database v5.0.2 (Huerta-Cepas et al., 2017) with standard parameters. The presence of virulence and antimicrobial resistance genes were verified using PanViTa with VFDB and CARD databases.

Quality assessment and taxonomy

All genomes were classified as C. pseudotuberculosis. Completeness and contamination ranged between 97.37% and 100% and 0.14% and 6.43%, respectively. No evidence of chimerism was detected by GUNC (Data S1).

Analysis of Cp162 genome assemblies

The comparisons between the three assembly versions showed increased genome size and the number of coding sequences (CDS) (Table 1). Across all three versions, we identified 2,128 genes, 2,010 of them shared. About virulence genes, the first assembly version lacks the virulence genes nanH while spaC is fragmented as CP162_RS09080 and CP162_RS09085). The synteny analysis showed two inversions and one rearrangement in the first version, which were later corrected in the subsequent versions using the optical map (Fig. 1). Since the first assembly, the genome had an increase of 88 Kbp and 121 CDSs and a decrease of pseudogenes from 139 to 53. The third version has 18 exclusively present genes (Table 1, Data S4).

In the third assembly version, we predicted one incomplete prophage with 8.9 Kb and 18 CDSs (Data S5, S6 and S7), 13 GEIs, five of them pathogenicity islands (Fig. 2, Data S7), three CRISPR arrays, and a Type I-E CRISPR-Cas system (Data S7). PanViTa identified 30 virulence genes. The first and second assembly versions had strB and sapA as exclusive virulence genes, respectively (Table 1, Data S8). The analysis of Mauve alignment revealed that the CDSs classified as srtB (v1 locus_tag: CP162_RS09100, old_locus_tag: Cp162_1849) and sapA (v2 locus_tag: CP162_RS09670, old_locus_tag: Cp162_1968) were present in the three assembly versions, but with variations in size. PanViTa identified two antimicrobial resistance proteins: Rifampin-resistant beta-subunit of RNA polymerase (rpoB2, WP_041481489.1) and RbpA bacterial RNA polymerase-binding protein (rpbA, WP_014800420.1) (Table 1, Data S8).

Phylogeny and pangenome

A tree was generated by IQ-TREE2 using core genome alignment from Panaroo and MAFFT. The species tree in Fig. 3 shows two main clades. The first one contains a subclade composed of strains Cp162 (camel), G1 (alpaca), and I37 (cow, Israel) and another subclade containing strains isolated from horses and buffalo, separated by a host. The second one contains strain 262 (cow, Belgium) and all biovar ovis strains (collapsed) as its sister group. In the pangenome analysis, 2,332 genes were identified in the pangenome (α = 1.27), 1,877 as core, 68 as softcore, 173 as shell, and 214 as cloud. Of 2,332 genes, 2,181 were scanned by eggNOG-mapper, and 1,953 had a functional annotation (Data S9). The PanViTa analysis showed the presence of the virulence genes tufA, DIP_RS14950, mprA in most genomes and the diphtheria toxin gene (tox) only in genomes from strains isolated from buffalo (strains 31, 32, 33, 34, 35, 36, 38, 39 and 48). The antimicrobial resistance genes rpoB2 and rbpA were present in all genomes (Data S8).

Exclusive genes of Cp162

From the pangenome analysis, we also identified four proteins exclusively present in Cp162 and three exclusively absent in this lineage (Table 2). In the exclusively present group, two were predicted as the same transposase (WP_048653436.1). The other two are truncated (41 and 45 aa), none showed conserved domains, and one is in GEI 6. Analysis against the GenBank database using BLASTp and WP_275060758.1 as a query showed 92% of coverage with 52% of identity to an ATP-dependent helicase of Streptomyces spp. The same analysis with WP_231131458.1 showed 97.8% of identity with other truncated proteins from CpE19_1664 (AKS14002.1).

From the group of proteins exclusively absent in Cp162, one was recognized by eggNOG-mapper as a Transcriptional Regulator named lysG (COG category: K), another as an enzyme from NUDIX superfamily (COG category: L), and the last one as a hypothetical protein with no domains.