Transcriptome, metabolome, and inflammatory and oxidative properties of Clostridium butyricum CB1002 clinical strain and its isogenic mutant Δhbd

Article type: Research Article

Keywords: interferon gamma, metabolome, transcriptome, redox, necrotizing enterocolitis, NEC

Authors: Lancelot Nicolas, Victoria Mesa, Laurent Ferraris, Johanne Delannoy, Thomas Cokelaer, Céline Chollet, Florence Castelli, Amal Zerrad-Saadi, Frédéric Barbut ⚬, Julio Aires ⚬

Affiliations: Université Paris Cité, Inserm U1139 (FPRM)555089https://ror.org/05f82e368, Paris, France; Plate-forme Technologique Biomics, Institut Pasteur, Université Paris Cité555089https://ror.org/05f82e368, Paris, France; Laboratoire d’Etude du Métabolisme des Médicaments, CEA, Université Paris Saclay27048https://ror.org/03xjwb503, Gif-sur-Yvette, France

License: Copyright © 2026 Nicolas et al. CC BY 4.0 This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license.

Article links: DOI: 10.1128/aem.01106-25  |  PubMed: 41729214  |  PMC: PMC12997760

Relevance: Moderate: mentioned 3+ times in text

Full text: PDF (1.5 MB)

INTRODUCTION

Necrotizing enterocolitis (NEC) is one of the most common neonatal intestinal diseases occurring in extremely preterm neonates (PNs) (ref. 1, ref. 2). In high-income countries, the incidence of NEC ranges from 2% to 13%, depending on the geographic area and the reporting period (ref. 3). NEC is associated with a high mortality rate up to 30% as well as long-term complications, such as short bowel syndrome, growth delays, and neurodevelopmental delays in infants who survive surgical intervention (ref. 4). NEC is a multifactorial disease, whose pathogenesis remains unclear. Data on microbial dysbiosis have highlighted the impact of abnormal microbiota of the PN on NEC (ref. 5); however, the mechanistic understanding of the relationship between bacteria and NEC development is incomplete. No specific microorganism has been definitively linked to NEC. However, the most commonly proposed bacterial etiological agents of NEC belong to the Enterobacteriaceae and the Clostridium genera (ref. 6, ref. 7). Concerning Clostridium, epidemiological studies, clinical signs, and animal models have reported the involvement of clostridia strains in the development of NEC, including Clostridium butyricum (ref. 7).

C. butyricum belongs to the sensu stricto cluster I of the genus Clostridium, a strictly anaerobic, Gram-positive, sporulating bacterium that is widely distributed in the environment and in the intestines of humans and animals (ref. 8). C. butyricum has been reported to be part of the gut microbiota in PN (ref. 9). However, intestinal colonization of PNs by C. butyricum has been associated with NEC (ref. 10ref. –ref. 12). In a preterm piglet model, C. butyricum was significantly overrepresented in colonic mucosal samples from piglets with NEC (ref. 13). Smith et al. reported a correlation between C. butyricum presence and pneumatosis intestinalis in tissue samples from neonates with NEC (ref. 14).

Pneumatosis intestinalis and portal venous gas are pathognomonic radiographic signs of NEC and have been proposed to result from gas production due to bacterial fermentation metabolism. Using germ-free quails fed with a lactose diet and mono-associated with C. butyricum, animals developed cecal NEC-like lesions (ref. 15). In contrast, animals colonized with Klebsiella pneumoniae or Clostridioides difficile did not develop NEC-like lesions (ref. 15). Furthermore, animals colonized with C. butyricum CB1002 clinical strain showed significantly more frequent and severe intestinal NEC-like lesions than those colonized with CB1002 isogenic strain deleted in the hbd gene, which encodes a β-hydroxybutyryl-CoA dehydrogenase involved in the initial steps of the central carbohydrate fermentation pathway, converting pyruvate to butyrate production (ref. 16). Of note, inflammatory and redox analysis of the host-cellular response has not been performed.

No specific phenotypic or molecular marker can distinguish asymptomatic carriage of C. butyricum isolates. Efforts to identify the genes responsible for encoding specific virulence factors in this non-toxigenic species have been unsuccessful. However, possible non-specific virulence genes have been suggested based on a comparative genomics analysis (ref. 17, ref. 18). To test the hypothesis that gene expression may contribute to virulence, we performed a transcriptomic comparison of the clinical C. butyricum 1002 and C. butyricum 1002 Δhbd strains. In addition, we investigated the pro-inflammatory and redox properties of both strains in vitro.

MATERIALS AND METHODS

Bacterial strains and culture conditions

Strain C. butyricum CB1002 was isolated from the fecal samples of a fatal NEC case (ref. 19). CB1002 Δhbd (HBD) is the isogenic strain of CB1002 WT, inactivated for butyrate production via homologous recombination of the hbd gene, which encodes a β-hydroxybutyryl-CoA dehydrogenase involved in butyrate fermentation (ref. 16). Strains were cultured on Columbia cysteine agar (Oxoid; Thermo Fisher Diagnostics SAS, Dardilly, France) supplemented with 5% sheep blood at 37°C for 72 h. Liquid cultures were grown overnight at 37°C in either 10 mL brain heart infusion (BHI) or Tryptone-glucose-yeast extract (TGY) broth. All cultures were incubated in an anaerobic chamber (Coy Laboratory Products Inc, Grass Lake, MI, USA) with a gas mixture of 90% N₂, 5% H₂, and 5% CO₂.

RNA extraction, rRNA depletion, and cDNA library preparation

For RNA extraction, CB1002 WT and CB1002 HBD cultures were grown overnight in 20 mL BHI broth at 37°C, then sub-cultured in 250 mL BHI broth to mid-exponential phase (OD600 = 1.1 ± 0.5). Cultures were harvested at 4°C and resuspended in 500 µL RNAlater (Invitrogen, Carlsbad, CA, USA), then stored at −80°C. RNA was extracted using the FastRNA Pro Blue kit (MP Biomedicals, Irvine, CA, USA) and bead-beater instrument, then purified with the ZYMO Direct-zol RNA Miniprep Plus kit (Zymo Research, Irvine, CA, USA). RNA concentration and purity were measured using the Qubit RNA HS assay kit and Qubit 2.0 fluorometer (Life Technologies, Grand Island, NY, USA). Integrity was confirmed using the Agilent 2100 Bioanalyzer (Santa Clara, CA, USA) based on an RNA integrity number (RIN) > 7. A total of 100 ng RNA was used for RNA-Seq library construction with the Illumina stranded Total RNA Ligation Ribo0+ kit, followed by ribosomal RNA depletion. Libraries were sequenced on an Illumina NextSeq 2000 (Illumina, Inc., San Diego, CA, USA) using single-read mode (50–75 base-length reads).

RNA-Seq analysis

RNA-Seq analysis was performed using the Sequana RNA-Seq pipeline (v0.92.2; https://github.com/sequana/rnaseq) with default parameters, unless stated otherwise (ref. 20). The pipeline first checked ribosomal content (<5%) and read quality with FastQC v0.11.9 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were then trimmed for adapters and low-quality bases using fastp v0.20.1 (ref. 21). The filtered reads were aligned to the CB1002 WT genome using Bowtie2 v2.4.5A D (ref. 22). Normalization corrected for technical biases, and differential analysis between HBD and WT groups was done using DESeq2 (ref. 23) with the rnadiff container from the Damona project (damona.readthedocs.io). Genes with an adjusted P-value < 0.05 were considered significant..

CB1002 genome sequencing, assembly, and annotation

DNA extraction was performed on 24-h bacterial liquid cultures using the DNA easy UltraClean microbial kit (Qiagen, Hilden, Germany), following the manufacturer’s instructions. Paired-end libraries were prepared using the TruSeq DNA PCR-free sample prep kit (Illumina, Inc., San Diego, CA, USA) and sequenced on a HiSeq 2500 system (Illumina, Inc., San Diego, CA, USA) according to the manufacturer’s instructions. Genome sequences were assembled as previously reported (ref. 24). The genome assembly was annotated automatically using the MicroScope pipeline platform v3.14.3 (ref. 25).

GO annotation and KEGG pathway analysis of differentially expressed genes

To predict metabolic pathways, differentially expressed genes (DEGs) were compared with the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Clusters of Orthologous Groups (COG) databases for annotation and classification. Pathways were annotated using the KAAS-KEGG Automatic Annotation Server (https://www.genome.jp/kegg/kaas/), and gene families were characterized by COG functional category using eggNOG-mapper v5.0 software (ref. 26).

For GO and KEGG enrichment analysis, we used the “clusterProfiler” (version 4.10.1) package of R software to conduct GO function and KEGG pathway enrichment analysis on DEGs, particularly over-representation analysis (ORA) to evaluate associated pathways.

Protein-protein interaction network

To investigate the inter-relationship of the DEGs, protein-protein interaction (PPI) network of target genes was retrieved using the online Search Tool for the Retrieval of Interacting Genes (STRING) database version 12 (ref. 27). We set the STRING database to include both experimentally verified interactions and predicted functional connections between proteins. It assigns scores to protein pairs, indicating the strength of their functional relationships. These STRING scores are derived from the gene co-occurrence, co-expression, and previous knowledge in the database. MCL clustering with inflation parameters was applied to find out the top three clusters. The minimum required interaction score was set to 0.4.

Cell culture and isolation

All medium components were purchased from Gibco (Fisher Scientific Diagnostics SAS, Dardilly, France), unless otherwise stated. Cells were incubated for 24 h at 37°C, 95% air, 5% CO₂, and 100% humidity. Human colon adenocarcinoma Caco-2 cells were cultured and maintained in Dulbecco’s Modified Eagle Medium (DMEM)-GlutaMax-I, with [+]4.5 g/L D-glucose, supplemented with 1% (v/v) non-essential amino acids, 10% (v/v) fetal bovine serum (FBS), 10 units/mL penicillin, and 10 μg/mL streptomycin, as previously reported (ref. 28). Peripheral blood mononuclear cells (PBMCs) were isolated from healthy adult donors’ blood, without pooling, using an Histopaque solution (Sigma-Aldrich, Inc., St. Louis, USA) (centrifugation 400 × g, 20 min, at room temperature). Cells were washed twice in Dulbecco’s Phosphate Buffered Saline (DPBS)−1×, centrifuged at room temperature (400 × g, 10 min), and resuspended in Roswell Park Memorial Institute (RPMI)-GlutaMax-I, supplemented with FBS (10% v/v), HEPES buffer solution (1% v/v), and penicillin/streptomycin (1% v/v).

Bacteria and bacterial components for cellular assays

Bacteria and bacterial components were prepared using a 24-h TGY bacterial sub-culture (20 mL). After centrifugation (3,500 × rpm, 20 min, at 15°C), the culture supernatants were filtered (0.2-μm-pore-size filter; Millipore) and stored at –80°C for metabolomics analysis. The bacterial pellets were washed, resuspended in DPBS (1×), and divided into aliquots as follows: live bacteria, dead bacteria after boiling (100°C for 15 min); bacterial debris and bacterial cytoplasmic contents obtained after sonication (30 s, 15 s pause, 30 s) (Vibra-Cell, Bioblock Scientific, France) and centrifugation (13,000 rpm, 10 min, at 4°C). Bacteria and bacterial components were distributed into each well of the culture plates at a concentration of 1/20 (v/v).

Cytokine assays

Twenty-four-well culture plates were coated for 1 h at 37°C under CO₂ atmosphere (95% air, 5% CO₂, and 100% humidity) with Ultra-LEAF Purified anti-human CD3 Antibody (BioLegend, San Diego, CA, USA) (2 μg/mL in RPMI). Then, PBMC cells (1 × 10⁶ cells/mL) were added to the plates and activated with Ultra-LEAF Purified anti-human CD28 Antibody (BioLegend, San Diego, CA, USA) (4 μg/mL in RPMI), and incubated for 24 h with live bacteria and bacterial components. After centrifugation, supernatants were collected and frozen at −20°C. IL-10, IL-17, IL-22, and IFNγ secretion were quantified by ELISA (R&D Systems Inc., Minneapolis, MN, USA), following the manufacturer’s instructions. The control consisted of PBMC cells that were activated using an anti-human CD3/CD28 antibody complex. The assays included three biological experiments with technical duplicates.

Redox assay

Reactive oxygen species (ROS) production was assessed using the CM-H2DCFDA probe according to the manufacturer’s instructions (Invitrogen, Carlsbad, CA, USA). Briefly, CM-H2DCFDA probe was suspended in DMSO and diluted in Hanks’ Balanced Salt Solution (HBSS) at a final concentration of 5 µM. Caco-2 cells (1 × 10⁶ cells/well/plate) were incubated with the probe for 1 h (37°C, 5% CO₂, and 100% humidity). Then, live bacteria and bacterial components were added for 30 min co-incubation in dye-free medium. The medium was then replaced with fresh dye-free medium, and measurements were performed using a microplate reader at Ex/Em = 492/518 wavelength (Enspire 2300 Multilabel Reader—PerkinElmer, Waltham, MA, USA). The control consisted of Caco-2 cells incubated with the CM-H2DCFDA probe and an exogenous inducer of oxidative stress, tert-butyl hydroperoxide. O₂^•– production levels by NADPH oxidase were measured as previously described and adapted for bacterial samples (ref. 29). ONOO⁻ production levels were measured using Peroxynitrites Assay Kit (Cell-based) ab233468 (Abcam, Cambridge, UK) according to the manufacturer’s instructions. The controls consisted of Caco-2 cells incubated with Escherichia coli K12 lipopolysaccharide at a final concentration of 1,000 ng/mL. Assays included three to five biological experiments with technical duplicates.

Untargeted metabolomics

Culture supernatants were processed for metabolite extraction and analysis by the LEM facility (CEA Paris-Saclay, France), as previously described (ref. 30). Metabolites were extracted using methanol-assisted protein precipitation and analyzed via liquid chromatography-high-resolution mass spectrometry (LC-HRMS), combining reversed-phase (C18) and HILIC methods (excluding fatty acids and complex lipids). LC-HRMS was performed on an Ultimate 3000 system coupled to Exactive or Q-Exactive mass spectrometers (Thermo Fisher Diagnostics SAS, Dardilly, France), using electrospray ionization in both positive (C18) and negative (HILIC) modes. Metabolite annotation was done using an internal spectral database, and peak integration was manually reviewed. Peak areas were normalized to QC samples using the LOESS algorithm to correct for batch effects (ref. 31). Final data were expressed in arbitrary units, representing chromatographic peak areas per metabolite.

Short-chain fatty acid analyses

Gas chromatography was used to measure acetate, propionate, and butyrate levels in the cell-free culture supernatant samples of the wild-type (WT) and Δhbd mutant strains, as described previously (ref. 16).

Statistical analysis

Experimental data were analyzed using GraphPad Prism 8.0.1 (La Jolla, CA, USA). Comparisons between CB1002 HBD and WT strains were performed using the unpaired, Welch-corrected Student’s t-test. Metabolite data were analyzed using both univariate and multivariate approaches in R (R Core Team, 2022; https://www.R-project.org). The relationship between gene expression and metabolites was tested using Spearman’s rank correlation coefficient. The R software package “corrplot” (Version 0.95) was used for correlation visualization (https://github.com/taiyun/corrplot). Statistical significance was set at P < 0.05.

RESULTS AND DISCUSSION

De novo transcriptome assembly

RNA-Seq produced high-quality data sets averaging 24 million reads per sample, with 93% of bases scoring above Q30, and 99% mapping to the CB1002 WT genome. Read normalization minimized false positives (Fig. S1), and a clear separation was observed between the CB1002 HBD and WT samples (Fig. S2). To analyze transcriptomic changes, the CB1002 WT genome was sequenced and annotated, revealing a genome size of 4.57 Mb with 28.52% G+C content, which is consistent with the reference strain VPI3266^T. A total of 4,198 coding sequences were annotated.

Identification of DEGs

Of the annotated genes, 2,054 exhibited differential expression, 670 of which were significantly different (P < 0.05). This included 328 upregulated genes and 342 downregulated genes. KEGG enrichment analysis revealed that these genes were involved in metabolic and biosynthetic pathways (40%), ABC transporters (3%), two-component systems (2%), quorum sensing (1%), and chemotaxis and motility (1%) (Fig. S3). According to the COG enrichment analysis, most of the matched genes were associated with unknown functions (29%), replication and transcription and translation (21%), carbohydrate metabolism (7%), energy production (5%), and ion transport (6%) (Fig. S4).

To investigate the observed loss of enteropathogenicity observed in the quail NEC model following the hbd mutation (ref. 16), we specifically focused on potential virulence genes that were differentially expressed between CB1002 HBD and WT. These included genes associated with bacterial chemotaxis, flagellar and pilus assembly, quorum sensing, metabolism, cell wall synthesis, and oxidative stress response. Table 1 lists the 61 genes of interest that we identified.

Bacterial chemotaxis, flagellar, and pilus assembly-associated DEGs

Of the 61 DEGs, 14 belong to the chemotaxis-related KEGG classes, including cheY, cheX, cheV, nine mcpB genes, pilT, and pilB (Table 1). In CB1002 HBD, cheY was downregulated, while cheX and cheV were upregulated. The methyl-accepting chemotaxis McpB proteins interact with the CheA/CheY histidine protein kinases two-component signaling system, which controls flagellar movement in response to repellents. cheX and cheV, which encode coupling proteins, modulate CheY activity and affect motility (ref. 32). The flagellar assembly genes motA and fliB were also repressed in CB1002 HBD (Table 1). These data are consistent with the reduced swimming and swarming motility previously observed in CB1002 HBD compared with WT (ref. 14). Although not statistically significant (P > 0.05), 32 other flagellar assembly genes showed differential expression (Table S1). Additionally, pilB and pilT were repressed in CB1002 HBD. PilB is an ATPase involved in type IV pilus assembly, and PilT provides the power for pilus retraction, which is required for twitching motility (ref. 33). These findings are consistent with earlier observations of reduced twitching motility in CB1002 HBD (ref. 16). Chemotaxis and motility (flagella and pili) are tightly linked and play multiple roles in bacterial behavior, including adhesion, biofilm formation, autoaggregation, and host interactions. All of these processes are linked to virulence (ref. 34).

Quorum sensing-associated DEGs

As shown in Table 1, we found that oppA, oppB, oppC, and oppD were downregulated in CB1002 HBD. These genes form an operon that encodes membrane-associated proteins of an ATP-binding cassette transporter (ref. 35). A variety of roles have been described for the bacterial Opp systems. In particular, they are involved in short peptide transport, including amino acid uptake and peptide-mediated signaling, through RRNPP type quorum sensing systems that regulate virulence gene expression (ref. 35). In addition, expression of the argCJBDF operon was repressed in CB1002 HBD (Table 1). agrD encodes a precursor peptide that is processed by AgrB into an autoinducing peptide (AIP), which activates the sensor kinase AgrC and initiates a response regulator cascade for gene expression (ref. 35). The Agr system also plays a fundamental role in quorum sensing, regulating the spatiotemporal expression of virulence factors, metabolic flexibility, and survival (ref. 36, ref. 37). In C. difficile, the activation of argD and the downstream genes argB and argC may participate in toxin production, flagellar biosynthesis, and cyclic diguanylate signaling (ref. 36, ref. 37). Several other quorum sensing-related genes (dppD, ciaR, ciaH, and nisin biosynthesis-related genes) were differentially expressed but did not reach statistical significance (Table S1).

Metabolism-associated DEGs

We confirmed that the hbd gene was significantly repressed in CB1002 HBD, as were the ptb and buk genes. The ptb gene encodes a phosphotransbutyrylase, and the buk gene encodes a butyrate kinase. These are the enzymes involved in the production of butyrate (ref. 38) (Table 1). Several sugar metabolism genes were also downregulated, including lacZ, cbgA, agaA, ilvC, ilvB, and lctP. LacZ, CbgA, and AgaA are glucosidases involved in sugar catabolism. ilvB and ilvC encode acetohydroxy acid synthase (IlvB) and ketoacid reductoisomerase (IlvC), respectively, which are involved in the biosynthesis of branched-chain amino acids (BCAAs). The hbd mutant may shift to use BCAAs as a redox-neutral overflow or for biosynthesis due to carbon/nitrogen imbalance. lctP encodes a lactate permease, whose repression suggests a limitation in lactate uptake. It is likely that CB1002 HBD downregulates certain sugar utilization and biosynthesis genes to avoid overloading its compromised fermentation network and maintain energy/redox balance. Virulence-related functions, such as galactosidase activity, lactate, and amino acid metabolism, have been directly and indirectly linked to pathogenicity (ref. 39ref. –ref. 42).

In contrast, galR, galTKE, glpKF, and ldh-1 were upregulated. Due to the hbd mutation, which blocks the butyrate pathway and impairs NADH oxidation, the upregulation of ldh-1 may represent a compensatory response to restore redox balance. The galTKE operon, which is regulated by the GalR regulator, encodes the enzymes that are responsible for the catabolizing galactose via the Leloir pathway. This allows glucose-1-phosphate to enter glycolysis and be oxidized to pyruvate, enabling bacteria to survive and proliferate using alternative energy sources. This is consistent with the overexpression of the glycerol kinase GlpK and the glycerol uptake facilitator GlpF, which suggests a potential shift toward 1,3-propanediol production. Taken together, these data suggest a redirection of the bacterial metabolism affecting its biological properties and virulence.

Bacterial cell wall-associated DEGs

The genes involved in the peptidoglycan biosynthesis were significantly differentially expressed in the CB1002 HBD strain (Table 1). glmU and murA were repressed in CB1002 HBD. These genes encode an N-acetylglucosamine-1-phosphate uridyltransferase and a UDP-N-acetylglucosamine-1-carboxyvinyltransferase, respectively. They are involved in the cytoplasmic steps of peptidoglycan UDP-N-acetylglucosamine biosynthesis (ref. 43). The ddlB gene, encoding a D-alanine ligase that catalyzes the incorporation of D-alanine into the bacterial cell wall lipoteichoic acids, was overexpressed in the CB1002 HBD strain (ref. 43). In C. butyricum, D-alanylation increases the hydrophobicity of the cell wall (ref. 44). Previously, we have reported that CB1002 HBD has higher surface hydrophobicity than CB1002 WT (ref. 16). The genes penP and vanR, encoding for a putative penicillin-binding protein and a putative cytoplasmic response regulator of the two-component system VanSR, were also differentially expressed in CB1002 HBD (Table 1). PenP and VanSR are associated with peptidoglycan modification (ref. 45). Interestingly, comEC and comB, which encode proteins involved in bacterial competence, were overexpressed in the mutant, suggesting possible cell wall modifications (ref. 46, ref. 47). These modifications may contribute to changes in bacterial virulence through adhesion, biofilm formation, stress response, and resistance to cationic antimicrobial peptides (ref. 43). Of note, the dltA, dltC, dltD, and vanY were differentially expressed, though not significantly (Table S1).

Oxidative stress-associated DEGs

Oxidative stress-related genes (sodA, nirA_1, and ldh) were significantly overexpressed in the CB1002 HBD strain compared with the WT strain (Table 1). sodA encodes a superoxide dismutase that detoxifies superoxide and works with catalase to scavenge ROS (ref. 48). nirA enables denitrification under nitrosative stress (ref. 49), and ldh supports redox balance (ref. 39). The hbd mutant, which is unable to oxidize NADH via butyrate synthesis, may shift to a redox compensation mode. These genes have been associated with increased survival and virulence under oxidative stress (ref. 50). Conversely, trxB and bcp were upregulated, while trxA was downregulated. trxA encodes an electron carrier of the thioredoxin system that donates electrons to reduce disulfide bonds in proteins. trxB encodes a thioredoxin reductase that regenerates reduced thioredoxin using NADPH (ref. 51). This expression pattern of the thioredoxin system may indicate that the bacteria are redirecting NADPH and defending against oxidative stress. This is consistent with the upregulation of bcp, which encodes a peroxiredoxin that detoxifies peroxides. It is also consistent with sodA expression, indicating that the bacteria is defending against both superoxide and peroxides. gloA, which encodes glyoxalase I for detoxifying the byproduct of central carbon metabolism, methylglyoxal, was also overexpressed (ref. 52). In the CB1002 HBD mutant strain, where butyrate production is blocked, causing pyruvate and upstream intermediates to accumulate. Overexpression of gloA helps prevent carbon-induced bacterial viability loss by detoxifying methylglyoxal. These systems are also involved in resistance to stress (ref. 51).

Protein-protein network analysis

We analyzed the interactions among the 61 DEGs of interest. The interaction network for C. butyricum comprised 29 nodes and revealed eight clusters (Fig. 1A), including arginine, butyrate, galactose metabolism, bacterial chemotaxis, amino sugar biosynthesis (valine, leucine, and isoleucine biosynthesis), glycerol dehydrogenase and carbohydrate kinase, and oxidative stress-related proteins. We also compared the PPI networks of C. butyricum and Clostridium neonatale. C. neonatale is the closest species relative to C. butyrium within the genus Clostridium sensu stricto (cluster I) (ref. 53, ref. 54). Additionally, C. neonatale has previously been isolated from NEC cases (ref. 55). Similar to C. butyricum, deletion of the hbd gene in the C. neonatale 250.09 clinical strain abolished enteropathogenicity in the quail NEC animal model (ref. 16). Analysis of the C. neonatale PPI network identified 19 nodes with clusters similar to those of C. butyricum, including valine, leucine, and isoleucine biosynthesis, galactose and butyrate metabolism, and oxidative stress-related proteins (Fig. 1B). Thirty-three proteins overlapped between the two species. The smaller C. neonatale network compared with C. butyricum suggests species-specific differences. However, the shared protein networks point to conserved pathways, which aligns with their similar enteropathogenic behavior in the quail NEC model (ref. 16).

Gene expression correlations

Figure 2 highlights distinct gene clusters whose expression was significantly altered following hbd deletion, suggesting disruption of normal metabolism. Overexpression of galE, galR, galK, galT, agaA, sodA, ldh, bcp, mcpB, fliB, comB, grxC, and nirA in CB1002 HBD suggests a shift in carbon flux, changes in cell wall composition, and increased use of alternative sugars and electron acceptors. Upregulation of nirA indicates an adjustment in nitrogen metabolism, possibly involving nitrate. This metabolic reprogramming may increase toxic byproduct formation, as reflected by coexpression of sodA, bcp, grxC, and trxB, which in turn may trigger gloA-mediated detoxification. Correlated expression of the opp operon, arg operon, motA, pilT, and mcpB points to a coordinated response to nutrient stress, enhancing nutrient uptake and environmental sensing involving quorum sensing. Upregulation of mcpB, fliB, and che genes supports a modification of chemotactic response to unfavorable conditions caused by the loss of the β-hydroxybutyryl-CoA dehydrogenase function. Notably, gloA, nirA, and cheV are also regulated by general stress-response pathways (ref. 51).

Gene expression and metabolite correlations

We confirmed that the culture supernatant of CB1002 HBD contained no detectable butyrate and produced acetate and propionate at levels similar to those previously reported for CB1002 WT (ref. 14). A total of 28 metabolites were identified through non-targeted metabolomics analysis of the same culture supernatant used for RNA-Seq analysis (Table S2). Figure 3 illustrates some of the correlations identified between metabolites and gene expression. Consistent with the overexpression of galE, galR, galK, galT, agaA, sodA, ldh, bcp, mcpB, fliB, comB, grxC, and nirA, we identified metabolites associated with a metabolic shift in CB1002 HBD. CB1002 HBD had higher levels of lactic acid, acetoin, oxobutyric acid, and sugars, which are associated with fermentation flux products. This is indicative of pyruvate overflow and shunting to multiple pathways (Table S2). We also identified indole lactic acid and phenyl lactic acid, derivatives of tryptophan and phenylalanine, respectively. These are associated with the degradation of amino acids during overflow metabolism under anaerobic and fermentative conditions. Notably, indole-lactic acid and phenyl-lactic acid have been reported to act on intestinal epithelial and immune cells, where they can signal through receptors, such as the aryl hydrocarbon receptor to dampen NF-κB-dependent pro-inflammatory cytokine production, enhance epithelial barrier function, and promote anti-inflammatory immune phenotypes (ref. 56). Thus, the accumulation of these metabolites in CB1002 HBD supernatant suggests that, in vivo, this strain may expose intestinal tissues to a qualitatively different metabolite profile with the potential to modulate epithelial responses and the local mucosal immune environment. The presence of the fatty acid derivative hexadecanedioic acid suggests a shift in lipid metabolism that may be associated with membrane remodeling, oxidative stress, or carbon overflow in fatty acid oxidation processes. The presence of glutamine, N-acetylserine, indole lactic acid, and phenyl lactic acid suggests rewiring of sulfur metabolism or oxidative stress. Therefore, in the context of a metabolic shift, intermediate byproducts accumulate, leading to stress responses. This is consistent with the identification of metabolites associated with nucleotide and polyamine-related metabolites, such as deoxyadenosine, cAMP, N8-acetylspermidine, and N6-acetyllysine, which are associated with signaling and cell stress response. This agrees with the overexpression of sodA, bcp, grxC, and trxB, and gloA, as well as increased oxidative stress and redox imbalance. Meanwhile, the identification of the following flavonoids—apigenin, genistein, pelargonidin, naringenin, and glycitein—may reflect mechanisms that help bacteria adapt to bacterial metabolic and environmental changes. This is consistent with the differential expression of genes that reflect chemotaxis, flagellar, and pilus processes.

In terms of gene-metabolite function associations, based on the KEGG MAPPER database and the KEGG compound ID, we confirmed that the identified metabolites were associated with map01100_Metabolic pathways (n = 23 metabolites), map01110_Biosynthesis of secondary metabolites (n = 14), map01120_Microbial metabolism in diverse environments (n = 8), map01230_Biosynthesis of amino acids (n = 6), map00941_Flavonoid biosynthesis (n = 3) and map00942_Isoflavonoid biosynthesis (n = 4), map01240_Biosynthesis of cofactors (n = 3), and map00230_Purine metabolism (n = 3).

Inflammatory property analysis

Intestinal inflammation has been identified as a significant contributing factor to the development of NEC (ref. 1, ref. 2). Furthermore, bacterial involvement in local inflammation has been observed in NEC. In this study, we investigated the in vitro immunostimulatory properties of CB1002 strains on PBMCs. The levels of various cytokines, including IL-10, IL-17, IL-22, and IFNγ, in the culture supernatant of PBMCs after co-incubation with CB1002 HBD and WT live bacteria and bacterial components were not statistically significant, except for IL-10 production, which was significantly lower in the presence of bacterial debris from CB1002 HBD (Fig. 4). This finding suggests that the hbd deletion may alter bacterial cellular stimulation. PBMC cells from the control condition produce a higher concentration of IL-10 compared with the other experiment conditions, without statistical difference.

Redox stress analysis

Oxidative stress has also been described as a contributing factor to NEC development (ref. 1, ref. 2). Therefore, we analyzed the production levels of ROS, O₂^•–, and ONOO^– production in Caco-2 cells after stimulation with the different bacteria samples. No statistically significant differences in total ROS, O₂^•–, or ONOO⁻ levels were observed between CB1002 HBD and CB1002 WT under any of the conditions tested (Fig. 5). Thus, the butyrate-deficient mutant does not appear to differentially contribute to epithelial redox activation relative to the WT strain in the Caco-2 model. When compared with their respective unstimulated controls, however, all bacterial conditions tended to show a reduction in bulk ROS together with a concomitant increase in ONOO⁻, with significant ONOO⁻ increases detected for the WT DB, WT S, and HBD S conditions (Fig. 5A and C). These findings indicate that exposure to CB1002 preparations is associated with a redistribution of reactive species from classical ROS toward reactive nitrogen species, while the overall pattern of redox modulation remains similar between the WT and HBD strains in this epithelial model. Although the ONOO⁻ increase observed here was modest and comparable between CB1002 WT and HBD, a similar shift from ROS to peroxynitrite in an immature, inflamed intestine could still amplify nitrosative injury, promote enterocyte apoptosis, and impair barrier function even without large changes in bulk ROS. In this sense, CB1002 strains may not differ in how they modulate epithelial redox pathways relative to each other, but they do engage host redox signaling in ways that, in a susceptible neonatal gut, could intersect with established NO/iNOS–peroxynitrite-driven mechanisms of NEC (ref. 57).

Limitations

One limitation of our study is that the RNA-Seq analysis identified less than 50% of the expressed genes from the 4,192 annotated genes in the sequenced CB1002 WT genome. This is due to the automatic annotation of the CB1002 WT genome, which identified only 34% of the gene products and assigned an unknown function to the remaining 66% of the genes. Manually annotating and curating the CB1002 genome may allow for a better analysis of our data. Another limitation is that the experimental data and analysis are limited to the exponential growth phase of bacterial cultures in a rich liquid medium. Although we confirmed previous experimental data, this descriptive in vitro work does not account for interactions with the gut microbiota and requires further functional validation through additional experimental studies. One possible way to overcome this limitation would be to culture human gut microbiota in vitro using devices such as the Shime or Mipro systems and then analyze how the composition and functions of the microbiota change upon exposure to C. butyricum strains.

Conclusion

C. butyricum is an intestinal commensal bacterium that has been observed to behave similarly to opportunistic pathogens that cause NEC. However, no classical repertoire of virulence genes has yet been identified that would allow for the unambiguous identification of pathogenic strains. In this study, we identified differential gene expression in an avirulent hbd mutant strain, which suggests that the bacterial cellular processes are reprogrammed to compensate for the consequences of the hbd deletion. These changes were consistent with different levels of various metabolites correlated with gene expression and protein interaction networks. If these cellular processes are indeed associated with the biological activities of C. butyricum, they may be implicated in its pathogenesis. Experimental data showed that the hbd deletion had no effect on in vitro inflammatory or oxidative activity when cellular models were exposed to the mutant strain. These results suggest that these processes may not contribute to the development of NEC-like lesions observed in an animal model of NEC. In conclusion, we identified several genes in C. butyricum that are promising targets for further study, particularly in the context of NEC pathogenesis. Generating mutants and in-depth characterizing those in vitro and in vivo settings will facilitate identifying the critical bottlenecks hindering invasion and survival. This will allow us to assess immune susceptibility and evaluate the enteropathogenicity of selected mutants in an NEC animal model. This will pave the way for validating new hypotheses about the role of C. butyricum and other related species as potential pathogens.

References

1. Necrotizing enterocolitis.. N Engl J Med, 2011. [DOI | PubMed]
2. Necrotizing enterocolitis: long term complications.. Curr Pediatr Rev, 2019. [DOI | PubMed]
3. Global incidence of necrotizing enterocolitis: a systematic review and meta-analysis.. BMC Pediatr, 2020. [DOI | PubMed]
4. Is intestinal cell death in necrotising enterocolitis assorted and multifarious? A special focus on risk factors and their pathogenic mechanisms.. EMJ Gastroenterol, 2023. [DOI]
5. Intestinal dysbiosis in preterm infants preceding necrotizing enterocolitis: a systematic review and meta-analysis.. Microbiome, 2017. [DOI | PubMed]
6. Infectious causes of necrotizing enterocolitis.. Clin Perinatol, 2015. [DOI | PubMed]
7. Clostridia and necrotizing enterocolitis in preterm neonates.. Anaerobe, 2019. [DOI | PubMed]
8. 8 Ruan J. 2013. Bergey’s manual of systematic bacteriology, p 521–530. In The study of Actinomycetes systematic in China, second edition. Vol. 53.
9. Clostridia in premature neonates’ gut: incidence, antibiotic susceptibility, and perinatal determinants influencing colonization.. PLoS One, 2012. [DOI | PubMed]
10. Clostridium butyricum strains and dysbiosis linked to necrotizing enterocolitis in preterm neonates.. Clin Infect Dis, 2015. [DOI | PubMed]
11. Nutritional strategies and gut microbiota composition as risk factors for necrotizing enterocolitis in very-preterm infants.. Am J Clin Nutr, 2017. [DOI | PubMed]
12. Multidisciplinary evaluation of Clostridium butyricum clonality isolated from preterm neonates with necrotizing enterocolitis in South France between 2009 and 2017.. Sci Rep, 2019. [DOI | PubMed]
13. Acute necrotizing enterocolitis of preterm piglets is characterized by dysbiosis of ileal mucosa-associated bacteria.. Gut Microbes, 2011. [DOI | PubMed]
14. Community analysis of bacteria colonizing intestinal tissue of neonates with necrotizing enterocolitis.. BMC Microbiol, 2011. [DOI | PubMed]
15. Evidence for clostridial implication in necrotizing enterocolitis through bacterial fermentation in a gnotobiotic quail model.. Pediatr Res, 2005. [DOI | PubMed]
16. Neonatal necrotizing enterocolitis: Clostridium butyricum and Clostridium neonatale fermentation metabolism and enteropathogenicity.. Gut Microbes, 2023. [DOI | PubMed]
17. Analysis of the core genome and pangenome of Clostridium butyricum. Genome, 2021. [DOI | PubMed]
18. Pangenome analyses of Clostridium butyricum provide insights into its genetic characteristics and industrial application.. Genomics, 2024. [DOI | PubMed]
19. Experimental cecitis in gnotobiotic quails monoassociated with Clostridium butyricum strains isolated from patients with neonatal necrotizing enterocolitis and from healthy newborns.. Infect Immun, 1989. [DOI | PubMed]
20. “Sequana”: a set of snakemake NGS pipelines.. JOSS, 2017. [DOI]
21. Fastp: an ultra-fast all-in-one FASTQ preprocessor.. Bioinformatics, 2018. [DOI | PubMed]
22. Fast gapped-read alignment with Bowtie 2.. Nat Methods, 2012. [DOI | PubMed]
23. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2.. Genome Biol, 2014. [DOI | PubMed]
24. Core-, pan- and accessory genome analyses of Clostridium neonatale: insights into genetic diversity.. MicrobGenom, 2022. [DOI]
25. MicroScope: an integrated platform for the annotation and exploration of microbial gene functions through genomic, pangenomic and metabolic comparative analysis.. Nucleic Acids Res, 2019. [DOI]
26. eggNOG-mapper v2: functional annotation, orthology assignments, and domain prediction at the metagenomic scale.. Mol Biol Evol, 2021. [DOI | PubMed]
27. The STRING database in 2023: protein-protein association networks and functional enrichment analyses for any sequenced genome of interest.. Nucleic Acids Res, 2023. [DOI | PubMed]
28. Clostridial strain-specific characteristics associated with necrotizing enterocolitis.. Appl Environ Microbiol, 2018. [DOI | PubMed]
29. NADPH oxidase is the major source of placental superoxide in early pregnancy: association with MAPK pathway activation.. Sci Rep, 2019. [DOI | PubMed]
30. Annotation of the human serum metabolome by coupling three liquid chromatography methods to high-resolution mass spectrometry.. Journal of Chromatography B, 2014. [DOI]
31. Procedures for large-scale metabolic profiling of serum and plasma using gas chromatography and liquid chromatography coupled to mass spectrometry.. Nat Protoc, 2011. [DOI | PubMed]
32. The effect of bacterial chemotaxis on host infection and pathogenicity.. FEMS Microbiol Rev, 2018. [DOI]
33. A review on pilus assembly mechanisms in Gram-positive and Gram-negative bacteria.. Cell Surf, 2022. [DOI | PubMed]
34. Multiple functions of flagellar motility and chemotaxis in bacterial physiology.. FEMS Microbiol Rev, 2021. [DOI | PubMed]
35. The oligopeptide ABC-importers are essential communication channels in Gram-positive bacteria.. Res Microbiol, 2019. [DOI | PubMed]
36. Editorial: quorum-sensing in Gram-positive pathogens – mechanisms, role in infection, and potential as a therapeutic target.. Front Cell Infect Microbiol, 2023. [DOI | PubMed]
37. Combined and distinct roles of Agr proteins in Clostridioides difficile 630 sporulation, motility, and toxin production.. mBio, 2020. [DOI | PubMed]
38. Regulation of carbon and electron flow in Clostridium butyricum VPI 3266 grown on glucose-glycerol mixtures.. J Bacteriol, 2001. [DOI | PubMed]
39. Microbial lactate utilization: enzymes, pathogenesis, and regulation.. Trends Microbiol, 2014. [DOI | PubMed]
40. Characterization of novel beta-galactosidase activity that contributes to glycoprotein degradation and virulence in Streptococcus pneumoniae.. Infect Immun, 2010. [DOI | PubMed]
41. Host glycan sugar-specific pathways in Streptococcus pneumoniae: galactose as a key sugar in colonisation and infection [corrected].. PLoS One, 2015. [DOI | PubMed]
42. Branching out: alterations in bacterial physiology and virulence due to branched-chain amino acid deprivation.. mbio, 2018. [DOI | PubMed]
43. Lipoteichoic acid synthesis and function in gram-positive bacteria.. J Bacteriol, 2014. [DOI | PubMed]
44. The functional dlt operon of Clostridium butyricum controls the D-alanylation of cell wall components and influences cell septation and vancomycin-induced lysis.. Anaerobe, 2015. [DOI | PubMed]
45. Genetic mechanisms of vancomycin resistance in Clostridioides difficile: A Systematic Review.. Antibiotics (Basel), 2022. [DOI | PubMed]
46. Genetic competence in Helicobacter pylori: mechanisms and biological implications.. Res Microbiol, 2000. [DOI | PubMed]
47. Topology and membrane interaction of Helicobacter pylori ComB proteins involved in natural transformation competence.. Int J Med Microbiol, 2003. [DOI | PubMed]
48. SOD enzymes and microbial pathogens: surviving the oxidative storm of infection.. PLoS Pathog, 2016. [DOI | PubMed]
49. NirA is an alternative nitrite reductase from Pseudomonas aeruginosa with potential as an antivirulence target.. mbio, 2021. [DOI | PubMed]
50. Bacterial Stress responses during host infection.. Cell Host Microbe, 2016. [DOI | PubMed]
51. Oxidative stress, protein damage and repair in bacteria.. Nat Rev Microbiol, 2017. [DOI | PubMed]
52. Bacterial responses to glyoxal and methylglyoxal: reactive electrophilic species.. Int J Mol Sci, 2017. [DOI | PubMed]
53. Clostridium neonatale sp. nov. linked to necrotizing enterocolitis in neonates and a clarification of species assignable to the genus Clostridium (Prazmowski 1880) emend. Lawson and Rainey 2016.. Int J Syst Evol Microbiol, 2018. [DOI | PubMed]
54. An atypical Clostridium strain related to the Clostridium botulinum group III strain isolated from a human blood culture.. J Clin Microbiol, 2014. [DOI | PubMed]
55. An outbreak of necrotizing enterocolitis associated with a novel Clostridium species in a neonatal intensive care unit.. Clin Infect Dis, 2002. [DOI | PubMed]
56. Indole-3-lactic acid, a metabolite of tryptophan, secreted by Bifidobacterium longum subspecies infantis is anti-inflammatory in the immature intestine.. Pediatr Res, 2020. [DOI | PubMed]
57. Redox chemistry: implications for necrotizing enterocolitis.. Int J Mol Sci, 2024. [DOI | PubMed]