Patients who have been followed-up for at least 2 years after index CBD stone removal and had at least one of the high-risk factors for CBD stone recurrence were included in this study. High-risk factors for CBD stone recurrence were as follows: maximum CBD diameter > 15 mm, CBD angle ≤ 145°, and the presence of periampullary diverticulum (PAD) [12]. Patients included in this study had undergone cholecystectomy at the time of index ERCP since they presented with both GB and CBD stones. These patients were divided up into two groups according to the history of CBD stone recurrence. Patients who belonged to the recurrence group were those who had at least 2 episodes of recurrences during the follow-up period. The non-recurrence group was defined as those who had no recurrence event during the follow-up period. A recurred CBD stone was defined as a particle with a diameter greater than 3 mm which took on the shape of a stone. CBD stones that were found and removed within 6 months after index ERCP were not considered as recurrent stones but as residual stones due to incomplete CBD stone clearance. Only those detected after 6 months following endoscopic CBD stone removal were regarded as recurrent stones [3]. Of the aforementioned patients, 20 who agreed to undergo surveillance endoscopic retrograde cholangiopancreatography (ERCP) procedures to check for CBD stone recurrence were enrolled. All patients were asymptomatic at the time of surveillance ERCP without clinical evidence of biliary tract infection, and had not taken antibiotics nor proton pump inhibitors for at least 3 months prior to bile sample collection.

From these 20 subjects, a total of 20 bile samples (11 from the CBD stone recurrence group and 9 from the non-recurrence group) were collected during scheduled surveillance ERCP from September 2018 to August 2019. After advancing the side-viewing endoscopes (TJF240/JF-260V; Olympus Optical, Tokyo, Japan) into the duodenum 2nd portion and facing the papilla enface, the cannulation catheter was passed through the working channel and inserted into the CBD through the manipulated ampulla. Two to five milliliters of the patients’ bile aspirates from bile duct were sucked out via the catheter before the injection of contrast agent for ERCP procedure. A portion of the sample was stored at −80 °C for NGS analysis and the remainder was stored at room temperature for analysis of bile juice components. Afterwards, dye was injected to evaluate for the presence of CBD stones and removed when present. CBD stones were broadly classified as being cholesterol or pigment stones (black or brown) based on their characteristic external appearance by gross inspection [13], and all recurred CBD stones were brown pigment stones.

This study was approved by the Institutional Review Board of Korea University Ansan Hospital (No. 2018AS0208) and the research was conducted in accordance with the Declaration of Helsinki. All patients provided written informed consent upon enrollment.

To analyze the bacterial community, the V3/V4 regions of the 16S rRNA gene of DNA extracted from the bile samples were amplified, sequenced, and analyzed. The extraction method for bacterial DNA was performed by using a PowerMax Soil DNA Isolation Kit (MO BIO). Each sequenced sample was prepared according to the Illumina 16S Metagenomic Sequencing Library protocols to amplify the V3 and V4 regions (519F-806R). The DNA quality was measured by PicoGreen and NanoDrop. Input genomic DNA (10 ng) was PCR amplified. The barcoded fusion primer sequences used for amplifications were as follows: 341F: 5′ CCTACGGGNGGCWGCAG 3′, and 806R: 5′ GACTACHVGGGTATCTAATCC 3′ containing forward overhang adapter pair. The final purified product was then quantified using qPCR according to the qPCR Quantification Protocol Guide (KAPA Library Quantification kits for Illumina Sequencing platforms) and qualified using the LabChip GX HT DNA High Sensitivity Kit (PerkinElmer, MA, USA). Then, the paired-end (2 × 300 bp) sequencing was performed by the Macrogen using the MiSeq™ platform (Illumina, San Diego, CA, USA).

The operational taxonomic units (OTUs) were generated to calculate the number of microbial species in the bile. In order to compare the bacterial community richness of bile between the two groups, observed OTU-based analyses (OTU richness and Chao 1) were performed. OTU richness is the number of bacterial species observed in each individual sample. Chao 1 represents species richness estimators of the expected OTUs present in a group. As Chao 1 measures OTUs expected in samples, given all the bacterial species that were identified in the samples, Chao1 gives more weight to the low abundance species. To determine α-diversity, the Shannon diversity index and the inverse Simpson index were used. The Shannon diversity index is a quantitative measure that shows how evenly the basic individuals are distributed by taking into account and reflecting the number of different species in each sample; a high Shannon diversity index signifies evenness, whereas a low Shannon diversity index denotes unevenness. The Simpson index measures the degree of concentration or dominance of certain species in a sample. Therefore, the inverse Simpson index is an indicator of how evenly the species are distributed; a high inverse Simpson index implies evenness or a lack of dominance, and a low inverse Simpson index suggests unevenness or the presence of dominance. Evaluation of β-diversity, which analyzes community similarity, was performed by calculating pairwise distances using the phylogenetic metric UniFrac.

The bile samples were analyzed for composition of bile, including electrolyte, total bilirubin, total bile acid, cholesterol, phosphorus and Ca2+. Total bilirubin concentration was measured in aliquot on a Cobas 8000 C702 (Roche Diagnostics System, Switzerland). pH, sodium, potassium, chloride, phosphorus, and calcium estimations were carried out on a cobas 8000 ISE (Roche Diagnostics System, Switzerland) using an ion-selective electrode method. For measuring cholesterol and triglycerides, bile was diluted one to five times with methanol and the protein precipitate was removed by centrifugation. The methanol was then evaporated and the residue was taken for analysis. The cholesterol and triglycerides were analyzed by photometric analysis on a Cobas 8000 C702 (Roche Diagnostics System, Switzerland). Total bile acids were also measured with an enzymatic method on a Roche cobas 8000 C702 analyser.

Statistical analyses were performed using IBM SPSS Statistics version 20.0 (IBM, Armonk, NY, USA). The continuous variables were compared using the Student’s t-test. Fisher’s exact test was used to analyze the categorical variables. The sequences generated from pyrosequencing were analyzed with variable software for pre-processing (quality-adjustment, barcode split), identification of OTUs, taxonomic assignment, community comparison, and statistical analysis. The FASTP program was used to remove the adapter sequence, and error-correction was performed for areas where the two reads overlap. The paired-end data, separated by each sample, was assembled in one sequence using FLASH (v. 1.2.11). To ensure that any subsequent analysis was highly accurate, sequences shorter than 400 bp were discarded. By using CD-HIT-OTU, after removing the low-quality sequence, ambiguous reads, chimera reads, etc., which are based on errors, the filtered reads were clustered by identity as OTUs at 97% similarity. The representative sequence of each OTU was performed by BALSTN (v. 2.4.0) to the NCBI 16S Microbial DB, performing a taxonomic assignment with the organism information of the most similar subject. The taxonomy is not assigned if the query coverage and identity score matched to the reference are less than 85%. Community diversity was estimated by using the Shannon and inversed Simpson indexes. The weighted UniFrac distance method was used to perform a principal coordinates analysis (PCoA), and trees were built by the unweighted-pair group method with arithmetic mean (UPGMA).

The mean age of the patients was 62.3 ± 16.6 years (range, 29 to 83 years) and 9 (45.0%) patients were female. The anatomical risk factors for recurrent CBD stones, including CBD diameter, presence of PAD, and angulation of the bile duct ≤ 145°, were not statistically different between the two groups (Table 1).

A total of 303 bacterial species were identified from 20 bile samples. The average OTU, obtained by calculating the mean values after counting the number of observed OTUs in each sample, was reduced in the 11 individuals with CBD stone recurrence compared to the 9 with non-recurrence (29.6 ± 3.6 vs. 60.4 ± 11.6, p = 0.0005, Wilcoxon rank-sum test; Figure 1a). The expected total number of OTUs for each group, as estimated by the Chao1 estimator, was also considerably lower in the recurrence group than in the non-recurrence group (29.8 ± 3.5 vs. 60.5 ± 11.6, p = 0.005, Wilcoxon rank-sum test; Figure 1b).

The α-diversities of bile samples, as measured by the Shannon diversity index and the inverse Simpson index, were significantly lower in the recurrence group than those of the non-recurrence group samples. The Shannon index, which reflects equal microbial proportional abundance in a species, was lower in the recurrence group compared to the non-recurrence group (0.65 ± 0.36 vs. 3.12 ± 0.45, p = 0.027, Figure 2a), signifying that microbial proportional abundance was uneven in the recurrence group. The inverse Simpson index, which is another way of quantifying the evenness of a species, was also lower in the recurrence group compared to the non-recurrence group (0.46 ± 0.10 vs. 0.75 ± 0.06, p = 0.036, Figure 2b), meaning that certain species were more dominant than other species in the recurrence group. According to the Shannon index and inverse Simpson index, the diversity index was statistically higher in the bile of the non-recurrence group than that obtained in the patients with recurrent CBD stones.

The unweighted UniFrac distance metric, which measures the similarity in the microbial communities of the bile samples (β-diversity), revealed an overall microbial composition difference between the two groups (unweighted UniFrac, PERMANOVAR, pseudo-F: 1.655, p = 0.05, Figure 3). This result indicates that the microbiome clustering in the recurrence group was different from that of the non-recurrence group.

At the phylum level, Bacteroidetes spp. were significantly decreased in the recurrence group, whereas Actinobacteria and Firmicutes spp. were over-represented in the recurrence group relative to the non-recurrence group (Table 2). At the genus level, three bacterial taxa displayed significantly different abundance between the recurrence group and the non-recurrence group. The three genera, Neisseria, Capnocytophaga, and Gemella, were enriched in the non-recurrence group. Conversely, E. coli, Enterococcus spp., Klebsiella spp., Acinetobacter spp., Streptococcus spp., and Staphylococcus, which exert β-glucuronidase activity, were more frequent in the recurrence group, even though there was no statistical difference (Table 3).

Bile acid concentration was higher in the non-recurrence group than in the recurrence group (254.51 ± 82.0 mmol/L vs. 147.1 ± 64.2 mmol/L, p < 0.01). Neither the pH values nor Ca²⁺ showed statistical differences between the bile samples of both groups (Table 4).