We prospectively and randomly enrolled 38 patients with HCC, including 18 with stage I–II (HCC12) and 20 with stage III–IV disease (HCC34), as well as 19 patients with liver cirrhosis (LC) and 19 healthy adults (CG), who presented to the First Hospital of Shanxi Medical University between September 2023 and July 2024 for diagnosis and/or clinical management. Fresh fecal samples from patients with liver cirrhosis and HCC were collected on the first day of admission, prior to any therapeutic intervention. Samples were rapidly frozen in liquid nitrogen for 15 min, and subsequently stored at −80 °C. Clinical features, including medical history, laboratory findings, clinical presentation, and disease classification, were retrieved from electronic medical records and the laboratory information system. The diagnosis of liver cirrhosis was established based on hematological tests combined with imaging or histopathology. HCC diagnosis was confirmed using AFP, AFP-L3, and Pivka-II in combination with at least two imaging modalities or histopathological examination. HCC staging followed the latest China liver cancer staging (CNLC) system (Zhou et al., 2025), which incorporates performance status scoring to comprehensively assess treatment tolerance. Staging information for enrolled HCC patients is summarized in Table 1.

All participants had no constipation, hematochezia, diarrhea, dyspepsia, or other gastrointestinal symptoms. None had taken antibiotics or acid suppressants within the preceding month. None of the patients with HCC had undergone surgical resection, local ablation, transarterial therapies, systemic therapy, or immunotherapy at the time of sample collection. No patient had a history of systemic malignancies other than HCC.

Genomic DNA was extracted from fecal samples via the Magnetic Soil and Stool DNA Kit (TianGen, China; Catalog No. DP712) following the manufacturer’s instructions. DNA quality was assessed using 1% agarose gel electrophoresis, and qualified samples were diluted with nuclease-free sterile water to a final concentration of 1 ng/μL. The V3–V4 hypervariable regions of the 16S rRNA gene were amplified using primers 341F (CCTAYGGGRBGCASCAG) and 806R (GGACTACNNGGGTATCTAAT). Each PCR reaction contained 15 μL of Phusion^® High-Fidelity PCR Master Mix (New England Biolabs), 0.2 μM of each primer, and 10 ng of genomic DNA. The thermal cycling protocol consisted of an initial denaturation at 98 °C for 1 min; 30 cycles of 98 °C for 10 s, 50 °C for 30 s, and 72 °C for 30 s; followed by a final extension at 72 °C for 5 min. PCR amplicons were purified using magnetic beads, and target fragments were recovered via the Universal DNA Purification Kit (TianGen, China; Catalog No. DP214). Library preparation was performed via the NEBNext^® Ultra^™ II FS DNA PCR-Free Library Prep Kit (New England Biolabs, United States; Catalog No. E7430L). Libraries were quantified via Qubit 2.0 fluorometry and qPCR prior to sequencing on the NovaSeq 6000 platform with a paired-end 250 bp strategy.

Raw sequencing data were processed by merging paired-end reads, performing stringent quality filtering, and removing chimeric sequences to obtain high-quality effective tags. The DADA2 module in quantitative insights into microbial ecology 2 (QIIME2, v2022.02) was used for denoising to generate amplicon sequence variants (ASVs) and a feature table. Taxonomic annotation was conducted via the classify-sklearn algorithm in QIIME2 against the SILVA 138.1 reference database.

Multiple sequence alignment of all ASV representative sequences was conducted in QIIME2 to infer phylogenetic relationships. To minimize sequencing depth bias, all samples were rarefied to the minimum sequencing depth observed across the dataset. Alpha diversity indices were calculated via QIIME2 and visualized with R (v4.0.3). Beta diversity was assessed based on weighted and unweighted UniFrac distances. Beta diversity heatmaps, and non-metric multi-dimensional scaling (NMDS) analyses were performed via QIIME2, R (v4.0.3), and Perl (v5.26.2). Microbial functional prediction was performed via Tax4Fun (v0.3.1).

Based on the top 15 microbial families and genera identified by the Kruskal–Wallis test in the microbiome analysis, together with clinical features from 38 HCC patients collected at our center, we randomly divided the dataset into a training set and a test set at a 7:3 ratio. Using the scikit-learn package in Python (v3.5.0), we constructed three ensemble machine learning models: random forest (RF), gradient boosting decision tree (GBDT), and extreme gradient boosting (XGB). Model performance was evaluated across multiple metrics, including accuracy, recall, F₁-score, Matthews correlation coefficient (MCC), and area under the ROC curve (AUC), to identify the optimal classifier. For the final selected model, we further assessed stability and generalizability using 5-fold cross-validation (CV) and 200 bootstrap resampling iterations. SHAP analysis was then applied to quantify the contribution of each feature to the model’s classification output and to identify key predictive features based on importance ranking. Finally, external validation was performed using 74 published HCC cases (Zhang et al., 2021), with additional evaluation using 10-fold CV and 200 bootstrap iterations. Because Enterococcus and Bacilli were not directly available in the external dataset, they were substituted with Enterococcaceae and Lachnospiraceae, respectively. This substitution does not imply strict taxonomic equivalence and was applied to enable approximate feature alignment across datasets. In addition, and missing values for three features (Pivka II, Pseudomonas, and Moraxellaceae) were imputed using the mean values from the training set.

Clinical data were analyzed via the t-test, Kolmogorov–Smirnov test, analysis of variance (ANOVA), Welch’s ANOVA, or the Kruskal–Wallis test, as appropriate. No clinical data were missing. Differences in gut microbial communities were evaluated via zero-inflated negative binomial regression (ZINB), Kruskal–Wallis rank-sum tests, and linear discriminant analysis effect size (LEfSe, v1.1.01). All analyses and visualizations were performed via QIIME2 (v2022.02), Perl (v5.26.2), Python, and R (v3.4.3). A p-value <0.05 was considered statistically significant, and multiple comparisons were adjusted via the false discovery rate (FDR) according to the Benjamini–Hochberg procedure.

A total of 76 participants were ultimately recruited for this study. The clinical features of all participants are summarized in Table 2. Males predominated in all groups, consistent with the known epidemiological features of cirrhosis and HCC. To ensure that the healthy participants were free of any disease, their mean age was significantly lower than that of the other groups, averaging 41.37 years. Alanine aminotransferase (ALT) levels in the HCC12 group were higher than those in the LC group, whereas aspartate aminotransferase (AST) levels did not differ among the HCC12, HCC34, and LC groups. Compared with the HCC12 group, the HCC34 group exhibited significantly lower albumin (ALB) and cholinesterase (ChE) levels, along with significantly higher total bilirubin (TBil) levels and Child–Pugh scores. Furthermore, the lg(Pivka II) values were significantly higher in the HCC34 group than in the HCC12 group.

Venn diagram analysis based on 16S rRNA amplicon sequencing revealed 624 and 727 unique ASVs in the HCC12 and HCC34 groups, respectively (Figure 1A). Gut microbial diversity appeared to play a role in distinguishing HCC patients. The α-diversity indices of both the HCC12 and HCC34 groups were significantly lower than those of healthy controls (Figures 1B–E). Specifically, analyses via the Chao1, Shannon, Simpson, and Pielou’s evenness indices indicated that both the diversity and abundance of the gut microbiota in HCC12 and HCC34 patients were reduced compared with the CG group. The Chao1 index further demonstrated that microbial richness in the HCC12 group was lower than that in the LC group. Pielou’s evenness index showed that species evenness in the LC group was decreased relative to the control group. Although the HCC34 group did not differ significantly from the LC group, all α-diversity indices were lower than those in the LC group.

The distance matrix heatmap showed that the HCC34 group exhibited the greatest dissimilarity with the CG group, whereas the HCC12 group showed the smallest dissimilarity with the LC group (Figure 1F). NMDS analysis indicated that samples from the CG group were tightly clustered, while the HCC12 and HCC34 groups were more dispersed, and the LC group overlapped with the other three groups (Figure 1G). Statistical analysis using unweighted and weighted UniFrac Kruskal–Wallis tests revealed that β-diversity differed significantly between the CG group and the other groups, as well as between the LC group and the other groups; however, no significant differences were observed between the HCC12 and HCC34 groups (Figures 1H,I).

Analysis of the top 15 taxa at the order, family, and genus levels showed that the relative abundances of Enterobacterales, Enterobacteriaceae, and Escherichia-Shigella (the old NCBI hierarchical classification) progressively increased from the LC to HCC12 and HCC34 groups, reaching their highest levels in the HCC34 group (Figures 2A–C). Furthermore, the relative abundances of Lactobacillales, Enterococcaceae, and Enterococcus were also elevated in the HCC34 group (Figures 2A–C). This trend indicates a stepwise enrichment of specific gut microbiota as liver disease progresses from cirrhosis to HCC.

To identify taxa with significant differences between groups, we applied multiple statistical approaches. At the genus level, Escherichia-Shigella, Enterococcus, and Succinivibrio, and at the family level, Enterobacteriaceae, Enterococcaceae, and Succinivibrionaceae were dominant in the HCC34 group, with Enterococcus and Enterococcaceae showing particularly pronounced enrichment (Figures 2D,E). LEfSe analysis further identified four significantly enriched taxa across the groups (LDA score >4, p < 0.05). Notably, taxa within the same evolutionary lineage, Enterobacterales, Enterobacteriaceae, and Escherichia-Shigella, were enriched in HCC34, as were Enterococcaceae and Enterococcus (Figures 2F,G). The HCC12 group exhibited significant enrichment of Ruminococcus, whereas the LC group was dominated by Streptococcaceae, Streptococcus, and Veillonella (Figure 2F). Kruskal–Wallis rank-sum tests of family- and genus-level relative abundances yielded consistent results: Enterobacteriaceae, Enterococcaceae, Lactobacillaceae, Escherichia-Shigella, and Enterococcus were enriched in HCC34 (Figures 2H,I; Supplementary Figures S1A–E). Additionally, Ruminococcus was enriched in HCC12 and the control group but markedly reduced in HCC34, whereas Veillonella was elevated in both HCC34 and LC groups (Figure 2I; Supplementary Figures S1F,G).

We performed LEfSe analysis specifically within the HCC groups, which revealed that Enterococcus, Bacilli, Lactobacillales, and Enterococcaceae were significantly associated with distinguishing the HCC34 group (LDA score >4, p < 0.05) (Supplementary Figures S2A,B).

We examined the correlations between the top 15 family- and genus-level taxa identified by Kruskal–Wallis rank-sum tests and clinical features. At the family level (Figure 3A), Enterococcaceae, Enterobacteriaceae, Streptococcaceae, and Lactobacillaceae were positively correlated with ALT, AST, prothrombin time (PT), alkaline phosphatase (ALP), TBil, and Child–Pugh scores, with Enterococcaceae showing the strongest associations. In contrast, Eubacterium coprostanoligenes group, Rikenellaceae, and Oscillospiraceae were negatively correlated with ALT, AST, TBil, and Child–Pugh scores. These findings indicate that liver function impairment is accompanied by a relative decrease in Eubacterium coprostanoligenes group, Rikenellaceae, and Oscillospiraceae, and a relative increase in Enterococcaceae, Enterobacteriaceae, Streptococcaceae, and Lactobacillaceae. Conversely, ChE levels were significantly negatively correlated with Enterococcaceae and Streptococcaceae, and positively correlated with Ruminococcaceae, Moraxellaceae, Erysipelatoclostridiaceae, Eubacterium coprostanoligenes group, Rikenellaceae, Oscillospiraceae, and Pseudomonadaceae. Body mass index (BMI) exhibited a trend of negative association with liver cell injury, consistent with disease progression and malnutrition-related wasting.

At the genus level (Figure 3B), Agathobacter, Alistipes, Ruminococcus, Blautia, Dialister, and Roseburia were negatively correlated with liver injury-related clinical parameters, including TBil, Child–Pugh scores, ALP, ALT, and AST. In contrast, Enterococcus, Veillonella, and Clostridioides showed positive correlations with these clinical features. These findings suggest that Enterococcus, Veillonella, and Clostridioides, together with their corresponding families (Enterococcaceae, Enterobacteriaceae, Streptococcaceae, and Lactobacillaceae), are associated with the severity of liver injury.

In addition, Erysipelatoclostridiaceae, Lachnospiraceae, and Blautia were negatively correlated with the HCC biomarkers Pivka II, AFP, and AFP-L3 (Figure 3C), which may be related to treatment responses in HCC. Conversely, Streptococcaceae and Streptococcus showed positive correlations with these biomarkers (Figure 3D), suggesting the potential diagnostic value of Streptococcus in HCC. Ruminococcaceae exhibited a notably negative association with Pivka II (Figure 3C), a biomarker with high specificity for diagnosing HCC, indicating that Ruminococcaceae may also possess diagnostic potential for hepatocellular carcinoma.

Among the 38 patients with HCC, we first constructed machine learning classification models for advanced HCC using clinical variables, including RF, GBDT, and XGB. The XGB model demonstrated the best performance (AUC = 0.889) (Figures 4A–H). We then developed models using the top 15 family- and genus-level microbial taxa identified by Kruskal–Wallis rank-sum tests. The XGB model based solely on microbial features showed strong discriminatory ability for identifying HCC34, with the highest performance (AUC = 0.926) (Figures 4I–P). Finally, integrating clinical features with key microbial taxa further improved the diagnostic performance for HCC34, achieving optimal discrimination (AUC = 1) (Figures 5A–D). To minimize overfitting, we applied both bootstrap resampling and k-fold CV for internal validation of each model. The results indicated that differential microbial taxa possess diagnostic value for HCC staging, and that combining clinical variables with microbial features yields the best performance (XGB with bootstrap: AUC = 0.943; XGB with k-fold CV: AUC = 0.766) (Figures 5E–H).

Feature importance analysis revealed that Enterococcus, Pivka II, Child–Pugh scores, Blautia, ALT, ALP, TBil, Moraxellaceae, and Pseudomonas were key contributors to distinguishing HCC34 from HCC12 (Figures 5I,J). Except for ALT, the remaining eight features were positively associated with HCC34. Notably, PIVKA-II and Child–Pugh scores are integral components of the CNLC staging system (Zhou et al., 2025). Enterococcus, Moraxellaceae, and Pseudomonas were specifically enriched in the HCC34 group, consistent with the earlier differential abundance analyses. Among them, Enterococcus showed the strongest discriminatory value for differentiating HCC34 from HCC12.

Using the clinical characteristics and key microbial taxa reported by Zhang et al. (2021) for 74 patients with HCC, we reclassified their early and intermediate groups as HCC12 and their terminal group as HCC34 (Supplementary Table S1). External validation using our optimized XGB model demonstrated that the top nine features selected by SHAP (Figure 5I) showed excellent discriminatory performance for advanced HCC (AUC = 1) (Figure 6A). Independent validation with bootstrap resampling and k-fold CV also yielded robust results (bootstrap: AUC = 0.924; k-fold CV: AUC = 1) (Figures 6B,C). Moreover, the feature importance ranking derived from the XGB model, which highlighting Child–Pugh scores, ALP, ALT, Enterococcaceae, and Lachnospiraceae, was consistent with the key features identified in our cohort (Figure 6D). Evaluation metrics for all models are provided in Supplementary Table S2.

Based on t test analyses of microbiome functional predictions, the HCC34 group exhibited enrichment of peptidases, glutathione metabolism, and K03564 (thioredoxin-dependent peroxiredoxin) compared with the CG group (Supplementary Table S3). The enhanced peptidase activity may reflect the hypermetabolic and catabolic state characteristic of advanced HCC, in which increased proteolysis accelerates the breakdown of luminal proteins. The upregulation of glutathione metabolism and thioredoxin-dependent peroxiredoxin suggests pronounced oxidative stress within the gut microenvironment. Glutathione represents a major endogenous antioxidant, while K03564 is crucial for detoxifying peroxides. Their concurrent elevation indicates increased levels of reactive oxygen species (ROS), triggering peroxiredoxin-mediated peroxide removal and compensatory activation of glutathione metabolism. This oxidative stress–driven feedback loop may further exacerbate disease progression in advanced HCC.