Biliary Microbial Community and Metabolic Potential in Patients with Multiple Common Bile Duct Stones

Introduction

Cholelithiasis is a global health problem. Common bile duct stones (CBDS) are primarily caused by the migration of gallstones and present in 10–20% of individuals with symptomatic gallstones.¹ Extensive studies revealed the crucial role of abnormal cholesterol and bile acid metabolism in the development of CBDS.² Different from the previous view on sterile bile, biliary system is considered to harbor a microbial community based on the experimental results from next-generation sequencing (NGS), particularly in patients with calcium bilirubinate stones.³ Niches of biliary microbiota were demonstrated in different population, including healthy donors⁴ and those with hepatobiliary diseases.

Emerging studies^5–9 highlighted the influence of biliary microbiota on CBDS formation. Bacteria, adhering to pigment solids, can alter the physical and chemical properties of bile via their enzymatic activities, thereby promoting the formation of pigment stones.⁵ The bacteria expressing β-glucuronidase and phospholipase have been shown to be involved in the formation of pigment stones.⁵ The β-glucuronidase can hydrolyze bilirubin diglucuronide into glucuronic acid and unconjugated bilirubin, which subsequently combines with calcium ions in bile to form pigment stones. The phospholipase can hydrolyze phosphatidylcholine, a protective component of biliary epithelial cells, into free fatty acids, which contribute to the formation of pigment stones. For another thing, the stones provide a protective environment where the flora can establish a distinct niche resistant to the antimicrobial effects of bile. Thus, the conditions of bile will become more favorable for the growth of stones.⁶

Nowadays, endoscopic retrograde cholangiopancreatography (ERCP) is the preferred treatment strategy for CBDS. Successful management is often hampered by a variety of conditions, including numerous of stones, unusual shape, incarcerated location and large diameter,^10,11 which also increases the morbidity of postoperative complications.¹⁰ The previous studies usually focused on the composition and function of biliary flora involved in secondary CBDS⁷ and recurrent CBDS.⁸ In the prior work, we focused on the characteristics of biliary microbiota in giant CBDS.⁹ The correlation between biliary microbiota and the number of CBDS remains unclear.

In this study, we employed 16S rRNA sequencing to reveal the distinct microbial community present in the bile from patients with multiple CBDS. And then we used bioinformatic analyses to explore the potential metabolic functions that affect CBDS formation. Furthermore, we conducted quantitative analysis on 15 bile acids using ultra-performance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS) system to evaluate the association among biliary microbiota, bile acid metabolism and CBDS.

Methods

Study Design and Patient Enrollment

The present study was designed as a case–control study. Consecutive patients were prospectively enrolled from December 2022 to October 2023 at the First Affiliated Hospital of Soochow University. This study complies with the Declaration of Helsinki. It was approved by the ethics committee, and all the patients provided written informed consent. The inclusion criteria were: (1) CBDS were confirmed through at least one of the following examinations, including abdominal ultrasound, computed tomography (CT), magnetic resonance cholangiopancreatography (MRCP) and endoscopic ultrasound (EUS); (2) no use of antibiotics within the past three months. The exclusion criteria were as follows: (1) presence of acute obstructive suppurative cholangitis; (2) history of cholecystectomy; (3) history of upper gastrointestinal tract surgery; (4) history of gastrointestinal or hepatobiliary-pancreatic neoplastic diseases. The number of CBDS was ultimately confirmed during ERCP procedure. Patients were divided into two groups based on the number of stones, single stone group (SS group, number of stones = 1) or multiple-stone group (MS group, number of stones ≥ 2).

Sample Collection

ERCP was performed on the patients who suffered from CBDS after the informed consent was obtained. The duodenoscope and all instruments were strictly disinfected to maintain the sterility of working channel. Based on standard ERCP protocol, a guide wire was firstly inserted into common bile duct under X-ray guidance. Bile samples were aseptically aspirated through a sterile catheter before the injection of contrast medium. Ten milliliters of bile was collected from each patient, and the samples were immediately stored in −80°C freezer for preservation.

DNA Extraction, PCR Amplification and 16S rRNA Sequencing

Microbial DNA was extracted from bile samples using the QIAamp^® DNA Mini Kit (250) (QIAGEN, Germany) in accordance with the protocols. The quality of total DNA was assessed using Thermo NanoDrop 2000 UV microspectrophotometer and 1% agarose gel electrophoresis. We used primers 341F (5’-CCTACGGGRSGCAGCAG-3’) and 806R (5’-GGACTACVVGGGTATCTAATC-3’) to perform PCR amplification targeting the V3-V4 region of 16S rRNA genes. PCR reactions were carried out in a 25μL mixture containing 12.5μL of KFX HiFi 2×PCR Master Mix, 1μL of each primer (10μM), 50ng of template DNA and ddH2O. Amplicons were extracted from 2% agarose gels and purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, U.S,). The concentration of the purified products was determined using the Qubit^®2.0 (Invitrogen, U.S.). Purified amplicons were prepared for library construction and subjected to 250bp paired-end sequencing using the Illumina NovaSeq platform (Illumina, Inc., CA, USA). 16S rRNA sequencing data are available from the National Center for Biotechnology Information BioProject database under accession no. PRJNA1121231.

UPLC–MS/MS Analysis of Bile Acids

The bile samples were thawed at 4 °C. A mixture of methanol and acetonitrile (1:1, v/v) containing internal standard was added to the samples. The mixture was vortexed for 10 min to precipitate proteins, followed by the centrifugation at 14,000 rpm for 10 min. Subsequently, the supernatant was dried under a stream of nitrogen gas. Next, 100 μL of an acetonitrile-water solution (acetonitrile: water =1:1, v/v) was added to reconstitute the samples. The samples were vortexed again and centrifuged at 14,000 rpm for 10 min. Finally, the supernatant was injected into UPLC-MS/MS system (Triple Quad 6500+, SCIEX, the Netherlands) for analysis.

The ion source operated in negative ion mode and used an electric spray ion source. The ion source parameters consisted of 60 psi for desolvent gas, 55 psi for heating gas, 500 °C for desolvent gas, 30 psi for Curtain Gas, 12 psi for Collision Gas and −4500 V for spray voltage. Multiple reaction monitoring mode was employed for scanning. Liquid phase separation was conducted using C18 chromatographic column (Phenomenex, 100×2.1 mm, 1.7 μm). The mobile phase consisted of water and methanol with the addition of 5 mm ammonium acetate. The gradient elution was carried out over a duration of 9 minutes. The gradient elution program was as follows: started at 20%; linearly increased from 20% to 80% and held at 80%; linearly increased from 80% to 100% and held at 100%; linearly decreased from 100% to 20% and finally held at 20%. Blank samples were injected before data acquisition to balance the instrument. The testing sequences were blanks, calibration curves, blanks, quality control (QC) samples, samples and QC samples.

Data Analysis

Demographic characteristics were recorded in detail, including age, sex, body mass index (BMI), stone number, laboratory test results and past medical history. Statistical analysis was conducted using SPSS Statistics 27.0. Continuous variables were presented as mean ± standard deviation (mean±SD) or median and analyzed using Student’s t-test or Mann–Whitney test. Categorical variables were presented as count and percentage and compared using Chi-square test or Fisher’s exact test. And p<0.05 was considered statistically significant.

Upon successful sequencing of the bile samples, the tags whose quality ranged from 250bp to 500bp were immediately examined and filtered. Phred score of bases was better than 30 (Q30) and less than 1 ambiguous N. The redundancy of repeated tags would be removed when the copy number was enumerated. Only the sequences whose frequency was greater than 1 could be clustered into Operational Taxonomic Units (OTUs). OTUs were clustered on a 97% similarity threshold using UPARSE (http://drive5.com/uparse/). Chimeric sequences were removed using Userach (version 7.0). The taxonomic classification of each representative sequence was identified by RDP Classifier (http://rdp.cme.msu.edu/) against the RDP database (http://rdp.cme.msu.edu/) with a confidence threshold of 80%.

OTU profiling table and alpha/beta diversity analyses were performed using QIIME (version 1.9.1). Sequencing depth was assessed through good’s coverage index. Alpha diversity analysis was used to calculate the community richness and diversity. Beta diversity analysis was employed to evaluate microbial composition and distribution between the two groups based on unweighted UniFrac distance. And the significance was determined by principal coordinates analysis (PCoA) and Anosim test. The core microbiome (Venn diagram) was generated using R/Perl SVG. We used linear discriminant analysis (LDA) to conduct LDA Effect Size (LEfSe) analysis, in order to estimate the influence of each component abundance. The LDA threshold was set at 2.0. Spearman correlation heatmap was figured out using corrplot package in R based on the top 21 different genera.

Metabolic pathways were analyzed by PICRUSt2.0 and matched to MetaCyc database to predict the metabolic functions. The data of bile acids were collected by Analyst 1.6.2 software, and the quantitative analysis was carried out using Multiquant software. To investigate the relationship between differential genera and bile acids, the correlation analysis was conducted using cor.test in R, and the heatmap was generated through heatmap.2 package.

Results

Demographics and Clinical Baseline Characteristics

A total of 31 patients with choledocholithiasis were enrolled in this study, 14 in SS group and 17 in MS group. No significant differences were observed in terms of sex distribution (p=0.456), age (SS group vs MS group, 69.00±13.42 vs 68.94±18.07, p=0.992) and BMI (SS group vs MS group, 22.34±2.57 vs 21.69±2.32, p=0.472) (Table 1). Recurrence of CBDS was similar between the two groups. In detail, 3 patients in SS group and 6 patients in MS group suffered from recurrent CBDS (p=0.456). The two groups shared similar past medical history (all p>0.05). The results of preoperative laboratory examinations were comparable between the two groups, including hemoglobin (Hb), white blood cell (WBC), platelet (PLT), total bilirubin (TBIL), direct bilirubin (DBIL), alanine aminotransferase (ALT), aspartate aminotransferase (AST), γ-glutamyl transpeptidase (γ-GGT), alkaline phosphatase (ALP), albumin (ALB) and CA19-9 (all p>0.05) (Table 1).

Table 1 Demographics and Preoperative Characteristics

Lower Richness and Different Composition of Biliary Microbiota in Patients with Multiple CBDS

The average sequence length ranged between 420bp and 440bp. The good’s coverage index was comparable between the two groups (p=0.89) (Figure 1A). And both of the two indexes were approaching to 1.00. The sequencing depth in this study was enough for further analyses.

Figure 1 Diversity analysis for biliary microbial richness and community composition: (A) good’s coverage index was comparable between the two groups (p=0.89), and both approached to 1.00; (B) richness of biliary microbiota in MS group was lower than SS group (p=0.013); (C) PCoA analysis using unweighted UniFrac distance (p=0.023 on PCoA1 axis); (D) the difference between groups was greater than the difference within groups based on Anosim analysis (R=0.139, p=0.013).

The richness of biliary microbial community in MS group was significantly lower than SS group according to alpha diversity analysis (p=0.013) (Figure 1B). For differentiating multiple CBDS from single one, beta diversity analysis was further carried out. The distribution of biliary microbiota was demonstrated to be significantly different between the two groups by PCoA analysis (p=0.023) (Figure 1C). The difference between groups was significantly greater than the difference within groups (p=0.013) based on Anosim analysis using unweighed unifrac distances (Figure 1D).

There were 528 OTUs and 312 OTUs identified in SS group and MS group, respectively, of which 231 OTUs were shared between the two groups (Figure S1). The distribution of biliary microbiota differed between the two groups (Figure S2).

Characteristics of Biliary Microbiota in Patients with Multiple CBDS

In order to describe the details of biliary microbiota, the relative abundance and composition of detected flora were calculated and compared. The numbers of OTUs assigned to different taxonomic levels are shown in Table 2. The distribution and composition of detected bacteria varied from each other at each taxonomic level (Figure 2A and B & Figure S3). Firmicutes and Proteobacteria occupied the dominant niches at phylum level in all the samples (Figure 2C), while Escherichia/Shigella, Enterococcus and Klebsiella were the predominant genera (Figure 2D). We present detailed quantitative data in Table S1.

Table 2 Assignment of OTUs to Different Levels

Figure 2 Species annotation analysis showed the relative abundance and distribution of biliary microbiota at different taxonomic levels: (A) comparison at phylum level; (B) comparison at genus level; (C) microbial community at phylum level in each sample; (D) microbial community at genus level in each sample. (S represents SS group, M represents MS group).

Firmicutes (54.11%) played a dominant role in SS group at phylum level. And Proteobacteria, Bacteroidetes, and Fusobacteria accounted for 32.27%, 5.25% and 4.52%, respectively. In contrast, Proteobacteria (50.41%) was the predominant phylum in MS group, followed by Firmicutes (45.05%), Actinobacteria (2.34%) and Fusobacteria (2.04%). At genus level, the top five genera in MS group were Enterococcus (29.94%), Klebsiella (27.32%), Escherichia/Shigella (20.80%), Veillonella (4.65%), and Abiotrophia (4.33%), while Escherichia/Shigella (27.77%) accounted for the highest proportion in SS group, followed by Enterococcus (15.88%), Streptococcus (15.79%), other unclassified bacteria (13.44%) and Staphylococcus (7.31%). The significant reduction in Streptococcus between the two groups approached to 83%.

Further analyses were conducted to illustrate the differences in biliary microbial structure between the two groups. We found the distinct composition of biliary microbiota and the unique flora in the two groups (Figure 3A). The relative abundance of 54 genera were identified to be significantly different, of which 50 genera were abundant in SS group. While the other 4 genera were enriched in MS group (Figure 3B). In detail, Klebsiella, Aquabacterium, Morganella and Diaphorobacter were illustrated to be significantly enriched in MS group at genus level (Figure 3B). Both Morganella and Aeromonas were exclusively found in MS group based on annotation analysis, along with the absence of Metaprevotella (Table 3). The genera abundant in each group were positively correlated with each other, however, negatively correlated with the genera that enriched in the other group (Figure 3C).

Table 3 Relative Abundance of the Top 20 Genera

Figure 3 Characteristics of biliary microbiota in patients with multiple CBDS: (A) LEfSe analysis revealed the different composition and the unique predominant bacteria; (B) the relative abundance of the 54 significantly different bacteria (LDA score threshold was 2.0); (C) the genera abundant in MS group were positively correlated with each other, however, negatively correlated with the genera that enriched in SS group.

Changes of Biliary Metabolism

Bile acids were extracted from all the bile samples. We conducted quantitative calculation to compare the contents of each bile acid (Table 4). Chenodeoxycholic acid (CDCA) was significantly abundant in MS group (MS group vs SS group, 14.77μmol/mL vs 2.69μmol/mL, p=0.048), while no significant differences were found between the two groups concerning the other 14 bile acids (all p>0.05). An obvious decreasing tendency of several conjugated bile acids was also observed in MS group, including glycodeoxycholic acid (GDCA), glycolithocholic acid (GLCA), taurodeoxycholic acid (TDCA), taurochenodeoxycholic acid (TCDCA) and tauroursodeoxycholic acid (TUDCA). CDCA was illustrated to have negative correlation with Enhydrobacter, Massilia and Neglecta (Figure 4A), which were significantly enriched in SS group (Figure 3B). On the other hand, GDCA and GLCA were both positively correlated with Mogibacterium, Dialister, Prevotella, Leptotrichia, Lancefieldella, Centipeda, Selenomonas and Catonella (Figure 4A), which were all abundant in SS group (Figure 3B).

Table 4 UPLC-MS/MS† Analysis of Bile Acids

Figure 4 Changes of biliary metabolism: (A) correlation analysis between biliary microbiota and bile acids; (B) metabolic pathways were calculated by PICRUSt 2.0 based on the abundance of biliary microbiota to predict functions; (C) PCA analysis indicated that metabolic functions of biliary microbiota in MS group differed from SS group.

We performed metabolic prediction using PICRUSt 2.0 method based on MetaCyc database, to illustrate the potential metabolic functions which increased stone number. The MS and SS groups exhibited enrichment of 21 and 14 metabolic pathways, respectively (Figure 4B). Metabolic functions of biliary microbiota in MS group showed obvious differences from SS group (Figure 4C). The metabolic superpathways enriched in MS group included synthesis and metabolism of cofactors and vitamins, sugar biosynthesis and metabolism, energy metabolism and biosynthesis of secondary metabolites.

Discussion

The development of endoscopic technology has facilitated the microinvasive treatment for cholelithiasis. CBDS can be removed by ERCP for nearly 97% of patients.¹² But the successful extraction is commonly limited by various factors, including unusual shape, incarcerated location, large diameter and multiple stones,⁹ of which the number of stones plays a vital role. It is more time-consuming for endoscopists to manage multiple CBDS through ERCP. And the success rate is lower than single CBDS.^10,11 Based on the advancements in sequencing technology, emerging studies have focused on the microbial community in biliary tract, which was regarded as a sterile space. The relative abundance of biliary microbiota potentially influenced biological and medical characteristics. It is deserved to be researched on the association between biliary microbiota and hepatobiliary diseases, especially cholelithiasis.

Recurrent cholelithiasis¹³ is the most popular topic of studies concerning biliary microbiota. Firmicutes, Bacteroidetes, Proteobacteria and Actinobacteria were dominant phyla in bile, either non-diseases or with cholelithiasis. The biliary microbial distribution of cholelithiasis was significantly distinct from healthy volunteers.⁴ Furthermore, species diversity and composition of microbiota in bile differed from gut. The stones shared 85% of OTUs with bile.¹⁴ Cholesterol aggregation and mucin secretion were considered as key points between biliary microbiota and the formation of cholesterol stones.¹⁵ To the best of our knowledge, no studies focus on how biliary microorganism influences the number of CBDS. Not only the species richness (Figure 1B), but also the composition of biliary microbiota (Figure 1C and D) from the patients with multiple CBDS was significantly different from the patients with single CBDS. Like the previous findings,^4,16 Firmicutes and Proteobacteria were the predominant phyla in all samples (Figure 2C), while the distribution was totally different in the two groups. Compared with our prior work,⁹ Firmicutes was the dominant phylum in the giant CBDS, while Proteobacteria occupied the predominant niche in the multiple CBDS. These two phyla are capable of promoting the formation of CBDS. It seems that Firmicutes could enhance the diameter of stones, while Proteobacteria could increase the number of stones. It deserves to conduct further research to reveal the effect of Firmicutes and Proteobacteria.

In the present study, we found that Enterococcus, Klebsiella, Escherichia/Shigella, Veillonella and Abiotrophia accounted for the majority at genus level in MS group. Taking the prior work⁹ into consideration, four of the top five genera were shared by giant CBDS and multiple CBDS, including Enterococcus, Klebsiella, Escherichia/Shigella and Abiotrophia. Ye et al found a potential link between the invasion of Escherichia coli during ERCP procedure and the formation of stones.¹⁷ Escherichia coli originated from gut was identified to have the ability to cause acute cholecystitis.¹⁸ According to our findings, Escherichia/Shigella occupied a high relative abundance, no matter in which group (Table 3). Enterococcus was illustrated to be involved in the development of pigment stones due to the higher abundance.¹⁹ Both of Enterococcus and Klebsiella have the natural resistance to bile and contribute to the biofilm formation of stones.²⁰ The relative abundance of Enterococcus and Klebsiella were the top two in MS group (Figure 2B). Klebsiella was illustrated to be significantly enriched in MS group, followed by Aquabacterium, Morganella and Diaphorobacter (Figure 3B). Moreover, Klebsiella is the only genus that exists in giant CBDS, multiple CBDS and control individuals. Lactobacillus is able to absorb cholesterol and reduce total serum cholesterol levels.^6,21 The deficiency of Lactobacillus served as a risk factor for cholelithiasis.²² In the present study, the presence of Lactobacillus was only found in four samples, 2 in SS group and 2 in MS group, respectively. Additionally, the significant reduction of Streptococcus may play a critical role in increasing the number of CBDS. The findings indicate the potential role of these species in the formation of CBDS. Further studies are going to focus on the mechanism by which these bacteria contribute to the stone formation and the increase of stone diameter and number.

The risk of cholestasis can be increased by the accumulation of hydrophobic bile acids in hepatobiliary system, which contribute to the formation of CBDS.²³ CDCA is a typical hydrophobic bile acid. In our prior work, the content of CDCA was found significantly higher in the bile of giant CBDS.⁹ In the present study, CDCA was negatively correlated with the genera enriched in SS group. And it was the only bile acid that was significantly abundant in the bile of patients with multiple CBDS (Table 4), while no more significant differences about other bile acids could be observed between the two groups. The change of biliary microbiota altered bile acid metabolism. The increase of CDCA influenced bile cholestasis, leading to the increase in the diameter and number of CBDS.

Based on the calculation using MetaCyc database, the number of CBDS could be influenced by various metabolic functions, including synthesis and metabolism of cofactors and vitamins, sugar biosynthesis and metabolism, energy metabolism and biosynthesis of secondary metabolites. Methionine restriction (MR) is valuable for keeping balance of bile acids circulation, to maintain the stability of bile physiology and reduce the incidence of cholelithiasis.^24,25 In our study, the superpathway of L-methionine biosynthesis (PWY-5347) was found to be significantly enriched in MS group. The increase in the number of CBDS may be related to the enhanced methionine biosynthesis that is associated with biliary microbiota. Moreover, several metabolic pathways concerning cholelithiasis were also significantly enriched in MS group, including aspartate superpathway (PWY0-781), glucose and glucose-1-phosphate degradation (GLUCOSE1PMETAB-PWY) and superpathway of glycolysis and the Entner-Doudoroff pathway (GLYCOLYSIS-E-D). The distinct composition of biliary microbiota altered the metabolism in bile, leading to the increased number of CBDS. It is meaningful to design further experiments to elucidate the mechanisms.

To avoid the potential bias from different regions on microbiota community, the present study was conducted in our single hospital. A large number of patients were initially admitted to primary hospitals due to cholangitis. Antibiotic treatment would be scheduled for patients once the diagnosis was confirmed. However, these patients would be ruled out when they were transferred to our hospital for ERCP. As a result, the sample size was relatively small. Further studies are going to improve referral process for controlling the use of antibiotics. In this study, the metabolic pathways were analyzed by PICRUSt 2.0 and matched to MetaCyc database. Maybe some pathways would be missed because of the focus on MetaCyc. More calculations based on other databases are needed to avoid the limitation in the further studies.

Conclusions

This study revealed the characteristics of biliary microbiota and associated metabolic functions in patients with multiple CBDS. Microbial niche change in biliary tract increased the cholelithiasis burden. Distinct microbial richness and composition can be used to evaluate medical condition of multiple CBDS. The metabolic functions abundant in MS group indicated the underlying mechanism of multiple CBDS formation. The findings innovatively highlighted the importance of biliary microbiota in the etiology of multiple CBDS.