Bacterial strains used in the study were maintained in the culture collection of the Institute of Bacterial Infections and Zoonoses (IBIZ), Friedrich-Loeffler-Institut (FLI), Jena, Germany. C. chauvoei strains (60 strains) were mainly of European origin representing Germany (38 strains), Austria (ten strains), Switzerland (two strains), Italy (one strain), and few strains of unknown origins (nine strains). Within Germany, blackleg outbreaks are reported from six different states and the strains involved in the study were representative from all these reported states (Annual Animal Health Report, 2017¹). Schematic representation of blackleg outbreaks in Germany from 1995 to 2016 is shown in a geographical map in Figure 1. One strain was from Canada and another strain (NCTC 08361) from a sheep in South Africa. In addition, four previously published strains of C. chauvoei (genome sequence/reads) were included in the analysis (Falquet et al., 2013; Thomas et al., 2017). Details of strains involved in the current study with respect to epidemiological data (host, year, and country of origin) and Sequence Read Archive (SRA)/Genome assembly accession numbers are provided in Supplementary Table 1. These strains were mostly isolated from various tissues such as muscles, spleen, and liver following blackleg outbreaks. All strains were cultured on sheep blood agar plates (Yeast Extract Cysteine Medium with Sheep Blood (Beerens Formulation) Thermo Scientific, Oxoid, Germany) at 37°C for 24–48 h under anaerobic conditions. The culture material scraped off from one to two agar plates was used for DNA extraction. The genomic DNA was isolated using DNeasy Blood and Tissue Kit (Qiagen, Germany) with slight modification, as 40U of achromopeptidase (Frey et al., 2012) was included in the enzymatic buffer along with lysozyme. A PCR for species confirmation was based on published primers specific region for 16–23 S rDNA spacer regions for differentiating C. chauvoei and C. septicum (Sasaki et al., 2000). Library preparation from genomic DNA was carried out using Nextera^TM (Illumina, Netherlands) library preparation method. Genome sequencing was carried out using MiSeq^TM System (Illumina, United States) paired-end sequencing technology (2 × 300bp) at the Institute of Microbiology and Epizootics (IMT), Freie Universität Berlin, Berlin.

Paired-end Illumina reads were first checked for quality using the FastQC version 0.11.6²). Kraken2 was then used for taxonomic classification (Wood and Salzberg, 2014). For the assembly, we used Shovill version 1.0.4³ using the option “-trim” enabling adapter and quality trimming using Trimmomatic version 0.39 (Bolger et al., 2014). The two complete genome sequences (DSM 7528^T and 12S0467) previously published (Thomas et al., 2017) based on PacBio assembly, were downloaded and polished for four rounds with Illumina reads using Pilon tool version 1.22 (Walker et al., 2014). Assembly statistics were generated using QUAST (Gurevich et al., 2013). Genome annotations were carried out using Prokka version 1.11 (Seemann, 2014). The presence of previously reported plasmid (4.1Kb; NZ_CP018631.1) was searched among the strain assembly files using ABRicate version 1.0.1⁴.

CRISPR elements and spacers were predicted using the CRISPRDetect program (Biswas et al., 2016) followed by CRISPR array visualization using CRISPRStudio (Dion et al., 2018). When CRISPR regions were represented in more than one contig (six strains), the orientation and directions were assigned as per DSM 7538^T strain. Prophage elements and genomic islands were predicted using Prophage Hunter web tool (Song et al., 2019) and IslandViewer 4 (Bertelli et al., 2017), respectively. BLAST atlas of all the genomes with locations of prophages and genomic islands was created using BLAST Ring Image Generator (BRIG) (Alikhan et al., 2011).

Pangenome analysis was carried out with Panaroo pipeline version 1.2.7 (–clean-mode strict; –remove-invalid-genes) (Tonkin-Hill et al., 2020) using the Prokka (Seemann, 2014) annotated genomes. Curves for the core and pangenome were calculated and plotted using PanGP version 1.0.1 (Zhao et al., 2014). New gene discovery was plotted as a function of the number of genomes added sequentially using a distance guide algorithm with 100 replicates and 3000 permutations of genome order. PanGP uses power-law regression (y = A^xB + C) to model the pangenomes generated from all permutations, where y is the total number of gene families in the pangenome, x is the number of genomes considered, A, B, and C are fitting parameters. When 0 < B < 1, the pangenome should be considered open. The exponential curve fit model, y = Ae^Bx + C was used to fit the core genome. Here, y denotes core genome size, x denotes the number of genomes; A, B and C are fitting parameters. The curve for new gene discovery represents the least-squares fit for the function y = A^xB, where y denotes the new genes, x denotes the number of genomes and A and B are the fitting parameters.

Single nucleotide polymorphism analysis was carried out using Snippy version 3.2⁵. The pipeline maps read-pairs to reference strain (DSM 7528^T) using Burrows-Wheeler Aligner version 0.7.12 (Li and Durbin, 2009). Average read depths were calculated with SAMtools version 1.2 (Li et al., 2009) and SNPs were identified using FreeBayes version 0.9.21 (Garrison and Marth, 2012) with a minimum depth of 5x and a minimum variant allele proportion of 0.9. Snippy was then used to pool all identified SNP positions (excluding indels) called in at least one isolate, and a multiple sequence alignment of core SNPs was generated. Recombination analysis of the strains based on the core genome was carried out using Gubbins version 2.2.1 (Croucher et al., 2015) with default parameters. Pairwise SNP differences between groups of strains based on geographical origin were calculated and summarized for strain groups using pairwise_snp_differences⁶. The BCFtools⁷ was used to extract strain-specific and group-specific SNPs from outbreak strains that were recovered from the same animal or different animals on the same farm. Protein coding genes that harbor any of the identified SNPs were extracted and functional annotation was carried out using OmicsBox version 1.4.11 (Götz et al., 2008; BioBam Bioinformatics., 2019).

Pangenome SNP analysis independent of multiple alignments and a reference genome was carried out using kSNP3.0 (Gardner et al., 2015). kSNP3 identifies SNPs based on unique stretches of nucleotides present in all genomes having SNPs in their middle position. The optimal size of the nucleotide regions flanking the SNPs (kmer) was identified using the program Kchooser available with the package. The SNP matrix file (SNPs_all_matrix.fasta) generated by the tool was used to produce a parsimony tree, which is a consensus of up to 100 equally parsimonious trees. The phylogenetic tree was visualized using Interactive Tree Of Life (iTOL) version 4 (Letunic and Bork, 2019).

A core genome-based multilocus sequence typing (cgMLST) analysis was performed for C. chauvoei using Ridom SeqSphere version 7.1.0 (Junemann et al., 2013). The tool “cgMLST Target Definer” extracted genes from the reference genome (DSM 7528^T; NZ_CP018624.1) and compared these genes through BLAST (Altschul et al., 1990) against multiple query genome sequences with gene identity and query coverage set to 90 and 100%, respectively. The 63 assembled genomes in this study were used as query genomes. BLAST (blastn) options included match reward 1, word size 11, gap open costs 5, mismatch penalty −1, gap extension cost 2. The analysis identified cgMLST targets and accessory targets. SeqSphere software assigned alleles for each gene and generated an allelic profile for all strains. Allelic profile data were used to generate minimum spanning trees (MST) using the parameter ‘‘pairwise ignore missing values’’ during distance calculation. The developed C. chauvoei cgMLST scheme is publicly available at PubMLST.org website⁸ (Jolley et al., 2018).

Variations in primary virulence factors were studied based on sequence searches using ABRicate version 1.0.1⁹ with custom sequence databases followed by alignment of individual genes from all strains. Antivirulence loci (AVL) are genetic elements present in genomes of an ancestor but absent or inactive in the pathogen since AVL expression is detrimental to the expression of some virulence phenotype. Besides, variation concerning primary virulence factors was studied based on protein alignments of individual genes created using all strains followed by visualization in Geneious Prime 2019 (Kearse et al., 2012). The presence or absence of antivirulence loci in C. chauvoei genome sequences was determined based on ABRicate [32] searches using identified homolog AVL loci (nadA, nadB, and nadC; 3.1Kb) from Clostridium septicum DSM 7534^T strain (NZ_CP023671). Nicotinate and nicotinamide metabolism pathways for C. chauvoei and C. septicum were accessed through KEGG database (Kanehisa et al., 2017).

Taxonomic classification of the sequence reads revealed 98% read homology for all strains with the NCBI Clostridium databases except one strain (strain 16S0574) whereas only 92% of the reads have matched. Further analysis of this strain showed read contamination with exogenous genetic material; hence we first mapped the sequence data of this strain to reference genomes of C. chauvoei previously described (Thomas et al., 2017), and only the mapped reads were filtered out and used for subsequent analysis. The GC content averaged 27.9% for all strains. The de novo genome assembly and annotation of the strains reflected similar genome sizes and number of encoded proteins respectively. The genome size of the assembled genomes varied between 2.76 and 2.79 Mb for all strains, with N50 values between 63 and 90Kb and a total contig number between 63 and 81 contigs per strain (Supplementary Table 2). Genome annotation revealed 2560 to 2581 protein coding genes among the strains (Supplementary Table 2).

The genomes harbored two CRISPR repeat sites (repeat region 1 and 2) separated by a protein coding gene except for one strain of unknown origin (S0132-09) where the CRISPR repeat region was continuous. The strains showed variations in the composition of CRISPR spacer sequences among the isolates. Variations were observable at the CRISPR site region which was close to the Cas proteins. The total number of CRISPR spacers among all the strains was 2,608 (1,988 for repeat region 1 and 620 for repeat region 2, respectively) whereas the number of unique spacers identified was 45 (35 for repeat region 1 and 10 for repeat region 2, respectively). The maximum and minimum number of spacers ranged from 35 to 17 and 10 to 8 for repeat regions 1 and 2, respectively. Most of the strains from Europe and the type strain had 43 spacers (33 and 10 spacers at repeat regions 1 and 2, respectively). The number of spacers identified for each strain at repeat regions, positions of spacers in form of an array is detailed in Supplementary Table 3. Two strains (12S0470; spacer 44 and S0132-09; spacer 45) each had an additional unique spacer that was not present in the type strain (Supplementary Table 3). Recent studies have pointed out the applicability of CRISPR spacers for differentiation of strains of C. chauvoei (Rychener et al., 2017). CRISPR spacers for all 64 strains were visualized as an array and were representing 18 unique patterns. The CRISPR spacers, array pattern, strain name, and region wise and country information are represented in Supplementary Figure 1. The array for spacers showed marked variations with respect to strains from Austria, Switzerland, 15S0023 (South Africa), strains from North Rhine-Westphalia, Bavaria and few strains from Lower Saxony and of unknown origin. Prophage Hunter tool predicted prophages categorized as active/ambiguous and with a prediction score above 0.5 were considered (Song et al., 2019). Four prophages and nine genomic islands were predicted for C. chauvoei strain DSM 7528^T, respectively. The sizes and protein coding genes encoded by respective prophages and genomic islands are shown in Supplementary Table 4. BRIG blast atlas showing sequence variability among strains and with respect to the prophages and genomic islands are depicted in Supplementary Figure 2.

Panaroo predicted 2,494 core genes strictly present in the 64 investigated strains comprising more than 95% of the pangenome of the strains (total 2,624 pangenes), while the accessory genome comprised less than 5% (148 out of 2,624 genes). The curve for the pangenome represents the least-squares fit for the function y = Ax^B + C with the best fit obtained with a correlation r² = 0.986 for A = 140.68, B = 0.08, C = 2425.06 (Figure 2A). The extrapolated C. chauvoei pangenome size expected after the number of genomes increases to 100 was 2,628. The fitting parameter B value was 0.08 i.e., between 0 and 1. The number of core genes after the addition of each new genome was plotted as a function of the number of genomes to create a core genome plot. The curve for the core genome represents the least-squares fit for the function y = Ae^Bx + C with the best fit obtained with a correlation r² = 0.997 for A = 90.4, B = −0.02, C = 2,475.36 (Figure 2A). The extrapolated C. chauvoei core genome size for 100 genomes was estimated as 2,487 genes. The curve for new genes represents the least-squares fit for the function y = Ax^B with the best fit obtained with a correlation r² = 0.835 for A = 39.854, B = −1.37 (Figure 2B). The extrapolated number of new genes expected after the number of genomes increases to 100 was zero indicating a closed pangenome structure. Accessory genes were majorly comprised of genes related to prophages, mobile genetic elements and insertional sequences. The variations attributed by the strains even involving different geographical locations were limited as there was no clear pattern of accessory genes sharing among the strains.

A total of 1,677 core genome SNPs was identified among all strains based on the sequence read mapping to the DSM 7528^T reference genome by Snippy analysis. Core SNP density showed an almost even distribution across genomes and ranged from 31 to 81 SNPs within every 100Kb genomic regions. Pairwise SNP variations ranged from 0 to a maximum of 331 SNPs among any pair of isolates. Pairwise SNP differences among strains are depicted in Supplementary Table 5. Pairwise SNP difference analysis within groups based on geographical origin showed a high median pairwise SNP difference value of 89 SNPs (10 strains) for Austria as compared to 41 SNPs (38 strains) from Germany (Supplementary Figure 3A). Similarly, the strains from the mountain pastures (Bavaria, Austria and Switzerland) showed a higher pairwise SNP difference value of 77 (17 strains) as compared to 38 (34 strains) from plain regions (Lower Saxony, North Rhine-Westphalia, Mecklenburg-Western Pomerania, and Schleswig-Holstein) (Supplementary Figure 3B). This is indicative of higher strain diversity in seasonal mountain pastures as compared to other regions.

The dataset investigated herein included strains obtained from six outbreaks in Germany. Each of these outbreaks was represented by two to three C. chauvoei strains. Though phylogenetic analysis revealed that the outbreak strains were closely related, the pairwise SNP distances (Supplementary Table 5) between these strains revealed a genetic distance between 11 to 112 SNPs. These strains were however isolated either from the same animal (one outbreak) or from two different animals on the same farm in the event of the outbreak (five outbreaks). To understand the genetic variations associated with these outbreak strains, we parsed the Snippy output files using BCFtools and extracted SNP variations among the outbreak strains. To define SNPs exclusively associated with outbreak strains, closely related strains that were not isolated during the outbreak were not included. For example, we excluded strain S0121-09 that was isolated in 2009 but was phylogenetically grouped with the 2010 outbreak strains (S0121-10 and S0122-10) from North Rhine-Westphalia. Similarly, we excluded strain S0260-09 isolated in 2009 but grouped with the 2012 outbreak strains (12S0468 and 12S0467) in Lower Saxony. The number of SNPs identified from the same host and outbreak strains is shown in Table 1. In total, 240 strain-specific SNPs were identified among the same animal and outbreak strains, 194 SNPs belonged to protein-coding genes and were predominately represented by non-synonymous (NS) substitutions indicating diversifying selection occurring inside the host. The number of SNP variations uniquely represented by strains recovered in the same animal and from five outbreaks is shown in Supplementary Table 6. SNPs in different positions for the same protein-coding genes observed for the outbreak strains were coding for HTH-type transcriptional regulator CymR and Nitrate reductase (three SNPs) and Sensor histidine kinase RcsC, Putrescine transporter PotE, putative M18 family aminopeptidase 1, Magnesium transporter MgtE, Putative competence-damage inducible protein, Flagellar biosynthesis protein FlhA, Leucine–tRNA ligase, UDP-3-O-(3-hydroxymyristoyl) glucosamine N-acyltransferase and for three proteins coding genes predicted as hypothetical proteins (two SNPs). Based on Blast2GO, these genes were involved in various biological process and predominantly for cellular macromolecule biosynthetic process (GO:0009059, GO:0044260, GO:0044249) (Supplementary Figure 4). Based on KEGG pathway analysis, 21 CDS were involved in purine metabolism followed by thiamine metabolism (13 CDS), amino sugar and nucleotide sugar metabolism (7 CDS).

Recombination analysis of the strains under study showed limited recombination events based on Gubbins predictions. Eight recombinant regions were predicted involving a minimum number of 4 to a maximum of 14 SNPs. Among eight, three regions involved protein coding genes representing capsular polysaccharide biosynthesis protein (10 SNPs), flagellin C (14 SNPs) and metal ABC transporter permease (6 SNPs).

Pangenome phylogenetic analysis was carried out using kSNP3. The optimal kmer value was 21. A total of 2,389 SNPs were identified in the 64 C. chauvoei genomes, of which core SNPs were 2,125 and non-core SNPs were 264. Among the SNPs occurring in protein-coding genes, 1,451 corresponds to non-synonymous (NS) substitutions and 178 were synonymous (S) type resulting in a high NS/S ratio of 8.151 indicating diversifying selection. A total of 45 SNPs were homoplastic which corresponds to SNPs shared by groups of genomes that are not clustered in the consensus maximum parsimony tree used for clustering the strains. The phylogenetic tree indicating the clusters with branch lengths expressed in terms of changes per number of SNPs and node labels showing shared SNPs (>5) among strains are shown in Figure 3.

The parsimony tree based on the pangenome SNPs indicated clusters defined based on the presence of SNPs shared. Three strains represented type strain ATCC 100092 (S0136-09, BC97-06 and DSM 7528^T) received from various sources showed relatedness with respect to shared SNPs. As mentioned earlier, strains recovered from the same farm or animal/outbreak showed high genetic relatedness. Also strains originating from the same farm or region but in longitudinal time gaps of one (S0121-09, S0021-10 and S0022-10), two (12S0467, 12S0468 and S0260-09), and 4 years (BS79-01, BS80-01 and 11S0318) showed striking relatedness. This indicates persistent lineages can be involved in repeated outbreaks in the same farm (Figure 3). Exceptions were also observed when one strain from Lower Saxony (12S0464) isolated a year before (2010) was not forming a cluster with three strains isolated from the same farm a year later (11S0315, 11S0316) and one isolated from the same animal (12S0471). Similarly, strains involved in the same cluster were originating from two farms from the Lower Saxony region (BS90-02 and BS106-02); indicating diverse genotypes can be prevalent on farm/region levels.

Clusters were mostly based on the strains’ geographical origins. The SNPs shared by a cluster had a maximum number of 76 (three strains from North Rhine-Westphalia). Shared SNPs identified in the inner nodes of the clusters ranged from 0 to 42. One exceptional cluster with divergent strains was observed where the type strain (DSM 7528^T) which was isolated before 1950 was found to be genetically related to the three strains from Austria (S0013-08, S0014-08 and S0014-08) isolated in 2008. The same cluster also involved a Canadian strain (16S0578) and was showing more relatedness to a strain of unknown origin (S0134-09). Strains from Germany were also sharing relatedness with strains from Austria. These strains were from Bavaria which shares geographical closeness to Austria. Similarly, few strains from Austria clustered with strains from Switzerland and indicate geographical proximity. Other clusters mostly involved strains from different regions of Germany viz., Lower Saxony, Lower Saxony and Mecklenburg-Western Pomerania, North Rhine-Westphalia, Schleswig-Holstein and Baden-Württemberg. Strain-specific SNPs ranged from 1 to 207. The strains harboring specific SNPs above 100 were S0132-08 (unknown origin), 1520023 (South Africa) and DRS0104052 (unknown origin) (Figure 3).

The genome sequence analysis of strain JF4335 and previous studies have characterized primary virulence factors of C. chauvoei (Frey and Falquet, 2015). The identified virulence factors included C. chauvoei toxin A (CctA), hemolysins, sialidases, hyaluronidases and internalin A proteins (Frey and Falquet, 2015). Most of the known primary virulence factors for C. chauvoei showed limited variations at the amino acid level. The CctA belongs to the leucocidin superfamily of bacterial toxins and is considered the major virulence factor of the pathogen (Frey et al., 2012). The CctA showed high conservation at the amino acid level among all strains under study. Sialidases and hyaluronidases are enzymes hydrolyzing the glycoside linkage between sialic acid molecules and break hyaluronate (a carbohydrate polymer that is part of the extracellular matrix), respectively. Both virulence factors are thought to be involved in the rapid spread of C. chauvoei through tissues (Frey and Falquet, 2015; Abreu et al., 2017). NanA sialidase showed one amino acid substitution for strain 15S0023 (South African strain) and another one for two strains originating from Switzerland (BS169-00, BS171-00). Hyaluronidase (nagH) showed multiple deletions of amino acids for strain S0131-09 and three single substitutions were observed for three strains. Except for CctA, the role of other hemolysins in blackleg pathogenesis is unclear. There was also no variation observed for the hemolysin III and hemolysin XhlA. The amino acid substitutions and deletions identified in the virulence factors and respective positions are depicted in Table 2.

Antivirulence loci (AVL) are described as genes whose expression is incompatible with the pathogenic lifestyle. Earlier studies for Shigella have described two genes nadA and nadB of Shigella as antivirulence genes (Prunier et al., 2007). There was a notable absence of antivirulence loci encoding nadA (quinolinate synthase [EC:2.5.1.72]), nadB (L-aspartate oxidase [EC:1.4.3.16]), and nadC (nicotinate-nucleotide pyrophosphorylase (carboxylating) [EC:2.4.2.19]) genes in C. chauvoei in comparison to C. septicum. Both genes coding enzymes are required for the conversion of L-aspartate to quinolinate, a precursor to NAD synthesis. In contrast, all the four strains of C. septicum genomes (strains DSM 7534 (NZ_CP023671/NZ_CP023672.1), VAT12 (NZ_CP034358), P1044 (NZ_CABMIZ000000000.1) and MGYG-HGUT-02373 (NZ_FLTT00000000.1) were predicted for an intact coding sequence for nadA (100% consensus at nucleotide and protein levels) whereas both nadB and nadC coding sequences were disrupted by stop codons.

The C. chauvoei cgMLST scheme consisted of 2,223 core genes for assessing the genetic diversity among isolates. The minimum spanning tree for the strains highlighting the origin information of isolates is shown in Figure 4. Similar to kSNP3 clustering, the three isolates representing the type strain (ATCC 100092) formed a group with three strains form Austria and another strain of unknown origin in MLST tree as well. This may further indicate the possibility of the type strain being initially isolated from Austria (Tyrol region) as they also represent similar CRISPR spacer patterns and share similar genetic variations for one virulence factor (Table 2). The limited genetic variations observed among the strains may be due to laboratory passages/experiments incurred SNPs and/or biases from different sequencing technologies applied in the study. The tree shows that the SNP distance for strains outside Europe is greater when compared to strains within Germany. Within Germany, the strains from North Rhine-Westphalia are the most divergent strains. Interestingly, the core genome MLST phylogeny showed collinearity with most of the outliner CRISPR arrays.