Analysis of the relationship between bile duct and duodenal microbiota reveals that potential dysbacteriosis is the main cause of primary common bile duct stones.

Bacteria play an important role in the formation of primary Common Bile Duct (CBD) stones. However, the composition and function of the microbiota of bile duct in patients with primary CBD stones remained to be explored. We utilized the 16S rRNA gene high-throughput sequencing technology to analyze the microbial diversity and community composition of biliary and duodenal microbiota in 15 patients with primary CBD stones and 4 patients without biliary tract diseases. Alpha diversity analysis showed that the microbiota richness was similar in bile and intestinal fluid; Beta diversity analysis showed that there were differences in the composition between biliary microbiota and the duodenal microbiota, but the abundance of the main groups showed similarities. The composition of the biliary microbiota from gallstone patients was more complex, as was the duodenal microbiota. Proteobacteria and Firmicutes were the dominant bacteria at phylum level, accounting for at least 75% of the total reads in each subgroup. Pseudomonas and Escherichia-Shigella were the major genus among subgroups, but Escherichia-Shigella had increased abundance in duodenal microbiota with primary choledocholithiasis, which may play an important role in stone formation. It is noteworthy that Clostridiumsensu_stricto, Lachnospiraceae _UCG-008, Butyrivibrio and Roseburia which could produce short chain fatty acids (SCFAs), were significantly decreased in biliary microbiota with primary CBD stones (p < 0.05). Our study provided new insights into the compositional of normal biliary microbiota. The micro-ecology of biliary and duodenal in patients with stones is complex and closely related, and there is a potential for dysbacteriosis. The decrease in abundance of certain major acid-producing bacteria affects the health of the biliary tract and thus leads to the formation of stones.

Introduction

Even being a benign disease, primary choledocholithiasis can cause high morbidity and serious complications in China. The morbidity of primary CBD stones accounts for 5%–30% of patients with cholelithiasis [1], and its main component is brown pigment stones in China. The pathogenesis of primary CBD stones is not clear yet. Therefore, the study of the causes and mechanisms of primary CBD stones has always been the focus and hotspot in order to prevent and treat it better. Among the many risk factors, bacteria are closely related to the formation of primary CBD stones, and it had long been confirmed that bacteria are present in bile with CBD stones [[2], [3], [4], [5], [6]]. In addition, the intestinal microbiota has been considered as the main source of biliary bacteria, and mainly from the upper digestive tract [7]. Intestinal dysbacteriosis is associated with the development of various diseases such as metabolic syndrome, Crohn's disease, colon cancer, and cholecystolithiasis [[8], [9], [10], [11]]. We assumed that, the formation of CBD stones may be closely related to the bile duct and duodenal microbiota alterations.

The previous studies could not provide an overall explanation of the composition and function of the biliary microbiota, especially the normal biliary microbiota. What's the relationship between biliary and gut microbiota and how do they affect the formation of primary CBD stones have rarely been reported. With the development of science and technology, high-throughput sequencing based on 16S rRNA gene PCR amplification has been widely used in microbial research because of its high throughput, high sensitivity, and rapid sequencing [12].

Based on the 16S rRNA gene high throughput sequencing technique, the bile and duodenal juice of the patients with primary CBD stones and without biliary tract disease (control group) in our hospital are sequenced using the Illumina Miseq sequencing platform, the composition and diversity of the microbiota were analyzed, and the differences between the bile duct bacteria and intestinal bacteria were compared, to explore its significance for the formation of primary CBD stones.

Material and methods

Study design and sample collection

We selected the inpatients at the Second Hospital of Hebei Medical University, through clinical symptoms, signs, ultrasonography, abdominal CT and other imaging examinations, and finally confirmed by ERCP as primary CBD stones patients (EG) without biliary tract infections. In EG, there were 15 cases, including 10 males and 5 females, aged 68.13 ± 13.53 years, and 4 patients without hepatobiliary diseases (CK), all were male and aged 39.25 ± 12.09 years old. Inclusion criteria were as follows: (1) Patients with related clinical manifestations of choledocholithiasis and diagnosed according to the relevant diagnostic methods, or without any hepatobiliary diseases. (2) No major diseases such as severe heart, brain, kidney, etc. (3) Good condition, could undergo ERCP. (4) signed informed consent. Exclusion criteria: (1) antibiotics or probiotics were applied within 3 months before operation. (2) coexisting cholelithiasis. (3) common bile duct obstruction caused by malignant tumors.

The four subgroups, namely EG group bile (B subgroup) and duodenal juice (D subgroup) and CK group bile (BCK subgroup) and duodenal juice (DCK subgroup) were respectively extracted from two groups of patients under ERCP.

The entire ERCP process was strictly aseptic and all experimental materials were sterilized according to the procedure. During the entry of the duodenoscope (TJF240/JF260V; Olympus), any inhalation was strictly prohibited to avoid contamination. Before duodenal fluid was drawn, the duodenoscope channel was washed with 10 ml acidified water (pH = 2.5 ± 0.2) and 20 ml physiological saline at the level of the duodenum, and used balloon extraction catheter to draw at least 2 ml of intestinal fluid within 3 cm of the duodenal papilla. Then, the channels and papilla were flushed again with 20 ml of normal saline, and 5–10 ml of bile was extracted from the common bile duct using a triple lumen sphincterotome. The samples were placed in sterile sputum cup and immediately stored in the −80 °C refrigerator for later bacterial DNA extraction and high throughput sequencing analysis.

All patients in this study expressed their understanding and signed informed consent and obtained approval from the Ethics Committee of the Second Hospital of Hebei Medical University.

16S rRNA gene amplification and sequencing of amplicons

The genomic DNA of bile and duodenal fluid samples were extracted by means of Soil DNA Kit (Omega). The obtained DNA was subjected to 0.8% agarose gel electrophoresis for molecular size determination, and the DNA was quantified by ultraviolet spectrophotometer. The DNA is stored in the −20 °C refrigerator for subsequent analysis after passing the test. In order to ensure the quality of sequencing, the insertion fragment range of the best sequencing was 200–450bp by Miseq sequence reading length. In this experiment, the V3∼V4 region of 16S rRNA gene was selected for sequencing. The PCR primers used were 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′). Conventional PCR thermal cycling program: initial denaturation at 98 °C for 30s, denaturation at 98 °C for 15s, annealing at 50 °C for 30s, extension at 72 °C for 30s. The above three processes were carried out for 25 cycles and lastly extended for 5 min at 72 °C, finally stored at 4–10 °C conventionally. The amplification results were performed on 2% agarose gel electrophoresis. After passing the test, the target fragment was excised and the target fragment was recovered using the Axygen Gel Recovery Kit. The PCR products were quantified on a Microplate reader (BioTek, FLx800) using the Quant-iTPico Green dsDNA Assay Kit, and then mixed according to the amount of data required for each sample. The obtained PCR products above were used to construct a sequencing library by using Illumina truseq. The Qualified sequencing library (index sequence cannot be repeated) was diluted according to the gradient, then mixed according to the required amount of sequencing, and denatured by NaOH into single strands for on-machine sequencing. Barcoded V3–V4 amplicons were sequenced using the pair-end method by Illumina Miseq [13]. The optimal sequencing length was chosen to be 200–450 bp. Using FLASH software [14] to pair pair-end sequences based on the initial screening according to the overlapping bases: the overlapped bases of both Read 1 and Read 2 sequences are required to be ≥ 10 bp in length, and base mismatches are not allowed. Afterwards, the valid sequences were obtained according to the Index information corresponding to each sample. Finally, use QIIME software (Quantitative Insights Into Microbial Ecology, v1.8.0) [15] to identify the error sequence and remove it. Amplification and sequencing of 16S rRNA gene was completed by Personal Biotechnology Co., Ltd. (Shanghai, China).

OTU clustering and statistical analysis

Operational Taxonomic Unit (OTU) [16]usually refers to a set of human-defined sequence similarity thresholds. The sequences from one or more samples are merged, and the sequences whose similarity is higher than the threshold will be merged into one OTU. Using the UCLUST sequence comparison tool in QIIME software [17], the sequences obtained above were merged and OTU-divided with a sequence similarity of 97%, and selected the highest abundance sequence in each OTU as the representative sequence of the OTU. The OTU with an abundance value less than 0.001% of the all sample total sequencing was removed [18], then used Silva database (Release 132,http://www.arb-silva.de) as the template sequence for the identification of the OTU taxonomic status.

Statistical analyses were performed using the SPSS 23.0 software. The Fisher's exact test was used to enumeration data, and the Kruskal-wall is test to measurement data. The OTU was classified by QIIME software using the Silva database and the bacterial annotation was carried out. Differences within or between groups were compared by Alpha diversity, Beta diversity, Metastats analysis, PLS-DA analysis, and Spearman correlation network analysis. PICRUSt was used to predict bacterial function through the KEGG database.

Results

Sequencing data and bacterial diversity analysis

Through the metagenomic DNA extraction, 16S rRNA gene amplification and amplicon sequencing of the 38 biliary bile and duodenal juice samples from 19 persons including 15 patients and 4 controls, a total of 1,811,427good reads were obtained and 2801 OTUs were annotated. Sequencing data statistics and diversity analysis each sample yielded 47669.13 sequences (range = 33338–76885, SD = 10186.88) on average. The Goods coverage (%) of each group was above 99% and was no statistically significant different (p > 0.05), that is, the amount of sequencing of each sample is nearly saturated. Thus, the study has a very high reliability and reference value. After statistical analysis (Table S1), it was found that the bacterial taxa in bile were less than that of duodenal fluid at the level of phylum, genus, and species in the experimental group (p < 0.05).

The diversity indices (Chao 1, ACE, Simpson and Shannon) were applied to reflect its alpha diversity. The species accumulation curve in each sample reached plateau phase, which confirmed the sequence numbers is large enough to cover most of the microbial information in the sample, and abundance and evenness of the microbial species in the current study meet the requirements for sequencing and analysis. The Shannon, Simpson, Chao 1 and ACE index demonstrated there was no striking difference in the microbial community abundance between groups (p > 0.05) (Fig. 1). This result indicated that the microbiota richness is similar in bile and intestinal fluid.

Fig. 1

Microbial alpha diversity with a box plot exhibiting the community diversity.

Beta diversity analyses were used by employing Principal component analysis (PCA) and Multidimensional scaling (MDS) of the similarity of microbial community structure in different groups. Mapping data were obtained from UniFracPCoA master coordinate analysis via QIIME [19,20]. Based on the Unweighted UniFrac distance, almost all bile samples were on the left side of the vertical axis and the duodenal fluid samples were on the right side (Fig. 2a), and Non-metric Multidimensional scaling (NMDS) analysis based on Unweighted UniFrac distance obtained similar mapping data (Fig. 2c). On the contrary, based on the Weighted UniFrac distance, the distribution of all samples was more concentrated and some of them were clustered and overlapped (Fig. 2b), and NMDS analysis based on Weighted UniFrac distance also showed similar results (Fig. 2d). It indicated that there were significant divergence in the composition between biliary microbiota and the duodenal microbiota. ANOSIM testing (Table S2)confirmed that significant separation occurred between B subgroup and D subgroup (Unweighted UniFrac, p = 0.001; weighted UniFrac, p = 0.001), as well as between BCK subgroup and DCK subgroup (Unweighted UniFrac, p = 0.028; weighted UniFrac, p = 0.021).

Fig. 2

Multidimensional scaling (MDS) of the similarity of microbial community structure in different groups. (a. Unweighted UniFracPCoA; b. weighted UniFracPCoA; c. Unweighted UniFrac NMDS; d. weighted UniFrac NMDS.)

We presented the proportion of shared and unique OTUs in each subgroup by using the Venn diagram visually (Fig. 3). There were 1205 shared OTUs among the four subgroups, occupying 43% of all OTUs (total 2801 OTUs). More than 70% OTUs were shared in B subgroup and D subgroup, but only about 60% shared in BCK subgroup and DCK subgroup. The most dominant shared phylum and genus were Proteobacteria (95%) and Pseudomonas (83%), respectively. In the Venn diagram, we could also observe the unique OTUs in each subgroup. The unique OTUs of genus in the B subgroup were mostly Pseudomonas, which corresponding species were unidentified. Only 6 unique OTUs of genus were observed in BCK subgroup: Dechlorobacter, Halomonas, Pseudomonas and other unidentified genera.

Fig. 3

A Venn diagram showing shared and unique OTUs among the four subgroups.

Analyses the microbial community structure in different groups

According to the OTU division and classification status, the specific microbial composition of samples at each classification level can be obtained. A total of 23 phyla were identified in this work. The microbial composition varied among subgroups at the phylum and genus level, but Proteobacteria and Firmicutes accounted for at least 75% of all phyla (Fig. 4(a)), and Pseudomonas and Escherichia−Shigella were the major genus (Fig. 4(b)).

Fig. 4

The microbial community structure at the phylum level (a) and genus level (b) in different subgroups.

Comparison of the microbial community structure in B subgroup and BCK subgroup

At the Phylum level, a total of 18 phyla were identified in the two groups. Proteobacteria was the most abundant phylum, accounting for more than 90% in each group. Metastats (http://metastats.cbcb.umd.edu/) was utilized to analysis the microbial community difference between groups [21]. Four phyla (Actinobacteria (p = 0.007867,q = 0.068733), Saccharibacteria (p = 0.005733,q = 0.068733), Deinococcus-Thermus (p = 0.023533,q = 0.154212) and SR1 (p = 0.0068,q = 0.68733)) were discovered to have significant difference between B subgroup and BCK subgroup (Fig. 5a).

Fig. 5

Comparison of microbial community structure at Phylum and Genus level in different groups by Metastats. a. B-BCK (Phylum); b. B-BCK (Genus); c. D-DCK (Genus); d. B-D (Phylum); e. B-D (Genus); f. BCK-DCK (Phylum); g. BCK-DCK (Genus).

A total of 232 genera were identified in the two groups, of which the most dominant genus Pseudomonas and sub-dominant Escherichia-Shigella accounted for more than 80% of the two groups, both of which belonged to the Gammaproteobacteria of Proteobacteria. Among 232 genera, Actinomyces (p = 0.000143,q = 0.011615), Proteus (p = 0.000558,q = 0.022678), Clostridiumsensu_stricto (p = 0.000721,q = 0.023452), Klebsiella (p = 0.001571,q = 0.028393), Actinobacillus (p = 0.000143,q = 0.011615), Lachnospiraceae_UCG-008 (p = 0.002388,q = 0.029789), Butyrivibrio (p = 0.002565,q = 0.029789), Roseburia (p = 0.003211,q = 0.031885), Porphyromonas (p = 0.010422,q = 0.080702), Streptococcus (p = 0.012286,q = 0.086862), Helicobacter (p = 0.018143,q = 0.118012), Prevotella 1(p = 0.037156,q = 0.163601), 7(p = 0.044231,q = 0.174967), Enterobacter (p = 0.039463,q = 0.168874) et al. showed significant difference with p < 0.05. It is notable that except for Brevundimonas and Prevotella_1, the abundance of all the different genera in B subgroup was lower than that in BCK subgroup (Fig. 5b).

Comparison of the microbial community structure in D subgroup and DCK subgroup

A total of 22 phyla were identified in the two groups, of which Proteobacteria was the most dominant phylum, accounting for more than 70% of OTUs in both groups, and Firmicutes was the sub-dominant phylum. No phylum showed significant difference.

A total of 291 genera were identified in the two groups, and the dominant genera were Pseudomonas, Escherichia-Shigella and Streptococcus in sequence. 14 genera were significant different with p < 0.05 between two groups. The abundance of Escherichia-Shigella (p = 0.0336,q = 0.698922), Geobacillus (p = 0.008794,q = 0.418704), Stomatobaculum (p = 0.027543,q = 0.698922), Campylobacter (p = 0.022549,q = 0.698922) and Solibacillus (p = 0.04092,q = 0.698922) in D subgroup was significantly increased, the remaining 9 genera such as Ochrobactrum (p = 0.004337,q = 0.418704), Moraxella (p = 0.0076,q = 0.418704), etc. increased significantly in DCK subgroup (Fig. 5c).

Comparison of the microbial community structure in B subgroup and D subgroup

Similarly, Proteobacteria and Firmicutes were still the most abundant phyla. The bacterial phyla in bile were contained in the duodenal fluid, however the duodenal fluid additionally contained Nitrospirae and Gracilibacteria. Eight phyla were significant different with p < 0.05 between B subgroup and D subgroup. In addition to the abundance of Proteobacteria (p = 0.000999,q = 0.003959) and Cyanobacteria (p = 0.036963,q = 0.109868) increased, the others decreased significantly in bile (Fig. 5d).

There were 295 genera were identified in the two groups, of which Pseudomonas, Escherichia-Shigella and Clostridiumsensu_stricto_1were the dominant genera in sequence. There are 15 unique genera in the bile, but there are up to 40 unique genera in the duodenal fluid. Even more striking is that 69 genera were significant different with p < 0.05 between two groups (Fig. 5e). Among 69 genera, the abundance of Pseudomonas (p = 0.007992,q = 0.045899)was significantly increased in the bile, others such as Streptococcus (p = 0.000999,q = 0.009562), Veillonella (p = 0.000999,q = 0.009562), Lactobacillus (p = 0.001998,q = 0.015187), Prevotella_7 (p = 0.000999,q = 0.009562), Prevotella_6 (p = 0.026973,q = 0.126743), etc. increased significantly in duodenal fluid.

Comparison of the microbial community structure in BCK subgroup and DCK subgroup

A total of 18 phyla were identified in the two groups, of which Proteobacteria, Firmicutes and Actinobacteria were the dominant phyla in sequence. Only phylum Bacteroidetes (p = 0.004263,q = 0.052201) was found obvious difference between bile and duodenal fluid in CK, of which the abundance was significantly increased in the duodenal fluid (p < 0.05) (Fig. 5f).

At the genus level, we identified 234 genera from the two groups. Similarly, Pseudomonas and Escherichia-Shigella in the phylum Proteobacteria are in turn the two most dominant. 17 genera were found obvious significant difference between two groups (Fig. 5g), of which only the abundance of Aquabacterium (p = 0.018142,q = 0.406811) and Enhydrobacter (p = 0.048049,q = 0.442801)was significantly increased in the duodenal fluid.

Key genera screening and association analysis

Partial least squares discriminant analysis (PLS-DA)

The PLS-DA discriminant model was constructed based on the species abundance matrix and sample grouping data. Just like Fig. 6 showed, both bile samples and duodenal fluid samples received relatively better grouping models. At the same time, at the genus level, we calculated the VIP (Variable importance in projection) coefficient for each genus (VIP>1, the greater the value, the greater the contribution of the species to differences between groups). In bile samples, there were 75 genera that contributed to the bile grouping (Table S3). Correspondingly, 87 contributing genera were discriminated for duodenal fluid grouping in the duodenal fluid (Table S4).

Fig. 6

The PLS-DA discriminant model was constructed in bile group (a) and duodenal fluid group (b).

Spearman Association Network Analysis of Dominant genera

We calculated the Spearman rank correlation coefficient between the dominant Genera, whose abundance was located in the top 50, and constructed the associated network for the related advantage where rho >0.6 and P value < 0.01 [22]. The results showed (Fig. 7) that there were abundant correlated networks between EG and CK. Most of the bacterial genera existed in the same environment with collaborative manners. However, the associations among bacterial genera in EG were more complicated. In CK, the relationship between the genera Ralstonia and Ochrobactrum, Streptococcus and Granulicatella, Acinetobacter and Alloprevotellawere antagonistic to each other (Fig. 7a), but they showed cooperative relationships in EG (Fig. 7b).

Fig. 7

Network plots highlight correlations among bacterial genera in EG (a) and in CK (b). (Red lines refer to cooperative correlation and green lines refer to antagonistic correlation).

4 Microbial function prediction

In order to further understand the differences in the functions of biliary and intestinal microbiota, we performed functional prediction of the 16S rRNA gene sequence by PICRUSt in the KEGG database [23]. At KEGG level 3, we obtained a total of 304 pathways. 14 pathways were discovered to have significant difference in B subgroup and BCK subgroup (p < 0.05), and they were down-regulated in B subgroup than those in BCK subgroup. 8 pathways were significant different with p < 0.05 in duodenal fluid samples between D subgroup and DCK subgroup, including 6 up-regulated pathways and 2 down-regulated pathways.

The 6 major metabolic pathways, Metabolism, Genetic Information Processing, Environmental Information Processing, Cellular Processes, Organismal Systems and Human Diseases, all showed differences in some pathways at KEGG level 3 in bile samples between EG and CK. Microbiota in B subgroup had lower abundance of functions involved in Metabolism pathways such as Cysteine and methionine metabolism, Glycerophospholipid metabolism and Tetracycline biosynthesis than those in BCK subgroup (p < 0.05). In addition, the abundance of cellular processes (Apoptosis and Endocytosis), environmental information processing (Two-component system and CAM ligands) and the immune system (RIG-I-like receptor signaling pathway) decreased significantly in B subgroup (p < 0.05) (Fig. 8a).

Fig. 8

Difference of functional pathway at KEGG level 3 in bile samples (a), duodenal fluid samples (b) between EG and CK.

In duodenal fluid samples (Fig. 8b), microbiota in D subgroup had higher abundance (p < 0.05) of functions involved in Metabolism pathways (Purine metabolism, Metabolism of xenobiotics by cytochrome P450 and Naphthalene degradation) and human diseases (Bacterial invasion of epithelial cells, Measles). But the abundance of Digestive system (Pancreatic secretion) and nervous system (Glutamatergic synapse) decreased significantly in D subgroup (p < 0.05).

Discussion

Due to differences in lifestyle, dietary structure, and living environment, primary CBD stones maintain a high incidence in China, which has brought about property damage and health hazards for our country and people. As an important condition for the occurrence of the disease, bacteria have been studied by scholars. Most previous experiments were based on cholecystolithiasis or in vitro model to infer the relationship between the two. In addition, the intestinal microbiota plays a role in human diseases with its rich functions. There are a number of diseases such as Crohn's disease, colon cancer, and metabolic diseases that have been shown to be related to the intestinal flora [[8], [9], [10]]. The duodenal microbiota, as the nearest neighbor to the biliary tract, was considered to be the main source of possible biliary bacterial infections [24]. In this study, bile and duodenal fluid were directly obtained through ERCP techniques and diversity analysis was performed. Through the rigorous screening of the subjects, we obtained a total of 19 samples of patients. It is worth noting that 4 of them as a control group were pancreatic-related diseases patients without hepatobiliary system disease. Although we were unable to obtain samples from healthy people, these 4 patients still have research significance, at least for their normal bile. The study was ethical and informed consents were obtained from the patients. To ensure that samples were not contaminated with each other was a difficult point in the sampling process. Thus, the duodenal fluid was first taken and then bile was extracted under the premise of strict aseptic processing to ensure the fidelity of this study.

In the field of microbial research, the widespread application of high-throughput sequencing technology has broken the limitations of traditional culture techniques. In our study, the 16S rRNA gene high-throughput sequencing technology was used for biliary and intestinal (duodenal) microbiota in patients with CBD stones and in patients without hepatobiliary system diseases. 748,340 high-quality gene sequences were obtained in B subgroup, and 2062 bacterial OTUs were identified, which is similar to previous studies by Shen et al. [25]. Interestingly, at the same time, we obtained a large number of effective high-quality sequences of 172,025 and 1597 OTUs in BCK subgroup (without hepatobiliary system disease). In the past, based on bacterial culture, it was found that there was no bacteria in normal bile which without biliary tract infections, because of its potent antimicrobial properties [26]. Some scholars have also demonstrated the existence of bacteria in bile without common bile duct stones. For example, Pereira P et al. found the presence of bile bacterial communities in PSC patients [27]. Liu et al. found bacterial structures in acute liver failure and normal mice [28]. However, none of the patients enrolled in this study had the above-mentioned liver and biliary tract diseases, which provided valuable resources for the study of normal biliary microbiota.

The sequence quantity and the number of OTUs identified indicated that the bile duct bile and duodenal juice are rich in bacteria, regardless of whether the patients have hepatobiliary disease. At the same time, we observed that in the same group, the Alpha diversity of bile and duodenal microbiota was similar. But there was also study showing that the Alpha diversity of biliary microbiota was higher than that of gut microbiota and had a higher richness [11]. The factors for this difference may be related to the sites of the bile samples (gallbladder or common bile duct) and the sites of the intestinal microbiota samples (fecal or duodenal fluid).

We found that the most dominant phylum in bile and duodenal fluid was Proteobacteria in obtained 23 phyla. Other phyla with higher abundances included Firmicutes, Actinobacteria, Bacteroidetes, etc., this is similar to previous studies [6,29,30]. The difference is that Synergistetes (0.63%) with a high abundance in B subgroup were not found in BCK subgroup. Synergistetes, the main pathogens of oral diseases such as periodontitis [31], may also participate in the formation of CBD stones. The analysis of the composition of the biliary and duodenal microbiota showed that phyla contained in the bile could be found in the duodenal fluid almost. Similar result was also found in previous work of Ye et al. [7], their work showed the biliary microbiota had a comparatively higher similarity with the duodenal microbiota. However, the abundance of Proteobacteria and Cyanobacteria in the bile of the experiment group was significantly higher than that in the duodenal fluid. It shows that although the composition of biliary and duodenal bacterial community is similar, but the composition ratio is different.

At the genus level, the most dominant genus was Pseudomonas of the phylum Proteobacteria, followed by Escherichia-Shigella. Although they are classic contaminants [32], this does not invalidate the results as it has been shown that they are common in the small intestine and bile [33,34]. This is different from previous research results, Ye et al. found three Enterobacteriaceae genera (Escherichia, Klebsiella, and an unclassified genus) and Pyramidobacter were abundant in bile of six gallstone patients, while Li et al. found Acinetobacter and Prevotella were dominant genera in the duodenum of nine healthy volunteers [7,35]. According to Metastats and PLS-DA, Pseudomonas (VIP<1) and Escherichia-Shigella (VIP<1), the two most abundant genera may only exist as normal flora in bile, not as a risk factor for CBD stones. But Escherichia-Shigella plays an important role in the duodenum of patients with CBD stones. The bile or duodenal fluid in EG contained more genera than CK. The results showed that although the floras of the samples were relatively conservative at phylum level, there was dysbacteriosis and the imbalance of the composition ratio in the bile or duodenal fluid on genus level in EG. Meanwhile, the associations among bacterial genera in EG were more complicated by Spearman Association Network Analysis. It further proves that dysbacteriosis occurs in both the biliary tract and the duodenum in the relationship among bacterial genera.

Studies have shown that changes in the bacterial community of the duodenum were associated with various diseases such as small intestinal bacterial overgrowth (SIBO), irritable bowel syndrome (IBS) and celiac disease (CD) [36]. It was found that Clostridium sensu_stricto (VIP>1), Lachnospiraceae_UCG-008 (VIP>1), Butyrivibrio (VIP>1), Roseburia (VIP>1), which can produce short chain fatty acids (SCFA), were significantly reduced in bile from patients with CBD stones. SCFA has the functions of sterilization, bacteriostasis, and provide energy for host antibodies to improve immunity [37,38]. SCFA not only affects the intestinal epithelium, but has a certain effect on airway epithelium, enteric nervous system (ENS), central nervous system (CNS), and peripheral neurons [39]. The occurrence of multiple diseases is associated with a decrease in the proportion of “good bacteria” that can produce SCFA in the intestine, such as inflammatory bowel disease [40], type 2 diabetes [41], Asthma [42] and so on. So we can conclude that the low expression of “good bacteria” in bile has a negative impact on biliary health in patients with common bile duct stones, which ultimately leads to the development of stones. In addition, Klebsiella (VIP>1) and Helicobacter (VIP>1) have been reported in previous studies to be associated with CBD stone formation. For example, Klebsiellapneumoniae and Helicobacter pylori have been cultured by bacteria or identified by PCR technique in bile from patients with CBD stones [25,43]. Similar results were obtained in this study. But the two genera were more abundant in BCK subgroup. This phenomenon suggests that the disorder of the biliary bacteria community may be more important than the number of related pathogenic bacteria in the formation of CBD stones.

The comparison of biliary microbiota and duodenum microbiota indicated that gut microbiota showed a higher diversity of genera, both in EG and in CK. At the same time, the specificity and complexity of the two bacterial communities could be seen from as many as 60 differential genera between biliary microbiota and duodenum microbiota in EG. According to the Spearman Association Network Analysis of Dominant genera in EG and the normal group, we found that the relationship between Ralstonia and Ochrobactrum, Streptococcus and Granulicatella, Acinetobacter and Alloprevotellawere antagonistic to each other in CK, but they showed cooperative relationships in EG. This shift in bacterial relationship may be the biological basis for the occurrence and development of biliary tract diseases (e g. primary CBD stones). The reasons for this shift are not referenced in the relevant literature, may be related to the elimination or reduction of specific genus, or may be affected by other genera. These speculations need further research to explore.

The relationship between biliary and intestinal microbial communities has been a subject of continuous research and discussion. Bile bacterial infection was an important factor in the formation of CBD stones. Based on human anatomy and related studies, biliary microbiota in patients with CBD stones originated from the intestinal tract [44,45]. In recent years, research has also tried to confirm this view. Ye et al. provided evidence that the biliary flora in patients with CBD stones originated in the intestine and mainly from the upper gastrointestinal tract by the 16S rRNA gene amplicon sequence [7]. In this study, the beta diversity analysis of the sample community structure showed that the biliary and the duodenum microbiota showed some similarities in both EG and CK, especially in terms of bacterial population abundance (slightly different in bacterial composition), moreover this similarity was more evident in the EG. At the same time, a Venn diagram showed “core microbiome” (1205 shared OTUs), which was first proposed by Turnbaugh et al. [46]. Thus we can conclude that the most of biliary microbiota and the duodenal microbiota have a common source. More importantly, this study confirmed the existence of self-microbiota in the normal biliary tract, which constituted a microenvironment of bile together with bile salts, bile acids, and so on. In addition, though the abundance of unique OTUs was low in four subgroups, we still have the possibility to find biomarkers related CBD stones among them.

In order to gain more understanding of the biliary and duodenal microbiota, we performed microbial function prediction by PICTUSt. Both the duodenal and the biliary microbiota have gained rich function pathways. There are differences in some pathways at KEGG level 3 in bile samples between EG and CK, such as cellular processes (Apoptosis, Endocytosis), environmental information processing (Two-component system, CAM ligands) and the immune system (RIG-I-like receptor signaling pathway), of which the abundance in EG decreased. We believe that some genera appeared to have a decline of self-repair ability and environmental adaptability in bile in EG and were in a “pathological state.” In duodenal fluid samples, microbiota in EG had much higher abundance of functions involved in Naphthalene degradation and Bacterial invasion of epithelial cells. Studies have shown that Pseudomonas sp. JLR11 and E. coli were related to pathway of Naphthalene degradation [47,48]. Some bacteria such as Shigella, Streptococcus and Listeria could enter the epithelial cells through pathway of Bacterial invasion of epithelial cells [[49], [50], [51]]. It is indicated that the increase of bacteria with invasive ability of duodenal microbiota may indirectly lead to the formation of the stone, which may also explain the increase in the abundance of Escherichia-Shigella in the duodenal fluid of EG. However, in the Infectious Diseases pathway, Shigellosis pathways and pathogenic Escherichia coli infection pathways had no difference in duodenal microbiota between two groups. This is also a good illustration of why the high abundance of Escherichia-Shigellain duodenal flora does not cause related diseases.

This study further explored the composition of the biliary and intestinal microbiota in patients with primary CBD stones and the differences in key genera. And we had a completely new understanding of the biological communities in the normal biliary tract. However, this experiment still has its limitations. The small sample size may lead to certain deviations. Controls are not healthy volunteers, they could suffer from dysbacteriosis as well. As far as research methods are concerned, in the future, animal experiments and other methods may be needed to study and validate the biliary microbiota, so as to further clarify the role of bacteria on the formation of primary CBD stones.

Conclusion

In conclusion, our results give a more comprehensive disclosure of the composition of the biliary microbiota for the first time. The microecology of biliary and duodenal in patients with stones is complex, potential dysbacteriosis of them is of great importance in the formation of CBD stones, which is mainly reflected in the reduction of some dominant bacterium and the change of the association network of dominant genera. These are new understandings of the biliary microbiota and pathogenic mechanisms. Equally important is the discovery of reduction of bacteria which could produce short chain fatty acids (SCFA), the decrease may affect the health of the biliary tract and lead to the formation of stones. Maybe this discovery will guide us in the prevention and treatment of primary CBD stones in the future.

Availability of data and materials

16S rRNA sequence reads of 38 samples have been submitted to the NCBI SRA metadata under accession code SRP155620.

Ethical approval

This study was approved by the ethics committee of the second hospital of Hebei Medical University, Shijiazhuang, P. R. China.

Informed consent

Informed consent was obtained from all individual participants for the information included in this article.

CRediT authorship contribution statement

Zhitang Lyu: Managed the project, designed the analysis, wrote the paper. Tingting Yu: Performed sample collection and clinical study, performed DNA experiments, performed the data analysis, wrote the paper. Lichao Zhang: Performed sample collection and clinical study. Xiaona Xu: Performed DNA experiments. Yijun Zhang: Performed sample collection and clinical study, performed the data analysis. Jihong Li: Performed the data analysis. Zhirong Li: Performed DNA experiments. Wei Zhang: Performed sample collection and clinical study. Senlin Hou: The project idea was conceived and the project was designed by, designed the analysis. All authors have read and approved the manuscript.

Declaration of competing interest

The authors have no conflict of interest and competing Interests to declare.