Four different sample types which included meat processing equipment (n = 90), bovine fecal samples (n = 3), and meat juice samples of unspoiled meat (n = 153) or BPS meat (n = 3) were screened for the presence of CEC by quantitative real-time PCR (qPCR) as previously described (Wambui et al., 2020c). Briefly, sponge samples obtained from the meat processing equipment were homogenized with Reinforced Clostridia Media (RCM; 25 ml) for 60 s using a stomacher. Hungate tubes containing 10 ml of the homogenized meat processing equipment sample were heated at 80°C for 10 min. The meat processing equipment and meat juice samples of unspoiled meat were incubated anaerobically for 10–14 days at 8°C prior to screening while meat juice samples from BPS samples were screened immediately after sampling. Current and subsequent anaerobic incubations were carried out in rectangular anerobic boxes (7.0 L; bioMérieux, Inc., Marcy l’Etoile, France), and the anaerobic conditions were generated by three 2.5 L AnaeroGen Sachets (Thermo Fisher Scientific, Waltham, MA, United States) per box. DNA was extracted from 100 μl of each sample using the MagNa Pure LC DNA Isolation Kit III (Roche, Rotkreuz, Switzerland) by the MagNa Pure LC robotic workstation (Roche). Each sample was mixed with 10 μl lysozyme (20 mg/ml in phosphate-buffered saline) solution and incubated at 37°C for 30 min. Subsequent steps from lysis with proteinase K to DNA elution were carried out according to the instruction manual of the DNA extraction kit. The qPCR primers, probes and protocol used herein were as previously described (Brightwell and Clemens, 2012). qPCR-positive meat processing equipment and meat juice samples of unspoiled meat were enriched anaerobically 4 weeks at 8°C. The bovine fecal samples were obtained from a previous study (Wambui et al., 2021a) where attempts to isolate CEC strains were unsuccessful and thus had been stored anaerobically at 4°C.

Isolation of CEC from meat juice samples of BPS samples was carried out by direct plating. Briefly, the samples were serially diluted 10-fold in saline solution and plated on Columbia agar supplemented with 5% defibrinated sheep blood (CBA). Isolation of CEC from the enriched meat processing equipment, bovine fecal samples, and meat juice samples of unspoiled meat was carried out after elimination of competitive microflora and activation of CEC spores as previously described without modification (Wambui et al., 2021b). Colonies on CBA displaying clostridial characteristics described previously (Húngaro et al., 2016) were selected and purified anaerobically on CBA for 3 weeks at 8°C. The isolates were confirmed as members of CEC by qPCR as described above.

The extraction and whole-genome sequencing of CEC isolates was carried out as previously described (Wambui et al., 2020a). Briefly, the genomic DNA was isolated using the DNA blood and tissue kit (Qiagen, Hombrechtikon, Switzerland). The sequencing outputs, which were 150–300 bp pair-ended reads, were prepared using Nextera DNA Flex chemistry using the Nextera DNA Flex Sample Preparation Kit (Ill) as per manufacturer’s guidelines (Illumina, San Diego, CA, United States). The resulting transposon-based libraries were sequenced on a MiniSeq sequencer (Illumina) with a minimal coverage of 50-fold. The MiniSeq MidOutput Reagent Cartridge (300 cycles) was used. Demultiplexing and adapter trimming was done using the Miniseq local run manager version 2.4.1 using standard settings. The reads were checked for quality using FastQC (Andrews, 2010) then assembled with SPAdes v. 3.12.0 (Bankevich et al., 2012) using Shovill 1.0.9.¹ Since Clostridium algoriphilum DSM 16153^T has not been sequenced yet, we purchased the strain from the Leibniz Institute DSMZ-German Collection of Microorganisms and Cell Cultures, revived it as previously described, and sequenced it as described above. The quality of the genomes was checked using ContEST16S (Lee et al., 2017) and CheckM (Parks et al., 2015). The complete rpoB gene sequences were identified using blastn and extracted from the genomes with samtools using the faidx option and the blast results as input (Altschul et al., 1990; Li et al., 2009). These and the rpoB sequences of 34 publicly available CEC genomes (Wambui et al., 2021a), were aligned in CLC Workbench Genomics v. 8.1 (Qiagen, Aarhus, Denmark) using default settings while the phylogenetic tree was created with the same software using the Maximum likelihood Phylogeny method whereby the neighborhood joining method using the Jukes Cantor model were applied, respectively. Bootstraps were based on 1,000 replicates. The correct species assignment of the CEC isolates was carried out in silico from the pairwise nucleotide comparison of Average Nucleotide Identity (ANI), which was determined with pyANI (Pritchard et al., 2016) using the BLAST algorithm. The rpoB-based phylogeny consisting of the Type/Representative strains of CEC, Clostridium algidicarnis, Clostridium frigidicarnis, Clostridium gasigenes, Clostridium putrefaciens, and Clostridium argentinense was used to validate a divergent isolate as a member of CEC. Where necessary, resulting newick files were visualized in iTOL (Letunic and Bork, 2019). Core genome MLST was performed in Seqsphere+8.2.0 (Ridom, Münster, Germany). The C._estertheticum_cgMLST scheme was constructed using the type strain DSM 14864^T as seed and all other available C._estertheticum_genomes as testers. Standard settings were applied for target selection.

The 16 and 34 CEC genomes from the current and previous studies, respectively (Supplementary Tables 1A,B) were annotated using Prokka (Seemann, 2014) and RAST (Brettin et al., 2015). The GFF3 output files from Prokka were used for pan-genome analysis of the 50 CEC genomes in Roary v3.11.2 (Page et al., 2015) with the minimum percentage identity for BLASTP set at 60%, 70%, 80%, 90%, and 95%. For each setting, four different gene classes grouped into “core” (99% ≤ strains ≤100%), “soft-core” (95% ≤ strains ≤100%), “shell” (15% ≤ strains <95%), and “cloud” (0% ≤ strains <15%) were obtained with the -s setting for paralogues clustering. Pan-genome analysis were further carried out for the 25 C. estertheticum genomes in Roary with the minimum percentage identity for BLASTP set at 80%. Recombination events in CEC were estimated using Gubbins version 2.4.1. (settings: -f 30) using coreSNP alignment and tree as input (Croucher et al., 2015). For this purpose, the core genome alignment was created using Parsnp v1.5.3 (Treangen et al., 2014). The maximum likelihood phylogeny of C. estertheticum was determined using RAxML (Randomized Axelerated Maximum Likelihood; Stamatakis, 2014). The phylogenetic trees were visualized and annotated with metadata in iTOL (Letunic and Bork, 2019). Pan-genome characteristics were calculated using the package micropan in R using 500 the random permutations of genome ordering (Snipen and Liland, 2015).

Function annotation and classification of proteins were performed by sequence comparison against the eggNOG database 5.0 (Huerta-Cepas et al., 2019), using the eggNOG-Mapper v2 (Cantalapiedra et al., 2021). Putative enzymes involved in the breakdown, biosynthesis, or modification of carbohydrates were identified using CAZy database using dbCAN using default settings (Drula et al., 2022). The genes involved in the acetone–butyrate–ethanol fermentation were identified from the previous work on CEC (Palevich et al., 2021) and other solventogenic clostridia (Li et al., 2020). In addition, the genes annotated as hydrogenases were manually identified in the RAST annotated CEC genomes. In total, a comprehensive list of 48 genes (Supplementary Table 2) was established and each gene was manually identified and counted in each of the RAST annotated CEC genome. Where applicable, the genes and pathways were validated using the KEGG (Kyoto Encyclopedia of Genes and Genomes) database (Kanehisa et al., 2016).

For phenotypic characterization, respective cultures were initially grown anaerobically at 8°C on CBA then standardized to McFarrand 3.0. Utilization of fermentable substrates by C. estertheticum and C. tagluense was determined using the standardized cultures as previously described (Wambui et al., 2020a) using the API20A kit (bioMeriéux, Marcy l’Etoile, France) and following manufacturer’s specifications. The API20A strips were incubated at 8°C anaerobically for 2 weeks. For motility test and gas production in meat, standardized cultures were first inoculated 1:100 into RCM broth and incubated at 8°C for 10 days. The motility test was performed in liquid media which essentially determines swimming motility, requires shorter incubation time than growth on agar, and is easily reproducible. For this test, 200 μl of each culture was slowly injected into Hungate tubes containing fresh RCM broth using a syringe and the tubes incubated for 2 weeks at 8°C without shaking. Motility was determined as the movement of the cultures from the surface toward the bottom of the tubes. Gas production in meat was carried out as previously described with slight modifications (Dorn-In et al., 2022). Five hundred microliter of each culture was spread on the surface of approximately 500 g of beef sample. The samples were vacuum-packed and stored for 28 or 56 days at 8°C. Uninoculated vacuum-packed beef samples were included as negative controls. At the end of the storage period, meat spoilage was scored from 0 to 5 according to a previous description (Boerema et al., 2007). Loss of flagellar genes has been shown to increase exopolysaccharide synthesis (Watnick et al., 2001). Therefore, CEC colony morphology that could be related loss of flagellar genes, including sliminess, were determined on cultures grown on CBA for 8 weeks at 8°C. Each experiment was carried out in three biological replicates apart from the gas production in meat with two biological replicates.

Data for substrate utilization, motility test, and gas production were described qualitatively. The distribution of carbohydrate-active enzymes was described as means and SD and the means compared using the Mann–Whitney U (p = 0.05) with Bonferroni correction. Data analysis and visualization were carried in R Studio Version 1.1.463 (RStudio, Inc., Boston, United States).

In an effort to expand the existing collection of sequenced CEC strains, 15 new strains, which were confirmed by qPCR, were isolated from meat processing equipment (n = 2), bovine fecal samples (n = 3), and meat juice samples of unspoiled meat (n = 7) and BPS meat (n = 3). Based on the rpoB gene phylogeny (Figure 1A), 14 strains were identified as C. tagluense (n = 8) and C. estertheticum (n = 6). Their correct species identification was confirmed and validated by ANIb analysis (Supplementary Figure 1). The ANIb between the C. tagluense strains and C. tagluense DSM 17763^T ranged between 97.27% and 97.32% while the ANIb values between the C. estertheticum strains and C. estertheticum DSM 8809^T ranged from 95.80% to 98.18% confirming the correct species identification of the strains. Strain, CS001, formed an outgroup compared to other known 49 CEC strains (Figure 1A) suggesting a novel species within the CEC. Interspecies ANIb between CS001 and the other CEC strains ranged from 78.60% to 81.36% for C. psychrophilum DSM 14207^T and C. tagluense DSM 17763^T, respectively. Furthermore, strain CS001 had an ANIb of <81.7% with any of the 122 representative genomes of the genus Clostridium in the NBCI database, confirming it is a novel CEC species. The new species is herein referred to as genomospecies4. rpoB-based phylogenetic reconstruction with 16 type/representative strains from 11 CEC species, four non-CEC Clostridium species that are commonly isolated from meat (C. algidicarnis, C. frigidicarnis, C. gasigenes, and C. putrefaciens) and C. argentinense, which is the known closest Clostridium species to CEC, further provided evidence that genomospecies4 is a novel species within CEC despite its high divergence (Figure 1B). The novel species shows CEC is both expansive and diverse.

The strains, CM032, CM033, and CM034, isolated from BPS cases (Figure 1C; Supplementary Figure 2) were identified as C. estertheticum, and clustered together (Figure 1A). cgMLST of species C. estertheticum demonstrated the three strains had SNP differences of two to three, indicating that they are clonal (Supplementary Figure 3). Interestingly, strains CM033 and CM034 were isolated from BPS involving horse meat (Figure 1C), while their involvement in BPS was validated with CM034, which caused BPS in vacuum-packed beef (Figure 1D; Table 1). To the best of our knowledge, this is the first reported case of BPS in horse meat caused by C. estertheticum. Phenotypic analysis revealed CS001 from the novel species caused meat spoilage with a score of 3 which corresponds with blown puffy packs (Table 1). We further compared the new species with the other recently described species, genomospecies2 and genomspecies3 whose representative strains, CM028 and CF012, caused spoilage with a score of 3 and 2, respectively with the latter score corresponding to loss of vacuum without evidence of substantial gas production. This further demonstrates the significant involvement in meat spoilage of not only the known CEC species, but also other newly described species.

Clostridium algoriphilum DSM 16153^T was also sequenced within this study. On the rpoB phylogeny (Figure 1A), the strain was closely related to C. psychrophilum DSM 14207^T. This was further confirmed by the ANIb of 84.63% between the two genomes. Clostridium algoriphilum DSM 16153^T was distantly related to the newly identified genomospecies4 CS001 (ANIb = 79.53%).

Although CEC forms a single clade within the Clostridium genus, their genomes demonstrated high variability with respect to size, number of coding sequences (CDS), and G + C content (Supplementary Table 1). The size varied from 3.662 Mbp in genomospecies2 CM028 to 5.555 Mbp in genomospecies1 FP3 while the CDS ranged from 3.337 in genomosspecies2 CM028 to 5.271 in C. tagluense FP2. The G + C content ranged from 29.9% in C. psychrophilum DSM 14207 to 31.7% in genomospecies2 CM027 and genomospecies3 CF012. The size, CDS and G + C of CS001 from the newly identified species, genomospecies4, were 4.139 Mbp, 3,831 and 31.6%, respectively, and 4.520 Mbp, 4,524, and 30.7%, respectively, in the newly sequenced C. algoriphilum DSM 16153^T. Overall, genomospecies2 has the smallest genome size within CEC (Supplementary Figure 4).

The variability of CEC genomes may have contributed to interspecies or intraspecies diversity, population genetics, and evolution. We therefore explored the entire genomic repertoire of the CEC by systematically determining its pan-genome with the BLAST identity cutoff at 60%, 70%, 80%, 90%, and 95%. The cutoff had a significant effect on the CEC pan-genome whereby its increase corresponded with a decrease in core- and soft-core genomes and an increase in shell- and cloud-genomes (Supplementary Table 3). Because of these effects, further pan-genome analyses were performed on the data set from 80% BLAST identity cutoff whereby the core- and soft-core genomes comprised of 748 and 942 genes, respectively, while the shell- and cloud-genomes had 7,296 and 28,741 genes, respectively (Figure 2A). The core genome represented 14.19%–22.40% of the gene content of each strain, while the cloud genes represented 77.72% of the pan-genome. The large accessory genome highly suggests a significant degree of genomic diversity within CEC. The pan-genome frequency showed a proportional relationship between the number of genes and genomes and the curve analysis indicates the CEC pan-genome is open (Figure 2B). The pan-genome of CEC had an alpha of 0.456. Similarly, the pan-genome of C. estertheticum had an alpha of 0.392. Given that the pan-genome is closed if the estimated alpha is above 1.0, both pan-genomes can be considered to be open.

In a next step, we evaluated the cluster of orthologous groups (COGs) functional classification to define possible differences in the functions encoded by the core, soft-core, cloud-, and shell-genomes of CEC (Figures 2C–E). In all four cases, the highest proportion of genes were COG category S (Function unknown). Among the known functions, both core and soft-core gnomes consisted mostly of COG categories J (Translation), L (Replication and repair), K (Transcription), E (Amino Acid metabolism and transport), and F (Nucleotide metabolism and transport), which suggest that the processes controlling the flow of genetic information in CEC are conserved. On the other hand, the shell-genome was biased toward K (Transcription), T (Signal Transduction), E (Amino Acid metabolism and transport), C (Energy production and conversion), G (Carbohydrate metabolism and transport) while the cloud genome consisted mostly of L (Replication and repair), K (Transcription), M (Cell wall/membrane/envelop biogenesis), G (Carbohydrate metabolism and transport) and T (Signal Transduction) categories. The identification of carbohydrate metabolism and transport category in both cloud- and shell-genomes suggests that these processes are among the most genetically variable within CEC.

To determine whether recombination events have influenced the population structure of CEC, we first determined the occurrence of recombination within the core genome of CEC. The core-genome-based phylogenetic tree (Figure 3A) confirmed all species identified by the rpoB gene (Figure 1A), but the population structure varied slightly whereby it was divided into Clade 1 consisting of C. estertheticum, C. frigoris, and C. psychrophilum, Clade 2 consisting of C. bowmanii C. algoriphilum C. lacusfryxellense genomospecies2 and genomospecies4 and Clade 3 consisting of C. tagluense, genomospecies1, and genomospecies3. Within the core genome, we identified 252 recombination events, out of which 70 were ancestral, that is, they occurred between or among strains, while 182 were strain-specific (Supplementary Figure 5). Among, the ancestral recombination events, the majority occurred within the species and more specifically between strains that were more closely related than distantly related. The reconstruction of the core phylogenetic tree without the recombined regions (Figure 3B) revealed the influence of recombination on the population structure of CEC. While three clades were identified, the species within these clades varied slightly with those before the removal of recombination events. Clade 1 consisted of C. estertheticum and C. frigoris, Clade 2 consisted of C. algoriphilum, C. lacusfryxellense, and C. psychrophilum and genomspecies2 while Clade 3 consisted of C. tagluense, C. bowmanii, genomospecies1, genomospecies3, and genomospecies4. The data show that recombination events within CEC occur regularly and have also contributed to the current population structure of the complex.

The variable interspecies genome sizes within CEC shows genes have been acquired or lost within the complex, which might contribute to speciation. This prompted us to determine signature genes that have been acquired or lost within CEC and determine whether the loss or gain may have contributed to speciation events within the complex. In this case, we investigated genomospecies2, which has the least genome size of less than 4.0 Mbp (Supplementary Figure 4). An analysis of the CEC accessory genome led to the identification of nine genes that were exclusively absent in genomospecies2. These included chemotaxis-related genes (cheA, cheD, and cheW), RNA polymerase sigma-D factor, and two putative type III secretion system genes that were annotated as putative ATP synthase YscN and Yop proteins translocation protein U. Further analysis of the genes showed each CEC strain had 23–28 genes related to chemotaxis and flagella assembly while genomospecies2 strains had one gene each (data not shown). A whole-genome sequence-based comparison of genomospecies2 and the two closely related species, C. bowmanii and C. lacusfryxellense using RAST annotation server showed up to 55 related genes were absent in genomospecies2 (Figure 4A) further confirming the absence of these genes. The lack of these genes was consistent with reduced motility of the genomospecies2 CF011 when compared with C. bowmanii DSM 14206 in RCM media (Figure 4B). An investigation of genomospecies2 strains CF011, CM027, and CM028 colony morphology showed they form significantly larger and slimy colonies after extended incubation periods (>8 weeks) compared with all other CEC species (Figure 4C). Further analysis of gene presence/absence in the pan-genome matrix did not reveal any exopolysaccharide genes that were specific to genomospecies2. Therefore, genomospecies2 is currently the only known CEC species lacking genes related to flagella signaling and assembly and this loss influences the motility and may also be linked with the unique colony morphology.

Although CEC causes meat spoilage, the spoilage is distinguished as either corresponding with formation of copious or low amounts of gas (Table 1). This prompted us to investigate the presence/absence of conserved metabolic genes that can explain these differences using C. estertheticum and C. tagluense as representative species. The analysis of the CEC accessory genome revealed 166 and 412 genes that were only present in C. estertheticum and C. tagluense, but no obvious genes could be linked to meat spoilage that is characteristic of the two species or CEC in general (data not shown). Therefore, we hypothesized that these metabolic genes are also conserved with other species given the high relatedness of closely related CEC species. Therefore, the CEC species were divided into two major phylogroups, phylogroup 1 (PG1) and Phylogroup 2 (PG2), which is consistent with the rpoB gene phylogenetic tree. PG1 consisted of C. estertheticum, C. frigoris, C. psychrophilum, and C. algoriphilum while PG2 consisted of C. tagluense, C. lacusfryxellense, C. bowmanii, genomospecies1, genomospecies2, genomospecies3, and genomospecies4. The analysis identified 63 and 27 genes that were conserved in PG1 and PG2, respectively. Elimination of seven genes with similar annotations between the two groups and 12 hypothetical proteins allowed 71 (49 and 14 genes in PG1 and PG2, respectively) genes to be compared. Interestingly, genes linked to hydrogen and carbon dioxide production could be identified among the 71 genes. Among these included formate hydrogenlyase subunits 4 and 7 and Hydrogenase-4 components B and G that were specific to PG1 while hydrogenase isoenzymes formation protein HypE was specific to PG2. This suggested that PG1 and PG2 have evolved different fermentation pathways that might influence gas production. Further comparison of the 71 genes revealed pentose metabolic genes (arabinose metabolism transcriptional repressor, arabinose operon regulatory protein, arabinose-proton symporter L-arabinose isomerase, L-ribulose-5-phosphate 4-epimerase, ribulokinase, and transaldolase) and lactose and galactose metabolic genes (lactose permease and tagatose-6-phosphate kinase) were identified in PG1 only. This supports the COG analysis and further suggests PG1 and PG2 have evolved different carbohydrate metabolic pathways which might not only influence their habitat adaptation, but also meat substrate utilization and their respective meat spoilage mechanisms.

To determine the differences in carbohydrate metabolic pathways between PG1 and PG2, we specifically investigated the distribution of six groups of carbohydrate-active enzymes (CAZy), namely, carbohydrate esterases (CE), carbohydrate-binding modules (CBM), glycoside hydrolases (GH) glycosyltransferases (GT), polysaccharide lyases (PL), and auxiliary activities (AA), in C. estertheticum and C. tagluense genomes (Supplementary Figure 6), which represented PG1 and PG2, respectively. In C. estertheticum, GHs were the most frequent CAZy genes and ranged from 40 to 69 per genome. In contrast, GTs, which ranged from 30 to 42 per genome, were the most frequent in C. tagluense. Overall, the frequencies of AAs, PLs, CBMs, GTs, and GHs were higher (p < 0.05) in C. estertheticum than in C. tagluense. However, GTs were more frequent (p < 0.05) in C. tagluense than in C. estertheticum. An in-depth analysis of the GHs revealed the two species encoded 53 different types of enzymes (Figure 5), with C. estertheticum having the highest diversity of the enzymes. Specifically, 36 GHs were only identified in C. estertheticum and eight of these GHs (GH3, GH20, GH32, GH36, GH42, GH43, GH51, and GH115) were present in all 25 C. estertheticum genomes while two GHs (GH28 and GH105) were present in 24 C. estertheticum genomes. In contrast, only three GHs were exclusively identified in C. tagluense. Among these, GH171 was present in all C. tagluense genomes while GH129 was absent in only one genome. Despite the differences, nine GHs (GH13, GH18, GH23, GH39, GH65, GH73, GH94, GH109, and GH156) were present in all C. estertheticum and C. tagluense genomes. Among these, GH13, GH94, and GH109 were significantly higher (p < 0.05) in C. estertheticum than in C. tagluense while GH73, and GH156 were significantly higher (p < 0.05) in C. tagluense than in C. estertheticum (Supplementary Figure 7). Therefore, this confirmand the variable distribution of carbohydrate utilization genes between C. estertheticum and C. tagluense.

The variable distribution of carbohydrate utilization genes between the two species suggested the range of substrates they can utilize also differs. To validate this, we determined the utilization of 16 substrates (Table 2). Clostridium estertheticum strains utilized D-glucose, D-lactose, D-sucrose, D-maltose, D-salicin, D-xylose, L-arabinose, D-cellobiose, and D-raffinose. The strains showed variable utilization of six other substrates and did not utilize D-trehalose. Comparatively, C. tagluense strains utilized only three substrates, D-glucose, D-maltose, and D-trehalose while three strains utilized glycerol. Therefore, the high number of GHs in C. estertheticum corresponded to utilization of diverse fermentable substrates compared to C. tagluense.

Carbohydrate utilization is liked with gas production in Clostridium genus through the acetone–butanol–ethanol fermentation pathway. To determine whether the differences in gas evolution during meat spoilage between the C. estertheticum and C. tagluense is solely linked with differences in the carbohydrate utilization genes or there exists other underlaying genetic influences, we first reconstructed the acetone–butyrate–ethanol fermentation pathway to include steps for the production of lactate, formate, ethanol, acetate, butyrate, acetone, acetoin, and butanol (Figure 6). The genes involved in the pathway were present in both species with minor variations including copy numbers (Figure 7). This suggested the central acetone–butyrate–ethanol fermentation pathway was conserved between the two species at the genetic level and cannot fully explain the differences in gas production. We therefore carried in-depth analysis whereby we identified and compared two different processes that produce H₂ and CO₂ through the oxidation of pyruvate in the acetone–butyrate–ethanol fermentation pathway (Figure 6). The first process involves the oxidative decarboxylation of pyruvate to acetyl-coA and CO₂ by pyruvate ferrodoxin oxidoreductase (PFOR). This process requires the reduction of ferredoxin, which is further oxidized simultaneously by [FeFe] hydrogenase, HydA, and ferredoxin NADH oxidoreductase (FNOR) resulting in the formation of H₂ (Figure 6). All genes involved in this process including pfor, fnor, and hydA were present in both species although the copy number of hydA in C. tagluense was higher than in C. estertheticum (Figure 7). The genes involved in the maturation of HydA, that is hydE, hydF, and hydG, were also identified (Figure 7). The second process involves the decarboxylation of pyruvate to acetyl-coA and formate by pyruvate formate lyase (PFL; Figure 6). The formate is then oxidized to CO₂ by a formate dehydrogenase (Fdh) or to CO₂ and H₂ by formate-hydrogenase lyase complex (FHL; comprising of Fdh and a hydrogenase; Figure 6). The presence of the pfl gene indicated the process of formate production is conserved in both species (Figure 7). However, the fdh gene was present in all C. tagluense genomes and six C. estertheticum genomes (Figure 7).

Interestingly, we identified other gene clusters, which likely encode putative hydrogenases that can play an important role in production of molecular hydrogen. A gene cluster for bifurcating [FeFe] hydrogenase encoding genes for alpha, beta, and gamma subunits was identified in both species. The bifurcating [FeFe] hydrogenase is herein referred to as Hya and is encoded by the hya cluster consisting of hyaABC genes, encoding the alpha, beta, and gamma subunits, respectively (Figure 7). The presence of hyaABC genes suggests the metabolic process of its product is conserved in both species and together with HydA hydrogenase may also not fully explain the differences observed in gas production at the genetic level. Additionally, we found three different hydrogenase gene clusters in both species (Figure 7). One cluster herein referred to as hyc cluster consists of six genes hycBCDEFG was exclusive to C. estertheticum. The second cluster, herein referred to as hyp cluster consisting of nine [NiFe] hydrogenase genes in the order hypABCDEFGHB was present in all C. tagluense strains and nine C. estertheticum strains. The strains with the hyp cluster formed two distinct clusters within C. estertheticum suggesting the intraspecies presence or absence of the cluster was evolutionary. The third [FeFe] hydrogenase gene cluster, here in referred to as hyb cluster consisted of four genes hybABCD, and was exclusive to C. tagluense. We speculated the variable distribution of the hyb, hyc, and hyp hydrogenase gene clusters could explain the phenotypic differences in gas production in meat between C. estertheticum and C. tagluense.

To link the genome and phenotype, we tested gas production by both species in meat. Three out of the four C. estertheticum strains DSM 14864, CM020, and CM034 caused pack distention while one C. estertheticum strain, CEST001, and all four C. tagluense strains CM008, CM023, CM024, and CS002, did not cause pack distention (Figures 8A,B). This was surprising since variable differences in pack distention by C. estertheticum have not been reported before. At the genomic level, the four C. estertheticum strains had hyd, hya, and hyc gene clusters while all four C. tagluense strains had hyd, hya, hyb, and hyp genes. Clostridium estertheticum CEST001, which did not cause pack distention also had the hyp genes, which led to the suggestion that the presence of the hyp gene cluster was linked with absence of pack distention. Further analysis of selected C. estertheticum strains with or without the hyp genes consistently showed C. estertheticum CM035 (hyp −) caused pack distention while C. estertheticum CF003 and CF007 (hyp +) did not cause pack distention despite extended incubation periods of 56 days (Figure 8C). Therefore, the differences in the intra- and interspecies distribution of hydrogenases gene clusters, and specifically the hyp cluster, may be associated with meat spoilage with or without gas production by C. estertheticum and C. tagluense species.

We finally took advantage of the large repertoire of C. estertheticum genomes (n = 25) to perform an intraspecies comparative genomics. Using a 80% BLAST identity cutoff, we showed C. estertheticum has a pan-genome of 20,111 genes. The core- and soft-core genomes comprised of 1,771 and 2,036 genes, respectively, while the shell- and cloud-genomes had 4,354 and 13,721 genes, respectively. The large accessory genome highly suggests a significant degree of genomic diversity within C. estertheticum.

We further created a maximum likelihood phylogenetic tree which revealed six C. estertheticum putative clades (Supplementary Figure 8). Interestingly, the C. estertheticum reference strains known to cause BPS, DSM 8809, DSM 14864, MA19, and MA41 and the BPS causing strains in the present study, CM020 and CM035 clustered together in Clade 5 (Supplementary Figure 8), which is hyp − (Figure 7). The three strains causing BPS cases in Switzerland CM032, CM033, and CM034 clustered in Clade 3, which is also hyp − (Figure 7). Because of the putative clades’ relevance to BPS, we determined the presence of other important genes that are uniquely specific to these clades. Four genes annotated as “putative response regulatory protein,” “Spore germination protein B2,” “Spore germination protein B1,” and “DegV domain-containing protein” were identified. The identification of the two spore germination proteins genes that are specific to these two putative clades suggests they have evolved different spore germination pathways.