Specific gut microbiome signature predicts hepatitis B virus-related hepatocellular carcinoma patients with microvascular invasion

Abstract

Background

We aimed to assess differences in intestinal microflora between patients with operable hepatitis B virus-related hepatocellular carcinoma (HBV-HCC) with microvascular invasion (MVI) and those without MVI. Additionally, we investigated the potential of the microbiome as a non-invasive biomarker for patients with MVI.

Methods

We analyzed the preoperative gut microbiomes (GMs) of two groups, the MVI (n = 46) and non-MVI (n = 56) groups, using 16S ribosomal RNA gene sequencing data. At the operational taxonomic unit level, we employed random forest models to predict MVI risk and validated the results in independent validation cohorts [MVI group (n = 17) and non-MVI group (n = 15)].

Results

β diversity analysis, utilizing weighted UniFrac distances, revealed a significant difference between the MVI and non-MVI groups, as indicated by non-metric multidimensional scaling and principal coordinate analysis. We also observed a significant correlation between the characteristic intestinal microbial communities at the genus level and their main functions. Nine optimal microbial markers were identified, with an area under the curve of 79.76% between 46 MVI and 56 non-MVI samples and 79.80% in the independent verification group.

Conclusion

This pioneering analysis of the GM in patients with operable HBV-HCC with and without MVI opens new avenues for treating HBV-HCC with MVI. We successfully established a diagnostic model and independently verified microbial markers for patients with MVI. As preoperative targeted biomarkers, GM holds potential as a non-invasive tool for patients with HBV-HCC with MVI.

Keywords:

• Hepatitis B virus-related hepatocellular carcinoma
• microvascular invasion
• gut microbiomes
• microbiome functions
• prediction model

Introduction

Hepatocellular carcinoma (HCC) is a major health concern, ranking as the third leading cause of cancer-related deaths worldwide. Its incidence, notably in China [1], has been steadily increasing. HCC primarily arises from chronic liver diseases and is closely associated with hepatitis B or C viral infections and cirrhosis, particularly hepatitis B virus (HBV) infection in China [2]. Various treatment modalities, including hepatectomy, liver transplantation, transcatheter arterial chemoembolization, radiotherapy, and ablation, have been deployed to address HCC at different stages, aiming for curative or palliative care [3]. Nevertheless, advanced cases often display limited responsiveness to treatments such as tyrosine kinase and immune checkpoint inhibitors [4]. Despite the variety of treatment options available for liver cancer, hepatectomy remains the most commonly used and effective option [5], Owing to advances in perioperative management and surgical techniques, the mortality rate and complications after post-hepatectomy liver cancer patients have greatly decreased [6]. However, early postoperative recurrence and metastasis rates hamper treatment outcomes [7]. Recent insights into microvascular invasion (MVI) have provided a deeper understanding of local invasion, distant metastasis, and postoperative recurrence in liver cancer. MVI is an independent risk factor for tumor recurrence following hepatectomy in patients with HCC [8]. Patients exhibiting recurrent HCC with MVI at initial first resection typically experience poorer survival rates [9], necessitating wider resection margins for curative surgery when MVI is present [10]. MVI is conventionally identified through pathological examination of surgical specimens, occurring in 15–60% of resected HCC samples [11]. Preoperatively predicting MVI could facilitate more tailored surgical strategies and reduce postoperative tumor recurrence [10]. Recently, there has been increasing evidence indicating the possibility of accurate and specific preoperative prediction of MVI. Despite ongoing efforts, current clinical parameters, serological markers, and predictive models fall short of a good assessing the preoperative risk of MVI in patients with HCC [12,13]. However, predicting MVI to establish a precise hepatectomy range before surgery remains challenging, and new methods are needed to overcome this.

The gut microbiome (GM), often dubbed the ‘hidden organ’, plays a crucial role in human health and has garnered considerable interest in recent years [14]. The GM strongly influences the metabolic, endocrine, and immune systems, both locally at the mucosal level and systemically at the host level [15]. The dysregulation of GM has been implicated in various diseases, including mental illness [16], type 2 diabetes [17], lung cancer [18], and colorectal cancer [19], while GM signatures have exciting diagnostic potential [18]. The liver’s role in metabolism has prompted research into the interaction between GM and liver diseases. The GM plays a significant role in liver damage initiation and progression [20]. For example, the step-up abundance of Bacteroides, Prevotella, and ethanol-producing GM and their detrimental metabolites in nonalcoholic fatty liver disease (NAFLD) can induce liver inflammation by furnishing toll-like receptor ligands and changing lipid metabolic pathways [21]. Chronic hepatitis B can cause advanced liver fibrosis, cirrhosis, and HCC, concurrently dysregulating the GM and reducing its diversity [22]. Bacteroides, Proteus, and Sclerotinia are the characteristic GM genera at the onset of hepatitis B. With disease progression (from hepatitis B to liver fibrosis/cirrhosis), the abundance of butyrate-producing bacterial families (Lachnospiraceae and Ruminococcaceae) decreases, which facilitates the permeability of the intestinal walls and the gut-liver axis, converting some bacteria from beneficial to pathogenic [23,24]. These facts further confirm a close correlation between the GM and the liver. Oh et al. conducted an external validation study on a characteristic GM using a cohort of three relatives with NAFLD-cirrhosis. They employed a random forest model to predict the diagnosis of liver cirrhosis and obtained favorable results [25]. Furthermore, Zheng et al. discovered that the intestinal microbiota can influence various metabolic pathways, potentially expediting the progression of the disease and serving as superior non-invasive biological diagnostic markers for HCC [26]. Nonetheless, to date, there have been no reports on the gut microbial characteristics in clinical patients with HCC and MVI.

In this study, we aimed to compare the intestinal microflora of patients with operable hepatitis B virus-related hepatocellular carcinoma (HBV-HCC) with and without MVI. Additionally, we aimed to examine the characteristics and functions of the intestinal microflora and evaluate GM’s potential as a non-invasive biomarker for patients with MVI. This study involved a discovery cohort comprising 46 HCC patients with MVI and 56 without. Using 16S rRNA gene sequencing analysis, we delineated the microbial spectrum of HCC patients with MVI, constructed a prediction model, and identified a specific microbiome signature with potential as a non-invasive biomarker for HCC patients with MVI.

Materials and methods

Study design

Patients diagnosed with HBV-HCC, both with and without MVI, underwent hepatectomy at the Department of Hepatobiliary Surgery, Guangxi Medical University Cancer Hospital (Nanning, China). Pre-operative faecal stool samples were collected during two periods: from September 2020 to August 2021 for the development cohort (MVI group = 46, non-MVI group = 56) and from September 2021 to December 2021 for the validation cohort (MVI group = 17, non-MVI group = 15). The diagnosis of MVI was based on postoperative pathological sections. Our study included patients who were diagnosed with primary liver cancer for the first time and who had not received any antitumor treatment before surgery (specifically no immunosuppressive regimen, These patients were only surgically treated with their tumors before this outcome). None of the participants had a history of drug allergies, asthma, alcohol addiction, nonalcoholic fatty liver disease, chronic infectious diseases, gastrointestinal disease, inflammatory bowel diseases, diabetes, or malignant tumors other than those in the liver. Additionally, participants refrained from consuming acid suppressants, pro-gastrointestinal prokinetic agents, probiotics, or antibiotics for at least four weeks before surgery. Prophylactic antibiotics were administered after surgery for a duration of 48 h and based on the same initial timing. All aspects of the study were conducted consistently. Clinical data, including age, sex, BMI, WBC count, platelet count, total bilirubin level, albumin level, prealbumin level, alanine aminotransferase level, aspartate aminotransferase level, HbsAg status, hepatitis B virus DNA level, INR, serum AFP level, γ-glutamyl transpeptidase level, neutrophil-to-lymphocyte ratio, liver cirrhosis Child–Pugh grade, 15-min ICG retention rate (ICG-R15), China liver cancer (CNLC) stage, Barcelona Clinic Liver Cancer (BCLC) stage, CT tumor size before the operation, tumor number, presence of portal hypertension, presence of an envelope or intact envelope, and Ki67 cellular index, were collected. Fecal samples were collected 1 day before surgery and immediately stored at −80 °C.

DNA extraction and 16S rRNA gene amplicon sequencing

DNA from fecal samples was extracted utilizing cetyltrimethylammonium bromide for bacterial DNA detection and library construction. The V3-V4 region of the 16S rRNA gene of the small bacterial subunit was amplified using the primer pair 341 F (5′-CCTAYGGGRBGCASCAG-3′)/806R (5′-GACTACHVGGGTATCTAATCC-3′) and the universal sequence. The amplicon libraries were constructed using TruSeq DNA amplicon libraries (Thermo Fisher Scientific, Waltham, MA, USA) and the TruSeq DNA PCR-Free Sample Preparation kit (Illumina, San Diego, CA, USA). Library quality was assessed using a Qubit 2.0 Fluorometer (Thermo Fisher Scientific) and an Agilent Bioanalyzer 2100 system (Agilent, Santa Clara, CA, USA). Sequencing was performed using Illumina Novaseq 6000.

16S rRNA data analysis

Analysis was conducted on a quality-filtered dataset using FLASH (version 1.2.7) and QIIME (version 1.9.1) with the original FASTQ files [27,28]. Operational taxonomic units (OTUs) were defined at 97% sequence similarity using an algorithm and taxonomic annotation at a confidence threshold of 80% [29]. Normalized abundance data of OTUs at various taxonomic levels (phylum, class, order, family, genus, and species) were used for α- and β-diversity and abundance analyses. Analysis of similarity (AMOVA) with weighted UniFrac and unweighted UniFrac formulae was used to assess differences in bacterial composition. Linear discriminant analysis effect size (LEfSe) identified differentially represented biomarkers between the MVI and non-MVI groups. PICRUSt2 was used to predict metabolic pathways based on the Kyoto Encyclopedia of Genes and Genomes database from Integrated Microbial Genome (IMG) data [30].

OTU biomarker identification and model construction

In this phase, we applied a standardized procedure known as counts per million (CPM) [31] to preprocess the absolute OTU matrix, akin to the approach used in single-cell expression profiling. The process involved the following steps: i) aggregating the expression levels of each OTU in the original OTU count matrix, followed by dividing 10,000 by the sum of the expression of each OTU to derive a scale factor vector for each OTU ii) employing these scaling factors to smooth the original expression spectrum using the log1p function, rendering the data more compliant with a Gaussian distribution; iii) After data preprocessing, the standardized OTU expression profiles obtained were used for subsequent feature screening. We conducted three feature screenings using the Boruta method [32] and identified the feature OTUs considered effective across all three screenings. Based on the feature OTU expression profiles from the training set, we employed 10-fold cross-validation for model construction and training evaluation. We applied the SMOTE [33] algorithm to address sample distribution imbalances during training. The performance of the constructed biomarkers and effectiveness of the models were evaluated using receiver operating characteristic (ROC) curves.

Statistical analysis

The preoperative routine demographic and clinical characteristics were compared between the MVI and non-MVI groups using R software version 4.1.2. Measurement data adhering to a normal distribution were presented as median (P25, P75) and compared using the Wilcoxon rank-sum test. Categorical data were presented as case numbers (percentages) and compared using the Chi-square test or Fisher’s exact probability method. Statistical significance was defined as p < 0.05 (two-sided).

α-diversity was assessed using QIIME (version 1.9), while β-diversity non-metric multidimensional scaling (NMDS) and principal coordinate analysis (PCoA) were performed using R software (version 2.15.3) with the ‘WGCNA’, ‘PCoA’, ‘stats’, ‘ggplot2’, and ‘vegan’ packages. LEfSe software was employed with an LDA score threshold of 3. The Wilcoxon rank-sum test was used to identify differences in microflora. Heat map analysis was performed to reveal correlations between differential GM and function. The area under the ROC curve (AUC) was used to represent the ROC effect. The evaluation indicators were calculated using the MLeval (v0.3) R package.

Results

Baseline clinical characteristics

To investigate the potential association between GM changes and MVI, we conducted 16S rRNA gene sequencing using samples from a Chinese cohort, serving as the ‘discovery cohort’, consisting of 46 patients with MVI and 56 non-MVI participants. For validation, we recruited an independent cohort of 17 patients with MVI and 15 patients without MVI (the ‘validation cohort’). The clinical manifestations and demographic characteristics of the discovery cohort are detailed in Table 1. Among all the parameters between the MVI and non-MVI groups, statistically significant differences were observed in albumin level, platelet-to-lymphocyte ratio, γ-glutamyltransferase level, T-helper/T-suppressor ratio, CT tumor size, and Ki67 cellular index. Except for the aspartate aminotransferase level, none of the other indices in the validation cohort showed significant differences (Table S1).

Table 1. Baseline characteristics of the discovery cohort.

Characteristic	MVI Group(n = 46)	Non-MVI Group(n = 56)	p Value
Age(year)			0.301
<45	31 (67.4%)	31 (55.4%)
≥45	15 (32.6%)	25 (44.6%)
Sex			0.573
Female	4 (8.7%)	8 (14.3%)
Male	42 (91.3%)	48 (85.7%)
BMI	21.95 (20.59–23.82)	22.50 (20.17–24.04)	0.734
PLT			0.477
100–300	41 (50.0%)	46 (50.0%)
<100 or >300	5 (10.9%)	10 (17.9%)
RBP	25.35 (20.80–34.30)	26.30 (18.80–31.02)	0.598
TBIL			0.062
≤17.1	36 (78.3%)	33 (58.9%)
>17.1	10 (21.7%)	23 (41.1%)
ALB			0.004*
<35	13 (28.3%)	3 (5.4%)
≥35	33 (71.7%)	53 (94.6%)
PA	200.10 (160.10–242.93)	213.10 (185.75–246.02)	0.723
ALT			>0.999
≤40	32 (69.6%)	38 (67.9%)
>40	14 (30.4%)	18 (32.1%)
AST			0.175
≤40	20 (43.5%)	33 (58.9%)
>40	26 (56.5%)	23 (41.1%)
HBsAg(ng/ml)	626.60 (356.57–792.18)	606.56 (431.30–744.91)	0.888
PLR	122.71 (86.05–159.97)	100.37 (72.30–131.74)	0.037*
NLR	2.10 (1.44–2.65)	1.85 (1.39–2.41)	0.220
MLR	0.27 (0.21–0.35)	0.24 (0.19–0.32)	0.397
PT			0.331
≤13	26 (56.5%)	38 (67.9%)
>13	20 (43.5%)	18 (32.1%)
INR			0.475
≤1.2	43 (93.5%)	55 (98.2%)
>1.2	3 (6.5%)	1 (1.8%)
AFP			0.271
≤400	22 (47.8%)	34 (60.7%)
>400	24 (52.2%)	22 (39.3%)
γ-GGT	71.00 (41.00–119.50)	42.50 (30.75–63.25)	0.005
TCH	4.56 (4.02–5.06)	4.63 (4.28–5.45)	0.382
Ascites			0.757
No	41 (89.1%)	52 (92.9%)
Yes	5 (10.9%)	4 (7.1%)
Th/Ts	1.60 (1.30–2.28)	2.20 (1.60–2.52)	0.037*
Smoking			>0.999
No	30(65.2%)	36 (64.3%)
Yes	16 (34.8%)	20 (35.7%)
Drinking			0.471
No	36 (78.3%)	48 (85.7%)
Yes	10 (21.7%)	8 (14.3%)
Portal hypertension			0.768
No	42 (91.3%)	49 (87.5%)
Yes	4 (8.7%)	7 (12.5%)
With or without envelope			0.533
No	6 (13.0%)	11 (19.6%)
Yes	40 (87.0%)	45 (80.4%)
Is the envelope intact			0.792
No	21 (45.7%)	23 (41.1%)
Yes	25 (54.3%)	33 (58.9%)
Tumor size(CT)			0.001*
1	10 (21.7%)	31 (55.4%)
2	21 (45.7%)	19 (33.9%)
3	15 (32.6%)	6 (10.7%)
Liver Cirrhosis			0.736
No	28 (60.9%)	37 (66.1%)
Yes	18 (39.1%)	19 (33.9%)
Tumor Number (CT)			0.583
1	37 (80.4%)	43 (76.8%)
2	6 (13.0%)	6 (10.7%)
3	3 (6.5%)	7 (12.5%)
BCLC stage			0.728
0	2 (4.3%)	4 (7.1%)
A	30 (65.2%)	37 (66.1%)
B	8 (17.4%)	11 (19.6%)
C	6 (13.0%)	4 (7.1%)
Ki67	0.40 (0.20–0.49)	0.20 (0.10–0.30)	0.001*

The p-value was based on Wilcoxon Rank sum test, Fisher Exact test and Chi-square test. MVI Group: microvascular invasion Group; non-MVI Group; non-microvascular invasion Group; BMI: body mass index; PLT: platelet; RBP: retinol blinding protein; TBIL: total bilirubin; ALB: albumin; PA: prealbumin; ALT: alanine aminotransferase; AST: aspartate aminotransferase; HBsAg: hepatitis B surface antigen; PLR: platelet-to-lymphocyte ratio; NLR: neutrophil-to-lymphocyte ratio; MLR: monocyte-to-lymphocyte ratio; PT: prothrombin time; INR: international normalized ratio; AFP: alpha fetoprotein; γ-GGT: gamma-glutamyl transferase; TCH: total cholesterol; Th/Ts: helper T cell/suppressor T cell; BCLC stage: Barcelona Clinic Liver Cancer stage; Ki67: Ki-67. *pvalue < 0.05.

OTUs and α- and β-diversity

OTUs represent different microbial taxa, with a 97% identity in rRNA sequences used as a cutoff for defining OTUs. Overall, 5,747 OTUs were identified from the discovery cohort, of which 2,024 OTUs were shared between the MVI and non-MVI groups (Figure 1A). Furthermore, both the Shannon and Simpson indices (Figure 1B) and the Ace and Chao1 indices (Figure S1) were comparable between the two groups, indicating similar levels of α-diversity between the groups. Using weighted UniFrac distances and unweighted UniFrac distances, we employed β-diversity for NMDS and PCoA (Figure 2A). The NMDS (weighted UniFrac distance stress = 0.105 and unweighted UniFrac distance stress = 0.140) revealed significant compositional and abundance differences in the microbiomes between the MVI and non-MVI groups based on the PCoA and AMOVA tests using weighted UniFrac (p = 0.003). However, unweighted UniFrac (p = 0.219) showed no statistical differences. Our results demonstrate that the composition of GM differs between HBV-HCC patients with and without MVI.

Figure 1. Identification of gut microbe using metagenomics analysis. (A) Venn diagrams show the common OTUs among MVI Group and non-MVI Group. (B) MVI Group and non-MVI Group comparison of α-diversity using shannon and simpson. *, 0.01<p ≤ 0.05. OTUs: operational taxonomic units; MVI Group: microvascular invasion Group; non-MVI Group: non-microvascular invasion Group.

Figure 2. β-diversity analysis and relative abundance of the top 10 gut microbiome among groups MVI and non-MVI. (A) β-diversity analysis.(a)NMDS plot based on weighted UniFrac distances (b)NMDS plot based on unweighted UniFrac distances.(c)PCoA plot based on weighted UniFrac distances, (d)PCoA based on unweighted UniFrac distances. NMDS, non-metric multidimensional scaling; PCoA, primary coordinate analysis; stress values <0.2 were considered good after dimensionality reduction. p < 0.05 is statistically significant.(B)Relative abundance of the top 10 gut microbiome among groups MVI and non-MVI at genus level. MVI Group: microvascular invasion Group; non-MVI Group: non-microvascular invasion Group.

Specific gut microbiome signature in MVI and non-MVI patients

We examined the abundance of the top 10 dominant core taxa at the phylum, class, order, family, genus, and species levels in the discovery cohort between the MVI and non-MVI groups (Figure S2 and Figure 2B). In terms of phylum abundance, Firmicutes (62.83% and 48.09%, respectively) was significantly less abundant in the MVI group than in the non-MVI group, while Proteobacteria (26.22%) and Bacteroidota (19.17%) were more abundant in the MVI group. At the class level, the dominant bacteria in the MVI and non-MVI groups were Clostridia (37.23% and 52.69%, respectively), Gammaproteobacteria (25.99% and 19.87%, respectively), and Bacteroidia (19.16% and 12.58%, respectively). At the order level, a downward trend was observed in the abundance of Lachnospirales (20.18%) and Oscillospirales (12.27%). In contrast, the abundance of Pseudomonadales (11.34%), Bacteroidales (19.12%), and Enterobacterales (14.26%) showed an upward trend in the MVI group. At the order and family levels, the dominant bacteria in the MVI and non-MVI groups were identical, with some proportional differences. Moreover, the 10 most abundant genera in the MVI and non-MVI groups were Pseudomonas (11.32% and 6.63%, respectively), Prevotella_9 (4.25% and 1.76%, respectively), Escherichia-Shigella (10.30% and 9.54%, respectively), Bacteroides (12.01% and 8.08%, respectively), Enterococcus (0.98% and 1.50%, respectively), Blautia (7.01% and 11.56%, respectively), Faecalibacterium (5.72% and 10.40%, respectively), Megamonas (1.25% and 0.18%, respectively), Subdoligranulum (2.78% and 2.31%, respectively), and Fusobacterium (1.60% and 0.56%, respectively). At the species level, the abundance of Faecalibacterium prausnitzii was considerably decreased; however, the levels of Bacteroides vulgatus and Fusobacterium mortiferum increased.

The analysis of the LEfSe data revealed a significant difference between the gut microbiota of the MVI and non-MVI groups. In total, one phylum, eight genera, and nine species were significantly enriched in patients with MVI, while one phylum, ten genera, and three species were significantly enriched in non-MVI patients. At the genus level, the relative abundance of Blautia, Faecalibacterium, Ruminococcus, Agathobacter, [Eubacterium]_hallii_group, Fusicatenibacter, Erysipelotrichaceae_UCG-003, Anaerostipes, Dorea, and Peptoniphilus was highly enriched in the non-MVI group compared with that in the MVI group. Conversely, the MVI group showed higher relative abundances of Coprobacter, Veillonella, Megasphaera, Dialister, Fusobacterium, Subdoligranulum, Prevotella_9, and Bacteroides (abundance < 0.0001) than did the non-MVI group (Figure 3).

Figure 3. Bacterial taxa that best characterize the groups were identified by using LDA of effect size (LEfSe) on OTU tables among MVI Group vs. non-MVI Group. (A) The bar plot based on the LDA selection. (B) Cladogram representing the taxonomic hierarchical structure. LDA: linear discriminant analysis; OTUs: operational taxonomic units; MVI Group: microvascular invasion Group; non-MVI Group: non-microvascular invasion Group.

Gut microbiota functional signatures

We performed PICRUST2 analysis to predict and test for differences in the KEGG pathways between the MVI and non-MVI groups. Our analysis identified 142 pathways related to the differential gut microbiota at the genus level (Table 2). We observed that 65 pathways increased, while 77 decreased in the MVI group (Figure 4). Branched-chain amino acid (BCAA) metabolism, essential and non-essential amino acid metabolism generally decreased in the MVI group. Additionally, several pathways such as the guanosine and adenosine deoxyribonucleotide de novo biosynthesis pathway, purine and pyrimidine deoxyribonucleoside degradation pathways, 5-aminoimidazole ribonucleoside degradation pathway, purine ribonucleoside degradation, adenosine nucleotides de novo biosynthesis, glycogen degradation pathway, rhamnose biosynthesis pathway, galactose degradation pathway, peptidoglycan biosynthesis pathway, and lactate relative pathway were the main pathways that decreased in the MVI group. However, there was mainly a negative correlation between the characteristic flora of the MVI group and these pathways and a significant positive correlation between these pathways and the characteristic flora of the non-MVI group.

Figure 4. The function of gut microbiome. (A) Functional alteration caused by gut microbiome change through PICRUST2 prediction in MVI Group versus non-MVI Group.(B)correlation between genera and functional modules in MVI Group versus non-MVI Group. The heatmap panel shows the spearman correlation coefficient between the genera and functional modules. Significance levels are expressed as follows: *, 0.01<p ≤ 0.05; **, 0.001<p ≤ 0.01; ***, p ≤ 0.001. MVI Group: microvascular invasion Group; non-MVI Group: non-microvascular invasion Group.

Table 2. Functional alteration caused by microbial change through PICRUSt2 analysis.

Pathway ID	Altered pathway	p Value
NONOXIPENT-PWY	pentose phosphate pathway (non-oxidative branch)	0.007
PWY-7111	pyruvate fermentation to isobutanol (engineered)	<0.001
PWY0-1319	CDP-diacylglycerol biosynthesis II	<0.001
PWY-5667	CDP-diacylglycerol biosynthesis I	<0.001
PWY-6737	starch degradation V	0.013
PWY-5101	L-isoleucine biosynthesis II	0.003
PWY-7222	guanosine deoxyribonucleotides de novo biosynthesis II	0.016
PWY-7220	adenosine deoxyribonucleotides de novo biosynthesis II	0.016
GLYCOCAT-PWY	glycogen degradation I (bacterial)	0.001
PWY-7229	superpathway of adenosine nucleotides de novo biosynthesis I	0.005
ANAGLYCOLYSIS-PWY	glycolysis III (from glucose)	0.007
VALSYN-PWY	L-valine biosynthesis	0.003
ILEUSYN-PWY	L-isoleucine biosynthesis I (from threonine)	0.003
GLYCOGENSYNTH-PWY	glycogen biosynthesis I (from ADP-D-Glucose)	<0.001
PWY-5104	L-isoleucine biosynthesis IV	0.005
PHOSLIPSYN-PWY	superpathway of phospholipid biosynthesis I (bacteria)	<0.001
PWY-2942	L-lysine biosynthesis III	0.009
PWY-7219	adenosine ribonucleotides de novo biosynthesis	0.010
COMPLETE-ARO-PWY	superpathway of aromatic amino acid biosynthesis	0.001
PWY-6126	superpathway of adenosine nucleotides de novo biosynthesis II	0.008
PWY-621	sucrose degradation III (sucrose invertase)	0.002
PWY-7208	superpathway of pyrimidine nucleobases salvage	0.022
PWY-5097	L-lysine biosynthesis VI	0.007
PWY-5686	UMP biosynthesis	0.010
BRANCHED-CHAIN-AA-SYN-PWY	superpathway of branched amino acid biosynthesis	0.004
ARO-PWY	chorismate biosynthesis I	0.001
PWY-5103	L-isoleucine biosynthesis III	0.003
PWY4FS-8	phosphatidylglycerol biosynthesis II (non-plastidic)	<0.001
PWY4FS-7	phosphatidylglycerol biosynthesis I (plastidic)	<0.001
PWY-3001	superpathway of L-isoleucine biosynthesis I	0.013
PWY-6163	chorismate biosynthesis from 3-dehydroquinate	0.003
PWY-6121	5-aminoimidazole ribonucleotide biosynthesis I	0.030
PWY-6609	adenine and adenosine salvage III	0.017
PWY-5100	pyruvate fermentation to acetate and lactate II	0.003
PWY-6122	5-aminoimidazole ribonucleotide biosynthesis II	0.018
PWY-6277	superpathway of 5-aminoimidazole ribonucleotide biosynthesis	0.018
PWY-6317	galactose degradation I (Leloir pathway)	0.001
PWY-6386	UDP-N-acetylmuramoyl-pentapeptide biosynthesis II (lysine-containing)	0.011
PWY-6387	UDP-N-acetylmuramoyl-pentapeptide biosynthesis I (meso-diaminopimelate containing)	0.014
TRNA-CHARGING-PWY	tRNA charging	0.010
PWY-7221	guanosine ribonucleotides de novo biosynthesis	0.030
PEPTIDOGLYCANSYN-PWY	peptidoglycan biosynthesis I (meso-diaminopimelate containing)	0.014
PWY-6385	peptidoglycan biosynthesis III (mycobacteria)	0.014
P161-PWY	acetylene degradation	0.002
SER-GLYSYN-PWY	superpathway of L-serine and glycine biosynthesis I	0.004
DTDPRHAMSYN-PWY	dTDP-L-rhamnose biosynthesis I	0.008
DAPLYSINESYN-PWY	L-lysine biosynthesis I	<0.001
PWY0-1586	peptidoglycan maturation (meso-diaminopimelate containing)	0.006
COA-PWY	coenzyme A biosynthesis I	0.034
1CMET2-PWY	N10-formyl-tetrahydrofolate biosynthesis	0.028
HISTSYN-PWY	L-histidine biosynthesis	0.003
COBALSYN-PWY	adenosylcobalamin salvage from cobinamide I	0.046
PWY-5695	urate biosynthesis/inosine 5′-phosphate degradation	0.017
PWY-6151	S-adenosyl-L-methionine cycle I	0.001
PWY0-1296	purine ribonucleosides degradation	0.001
OANTIGEN-PWY	O-antigen building blocks biosynthesis (E. coli)	<0.001
ARGSYNBSUB-PWY	L-arginine biosynthesis II (acetyl cycle)	0.001
PWY0-1297	superpathway of purine deoxyribonucleosides degradation	0.001
UDPNAGSYN-PWY	UDP-N-acetyl-D-glucosamine biosynthesis I	0.003
ARGSYN-PWY	L-arginine biosynthesis I (via L-ornithine)	0.001
PWY-7400	L-arginine biosynthesis IV (archaebacteria)	0.001
PWY0-1298	superpathway of pyrimidine deoxyribonucleosides degradation	0.003
FOLSYN-PWY	superpathway of tetrahydrofolate biosynthesis and salvage	0.012
PYRIDNUCSAL-PWY	NAD salvage pathway I	0.021
PWY0-1061	superpathway of L-alanine biosynthesis	0.024
FASYN-ELONG-PWY	fatty acid elongation -- saturated	0.001
DENOVOPURINE2-PWY	superpathway of purine nucleotides de novo biosynthesis II	<0.001
GLUTORN-PWY	L-ornithine biosynthesis	0.016
PWY-841	superpathway of purine nucleotides de novo biosynthesis I	<0.001
PWY0-162	superpathway of pyrimidine ribonucleotides de novo biosynthesis	<0.001
PWY-5659	GDP-mannose biosynthesis	0.003
PWY-6612	superpathway of tetrahydrofolate biosynthesis	0.003
PWY-5505	L-glutamate and L-glutamine biosynthesis	<0.001
P4-PWY	superpathway of L-lysine, L-threonine and L-methionine biosynthesis I	0.019
PWY-6969	TCA cycle V (2-oxoglutarate:ferredoxin oxidoreductase)	0.002
FERMENTATION-PWY	mixed acid fermentation	0.031
PWYG-321	mycolate biosynthesis	<0.001
P42-PWY	incomplete reductive TCA cycle	0.001
PWY-7664	oleate biosynthesis IV (anaerobic)	<0.001
PWY-7539	6-hydroxymethyl-dihydropterin diphosphate biosynthesis III (Chlamydia)	0.005
PWY-6147	6-hydroxymethyl-dihydropterin diphosphate biosynthesis I	0.016
PWY-7196	superpathway of pyrimidine ribonucleosides salvage	<0.001
PWY-6125	superpathway of guanosine nucleotides de novo biosynthesis II	<0.001
TCA	TCA cycle I (prokaryotic)	0.005
GLCMANNANAUT-PWY	superpathway of N-acetylglucosamine, N-acetylmannosamine and N-acetylneuraminate degradation	0.033
PWY-6901	superpathway of glucose and xylose degradation	0.009
PWY0-862	(5Z)-dodec-5-enoate biosynthesis	<0.001
PWY-6282	palmitoleate biosynthesis I (from (5Z)-dodec-5-enoate)	<0.001
PWY-5989	stearate biosynthesis II (bacteria and plants)	<0.001
PWY-7228	superpathway of guanosine nucleotides de novo biosynthesis I	<0.001
POLYISOPRENSYN-PWY	polyisoprenoid biosynthesis (E. coli)	0.002
PRPP-PWY	superpathway of histidine, purine, and pyrimidine biosynthesis	0.018
PWY-5384	sucrose degradation IV (sucrose phosphorylase)	0.001
PWY-7187	pyrimidine deoxyribonucleotides de novo biosynthesis II	0.039
PWY0-166	superpathway of pyrimidine deoxyribonucleotides de novo biosynthesis (E. coli)	0.003
PENTOSE-P-PWY	pentose phosphate pathway	0.009
PWY-6471	peptidoglycan biosynthesis IV (Enterococcus faecium)	0.002
PWY-6700	queuosine biosynthesis	<0.001
THISYN-PWY	superpathway of thiamin diphosphate biosynthesis I	<0.001
PWY-6519	8-amino-7-oxononanoate biosynthesis I	<0.001
BIOTIN-BIOSYNTHESIS-PWY	biotin biosynthesis I	<0.001
PWY-5971	palmitate biosynthesis II (bacteria and plants)	0.012
FASYN-INITIAL-PWY	superpathway of fatty acid biosynthesis initiation (E. coli)	<0.001
COLANSYN-PWY	colanic acid building blocks biosynthesis	0.006
PWY-7211	superpathway of pyrimidine deoxyribonucleotides de novo biosynthesis	0.003
PWY-7200	superpathway of pyrimidine deoxyribonucleoside salvage	<0.001
PWY-6703	preQ0 biosynthesis	<0.001
PWY-7184	pyrimidine deoxyribonucleotides de novo biosynthesis I	<0.001
PWY-5154	L-arginine biosynthesis III (via N-acetyl-L-citrulline)	0.002
PWY-6507	4-deoxy-L-threo-hex-4-enopyranuronate degradation	0.006
PWY-5913	TCA cycle VI (obligate autotrophs)	0.003
PWY-7197	pyrimidine deoxyribonucleotide phosphorylation	<0.001
PWY-7323	superpathway of GDP-mannose-derived O-antigen building blocks biosynthesis	0.003
P105-PWY	TCA cycle IV (2-oxoglutarate decarboxylase)	0.016
PWY-5345	superpathway of L-methionine biosynthesis (by sulfhydrylation)	0.041
GALACTUROCAT-PWY	D-galacturonate degradation I	0.034
P23-PWY	reductive TCA cycle I	0.035
PWY-6895	superpathway of thiamin diphosphate biosynthesis II	<0.001
RHAMCAT-PWY	L-rhamnose degradation I	0.004
PWY-1269	CMP-3-deoxy-D-manno-octulosonate biosynthesis I	<0.001
NAGLIPASYN-PWY	lipid IVA biosynthesis	<0.001
PWY0-845	superpathway of pyridoxal 5′-phosphate biosynthesis and salvage	0.030
METH-ACETATE-PWY	methanogenesis from acetate	0.006
P108-PWY	pyruvate fermentation to propanoate I	0.004
PWY-6467	Kdo transfer to lipid IVA III (Chlamydia)	<0.001
PWY490-3	nitrate reduction VI (assimilatory)	0.038
PYRIDOXSYN-PWY	pyridoxal 5′-phosphate biosynthesis I	0.022
PWY-5840	superpathway of menaquinol-7 biosynthesis	0.016
PWY-5838	superpathway of menaquinol-8 biosynthesis	0.013
PWY-5899	superpathway of menaquinol-13 biosynthesis	0.021
PWY-5897	superpathway of menaquinol-11 biosynthesis	0.021
PWY-5898	superpathway of menaquinol-12 biosynthesis	0.021
PWY-7456	mannan degradation	0.019
PWY-6749	CMP-legionaminate biosynthesis I	0.036
KDO-NAGLIPASYN-PWY	superpathway of (Kdo)2-lipid A biosynthesis	0.017
HISDEG-PWY	superpathway of ubiquinol-8 biosynthesis (prokaryotic)	0.006
PWY-7210	pyrimidine deoxyribonucleotides biosynthesis from CTP	0.008
PWY-5304	superpathway of sulfur oxidation (Acidianus ambivalens)	0.001
PWY0-1241	ADP-L-glycero-β-D-manno-heptose biosynthesis	0.033
POLYAMINSYN3-PWY	superpathway of polyamine biosynthesis II	0.036
TEICHOICACID-PWY	teichoic acid (poly-glycerol) biosynthesis	0.001
PWY-7198	pyrimidine deoxyribonucleotides de novo biosynthesis IV	0.009

The pathways related to pyrimidine ribonucleotides and pyrimidine deoxyribonucleoside de novo synthesis and salvage, guanosine and purine nucleotides de novo biosynthesis, tetrahydrofolate biosynthesis pathway, fatty acid metabolism, oleate biosynthesis, palmitoleate biosynthesis, stearate biosynthesis, lipid biosynthesis, and L-rhamnose degradation pathways were markedly increased in the MVI group but decreased in the non-MVI group. These results showed that the significantly enriched characteristic flora of the MVI group was mostly positively linked to the pathways described above. In contrast, the characteristic flora in the non-MVI group showed a negative correlation with these pathways. The data from this study suggest that MVI development may be accompanied by potentially interesting metabolic reprogramming in GM.

Identification of GM signature for MVI prediction

We developed an approach based on reference biomarkers to construct a clinical index for early diagnosis using either sequencing or real-time quantitative PCR. To accomplish this, we first identified the characteristic reference OTUs associated with the MVI (n = 46) and non-MVI (n = 56) groups (OTU markers shown in Table 3). Using preprocessed data and the Boruta method, we obtained nine characteristic OTUs. Subsequently, we used 10-fold cross-validation for performance assessment and trained the optimal random forest machine (RF) model to predict the MVI status (Table 4). The consistency rate (classification accuracy) and false negative rate (FNR) were both considered during the selection process to keep error rates low for predicting diseased individuals. We emphasized the confidence interval of AUC relative minimization, which more accurately predicted the MVI status in the discovery cohort (sensitivity = 76.09%, specificity = 76.79%, consistency rate = 76.47%, FNR = 23.91%, AUC = 79.76%).

Table 3. Nine OTU-based markers.

OTU	Level name
OTU_56	k__Bacteria;p__Firmicutes;c__Clostridia;o__Lachnospirales;f__Lachnospiraceae;g__Lachnospiraceae_ND3007_group;
OTU_57	k__Bacteria;p__Firmicutes;c__Negativicutes;o__Acidaminococcales;f__Acidaminococcaceae;g__Phascolarctobacterium;s__Phascolarctobacterium_faecium
OTU_58	k__Bacteria;p__Firmicutes;c__Clostridia;o__Lachnospirales;f__Lachnospiraceae;g__Lachnospiraceae_NK4A136_group;
OTU_596	k__Bacteria;p__Firmicutes;c__Negativicutes;o__Veillonellales-Selenomonadales;f__Veillonellaceae;g__Veillonella;s__Veillonella_atypica
OTU_983	k__Bacteria;p__Firmicutes;c__Negativicutes;o__Acidaminococcales;f__Acidaminococcaceae;g__Acidaminococcus;
OTU_1590	k__Bacteria;p__Firmicutes;c__Clostridia;o__Lachnospirales;f__Lachnospiraceae
OTU_1690	k__Bacteria;p__Firmicutes;c__Clostridia;o__Lachnospirales;f__Lachnospiraceae;g__Blautia;s__Blautia_sp_YL58
OTU_2011	k__Bacteria;p__Firmicutes;c__Clostridia;o__Lachnospirales;f__Lachnospiraceae;g__Blautia;
OTU_4668	k__Bacteria;p__Firmicutes;c__Negativicutes;o__Veillonellales-Selenomonadales;f__Veillonellaceae;g__Veillonella

Table 4A. Table showing the prediction performance of random forest machine models using characteristic OTUs (in columns) methods in the discovery cohort.

Level	Model		Sensitivity (%)	Specificity (%)	False negative rate(%)	False positive Rate(%)	Consistency rate	AUC
OTUs	RF	76.09		76.79	23.91	23.21	76.47	79.76

Table 4B. Table showing the prediction performance of random forest machine models using characteristic OTUs (in columns) methods in the validation cohort.

Level	Model		Sensitivity (%)	Specificity (%)	False negative rate(%)	False positive Rate(%)	Consistency rate	AUC
OTUs	RF	82.35		73.33	17.65	26.67	78.12	79.80

The validation cohort was used as an independent test to determine whether the markers of the discovery cohort were generalizable to other samples with high predictive power (sensitivity = 76.09%, specificity = 76.79%, consistency rate = 76.47%, FNR = 17.65%, and AUC = 79.80%) (Figure 5). These results provide evidence that our approach can distinguish between MVI and non-MVI status in individuals with HBV-HCC.

Figure 5. OTU-based diagnostic biomarkers for HBV-HCC with MVI. Receiver operating characteristic (ROC) curves of prediction efficacy for the OTU-based predictors using RF. OTUs: operational taxonomic units; HBV-HCC: hepatitis B virus -related hepatocellular carcinoma; MVI: microvascular invasion; red line and green lines indicated the discovery (MVI Group = 46, non-MVI Group = 56) and validation (MVI Group = 17, non-MVI Group = 15) cohorts, respectively.

Discussion

The presence of MVI, a form of intrahepatic metastasis in HCC, significantly contributes to the poor prognosis of this disease [34]. It is well-established that intrahepatic metastasis is closely related to postoperative recurrence after HCC resection, which is a major obstacle to long-term survival [35]. While a growing body of evidence highlights the link between GM and various liver conditions, the specific differences in GM between the MVI and non-MVI patients with HCC and their implications, particularly in predicting the risk of microtumor thrombus before surgery, have remained largely unexplored.

Our study illuminates distinct GM alterations in HBV-HCC patients with and without MVI, shedding light on potential microbial candidates that may be attributed to the development of MVI. Importantly, we introduce an innovative diagnostic model for HBV-HCC complicated by MVI, leveraging microbial OTU markers, which were independently verified.

Intriguingly, our α- and β-diversity analyses revealed that the gut microbial species in HBV-HCC patients with MVI did not differ significantly from those in patients without MVI. Instead, variations in microbial abundances at different levels were pronounced, indicating an imbalance in the GM between MVI status and non-MVI status. These differences primarily impacted microbial relative abundance rather than overall community diversity. Weighted UniFrac quantitatively detected further evolved variations between samples with different lineages [36]^. Studies have shown that the variation in GM abundance is related to host genetic, etiological, and environmental factors [37]. Additionally, due to continuous evolution, primary tumors exhibit different intra-tumor heterogeneities through different subclonal manifestations. These lineages may have differentiated many generations ago. There is evidence that the evolution of different lineages is related to tumor heterogeneity and that different subclones may lead to different outcomes [38,39]. Tumors are populations that have undergone continuous evolution of cellular heterogeneity and interrelated genetic factors, and GM has also experienced little-known evolutionary trees, pedigree changes, genetic factors, and other significant influences [40]. Growing evidence also suggests that differences in preoperative flora abundance between the MVI and non-MVI states indicate that patients with different tumor subclones may have varying effects on flora, potentially influencing the prognosis of patients undergoing surgery.

Our study identified specific alterations among the top 10 abundant GM taxa, primarily at the genus level in pre-surgery MVI and non-MVI patients. This level of analysis provides a more precise reflection of species differences within the 16S rRNA gene [41]. Employing LEFSe, our findings pinpointed 17 characteristic differences in genus abundance between MVI and non-MVI patients (including Blautia, Faecalibacterium, Ruminococcus, Agathobacter, [Eubacterium]_hallii_group, Fusicatenibacter, Erysipelotrichaceae_UCG-003, Anaerostipes, Dorea, and Peptoniphilus in non-MVI, and Coprobacter, Veillonella, Megasphaera, Dialister, Fusobacterium, Subdoligranulum, Prevotella_9, and Bacteroides in MVI).

Zhang et al. [42] identified the enrichment of Prevotella_9, Bacteroides, Megasphaera, and Veillonella in lung tumor tissues, suggesting their potential as invasion indicators. The abundance of Dialister in the gastric mucosa of patients with gastric cancer is higher than that of precancerous individuals, which is related to the progression of the disease [43]. Fusobacterium and Subdoligranulum promote inflammation [44,45]. Additionally, Coprobacter is associated with intestinal barrier and inflammation promotion [46], while Blautia induces regulatory T cell proliferation, which may help limit inflammation [47]. Faecalibacterium and Ruminococcus can induce the synthesis of anti-inflammatory regulatory T-cells [48]. Agathobacter is more associated with lymph node metastasis in colorectal cancer [49]. Microbes related to butyrate production ([Eubacterium]_hallii_group, Dorea, and Anaerostipes) are enriched in patients with rectal cancer who benefited from neoadjuvant chemotherapy and play a protective role in the host body [50]. The prognosis of older patients with liver cancer is significantly better than that of young patients, while the significant enrichment of Fusicatenibacter and Erysipelotrichaceae_UCG-003 in the intestinal bacteria of older patients is negatively correlated with strong tumor proliferation and metastasis [51]. The abundance of Peptoniphilus is higher in the advanced stages of endometrial carcinoma than in the early stages [52]. Thus, we argue that certain characteristic bacterial genera may play key roles in the MVI and non-MVI status of HBV-HCC.

Our PICRUST2 analysis showed significant differences in the microbiome functional abundance between the MVI and non-MVI groups, particularly in core functions. The relationship between these differential functions and signature bacteria in both states underscores their potential significance. Interestingly, many metabolic function modules exhibited opposing trends in the MVI and non-MVI patients. This reciprocal regulation among metabolic reprogramming-induced cancer cell metabolism, the host immune system, and GM is vital for regulating cancer progression, including malignant progression and the proliferation change of cancer cells [53,54]. A previous study reported that BCAAs induce cell cycle arrest, leading to apoptosis in HCC cells [55]. Furthermore, the inhibition of tumor cell proliferation can also be accomplished by inducing apoptosis in tumor cells [56].

Research has shown that asparagine and aspartic acid can be interconverted, similar to glutamic acid and glutamine. Essential amino acids, including threonine, lysine, and histidine, are synthesized from aspartate. Additionally, aspartic acid-related enzymes regulate amino acid metabolism by coordinating the uptake of non-essential amino acids such as arginine and serine [57]. Therefore, aspartic acid may act as a bridge in amino acid metabolism. In contrast, asparagine inhibits apoptosis and has potential anti-proliferative activity in tumor cells [58]. However, our study found that the metabolism of these amino acids decreased in the MVI group and became active in the non-MVI group. The characteristic bacterial genera in the MVI group showed an almost negative correlation with the metabolism of these amino acids, whereas the non-MVI group showed the opposite trend. These findings suggest that multiple amino acids are needed for cancer cell growth and distant metastasis, and depleting certain amino acid sources may be a possible therapeutic method for preventing recurrence or metastasis [57]. Deoxyribo­nucleotides and ribonucleotides can interconvert via a central pathway, and an increase in the de novo anabolic metabolism of pyrimidine DNA can result in significant functions of genes that promote cancer cell proliferation [59]. Additionally, tumor cells rely on the pyrimidine deoxyribonucleic acid or pyrimidine ribonucleoside salvage pathway for proliferation; thus, inhibiting this pathway can reduce tumor growth [60]. Tumor activation can fulfill the needs of highly proliferative cells by increasing the contribution of nucleotides to the purine biosynthesis pathway [61]. Research has shown that accelerated tumor growth exhibits an increased dependency on de novo guanosine nucleotide biosynthesis [62]. A study revealed that serine is the main precursor of purine nucleotide biosynthesis, and the depletion of purine deoxyribonucleotides may be the effect of serine reduction, which is negatively related to tumor proliferation to some extent. Moreover, purine deoxynucleosides have been found to inhibit tumor proliferation [57,63]. Similarly, most studies suggest that tumor cell proliferation is inhibited by 5-aminoimidazole ribonucleotides, which can inhibit the growth of various cancer cells (such as prostate cancer) in vitro through AMPK activation [64,65]. The role of adenosine in regulating human inflammation is complex and involves both pro-inflammatory and anti-inflammatory functions [66]. One of the hallmarks of cancer is inflammation, which results in the exhaustion of immune cells, leading to tumor immune escape and progression [67]. Under certain conditions, adenosine can protect cells against injury and inhibit inflammation [68]. Based on these findings, purine and pyrimidine nucleotides can be considered critical regulators of both the direct and indirect pro-metastatic potential of tumor cells, which is consistent with the results of our study. High proliferation of tumor cells can be maintained through fat metabolism and accompanied by a hypoxic environment, which can also promote tumor progression. In hepatoma cells, disruption of the Akt/mTOR signaling pathway reduces fatty acid synthesis and upregulates peroxisome proliferator-activated receptor α, thereby increasing hepatoma cell proliferation and metastatic tumor formation in cell lines and mouse models [69]. Palmitic acid causes epigenetic changes in cancer cells, leading to tumor recurrence. Palmitic acid has been shown to promote tumor metastasis in cancer studies, and it does not need to be ingested over a long period to have this adverse effect [70]. Palmitic acid can be further converted to palmitoleic acid [71], which can elongate into oleate acid and stearic acid [72]. Our results suggest that several metabolic enhancements related to fat metabolism can promote tumor migration.

Research has indicated that the inhibition of dihydrofolate reductase results in the depletion of intracellular tetrahydrofolate, limiting thymidylate synthesis and resulting in apoptosis or autophagy of highly proliferative cells [73]. The metabolic reprogramming of tumor initiation involves cancer cells absorbing glucose and storing it intracellularly as energy in the form of glycogen, which accumulates and inhibits the Hippo signaling pathway through liquid-liquid phase separation, activating downstream proto-oncoprotein YAP and driving tumor initiation and invasion [74]. Zhuang et al. found that rhamnose lectins exist on the surface of some tumor cells, which can bind rhamnose specifically, and glycation modification of antitumor drugs can enhance the targeting of drugs and reduce the damage to normal cells [75]. Galectins also play a key role in tumor progression and metastasis by interfering with certain mechanisms [76]. An increase in lysine-synthesizing proteins leads to an overall increase in peptidoglycan synthesis, which can enhance metastasis [77]. The reduction in tumor growth was accompanied by a significant increase in intratumor lactate, consistent with our observations [78]. At the genus level, our study found that the tumor growth mechanism shown by the core functions of these characteristic GM may potentially impact patients with MVI (Figure S3). Therefore, it is possible to prevent ineffective interventions and improve prognosis.

In addition to characterizing the GM of patients in both groups, our study proposes preoperatively predicting MVI status based on species OTU levels. We found that OTU-level-based predictions were a good effect and successfully achieved independent validation with good reproducibility, possibly due to the use of detailed 16S rRNA gene sequences. Among the nine most promising OTU-based microbial markers, some coincide with our signature bacteria (Blautia and Veillonella), which differentiated the MVI group from the non-MVI group at the genus level and are linked to specific host functions. We suggest that the signature of MVI-associated microbial dysbiosis may be generalized in a manner similar to the characterization of GM in patients with early-stage HCC [26]. These results suggest that GM-targeted biomarkers could be potential non-invasive tools for preoperative MVI diagnosis.

Our study had some limitations. While 16S rRNA gene sequencing provides valuable insights into the GM, metagenomic sequencing offers a more comprehensive genetic view. Additionally, we observed an association between the MVI and non-MVI groups and the intestinal microbiome but did not establish direct causality. Finally, our findings achieved good universality in the independent verification group; nevertheless, further confirmation in a larger cohort is needed. Future research employing metagenomic sequencing, strain isolation and culture, and animal models will enhance our understanding of the mechanisms of intestinal bacteria in the state of MVI while expanding the cohort size promises to improve predictive accuracy.

Conclusions

We have delineated the systemic alterations in GM between MVI and non-MVI states in HBV-HCC patients. We identified microbial signatures associated with each status and developed a highly accurate preoperative diagnostic tool based on OTU levels for determining MVI status. Overall, our study reveals critical insights into the role of GM in HBV-HCC, particularly in relation to MVI. These findings not only enhance our understanding of the disease but also pave the way for developing innovative diagnostic and therapeutic strategies. Harnessing the power of microbial signatures for early MVI detection holds promise for improving patient outcomes and warrants further investigation in larger clinical cohorts.

Ethical approval and consent to participate

This study was conducted in accordance with the 2013 revision of the Declaration of Helsinki. The Guangxi Medical University Cancer Hospital Research Ethics Committee approved this study (No. KY2019009), and all participants provided informed consent.

Author contributions

(I) Conception and design: Yu-Chong Peng, Lu-Nan Qi; (II) Administrative support: Lu-Nan Qi; (III) Provision of study materials or patients: Yu-Chong Peng, Jing-Xuan Xu, Xue-Mei You, Lu-Nan Qi; (IV) Collection and assembly of data: Yu-Chong Peng, Jing-Xuan Xu, Xue-Mei You, Le-Qun Li, Lu-Nan Qi; (V) Data analysis and interpretation: Yu-Chong Peng, Jing-Xuan Xu, Xue-Mei You, Yi-Yue Huang, Lu-Nan Qi; (VI) Manuscript writing: All authors; (VII) Final approval of manuscript: All authors.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [81972306 & 82273405].