Freshly emitted fecal samples of healthy individuals were collected as described in our recent study (14). After sample homogenizing and 10-folds serial dilutions, microbial cells in supernatant were plated onto brain heart infusion agar (BHI Agar, BD Difco/BBL, USA) supplemented with 5% defibrinated sheep blood and incubated at 37°C under anaerobic condition. Isolated colonies were subjected to 16S rRNA gene sequencing and the resultant sequences were identified via the BLAST program of the EzBioCloud databases. Based on the species demarcation threshold of 98.7% 16S rRNA gene sequence similarity (16), the isolated strain B. vulgatus 46 (CGMCC NO.17140) was sub-cultured, harvested and suspended in phosphate buffered solution (PBS) to a concentration of 5×10⁹ colony forming units (CFU)/mL. The B. vulgatus Bv46 (0.2 mL/mouse) was administered daily by gavage within 1 h after harvest. OD₆₀₀ measurements and viable counts were used to confirm the concentration of bacteria.

The genomic DNA of B. vulgatus Bv46 extracted using the Wizard Genomic DNA Purification kit (Promega) was sequenced on the Illumina HiSeq TM2000 platform, and the draft genome was assembled as described previously (17). Genome sequence of fourteen other B. vulgatus strains (12) were retrieved from NCBI. All the genomes were annotated using Prokka v1.14.6 (18).

To investigate the phylogenetic relationships and evaluate the genetic distance among B.vulgatus strains, orthoFinder v2.5.4 was used to cluster orthologous genes (19), and MAFFT v7.453 was used for constructing multi-genome alignments (20), a neighbor-joining tree was then constructed with the maximum composite likelihood method with GAMMA distributed correction on rates among sites using MEGA X (21). The alignment gap or missing data were treated as pairwise deletion. The branches were evaluated using 1,000 times bootstrap resampling. Pan-genome computation of B.vulgatus strains Bv46, ATCC 8482, mpk and FTJS7K1 was performed using Roary (22), unique genes of B.vulgatus Bv46 were assigned a functional Clusters of Orthologous Group (COG). To predict the carbohydrate unitization, carbohydrate-active enzymes (CAZymes) were identified by searching for matches in the CAZy database using the standalone tool dbCAN v3.0.7 (23).

Specific-pathogen free (SPF) female C57BL/6J mice (6 weeks old) weighting 16-18 g were obtained from Beijing Vital River Laboratory Animal Technology. After a one-week acclimation, the mice were randomly divided into three groups (n = 8/group): control, PBS-treated colitis (DSS group) and B. vulgatus Bv46-treated colitis group (BV group). Both the DSS and BV groups were fed 3% DSS (36-50kDa, S0798, MP Biomedicals Canada) for 7 consecutive days to induce colitis. Meanwhile, B. vulgatus Bv46 (1×10⁹ CFU/mouse) and PBS were given orally in the BV and DSS groups, respectively. Control group was given plain water. The experiment mice were monitored daily for body weight, stool consistency and rectal bleeding, and the disease activity index (DAI) was calculated as previously described (24). At the end of the experiment, fresh feces were collected, frozen in liquid nitrogen, and then stored at -80 °C refrigerator. After the mice were euthanized, the colons were collected for histopathological observation, and the cecal contents were collected for microbiota analyses. The colon lengths from the ileocecal junction to the anus were measured (25).

The distal colon segment of each mouse was fixed in 4% paraformaldehyde (PFA) (BL539A, biosharp, China) for histological analysis. Cross sections were stained with hematoxylin-eosin (H&E) and the pathological features were analyzed by two qualified histopathologists without knowledge of classification. Histopathological analysis was performed according to the colitis scoring system (2, 9) with scores ranging from 0 (no pathology) to 4 (highest pathology). The average percentage of field was calculated at a given score and compared between treatment groups.

Bacterial genomic DNA was extracted from mouse cecal contents using MagPure Stool DNA KF kit (Magen, China), and the V3-V4 hypervariable regions of prokaryotic 16S rDNA were magnified with universal primers 341F (5’ -ACTCCTACGGGAGGCAGCAG-3’) and 806R (5’-GGACTACHVGGGTWTCTAAT-3’). After quality assurance, equimolar amounts of validated amplicon libraries were sequenced on Illumina HiSeq 2500 platform by BGI Tech Solutions Co., Ltd (Beijing, China) following the manufacturer’s instructions. The resulting paired-end reads were analyzed using scripts of EasyAmplicon v1.14 (26). For β-diversity analysis, Bray-Curtis dissimilarity matrice was measured and used in a non-metric multidimensional scaling (NMDS) analysis.

To determine potential microbial biomarkers in three groups, Linear Discriminant Analysis Effect Size (LEfSe) was calculated using the online tool (27). For the quantification of B. vulgatus in feces, species specific real-time PCR was performed as previously described (28, 29). Microbiota sequencing data used in this study are available at https://www.ncbi.nlm.nih.gov/bioproject/PRJNA872285.

Total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA) according to the manufacturer’s instructions and genomic DNA was removed using DNase I (Qiagen, Valencia, CA). The RNA size, integrity, and total amount were measured using a Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA).

As described previously (30–33), RNA sequencing (RNA-seq) transcriptome libraries were prepared using TruSeqTM RNA sample preparation Kit (Illumina, San Diego, CA) according to the manufacturer’s protocol and sequenced on the Illumina sequencing platform (HiSeq 4000). RNA-Seq was performed by BGI Genomics Co., Ltd.

The raw paired end reads were trimmed and quality controlled by Trimmomatic with default parameters (34). Clean reads were then subjected to Salmon (35) in mapping-based mode for transcript quantification, and the current version of Mus_musculus (http://ftp.ensembl.org/pub/release-106/) was served as a reference. Differentially expressed genes (DEGs) between the BV and DSS groups were subsequently identified using the DESeq2 package (36), the thresholds were set as: adjusted P value < 0.05, log2(fold change) > 1.0 or log2(fold change) < -1.0. To further understand the biological roles of colonic DEGs, Gene Ontology (GO) enrichment analysis was conducted using clusterProfiler package (37), the enriched maps of GO biological processes were visualized using Cytoscape v3.7.2. For Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, the gene set enrichment analysis (GSEA) (38) was applied to pre-ranked colon genes based on the log2 transformed fold change in expression, gene sets with the false discovery rate (FDR) < 0.01 and the absolute value of normalized enrichment score (|NES|) > 1.0 were defined as statistically significant. RNA-seq data are publicly available at NCBI under accession number PRJNA872866.

Reverse transcription quantitative real-time PCR (qRT-PCR) was performed to validate five selected DEGs (Cd40, Cd19, Cd22, Ccl19 and Cxcr5) identified from transcriptome sequencing. qRT-PCR primers were listed in Table S1 in the Supplemental Material .

Total RNA extraction was performed using TRIzol reagent (Invitrogen, Carlsbad, CA) following the manufacturer’s instructions, and cDNA synthesis was conducted using the PrimeScript™ RT Reagent Kit (Perfect Real Time, TaKaRa). qRT-PCR was carried out using SYBR Premix Ex Taq II (Perfect Real Time; TaKaRa) in a Rotor-Gene Q thermal cycler (Qiagen, Valencia, CA). Data were analyzed with Rotor-Gene Q series software version 1.7 (Qiagen, Valencia, CA). The data were normalized to the endogenous reference gene β-actin and analyzed by the cycle threshold method (2^-ΔΔCt) (39). Three independent replicates were carried out for each target.

Colon samples weighing 0.1g were homogenized in PBS and then centrifuged at 12000 g (10 min at 4°C). The amount of tumor necrosis factor (TNF)-α, IL-6 and IL-1β in the supernatant of the colon homogenates were determined by enzyme-linked immunosorbent assay (ELISA).

The anti-inflammatory effect of B. vulgatus Bv46 on the lipopolysaccharide (LPS, L2880, Sigma-Aldrich) activated RAW 264_·7 macrophage cell line was tested as described previously (40). After 48 h of anaerobic culture in BHI broth at 37°C, B. vulgatus Bv46 was harvested and suspended in PBS at a concentration of 1×10⁸ CFU/mL. Then the B. vulgatus Bv46 live bacteria were added to RAW 264_·7 cells (multiplicity of infection, MOI = 20) for 6 or 20 h in the presence of 1 μg/mL LPS. At the same time, cells treated with 1 μg/mL LPS were used as positive control, and the non-treated cells as negative control. TNF-α, IL-1β and IL-6 levels of the cell supernatants were determined by ELISA.

After 48 h anaerobic incubation in BHI broth at 37°C, 0.1 mL supernatant of pure culture B. vulgatus Bv46 was collected, and then added to a mixture containing 0.05 mL of 50% H₂SO₄ and 0.2 mL of methylvaleric acid. The solution was swirled for 1 min and ultrasonically treated for 10 min (incubated in ice water). After centrifugation for 15 min at 9600 g, 4°C, the supernatant was transferred into a fresh vial for gas-chromatography tandem mass spectrometry (GC-MS) analysis.

For fecal samples, 150 mg feces was added to 1 mL H₂O, homogenized in a ball mill for 4 min at 45 Hz, and then ultrasonically treated for 5 mins (incubated in ice water). After centrifugation for 20 min at 5000 rpm, 4°C, 0.8 mL of the supernatant was collected and transfered to the mixture containing 0.1 mL of 50% H₂SO₄ and 0.8 mL of methylvaleric acid. Subsequently, the solution was vortexed (10 s), sonicated (10 min, ice water bath), and centrifuged (12,000 g, 15 min, 4°C), then the obtained supernatant was transferred to a fresh glass vial for GC-MS analysis. A Shimadzu GC2030-QP2020 NX gas chromatography-mass spectrometer (Kyoto, Japan) fitted with a HP-FFAP capillary column (30m × 250 μm × 0.25 μm, J&W Scientific, Folsom, CA, USA) was used for short chain fatty acids (SCFAs) quantification as described previously (41).

Statistical analysis was conducted using GraphPad Prism 9.0. Data presented as mean ± standard deviation (SD). The unpaired Student’s t-test was employed to compare two groups. One-way analysis of variance (ANOVA) followed by the Tukey multiple comparison test was applied to datasets comparing more than two groups. P < 0.05 was considered statistically significant.

To investigate the genomic differences between B. vulgatus Bv46 and other B. vulgatus strains, we conducted a comparative genome analysis. A neighbor-joining phylogenetic tree based on 1796 core genes was constructed, which showed that B. vulgatus Bv46 was placed in its own lineage that was sister to strain FJSWX62K35, and separated from B. vulgatus strains ATCC 8482, mpk and FTJS7K1, all of which have been reported to attenuate colitis ( Figure S1A ). Comparison of B. vulgatus strain Bv46 with ATCC 8482, mpk and FTJS7K1 showed that B. vulgatus Bv46 had 872 unique genes. Functional annotation assigned the B. vulgatus Bv46 strain-specific genes to 21 COG classes, with the majority in replication, recombination and repair (L, 44 genes), cell wall/membrane/envelope biogenesis (M, 40 genes), transcription (K, 23 genes), and carbohydrate transport and metabolism (G, 19 genes) apart from genes with unknown function ( Figures S1B, Table S2 ). Further carbohydrate-active enzymes analysis predicted that the unique genes of B.vulgatus Bv46 contained 25 glycosyltransferases (GTs) genes, 12 glycoside hydrolases (GHs) genes, 1 carbohydrate-binding module gene and 1 polysaccharide lyase gene ( Table S3 ).

To test whether B. vulgatus Bv46 is protective against DSS-induced colitis in mice, we determined the severity of colitis by monitoring the DAI index during the 7-day induction period and the colonic length at the end of experiment. As shown in Figures 1 and S2 , the DSS group displayed notable weight loss (days 4-7 of DSS exposure; P < 0.01 versus control group), high DAI score (days 4-6 of DSS exposure; P < 0.01 versus control group) and shortened colon length (at day 7 of DSS exposure; P < 0.01 versus control group). Comparing with the DSS group, the DAI score was significantly reduced in the BV group (P < 0.05 at days 4-6; Figure 1A ), and the colon length was significantly longer in the BV group at day 7 of DSS exposure (P < 0.01; Figures 1B, C ). However, the weight losses were comparable between the DSS and BV groups (P > 0.05; Figure S2 ).

Microscopic examination of colon histopathological sections showed that DSS-induced histopathological damage was generally moderate to severe, and mucosal lesions and submucosal inflammatory cell infiltration were observed in some visual fields. In the DSS group, 23.3% visual fields showed moderate scores (grade 2), 75% visual fields showed high scores (grade 4) and only 1.7% visual fields with no pathology (grade 0) ( Figure 2A ). In the control group, 100% visual fields had no pathology (grade 0). There were significant differences in pathological grades 0, 2 and 4 between the control and DSS group (P < 0.01). In the BV group, 36.7% visual fields had no pathology (grade 0), 28.3% visual fields had moderate pathology (grade 2) and 33.3% visual fields had high pathology (grade 4). Compared with the DSS group, histopathological scores in grades 0 and 4 were significantly reduced (P < 0.01) in the BV group ( Figure 2A ). Inflammatory cell infiltration, mucin depletion and epithelial erosion were all less evident in the BV group ( Figure 2B ).

The V3-V4 region of the 16S rDNA was sequenced to assess the modulatory effects of B. vulgatus Bv46 on gut microbiota of the DSS-induced mice. The rarefaction curves of individual groups indicated that the sequencing depth was adequate to reliably describe the bacterial microbiome associated with these three groups ( Figure S3A ). NMDS analysis showed that the gut microbiota profiles (β diversity) were grouped according to diet ( Figure 3A ). The top abundant taxa at the phylum and genus levels were shown in Figures 3B and S3B . Similar to previous studies (42), at the phylum level, Firmicutes and Bacteroidetes were the predominant constituents in mouse gut microbiota ( Figure S3B ). Figure 3B showed the top 15 abundant genera among the control, DSS and BV groups. In the control group, Ligilactobacillus had the highest relative abundance. However, in the DSS group, the relative abundance of Ligilactobacillus was obviously lower, and Kineothrix was the dominant genus. In the BV group, Phocaeicola emerged as a dominant genus.

To identify bacterial taxa that have changed significantly among the control, DSS and BV groups, colonic microbial composition was analyzed the LEfSe based on nonparametric factorial Kruskal-Wallis and Wilcoxon tests. Figures 3C, D showed the taxa with significant differences, indicated by an LDA score greater than 3.0, which reflected the degree of influence of a taxa with a significant difference among the three groups. At the phylum and class levels, Bacteroidetes and Bacteroidia were enriched in the BV group. Moreover, B. vulgatus Bv46 treatment significantly increased the abundance of Parabacteroides, Bacteroides, Anaerotignum and Alistipes at the genus level. For the DSS group, genera Kineothrix, Agathobaculum, Faecalimonas and Flintibacter were enriched. However, genus Ligilactobacillus was enriched in the control group. With respect to B. vulgatus, a higher biomass was observed in the feces of mice fed with B. vulgatus Bv46, compared with that of the DSS group ( Figure S3C ).

Colonic RNA-seq was performed to analyze the differential expression of colonic genes between the DSS and BV groups. There were 25 upregulated and 83 downregulated genes in the BV group compared to those in the DSS group ( Figure 4A ). Through GO-BP (Biological Process) enrichment analysis, 19 downregulated DEGs in the BV group were significantly enriched in the regulatory B cell response, such as B cell activation, proliferation and differentiation pathway, B cell receptor signaling pathway and regulation of neutrophil chemotaxis pathway; two upregulated DEGs (Hba-a1 and Hbb-bt) in the BV group were significantly enriched in hydrogen peroxide catabolic process ( Figure 4B , Table S4 ). According to the BP pathway analysis, the downregulated genes, Cd40, Cd79a, Cd19, Cd22, Siglecg, Bst1 and Tnfrsf13c participated in over six pathways ( Figure 4C , Table S5 ).

Further GSEA analysis of mice colon transcripts found that 12 upregulated and 34 downregulated KEGG pathway terms were enriched in the BV group ( Table S6 ). Among 34 downregulated KEGG pathway terms, 11 were immune-related signaling pathways, including T cell receptor signaling pathway, Toll like receptor signaling pathway, cell adhesion molecules, B cell receptor signaling pathway and cytokine-cytokine receptor interaction ( Figure 4D ). Five of 19 downregulated DEGs in the BV group, including Cd40, Cd19, Cd22, Ccl19 and Cxcr5 also partake in two or more immune-related GSEA KEGG signaling pathways ( Table S5 ). To further validate the RNA-Seq data, these five genes were confirmed by qRT-PCR, which exhibited a concordant direction in both RNA-Seq and qRT-PCR ( Figure 5A ).

The expression of colonic inflammatory cytokines TNF-α, IL-6 and IL-1β in the DSS and BV groups were determined. The expression of these three colonic proinflammatory cytokines in the BV group were lower than those in the DSS group (P < 0.05; Figure S4 ).

The macrophage RAW 264.7 was used to evaluate the potential anti-inflammatory effect of B. vulgatus Bv46 (40, 43, 44). As showed in Figure 5B , B. vulgatus Bv46 markedly reduced the secretion of TNF-α, IL-1β and IL-6 in LPS-stimulated RAW 264.7 cells (P < 0.05).

Five SCFAs, including acetic acid, propionic acid, butyric acid, isobutyric acid and isovaleric acid were detected in the culture supernatant of B. vulgatus Bv46 with concentrations of 1486.21, 173.75, 0.95, 6.07 and 8.15 μg/ml, respectively.

As shown in Figure 6 , the fecal level of acetic acid, butyric acid, propionic acid, isobutyric acid, valeric acid and isovaleric acid were determined for all groups. Excepting for butyric acid, there was no significant difference in fecal SCFAs between the CT and DSS groups. The BV group showed significantly higher concentration of five of the six fecal SCFAs (except for acetic acid) compared with the DSS group (P < 0.05), indicating that B. vulgatus Bv46 supplementation contributed to the increase of these SCFAs.