Introduction

Advanced hepatocellular carcinoma (HCC) is an intractable global health concern associated with a poor prognosis, and it is the sixth leading cause of cancer-related mortality worldwide.1 Recent noteworthy advancements in immune checkpoint inhibitor (ICI) therapies targeting the programmed death 1 (PD-1) and programmed death ligand 1 (PD-L1) axes have demonstrated clinical efficacy in treating advanced HCC, both as first-line and second-line treatments.2 However, the response rates are heterogeneous and non-robust and a proportion of patients with advanced HCC fail to respond to this therapy or experience side effects.2 Biomarkers such as PD-L1 expression, high tumor mutation burden, high microsatellite instability, microbiota, and metabolic profiles are used to select patients who are likely to show a positive response to ICI treatment in HCC and other cancers, and studies suggest that gut microbiome and metabolite properties are of great value.3–7

Recent studies have revealed that the gut microbiota mediates diverse responses to immunotherapy and may serve as a potential predictive biomarker.8 Both preclinical and clinical research studies have provided growing evidence that gut microbiota and metabolites can modulate the oncogenesis of cancer and shape antitumor immune responses, affecting the efficacy and toxicity of immunotherapies, particularly ICIs.8 9 The human gut harbors more than 103 microbes, including the multi-kingdom microbes, prokaryotes, viruses, protists, and fungi, which play vital roles in intestinal ecology and human health.3 10 However, current research in the field of immunotherapy primarily focuses on the bacterial components of gut ecology.3 The intricate associations between other microorganisms, particularly fungi or the ‘‘mycobiome”, and ICI, remain unexplored. Despite constituting less than 0.1% of the gut flora, the mycobiome plays an indispensable role in the human gut microbiota and is tightly integrated into the bacterial community. Extensive evidence has shown its involvement in diverse disease pathogenesis and its profound influence on the host immune system, including antitumor responses.11–14 Gut fungal and bacterial communities are intrinsically interconnected and likely impact each other in various ways and are intertwined with host metabolic pathways.10 14 Growing evidence supports a close relationship between gut microbiota and various metabolic pathways, ultimately modulating the efficacy of ICIs through the regulation of gut bacterial functions and metabolites across different tumor types. For instance, microbiota-derived short-chain fatty acids (SCFAs) and bile acids (BA) are key mediators of this process.3 15

Previous studies, including ours, have provided insights into the association between specific gut bacteria and metabolites, and the efficacy of anti-PD-1 therapy in HCC. However, the findings vary significantly.4–7 Therefore, it remains imperative to further investigate the characteristics of the gut microbiota in advanced HCC and to elucidate the potential roles of the gut microbiota and metabolites in predicting responses to ICI treatment. In this study, we performed a comprehensive analysis of metagenomic, ITS2, and metabolomic data sets to assess the characteristics of multi-kingdom microbiota, including the gut microbiome, mycobiome, and metabolome in ICI-treated HCC (online supplemental figure S1 and table S1).

MethodsStudy cohort and fecal sample collection

We prospectively included 188 patients with advanced HCC from three registered clinical trials (NCT03892577, NCT04010071, and NCT03895970) conducted at Peking Union Medical College Hospital (PUMCH), from November 2018 to April 2022. The patients included in the study met the following criteria: pathological or clinical diagnosis of HCC; treatment with a PD-1/PD-L1 inhibitor, with or without targeted therapy; and other criteria, including adequate function of major organs, an Eastern Cooperative Oncology Group (ECOG) performance status (PS) of 0–2, and radiologically-evaluable lesions. Patients who had used other experimental drugs within 4 weeks prior to the start of the study, received antibiotics before or during treatment, underwent invasive biliary tract procedures, had organ transplants, or had autoimmune diseases were excluded. A total of 108 patients were excluded for these reasons and 80 patients were enrolled in this study (online supplemental figures S1 and S2).

The PD-1/PD-L1 inhibitor was administered intravenously every 3 weeks at the indicated dose for each drug. The treatment regimen was terminated when disease progression or intolerable toxicity was observed. Clinical response was assessed every 6–8 weeks according to Response Evaluation Criteria in Solid Tumours (RECIST) V.1.1 criteria.16 Patients with efficacy assessed as continuous partial response (PR), complete response (CR), or stable disease (SD) ≥6 months, were defined as having achieved durable clinical benefit (DCB), while disease progression (PD) or continuous SD <6 months was defined as a non-durable clinical benefit (NDB).17 Immune-related adverse events (irAEs) were assessed and monitored for each patient from the start of ICI treatment to the final follow-up time. IrAEs were diagnosed according to the National Comprehensive Cancer Network guidelines on ICI-related toxicities.18 All irAEs were co-confirmed by at least two or more oncologist experts and excluded the possibility of other adverse events caused by targeted therapy. Fecal samples were dynamically collected at PUMCH within 2 weeks before the first treatment and before every second cycle of PD-1/PD-L1 inhibitor administration. All fresh fecal samples were stored in sterile containers and immediately stored in a −80°C refrigerator.

Metagenomic sequencing and bioinformatic analysis

Genomic DNA was extracted using a Feces Genomic DNA Purification Kit (BIOER, Binjiang, China). DNA concentration and purity were measured simultaneously using a Qubit V.3.0 (Thermo Fisher Scientific, Waltham, Massachusetts, USA) and NanoDrop One (Thermo Fisher Scientific), respectively. Sequencing libraries were generated using NEB Next Ultra DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, Massachusetts, USA) following the manufacturer’s recommendations and index codes were added. The library was sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego, California, USA) and 150 bp paired-end reads were generated. Raw data were processed using fastp (V.0.19.7) to acquire clean data for subsequent analysis.19 Trimmed reads aligned to the Homo sapien genome assembly hg3720 were removed using KneadData integrated Bowtie (https://huttenhower.sph.harvard.edu/kneaddata/) to obtain metagenomic DNA sequences.

The MetaPhlAn tool (V.4.0.3) was used to quantitatively profile the taxonomic composition of the metagenome,21 whereas HUMAnN (V.3.0.0. alpha.4) was used to estimate the microbial metabolic and functional pathways.22 Functional pathway analysis included metabolic pathway (MetaCyc) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functions. Alpha diversity was used to evaluate the complexity of the species diversity for each sample using the Shannon index. Beta diversity calculations were applied to analyze the diversity of samples for species complexity using principal coordinate analysis (PCoA). The p value for the analysis of similarities was obtained using a permutation test. Linear discriminant analysis effect size was determined using the non-parametric factorial Kruskal-Wallis rank sum test and linear discriminant analysis (LDA) to determine the differential taxa between groups. The parameters were set as follows: the default logarithmic LDA score was 2.0 and the default p value was 0.05. Fivefold cross-validation was performed using the least absolute shrinkage and selection operator (LASSO) regression model (R V.4.2.0, glmnet package) using the 25-species which differed significantly between the DCB group (n=25) and NDB group (n=24) with a p value<0.05 in the training cohort. 18 species with minimum errors were chosen as the optimal set. The probability of the test cohort, including the 17 DCB and 14 NDB groups, was calculated using this set of species, and a receiver operating characteristic curve was drawn (R V.4.2.0, pROC package).

Metabolome sequencing and bioinformatic analysis

Fecal samples were analyzed using a liquid chromatography-tandem mass spectrometry (LC-MS/MS) system. LC-MS/MS analyses were performed using a UHPLC system (Vanquish, Thermo Fisher Scientific) with a UPLC BEH (2.1 mm × 100 mm, 1.7 µm) amide column coupled to Q Exactive HFX mass spectrometer (Orbitrap MS, Thermo Fisher Scientific). Metabolites were detected in both the negative and positive ion models. Raw data were converted to the mzXML format using ProteoWizard,23 and an internal standard was used for data normalization. Principal component analysis was performed using SIMCA software. To further investigate the differences between the groups, orthogonal partial least squares-discriminant analysis (OPLS-DA) modeling was performed, and the quality of the model was tested using sevenfold cross-validation. Student’s t-test p<0.05) was used to screen differential metabolites, and the variable importance in the projection of the first principal component of the OPLS-DA model was greater than 1. KEGG (http://www.genome.jp/kegg/) and PubChem (https://pubchem.ncbi.nlm.nih.gov/) were used to annotate the metabolites.

ITS2 data analysis

The DNA was extracted using a DNA extraction kit for the corresponding sample. The concentration and purity were measured using a NanoDrop One (Thermo Fisher Scientific). The ITS2 hypervariable regions of the fungal ITS rRNA gene were amplified with the specific primers ITS3F: GCATCGATGAAGAACGCAGC and ITS4R: TCCTCCGCTTATTGATATGC using a Bio-Rad S1000 PCR system (Bio-Rad Laboratories, Hercules, California, USA).24 Raw data were demultiplexed using the MiSeq Controller Software (Illumina). QIIME2 (V.2020.2) was used for downstream analysis.25 Demultiplexed ITS2 sequences were denoised and grouped into amplicon sequence variants (ASVs) using DADA2.26 ASV features present in only one sample were excluded from all cohorts. Individual ASVs were taxonomically classified based on the UNITE (V.8.2) database using the VSEARCH tool wrapped in QIIME2. Alpha diversity was estimated using the Shannon index at the operational taxonomic unit (OTU) level. PCoA was performed based on the unweighted UniFrac distance. FunGuild was employed for the functional prediction of fungi.27

Statistical analysis and correlation network

Differences in clinical characteristics were evaluated using the Pearson’s χ2 test or Fisher’s exact test. Permutational multivariate analysis of variance (permutations=9,999, R V.4.2.0, vegan package) was performed to investigate the effect of each clinical characteristic based on the Bray-Curtis distance matrix of the species abundance profiles of the samples. A canonical correlations analysis (CCA) (R V.4.2.0, vegan package) was used to assess the effects of each clinical characteristic on the species abundance profile of the samples. Spearman’s correlation coefficients between differentially enriched features were determined. Patients were categorized into high-ambudance and low-abundance groups based on the mean values of the relative abundances of different taxa. Survival-related taxa can be categorized as either beneficial for survival (high-abundance in DCB group, or low-abundance in NDB group) or not conducive to survival (low-abundance in DCB group, or high-abundance in NDB group). Finally, the survival benefits were categorized based on the combined expression of gut bacterial and metabolite profiles. A group exhibiting both benefits was classified as the good group, a group with no benefits was categorized as the poor group, and a group displaying one benefit and one lack of benefit was classified as the moderate group.7 Survival curves were estimated using the Kaplan-Meier method and compared using the log-rank test. The follow-up time was estimated using the reverse Kaplan-Meier method. Progression-free survival (PFS) was defined as the time from the initiation of PD-1/PD-L1 inhibitor treatment to tumor progression or death. Overall survival (OS) was defined as the total time from the initiation of PD-1/PD-L1 inhibitor treatment until death from any cause. Univariate and multivariate analyses were performed using a Cox proportional hazards model. Unless otherwise stated, all statistical analyses were performed using R software (V.4.2.0), and p values<0.05 were considered statistically significant.

ResultsCharacteristics of study cohort

80 patients with HCC (mean (SD) age, 57.4 (10.4) years; 72 (90.0%) men) were enrolled in this study (online supplemental figure S2). The study was followed-up until August 2023, and the median follow-up time was 27.7 (95% CI: 23.5 to 33.4) months. Patient characteristics of the training and test cohorts are shown in online supplemental tables S2 and S3, respectively. Hepatitis, including chronic hepatitis B and C, was not significantly different between the DCB and NDB groups in both train and test cohorts. 25 (49.0%) patients in the train cohort developed an irAE, and there were only four Grade ≥3 irAEs, of which three were dermatological irAE and one was pneumonitis irAE (online supplemental table S2). All patients received ICI plus molecular targeted therapy, with the highest use of ICI for toripalimab (N=36, 45.0%) and camrelizumab (N=31, 38.7%), and targeted drugs for lenvatinib (N=62, 77.5%) (online supplemental table S4).

Diversity and composition of gut microbiome and mycobiome

At baseline, there were no significant differences in the alpha diversity of bacteria (species level) or fungi (OTU level) between the DCB and NDB groups (online supplemental figure S3A,D). Beta diversity in the DCB group was significantly higher in the DCB group than that in the NDB group for both bacteria (species level) and fungi (OTU level) (online supplemental figure S3B,E). The overall difference between the DCB and NDB groups was significant at the bacterial level, whereas the difference at the fungal level was relatively insignificant (online supplemental figure S3C,F). At the bacterial species level, no significant changes in alpha diversity were observed following ICI treatment (online supplemental figure S4A). However, after 6–8 weeks (point 2), the difference in bacterial beta diversity (species level) gradually disappeared between the two groups, suggesting that immunotherapy had a greater effect on the distribution of intestinal bacteria (online supplemental figure S4B). At the fungal level, there was an overall decreasing trend and then a gradual recovery of alpha diversity with the use of ICI. Differences between the DCB and NDB groups appeared after 6–8 weeks (point 2), after which they disappeared, and beta diversity was relatively stable (online supplemental figure S4C,D).

Based on multiomics data, 894 bacteria (based on species), 1,906 fungi (based on OTU), and 620 metabolites were detected (online supplemental figure S5A,E and table S5). The relative abundances of the DCB and NDB groups were compared at different phylogenetic levels to investigate the relationship between the composition of gut bacteria and fungi and ICI treatment. At the phylum level, the predominant bacterial (Firmicutes and Bacteroidetes) and fungal compositions (k_Fungi; p_Ascomycota and k_Fungi; p_Basidiomycota) were consistent in both the DCB and NDB groups (online supplemental figure S5B,F). At the species level, the bacteria with the highest percentages in the DCB and NDB groups were s_Escherichia_coli and s_Faecalibacterium_prausnitzii, respectively, and the fungi with the highest percentages in both groups were k_Fungi; s_Candida_albicans (online supplemental figure S5D,H). The gut microbiome and fungal composition at the phylum, genus, and species levels in the DCB and NDB groups are shown in online supplemental figure S5. After immunotherapy, in both the DCB and NDB groups, the composition of bacteria and fungi was relatively stable at dynamic points at the phylum level, whereas at the species level, the composition and percentage of bacteria and fungi were relatively more variable at dynamic points (online supplemental figure S6).

Significantly different gut bacteria, fungi, and metabolites in ICI-treated HCC

We identified 25 enriched species, 18 enriched fungi, and 37 enriched metabolites in the DCB and NDB groups (figure 1A,B; online supplemental table S5). We observed that certain gut bacteria or fungi, which were enriched in the DCB group, had a lower frequency of detection or were even undetectable in the NDB group, such as s_Phascolarctobacterium_faecium (76.0% in DCB and 37.5% in NDB), s_Candidatus_Avimonas_narfia (20.0% in DCB and not detected in NDB), and k_Fungi; g_Cladosporium (27.3% in NDB and 4.3% in DCB). Similarly, certain gut bacteria or fungi enriched in the NDB group were detected less frequently or even undetectable in the DCB group, such as s_Actinomyces_sp_ICM47 (75.0% in NDB and 48.0% in DCB), s_Senegalimassilia_anaerobia (25.0% in NDB and 4.0% in DCB), s_Faecalibacillus_faecis (16.7% in NDB and not detected in DCB) and k_Fungi; g_Trichoderma (34.8% in NDB and 9.1% in DCB) (figure 1A; online supplemental figure S7A,B). Most of the differentially expressed metabolites were enriched in the DCB group (figure 1B; online supplemental figure S8). Furthermore, a phylogenetic tree was constructed based on the abundance of bacteria, which was mainly distributed in Candidatus_Saccharibacteria in the DCB group and Firmicutes in the NDB group (figure 1C). Based on the fungi, the phylogenetic tree demonstrated that k_Fungi; p_Ascomycota was distributed in the DCB group, and k_Fungi; p_Basidiomycota in the NDB group (online supplemental figure S7D). The distribution differences of bacteria were different between clinical assessment responses, PR, SD, and PD, with PD predominantly Actinobacteria at the phylum level and Actinomyces_sp_ICM47, Senegalimassilia_anaerobia, and Faecalibacillus_faecis at the species level (figure 1D; online supplemental figure S7C).

Figure 1

Significantly different gut bacteria, fungi, and metabolites in immune checkpoint inhibitor-treated hepatocellular carcinoma. (A) Significantly different abundant gut bacteria (species level) and fungi, and detection frequency in the DCB and NDB groups using LEfSe (LDA>2.0, p<0.05). **Represents the species incorporated in the predictive model. (B) Significantly different metabolites and variable importance in the projection between DCB and NDB groups. (C) Taxonomic cladogram from LEfSe showed different bacteria enriched in the DCB and NDB groups (LDA>2.0, p<0.05). (D) Ternary phase diagram for the prominent significantly different abundant species. The different colored dots show the phylum corresponding to species. DCB durable clinical benefit; LDA, linear discriminant analysis; LEfSe, Linear discriminant analysis effect size; NDB, non-durable clinical benefit; PD, disease progression; PR, partial response; SD, stable disease.

To gain insight into the potential interplay between multi-kingdom taxa and their potential roles in the response to HCC immunotherapy, we performed alterations in the multi-kingdom co-abundance network analysis based on the abundance of differential bacterial species, fungi, and metabolites. Generally, the ecological network of patients in the NDB group (75 taxa and 2,775 associations) was more complex than that of patients in the DCB group (69 taxa and 2,346 associations). Apart from strong correlations between intrakingdom taxa, we found substantial associations between interkingdom species, especially between bacteria and metabolites, in both the DCB and NDB groups (online supplemental table S6). After adjusting for association coefficients (r>0.4), we found a higher number of associations in the DCB group than in the NDB group, suggesting that moderately associated taxa (such as bacteria and metabolites) may play a greater role in the immune-beneficial population (figure 2). In particular, emerging interkingdom interactions were discovered in the DCB group; for example, correlations between the differential bacterial species, such as s_Butyricimonas_virosa and s_Candidatus_Avimonas_narfia, and metabolites, such as syringol, sphingosine, and L-glutamine (figure 2). The low number of fungi highly associated with other multi-kingdom taxa in the DCB and NDB groups suggests that fungi play a relatively stable role throughout immunotherapy.

Figure 2

Coabundance correlations among multiomics taxa in patients with DCB and NDB groups. Moderate coabundance networks in DCB and NDB groups, with absolute correlations above 0.4 and with a significance cut-off of p value<0.05. The color of the line represents the size of the correlation coefficient, lighter means a smaller correlation coefficient and darker means a larger correlation coefficient. DCB, durable clinical benefit; NDB, non-durable clinical benefit.

Microbial functional alterations in ICI-treated HCC

Owing to the great inter-individual heterogeneity of the microbiota, it seems possible that different strains in different individuals may use common pathways to trigger similar pathological changes. Therefore, functional annotation by multiomics was used to explore some of the possible variations in micro-organisms in HCC immunotherapy and to provide possible research strategies. Nine KEGG pathways, 74 KEGG Orthology (KO) genes, and 18 metacyclic pathways for the metagenome and 35 KEGG pathways (based on the strength of the correlation, 15 of which were closely related to metabolic pathways) for LC-MS/MS were analyzed (figure 3A–C; online supplemental tables S5,S7,S8). Based on the level of KEGG pathways for metagenomic sequencing, alpha diversity did not differ between the two groups; however, beta diversity was significantly higher in the DCB group than in the NDB group (online supplemental figure S9A,B). Alpha diversity was relatively stable with the use of ICI, while whole beta diversity first showed a declining trend and then a gradual recovery (online supplemental figure S4E,F). Based on the KEGG pathway database and metabolite annotation analysis, KO00250 (alanine, aspartate, and glutamate metabolism), KO00910 (nitrogen metabolism), and KO00750 (vitamin B6 metabolism) were the pathways of interest common to both bacteria and metabolites, with KO00250 being the most important pathway for both bacteria and metabolites (figure 3B,C). Both KO00250 and KO00910 were detected in the NDB group (figure 3D). Bacteria (KO genes) and metabolites interacted through the KO00250 pathway to modulate the response of patients with HCC to immunotherapy, which may be an important way for microbes to regulate immune mechanisms. The metabolite L-glutamine played an important role in three of these pathways (online supplemental figure S9C). Functional annotation and prediction of fungi using FUNGuild revealed that there was no difference between the DCB and NDB groups based on trophic mode, trait, and guild classifications, except for differences in growth classification annotated as secotioid and gasteroid (online supplemental figure S10 and table S9).

Figure 3

Functional annotation for metagenomic and metabolomic analysis. (A) Differential KEGG pathways for metagenomic sequencing. (B) Sankey bubble plot shows differences in the KEGG pathway, KO gene, and KO gene’s contribution to the KEGG pathway. (C) Enrichment metabolic pathways for differential metabolite annotation. (D) Relative abundance comparison of KO00250, KO00910, and KO00750 pathways between the DCB and NDB groups. DCB, durable clinical benefit; KEGG, Kyoto Encyclopedia of Genes and Genomes; KO, KEGG Orthology; LDA, linear discriminant analysis; NDB, non-durable clinical benefit.

Multi-kingdom intercorrelation and correlation with KEGG pathways

To further explore the ecological associations of bacteria, fungi, and metabolites, we analyzed the intracorrelations and intercorrelations among the 25 species, 18 fungi, and 37 metabolites, as well as their associations with nine KEGG pathways (figure 4). Overall, there were significantly more intracorrelations than intercorrelations, and intercorrelations were more frequent for bacterial-metabolite associations, whereas fungal-metabolite and bacterial associations were relatively rare. Regarding the correlation between multi-kingdom micro-organisms and KEGG pathways, s_Faecalibacterium_prausnitzii in bacteria was significantly associated with all nine KEGG pathways. In addition, KO00250 was positively correlated with s_Phascolarctobacteriumsuccinatutens, s_Roseburia_inulinivorans, k_Fungi; g_Papiliotrema, and k_Fungi; f_Rhynchogastremataceae. KO00910 was negatively correlated with s_Catabacter_hongkongensis, s_Massilimaliae_massiliensis, and s_Sellimonas_intestinalis and positively correlated with k_Fungi; g_Trichosporon, and k_Fungi; s_Papiliotrema_aurea. Metabolites and fungi showed a relatively weaker association with KEGG, whereas bacteria showed a stronger association with KEGG and were mostly positively correlated with metabolites expressed or enriched in the DCB group. The overall multi-kingdom microbial correlations suggest that bacteria are more closely associated with metabolites and KEGG pathways, whereas fungi have relatively fewer associations, suggesting that bacteria may play a greater role in HCC immunotherapy and influence KEGG pathways through multiple metabolites.

Figure 4

Multi-kingdom intercorrelation analysis among bacteria, fungi, metabolites, and KEGG pathways. DCB, durable clinical benefit; KEGG, Kyoto Encyclopedia of Genes and Genomes; NDB, non-durable clinical benefit.

Predictive model and the association between gut microbiome, microbe-derived metabolites, and survival

To investigate potential biomarkers for predicting response status, we developed a LASSO regression model based on the species in the DCB and NDB groups. In the training set, 18 key microbial biomarkers were identified and a LASSO classifier was constructed with an area under the curve (AUC) value of 100% (figure 1A; online supplemental figure S7A). It could also predict the response status in the test set, with an AUC value of 75.63% (95% CI: 58.34% to 92.93%; figure 5A). These results suggest that the model of 18 bacterial associations is a good predictor and differentiator of whether patients with HCC will derive clinical benefit from immunotherapy.

Figure 5

Predictive model and the association between gut microbiome, microbe-derived metabolites and survival. (A) Receiver operating characteristic curve of least absolute shrinkage and selection operator regression model for the test set. (B) Progression-free survival (PFS) and overall survival (OS) depended on the relative abundance of s_Actinomyces_sp_ICM47, s_Senegalimassilia_anaerobia, and galanthaminone. PFS and OS of s_Actinomyces_sp_ICM47 combined with galanthaminone (C) and s_Senegalimassilia_anaerobia combined with galanthaminone (D). Good signature: coexistence of depleted s_Actinomyces_sp_ICM47/s_Senegalimassilia_anaerobia and enriched galanthaminone. Poor signature: coexistence of enriched s_Actinomyces_sp_ICM47/s_Senegalimassilia_anaerobia and depleted galanthaminone. Moderate signature: coexistence of depleted both or enriched both of gut microbiome and metabolites.

To further explore the microbial taxa and metabolites associated with survival, we stratified patients with HCC in the training set into high-abundance and low-abundance groups based on the mean abundance of these differentially enriched taxa in the DCB and NDB groups. Patients with enriched s_Actinomyces_sp_ICM47 had significantly worse PFS (median PFS: 7.45 vs 12.0 months, p=0.013) and OS (median OS: 17.3 vs NA months, p=0.021; figure 5B). s_Senegalimassilia_anaerobia was also associated with worse PFS (median PFS: 6.37 vs 12.0 months, p=0.025) and OS (median OS: 14.4 vs NA months, p=0.043; figure 5B). A survival benefit was observed in patients with enriched galanthaminone (median PFS: 15.18 vs 9.33 months; p=0.029; median OS: NA vs 24.2 months, p=0.028; figure 5B). Patients with HCC were further categorized into good, moderate, and poor survival groups based on whether the bacteria or metabolites were simultaneously beneficial for survival. The results of univariate and multivariate analyses suggested that age (<60 years), ECOG PS (score 1), and maximum tumor diameter (≥5 cm) were risk factors, while s_Actinomyces_sp_ICM47 or s_Senegalimassilia_anaerobia combined with galanthaminone (good group) was a positive factor for OS in patients with HCC (online supplemental table S10). Furthermore, combined bacterial and metabolite survival analyses revealed that low expression of s_Actinomyces_sp_ICM47 or s_Senegalimassilia_anaerobia combined with galanthaminone predicted better survival (figure 5C,D). Our findings suggest that the enrichment of specific taxa or metabolites in gut microbes may enhance the clinical response to anti-PD-1 therapy and prolong survival, thereby predicting survival benefits in patients with HCC receiving immunotherapy. The combination of bacteria and metabolites predicted even higher survival, suggesting that the co-presence of these highly enriched bacteria or metabolites may potentiate the effects of immunotherapy in patients via some mechanism. An analysis of the patients, grouped into good, moderate, and poor survival categories based on the combination of Actinomyces_sp_ICM47 and galanthaminone, did not reveal any statistically significant differences in the incidence of overall irAEs or organ-specific irAEs. Similarly, when the patients were grouped based on Senegalimassilia_anaerobia and galanthaminone, comparable results were observed (online supplemental table S11).

Gut bacteria, fungi, and metabolites are influenced by clinical factors

Gut microbiology, including bacteria, fungi, and metabolites, is a complex, dynamically balanced system influenced by multiple factors. To clarify the possible influence of clinical phenotypes on gut microbes, we performed a correlation analysis of differential bacteria, fungi, metabolites, and KEGG pathways with clinical phenotypes and explored the influence of clinical phenotypes on the distribution of bacteria and metabolites using CCA analysis. Association analyses suggested that differentially expressed bacteria, such as s_Clostridium_scindens and s_Clostridia_bacterium_UC5_1_1D1, were strongly associated with clinical phenotypes such as total bilirubin (Tbil) and BA. In contrast, bacteria predicted to be survival associated, such as s_Actinomyces_sp_ICM47 and s_Senegalimassilia_anaerobia, were less influenced by clinical factors (figure 6). Tumor number and intrahepatic metastasis (IHM) were strongly associated with differential metabolites, whereas the survival-related metabolite, galanthaminone, was not affected by clinical factors (figure 6). Most of the differential fungi and KEGG pathways were less influenced by the clinical phenotype, whereas the KO00250 pathway was not significantly associated with the clinical phenotype, revealing the relative stability of the fungi and the importance of the role played by KO00250 pathway throughout immunotherapy (figure 6). The results of the CCA analysis suggested that the overall effect of the clinical phenotype on bacterial distribution was not significant, whereas it was relatively significant for the overall distribution of metabolites between the DCB and NDB groups (figure 7A,B,D,E). Multiple tumor-associated factors, such as aspartate aminotransferase, gamma glutamyl transferase (GGT), alanine aminotransferase, size, and vascular invasion significantly influenced the distribution of gut bacteria, such as s_Candidatus_Avimonas_narfia, whereas GGT, carbohydrate antigen 19–9, IHM, tumor number, and TBil significantly influenced the distribution of metabolites, such as galanthaminone, and sphingosine (figure 7B,C,E,F). The effect of clinical factors on the heterogeneity of gut microbial enrichment in the DCB and NDB groups suggests that the clinical response and survival benefit of ICI treatment depend on the diversity of microbes and differentially enriched taxonomic communities throughout the gut.

Figure 6

Correlations between gut bacteria, fungi, metabolites, and KEGG pathways with clinical factors. Spearman’s rank correlation analysis with two-tailed p values. *p<0.05; **p<0.01; ***p<0.001. DCB, durable clinical benefit; KEGG, Kyoto Encyclopedia of Genes and Genomes; NDB, non-durable clinical benefit.

Figure 7

The distribution of gut microbiome and microbe-derived metabolites was affected by clinical factors. Canonical correspondence analysis (CCA) with permutation test showed the clinical factors associated with the distribution and contribution for gut microbiome (A) and microbe-derived metabolites (B) of patients with HCC in DCB and NDB groups. *p<0.05. ALP, alkaline phosphatase; ALT, alanine aminotransferase; AST, aspartate aminotransferase; BA, bile acids; CA19-9, carbohydrate antigen 19–9; CEA, carcinoembryonic antigen; DCB, durable clinical benefit; ECOG, Eastern Cooperative Oncology Group; GGT, gamma glutamyl transferase; IHM, intrahepatic metastasis; NDB, non-durable clinical benefit; TBil, total bilirubin; TNM, tumor, node, metastasis; VI, vascular invasion.

Discussion

Immunotherapy has provided new hope for treating advanced HCC. The stability of gut micro-organisms plays a regulatory role in many diseases and states in the human body. In this study, the status and characterization of multi-kingdom microbiota, including gut bacteria, fungi, and their metabolites, in patients with HCC treated with ICI were described by multiomics sequencing for the first time. We also explored the alterations in gut microbes in the DCB and NDB groups, established 18 species of bacterial models as predictive biomarkers for predicting whether immunotherapy is of sustained benefit, and screened 2 bacterial species and 1 metabolite as prognostic biomarkers for predicting survival in patients with HCC treated with ICI. These results provide a foundation for future efforts to maximize the benefits of immunotherapy in patients with advanced HCC by determining and modifying gut microbial status.

Several studies, including those on HCC, have confirmed that bacterial diversity is correlated with the efficacy of tumor immunotherapy; however, the results are not entirely consistent.3 5 28 In our study, alpha diversity based on bacteria and fungi at baseline did not differ between the DCB and NDB groups; however, beta diversity was higher in the DCB group than in the NDB group (significant bacteria-based and non-significant fungal-based differences). Bacteria-based differences in beta diversity disappeared after 6–8 weeks, whereas fungi-based dynamic alpha diversity showed an overall declining trend followed by recovery. The time point before the third course of immunotherapy (around 6–8 weeks after initiation of immunotherapy) is a watershed point in terms of the greater impact on bacterial and fungal diversity, and it is possible that this point is the cut-off point at which differences in the ability of patients with HCC to have a sustained benefit from immunotherapy begin to emerge. In future studies, more attention should be paid to changes in diversity after 6–8 weeks of immunotherapy, in addition to baseline gut bacterial and fungal diversity.

Studies have suggested that gut bacteria are good predictive biomarkers for tumor immunotherapy.29 30 Fungi may be relatively more stable across disease states than bacteria, but can exert their influence on diseases through interactions with bacteria.10 24 31–33 A cohort study of four types of tumors (gastrointestinal cancer, melanoma, non-small cell lung cancer, and renal cell carcinoma) suggested that fungal and fungal-bacterial combinations are good predictors of efficacy in tumor immunotherapy.32 In our study, based on the fact that the 18 key microbes were a good class of taxa for predicting the efficacy of HCC immunotherapy, fungal differences were relatively less pronounced than those of bacteria in the DCB and NDB groups. Two bacteria, a metabolite, and a bacteria-metabolite combination, were prognostic predictors of HCC using immunotherapy; however, no suitable markers were screened at the fungal level. Based on the above findings, bacteria are good predictive and prognostic markers for tumors; however, more studies are needed to confirm whether fungi can be used as stable markers.

The bacteria and metabolites differed significantly between the DCB and NDB groups, whereas the differences were smaller for fungi. Faecalibacterium prausnitzii was mainly enriched in the NDB group and had the highest abundance and LDA scores among the species in the DCB and NDB groups. F. prausnitzii has been consistently reported as one of the main butyrate producers in the intestine.34 Butyrate plays a crucial role in the gut physiology and host well-being. Interestingly, higher basal fecal levels of butyric acid have also been associated with decreased survival.35 F. prausnitzii may reduce the efficacy of butyric acid immunotherapy; however, the exact mechanism remains unknown. In our study, both Actinomyces and Senegalimassilia were potential prognostic markers that were negatively associated with immunotherapy. Actinomyces were enriched in non-responders to advanced melanoma treated with immunotherapy, a non-response-associated module that linked Actinomyces to butyric acid (and derivatives) and cyclohexanecarboxylic acid ethyl ester, as well as to higher levels of neutrophils.35 Actinomyces is more abundant in the high-risk group with diffuse large B-cell lymphoma36 and lung cancer.37 The presence of Senegalimassilia, which was detected almost exclusively in the NDB group in the current study (25% in NDB and 4% in DCB), may have a poorer prognosis, and a larger cohort is needed to confirm its feasibility and stability. There are fewer studies on Senegalimassilia and the results are inconsistent across different tumors. For example, Yin et al38 found that Senegalimassilia is a risk factor for bladder cancer, while Jiang et al39 reported that Senegalimassilia is indeed a protective factor for pancreatic cancer. No mechanistic studies of Senegalimassilia in tumorigenesis or tumor immunotherapy have been reported. The cornerstone of ICI treatment for tumors lies in the augmentation of T cell-mediated immune cytotoxicity by targeting coinhibitory molecules such as PD-1/PD-L1, thereby preventing immune evasion by tumor cells. The composition of gut bacteria and their metabolic products exert a significant influence on the efficacy of immunotherapy by modulating various types of T cells, either enhancing or inhibiting the therapeutic response.40 For instance, specific bacteria, including Lactobacillus rhamnosus, have been shown to enhance the antitumor activity of PD-1 immunotherapy through the cyclic GMP-AMP synthase (cGAS)/stimulator of IFN genes (STING) signaling pathway, leading to the activation of interferon-α/β signaling and cytotoxic CD8+T cells.41 42 SCFAs have been implicated in limiting the antitumor effects of cytotoxic T-lymphocyte-associated antigen-4 (CTLA-4) blockade by alleviating the suppressive activity of regulatory T cells. Moreover, a higher concentration of butyrate has been found to decrease the anticancer activity of ipilimumab by inhibiting the accumulation of CD4+T cells involved in the response.43 Additionally, gut microbes play a role in nutrient metabolism and produce metabolites that can influence the body’s immune response, ultimately impacting the effectiveness of ICIs.40

The mechanisms of action of bacteria, fungi, and metabolites in tumor immunotherapy may be multifaceted.44 45 In our study, the KO00250 and KO00910 pathways were mainly upregulated in the NDB group, which may be associated with negative immunotherapeutic effects, and the metabolite L-glutamine may play an important role. Glutamine is an important non-essential amino acid. In tumor cells, “glutamine addiction” is evident. Studies have shown that an increase in glutamine levels promotes purine metabolism to support cancer cell growth and metabolism.46 47 Additionally, altered glutamine metabolism contributes to a hypoxic, acidic, nutrient-depleted tumor microenvironment (TME), which is detrimental to the antitumor immune response.47 In the KO00250 path