This cross-sectional study was conducted between July 2022 and September 2023 at Bengbu Third People's Hospital following institutional review board approval. A total of 28 adult patients (aged 29–69 years) scheduled for elective cholecystectomy with radiologically confirmed gallbladder stones were prospectively enrolled. The diagnosis of GSD was established through ultrasonographic examination demonstrating hyperechoic structures with acoustic shadowing within the gallbladder lumen, subsequently confirmed by pathological examination of surgical specimens. The exclusion criteria were as follows: (1) current or recent antibiotic therapy within 30 days prior to surgery; (2) presence of acute cholangitis or severe cholecystitis requiring emergency intervention; (3) documented inflammatory bowel disease or other chronic gastrointestinal disorders; (4) active malignancy or immunosuppressive therapy; (5) pregnancy or lactation; or (6) inability to provide informed consent.

Intraoperative bile collection was performed using standardized aseptic techniques to minimize contamination. Following gallbladder exposure, but before dissection from the hepatic bed, the surgeon performed direct aspiration of bile using an 18-gauge needle attached to a sterile syringe. Approximately 10–15 ml of bile was collected from each patient, immediately transferred to sterile cryovials, placed on ice and transported to the laboratory within 30 min of collection. Samples were aliquoted to avoid repeated freeze–thaw cycles and stored at −80 °C until deoxyribonucleic acid (DNA) extraction. Comprehensive clinical data were collected through structured questionnaires and electronic medical record review, including demographic characteristics [age, sex, and body mass index (BMI)], medical history (comorbidities, medication use), laboratory parameters [complete blood count and liver function tests, including alanine aminotransferase (ALT), aspartate aminotransferase (AST), gamma-glutamyl transferase (GGT), bilirubin levels, fasting glucose, C-reactive protein (CRP), total bile acids], imaging findings (gallbladder wall thickness, stone number and size, presence of fatty liver). The pathological examination of resected gallbladders was performed by experienced pathologists blinded to the microbiome analysis results.

Fasting venous blood samples were collected within 24 h prior to surgery following an overnight fast of ≥12 h. Serum lipid profiles were determined using automated biochemical analysers (Roche Cobas 8000 series; Roche Diagnostics GmbH, Mannheim, Baden–Württemberg, Germany; authorized local distribution in China: Roche Diagnostics Products (Shanghai) Co., Ltd., Shanghai, Shanghai Municipality, China) with standardized enzymatic methods. Total cholesterol (TC) was measured using the cholesterol oxidase–peroxidase method, TGs by the glycerol phosphate oxidase–peroxidase method, high-density lipoprotein cholesterol (HDL-C) by direct homogeneous enzymatic assay and low-density lipoprotein cholesterol (LDL-C) was calculated using the Friedewald equation when TGs were < 400 mg/dl or measured directly when TGs exceeded this threshold.

For the primary analysis, patients were stratified into two groups based on TG levels according to contemporary dyslipidaemia guidelines: (1) a hypertriglyceridaemia (HTG) group, defined as fasting TGs ≥1.7 mmol/L (150 mg/dl); and (2) a non-HTG (NTG) group, with TGs < 1.7 mmol/L. This binary classification was chosen to focus on the relationship between elevated TGs—a key component of metabolic syndrome strongly associated with gallstone formation—and bile microbiome composition. The TG threshold of 1.7 mmol/L represents the internationally accepted cutoff for defining HTG and increased cardiovascular risk. This stratification approach served as the primary grouping variable for all subsequent microbiome comparative analyses.

Microbial DNA extraction from bile samples was performed using a modified protocol optimized for low-biomass samples with high concentrations of bile salts and other inhibitory compounds. Extraction blanks (no sample input) were included as negative controls to monitor potential contamination during DNA extraction and library preparation. Bile samples (200 μl) were first treated with mutanolysin (25 U/ml) and lysozyme (20 mg/ml) for 1 h at 37 °C to enhance bacterial cell wall disruption. Subsequently, mechanical lysis was achieved through bead-beating using 0.1 mm zirconia/silica beads (BioSpec Products) in a FastPrep-24 homogeniser (MP Biomedicals) at 6.0 m/s for 60 s, repeated twice with 5-min cooling intervals on ice. Extraction of DNA then proceeded using the QIAamp DNA Mini Kit (Qiagen) with modifications including extended proteinase K digestion (2 h at 56 °C) and addition of carrier RNA to improve DNA recovery from low-biomass samples.

The V3–V4 hypervariable region of the bacterial 16S rRNA gene was targeted for amplification using universal primers 341F 5′-CCTAYGGGRBGCASCAG3′) and 806R 5′-GGACTACNNGGGTATCTAAT3′) containing Illumina adapter sequences and sample-specific barcodes. Polymerase chain reaction (PCR) amplification was performed in triplicate 25 μl reactions containing 12.5 μl 2 × KAPA HiFi HotStart ReadyMix, 0.5 μM of each primer and 5 ng template DNA. Thermal cycling conditions consisted of initial denaturation at 95 °C for 3 min, followed by 30 cycles of 95 °C for 30 s, 55 °C for 30 s and 72 °C for 30 s, with final extension at 72 °C for 5 min. Triplicate PCR products were pooled, purified using AMPure XP beads (Beckman Coulter) and quantified using Qubit fluorometry and quantitative PCR (qPCR).

Library preparation followed the NEBNext Ultra II DNA Library Prep Kit protocol with size selection targeting 450–550 bp fragments. Final libraries were validated using the Agilent 2100 Bioanalyzer for size distribution and molarity calculation. Equimolar pooling of indexed libraries was performed based on qPCR quantification, with a 20% PhiX control spike-in to increase sequence diversity. Paired-end sequencing (2 × 250 bp) was conducted on the Illumina MiSeq platform following the manufacturer's protocols, generating approximately 100,000 reads per sample to ensure adequate coverage for rare taxa detection.

The bioinformatics pipeline integrated multiple quality control and analysis steps to ensure robust taxonomic profiling and diversity assessment. (1) Quality control and read processing began with raw sequence demultiplexing using bcl2fastq, followed by adapter and primer removal using cutadapt v3.5 with error tolerance of 0.15. Quality filtering was performed using fastp v0.23 with parameters including minimum quality score of Q20, minimum read length of 200 bp and sliding window trimming (window size = 4, quality threshold = 20). Paired-end reads were merged using USEARCH v11 with a minimum overlap of 16 bp and a maximum mismatch ratio of 0.1. (2) Amplicon sequence variant (ASV) inference employed the UNOISE3 algorithm for denoising, which simultaneously performs error correction and chimera removal. This approach generates exact sequence variants providing single-nucleotide resolution superior to traditional operational taxonomic unit (OTU) clustering. Additional chimera filtering was performed using UCHIME in reference mode against the SILVA database. Sequences identified as chloroplast, mitochondrial or unclassified at the kingdom level were removed from further analysis. Data preprocessing steps were implemented in accordance with established best practices to minimize bias and ensure data quality, as highlighted by He et al. (2020). (3) Taxonomic assignment utilized the SILVA v138.1 reference database with the Ribosomal Database Project Classifier algorithm at a 0.8 confidence threshold. This threshold balanced sensitivity and specificity for bile samples, which often contain taxa absent from reference databases. The phylogenetic tree construction for the unique fraction metric (UniFrac) analyses employed MUSCLE alignment followed by FastTree. (4) Diversity metrics were calculated after rarefaction to 10,000 sequences per sample to standardize the sequencing depth. All 28 samples were retained after rarefaction, with no samples excluded due to insufficient read depth. The alpha diversity indices included observed ASVs (richness), Chao1 (estimated richness), Shannon entropy (evenness and richness), Simpson index (dominance) and Pielou's evenness. Beta diversity was assessed using multiple distance metrics, including Bray–Curtis (abundance-based), Jaccard (presence/absence) and weighted/unweighted UniFrac (phylogenetically-informed). Ordination techniques included non-metric multidimensional scaling (NMDS) with stress assessment, principal component analysis (PCA) and principal coordinates analysis (PCoA). (5) Differential abundance testing employed multiple complementary approaches. Non-parametric tests (Wilcoxon rank-sum for pairwise comparisons, Kruskal–Wallis for multi-group) were applied with the Benjamini–Hochberg false discovery rate (FDR) correction at q < 0.05. The linear discriminant analysis effect size (LEfSe) technique was performed using the default parameters in the Galaxy platform. The LEfSe identified biomarkers with a linear discriminant analysis score threshold of 2.0, incorporating both statistical significance (Kruskal–Wallis, p < 0.05) and biological relevance (effect size). (6) Functional prediction utilized PICRUSt2 for inferring metabolic pathways from 16S data, mapping to the Kyoto Encyclopedia of Genes and Genomes and MetaCyc databases. Random forest classification was implemented using 1,000 trees with fivefold cross-validation for feature importance ranking. The random forest models were constructed using the randomForest package in R, with model performance evaluated using precision, recall and area under the receiver operating characteristic curve (AUC-ROC) in addition to the out-of-bag (OOB) error rate. Feature importance was quantified using two metrics: mean decrease accuracy (reduction in classification accuracy when the feature is permuted) and mean decrease Gini (reduction in Gini impurity when the feature is used for splitting nodes in the decision trees), both of which reflect the contribution of each OTU in distinguishing between the HTG and NTG groups.

Statistical analyses were performed using R v4.2.0 with specialized packages for microbiome analysis, including the phyloseq, vegan and microbiome packages. Continuous variables were assessed for normality using Shapiro–Wilk tests and presented as mean ± standard deviation for normally distributed data or median with interquartile range (IQR) for non-parametric data. Categorical variables were expressed as frequencies and percentages. Between-group comparisons employed Student's t-test or Mann–Whitney U test for continuous variables and the chi-squared or Fisher's exact test for categorical variables, as appropriate.

For microbiome-specific analyses, community structure differences were evaluated using permutational multivariate analysis of variance (PERMANOVA/Adonis) with 9,999 permutations and analysis of similarities (ANOSIM) to test for significant clustering by lipid phenotype. The PERMANOVA was performed using the adonis2 function in the vegan package, with Bray–Curtis distances as the input. The dispersion homogeneity was verified using PERMDISP to ensure the PERMANOVA assumptions were met. Multiple testing correction was applied throughout using the Benjamini–Hochberg method with an FDR threshold of 0.05. Effect sizes were calculated using Cohen's d for parametric comparisons and Cliff's delta for non-parametric tests.

A post-hoc power assessment indicated that the sample size of 28 patients provided 75% power to detect a moderate effect size (Cohen's d = 0.8) in the Simpson diversity index (the primary diversity outcome) between HTG and NTG groups at alpha = 0.05.

Given the small sample size (n = 28) and potential for overfitting, the model performance was evaluated using fivefold cross-validation, and the overfitting risk was assessed by comparing training and cross-validation error rates. No evidence of significant overfitting was observed (training error rate: 21.4%; cross-validation error rate: 32.1%). Additionally, the study population size is relatively modest, which may limit generalisability to larger cohorts; however, this sample size is consistent with similar low-biomass microbiome studies in biliary samples (Grigor'eva and Romanova, 2020; Fang et al., 2024).

The study cohort exhibited a significant gender imbalance (85.7% women), which was acknowledged as a potential source of bias. This gender distribution aligns with epidemiological patterns of GSD (women have 2–3 times higher prevalence than men) but may limit the generalisability of findings to male populations. Sensitivity analyses excluding male patients (n = 4) did not substantially alter the primary findings, suggesting that gender imbalance did not majorly confound the observed microbiome–lipid associations.

This study was conducted in accordance with the Declaration of Helsinki and received approval from the institutional review boards of all participating centers. Written informed consent was obtained from all participants prior to enrolment, with specific consent for bile sample collection and genetic/microbiome analyses. Patient confidentiality was maintained through de-identification of all samples and data, with unique study identifiers assigned to each participant. Biospecimens were handled according to international guidelines for human biological material, with appropriate biosafety measures for potentially infectious samples.

The study cohort comprised 28 patients with confirmed GSD and demonstrated a predominance of women (85.7%, n = 24); this is consistent with established epidemiological patterns, where women have 2–3 times higher gallstone prevalence than men, attributed to estrogen effects on cholesterol metabolism and gallbladder motility. The median age was 53.5 years (IQR 40.75–58.75, range 29–69), with a mean BMI of 24.98 ± 3.24 kg/m² (median 25.05, IQR 22.39–27.22, range 20.0–30.86), indicating a predominantly overweight population according to World Health Organization criteria (normal: 18.5–24.9 kg/m²; overweight: 25.0–29.9 kg/m²; obese: ≥30.0 kg/m²). Based on TG levels, 11 patients (39.3%) were classified into the HTG group with TGs ≥1.7 mmol/L, and 17 patients (60.7%) comprised the NTG group. The HTG group showed significantly higher median TGs (2.53 mmol/L, IQR 2.00–3.85) than the NTG group (1.09 mmol/L, IQR 0.85–1.24, p < 0.001). Other lipid parameters, including TC (HTG: 5.15 ± 1.02 vs. NTG: 4.87 ± 0.91 mmol/L, p = 0.472) and LDL-C (HTG: 2.98 ± 0.99 vs. NTG: 2.95 ± 0.94 mmol/L, p = 0.938) were comparable between groups, whereas HDL-C was significantly lower in the HTG group (0.96 ± 0.21 vs. 1.31 ± 0.27 mmol/L, p = 0.001).

Inflammatory markers showed modest elevation with a median CRP of 4.75 mg/L (IQR 2.08–5.00, range 0.90–38.20), with sample BM14 from the HTG group showing the highest CRP level at 38.20 mg/L, suggesting predominantly chronic rather than acute inflammatory processes. Liver function tests demonstrated mild abnormalities, with a mean GGT of 63.86 ± 84.15 U/L (median 24.50, range 9.0–336.0), whereas transaminases remained within normal ranges (ALT 32.38 ± 25.95 U/L, AST 26.73 ± 11.69 U/L). Ultrasonographic findings most commonly reported gallbladder stones (25.0%), multiple gallstones (14.3%) and gallstones with fatty liver (10.7%). Pathological examination confirmed chronic cholecystitis with gallstones in 71.4% of cases, chronic cholecystitis with stone formation in 14.3% and chronic cholecystitis with stones and mild glandular atypical hyperplasia in 3.6%. The complete baseline characteristics are presented in Table 1.

Taxonomic profiling revealed distinct bile microbiome compositions between the HTG and NTG groups. The stacked bar plot analysis (Figure 1) demonstrated compositional differences between groups, with visible variations in bacterial taxa distribution. A detailed OTU analysis identified 981 total OTUs, with the Venn diagram analysis (Figure 2) revealing distinct microbiome signatures between groups; the HTG group harbored 137 unique OTUs, the NTG group contained 467 unique OTUs and 357 OTUs were shared between groups. This distribution suggests that HTG is associated with a more restricted bile microbiome, with total rare OTUs representing 743 or 75.7% of all OTUs, including 89 HTG-specific rare OTUs, 412 NTG-specific rare OTUs and 242 shared rare OTUs.

At the genus level, the dominant taxa showed distinct abundance patterns. Pseudoalteromonas (OTU_1) represented 17.79% in HTG vs. 26.22% in NTG, and Vibrio (OTU_2) comprised 21.35% in HTG vs. 17.79% in NTG. Burkholderia–Caballeronia–Paraburkholderia (OTU_3) showed 9.58% in HTG vs. 12.58% in NTG, and Alcaligenaceae (OTU_5) demonstrated 6.21% in HTG vs. 3.07% in NTG. Notably, Lysinibacillus (OTU_6) exhibited a 4.9-fold enrichment in HTG with 5.48% vs. 1.12% in NTG. Additional taxa showed trends toward differential abundance, including Allorhizobium–Neorhizobium-Pararhizobium–Rhizobium (OTU_34) with 0.417% in HTG vs. 0.008% in NTG (a 52-fold increase from a low baseline), Streptococcus (OTU_37) showing 0.384% in HTG vs. 0.004% in NTG (a 96-fold increase from a low baseline) and Brevundimonas (OTU_35) demonstrating 0.343% in HTG vs. 0.018% in NTG (a 19-fold increase from a low baseline). Conversely, NTG-enriched taxa included Bacillaceae (OTU_150) with 0.024% in HTG vs. 0.281% in NTG (an 11.7-fold increase), Blautia (OTU_32) showing 0.006% in HTG vs. 0.297% in NTG (a 49.5-fold increase) and Gramella (OTU_40) demonstrating 0.0001% in HTG vs. 0.201% in NTG (a 2,010-fold increase from an extremely low baseline; Figure 3).

Rarefaction curve analysis (Figure 4) demonstrated adequate sequencing depth for capturing microbial diversity, with curves reaching asymptotic plateaus at approximately 10,000 sequences per sample. Alpha diversity metrics revealed largely overlapping distributions between groups, with Simpson diversity index analysis showing median values of 0.209 (IQR 0.168–0.311) for the HTG group and 0.312 (IQR 0.307–0.332) for the NTG group. Individual sample Simpson indices ranged from 0.166 to 0.355, with sample BM15 showing the highest diversity at 0.355 and BM7 showing the lowest at 0.166. Shannon diversity revealed an HTG median of 3.21 (range 2.98–3.45) vs. an NTG median of 3.28 (range 3.05–3.51; Figure 5). The observed OTUs at plateau showed an HTG mean of 287 ± 45 vs. an NTG mean of 342 ± 68, and the Chao1 richness estimator indicated an HTG of 312 ± 52 vs. an NTG of 378 ± 71. Pielou's evenness demonstrated an HTG of 0.71 ± 0.08 vs. an NTG of 0.74 ± 0.06, with sample BM15 showing maximum observed richness at 412 OTUs.

Beta diversity analysis revealed patterns of microbial community differentiation between the HTG and NTG groups. The NMDS ordination (Figure 6) achieved a stress value of < 0.2, indicating acceptable representation of community distances, with some separation between HTG and NTG samples, though with overlap. The PCA (Figure 7) showed PC1 accounting for 11.8% and PC2 for 10.9% of total variance, with partial separation between HTG and NTG samples along the principal components. The PCoA using weighted UniFrac distances (Figure 8) revealed that PC1 explained 71.3% of variance and PC2 explained 17.2% of variance, with cumulative variance of 88.5% in the first two components, demonstrating phylogenetically-informed clustering, with HTG and NTG samples showing distinct but overlapping distributions.

Detailed distance matrix analysis revealed maximum dissimilarity between any two samples of 1.414 between BM1 and BM5, minimum dissimilarity within the HTG group of 0.022 between BM19 and BM20 and minimum dissimilarity within the NTG group of 0.012 between BM21 and BM23. The average between-group distance was 0.812 ± 0.094, the average within-HTG distance was 0.456 ± 0.187 and the average within-NTG distance was 0.398 ± 0.156. The Bray–Curtis dissimilarity analysis revealed inter-group distances ranging from 0.776 to 0.836 between the HTG and NTG samples, whereas intra-group distances showed tighter clustering with HTG ranging from 0.783 to 0.827 and NTG ranging from 0.789 to 0.823. The distance heatmap (Figure 9) demonstrated clustering patterns, with red indicating high similarity and blue indicating dissimilarity between samples. Hierarchical clustering analysis (Figure 10) revealed sample groupings that partially corresponded to the HTG and NTG groups. The PERMANOVA results showed that R² = 0.068 with p = 0.083, indicating marginal significance, and ANOSIM demonstrated R² = 0.041 with p = 0.091.

Two-group comparisons (Figure 11) revealed OTUs with varying abundance levels between groups. Statistical analysis revealed that although most OTUs showed non-significant differences after FDR correction with q-values >0.05, several taxa showed strong trends with uncorrected p-values < 0.05 for differences between groups. Inter-OTU correlations revealed ecological networks, including a strong positive correlation >0.7 between OTU_1 and OTU_814 with r = 0.82, a strong negative correlation between OTU_6 and OTU_11 with r = −0.71 and hub taxa with >10 significant correlations, including OTU_1, OTU_2 and OTU_3.

The random forest classification (Figure 12) achieved moderate accuracy, with an OOB error rate of 32.1%, class-specific errors of 36.4% for HTG and 29.4% for NTG, precision of 0.68, recall of 0.64 and an AUC-ROC of 0.72. The model identified top predictor OTUs with feature importance rankings shown by both mean decrease accuracy and mean decrease Gini metrics. The top five predictor OTUs by mean decrease Gini were OTU_37 with Gini importance 0.171, OTU_43 with 0.101, OTU_18 with 0.097, OTU_56 with 0.094 and OTU_38 with 0.094. By mean decrease accuracy, the top predictors included OTU_38 with score 3.20, OTU_43 with 2.68, OTU_69 with 2.44, OTU_170 with 2.35 and OTU_37 with 2.33. The NTG-enriched predictors included OTU_11 (Maribacter) with 1.60% in NTG vs. 0.30% in HTG, OTU_17 (Wolbachia) with 0.92% in NTG vs. 0.02% in HTG and OTU_24 (Sphingobacterium) with 0.51% in NTG vs. 0.003% in HTG.

Sample-specific patterns revealed distinct heterogeneity, with sample BM6 showing the highest Proteobacteria at 84.2%, sample BM13 demonstrating the highest Firmicutes at 32.1%, sample BM5 exhibiting a unique OTU profile with 47 sample-specific OTUs and sample BM15 showing maximum observed richness with 412 OTUs.