All experiments were conducted in strict accordance with the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health. The animal experiments were reviewed and approved by the Institutional Animal Care and Use Committee of Sun Yat-sen University (Permit No: 2016-104) and the Medical Research Ethics Committee of Sun Yat-sen University (SYSU-IACUC-2019-B517).

B. subtilis CMCC(B) 63501 was purchased from Solarbio Science and Technology Co., Ltd (Beijing). B. subtilis was cultured in a lysogeny broth (LB) medium to spawn for 24 h. Then, bacteria were successively washed with 1 M NaCl and 1 M KCl and washed two times with distilled water. The LB agar plates were incubated at 27°C for 12–24 h, and the number of CFUs per plate was counted. Finally, the concentration of bacterial suspension was adjusted to 3 × 10⁸ CFU/ ml, 3 × 10⁹ CFU/ ml, or 3 × 10¹⁰ CFU/ ml. Each animal was given 0.3 ml of the final suspension every 3 days. The bacteria suspension was refreshed every week.

A total of 36 male pathogen-free BALB/c mice, approximately 6 weeks old (body weight: 18 ± 2 g), were purchased from the Experimental Animal Center of Southern Medical University. They were reared in plastic cages with free access to autoclaved chow and water in the Biosafety Level-2 (BSL-2) laboratory of Sun Yat-sen University under controlled temperature and humidity and a 12-h light and 12-h dark cycle. The animals were randomly divided into groups for further experiments. Oncomelania hupehensis was purchased from the Chinese Center for Disease Control and Prevention (Shanghai). After acclimating to the laboratory environment, each mouse was infected with 20 ± 2 S. japonicum cercariae via shaved abdominal skin. The probiotic supplement for mice was performed 1 week before the first day of infection, and treatment continued for 7 weeks. To conduct experiments, we randomly divided mice into four groups: normal group (NG), normal mice with OA B. subtilis group (NBS), S. japonicum–infected group (SI), and S. japonicum–infected mice with OA B. subtilis group (SIBS).

The mice were sacrificed after chloral hydrate asphyxiation and cervical dislocation of S. japonicum at 56 days postinfection (dpi). Left liver lobes and colons were collected and immediately fixed in 4% paraformaldehyde for histopathological analysis. Liver samples were also collected for hydroxyproline content measurement according to the protocol of the hydroxyproline assay kit (Nanjing, China). Stool samples were collected the day before sacrifice. Blood samples were drawn from orbital veins and centrifuged at 1,500 × g for 15 min, and then the serum was collected after clotting. In addition, both male and female S. japonicum worms were collected from the portal vein. The detections of worm length, worm burden, and egg burden were determined as described in a previous study (Shen et al., 2017).

For histopathological analysis, fixed fragments from the intestine and liver of mice were sliced into sections (5 μm thick). These slices were subjected to H&E staining and Masson’s trichrome staining. Images were captured under an inverted microscope (Olympus, Japan). The percentage of the fibrotic area was detected using a ZEISS Axio Scan. Z1 automated slide scanner microscope (Germany). A full view of the whole tissue was also obtained. We analyzed the whole tissue area and the blue-positive region using Image-Pro Plus 6.0 software (Media Cybernetics, USA). We calculated the percentage of the fibrotic area based on the area of the blue-labeled region/the total area of the whole liver tissue. Granulomatous responses in the colon were estimated by calculating the ratio of the integrated optical density (IOD) to the intestinal tissue area.

The liver and intestinal tissues were collected and stored in TRIzol reagent (Invitrogen, USA) at −80°C until processing. Total RNA was extracted as described in a previous study (Lin et al., 2020). We quantified the total RNA using a NanoDrop 2000 spectrophotometer (Thermo Scientific, America).

RNA integrity was assessed using the RNA Nano 6000 Assay Kit protocol in the Agilent Bioanalyzer 2100 system (USA). Total RNA was used as input material for the RNA sample preparations. PCR was performed with Phusion high-fidelity DNA polymerase, universal PCR primers, and index (X) primers. The PCR products were purified. The library quality was assessed on a Qubit2.0 Fluorometer and Agilent Bioanalyzer 2100 system. The prepared library was sequenced on an Illumina NovaSeq platform by Frasergen Company (Wuhan, China) and 150 bp paired-end reads were generated. After quality control of raw data, we carried out transcriptome assembly, gene functional annotation, differential expression analysis, KEGG enrichment analysis, and GO enrichment analysis. Genes with an adjusted fold change ≥2 were assigned as DEGs using the DESeq2 R package (1.20.0). The cluster profiler R package was used to test the statistical enrichment of DEGs in KEGG pathways.

All stool samples were frozen at −80°C for further study. Total DNA from stool samples was isolated under a sterile environment according to the Hipure Stool DNA Kit protocol (Magen, China). After extraction, the total DNA quality and quantity examination were conducted using a NanoDrop 2000 spectrophotometer (Thermo Scientific, America).

The V3–V4 region of the 16 S rRNA gene (approximately 500 bp) was amplified using a specific bacterial primer set (forward primer 5′-ACT​CCT​ACG​GGA​GGC​AGC​A-3′ and reverse primer 5′-GGACTACHVGGGTWTCTAAT-3′) and sequenced on an Illumina HiSeq 2500 platform. PCR amplification was performed using Takara PrimeStar DNA polymerase (China). The following PCR cycling conditions were used: denaturation at 95°C for 5 min, 25 cycles of 95°C for 30 s, 50°C for 30 s, 72°C for 40 s, and final extension at 72°C for 5 min. The PCR products were analyzed by 2% agarose gel electrophoresis. Finally, sequencing was performed by Biomarker Technologies (China).

After the base calling analysis, the raw data files from the sequencing platform were transformed into the original sequenced reads and stored in FASTQ format. QIIME (version 1.8.0) was used to cluster reads into operational taxonomic units (OTUs) and identified at 97% or more similarity (Caporaso et al., 2010). We rarified the OTU table and calculated the species abundance based on the OTUs and ACE indices using mothur (version 1.33.3) (Grice et al., 2009). To analyze the alpha diversity. As measures of beta diversity of the similarity, the nonmetric multidimensional scaling (NMDS) and partial least squares discrimination analysis (PLS-DA) plots with binary Jaccard distance were performed in R with the vegan package (Kambura et al., 2016). PERMANOVA was used to evaluate the beta diversity between samples or groups using mothur (version 1.33.3). We computed and explored the taxonomic content of the sequencing dataset using MEGAN (Huson et al., 2007). A ternary plot of the microbial community was generated as described in a previous study (Han et al., 2018). Based on the genus abundance of the gut microbiota of mice, random forest algorithms with 1,000 random permutations were conducted (Li et al., 2021). For range adjustment, all pairwise comparisons between two groups were tested using Student’s t-test.

We calculated the results using GraphPad Prism version 6.0 (USA). Data are expressed as the mean ± SEM. The differences between groups were analyzed by Student’s t-test using SPSS 19.0 software (USA). ^* p < 0.05, ^** p < 0.01, and ^*** p < 0.001 were considered statistically significant.

To examine the role of B. subtilis in schistosomiasis-related hepatic fibrosis in vivo, we infected the mice with a lethal dose of S. japonicum cercariae. The infected mice were gavaged with B. subtilis or PBS every 3 days. Experimental designs on schedule for parasite infection, OA B. subtilis or PBS, and sample withdrawal are shown (Figure 1A). The survival rate of mice is shown (Supplementary Figure S1A).

Hepatic pathological injuries and fibrosis are severe symptoms during schistosomiasis and may cause death in hosts. Therefore, we next investigated whether OA B. subtilis (3 × 10⁸ CFU/ ml) protection involved attenuation of hepatic injuries and fibrosis. Based on the results of Masson’s trichrome staining (Figures 1B–D) and hydroxyproline quantification (Figure 1E), we found that infected mice displayed a significant reduction in the total area of hepatic granulomas and fibrosis after OA B. subtilis. In addition, the worm burden, egg count in the liver, and intestine and worm length from the mice were not significantly different between the groups (Supplementary Figures S1B–D).

However, whether the dead B. subtilis or living B. subtilis plays a role in attenuating pathological injuries in S. japonicum–infected mice is unclear. Subsequently, the infected mice were administered living B. subtilis (SIBS) or heat-killed B. subtilis (SIDBS). We found no significant difference in the area of fibrosis between infected mice and infected mice with OA-dead B. subtilis (Supplementary Figures S2A, B), in addition to the differences in the hydroxyproline quantification (Supplementary Figures S2A, C). However, it was significantly different in infected mice with and without OA living B. subtilis (Supplementary Figures S2). These findings indicated that OA living B. subtilis attenuates hepatic pathological injuries in S. japonicum–infected mice.

In addition, schistosome infection causes a wide range of clinical symptoms including gut inflammation, mainly induced by releasing eggs trapped in the intestinal wall. We next investigated whether OA B. subtilis protection involved attenuation of intestinal injuries. We found that S. japonicum eggs were distributed in the submucosa or even involved in the mucosa of infected mice (Figure 1F). These results indicated that parasitic eggs may affect the integrity of the intestinal tissue. Importantly, we found that infected mice treated with OA B. subtilis displayed a significant reduction in the total area of intestinal granulomas, as shown by HE staining (Figures 1F,G).

Evidence suggests that S. japonicum infection contributes to the dysbiosis of gut microbiota in mice, so we next investigated whether OA B. subtilis can affect the gut microbiota of mice during schistosomiasis using 16 S rRNA gene high-throughput sequencing. We analyzed the bacterial community structure and diversity of the intestinal contents of 36 individual mice. A total of 2,813,530 valid reads and 2,515,622 clean tags were retained from 36 stool samples, with an average of 69,878 clean tags per sample after filtering (Supplementary Table S1). A total of 416 OTUs were obtained at a 97% similarity level among all samples based on the Silva database. The sequenced samples were taxonomically clustered into 9 phyla, 15 classes, 21 orders, 37 families, 97 genera, and 102 species (Supplementary Table S2).

The mouse gut microbiota harbored 21 bacterial phyla, with two dominant phyla (Firmicutes and Bacteroidetes) accounting for over 73.4% of the total relative abundance (Supplementary Figure S3A). The phylum Firmicutes was the most frequent taxon in the mouse gut at the phylum level, followed by Bacteroidetes, Proteobacteria, and Epsilonbacteraeota. We found that Bacteroidia, Clostridia, and Bacilli were the most abundant gut bacteria at the class level, with total average relative abundances of over 72.5% in the six groups (Supplementary Figure S3B). was the most prevalent abundant bacterium in all groups at the family level, followed by Muribaculaceae, Rikenellaceae, Ruminococcaceae, Prevotellaceae, and Lactobacillaceae. (Supplementary Figure S3C). At the genus level, uncultured bacteria of Muribaculaceae and Alistipes were the most prevalent gut microbiota detected in BS and NG group mice, while Alloprevotella and Lactobacillus were the most prevalent gut microbiota detected in SIBI and SI groups (Supplementary Figure S3D).

To evaluate bacterial community differences, we further analyzed the beta diversity of gut microbiota among groups. We found that the bacterial communities of normal mice with and without OA B. subtilis were similar (quantified by NMDS analysis) (Figure 2A). After S. japonicum infection, significant differences in bacterial community structures between infected mice and normal mice were found (quantified by NMDS and PLS-DA analysis, permanova: p < 0.05) (Figures 2A–C). Moreover, the community structures of the gut microbiota of S. japonicum–infected mice with OA B. subtilis were significantly different from those of infected mice (quantified by NMDS and PLS-DA analysis, permanova: p < 0.05) (Figures 2A,B,D). Our findings suggested that OA B. subtilis can significantly modulate the beta diversity of gut microbiota in S. japonicum–infected mice.

To evaluate how OA B. subtilis affects the diversity of the gut microbiota of S. japonicum–infected mice, we next analyzed the diversity of the gut microbiota in the intestinal contents. We found that S. japonicum–infected mice showed significantly lower diversities of gut microbiota than the control mice (Figures 2A,B). Interestingly, the alpha diversity of the gut microbiome of both infected and normal mice increased significantly after OA B. subtilis (Figures 2A,B). Our findings suggested that OA B. subtilis can modulate the alpha diversity of gut microbiota in S. japonicum–infected mice.

To further examine how OA B. subtilis affects the composition of the gut microbiota of S. japonicum–infected mice, we next performed taxonomic analysis using MEGAN. Differences in the top 38 bacterial taxa were noted among groups (Figure 3). SIBS mice had an overabundance of most bacterial genera, including Bacteroides, Allporevotella, Enterococcus, Lactobacillus, Blautia, Lachnoclostridium, Anaerotruncus, Ruminiclostridium, Rhodospirillales, and Enterobacter (Figure 3). Interestingly, most bacterial genera belonging to Firmicutes, which are considered beneficial bacteria, were overabundant in the SIBS population. Moreover, the ternary plot of microbial communities in different groups also revealed that the relative number of microbes belonging to the phylum Firmicutes was closer to the SIBS population than to the SI group (Supplementary Figure S4).

Since OA B. subtilis significantly changed the gut microbiota composition and diversity of infected mice, we asked whether OA B. subtilis could also modulate the host response to schistosome infection. To this end, we performed small intestine and colon RNA-Seq in NG, SI, and SIBS mice after 8 weeks of infection. We found that the mRNA expression profiles between small intestine samples were strongly correlated (r ≥ 0.894), in addition to those between colon samples (Figure 4A). Meanwhile, there was also a strong correlation (r ≥ 0.739) between the small intestine and colon samples (Figure 4A). Similar results between sequenced samples are shown in the heatmap (Figure 4B). In addition, 275 of the total OTUs in the small intestine were commonly shared between the SI and SIBS groups (Supplementary Figure S5A). There were 330 OTUs in the small intestine uniquely identified in the NG versus (vs.) SIBS groups compared with the NG vs. SI groups. The shared and unique OTUs of expression profiles in the colon between groups are also displayed in the Venn diagram (Supplementary Figure S5B). These findings suggested that host immune system’s abilities between small intestine and colon to respond to intestinal and liver diseases are different.

The DEGs (based on fold changes) in small intestinal tissues between the NG and SI groups are shown as heatmaps in Figure 5A. Based on the assessments of the mRNA expression profiles of various genes in the small intestine of the NG vs. SI groups, we found that 514 genes were upregulated and 397 genes were downregulated (Supplementary Figure S6A; Supplementary Table S3). To identify the KEGG biological pathways enriched between the NG and SI groups, we revealed the top 10 most change enriched pathways, including the Fc epsilon RI signaling pathway, B cell receptor signaling pathway, Fc gamma R-mediated phagocytosis, intestinal immune network for IgA production, hematopoietic cell lineage, NF-kappa B signaling pathway, phospholipase D signaling pathway, pancreatic secretion, and phagosome (Figure 5B). Based on the KEGG classifications of DEGs, transport and catabolism, signal transduction, translation, lipid metabolism, and immune system ranked at the top of five KEGG categories, which included cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems (Supplementary Figure S7A). Based on the GO classification, these DEGs were classified into the cellular process, cell, and binding categories, ranking at the top of the biological process, cellular component, and molecular function classes, respectively (Supplementary Figure S8A). In addition, the DEGs (based on fold changes) in small intestinal tissues between the NG and SI groups are shown as heatmaps in Figure 5C. The comparative DEGs of small intestinal mRNA expression in profiles between the SI and SIBS groups revealed 541 upregulated and 348 downregulated genes (Supplementary Figure S6B; Supplementary Table S4). The KEGG pathway enrichment results between the SI and SIBS groups showed that the top 10 most prevalent pathways included the intestinal immune network for IgA production, Fc epsilon RI signaling pathway, NF-kappa B signaling pathway, Fc gamma R-mediated phagocytosis, calcium signaling pathway, B cell receptor signaling pathway, hematopoietic cell lineage, phospholipase D signaling pathway, pancreatic secretion, and phagosome (Figure 5D). The KEGG classification (Supplementary Figure S7B) and GO classification (Supplementary Figure S8B) are shown in the additional figures.

In addition, the comparative mRNA expression profiles in the colon of the NG vs. SI groups and SI vs. SIBS groups revealed that 1,693 genes (598 upregulated and 1,095 downregulated) (Supplementary Figure S6C; Supplementary Table S5) and 2,358 genes (1,783 upregulated and 575 downregulated) (Supplementary Figure S6D; Supplementary Table S6) were differentially expressed. They are shown as heatmaps in Figures 6A,C, respectively. The KEGG pathway enrichment results showed that the intestinal immune network for IgA production, hematopoietic cell lineage, B cell receptor signaling pathway, natural killer cell–mediated cytotoxicity, NF-kappa B signaling pathway, CAMs, Fc epsilon RI signaling pathway, phagosome, Fc gamma R-mediated phagocytosis, and Th17 cell differentiation were the top 10 most frequently enriched pathways in the colon between the NG and SI groups (Figure 6B). The KEGG classification (Supplementary Figure S7C) and GO classification (Supplementary Figure S8C) between the NG and SI groups are shown. In addition, CAMs, intestinal immune network for IgA production, hematopoietic cell lineage, Th1 and Th2 cell differentiation, Th17 cell differentiation, calcium signaling pathway, chemokine signaling pathway, NF-kappa B signaling pathway, B cell receptor signaling pathway, and phagosome were the 10 most change enriched pathways in the colon between the SIBS and SI groups (Figure 6D). Based on the DEGs between the SIBS and SI groups, the KEGG classification (Supplementary Figure S7D) and GO classification (Supplementary Figure S8D) were performed.