A distinct microbiota signature precedes the clinical diagnosis of hepatocellular carcinoma

ABSTRACT

Oral, gut, and tumor microbiota have been implicated as important regulators in the carcinogenesis and progression of gastrointestinal malignancies. However, few studies focused on the existence and association of resident microbes within different body regions. Herein, we aim to reveal the durability of the oral-gut-tumor microbiome and its diagnostic performance in hepatocellular carcinoma (HCC). Our study included two cohorts: a retrospective discovery cohort of 364 HBV-HCC patients and 160 controls with oral or fecal samples, a prospective validation cohort of 91 cases, and 124 controls for matching samples, as well as 48 HBV, and 39 HBV-cirrhosis patients for gut microbial patterns examined by 16S rRNA gene sequencing. With the random forest analysis, 10 oral and 9 gut genera that could distinguish HCC from controls in the retrospective cohort were validated among the prospective matching participants, with area under the curve (AUC) values of 0.7971 and 0.8084, respectively. When influential taxa were merged, the AUC of the consistent classifier increased to 0.9405. The performance continued to improve to 0.9811 when combined with serum levels of alpha-fetoprotein (AFP). Specifically, microbial biomarkers represented by Streptococcus displayed a constantly increasing trend during the disease transition. Furthermore, the presence of several dominant microbiota species was confirmed in hepatic tumor and non-tumor tissues with fluorescence in situ hybridization (FISH) and 5 R 16S rRNA gene sequencing. Overall, our findings based on the oral-gut-tumor microbiota provide a reliable approach for the early detection of HCC.

KEYWORDS:

• Microbiota
• hepatocellular carcinoma
• 16S rRNA gene
• AFP
• fluorescence in situ hybridization
• early detection

Introduction

Hepatocellular carcinoma (HCC) remains one of the most fatal malignancies with a global increase in incidence and mortality, resulting in over 700,000 deaths annually¹. Due to the lack of precisely predictive biomarkers and the asymptomatic nature of the disease, a considerable proportion of patients are detected at an advanced stage to miss the optimal timing for surgical eradication. To date, HCC has been diagnosed with serum AFP levels and imaging examination². However, as the sole standard predictor for HCC, AFP levels can be elevated in other conditions, such as active hepatitis, gonad embryoma, secondary liver cancer, and pregnancy, restricting the disease specificity to some degree. Therefore, it is urgent to develop a brand-new indication with reasonably high accuracy to boost life expectancy.

Recently, there has been a surge of interest in the etiological role of gut microbiota in a variety of cancers. Referring to HCC, it has a complex etiology that comprises alcohol, hepatitis virus, and aflatoxin³, which are frequently accompanied by visible alterations in the microbiome. For instance, convergent evidence has confirmed that hepatitis virus infection can elicit structural changes in the intestinal barrier and increase microbial translocation via the gut-liver axis, resulting in impaired mucosal immunological function, and an excessive long-term immune response that contributes to the development of liver diseases, including fibrosis, cirrhosis, and even HCC^4,5. Simultaneously, the oral microbiota has been widely linked to digestive system diseases, such as oral squamous cell carcinoma (OSCC)⁶, inflammatory bowel disease (IBD),⁷ and liver cirrhosis⁸, indicating microbial migration of the intestine from the oral cavity. In particular, significant changes occur in the oral microbiome of patients with liver carcinoma and cirrhosis, which are expected to distinguish HCC from healthy conditions as noninvasive biomarkers⁹. Moreover, altered gut microbiome profiles have been implicated in the diagnosis of HCC, and Ren et al. established a diagnostic model with successful cross-region validation¹⁰. However, previous studies merely focused on the link between the oral or fecal microbiome and HCC, which remains as yet unexplored in the translation of combined improvements into a novel microbiome signature for early detection as a signal of imminent HCC development.

Nowadays, distinct intratumor microbiota compositions have been verified in several cancers with low biomasses that are mostly localized to cancer and immune cells, and are probably generated from circulating bacteria¹¹. In addition, it was identified that the dominant microbiome Proteobacteria in pancreatic cancer was similar to that in duodenal microbiome makeup, which reflected retrograde bacterial invasion from the gut to the pancreatic duct¹². Similarly, the presence of tumor-associated microbiota has been characterized in HCC patients with hepatitis virus¹³, which may assist in illustrating the potential pathological contribution of intratumoral microbiota from either the portal vein or bile duct to HCC. For this reason, the so-called oral-gut-liver axis has been proposed, in which dysregulated oral bacteria infiltrate the gut, disseminate negative events to the resident intestinal microbiome, and result in liver diseases⁸. Here, we prospectively evaluated a cohort of 91 newly diagnosed patients with HCC and 124 controls. Both oral and fecal samples, along with suitable hepatic tumor and normal tissues, were collected for microbial composition examination with 16S rRNA gene amplicon sequencing, IHC and FISH assays. Based on a set of 19 integrated genera, we identified a distinct microbial signature with high prediction accuracy. More importantly, HCC-enriched species were investigated in tumorous tissues, with a possible relationship with oral and gut residents, disclosing a potential involvement in HCC pathogenesis.

Methods

Patient enrollment and sample collection

Briefly, our study comprised two independent cohorts: a retrospective cohort of 364 patients with newly diagnosed HBV-HCC and 160 control individuals (healthy volunteers and in patients with hepatic hemangioma) without HCC-related risk factors were recruited from the First Affiliated Hospital of Wenzhou Medical University from May 2019 to September 2020. The 16S rRNA gene amplicons of oral or gut microbiome and clinical features were obtained. In the prospective population (ClinicalTrials.gov identifiers NCT04637048), we collected oral (tongue surface) and matching fecal samples from 124 patients with HCC and 91 controls between March and October 2021. Meanwhile, 16S rRNA gene amplicons of the fecal microbiome for the patients with HBV-Cirrhosis (n = 39) and HBV-Noncirrhosis (n = 48) were obtained. Moreover, 46 HCC tissues, 42 paracancerous tissues, and 11 healthy tissues from the control group (hepatic hemangioma) were analyzed by 5 R 16S rRNA gene sequencing, following a strictly sterile procedure with a blank control to avoid cross-contamination¹¹. Clinical data pertaining to age, gender, BMI, and follow-up information were acquired from the same hospital structured questionnaires. Serum AFP levels were examined with an electrochemiluminescence immunoassay according to the manufacturer’s instructions from the Medical Laboratory Center. Written informed consent and ethics approval were obtained from the study participants and the Ethics Committee in Clinical Research (ECCR) of the First Affiliated Hospital of Wenzhou Medical University. Prior to the tongue surface samples collection with the sterile cotton swab, the oral health condition was evaluated among the participants, excluding potential tongue-related diseases that might alter the innate oral microbiota, such as Sjogren’s syndrome, Meige’s syndrome, atrophic glossitis, oral leukoplakia, and tongue fungal infection, etc. In addition, participants were instructed to gargle before sampling during non-eating times. Following that, oral and fecal samples were collected in sterile tubes and temporarily moved to the −20°C freezer in each ward, and were permanently stored at −80°C for DNA extraction and taxonomic analysis. For a subset of individuals, tumor and paracancerous tissues were obtained during surgery, immediately flash-frozen in liquid nitrogen, and stored at −80°C. All samples were shipped with dry ice.

Inclusion criteria

Patients with HBV-related HCC (clinically and histologically diagnosed) and healthy participants aged over 18 were enrolled.

Exclusion criteria

(1) Alcoholic liver cancer, nonalcoholic fatty liver cancer, HCV-related liver cancer, etc. (2) intrahepatic cholangiocarcinoma and mixed liver cancer (3) history of other malignancies (4) antibiotics or probiotics usage within 1 month of sample collection (5) presence of hypertension, diabetes, and other metabolic diseases (6) Crohn’s disease, ulcerative colitis, and inflammatory bowel disease.

The inclusion and exclusion criteria were applied identically to each cohort.

Collection of sterile tissues

Tumor and normal hepatic tissues were dissected in the operating room using sterile scalpels that were changed before dissecting different samples to avoid cross-contamination. The cutting tissue was placed in a sterile and nucleus-free cryopreservation tube with a corresponding blank collection tube as a negative control for 5 R 16S rRNA gene sequencing, which was left in the collection area with the lid open for 30 s and immediately stored with the sample tubes.

DNA/RNA extraction, amplification, sequencing, and analysis

Genomic DNA of oral and fecal samples was extracted with the EZNA® Soil DNA kit in accordance with the manufacturer’s guidelines (D4015, Omega, Inc., USA). Genomic DNA from hepatic tumor and non-tumor samples was extracted with an EZNA® Tissue DNA kit (D3396, Omega, Inc., USA) in a standard protocol. The isolated DNA was immediately stored at −80°C for amplification. For oral and gut DNA, the 16S rRNA gene V3-V4 region was amplified with slightly modified microbiome primers (515F-805 R) in a two-step PCR protocol¹⁴. To characterize the intratumor microbiota, a modified 5 R 16S rRNA gene amplification was done after amplifying five regions on the multiplexed 16S rRNA genes¹¹. The amplification products were purified with AMPure XT beads (Beckman Coulter, USA) and quantified with Qubit (Invitrogen, USA). Subsequently, the amplicon pools were applied for sequencing, and the quantity and size of the amplified libraries were evaluated on Agilent 2100 Bioanalyzer (Agilent, USA). Samples were sequenced on the Illumina NovaSeq platform from the manufacturer’s recommendations, which were supported by Lc-Bio Technologies Co., Ltd (Hangzhou, China). Afterward, raw reads were analyzed with QIIME2 software, and quality filtration was performed in fqtrim software (V.0.9.4) to acquire high-quality clean tags. DABA2 software was employed to construct the sequence and feature tables containing the amplicon sequence variants (ASVs). A taxonomic analysis for species annotation was conducted using BLAST with the database of SILVA and NT-16S. MAFFT software was used to explore the predominant species in different groups and diverse sequence alignments. The alpha diversity of specimens was characterized by the Chao1 and Shannon indices, which were determined with QIIME2. Beta diversity was described by principal component analysis (PCA) and principal coordinates analysis (PCoA) with R ade4 and vegan packages, and the P value was analyzed by the Wilcoxon test.

Shotgun metagenomic sequencing

A number of backup gut samples from the HCC (n = 42) and control (n = 34) groups in the prospective cohort were further analyzed by the shotgun metagenomic sequencing. Total genomic DNA was obtained from backup gut samples with the EZNA® Stool DNA kit. After evaluation of the purity and concentration, genomic DNA was pooled for sequencing on an Illumina Hiseq 2500 platform with NovaSeq Reagent Kits. Raw data were sorted with assigned barcodes and filtered out before analysis. The Partek Flow and MetAMOS pipeline software¹⁵ were applied for the annotation and assembly of sequences. Afterward, the non-redundant Unigenes sets were obtained through sequencing clustering for relevant functional profiles. The Diamond program was used to annotate the functional profiling of KEGG pathway information by comparing it with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database¹⁶. A one-way ANOVA analysis was conducted to assess the statistical difference in functional items with Statistical Analysis of Metagenomic Profiles (STAMP)¹⁷.

Microbial data statistical analyses

All steps of data analyses were performed in the R Statistical Computing platform. The continuous variables were represented as mean ± standard error of mean (SEM) or median with interquartile range, and categorical data were described with frequency. Differences in alpha diversity were verified with an analysis of Wilcoxon. The beta diversity between differences in sample community composition was presented as PCA and PCoA with a weighted Unifrac index. The relationship between community composition and microbiome extrinsic element was analyzed using the PERMANOVA and VEGAN R packages. The FDR method of Bonferroni correction was applied to avoid a type I error (false positives) for multiple testing of multivariate microbiota data. The Bonferroni-adjusted P values were considered to compensate for the significant level [α = 0.05, (n = 812 for retrospective oral test, P = 6.2 × 10⁻⁵; n = 989 for the retrospective fecal test, P = 5.1 × 10⁻⁵; n = 766 for prospective oral test, P = 6.5 × 10⁻⁵; n = 782 for the prospective fecal test, P = 6.4 × 10⁻⁵)]. Wilcoxon-rank-sum analysis was utilized to evaluate the significant difference in microbial abundance with the STAMP software, and the statistically enriched microorganisms were determined using a linear discriminant analysis effect size (LEfSe) algorithm with an LDA score of more than three (Bonferroni-adjusted P < 0.05). The predominant data flow that occurred during the disease progression at phylum and genus levels was visualized with a bubble plot, heatmap and Sankey plot using the R pheatmap, stats, ggplot2, and ggalluvia packages, respectively. The SparCC algorithm¹⁸ was performed to identify the microbial correlation of biomarkers between HCC and control groups (correlation p < 0.05 and |rho| > 0.1). Redundancy Analysis (RDA) was conducted to select several relevant clinical variables from the fecal microbiota with the R vegan package.

Establishment of multivariable statistical models and evaluation

Classification models were established with the random forest algorithm in the RandomForest R package. Following that, taxa with poor richness and prevalence were excluded from the multivariable statistical training models for HCC prediction. Certain taxa identified between HCC and the control group were integrated separately for the model construction of the oral and fecal microbiome. The area under the receiver operating curve (ROC) with validation scores was used to evaluate the performance in the prediction models. Alternatively, features with a ROC<0.5 were eliminated after the cross-validation, therefore sorting out features enriched in HCC. The selected influential taxa in oral and fecal microbiota were merged to culminate in a new classifier. Likewise, the ROC value of the combined classifier was calculated by the pROC R package as well. The oral and fecal samples were collected for taxonomic analysis immediately after the participants were enrolled in the validation cohort, and genera-based signatures were applied for the identification of HCC patients based on the imaging examination plus AFP or pathological test. All data processing, model training, predictions, and evaluation were conducted by the SIAMCAT R package.

Evaluation of intestinal permeability

Intestinal permeability was evaluated by quantification of serum lipopolysaccharide (LPS), zonulin-1 (ZO-1) and fecal calprotectin. Accordingly, a competitive enzyme-link immunosorbent assay (ELISA) was conducted for the quantitative determinant of ZO-1 (ZO-1, Immundiagnostik AG, Bensheim, Germany), and calprotectin (Elabscience, Houston, USA). The corresponding antibody was pre-coated in the microplate with 100 μl standards, serum, and fecal suspensions (feces were dissolved in PBS and vortexed for 30 min) were incubated at 37°C for 90 min. The biotinylated detection antibodies for ZO-1, calprotectin, and avidin-HRP conjugate were converted to each microplate. The concentration of ZO-1 and calprotectin was determined in a microplate reader with an absorbance of 450 nm, and bicinchoninic acid (BCA, Thermo Fisher Scientific, USA) was used for the quantification of fecal protein concentrations. The bacterial LPS was detected and measured by the chromogenic LAL endotoxin assay (Piscataway, NJ, USA), which was based on a colorimetric reaction in which the endotoxin could activate a proenzyme in the Limulus Amebocyte Lysate (LAL). 100 μl standard, serum sample, and LAL reagent were dispensed into endotoxin-free vials, then the mixture was incubated with 100 μl reconstituted chromogenic substrate solution for 6 min at 37°C. The reaction was stopped with the 500 μl reconstituted stop solution (hydrochloric acid with LAL reagent water), and the absorbance at 545 nm was measured in a microplate reader. The endotoxin concentration was calculated from a standard curve.

Immunohistochemistry (IHC) assays

IHC staining was conducted in a standard procedure, which included deparaffinization and rehydration stages. The antigens were retrieved by slightly boiling the citrate buffer in a microwave for 15 min at low-to-medium power. Endogenous catalase activity was quenched with 3% hydrogen peroxide for 10 min. After blockage with 5% bovine serum albumin for 1 hour, the samples were incubated with primary LPS (Lipopolysaccharide Core, HM6011) and LTA (Lipoteichoic acid, HM2048-200UG) antibodies (1:50) at 4°C overnight, followed by incubation with secondary antibody for an hour at room temperature in a humid chamber. The slices were stained with DAB (1:1000 dissolved in 30% hydrogen peroxide) and terminated with washing water. The sections were then incubated with hematoxylin for 3 min and sealed with neutral resin. The images were captured and analyzed using a Leica DM6000 microscope with CytoVision software.

Fluorescence in situ hybridization (FISH)

The FISH test was carried out using a probe that targeted the 16S rRNA gene sequence of a specific bacterial taxon. The probes were displayed in Supplemental Table S2. HCC tumor specimens were obtained during the surgical resection and immediately stored in liquid nitrogen under sterile circumstances. After that, the samples were preserved in a −80°C freezer for long-term storage. The samples were transferred to an appropriate cutting temperature mold and cut on a cryotome into several sections of 4 μm. All materials were sanitized using ethanol after sample processing. Sample slices of 4 μm thickness were mounted on the slides. Tissues were fixed in fresh 4% paraformaldehyde, and bacterial wall permeabilization was enhanced by 1 mg/ml lysozyme at room temperature (RT) for 10 min, followed by treatment with 20 μg/ml proteinase K for 10 min at RT. The probes were diluted to 1 μM in hybridization buffer with 10% dextran sulfate, 25% formamide, 1 mg/mL E. coli tRNA, 0.02% BSA, and 2× SSC buffer. Tissue sections were hybridized with specific probes overnight at 42°C after pre-hybridization in the hybridization buffer at 37°C for an hour. The unbounded probes were washed with wash buffer for 30 min at RT. The sections were stained with 1 ng/ml DAPI for 5 min at RT in the dark. The germ-free (GF) section denoted a clinical HCC tumor sample that was determined to be sterile and devoid of colonized bacteria using 5 R 16S rRNA gene sequencing, and it was selected as a control for the hybridization protocol. The FISH images were captured with a Nikon Eclipse Ci-L microscope. The images were analyzed and scored on the basis of the fluorescence signal.

Results

Clinicopathological features of participants in the study

All the participants were recruited from our hospital with the same criteria. The characteristics of the two cohorts were depicted in Figure 1a. Clinicopathological variables including age, gender, BMI, alcohol consumption, smoking history, Child-Pugh classification, serum AFP level, and liver function indexes et al. were presented in Figure 1b (for the prospective cohort) and Supplemental Table S1 (for the retrospective cohort). Several basic variables, such as age, gender, and the Child-Pugh score showed no significant difference among the groups, whether in the retrospective or prospective study. In addition, microbial-related factors like obesity, cigarette smoking, alcohol consumption, and non-antibiotic drugs usage, such as chemotherapy agents, proton-pump inhibitors (PPIs), nonsteroidal anti-inflammatory drugs (NSAIDs), atypical antipsychotics (AAPs), osmotic laxative and Chinese medicine, were generally matched. Patients with diabetes or hypertensive who had previously taken metformin, angiotensin II receptor blocker (ARB), angiotensin-converting enzyme inhibitors (ACEI), thiazide, or similar medication were also excluded, removing potential confounders to affect microbiome profiles and group discrimination.

Figure 1. Community analysis of clinicopathological features in the longitudinal study. (a) Overview of the study population. Grey bands between bar plots represent samples of matching body regions within individuals. (b) Clinicopathological information and normalized abundance (log10 transformation and centralization) of influential genera (10 oral genera and 9 gut genera) shown as heatmap between HCC (n = 91) and the control groups (n = 124) in the prospective study with the statistical p values. The blank cell in the “oral” and “fecal” panels represented the unmatching data of corresponding individuals.

Abbreviations: BMI, body mass index; CCI, Charlson comorbidity index; MELD, model for end-stage liver disease; HBV, hepatitis B virus; INR, international normalized ratio; PT, prothrombin time; HDL, high-density lipoprotein; LDL, low-density lipoprotein; ALT, alanine transaminase; AST, aspartate transaminase; PLT, platelet; Tchol, total cholesterol; TG, Triglyceride; Tbil, Total bilirubin; Cre, creatinine; Alb, albumin; GGT, γ-Glutamyl transpeptidase; AFP, alpha-fetoprotein.

The microbial alterations discriminate HCC patients from controls with oral and fecal samples analyzed by the 16S rRNA gene sequencing

In the prospective study, we collected both oral and gut samples from a cohort between March and October 2021 for the 16S rRNA gene sequencing. As illustrated in Figure 2a–d, the HCC group showed separated clusters from the control group with statistical β-diversity (PCA and PCoA analyses, p < 0.05), whether in oral or gut samples. Besides, the statistical taxonomic analysis was conducted at the levels of phylum and genus to determine the prevalent microbiota composition. At the phylum level, Firmicutes, Proteobacteria, Bacteroidetes, Actinobacteria, and Fusobacteria dominated the two groups in terms of relative abundance taxa. In oral samples, Firmicutes and Actinobacteria showed a significant increase in the HCC group compared with the control group (Firmicutes: 26.81% versus 32.91%; Actinobacteria: 9.74% versus 13.38%), whereas the relative abundance of Proteobacteria, and Fusobacteria were significantly enriched in the subjects of the control group (Proteobacteria: 28.45% versus 21.35%; Fusobacteria: 11.11% versus 9.09%) (Figure 2e). In gut samples, from the stacked bar plot of four predominant phyla, it was identified that the HCC group had a considerably greater average Proteobacteria composition (20.22%) than the control group (16.15%) (Figure 2g). The top 30 genera were plotted at the genus level, while the rest were merged as Others (Figure 2f–h). The box plots generated certain statistical general that comprised 10 and 21 abundant taxa (abundance>1%, Bonferroni-adjusted p < 0.05) in the oral and fecal samples correspondingly (Figure 2i,j). Furthermore, LEfSe analysis was used to identify the bacterial differences in the cladogram (Figure 2k,l) with the criteria of LDA score>3 (Bonferroni-adjusted p < 0.05). Particularly, Streptococcus was more abundant in HCC than that in the control group in both samples. Similarly, the microbiota diversity and taxonomic profiling of the retrospective cohort were analyzed and manifested in Supplemental Figure S1. In short, the microbiota profiles alter in HCC and differ from those in healthy conditions, proposing a microbial transformation of the disease.

Figure 2. The microbial composition and difference between HCC and the control groups in the prospective study.

(a-d) The distributional difference of oral (a-b) and gut (c-d) microbiota was assessed using principal component analysis (PCA) and principal coordinates analysis (PCoA) with the weighted Unifrac index. P value was calculated by the analysis of similarities (ANOSIM). Red and blue plots represented the control and HCC samples.

(e-h) Stacked bar plot of mean proportions of oral (e-f) and gut (g-h) derived taxonomic composition between HCC and the control groups at phylum and genus levels.

(i-j) Statistically differential genera of oral and gut microbiota were evaluated with box plots (Bonferroni-adjusted P < 0.05, mean abundance>1%). The “*” before the taxa represented the unclassified genus.

(k-l) Predominant taxa distribution between groups in a phylogenetic tree with cladogram computed by linear discriminant analysis effect size (LEfSe) algorithm (LDA [log10 transformation] > 3, Bonferroni-adjusted P < 0.05). The circles radiating from inside to outside represented the classification level from Kingdom to Species. Each node on different levels represented the specific taxon, with the diameter corresponding to the relative abundance. Taxa with significant differences (LDA>3, Bonferroni-adjusted P < 0.05) were highlighted and labeled between HCC and the control groups with red and green nodes. Significant taxa were comprehensively described in another panel.

Gut microbial alterations during the transition from healthy status to cirrhosis and HCC

Well-known factors for HCC, HBV and cirrhosis play critical roles in the disruption of microbial composition. Therefore, HBV individuals with and without cirrhosis were enrolled in our prospective study to investigate the changes of intestinal microbiota in the progression from healthy condition to HCC with the taxonomic resolution of 16S rRNA gene amplicons. Among the control, non-cirrhosis, cirrhosis, and HCC groups, the bubble plot and heatmap revealed that the top 5 phyla and top 30 genera were shared across the four groups (Figure 3a,b). The Sankey plot, in particular, was constructed to exhibit the branch association between phylum and genus degree of taxa with different colors (Figure 3c). And the size of the bubble and the breadth of the branch represented the relative richness of each genus among the groups. The average composition of Streptococcus, Subdoligranulum, Prevotella_9, Agathobacter, Faecalibacterium, Bifidobacterium, Ruminococcus]_gnavus_group, Klebsiella, Escherichia-Shigella, and Bacteroides showed the predominant abundance at the genus level, with Firmicutes accounted for the majority at the phylum level. Statistically, the control individuals had the highest abundance of Faecalibacterium (9.12%), Bacteroides (5.65%), Prevotella_9 (3.16%), and Agathobacter (2.94%), and the lowest abundance of Streptococcus (5.67%) and Escherichia-Shigella (9.55%) among the groups. And the HCC group had the highest proportion of Ruminococcus]_gnavus_group (4.08%) and the lowest proportion of Subdoligranulum (1.73%). Of note, Streptococcus (14.14% in HCC, 14.95% in cirrhosis) and Escherichia-Shigella (11.5% in HCC, 13.21% in cirrhosis) were considerably more abundant in HCC and cirrhosis groups than those in the control group (Streptococcus 5.67%, Escherichia-Shigella 9.55%), accompanied by an ascending trend during disease progression. The average abundance of Agathobacter (1.58%) in HCC and Prevotella_9 (0.27%) in the cirrhosis group, on the other hand, was significantly decreased. Additionally, we attempted to characterize the microbiome phenotypes at the organism level with the BugBase tool that incorporated three potential characteristics, including anaerobic, contains mobile elements and facultatively anaerobic. Consequently, the scatter plot depicted an apparent decrease in anaerobic (healthy individuals: 87.23%, HBV-Noncirrhosis: 81.11%, HBV-cirrhosis: 71.51%, HBV-HCC: 73.65%; p = 0.0047), and an increase in mobile elements (healthy individuals: 2.36%, HBV-Noncirrhosis: 6.23%, HBV-cirrhosis: 16.40%, HBV-HCC: 16.28%; p = 0.0001), and facultatively anaerobic (healthy individuals: 5.74%, HBV-Noncirrhosis: 9.33%, HBV-cirrhosis: 18.44%, HBV-HCC: 20.54%; p = 1.425594e-05) throughout the disease evolution (Figure 3d–f), which might indirectly reflect the changes in the intestinal epithelial hypoxia as the disease progressed. Moreover, the altered intestinal permeability caused by epithelial hypoxia was evaluated between the control and HCC groups. As illuminated in Figure 3g–i, the HCC group showed enhanced intestinal permeability than the control group indicated with significantly elevated levels of serum LPS [(6.678 ± 0.78) EU/ml versus (4.025 ± 0.57) EU/ml, p = 0.0083], ZO-1 [(147.60 ± 26.26) pg/ml versus (38.67 ± 8.30) pg/ml, p = 0.0002] and fecal calprotectin [(1.84 ± 0.17) pg/mg versus (0.57 ± 0.05) pg/mg, p < 0.0001]. Taken together, these data showed the landscape of intestinal hypoxia and microbial alterations during hepatic carcinogenesis.

Figure 3. Dynamic microbial alteration in the HCC transition from healthy condition (fecal samples obtained from control group, non_cirrhosis group, cirrhosis group and HCC group) of the prospective cohort.

(a) Bubble plot of the gut genus affiliation to phylum with different colors and the abundance of genera with bubble size.

(b) Heatmap of different gut genera among the control, non-cirrhosis, cirrhosis, and HCC group. Color in the heatmap was utilized to describe specific general abundance within the four groups, with blue indicating lower abundance and red indicating higher abundance.

(c) Sankey plot of the taxonomic data changes with the breadth of the branch at genus (right side) and phylum (middle) levels during the disease progression (left side). The color and width of the branches represented the flow of specific genera within different phyla.

(d-f) Characteristic analysis of fecal microbiota phenotype among the four groups with Anaerobic (d), Contains mobile elements (e) and Facultatively anaerobic (f).

(g) Metagenomic sequencing validation of the mobile genetic element database between HCC (n = 42) and Control (n = 24) groups

(h-j) The levels of serum Zonulin-1 (g), LPS (h), and fecal calprotectin (i) between the control and HCC groups, with the data representing mean ± SEM, *P<0.05, **P<0.01, ***P<0.001, ****P<0.0001; ns, no significance.

Microbial relevance with clinical features and functional profiling between the groups

After recognition of the difference between HCC and the control group, SparCC network analysis was used to clarify the interactive correlation of the prevalent genera with each other (Figure 4A–B). Relevantly, Streptococcus was identified in both oral and intestinal correlation networks, where it exhibited positive and negative connections with a total of 11 genera. According to the 16s rRNA gene amplicons, the underlying relationship between clinical features and fecal microbiome was detected with the RDA plot in which the distribution of microbiome was positively related with MELD, Child-Pugh, Age, CCI, AFP, AST, ALT, GGT, Cre, LDL, Tchol, PT and INR, and negatively related with ALB, BMI, HDL, TG, and PLT in HCC patients. Streptococcus and Veillonella showed a consistent tendency whereas Faecalibacterium, Akkermansia, Ruminococcaceae_UCG−002, Rumin-ococcaceae_UCG−005, Ruminococcus_2, and Subdoligranulum showed the contrary (Figure 4c–e). More importantly, for direct identification of the functional profiles, a number of backup gut samples from the HCC (n = 42) and control (n = 34) groups in the prospective cohort were further analyzed by the shotgun metagenomic sequencing. Methodologically, the Diamond algorithm was constructed to annotate and calculate the statistical difference in the KEGG pathways between the two groups (Figure 5, adjusted p < 0.05). In the majority of the top 20 statistical KEGG items, the HCC group showed more activities than the control group, including Arginine and proline metabolism; Alanine, aspartate, and glutamate metabolism; Cysteine and methionine metabolism; Arginine biosynthesis; Valine, leucine, and isoleucine biosynthesis, which underscored the significance of secondary metabolites, namely the amino acid synthesis and metabolism in the HCC microbiome.

Figure 4. Microbial biomarkers (LDA>3, Bonferroni-adjusted P value<0.05) intercorrelation and relevance with clinical characteristics between HCC and the control groups in the prospective cohort.

(a-b) SparCC network of the abundance of genera with node size and the link between nodes with a correlation coefficient |rho| more than 0.1 by default in oral (a) and fecal (b) samples, with positive and negative correlations represented by solid and dashed lines, respectively (correlation p < 0.05).

(c-e) RDA analysis of clinical feature and fecal genus in HCC and control groups, with the acute and obtuse angle indicating the positive and negative association between specific genus and clinical indicator.

Figure 5. The major KEGG pathways between HCC and the control groups of the prospective cohort with the metagenomic sequencing data of fecal samples.

Differential shotgun metagenomic sequence-based KEGG pathways in gut microbiota between the two groups detected by Diamond software. The top 20 items are listed along with the appropriate 95% confidence intervals and adjusted-p values.

The performance of a combined oral and fecal microbiota-based model in the diagnosis of HCC

Having recognized the existence of a microbial signal for HCC at the coarse level of general community composition based on the 16S rRNA gene sequencing, the random forest method was employed to visualize specific taxa that contributed to the diagnostic potential between HCC and the control group. In the retrospective study, the top 20 influential taxa were identified in oral or fecal microbiota, as represented by the mean decrease accuracy of each taxon in Figure 6a–b. After qualification, 10 (Anderoglobus, Corynebacterium, Prevotellaceae_UCG-001, Pseudoramibacter, Atopobi-um, Lutispora, Cryptobacterium, Dialister, Lachnospiraceae_NK4A136_group, Sphaerochaeta) and 9 (Elizabethkingia, Burkholderia_Caballeron-ia_Paraburkholderia, Klebsiella, Delftia, Faecalib-aculum, Acetatifactor, Lactobacillus, Rumino-coccaceae_UCG-010, Stenotrophomonas) genera with disease-specific univariate (Figure 1b) were applied to establish oral and fecal models with favorable ROC values of 0.7971 and 0.8084 in the prospective cohort for validation, which was comparable to the accuracy of AFP signature (AUC = 0.8505, Figure 6c). Taken together, the oral and fecal influential taxa were merged to develop a combined prediction model. Afterward, the classification accuracy of such a prediction model between HCC and the controls was verified with an AUC of 0.9405, indicating its potent differential diagnosis capability in the prospective population. Combined with AFP, the performance of the diagnostic model further improved to 0.9811. Generally, our results explained the diagnostic potential of the microbiota and established related classifiers with favorable accuracy.

Figure 6. Establishment of classification models from the retrospective data and validation for distinguishing HCC from the control in the prospective cohort.

(a-b) The top 20 influential oral (a) and gut (b) genera for distinguishing HCC from the control were identified by random forest analysis from the retrospective data. Each genus is ranked with the mean decrease in accuracy. The mean relative abundance (%) of the predictive genus derived from the oral (n = 10) and fecal (n = 9) models was described as Control vs HCC between the two groups of the prospective data.

Abbreviations: B_C_P (Burkholderia_Caballeronia_Paraburkholderia), R_UCG_010 (Ruminococcaceae_UCG_010), L_N_group (Lachnospiraceae_NK4A136_group), P_UCG_001 (Prevotellaceae_UCG_001).

(c) External longitudinal-validation results of the AFP model, oral model, fecal model, oral-fecal model, and the combination model of AFP with microbial features in the prospective cohort were shown as ROC curves with corresponding colors.

Associations between the oral, intestinal, and hepatic microbiomes

The biliary duct anatomically interacts with the duodenum, which proposed an evident link for the microbiota to colonize the hepatic tissues and facilitate HCC development. To corroborate the presence of several influential taxa described in the prospective study, we taxonomically profiled the biopsies of HCC tumors (n = 46), normal hepatic tissues from paracancerous regions (n = 42), and hepatic hemangioma (n = 11) with 5 R 16S rRNA gene sequencing, following a strict procedure to eliminate putative bacterial contaminants. A relatively rich and diverse hepatic microbiome with the top 30 genera and corresponding 16 phyla were observed in over 50% of the samples (Figure 7a–b). Among these, a series of genera including Enterococcus, Escherichia-Shigella, Faecalibacterium, Streptococcus, Delftia, Bifidobacterium, Prevotella, Parabacteroides, and Stenotrophomonas were abundant in tumors (abundance>1%), along with their presence in the oral or gut samples of HCC patients. In particular, among the differential taxa in oral and fecal samples, the species Streptococcus parasanguinis was identified in tumors with a relatively low abundance. We then characterized the microbes using IHC assays with LPS (for gram-negative bacteria) and LTA (for gram-positive bacteria) staining, which showed the presence of microbiota as punctate dots (Figure 7c). FISH assays with specific primers were utilized to further verify the prevalence of Streptococcus parasanguinis, Streptococcus mitis, Streptococcus salivarius, Delftia acidovorans, Parabacteroides distasonis, and Stenotrophomonas maltophilia in tumor tissues. Owing to the decreased sensitivity of FISH, Stenotrophomonas maltophilia, Parabacteroides distasonis, Streptoco-ccus salivarius, and Streptococcus mitis were detectable in the tested samples (Figure 7d). In summary, the concordant amplicon, IHC, and FISH data suggested that distinct hepatic microbial residents may be derived from the oral and the gut.

Figure 7. Presence of microbiota in hepatic tissues with different conditions of the prospective cohort.

(a-b) Microbial composition in hepatic tumor and healthy tissue samples at phylum (a) and genus (b) levels, with top 16 phyla and 30 genera as determined by 5 R 16S amplicon data, with different colors representing the corresponding taxa. The matched samples were put next to each other (cancerous and paracancerous tissues) within the gray bar, along with benign liver tissues from the hepatic hemangioma inpatients. Two random regions of each large tissue were selected for microbial analysis, e.g. CA_6-1, CA_6-2. The blank bars of Germ-free samples were marked as GF.

(c) IHC of LPS and LTA staining in hepatic cancer, paracancerous tissues, and control samples as indicated by punctate dots (magnification x400). The expression of LPS and LTA between tumor and the paracancerous group showed no significant difference, *P < 0.05, **P < 0.01, **P < 0.01, ***P < 0.001, ****P < 0.0001; ns, no significance.

(d) Representative microscope images for Parabacteroides distasonis, Streptococcus mitis, Streptococcus salivarius, and Stenotrophomonas maltophilia with Fam and Cy3 fluorescent dyes, and DAPI was used for nucleus staining in tumor and control samples (GF tissues) (magnification x400). The expression of fluorescence density between tumor and control samples showed significant difference, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001; ns, no significance.

Discussion

As an insidious malignancy, HCC remains a formidable challenge that requires extensive efforts for early detection to mitigate its adverse burden. Herein, this was the first time we determined that the oral microbiota of HBV-HCC patients differed significantly from that of healthy individuals, as well as a similar gut microbial trend that was confirmed in the literature. These alterations thus prompted several distinct microbiota-based classifiers to predict HCC incidence. Importantly, microbial data of oral and fecal counterparts in the matching prospective cohort confirmed the diagnostic superiority of a combined classifier with an impressive AUC value. Furthermore, we identified the presence of hepatic microbiota and disclosed a potential relationship across microbes in the cancer ecosystem, mouth, and intestine.

Previous studies have explored associations between digestive tract diseases and oral¹⁹ or gut^20,21 microbiome with the taxonomic resolution of 16S rRNA gene amplicons, including early liver carcinoma with cirrhosis⁹. In particular, the gut microbiome has also been identified to effectively predict HCC incidence¹⁰, demonstrating its diagnostic potential for classification. In our retrospective study, we delineated the microbial composition of oral and fecal samples with 16S rRNA gene sequencing in HCC and healthy individuals. The taxa heterogeneity and distribution described with β diversity and the LEfSe algorithm (Supplemental Figure S1) significantly distinguished HCC from the control group. Subsequently, with the random forest algorithm, several influential oral (n = 10) and gut (n = 9) genera were recognized to improve the diagnostic efficiency of HCC patients from healthy individuals. For validation, both tongue and matching fecal samples were prospectively collected from 124 HCC patients and 91 control individuals, with taxonomic profiles showing consistent results. Respectively, we developed amplicon-derived classifiers that could accurately predict HCC based on characteristic oral and fecal microbial genera with AUC values of 0.7971 and 0.8084 (Figure 6c), which showed comparable performance to AFP (AUC = 0.8505), as the sole FDA-approved marker for HCC²². Moreover, the oral-gut-liver axis has been introduced to play a pivotal role in the maintenance of metabolism and the pathophysiology of hepatic diseases through blood circulation and enteral routine^23,24, which was confirmed with the intrinsic connection among the microbes and their correlation to the clinical markers in SparCC and RDA analyses (Figure 4). Reportedly, a combination of oral and fecal microbiota has been proposed to enhance the sensitivity of predictive models for pancreatic ductal adenocarcinoma and colorectal cancer^25,26. Therefore, the optimal 19 microbial biomarkers were merged to establish a distinct microbiota classification that considerably increased the accuracy to 0.9405 (Figure 6c). Importantly, when paired with the AFP level, the combined prediction model’s accuracy further improved to 0.9811 (Figure 6c), suggesting its complementary effect to the serum testable markers. These results indicated that oral and gut microbiota-targeted markers may become promising noninvasive tools for early diagnosis, demonstrating for the first time that an integrated analysis of the microbiome from subjects with HCC was conducted using oral and fecal samples.

To further explore when the alterations of microbiota occur in the participants, HBV and HBV-cirrhosis patients were enrolled and the relative abundance of potentially pathogenic bacteria was also evaluated to indirectly reflect the difference in microbial composition across HCC disease stages. Consequently, a couple of genera were identified to be crucial elements, some of which had been confirmed in the occurrence and progression of HCC, including Streptococcus ²⁷, Prevotella_9 ²⁸, Faecalibacterium ²⁹, and Bacteroides ³⁰. Interestingly, these potential pathogens represented by Streptococcus showed a continuous accumulation from the healthy condition to HCC, which was speculated to be related to the successive decline in the gastrointestinal mucosal barrier. Further microbial phenotype analyses showed a decreasing trend in anaerobic bacteria, and an increasing trend in mobile elements and facultatively anaerobic bacteria (Figure 3d–f), indicating that disruption of epithelial hypoxia and accumulation of harmful bacteria were linked to gut dysbiosis and inflammatory reactions that increased intestinal permeability³¹ for intestinal bacterial translocation during carcinogenesis. In addition, improved intestinal permeability was also confirmed with elevated levels of serum LPS, ZO-1 and fecal calprotectin in HCC patients compared with the control individuals (Figure 3g–i). Consistently, as a typical facultative anaerobic genus, Streptococcus displayed a continuous accumulation from healthy conditions to HCC in the Sankey plot (abundance: 5.67% to 14.14%, Figure 3c), signifying its potential relationship with HCC. Notably, emerging lines of evidence have linked microbiota metabolites and dysfunctional gut barriers to the occurrence of HCC. A variety of metabolites, including choline, bile acids (BAs), short-chain fatty acids (SCFAs) and amino acids have been reported to serve as critical signaling factors and metabolic substrates that affect liver function³². Accordingly, genera Streptococcus and Lactobacillus were enriched in HCC, whereas the Akkermansia, Prevotella_2, Subdoligranulum and Faecalibacterium decreased, all of which were associated with BAs synthesis, especially with lithocholic acid (LCA) and deoxycholic acid (DCA)^33–35. Intriguingly, Faecalibacterium was involved in the production of butyrate³⁶ and acted as a SCFAs-related bacterium that modulated inflammatory reactions. Meanwhile, KEGG pathway analyses discovered several biological metabolic pathways in HCC with the metagenomic sequencing (Figure 5), particularly the biosynthesis and metabolism of Arginine³⁷, proline³⁸, alanine, aspartate³⁹, glutamate⁴⁰, cysteine and methionine⁴¹, which may provide a promising prospect for the regulatory role of gut microbes in the biosynthesis and metabolism of amino acids during HCC development with metagenomic sequencing. With regard to the function of the gut barrier, of those exactly defined genera that were differentially distributed in control and HCC participants, alterations in the relative abundance of Akkermansia, Bacteroides and Lactobacillus were related to intestinal permeability^42,43 that contributed to the gut microbiota shifting for gut-liver reactions. Besides, dysbiosis of Bacteroides and Bifidobacterium in the tumor lesion microbiota was correlated with cirrhosis and HBV-related HCC^11,44. It is known that the leaky gut and bacterial translocation could facilitate microbial metabolites to the liver, followed by the impairment of BAs metabolism and inflammation via Toll-like receptors (TLRs) signaling, eventually leading to the progression of tumorigenesis⁴⁵. Additionally, the present study has also shown a positive relationship between the gut microbiome and clinical outcomes in patients with HBV-related HCC⁴⁶. Therefore, all this evidence underscored the importance of certain signature-derived microbiota in the development of liver disease and early detection of HCC, which reminded us that the intervention of gut microbiota in the precancerous stage may slowdown the deterioration among the patients.

Clinically, cohort studies have revealed that features of tissue-derived bacteria are correlated with cancer risks⁴⁷, treatment response,¹¹ and cancer prognosis⁴⁸ through various mechanisms, including increased mutation rate, modulation of oncogenic genes and pathways, and regulation of tumor immune microenvironment^49–51. In particular, a landmark work reported by Cai and his colleagues⁵² innovatively identified the intracellular microbiota that could protect circulating tumor cells from mechanical stress by reorganizing the actin cytoskeleton in breast cancer. Alternatively, bacteria had also been postulated to contribute to tumorigenesis through indirect mechanisms such as metabolites, biofilms, inflammation, and immunosuppression⁵³. In our study, after strictly aseptic collection and careful filtration of the 16S rRNA gene amplicon datasets with different controls to account for contamination during the procedures, several taxa could be traced from the gut and hepatic tissues with univariate enrichment in lesions, suggesting possible associations between HCC and intestinal microbiota. Based on specific antibodies and primers, IHC and FISH assays were used to describe the landscape of microbiome and certain abundant species in hepatic samples. Similarly, a significantly elevated level of oral-intestinal strain transmission was observed in patients with HCC, particularly of signature taxa, indicating that they originated intraindividually from the oral cavity. These findings illustrated that the oral, gut, and hepatic microbiomes may be intricately linked, which will be essential to disclose their respective roles in HCC etiology.

To our knowledge, it was the first longitudinal cohort study focusing on the diagnostic potential of the combined oral-gut microbiome in HBV-HCC patients, and the microbial relationship between liver tumors and the oral-gut axis had been systematically introduced through the prospectively matched population. A majority of HCC in our country was caused by HBV and developed into cirrhosis, which would influence the dynamic microbiome along with the predominant geographical variances^54–56. Hence, with the relatively large size of HBV-HCC patients, our co-territorial study achieved a few convincing results with internal validation and less geographical bias. However, there were certain limitations to our research. On the basis of differential abundance analysis, more robust methods including ANCOM-BC⁵⁷, LinDA,⁵⁸ and LOCOM⁵⁹ were required with a unique focus on the compositional nature of the original data and powerful performance in the FDR control. Sequential samples were obtained from a specific disease status without dynamic samples of the same patient over time, and the data were derived from the same race in a single center (Asia), therefore future study was required to verify the general application for the predictive signature in different races from other regions (e.g. North America, Africa). And it was innovative to demonstrate the visible presence and exact location of viable bacteria within tumor tissues under the corresponding medium, which may assist in clarifying its potential carcinogenic role in additional works. Nevertheless, the present study described a distinct microbiome signature that enabled robust performance for HCC detection with unique specificity, complementary to current biomarkers, and the possibility of effective HCC screening and monitoring. Beyond the use for early diagnosis, we believe that the panel of microbial species may be relevant to hepatic carcinogenesis, providing promising insights into proper HCC prevention and intervention.

List of abbreviations

HCC	=	hepatocellular carcinoma
AUC	=	area under the curve
AFP	=	alpha-fetoprotein
FISH	=	fluorescence in situ hybridization
OSCC	=	oral squamous cell carcinoma
IBD	=	inflammatory bowel disease
PCA	=	principal component analysis
PCoA	=	principal coordinates analysis
KEGG	=	Kyoto Encyclopedia of Genes and Genomes
STAMP	=	statistical analysis of metagenomic profiles
SEM	=	standard error of mean
LEfSe	=	linear discriminant analysis effect size
FDR	=	false discovery rate
RDA	=	redundancy analysis
ROC	=	receiver operating curve
LPS	=	lipopolysaccharide
ZO-1	=	zonulin-1
ELISA	=	enzyme-link immunosorbent assay
LAL	=	limulus amebocyte lysate
IHC	=	immunohistochemistry
LTA	=	lipoteichoic acid
RT	=	room temperature
GF	=	Germ-free
PPIs	=	proton-pump inhibitors
NSAIDs	=	nonsteroidal anti-inflammatory drugs
AAPs	=	atypical antipsychotics
ARB	=	angiotensin II receptor blocker
ACEI	=	angiotensin converting enzyme inhibitors
ANOSIM	=	analysis of similarities
BMI	=	body mass index
CCI	=	charlson comorbidity index
MELD	=	model for end-stage liver disease
HBV	=	hepatitis B virus
INR	=	international normalized ratio
PT	=	prothrombin time
HDL	=	high density lipoprotein
LDL	=	low density lipoprotein
ALT	=	alanine transaminase
AST	=	aspartate transaminase
PLT	=	platelet
Tchol	=	total cholesterol
TG	=	triglyceride
Tbil	=	total bilirubin
Cre	=	creatinine
Alb	=	albumin
GGT	=	γ-glutamyl transpeptidase
BAs	=	bile acids
SCFAs	=	short-chain fatty acids
LCA	=	lithocholic acid
DCA	=	deoxycholic acid
TLRs	=	toll-like receptors

Acknowledgments

The authors would like to thank Peng Wu, Qingqing Liang, Wei Wang, and Gaoyang He of Lc-Bio Technologies Co., Ltd (Hangzhou, China) for their technical support in this work. We thank the anonymous reviewers and the editors for their helpful remarks and useful feedback that improved this paper.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The data that support the findings of this study are available on request from the corresponding author, Gang Chen. The data are not publicly available due to their containing information that could compromise the privacy of research participants.

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/19490976.2023.2201159

Additional information

Funding

The work was supported by the National Natural Science Foundation of China [81772628, 81703310, 82072685]; the Research Foundation of the National Health Commission of China–Major Medical and Health Technology Project for Zhejiang Province [WKJ-ZJ-1706]