Skip to main content

An exosome mRNA-related gene risk model to evaluate the tumor microenvironment and predict prognosis in hepatocellular carcinoma

Abstract

Background

The interplay between exosomes and the tumor microenvironment (TME) remains unclear. We investigated the influence of exosomes on the TME in hepatocellular carcinoma (HCC), focusing on their mRNA expression profile.

Methods

mRNA expression profiles of exosomes were obtained from exoRBase. RNA sequencing data from HCC patients’ tumors were acquired from The Cancer Genome Atlas (TCGA) and the International Cancer Genome Consortium (ICGC). An exosome mRNA-related risk score model of prognostic value was established. The patients in the two databases were divided into high- and low-risk groups based on the median risk score value, and used to validate one another. Functional enrichment analysis was performed based on a differential gene prognosis model (DGPM). CIBERSORT was used to assess the abundance of immune cells in the TME. The correlation between the expression levels of immune checkpoint-related genes and DGPM was analyzed alongside the prediction value to drug sensitivity.

Results

A prognostic exosome mRNA-related 4-gene signature (DYNC1H1, PRKDC, CCDC88A, and ADAMTS5) was constructed and validated. A prognostic nomogram had prognostic ability for HCC. The genes for this model are involved in extracellular matrix, extracellular matrix (ECM)-receptor interaction, and the PI3K-Akt signaling pathway. Expression of genes here had a positive correlation with immune cell infiltration in the TME.

Conclusions

Our study results demonstrate that an exosome mRNA-related risk model can be established in HCC, highlighting the functional significance of the molecules in prognosis and risk stratification.

Peer Review reports

Introduction

Hepatocellular carcinoma (HCC) is one of the most fatal cancers worldwide [1]. In the past 20 years, treatment options of HCC have been primarily chemotherapy, radiotherapy, and surgery. Recently, anti-angiogenic drugs and immune checkpoint inhibitors (ICIs) have demonstrated promise here [2]. However, a large number of HCC patients inevitably experience relapse or disease progression after initial treatment [3]. To identify high-risk patients with a poor prognosis and inform treatment decisions are of vital importance in HCC.

Extracellular vesicles (EVs) are defined as lipid bilayer packages of biological materials, released from various type of cells into surrounding environment. EVs have been considered as ideal biomarkers for the diagnosis and prognosis of cancer [4]. EVs include particles such as exosomes, microvesicles, ectosomes, oncosomes, and apoptotic bodies. Among the main types of EVs, exosomes contain large amounts of RNAs, which can be transmitted among cells and modulate the gene expression of recipient cells [5]. Exosomes contain messenger RNA (mRNA), circular RNA (circRNA), long non-coding RNA (lncRNA), microRNA (miRNA), lipids, and proteins [6]. As intercellular messengers between cells, exosomes can regulate cell differentiation and tissue development [7]. Exosomal RNAs can interact with many types of cancers and are associated with several hallmarks features of cancer. The liquid biopsy approach of exosomes has been used as tumor markers [8]. Based on the ability of exosomes to carry biomolecules to different tissues, exosomes also have application potential in cancer therapy [9]. The potential of using exosomes to predict response to immunotherapy has also been investigated [10].

Exosomes play a vital role in the development, progression, and metastases of HCC [11]. Previous studies confirmed that exosomes could promote progression and metastasis of HCC by regulating multiple signaling pathways and modulating the TME [12]. Furthermore, as exosomes release inhibitory and stimulatory contents that facilitate the cross-talk of tumor cells and the TME, exosomes demonstrate potential for overcoming resistance mechanisms of anti-cancer drugs.

In this study, we aimed to explore potential functional mRNAs in progression and development of HCC and the immune microenvironment of HCC. The study results highlight that exosomal mRNAs which could act as prognostic biomarkers for HCC.

Materials and methods

The mRNA expression data collection

The exoRBase database (http://exorbase.org/exoRBaseV2/download/toIndex) is a depository of mRNAs, lncRNAs, and circRNAs from RNA sequencing (RNA-seq) data analyses in different human body fluids [13]. ExoRBase provides expression landscapes and a comprehensive annotation of extracellular vesicle long RNAs (exLRs), which will help discover novel exLR signatures and facilitate the identification of new circulating biomarkers for cancer therapy. In the current study, the mRNA expression profiles were obtained from the exoRBase database up to April 30, 2023, which included 112 HCC, 130 benign tumor and 118 healthy blood samples.

TCGA-LIHC cohort and ICGC (LIRI-JP) cohort

A total of 374 HCC patients from the TCGA-LIHC cohort were identified. Among them, the level 3 RNA-seq data, somatic mutations data and clinical data of 371 HCC patients with complete information were retrieved from the TCGA website (https://portal.gdc.cancer.gov/projects/TCGA-LIHC/). RNA-seq data, somatic mutations data and clinical information of 231 tumor samples (LIRI-JP cohort) from the ICGC database were downloaded from ICGC portal (https://dcc.icgc.org/projects/LIRI-JP). The samples from LIRI-JP cohort were primarily derived from Japanese HCC patients with HBV or HCV infection.

Differentially expressed mutant genes (DEMGs)

The 118 healthy blood samples were used as the normal group. The differentially expressed genes (DEGs) in blood samples from 130 benign and 112 HCC patients were identified of p < 0.05 and fold-change (FC) > 0 via “limma” R package, respectively. Integration analyses between normal vs. HCC and normal vs. benign comparisons were performed. Those mRNAs that were only significantly differentially expressed in normal vs. HCC comparison (not in normal vs. benign comparison) were screened out. Then, the DEMGs of HCC were obtained from the intersection genes that have mutations in TCGA, ICGC, and exoRBase.

Construction and validation of a differential gene prognosis model (DGPM)

The DEGs between tumor and adjacent tissues were identified with a false discovery rate (FDR) < 0.05 in the TCGA cohort. Univariate Cox analysis was carried out to screen DEMGs with prognostic values in terms of overall survival (OS). P values were adjusted by Benjamini & Hochberg (BH) correction. Visualization and comparison of gene alterations were conducted with cBioPortal. LASSO-penalized Cox regression analysis was used to construct a DGPM to minimize the risk of overfitting. LASSO algorithm was used for to select variables with “glmnet” R package. Penalty parameter (λ) was chosen by 10-fold cross-validation with the minimum criteria. The risk scores were calculated with the normalized expression level of each gene and corresponding regression coefficients. The risk score was calculated with the following formula:

$$The\;risk\;score=\sum_{i=1}^n{Coef}_i\times{Expr}_i$$

where Expri indicates the expression level of gene i, and coefi means the regression coefficient of gene i.

The patients were separated into high- and low-risk groups due to the median value of the risk score. Principal component analysis (PCA) was performed with “stats” R package. Besides, t-distributed stochastic neighbor embedding (t-SNE) was used to investigate the distribution of different groups using “Rtsne” R package. The optimal cut-off value for gene expression was determined by “survminer” R package. The time-dependent receiver operating characteristic (ROC) curve analyses were conducted with “survivalROC” R package. Survival analysis was performed to validate the prognostic performance independent from clinical parameters. For the validation studies, LIRI-JP cohort was employed. The risk score was calculated with the same formula used with TCGA cohort. The patients in the ICGC cohort were also divided into low- or high-risk subgroups by applying the median risk score from TCGA cohort, and these groups were then compared to confirm the gene model.

Functional enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed with “clusterProfiler” R package based on the DEMGs between high- and low-risk groups. The infiltrating score of immune cells and the activity of immune-related pathways were analysed with single-sample gene set enrichment analysis (ssGSEA) with “gsva” R package.

Construction of the prognostic nomogram

To comprehensively assess prognosis predictive ability of risk signature, tumor stage, gender, age, WHO grade, T category, N category and M category for 1-, 3-, and 5-year OS, time-dependent ROC curves was performed to calculate the area under the curve (AUC) values. Prognostic nomogram which containing DEMG-based risk model and clinical parameters was established.

Correlation of risk score with immune cells

The information for immune infiltration was downloaded from the tumor immune estimation resource (TIMER) (https://cistrome.shinyapps.io/timer/). A correlation between prognostic risk score and immune cell infiltration was performed. SsGSEA was performed to investigate the enrichment of the two subgroups in immune function-associated gene sets via “GSEAbase” R package. “ESTIMATE” R package was used to evaluate tumor purity and infiltrating cells, including immune cell and stromal cell. The fraction of 22 immune cell types was assessed by cell type identification by estimating relative subsets of RNA transcripts (CIBERSORT; https://cibersort.stanford.edu/).

Role of risk score in immune checkpoint blockade (ICB) treatment

Herein, 6 key genes of ICB treatment in HCC were extracted, including cytotoxic T-lymphocyte antigen 4 (CTLA-4), programmed death 1 (PD-1, or PDCD1), programmed death ligand 1 (PD-L1, or CD274), programmed death ligand 2 (PD-L2, or PDCD1LG2), T-cell immunoglobulin domain and mucin domain-containing molecule-3 (TIM-3, or HAVCR2), and indoleamine 2,3-dioxygenase 1 (IDO1). DEMG-based prognostic signature and expression level of 6 ICB key genes were correlated. Furthermore, the expression level of 47 ICB-related genes between low- and high-risk groups were also compared.

Statistical analysis

All statistical analyses were conducted by R software (version 4.1.1). Gene expression was analysed by student’s t-test. Differences in proportions were analysed with Chi-squared test. ssGSEA of immune cells or pathways between the groups were examined with Mann-Whitney test. Kaplan-Meier curves were employed to assess survival data. Pearson correlation test was used to analyze the correlation of risk score, clinical parameters, immune cell infiltration, and immune checkpoints. Independent prognostic performance of risk signature was evaluated with Cox regression models. p value < 0.05 was considered as statistically significant.

Results

Identification of prognostic related DEMGs from the exoRBase database

To explore potential mRNAs associated with development and progression of HCC, the blood samples of 112 HCC, 130 benign and 118 healthy people from the exoRBase database were investigated, after differential expression analysis between normal and benign or blood samples of HCC patients. Compared with normal blood samples, 59 and 42 mRNAs were significantly upregulated and downregulated in benign blood samples, respectively (Fig. 1A and Supplementary Table 1). In addition, a total of 132 significant DEGs (106 upregulated and 26 downregulated ones) were identified in normal vs. HCC (Fig. 1B and Supplementary Table 2). Integration analyses between normal vs. HCC and normal vs. benign were performed. The mRNAs only significantly differentially expressed in normal vs. HCC (not in normal vs. benign) were screened out. Consequently, 102 upregulated and 25 downregulated mRNAs were identified (Fig. 1C-D).

Fig. 1
figure 1

Candidate mRNAs in HCC. A Volcano plot of DEGs between 130 benign and 118 healthy blood samples of exosome. B Volcano plot of DEGs between 112 HCC and 118 healthy blood samples of exosome. C Interaction analysis of downregulated (C) and upregulated (D) DEGs in both compared groups

In addition, 374 HCC patients from TCGA-LIHC cohort and 273 HCC patients from ICGC cohort were enrolled. The frequently mutated genes in American HCC samples from TCGA cohort (Fig. 2A) and ICGC cohort (Fig. 2B) were identified. Of note, there were some frequently mutated genes in both American and Japanese patients. Comparative analysis of mutated genes between TCGA and ICGC cohorts were performed (Fig. 2C). Then, we analyzed the DEGs from the exoRBase database, and 24 HCC DEMGs were obtained (Fig. 2D).

Fig. 2
figure 2

Frequently mutated genes in HCC. A The frequently mutated genes in HCC from TCGA cohort were depicted with Oncoplot. B The frequently mutated genes in HCC from ICGC cohort were displayed by waterfall plot. C Venn diagram of genes covered by both TCGA and ICGC cohorts. D Interaction analysis for DEMGs in both compared groups

The expression levels of 24 DEMGs were compared in the pooled TCGA and Genotype-Tissue Expression (GTEx) data from 50 normal and 374 tumor tissues, and 13 DEMGs were identified (Fig. 3A). Among them, 4 genes (APOB, FGB, ALDOB, and ALB) were downregulated while 9 other genes (PPARG, PABPC3, MAPKAPK2, ADAMTS5, CCDC88A, PRRC2C, HNRNPK, DYNC1H1, and PRKDC) were enriched in the tumour group. All of the 13 genes were associated with OS in the univariate Cox regression analysis (Fig. 3B). Mutations of the 13 genes were analyzed by cBioPortal (Fig. 3C). The correlation of the genes is presented in Fig. 3D.

Fig. 3
figure 3

Candidate prognostic related DEMGs in the TCGA cohort. A Heatmap of the prognostic related DEMGs between normal and tumour tissues. B Forest plots showing the results of the univariate Cox regression analysis between gene expression and OS. C Landscape of prognostic related DEMGs alteration in HCC. D The correlation network of candidate genes

Development of a prognostic gene model in TCGA cohort

A prognostic model was established by LASSO Cox regression with the expression profile of the 13 genes identified above. As a result, the 4-gene signature based on the optimal value of λ was identified (DYNC1H1, PRKDC, CCDC88A, and ADAMTS5, Table 1). Survival analyses suggested that high expression of these genes correlated with a poor prognosis according to the optimal cut-off expression value of each gene (all adjusted p < 0.05). The patients were stratified into high- and low-risk group based on the median cut-off value (Fig. 4A). The patients in different risk groups were distributed in two directions as PCA and t-SNE analysis indicated (Fig. 4B-C). Patients in high-risk group had a higher probability of death earlier (Fig. 4D) and a significantly worse OS (Fig. 4E). The AUCs reached 0.733, 0.634, and 0.652 at 1-, 3-, and 5-year, respectively (Fig. 4F).

Table 1 Genes selected to construct the prognostic model
Fig. 4
figure 4

Prognostic value of the risk score in the TCGA cohort. A Distribution and median value of the risk scores. B PCA plot and C t-SNE analysis. D Distributions of OS status, OS and risk score. E Kaplan-Meier curves for the OS of patients in the high- and low-risk group. F AUC of time-dependent ROC curves confirmed the prognostic performance of the risk score. G Comparison of TMB between low- and high-risk groups. H Survival analysis based on the TMB. I Survival analysis for 4 groups by combining TMB and DEMG-based risk signature

Although that there was no significant difference in TMB between patients with high- and low-DEMG (Fig. 4G), low TMB was associated with better OS (Log-rank test, p < 0.001, Fig. 4H). When DEMG and TMB were e integrated to stratify the samples into TMBhigh/DEMGlow, TMBlow/DEMGlow, TMBhigh/DEMGhigh, and TMBlow/DEMGhigh groups, significant differences were found among all groups (Fig. 4I, Log-rank test, p < 0.001), and patients in the TMBlow/DEMGlow group exhibited the best OS.

Validation of the risk signature in ICGC cohort

To validate the robustness of the model constructed from TCGA cohort, a total of 273 HCC patients from ICGC cohort were utilized as the validation cohort. Based on the median risk score in TCGA cohort, the patients in ICGC cohort were also categorized into high- or low-risk groups (Fig. 5A). PCA and t-SNE analysis confirmed the discrete distribution of the patients in two subgroups (Fig. 5B-C). The survival outcomes of high-risk group were similar with those from TCGA cohort, Fig. 5D-E). ROC curve showed that this model had good predictive efficacy (AUC = 0.630, 0.607, and 0.638 for 1-, 3-, and 5-year survival, Fig. 5F). Although TMB is not significant between patients with DEMGhigh and DEMGlow subgroups (Fig. 5G), TMBlow (Fig. 5H) and TMBlow/DEMGlow (Fig. 5I) were also associated with good OS. These results collectively demonstrate that increased risk score was correlated with tumor progression.

Fig. 5
figure 5

Validation of the risk score in the ICGC cohort. A Distribution and median value of the risk scores. B PCA plot and C t-SNE analysis. D Distributions of OS status, OS and risk score. E Kaplan-Meier curves for the OS of patients in the high-risk group and low-risk group. F AUC of time-dependent ROC curves confirmed the prognostic performance of the risk score. G Comparison of TMB between low- and high-risk groups. H Survival analysis based on the TMB. I Survival analysis for 4 groups by combining TMB and DEMG-based risk signature

Independent prognostic value of the risk model

The univariate Cox regression analysis suggested that the risk score was an independent prognostic factor (HR = 2.613, 95% CI: 1.784−3.826 for TCGA and HR: 3.204, 95% CI: 1.730−5.932 for ICGC, Fig. 6A-B). The multivariate Cox regression analysis also indicated that the risk score was a prognostic factor, after adjusting for other confounding factors (HR = 2.405, 95% CI: 1.639−3.528 for TCGA and HR: 2.879, 95% CI: 1.499−5.529 for ICGC, Fig. 6C-D). The clinical features of TCGA cohort suggested that grade and survival status of the patients were diversely distributed between low- and high-risk subgroups (Fig. 6E, p< 0.01).

Fig. 6
figure 6

Univariate and multivariable Cox regression analyses of risk score. A Univariate and B multiple Cox regression analyses were performed in the TCGA cohort. C Univariate and D multiple Cox regression analyses were performed in the ICGC validation cohort. E Heatmap of clinical parameters for the TCGA cohort

Functional analyses

The DEGs were extracted by LASSO regression with “limma” R package using criteria FDR < 0.05 and |log2FC | ≥1 as indicated above. GO enrichment and KEGG pathway analysis were performed with these DEGs. DEGs in TCGA cohort were mainly correlated with the extracellular matrix, extracellular matrix (ECM)-receptor interaction, and PI3K-Akt signaling pathway (Fig. 7A-B). The biological processes, molecular functions and signaling pathways were validated by ICGC cohort (Fig. 7C-D).

Fig. 7
figure 7

Functional analysis of the DGPM in TCGA and ICGC cohort. A Bubble graph for GO enrichment and B barplot graph for KEGG pathways in the TCGA cohort. C Bubble graph for GO enrichment and D barplot graph for KEGG pathways in the ICGC cohort

Comparison of the immune activity between subgroups

The ratio of immune cell infiltration and correlation of immune cells in TCGA and ICGC databases were shown in Fig. 8. The enrichment scores of 16 types of immune cells and the activity of 13 immune-related pathways were compared between low- and high-risk groups in both TCGA and ICGC cohorts by employing ssGSEA. In TCGA cohort (Fig. 9A), higher levels of infiltration of immune cells in high-risk subgroup, especially activated dendritic cells (aDCs), dendritic cells (DCs), immature dendritic cells (iDCs), macrophages, T helper (Th) cells (Tfh, Th1, and Th2 cells) and regulatory T (Treg) cells. Except for the type-1 and type-2 IFN response pathway, other immune pathways exhibited higher activity in the high-risk group in the TCGA cohort (Fig. 9B). When the immune status in ICGC cohort was evaluated, similar conclusions were drawn (Fig. 9C-D).

Fig. 8
figure 8

Immune cell infiltrations of TCGA and ICGC cohorts. Relative proportion of immune cell infiltration in (A) TCGA and (B) ICGC (B). Correlation analysis of immune cells in (C) TCGA and (D) ICGC

Fig. 9.
figure 9

SsGSEA scores between different risk groups in TCGA and ICGC cohort. A The scores of 16 immune cells and B 13 immune-related functions in TCGA cohort are displayed in boxplots. C The scores of 16 immune cells and D 13 immune-related functions in ICGC cohort are displayed in boxplots. Adjusted P values were showed as: *, p < 0.05; **, p < 0.01; ***, p < 0.001

Construction of nomogram with prognostic signature with clinical features

The risk score increased significantly with tumor grade, stage, and T category (Fig. 10A-C). A nomogram was constructed based on the prognostic signature and clinical parameters (Fig. 10D). The nomogram-predicted survival closely matched with the optimal predictive performance. The AUCs for the 1-, 3-, and 5-year were 0.752, 0.687, and 0.710, respectively (Fig. 10E). The nomogram had similar performance to that of an ideal model (Fig. 10F).

Fig 10
figure 10

Correlation of DGPM with clinical features and construction of clinicopathological nomogram. A Correlation of risk score with A tumor grade, B clinical stage, and C T status. D Nomogram was constructed by grade, stage and risk signature for predicting survival. E AUCs for predicting 1-, 3-, and 5-year survival. F The 1-, 3-, and 5-year nomogram calibration curves

PRKDC independently affected prognosis

PRKDC was only gene whose expression level was upregulated among the four prognostic-related DEMGs (DYNC1H1, PRKDC, CCDC88A, and ADAMTS5) (Log FC >1). PRKDC expression level was lower in adjacent normal specimens compared that of tumor tissues (Fig. 11A). There are 4 analyses showing high expression of PRKDC from ONCOMINE website (Fig. 11B). The protein expression level of included gene of signature was verified by The Human Protein Atlas (Fig. 11C-D). The results found that PRKDC expressed statistical significantly among different pathological grades (Fig. 11E, most p < 0.05). The TIMER shows that the clinical outcome of age, gender, race, stage, and purity increased risk with the increase of PRKDC gene expression (Fig. 11F). Kaplan-Meier analysis showed lower PRKDC expression level was correlated with longer OS time (Fig. 11G).

Fig. 11
figure 11

The clinical significance of PRKDC in HCC. A PRKDC are overexpressed in HCC tumor tissue. B Four analyses from ONCOMINE platform showing high expression of PRKDC. C-D Protein expression level of PRKDC was shown by The Human Protein Altas by immunohistochemistry. E Correlation of risk score with tumor grade. F The TIMER shows that the clinical outcome increased risk with the increase of PRKDC gene expression. G Lower PRKDC level predicts longer OS

PRKDC correlated with ICB key genes

The correlation between six key ICB key targets and prognostic signature was analyzed to reveal the potential player of risk signature (Fig. 12A). Besides, expression levels of 36 out of 47 ICB-associated genes between low- and high-PRKDC groups were dysregulated in different subgroups (Fig. 12B). TIMER results adjusted by tumor purity showed PRKDC was positively correlated to CD274 (r = 0.433; p = 3.47e−17), CTLA4 (r = 0.214; p = 5.98e−05), HAVCR2 (r = 0.407; p = 3.43e−15), IDO1 (r = 0.181; p = 7.41e−04), PDCD1LG2 (r = 0.264; p = 6.34e−07) and PDCD1 (r = 0.209; p = 9.32e−05), suggesting PRKDC may exert a vital role in ICB treatment of HCC (Fig. 12C-H).

Fig. 12
figure 12

Association between PRKDC and immune checkpoint genes. A Correlation analysis between immune checkpoints CD274, PDCD1, PDCD1LG2, CTLA4, HAVCR2, and IDO1 and risk score. B Comparison of the expression levels of ICB-related genes between low- and high-PRKDC groups. Association between PRKDC and C CD274, D CTLA4, E HAVCR2, F IDO1, G PDCD1LG2, and H PDCD1

Role of PRKDC in TME

HCC patients were classified into high- and low- PRKDC expression groups based on the median PRKDC expression level. ESTIMATE results suggested that patients with higher PRKDC expression had a higher stromal score, higher ESTIMATE score and lower tumor purity (Fig. 13A). TIMER showed that low PRKDC was associated with good OS (Fig. 13B). Arm-level deletion was the main type of mutation (Fig. 13C-D). PRKDC expression was positively correlated with immune cell infiltration (Fig. 13E). The results of ssGSEA suggested that except for cytolytic activity and NK cells, the infiltration fraction of aDCs, antigen-presentation cell (APC) co-inhibition, iDCs, macrophages, MHC-class-I, para-inflammation and Treg expression were significantly decreased when risk score is declining (Fig. 13F).

Fig. 13
figure 13

The role of PRKDC in TME features. A Comparison of (A1) immune score, (A2) ESTIMATE score and (A3) tumor purity between low- and high-PRKDC groups. B TIMER showed that low PRKDC was associated with good OS. Copy number of C CD4+ T-cells and D CD8+ T-cells in HCC. E Relationship between PRKDC with (E1) CD4+ T-cells, (E2) neutrophils, (E3) macrophages and (E4) myeloid dendritic cells. F Comparison of ssGSEA enrichment between low- and high-PRKDC groups.

Discussion

Exosome-derived mRNAs have gained attention due to their potential role in cancer development, progression, and as a source of diagnostic and therapeutic targets. Exosomes released by cancer cells can carry various molecular cargoes, including mRNA, contributing to the communication between cancer cells and their microenvironment. ExLR, mainly mRNAs are potential biomarkers in HCC [14]. Compared to tissue mRNAs, circulating emRNAs may reflect their original tissues and the relative fraction of immune cell types [15].

Exosome mRNA-related gene expression is associated with cancer development and progression. However, the characteristics and landscape of exosomal mRNAs are not fully understood. In this study, four exosomal mRNA-related genes were included to construct the prognostic model by univariate Cox and LASSO Cox regression analysis. Patients in high- and low-risk group had significantly different survival outcomes. Furthermore, ROC curve demonstrated that the prognostic signature was robust in predicting the OS for HCC patients. Additionally, when combined the prognostic signature with clinical parameters, the nomogram had satisfactory predictive performance for HCC. The prognostic signature correlated with tumor microenvironment and the expression of immune checkpoints. The flow diagram of our study is shown in supplementary Figure S1.

A better understanding of the biology represented by the selected genes can be obtained via analysis of functional networks or pathways that these genes based on their biological functions. The functional analysis indicated the prognostic signature was enriched in extracellular matrix, ECM-receptor interaction, and PI3K-Akt signaling pathway. Tumor cells release exosomes to interact with ECM. In turn, ECM regulates exosome secretion and uptake. The biomolecules of exosomes can impact ECM remodeling, which is associated with cancer progression. Tumor-derived exosomes are capable of modulating ECM and TME by disruption of cell adhesion formation and stimulation of the extracellular receptor signaling. Exosomes reach distal sites where they bind to cell surfaces and experience endocytosis with specific mechanisms [16]. ECM reorganization by exosomes contributes to physiological and pathologic angiogenesis. The surface to volume ratio of EVs is relatively large, which enables efficient surface interactions of EVs with cells and molecules. The surface interactions determine the fate of EVs by orchestra them to certain tissues or cell membrane. Focal adhesion pathway has also been implicated in KEGG pathway, which could regulate ECM. Being dynamic integrin-based adhesion complexes, focal adhesions anchor actin cytoskeleton to ECM, which transfers environmental stimuli to the cells and to change cell motility, adhesion, and shape [17]. Taken together, ECM modulation of the host tissue by tumor-derived exosome is involved in the crosstalk between cancer and the premetastatic niche [18].

Tumor-derived exosomes carry immunostimulatory and immunosuppressive receptor or ligands, partially mimicking the profiles of the parent cells [19]. Exosomes are involved in immune responses for tumorigenesis [20]. Exosomes-mediated signaling is contextual, and tumor-derived exosomes mainly mediate suppression in TME. The role of exosomes in antigen specific immune responses have also been demonstrated. In our study, the enriched immune cell types in high-risk group were aDCs, DCs, iDCs, Th2 cells, and Tregs. The results also found that most of the 13 immune-related functions were highly activated in the high-risk group, especially APC co-inhibition, APC co-stimulation, check-point, HLA, T cell co-stimulation. A large body of evidences supports the potential of tumor-derived exosomes to promote antigen-processing and differentiation capabilities of DC in TME. Exosomes carrying tumor-associated antigens (TAAs) and costimulatory molecules reprogram APCs, which could promote immune responses. The effects of tumor-derived exosomes on T-cell subsets are complex [19]. The proteins carried by exosomes can inhibit cytotoxicity and regulate immune-related genes in T cells [19]. Previous studies have reported tumor-derived exosomes could promote Treg activity and expansion [21, 22]. These results collectively indicated the risk signature was associated with TME.

Inhibitory receptors have been identified in cancers, including but not limited to PD-1, CTLA-4, LAG3, and TIM3, etc. [23],. Many biochemical studies have revealed complex and delicate regulation of checkpoint expression on cell surface [24]. Upon ligand engagement, checkpoints follow distinct signaling mechanisms to inhibit antitumor immunity. Exosomes communicate between tumor cells of immune cells and stromal cells by transferring message, contributing to immune escape. For example, tumor-derived exosomes containing PD-L1 can mimic the function of PD-L1 on cell surface. The association between levels of circulating exosomal PD-L1 and response to anti-PD-1/PD-L1 antibody therapy has been documented [25]. The correlation of the prognostic signature of exosome mRNA-related genes with ICB key targets reveal the risk signature has a potential role in the ICB treatment of HCC.

There are several proposed exosome-mediated drug resistance mechanisms, including exosome-mediated transfer of miRNAs, neutralization of antibody-based drugs, and drug export via the exosome pathway [26]. Exosomes could also induce therapy-resistance by promoting anti-apoptotic pathways and alteration of signaling transduction [27]. The PI3K-Akt pathway has been linked to modulate the multidrug resistance of various cancers [28]. For example, the exosome-mediated PI3k/Akt/mTOR signaling pathway has been implicated in cervical cancer [29]. Similarly, gastric cancer-derived exosomes facilitate the proliferation of recipient cell via PI3K/Akt signaling pathway [30]. However, the role of exosome mRNA-related gene signature in the drug sensitivity prediction needs further validation in larger clinical samples. More studies are required to fully uncover the complexities of exosome-derived mRNA in cancer. Ongoing research aims to unravel the specific mRNA cargoes carried by exosomes, their functional implications, and their potential as therapeutic targets or tools for cancer management.

Among the four genes selected to construct the prognostic signature, PRKDC was recently identified as a new biomarker and potential target for immunotherapy [31]. PRKDC gene encodes catalytic subunit of a nuclear DNA-dependent serine/threonine protein kinase (DNA-PKcs) [32]. As a catalytic protein, DNA-PKcs together with Ku70 and Ku80 constitute a DNA-dependent protein kinase (DNA-PK) [33]. The degradation of DNA-PK is regulated by vasolin-containing protein (VCP), which is found in the exosomes from gliomas [34]. DNA-PK is involved in DNA double-strand break repair, immunocompetence, and genomic integrity [35]. Previous studies indicated DNA-PK is a candidate driver of hepatocarcinogenesis which can predict treatment response and patient survival [36]. DNA-PK has emerged as a potential therapeutic target in various cancer types. By inhibition of its kinase function, DNA-PK inhibitors could potentiate DNA damage [37]. The investigation of DNA-PK inhibitor employs with monotherapy and combination strategy [38].

Our study is the first to construct and validate an exosome mRNA-related gene risk model based on four exosome-related genes, which can serve as an independent prognostic factor in HCC patients. However, our study has some limitations. First, the validation of the prognostic signature was performed in public database, larger clinical cohorts are need to confirm the value of this signature. Second, the risk model was solely established on exosome mRNA-related genes, the information of lncRNA, miRNA, and circRNA were not included in our study. Third, the relationship of this signature with tumor microenvironment is preliminary and hypothesis-generating, therefore, experimental confirmation is required in the future.

Conclusions

In summary, an exosome mRNA-related prognostic risk model was established and validated to predict the prognosis of HCC patients and associated with immune infiltration. The risk model can serve as an independent prognostic factor for HCC patients and highlights the functional significance of mRNAs in exosomes.

Availability of data and materials

Publicly available datasets from TCGA (https://portal.gdc.cancer.gov/), ICGC (https://dcc.icgc.org/) and exoRBase (http://exorbase.org/exoRBaseV2/download/toIndex) were analyzed in this study. were analyzed in this study.

References

  1. Llovet JM, Kelley RK, Villanueva A, et al. Hepatocellular carcinoma. Nat Rev Dis Primers. 2021;7(1):6.

    Article  PubMed  Google Scholar 

  2. Rinaldi L, Vetrano E, Rinaldi B, et al. HCC and molecular targeting therapies: back to the future. Biomedicines. 2021;9(10):1345.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Craig AJ, von Felden J, Garcia-Lezana T, et al. Tumour evolution in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol. 2020;17(3):139–52.

    Article  PubMed  Google Scholar 

  4. Guo L, Guo N. Exosomes: Potent regulators of tumor malignancy and potential bio-tools in clinical application. Crit Rev Oncol Hematol. 2015;95(3):346–58.

    Article  PubMed  Google Scholar 

  5. Zhang L, Yu D. Exosomes in cancer development, metastasis, and immunity. Biochim Biophys Acta Rev Cancer. 2019;1871(2):455–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. van Niel G, D’Angelo G, Raposo G. Shedding light on the cell biology of extracellular vesicles. Nat Rev Mol Cell Biol. 2018;19(4):213–28.

    Article  PubMed  Google Scholar 

  7. Zhang Y, Liu Y, Liu H, et al. Exosomes: biogenesis, biologic function and clinical potential. Cell Biosci. 2019;9:19.

    Article  PubMed  PubMed Central  Google Scholar 

  8. Yu W, Hurley J, Roberts D, et al. Exosome-based liquid biopsies in cancer: opportunities and challenges. Ann Oncol. 2021;32(4):466–77.

    Article  CAS  PubMed  Google Scholar 

  9. Mathew M, Zade M, Mezghani N, et al. Extracellular vesicles as biomarkers in cancer immunotherapy. Cancers (Basel). 2020;12(10):2825.

    Article  CAS  PubMed  Google Scholar 

  10. Gilligan KE, Dwyer RM. Engineering exosomes for cancer therapy. Int J Mol Sci. 2017;18(6):1122.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Li LM, Liu ZX, Cheng QY. Exosome plays an important role in the development of hepatocellular carcinoma. Pathol Res Pract. 2019;215(8):152468.

    Article  CAS  PubMed  Google Scholar 

  12. Wang H, Lu Z, Zhao X. Tumorigenesis, diagnosis, and therapeutic potential of exosomes in liver cancer. J Hematol Oncol. 2019;12(1):133.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Lai H, Li Y, Zhang H, et al. exoRBase 2.0: an atlas of mRNA, lncRNA and circRNA in extracellular vesicles from human biofluids. Nucleic Acids Res. 2022;50(D1):D118-D28.

  14. Li Y, He X, Li Q, et al. EV-origin: enumerating the tissue-cellular origin of circulating extracellular vesicles using exLR profile. Comput Struct Biotechnol J. 2020;18:2851–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Li Y, Zhao J, Yu S, et al. Extracellular vesicles long rna sequencing reveals abundant mrna, circRNA, and lncRNA in human blood as potential biomarkers for cancer diagnosis. Clin Chem. 2019;65(6):798–808.

    Article  CAS  PubMed  Google Scholar 

  16. Gurung S, Perocheau D, Touramanidou L, et al. The exosome journey: from biogenesis to uptake and intracellular signalling. Cell Commun Signal. 2021;19(1):47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Wozniak MA, Modzelewska K, Kwong L, et al. Focal adhesion regulation of cell behavior. Biochim Biophys Acta. 2004;1692(2–3):103–19.

    Article  CAS  PubMed  Google Scholar 

  18. Mu W, Rana S, Zoller M. Host matrix modulation by tumor exosomes promotes motility and invasiveness. Neoplasia. 2013;15(8):875–87.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Whiteside TL. The effect of tumor-derived exosomes on immune regulation and cancer immunotherapy. Future Oncol. 2017;13(28):2583–92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Zhang X, Yuan X, Shi H, et al. Exosomes in cancer: small particle, big player. J Hematol Oncol. 2015;8:83.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Wieckowski EU, Visus C, Szajnik M, et al. Tumor-derived microvesicles promote regulatory T cell expansion and induce apoptosis in tumor-reactive activated CD8+ T lymphocytes. J Immunol. 2009;183(6):3720–30.

    Article  CAS  PubMed  Google Scholar 

  22. Szajnik M, Czystowska M, Szczepanski MJ, et al. Tumor-derived microvesicles induce, expand and up-regulate biological activities of human regulatory T cells (Treg). PLoS One. 2010;5(7):e11469.

    Article  PubMed  PubMed Central  Google Scholar 

  23. He X, Xu C. Immune checkpoint signaling and cancer immunotherapy. Cell Res. 2020;30(8):660–9.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Jung CY, Antonia SJ. Tumor immunology and immune checkpoint inhibitors in non-small cell lung cancer. Tuberc Respir Dis (Seoul). 2018;81(1):29–41.

    Article  PubMed  Google Scholar 

  25. Yin Z, Yu M, Ma T, et al. Mechanisms underlying low-clinical responses to PD-1/PD-L1 blocking antibodies in immunotherapy of cancer: a key role of exosomal PD-L1. J Immunother Cancer. 2021;9(1):e001698.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Giallombardo M, Taverna S, Alessandro R, et al. Exosome-mediated drug resistance in cancer: the near future is here. Ther Adv Med Oncol. 2016;8(5):320–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Mashouri L, Yousefi H, Aref AR, et al. Exosomes: composition, biogenesis, and mechanisms in cancer metastasis and drug resistance. Mol Cancer. 2019;18(1):75.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Liu R, Chen Y, Liu G, et al. PI3K/AKT pathway as a key link modulates the multidrug resistance of cancers. Cell Death Dis. 2020;11(9):797.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Zhang W, Zhou Q, Wei Y, et al. The exosome-mediated PI3k/Akt/mTOR signaling pathway in cervical cancer. Int J Clin Exp Pathol. 2019;12(7):2474–84.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. Zhang S, Zhang Y, Qu J, et al. Exosomes promote cetuximab resistance via the PTEN/Akt pathway in colon cancer cells. Braz J Med Biol Res. 2017;51(1):e6472.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Tan KT, Yeh CN, Chang YC, et al. PRKDC: new biomarker and drug target for checkpoint blockade immunotherapy. J Immunother Cancer. 2020;8(1):e000485.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Lee HS, Yang HK, Kim WH, et al. Loss of DNA-dependent protein kinase catalytic subunit (DNA-PKcs) expression in gastric cancers. Cancer Res Treat. 2005;37(2):98–102.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Goodwin JF, Knudsen KE. Beyond DNA repair: DNA-PK function in cancer. Cancer Discov. 2014;4(10):1126–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Jiang N, Shen Y, Fei X, et al. Valosin-containing protein regulates the proteasome-mediated degradation of DNA-PKcs in glioma cells. Cell Death Dis. 2013;4:e647.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Lieber MR. The mechanism of double-strand DNA break repair by the nonhomologous DNA end-joining pathway. Annu Rev Biochem. 2010;79:181–211.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Cornell L, Munck JM, Alsinet C, et al. DNA-PK-A candidate driver of hepatocarcinogenesis and tissue biomarker that predicts response to treatment and survival. Clin Cancer Res. 2015;21(4):925–33.

    Article  CAS  PubMed  Google Scholar 

  37. Mohiuddin IS, Kang MH. DNA-PK as an emerging therapeutic target in cancer. Front Oncol. 2019;9:635.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Hu S, Hui Z, Lirussi F, et al. Small molecule DNA-PK inhibitors as potential cancer therapy: a patent review (2010-present). Expert Opin Ther Pat. 2021;31(5):435–52.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

None.

Funding

This study is supported by grants from Medical Scientific Research Foundation of Zhejiang Province, China (2022KY545).

Author information

Authors and Affiliations

Authors

Contributions

ZW, LP and LZ designed and supervised the study. ZD, XH, LP, ZW, and LL analyzed the data and wrote the original draft. LP, XH, LC and JS edited the draft. All the authors have read and approved the final manuscript.

Corresponding authors

Correspondence to Ling Peng or Zhiqiang Wang.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

JS’ conflicts can be found at https://www.nature.com/onc/editors. None are relevant here. Other authors: none declared.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Du, Z., Han, X., Zhu, L. et al. An exosome mRNA-related gene risk model to evaluate the tumor microenvironment and predict prognosis in hepatocellular carcinoma. BMC Med Genomics 17, 86 (2024). https://doi.org/10.1186/s12920-024-01865-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-024-01865-z

Keywords