This study was conducted in May 2022. A total of 56 inhabitants living on Qingshan Island were engaged into the cross-sectional study, of whom 40 were finally enrolled for intestinal microbiota analysis and divided into two groups: patients with chronic schistosomiasis japonica-induced liver fibrosis group (CSJ group, n = 10) and a healthy control group (HC group, n = 30). CSJ participants were eligible for inclusion in the study if they met the following criteria: (a) with a clear history of S. japonicum infection; (b) history of treatment for S. japonicum; (c) have no prior history anti-helminthic or antibiotic treatment within 1 month (Schneeberger et al., 2018; Ramirez et al., 2020); (d) with liver fibrosis by ultrasound diagnosis, without ascites, portal hypertension and splenomegaly, etc., which referred to the previous studies (Abdel Wahab et al., 1992; Richter et al., 2003; Wu et al., 2015). Participants were excluded if they presented with cognitive or mental disabilities, gastrointestinal diseases, and/or other liver diseases, such as hepatitis, fatty liver and hepatocellular carcinoma. HC individuals were referred to reside on Qingshan Island, have no prior history of S. japonicum and test negative for S. japonicum. Fresh feces and blood were collected from all participants for analysis. Blood routine examination, including hemoglobin (HGB), neutrophils (NEU), eosinophils (EOS), lymphocyte (MON) and platelet (PLT), were tested by an automated hematologic analyzer (Sysmex XN-9000); Blood biochemical examination, including total triglycerides (TG), total cholesterol (TC), alanine aminotransferase (ALT), aspartate aminotransferase (AST), alkaline phosphatase (ALP), albumin (ALB), globulin (GLOB), and total protein (TP) were tested by an automatic biochemical analyzer (Hitachi 7600).

Total genomic DNA from fecal samples were extracted using OMEGA Soil DNA Kit (M5635-02, Omega Bio-Tek). DNA concentration and purity were monitored on 1% agarose gels. According to the concentration, DNA was diluted to 1 ng/μl using sterile water.

The V3–4 hypervariable region of the 16S rRNA genes was amplified using the specific primers with the barcode: 341F (5′-CCTAYGGGRBGCASCAG-3′) and 806R (5′-GGACTACNNGGGTATCTAAT-3′). All PCR reactions were carried out in 30 μl reactions with 15 μl of the Phusion^® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μM of the forward and reverse primers, and about 10 ng template DNA. Thermal cycling included an initial denaturation of 98°C for 1 min, followed by 30 cycles of denaturation at 98°C for 10 s, annealing at 50°C for 30 s, elongation at 72°C for 60 s, and finally heat preservation for 5 min at 72°C.

The same volume of 1× loading buffer (containing SYB green) was mixed with the PCR products and subjected to electrophoresis on a 2% agarose gel. Samples with a bright strip between 400 and 450 bp were chosen for further experiments. PCR products were mixed in equidensity ratios. Then, mixture PCR products was purified with Qiagen Gel Extraction Kit (Qiagen).

Sequencing libraries were generated using TruSeq^® DNA PCR-Free Sample Preparation Kit (Illumina) following manufacturer’s recommendations and index codes were added. The library quality was assessed on the Qubit^® 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. At last, the library was sequenced on an Illumina NovaSeq6000 platform and 250 bp paired-end reads were generated. Quality control on the raw reads was performed using Quantitative Insights into Microbial Ecology (QIIME 1.8.0) as previously described (Caporaso et al., 2010).

Paired-end reads from the original DNA fragments were merged using FLASH (Magoè and Salzberg, 2011), a very fast and accurate analysis tool designed to merge paired-end reads. When at least some of the reads overlap the read was generated from the opposite end of the same DNA fragment. Paired-end reads were assigned to each sample according to the unique barcodes. Sequence analysis was performed using the UPARSE software package with the UPARSE-Operational Taxonomic Units (OTU) and UPARSE-OTUref algorithms (Edgar, 2010, 2013). Sequences with ≥97% similarity were assigned to the same OTUs. A representative sequence was chosen for each OTU and the RDP classifier was used to annotate taxonomic information for each representative sequence (Silva 132 for 16S, UNITE for ITS). α (within sample) and β (among sample) diversities were estimated using MOTHUR software (v1.31.2) (Schloss et al., 2009) and the QIIME pipeline (v1.8.0), respectively. The Observed Species, Good’s coverage, Chao1, ACE, Shannon, and Simpson indexes were used to evaluate α-diversity. Shannon curves were generated based on these indexes. β-diversity analysis was calculated using the weighted UniFrac (Lozupone et al., 2011) and visualized by principal coordinates analysis (PCoA) downscaling. To measure differences in individual taxonomic abundance, linear discriminant analysis (LDA) effect size (LEfSe) was used to count species with significant differences in community structure between the samples. A LDA score >2 was considered a significant difference and used to screen for characteristic biomarkers of the intestinal microbiota between groups (Segata et al., 2011). Microbial functions were predicted by PICRUSt (phylogenetic investigation of communities by reconstruction of unobserved states) (Langille et al., 2013) in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Clusters of Orthologous Groups of proteins (COG) databases based on 16S rRNA gene sequence data. KEGG and COG analysis between two groups were conducted using LEfSe (Segata et al., 2011).

Statistical analyses were performed using SPSS software (v 26.0; SPSS, IL, USA). Continuous variables were expressed as the mean ± standard deviation (SD) or median [interquartile range (IQR)]. The Student’s t-test or Wilcoxon Rank-Sum and the χ² tests were used to evaluate differences between two groups for continuous and categorical variables, respectively. Receiver operation characteristic (ROC) analysis to predict the diagnostic efficiency of the different intestinal microbiota. The false discovery rate (FDR) was used to correct the calculated Q-values. A P-value < 0.05 was considered statistically significant.

A participant flow chart was shown in Figure 1. A total of 56 inhabitants of Qingshan Island participated in the cross-sectional study, of whom 40 were finally enrolled in the intestinal microbiota analysis based on the inclusion and exclusion criteria. These participants were divided into a CSJ group (n = 10) and a HC group (n = 30). The general characteristics of all participants were shown in Table 1, no significant differences were observed between the two groups.

In our study, 250 bp paired-end reads were generated by Illumina NovaSeq6000 platform. There were 3,060,101 sequences in total after demultiplexing and filtering for read quality, with a median of 77,322 sequences per experimental sample. The species accumulation curve of the samples reached a plateau, indicating that the study sample size was sufficient (Figure 2A). The Shannon curve indicated that the sequencing depth detected the most microbial information for each fecal sample (Figure 2B). To remove potential signal from contaminants, we removed from analysis taxa with >10% relative abundance in the blank swab samples and OTUs with <2 total occurrences among the samples. All reads clustered into 5,051 OTUs. Among them, 1,309 and 3,742 OTUs were detected in the CSJ and HC groups, respectively, and a Venn diagram indicated that 1,207 of these OTUs were found in both groups (Figure 2C). ANOSIM analysis was also performed and the between-group difference was greater than the within-group difference (R = 0.153, P = 0.001), indicating that the study groups were valid (Figure 2D).

Relative abundance of the intestinal microbiota at the phylum and family levels was also assessed. At the phylum level, the dominant microbiota of the two groups were similar, indicating that the disease status did not change the “Enterotype” of the patients. Firmicutes, Bacteroidetes, and Proteobacteria were the three most common phyla found in both groups, together accounting for 94.9 and 96.4% of the abundance in the CSJ group and HC groups, respectively. The abundance of Firmicutes and Bacteroidetes was higher in the CSJ group than in the HC group, while the abundance of Proteobacteria was significantly higher in the HC group than in the CSJ group (Figure 3A). At the family level, Lachnospiraceae, Prevotellaceae, and Erysipelotrichaceae were increased in CSJ group, while Enterobacteriaceae and Veillonellaceae were reduced (Figure 3B).

Linear discriminant analysis effect size was used to further explore the microbiota genera that differed significantly between the two groups. The abundance of more than 20 microbiota differed, however, some of these exhibited too low abundance at the family level to perform a clinical function. Combined with their relative abundance, Prevotella 7, Alloprevotella, and Holdemanella genera were significantly higher in the CSJ group than in the HC group (Figure 6A). We also did ROC analysis to explore the diagnostic efficiency of the Prevotella 7, Alloprevotella, and Holdemanella genera by SPSS. The results showed that the area under the curve (AUC) of Prevotella 7, Alloprevotella, and Holdemanella genera was 0.779, 0.769, and 0.840, respectively (Figures 6B–D), which indicated a potential value in non-invasive diagnosis of chronic schistosomiasis japonica-induced fibrosis.

Kyoto Encyclopedia of Genes and Genome and COG analyses were used to predict the function of the different intestinal microbiota in patients with chronic schistosomiasis japonica. KEGG analysis revealed that the Replication and Repair, and Translation pathways were significantly higher in the CSJ than in the HC groups (Figure 7A). COG analysis identified enrichment in the translation ribosomal structure and biogenesis and defense mechanisms after S. japonicum infection (Figure 7B).