Constructing a comprehensive Desulfovibrio database of 2658 genomes
For a comprehensive characterization of Desulfovibrio spp., we applied a MAG recovery approach. First, we analysed 11,017 samples from 10 diverse human metagenomic cohort studies, resulting in 1921 new Desulfovibrio high-quality MAGs (Supplementary Data 1+2). These cohorts span different diseases (including IBD and CRC), ages (young, elderly and centenarians) and geographical areas (i.e., Southern China, Japan, Nepal, Tanzania, the Netherlands, the US, Germany, France and Israel). The majority of the recovered MAGs (1726/1921; 89.85%) came from two Dutch cohorts. High Desulfovibrio prevalence was observed in both Dutch studies (LifeLineNext38 mother & infant study, prevalence in adults = 26.67%; Dutch microbiome cohort39, prevalence = 20.86%), while several Chinese populations showed a high degree of variation in prevalence (max prevalence: 60%, min prevalence: 5%, Supplementary Data 3). Interestingly, sequencing depth was not a major driver of Desulfovibrio spp. detection (Supplementary Figs. 1, 2). In addition, 578 Desulfovibrio genomes were recovered from MetaRefSGB40 (Supplementary Data 1), a resource that includes a collection of MAGs from > 100k metagenomes spanning different body sites, ages, countries and lifestyles. Furthermore, we identified 442 Desulfovibrio genomes in publicly available databases, including NCBI RefSeq, the Unified Human Gastrointestinal Genome (UHGG), HRGM, IMG/M and MGnify. Due to recent re-classification efforts, we re-annotated all of these genomes using GTDB-tk annotations41. This led to the re-assignment of 31 genomes to other genera within Desulfovibrionaceae (Suppl. Data 4), which were subsequently excluded from our analysis. Overall, this genome database represents Desulfovibrio diversity recovered from human, animal, and environmental samples (Fig. 1a).
Fig. 1: Genomic analysis of 2658 Desulfovibrio genomes.
a We reconstructed a Desulfovibrio database consisting of 1921 high-quality metagenome-assembled genomes (MAGs) identified from the analysis of 11,017 metagenomic samples from 10 human cohort studies with different health conditions (e.g. healthy individuals, inflammatory bowel diseases [IBD], colorectal cancer [CRC], and other diseases) and geographical areas, including Southern China, Japan, Nepal, Tanzania, the Netherlands, the US, Germany, France and Israel. In addition, 578 Desulfovibrio genomes were retrieved from MetaRefSGB, containing MAGS from > 100k metagenomes40, and major public genome databases. Furthermore, we isolated and obtained 24 Desulfovibrio strains from human cohort samples and public culture collections, which were subsequently sequenced. Additionally, 11 human isolate genomes were downloaded from NCBI RefSeq. All available genomes were merged and duplicate entries removed, resulting in 2658 Desulfovibrio genomes. HQ MAGs – high-quality metagenomic assembled genomes. b Phylogenetic tree of all Desulfovibrio genomes based on mash distance. Rings indicate host phenotype (diseased vs. healthy, “others” refers to genomes from animals, natural environments or unknown origin), sample source (human vs. animal vs. environment) and genome source (newly identified MAG vs. newly sequenced isolate or database entry). Stars indicate strains for which isolates are available. Species names follow GTDB taxonomy. c Overview of the number of genomes within each species. The species with the most genomes is D. sp900556755 (proposed name: Desulfovibrio elongatum sp. nov.). followed by D. piger_A, D. piger, D. sp900319575, and D. fairfieldensis. d Hierarchical clustering of 86 genomes from uncharacterized Desulfovibrio spp. based on mash distance (cutoff < 5%). Colors indicate Desulfovibrio clades: clade I (blue), clade II (orange), and clade III (green).
As isolates are crucial for any type of functional validation experiments, we next focused on culturing and procuring human Desulfovibrio isolates. We were able to isolate 2 strains from IBD and CRC patients and 7 strains from healthy individuals. In addition, 15 clinical isolates were acquired from microbial culture repositories, leading to a total of 24 newly sequenced human isolates representing 8 major Desulfovibrio species and one uncharacterized Desulfovibrio spp. (Fig. 1a, Suppl. Data 5).
After merging all available genomes, they were de-replicated (99.9% average nucleotide identity [ANI] cutoff42,43), resulting in a final curated database of 2,658 unique, high-quality Desulfovibrio genomes, spanning 37 species (including new candidate species, Fig. 1b, c), recovered from 174 metagenomic studies (including 101 human microbiome studies) from 32 countries, representing 90 different diseases and health conditions (Supplementary Data 1, 6). Most genomes belonged to Desulfovibrio elongatum sp. nov. (i.e., D. sp900556755 based on GTDB-tk) (n = 1528, 57.49%), which constituted the major species represented in this database. Other main species (containing > 50 genomes) were D. piger_A, D. piger, D. sp900319575, and D. fairfieldensis (Fig. 1b, c). In summary, this comprehensive Desulfovibrio database with its host-associated information will serve as the basis for exploring genomic diversity and strain- and gene-specific disease signals, which can then be further investigated using the representative clinical isolates.
Desulfovibrio spp. are found across different human diseases and health conditions
To investigate Desulfovibrio species diversity and phylogeny, we calculated pairwise mash distances between all genomes (Fig. 1b). While Desulfovibrio is often associated with disease1,2,3,4, they were also prevalent in healthy individuals with 686 Desulfovibrio genomes recovered from healthy samples versus 1503 recovered from patient samples. However, no obvious clustering by disease phenotype was observed, indicating that overall genome similarity is insufficient and detailed analyses are required for the identification of disease-specific functional traits. Interestingly, Desulfovibrio sp910575755 and Desulfovibrio sp910584445, two murine species, clustered separately, potentially reflecting host adaptation, while 14 environmental Desulfovibrio genomes (originating from soil, sand, sediment, freshwater, wastewater, bioreactor and agricultural waste), were interspersed with genomes from human samples.
Furthermore, our analysis identified 3 potential novel species. Among the 102 Desulfovibrio genomes that could not be assigned to any existing species, 76 genomes clustered into 3 clades (mash distance < 5%44,45,46, Fig. 1d). Desulfovibrio clade I contained 70 genomes, including 58 MAGs from the two Dutch cohorts, while the other two clades were smaller with 3 MAGs from the Dutch cohorts and 3 MAGs from a Hadza hunter population, respectively. The detection of these novel species highlights that Desulfovibrio spp. are highly diverse and that much of this diversity remains to be uncovered.
Desulfovibrio spp. differ in cell motility, ureases, NADH dehydrogenase and oxidoreductase genes
We next identified differential genetic features in different Desulfovibrio species. Several clusters of orthologous groups (COGs) emerged, in particular COGs related to cell motility (Fig. 2a), suggesting species-specific difference in motility. Especially a flagellin gene (COG1344) was present in all genomes for 28 Desulfovibrio species, including D. fairfieldensis (117/117), D. desulfuricans (42/42), D. legallii (42/42) and Desulfovibrio clade III (3/3, Fig. 1d), while completely absent in other species, including D. piger, Desulfovibrio elongatum sp. nov., and Desulfovibrio clade I and II. This indicates increased motility and potential microbe-host interactions47 of D. fairfieldensis, D. desulfuricans, D. legallii and Desulfovibrio clade III.
Fig. 2: Comparison of functional potential across different Desulfovibrio species.
a The heatmap shows that the functional capabilities of Desulfovibrio species differ in clusters of orthologous groups (COGs) related to cell motility, urease, NADH dehydrogenase, and Na+-translocating Fd/NADH oxidoreductase. For cell motility, only COGs present in at least one species with a minimum prevalence of > 10% were selected for the visualization. The flagellin COG is indicated in red and other flagellar-related COGs highlighted in brown. b, c The heatmap and barplots show the genomic diversity in (b) antimicrobial resistance genes and (c) virulence factors. Genes were identified using ABRicate with a cutoff of > 60% identify and > 80% coverage. Species names follow GTDB taxonomy and represent genome-based species-level clusters. Desulfovibrio spp. refers to multiple Desulfovibrio genomes assigned to the Desulfovibrio genus but cannot be resolved to a known species (based on GTDB).
In addition, we identified differences in genes encoding urease, NADH dehydrogenase complex, and Na+-translocating Rnf complex (Fig. 2a). Urease is broadly involved in nitrogen metabolism and environmental adaptation, while it also contributes to the virulence in some pathogenic bacteria48. For example, urease in Helicobacter pylori produces ammonia at levels toxic to host epithelial cells and can neutralize stomach acidity facilitating colonization49. More than 95% genomes of D. desulfuricans, D. legallii, Desulfovibrio sp910584445, and Desulfovibrio sp900547595 encoded the ureABCEFH operon required for urease synthesis and activity. In contrast, these genes were absent in D. piger and D. piger_A genomes. Furthermore, the nuo gene clusters (nuoABDHJKLMN) encoding proton-pumping type I NADH dehydrogenase was present in several species, including D. desulfuricans, D. legallii, and D. intestinalis, but absent in D. porci, D. sp900319575 and several other species (Fig. 2a). This NADH dehydrogenase complex was previously identified in Solidesulfovibrio magneticus (formerly Desulfovibrio magneticus)50. It connects NADH oxidation to proton translocation, involving energy transfer reactions during sulfate respiration51. Lastly, we also detected a rnfCDGEFB gene cluster encoding Na+-translocating Rnf complex in > 90% genomes of 18 Desulfovibrio species, including Desulfovibrio clade I, D. piger, and D. piger_A. The Rnf complex is homologous to the electrogenic NADH:ubiquinone reductase complex (Nqr), a sodium-pumping NADH dehydrogenase52. Interestingly, the nuo and rnf operon co-occurred in Desulfovibrio sp905214305. These results indicate that Desulfovibrio species encode different metabolic pathways related to energy transfer involving sulfate respiration and identify functional traits (i.e., ure) that have been shown to be involved in host-associated adaptation in pathogenic bacteria.
Antimicrobial resistance genes and virulence factors differ across Desulfovibrio species
To further explore disease-associated gene signals, we next investigated virulence factors and antimicrobial resistance (Fig. 2b, c). Virulence factors were widespread, including elongation factor tufA, chemotaxis protein cheY, catalase katA, and heat shock protein htpB, which were present in > 50% of Desulfovibrio species (Fig. 2b). Meanwhile, some virulence factors were species-specific, including type VI secretion system contractile sheath small subunit tssB, which was prevalent in D. faecigallinarum (27/31) and also detected in D. fairfieldensis (4/117), while a type III secretion system gene (pscR) was detected in D. desulfuricans_D (1/1). Antibiotic resistome profiling also revealed several prevalent but species-specific antibiotic resistance genes (Fig. 2c). Chloramphenicol resistance gene catB7 was present in D. piger (211/212) and highly prevalent in Desulfovibrio elongatum sp. nov. (1503/1528), while beta-lactamase resistance (blaDES-1) was specific to D. desulfuricans (42/42) and detected in D. sp009712225 (1/1). These analyses indicate that Desulfovibrio species are highly diverse in their encoded virulence and antibiotic resistance genes.
To further investigate the predicted antibiotic resistance of Desulfovibrio strains, we performed in vitro antimicrobial susceptibility testing for the newly isolated Desulfovibrio isolates (Supplementary Fig. 3). Isolates from 5 Desulfovibrio species (D. desulfuricans, Desulfovibrio paradesulfuricans sp. nov. [i.e. D. desulfuricnas_A based on GTDB-tk], D. legallii, D. porci and Desulfovibrio elongatum sp. nov.) were resistant to kanamycin (minimum inhibitory concentrations [MICs] value > 64 µg/ml), while only the D. porci strain was resistant to ceftazidime. Overall, the antibiotic resistance of different Desulfovibrio isolates was highly variable, suggesting species- and strain-level variation. Not all Desulfovibrio strains carrying putative antibiotic resistance genes were resistant to the corresponding antibiotic, where genes may not necessarily be expressed under laboratory conditions or the gene product may not be functional. Nevertheless, the detection of putative antibiotic resistance genes can still provide useful insights to identify potential antibiotic resistance reservoirs and monitor the risk of resistance emergence.
Distinct core and accessory genes of within-species clades associate with geography
To evaluate factors influencing within-species genomic diversity, we performed permutational analysis of variance based on accessory gene presence/absence for the major Desulfovibrio species ( > 50 available genomes). Desulfovibrio elongatum sp. nov., Desulfovibrio piger_A, Desulfovibrio piger, and Desulfovibrio clade I, showed highly similar core gene content with no distinct cluster based on a maximum-likelihood phylogenetic tree approach. Desulfovibrio sp900319575 differed in core and accessory genes with distinct clusters for genomes from East Asia and the Netherlands (Fig. 3a, b). This was further supported by mash distance-based phylogenetic analysis (Supplementary Fig. 4). To investigate how geography associates with a species’ functional capacity, we identified genes enriched in Dutch and East Asian D. sp900319575 (Fig. 3c, Supplementary Data 7). East Asian D. sp900319575 enriched genes (East Asia prevalence ≥ 30%, Netherlands prevalence ≤ 10%) included multidrug resistance genes acrB and mexA, transcriptional regulatory gene qseB and swarming motility regulation sensor gene rssA. The detected enrichment of multidrug resistance genes in D. sp900319575 within the East Asian population may be due to increased antibiotic usage53 leading to an increase in the prevalence of resistant lineages. Meanwhile, N,N’-diacetyllegionaminic acid synthase legI and antitoxin dinJ were enriched in D. sp900319575 within Dutch populations. Notably, legI is involved in legionaminic acid biosynthesis, which is incorporated into flagellin subunits54. Moreover, qseB (regulator of bacteria adherence) and dinJ (involved in bacterial biofilm formation) are differentially distributed within genomes from East Asia and the Netherlands, suggesting a distinct Desulfovibrio virulence gene repertoire. Lastly, several insertion sequence (IS) elements linked to bacterial adaptation and evolution were highly enriched in Dutch D. sp900319575 genomes, including ISSau1, ISAsp8, and ISAzs12 (Fig. 3c). The most prevalent IS element was ISAsp8 (56.04%), followed by ISAzs12 (42.86%) and ISSau1 (37.36%), with much lower prevalence in East Asian genomes (ISAsp8: 3.85%, ISAzs12: 5.77%, ISSau1: 7.69%). Interestingly, these IS elements were unique to D. sp900319575 genomes (only exception: ISSau1 was detected in one D. sp900546145 genome). To further explore the functional context of these IS elements, we profiled and annotated their neighboring genes (within 5 kbp). Functional annotations of neighboring genes were highly diverse, including multidrug resistance and toxin-antitoxin system. These results suggest increased genome plasticity in D. sp900319575 from East Asia population relative to Dutch populations.
Fig. 3: Associations of encoded gene content of Desulfovibrio sp900319575 with geography.
a Maximum-likelihood phylogenetic tree based on core gene alignment of Desulfovibrio sp900319575 genomes. Core genes were identified based on their presence in ≥ 99% of genomes of the species at ≥ 95% identity (Roary). Inner ring colors indicate host phenotypic characteristics. “Others” includes genomes from animals or unknown sources. Outer ring colors indicate geographic regions from which the genomes originate. Regions with <5 available genomes (Middle East, South America, Africa, Central Asia) were summarised as “Other regions”. b PCA analysis of accessory genes within Desulfovibrio sp900319575 genomes (Jaccard distance). Points represent Desulfovibrio sp900319575 genomes, colours indicate host regions. c Heatmap showing predicted genes (Prokka) that were enriched in the Netherlands (prevalence ≥ 30% in the Netherlands and prevalence ≤ 10% in East Asia) and East Asia (prevalence ≥ 30% in East Asia and prevalence ≤ 10% in the Netherlands), respectively. d Variance in accessory gene content of major Desulfovibrio species explained by disease, country, genome type (isolate vs. MAG) and cohort (or study) of origin were assessed by PERMANOVA (R: adonis test). Stars represent different significance level: *p ≤ 0.05; **p ≤ 0.01; ***p ≤ 0.001. N/A indicates that genomes are only from one category of the respective variable and excluded from the analysis.
To evaluate factors associated with within-species adaptation, we performed permutational analysis of variance for the major Desulfovibrio species ( > 50 available genomes, Fig. 3d). Country and cohort were associated with variation in accessory gene content for D. sp900319575 (country: R2 = 16.76%, p = 0.001; cohort: R2 = 24.08%, p = 0.001), Desulfovibrio elongatum sp. nov. (R2 = 4.23%, p = 0.001), D. piger_A (R2 = 8.46%, p = 0.001), D. piger (R2 = 7.70%, p = 0.001), D. fairfieldensis (R2 = 21.57%, p < 0.01) and Desulfovibrio clade I (R2 = 10.93%, p < 0.05). Furthermore, disease status showed a significant but weaker association with accessory gene content in Desulfovibrio elongatum sp. nov., D. piger_A, D. piger and D. fairfieldensis. In addition, genome type (i.e., MAG vs isolate) was weakly associated with gene content for Desulfovibrio elongatum sp. nov. (R2 = 0.098%, p < 0.05). Overall, country and cohort showed the largest effect on within-species accessory gene variation indicating species divergence.
Flagellated CD- and CRC-associated D. desulfuricans and D. legallii activate TLR5
Desulfovibrio are associated with diseases, in particular IBD and CRC1,2. We next identified high-level functional traits uniquely encoded by IBD- and CRC-associated Desulfovibrio genomes, including coenzyme and lipid transport and metabolism genes (UC, CD and CRC vs. healthy, Wilcoxon test, p < 0.05, Supplementary Fig. 5). To further investigate specific disease-associated Desulfovibrio functional traits, we analysed IBD and CRC-associated Desulfovibrio COGs (Fig. 4a, b). In CD, 78 COGs were enriched and 15 depleted compared to healthy individuals (Fisher’s exact test with Benjamini-Hochberg correction for multiple testing; adjusted p < 0.25, Fig. 4b, Supplementary Data 8). Toxin hipA (COG3550), which increases antibiotic resistance in E. coli55, was significantly enriched in UC (Fisher’s exact test, adj. p = 0.067, Supplementary Fig. 6). COGs associated with health included ribosome 30S subunit, while CD-enriched COGs were related to isoprenoid biosynthesis. Notably, urease genes ureABCEF were among the most significantly enriched Desulfovibrio operons in CD (Fig. 4b). Urease was associated with gut dysbiosis and is hypothesized to contribute to CD pathogenesis49,56 and our analyses identified a potential connection with D. desulfuricans and D. legallii. For CRC, 195 COGs showed increased prevalence in disease, including COG1111 containing putative ERCC4-related helicase MPH1, a DNA damage repair protein (Fig. 4a, Supplementary Data 9). The most significant gene in CRC (Fisher’s exact test with Benjamini-Hochberg correction for multiple testing; adjusted p = 0.00026) was a putative MPH1-like helicase gene, which was highly prevalent in CRC and present in 509 of 2658 Desulfovibrio genomes. Most (63.16%) of the identified MPH1-like genes possessed helicase core domains related to DNA repair helicases and genome stability, including PF00270 (DEAD/DEAH box helicase) and PF00271 (helicase conserved C-terminal domain). In addition, 77.78% of these genes also harbored the PF09369 domain (MrfA Zn-binding), previously reported as a DNA repair domain in Bacillus subtilis57,58. This suggests a potential association of this gene with bacterial survival under stress. CRC-depleted genes in Desulfovibrio included cell shape-determining gene mreC (COG1792) (Fisher’s exact test, adj. p = 0.00093), involved in organizing peptidoglycan synthesis and cell morphogenesis in bacteria. The depletion of mreC attenuates key virulence systems, i.e. SPI-1 Type III secretion system (T3SS) and flagella, in Salmonella59, establishing a connection between mreC and bacterial pathogenicity.
Fig. 4: Disease-associated functional trait across Desulfovibrio genomes.
a, b Volcano plots showing COGs differentially enriched in CRC (a) and Crohn’s disease (CD) (b). Each point represents a COG. X-axis indicates COG prevalence in CD (or CRC) patients minus prevalence in healthy individuals. Y-axis shows log-transformed adjusted p-values (Fisher’s exact test, Benjamini-Hochberg correction for multiple testing). Dashed horizontal lines indicate significance threshold (two-sided Fisher’s exact test, adjusted p-value < 0.25). c Electron microscopy image of D. desulfuricans (MS020) shows visible flagella. The experiment was performed once, with two slides examined and approximately 10 fields of view assessed per slide. Flagella were observed on nearly all bacterial cells, with only rare exceptions. d Flagellated D. desulfuricans (MS020) and D. legallii (MS029) activate TLR5, while non-flagellated Desulfovibrio elongatum sp. nov. (MS012) showed no TLR5 activation. Statistical significance was assessed using two-way ANOVA followed by Sidak multiple-comparisons for correction, star represent significance level: ****p < 0.0001. HEK-Blue TLR5 reporter cells were treated with bacterial culture supernatant filtrates. Six biological replicates were included, with 2 technical replicates per biological replicate. The barplot shows the individual values of each biological replicate, mean and standard deviation (SD). e Recombinant flagellin proteins confirm flagellin-specific TLR5 activation. All but one of the recombinantly expressed flagellins (Fla2DD_MS020) significantly activated TLR5 (HEK-Blue TLR5 reporter cell line, linear mixed-effects model with Sidak comparison, stars represent significance level: *p < 0.05; ****p < 0.0001). The concentrations for FlaBS and FlaST are 0.5 ng/ml and the concentrations for all recombinant Desulfovirbio flagellins are 1000 ng/ml. For each flagellin protein, 5 biological replicates were analyzed. The bar plot shows the values of each biological replicate, with error bars representing the standard error of the mean across biological replicates. Abbreviations of the flagellins on the x-axis are defined as follows: FlaBS, flagellin from Bacillus subtilis; FlaST: flagellin from Salmonella typhimurium; Fla1DD_MS021, flagellin 1 from D. desulfuricans (isolate MS021); Fla1DD_MS020, flagellin 1 from D. desulfuricans (isolate MS020); Fla2DD_MS020, flagellin 2 from D. desulfuricans (isolate MS020); Fla1DF: flagellin 1 from D. fairfieldensis (isolate MS026); Fla2DF: flagellin 2 from D. fairfieldensis (isolate MS026); Fla1DL: flagellin 1 from D. legallii (isolate MS029). f Fla2DD_MS020 flagellin induces distinct gene expression patterns in murine small intestinal organoids. PCA analysis based on bulk RNAseq of flagellin-treated organoid experiments. WT wild-type organoids, Tlr5KO = TLR5 knock-out organoids. Unstimulated organoids were used as negative controls. For each condition and each phenotype, 3 biological replicates were analyzed. g Differentially expressed genes induced by Fla2DD_MS020 include upregulation of TGF-beta signaling and downregulation of genes involved in oxygen response (Wald test implemented in R package DESeq2, Benjamini-Hochberg-adjusted p < 0.05).
Overall, 62 disease-enriched COGs were shared between CD and CRC, with 22 of these being flagellar-related, including flagellar basal body rod gene flgC, flagellar basal body-associated gene fliL and flagellin gene flgL (Fig. 4a, b). Flagellins have been linked to inflammation via their ability to stimulate innate immune receptor TLR560,61. To investigate immune system activation in connection with flagellated Desulfovibrio spp., we profiled flagellar-related genes in our clinical isolates. We identified disease-enriched flagellar-related genes in all D. desulfuricans, D. fairfieldensis and D. legallii isolates, while these were absent in D. piger isolates. Transmission electron microscopy (TEM) images confirmed visible flagella in D.desulfuricans (MS020), D. fairfieldensis (MS026) and D. legallii (MS029), while no flagellation was observed in Desulfovibrio elongatum sp. nov. (MS012), as expected (Fig. 4c, Supplementary Fig. 7). These results are consistent with a previously observed flagellar structure in a D. legallii strain28. N- and C-terminal functional domains responsible for TLR5 activation were present in flagellin genes flaB2, flaB3 and flagellar filament 33 kDa core gene in all D. desulfuricans and D. fairfieldensis and most D. legallii genomes (28/30), suggesting a central role of flaB subunits. Flagellin gene structure predictions support that the flaB subunits of D. desulfuricans isolates harbor orthologs of D0 and D1 domains of Salmonella-derived flagellin (FliC), a known TLR5 activator60 (Suppl. Fig. 8). We validated TLR5 activation by D. desulfuricans and D. legallii (Fig. 4d). For this, HEK-Blue TLR5 reporter cells were stimulated with bacterial culture supernatant filtrates, experimentally confirming TLR5 activation by flagellated Desulfovibrio. Overall, these results link flagellated Desulfovibrio to immune system activation in CD and CRC.
D. desulfuricans flagellin treated organoids suppress TGF-beta signaling
To further investigate the immune response induced by these flagellins, we recombinantly expressed and purified flaB2, flaB3 and flagellar filament 33 kDa core protein genes from D. desulfuricans (Fla1DD_MS020, Fla2DD_MS020, Fla2DD_MS021), D. legallii (Fla1DL) and D. fairfieldensis (Fla1DF and Fla2DF). All but one flagellin (Fla2DF) significantly activated TLR5 (HEK Blue reporter cell lines, Fig. 4e). We further profiled immune system interactions involving Desulfovibrio flagellins using murine small intestinal organoids (TLR5 knocked-out vs wild type organoids), which represent a well-established model for the investigation of epithelial TLR5 signaling60,62,63. Image-based analysis of all treated organoids confirmed no significant differences in growth between control flagellins (FlaBS and FlaST from Bacillus subtilis and Salmonella typhimurium, respectively) and Desulfovibrio-derived proteins (Suppl. Fig. 9). This indicates that these proteins do not induce acute toxic effects on epithelial cells and transcriptional differences were not due to skewed differentiation or growth of the organoids.
We subsequently exposed wild-type and TLR5 knockout mouse organoids to the two Desulfovibrio flagellins from the same isolate (Fla1DD_MS020 and Fla2DD_MS020), including FlaBS as a positive control, and performed bulk RNA-seq. Fla2DD_MS020–treated organoids exhibited distinct gene expression profiles compared to other flagellins (Fig. 4f). As expected, FlaBS induced an inflammatory response characterised by increased expression of TNF, NF-κB, and antimicrobial peptide genes (Supplementary Fig. 10)62,63. While both Fla1DD_MS020 and Fla2DD_MS020 lacked the same canonical TLR5 signature, Fla1DD_MS020 stimulation resulted in reduced expression of several genes via TLR5 signaling, including Trib3 (a regulator of stress and inflammatory signaling), and multiple solute carrier genes involved in cellular transport. Interestingly, Fla2DD_MS020 induced a robust TLR5-independent transcriptional signature with a significant downregulation of the TGF-beta signaling pathway, including Tgfbr1, Smad2/5, and Sox6/9 genes, while oxygen response pathways were upregulated, involving Nlrp6, Egr1 and Egln3 expression (Fig. 4g). No compensatory activation of other TLRs was detected (Supplementary Fig. 11). While flagellins can also activate intracellular receptors of the nucleotide binding oligomerization domain [NOD]-like receptors (NLR) family, only Nlrp6 showed strong upregulation by Fla2DD_MS020 (Wald test; wild type organoids: Fla2DD_MS020 vs. control, adjusted p = 5.17674e-43; TLR5 knock-out organoids: Fla2DD_MS020 vs. control, adjusted p = 6.69439e-47. Supplementary Fig. 12). NLRP6 is an intracellular inflammasome activator that is upregulated in response to TNF as well as microbial and viral triggers64, however, no expression changes of IL-1β, IL-18, and ASC/Pycard were observed, suggesting that increased Nlrp6 transcription did not lead to inflammasome activation. Overall, these results demonstrate that flagellated Desulfovibrio are active modulators of the innate immune response with distinct immune signatures compared to canonical flagellins from B. subtilis and S. typhimurium.
Comparison of sulfur metabolism capabilities across a wide range of Desulfovibrio spp
Sulfur metabolism is one of the key pathways in Desulfovibrio that has been hypothesised to play a role in disease. We next compared sulfur metabolism capabilities across Desulfovibrio genomes. For this, we curated a database of sulfur metabolism genes based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways65 and previous studies16,25,66,67. There are five major pathways involved in H2S production: dissimilatory and assimilatory sulfate reduction, cysteine degradation, tetrathionate metabolism and taurine metabolism (Supplementary Data 10). The majority of Desulfovibrio genomes (81.34%) encoded the complete dissimilatory sulfate reduction pathway (hmmscan, gene hits with e-value ≤ 1E-10). Further, > 99% of the genomes contained ≥ 80% of the pathway genes, where an incomplete pathway was associated with lower genome completeness (Supplementary Fig. 13). Overall, these findings confirm the established role of dissimilatory sulfate reduction as a conserved core metabolic feature across a wide range of Desulfovibrio species.
Furthermore, we identified several taurine metabolism gene clusters. One of these clusters was associated with taurine metabolism to isethionate, which is subsequently converted to sulfite (Supplementary Fig. 14) and was previously reported in D. desulfuricans DSM 642 and Oleidesulfovibrio alaskensis G20 (formerly Desulfovibrio alaskensis G20)67. Many taurine metabolism genes (e.g., putative taurine pyruvate aminotransferase [tpa] and sulfoacetaldehyde acetyltrasferase [xsc]) were detected in Desulfovibrio spp. (based on sequence similarity, Supplementary Fig. 15). However, D. desulfuricans DSM642 previously showed a lack of Tpa protein expression67, highlighting the importance of having culturable isolates available that can be tested in different culture conditions in combination with functional readouts to further investigate taurine metabolism signals in Desulfovibrio. We also identified species-specific differences in genes encoding alanine dehydrogenase (ald). All D. desulfuricans, D. legallii_A, D. intestinalis, and D. sp011039135 genomes encoded ald, while the gene was absent in D. piger, D. piger_A, Desulfovibrio clade I, Desulfovibrio clade II and clade III. These results suggest that while Desulfovibrio spp. shared core taurine metabolism pathway, variation in ald suggests species-specific difference in processing taurine-derived intermediates.
While dissimilatory sulfate reduction gene prevalence was comparable between individuals with IBD or CRC and healthy individuals (Fisher’s exact test, adjusted p > 0.25), other sulfur metabolism genes were enriched in disease, in particular ald gene for supporting taurine metabolism and ttrBCA genes responsible for tetrathionate respiration (enriched in CRC and CD, Fisher’s exact test, adjusted p < 0.25, Supplementary Fig. 16). Conversely, yhaO, involved in cysteine degradation, showed reduced prevalence in CD and CRC compared to healthy individuals. Overall, this suggests that other H2S pathways are encoded by IBD- and CRC-associated Desulfovibrio spp.
Investigating H2S production capacity of the IBD microbiome
Next, we investigated the overall IBD microbiome capacity for H2S production. We searched for sulfur metabolism genes in all MAGs in a deeply sequenced IBD cohort68 (hmmscan, gene hits with e-value ≤ 1E-10). Dissimilatory sulfate reduction was encoded by 40 of 54 (74.07%) Desulfovibrio MAGs spanning 7 known and one uncharacterized species (Supplementary Data 11, 12). In addition, 6 Mailhella (37.5%) MAGs encoded 80% of the pathway genes, where all Mailhella genomes lacked sulfate adenylyltransferase (sat/met3). The tetrathionate pathway was more widely distributed and identified in 41 MAGs from 9 genera: Desulfovibrio, Gordonibacter, Citrobacter, Proteus, Eggerthella, Escherichia, Rubneribacter, Morganella and Turicimonas (Fig. 5a, b). Overall, 443 MAGs from 16 genera, including Escherichia and Citrobacter, encoded seven genes involved in assimilatory sulfate reduction (cysC, cysD, cysH, cysI, cysJ, cysN, cysNC) and only lacked sulfate adenylyltransferase (sat) and sulfite reductase (ferredoxin) (sir) (Fig. 5a, b). As cysD and cysN together perform the same function as sat69,70, and cysIJ are alternatives for sir71,72, these bacteria still maintain the genomic capacity required for sulfide production through the assimilatory sulfate reduction pathway. While sulfide is typically directly incorporated into cysteine biosynthesis73,74 and does not result in free H2S75, cysteine can be further degraded, leading to an indirect contribution of free H2S. Overall, our genomic analysis revealed several bacteria that encode the complete tetrathionate metabolism, assimilatory and dissimilatory sulfate reduction pathway, identifying these pathways as potential H2S sources (including relevant precursors) in this IBD cohort.
Fig. 5: H2S pathway-encoding bacteria in IBD.
a Three pathways involved in H2S production (dissimilatory sulfate reduction, assimilatory sulfate reduction and tetrathionate metabolism) were encoded by different genera (FedericiS_2022 cohort68). Only genera, where 100% of the necessary pathway genes were present in at least one genome, are shown. Assimilatory sulfate reduction: cysD and cysN are functionally equivalent, sat and cysIJ are alternatives for sir, therefore strains lacking sat or sir but encoding cysD, cysN, and cysIJ were still classified as encoding the complete assimilatory sulfate reduction pathway. b The heatmap shows the presence of genes related to these three pathways involved in H2S production in species with > 3 MAGs encoding at least one complete pathway. *sat is involved in dissimilatory and assimilatory sulfate reduction according to KEGG modules M00596 and M00176. c Genes involved in tetrathionate metabolism were identified at community level with significantly increased abundance in UC (n = 188) and CD patients (n = 120) compared to controls (n = 229). Wilcoxon test: *p < 0.05; **p < 0.01; ***p < 0.001. Box plots show the distribution of values, with the centre line indicating the median, the lower and upper bounds of the box representing the 25th and 75th percentiles, respectively, and the box representing the interquartile range (IQR). Whiskers extend to the smallest and largest values within 1.5 × IQR of the lower and upper quartiles, respectively; values beyond the whiskers are shown as outliers. d Genomic context of the ttrBCA operon in the Proteus mirabilis MAGs. COG pathway information is indicated by different colors.
Due to the lack of a universally conserved gene sets and the presence of many alternative enzymes in bacterial taurine metabolism and cysteine degradation, we did not assess completeness of these pathways. Instead, we explored these pathways at the microbial community level to identify their associations with IBD. Several taurine metabolism genes were significantly increased in IBD patients compared to healthy individuals (adj. p < 0.05, Supplementary Fig. 17), including enzymes that catalyse the metabolism of taurine to isethionate with subsequent conversion to sulfite (tpa, sarD, islA and islB). In contrast, phosphotransacetylase (pta), which converts acetyl-phosphate to acetate, a short-chain fatty acid important for human epithelial integrity76, showed reduced abundance in IBD (adj. p < 0.05, Supplementary Fig. 17). Moreover, taurine dehydrogenase large subunit tauY abundances were increased in UC but decreased in CD compared to healthy individuals. These results indicated an association of changes in taurine metabolism gene abundances in IBD connected to acetate production. Similarly, most genes related to cysteine degradation are more abundant in IBD patients compared to healthy individuals (adj. p < 0.05, Supplementary Fig. 18). In particular, the major cysteine-degrading genes (dcyD, yhaM and sseA) were increased in IBD samples, where sseA has been found to be associated with bacterial virulence77. These results indicate an increase in bacteria encoding H2S-associated taurine and cysteine degradation pathways in IBD patients, which is in line with previously reported associations between organic sulfur-metabolizing bacteria and CRC risk66. Collectively, this supports the hypothesis that organic sulfur metabolism may represent an important H₂S source in the disease-associated microbiome.
The majority of MAGs (86.49%) encoding tetrathionate metabolism were recovered from IBD patients. Interestingly, these MAGs encoded the ttr genes, responsible for tetrathionate respiration, which has been shown to confer growth advantages in pathogenic bacteria during inflammation78. Indeed, the abundance of all genes involved in this pathway was significantly higher in IBD compared to healthy individuals (exception: phsC in UC vs. healthy, p = 0.074, Wilcoxon, all other adj. p < 0.05, Fig. 5c). To further investigate gene activity related to tetrathionate metabolism, we analysed IBD metatranscriptomic data79. In CD patients, all but one tetrathionate metabolism gene showed increased expression compared to non-IBD controls (exception: phsB in CD vs. non-IBD, p = 0.46, Wilcoxon, all other adj. p < 0.05, Supplementary Fig. 19). For UC patients, phsA, phsC and ttrB were expressed in most of the samples, while only ttrA showed a significant increase in UC vs. non-IBD, which may be due to smaller UC sample numbers. MAGs encoding the complete tetrathionate metabolism pathway included 5 Proteus MAGs from 3 species (P. mirabilis, P. vulgaris, P. sp003144505). Proteus mirabilis, a known CD-associated pathogen80, contained ttrRS and ttrBCA gene clusters (Fig. 5d), which can confer competitive growth advantages80. Overall, these results indicate that tetrathionate metabolism abundance is increased and transcriptionally active in CD and identifies other potential new players encoding disease-associated H2S production pathways.