Differentially expressed liver exosome-related genes as candidate prognostic biomarkers for hepatocellular carcinoma.

Background: Exosomes are involved in cell-to-cell communication, neovascularization, cancer metastasis, and drug resistance, which all play an important role in the occurrence and progression of hepatocellular carcinoma (HCC). Because there are few mechanistic studies about the function of exosomes in HCC, the goals of this study were to identify exosome-related genes in HCC, to establish a reliable prognostic model for HCC, and to explore underlying mechanisms.
Methods: The exoRBase and The Cancer Genome Atlas (TCGA) databases were used to analyze differentially expressed genes (DEGs). Cox regression and least absolute shrinkage and selection operator analyses were used to identify DEGs closely related to the overall survival of patients with HCC. An exosome-related prognostic model was then constructed in TCGA and validated in the International Cancer Genome Consortium database. A nomogram was developed to predict survival. CIBERSORT was used to estimate the abundance of different types of immune cells. Immunotherapy-related DEGs were used to predict the effect of immunotherapy.
Results: Forty-eight exosome-related DEGs were obtained; of them, five [exportin 1 (XPO1), lysosomal thiol reductase (IFI30), F-box protein 16 (FBXO16), calmodulin 1 (CALM1), MORC family CW-type zinc finger 3 (MORC3)] were selected to construct a predictive model. Patients with HCC were then divided into low- and high-risk groups using the best cut-off value, as determined by the X-tile software. Prognosis was significantly poorer in the high-risk than in the low-risk group (P=0.009; hazard ratio =2.65). Features related to exosomes were found to positively regulate immune response. Further analysis showed a higher risk score was associated with higher expression of immune checkpoint molecules, including programmed death ligand 1 (PD-L1), programmed death ligand 2 (PD-L2), T cell Ig and ITIM domain (TIGIT), and indoleamine-2,3-dioxygenase 1 (IDO1). Conclusions: This study has identified a novel signature based on exosome-related genes that has potential as a prognostic biomarker for HCC. Our research provides an immunological perspective for the development of precision treatment for HCC. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Hepatocellular carcinoma (HCC) accounts for 85–90% of primary liver cancers and is the fourth leading cause of cancer-related death worldwide. Approximately 841,000 new HCC cases are diagnosed worldwide each year (1,2). Despite new breakthroughs in targeted therapy, immune therapy, interventional therapy, surgical techniques, and liver transplantation, the prognosis of HCC remains poor due to its high rates of metastasis and recurrence (3-6). Traditional prognostic models for HCC use clinic-pathological parameters, such as stage, and other parameters; however, due to tumor heterogeneity among patients, these models no longer work satisfactorily. Therefore, it is necessary to investigate the underlying molecular mechanisms of HCC and to identify new diagnostic and prognostic markers for the disease (7,8).

Exosomes are tiny extracellular vesicles with a lipid bilayer membrane structure that were first discovered by Johnstone in 1989 in his study of reticulocytes (9). Later studies revealed that exosomes can transport biologically active molecules between cells to regulate the microenvironment between cells and the immune system, and are closely related to tumorigenesis and tumor development (3,10,11). Exosomes derived from tumor tissue contain numerous cancer-related biomarkers, such as microRNAs (miRNAs), which can be used to detect HCC at an early stage. Although some breakthroughs have been made in the field of exosome research, the specific biological functions of exosomes have yet to be fully elucidated (12-15). Therefore, the aim of this study was to comprehensively investigate exosome-related genes in HCC, which can help predict the outcome of HCC and discover potential new drug candidates for specific targeted therapy. We present the following article in accordance with the TRIPOD reporting checklist (available at https://atm.amegroups.com/article/view/10.21037/atm-21-4400/rc).

Methods

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Data collection

Genes related to liver-derived exosomes were downloaded from the exoRBase (http://www.exorbase.org/). Data on mRNA expression in 407 tumor and 58 normal tissue samples, and the corresponding clinical data, were obtained from The Cancer Genome Atlas (TCGA) (https://portal.gdc.cancer.gov/). For the validation dataset, data on gene expression in 240 tumor and 202 normal tissue samples were obtained from the International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/). The study flowchart is shown in .

Figure 1

Study flowchart. DEGs, differentially expressed genes; FC, fold change; ICGC, International Cancer Genome Consortium; LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas.

Identification of differentially expressed genes and competing endogenous RNA network construction

Differentially expressed genes (DEGs) related to exosomes and liver cancer were identified using the limma package in R software. Those DEGs with an adjusted P<0.05 and |log2(fold change)| >1 were selected for further analysis, as described below:

The miRcode database (http://mircode.org/) was used to search for interactions between long non-coding RNAs (lncRNAs) and miRNAs, and to match miRNA interactions;

The TargetScan (http://targetscan.org/) and miRanda (http://www.microrna.org/) databases were used to predict miRNA-targeted mRNAs;

The Encyclopedia of RNA Interactomes database (http://starbase.sysu.edu.cn/) was used to predict the relationships between circular RNAs (circRNA) and miRNAs.

A competing endogenous RNA (ceRNA) network comprising mRNAs obtained from these databases that overlapped with the exosome-related mRNAs was constructed. Results were graphed by Cytoscape 3.7.2. [Cytoscape is provided by the U.S. National Institute of General Medical Sciences (NIGMS)].

Development and validation of a prognostic model

First, we performed univariate Cox regression analysis in TCGA database and two-sided P<0.05 was set as statistically significant. The glmnet package in R software was used to narrow the genetic screening range, after which least absolute shrinkage and selection operator (LASSO) regression was performed (16). We performed 1,000 replicates on the dataset and selected markers with >500 times replicates. The regression coefficient (β) was obtained from the LASSO regression. The prognostic index (PI) was then calculated using the following formula:

The X-tile software was used to establish the optimal cut-off value to divide patients into high- and low-risk groups. The predictive power of the prognostic model was evaluated using Kaplan-Meier survival analysis, and the model was validated in the ICGC database.

Identification of independent prognostic parameters for HCC

Independent factors affecting clinical prognosis of patients with HCC were obtained by univariate and multivariate Cox proportional hazards regression analyses. Two-sided P<0.05 was considered significant. Hazard ratios (HRs) and 95% confidence intervals (CIs) were calculated for each factor.

Establishment and evaluation of a forecast nomogram

Independent prognostic factors obtained from the Cox analyses were used to construct a nomogram to predict the 1- , 2- and 3-year survival rates of patients with HCC. The nomogram was verified internally using a calibration chart, and its predictive ability was assessed through time-varying receiver operating characteristic (ROC) curve analysis (17,18).

Immune cell type scores and immunotherapy effect

CIBERSORT is a new analytical tool to characterize cell composition based on gene expression profiling which is frequently used to evaluate immune cell infiltration (19). In the present study, CIBERSORT was used to forecast differences in the ratios of immune cell types between the low- and high-risk patient groups. The Genefilter package of R was used to screen each sample, and P<0.05 was considered significant.

Immunotherapy is a very promising treatment method for patients with HCC, and its therapeutic effect is related to molecules including programmed death ligand 1 (PD-L1), programmed death ligand 2 (PD-L2), T cell Ig and ITIM domain (TIGIT), and indoleamine-2,3-dioxygenase 1 (IDO1). Differences in the expression of these molecules between the low- and high-risk groups were evaluated.

Statistical analysis

Cox regression and LASSO regression analyses were used to develop the prognostic model, and univariate and multivariate Cox proportional hazards regression analyses were used to identify independent prognostic parameters for HCC. Two-tailed P<0.05 was considered significant. Statistical analyses were performed with R 4.0.2.

Results

Construction of the ceRNA network

Using data downloaded from the exoRBase database, differences in mRNA, lncRNA, and circRNA expression between tumor and non-tumor tissues were identified. Relationships between mRNAs, lncRNAs, miRNAs, and circRNAs, and between lncRNA-miRNA, miRNA-mRNA, and circRNA-miRNA pairs were explored. Then, miRNAs were used to develop a ceRNA network with Cytoscape software. In all, 722 lncRNA-miRNA pairs, 680 miRNA-mRNA pairs, and 64 circRNA-miRNA pairs were obtained. Finally, 19 differentially expressed (DE) lncRNAs, 45 DEmiRNAs, 47 DEmRNAs, and 2 DEcircRNAs were used to construct a ceRNA network.

Construction of an exosome-related gene prognostic model in TCGA

Based on intersections between mRNAs found in the ceRNA network and TCGA genes, 48 DEGs were identified. Univariate Cox regression analysis revealed 18 DEGs that were significantly correlated with overall survival (OS) in patients with HCC (P<0.01). LASSO analysis was used to narrow down the range of genes. Five DEGS that appeared >500 times in a total of 1,000 replicates were selected, and these genes (XPO1, IFI30, FBXO16, CALM1, MORC3) were chosen to establish a prognostic model using the following equation:

X-tile indicated that the best cut-off value to discriminate between high- and low-risk patients was 2.06. The patients were divided into low- and high-risk groups, and Kaplan-Meier analysis revealed that OS was significantly poorer for the high-risk group than the low-risk group (P<0.001; HR 3.51; 95% CI: 1.68–7.31; ). The area under the time-dependent ROC curve for 3- and 5-year OS was 0.681 and 0.708, respectively, indicating the good predictive performance of the prognostic model ().

Figure 2

Establishment of a prognostic model of exosome-related genes in The Cancer Genome Atlas database and evaluation of the performance of the exosome-related prognostic model in the International Cancer Genome Consortium dataset. (A,C) Kaplan-Meier survival curves of the five prognostic genes in the testing (A) and validation (C) datasets. (B,D) Receiver operating characteristic (ROC) curve analysis of the predicted 3- and 5-year overall survival based on the five genes in the testing (B) and validation (D) datasets. AUC, area under the curve.

External validation of the exosome-related gene prognostic model in the ICGC

The predictive power of our prognostic model was validated using samples from the ICGC. As described above, patients were divided into low- and high-risk groups using a cut-off value of 2.06; the overall survival rate of patients in the low-risk group was higher than that of patients in the high-risk group (P=0.009; HR 2.65; 95% CI: 0.88–8; ). The area under the curve of the five-gene prognostic model for 3- and 5-year survival was 0.667 and 0.817, respectively (), indicating the reliability of the model.

Interpretation of the model genes

Gene set enrichment analysis functional annotation showed that the five HCC exosome-related genes in the prognostic model were enriched in the base excision repair, cell cycle, and glycolysis pathways (). As shown in , single deletions were found in XPO1 (2p15), FBXO16 (8p21.1), CALM1 (14q32.11), and MORC3 (21q22.13). Regression analysis between gene expression and methylation revealed that there were negative correlations between methylation and gene expression of CALM1, FBXO16, IFI30, and MORC3; however, no obvious correlation was found for XPO1 ().

Figure 3

Enrichment analysis, changes in copy numbers, and DNA methylation. (A) Gene set enrichment analysis of the five genes in the prognosis model. The graph (upper panel) shows the enriched pathways, whereas the lower panel shows expression levels. (B) Changes in the copy number of the five genes in patients with hepatocellular carcinoma. The outer circle shows the arrangement of chromosomes, whereas the inner markers are single deletions of four of the genes: XPO1 (2p15), FBXO16 (8p21.1), CALM1 (14q32.11), and MORC3 (21q22.13). (C-G) Regression analysis between gene expression and DNA methylation for the five genes in the prognostic model. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Establishing the predictive nomogram

Univariate Cox regression analysis revealed TNM stage and risk score to be significant predictors for survival (). Through multivariate Cox regression analysis, TNM stage and risk score were confirmed as independent prognostic factors ().

Figure 4

Construction of the predictive nomogram. (A,B) Forrest plots of the univariate (A) and multivariate (B) Cox regression analyses in HCC. (C) Building the nomogram for predicting the overall survival of patients with HCC in The Cancer Genome Atlas cohort. The nomogram plot was built based on two independent prognostic factors in HCC. (D-F) Calibration maps showing conformity between the predicted and observed results. AFP, α-fetoprotein; OS, overall survival; HCC, hepatocellular carcinoma.

Subsequently, a nomogram including TNM stage and the risk score was built. Using the nomogram, scores are assigned to independent factors at each level for individual patients, with the total score calculated by summing these scores. By transforming the relationship of the total score, the 1-, 2- and 3-year survival rates can be obtained (). A calibration map was used for internal verification of the nomogram, which showed that there was good conformity between the predicted and observed results ().

Computation of immune cell type fractions

CIBERSORT was used to directly analyze the types and expression of immune cells in tissues from patients with HCC, combined with the Leukocyte signature matrix (LM22). Differences in 22 types of immune cells were examined between the low- and high-risk HCC groups. Significant differences were found for two types of immune cells: regulatory T (Treg) cells (P=0.043) and M2 macrophages (P=0.048; ). Based on this result, we may conclude that high levels of Treg cells and M2 macrophages in patients with HCC are related to a poor prognosis. shows the proportion of each immune cell type in each sample, more intuitively showing the distribution of immune cells in different samples.

Figure 5

Relationship between the exosome-related prognostic signature and heterogeneity in immune cell infiltration. (A) Differences in immune cell infiltration between the high- and low-risk groups of patients with hepatocellular carcinoma. The data is shown as a violin plot, and the horizontal axis representing different immune cells and the vertical axis representing the degree of cell infiltration (*P<0.05; ns, not significant). (B) Histogram showing the relative infiltration of immune cell populations in tumor samples from The Cancer Genome Atlas dataset for which RNA-sequencing data were available.

Evaluation of immune-related DEGs

Blocking immune checkpoints is a promising method for treating many cancers. Therefore, in the present study we evaluated differences in the expression of key immune checkpoint molecules, including PD-L1, PD-L2, TIGIT, and IDO1, between low- and high-risk HCC patients. As shown in , the expression levels of immune checkpoint molecules were higher in the high-risk group than in the low-risk group (P<0.05). Based on this finding, we may conclude that high expression levels of exosome-related genes are related to a poor prognosis of HCC, which may account for patients’ heterogeneous responses to immunotherapy and help to identify specific patient populations that may benefit from immunotherapy.

Figure 6

Differences in the expression of the immunomodulatory molecules (A) PD-L1, (B) PD-L2, (C) TIGIT, and (D) IDO1 between the low- and high-risk groups of patients with hepatocellular carcinoma. And is shown as box plots, with the boxes indicating the interquartile range and the median value indicated by the horizontal line; whiskers show the range. Circles indicate individual values. PD-L1, programmed death ligand 1; PD-L2, programmed death ligand 2; TIGIT, T cell Ig and ITIM domain; IDO1, indoleamine-2,3-dioxygenase 1.

Discussion

Hepatocellular carcinoma is a fatal malignant tumor that has been a public health problem in many countries for many years. Although there has been considerable progress in the treatment of HCC, survival rates have not improved satisfactorily (19). Accumulating evidence shows that the malignant phenotype is affected by the cancer-related microenvironment. Hence, discovering immune biomarkers to predict the prognosis of patients with HCC is important and may also help with immunotherapy (20-22). In the present study, we used the TCGA database to develop a prognostic model based on exosome-related genes (XPO1, IFI30, FBXO16, CALM1, and MORC3) to predict the prognosis of patients with HCC, and then successfully validated the prediction ability of the model using an external database.

The genes used to build our prognostic model have been studied previously. Azizian et al. reported that XPO1 is frequently overexpressed in cancer cells and that its suppression leads to reductions in proteins such as MYC and epidermal growth factor receptor (20). XPO1 interacts with hundreds of proteins, which may affect its nuclear export activity. Originally identified as an autoantigen in inflammatory myopathies, MORC3 is recognized as a human ATPase and has been linked to cancer (21,22). CALM1 is a Ca2+ receptor protein which mediates a large number of signaling processes, including proliferation, motility, and differentiation. It had been found that CALM1 expression is strongly associated with many cancers (23). FBXO16 is a tumor suppressor that attenuates nuclear β‐catenin, and suppresses the growth, migration, and invasion of cancer cells (24). IFI30 has been found to be highly expressed in malignant glioma and to possibly have an influence on response to chemotherapy (25). Based on the results of the present study, IFI30 may affect the prognosis of HCC, but further studies are needed to illustrate the mechanism in detail. Furthermore, microenvironment analysis indicated that high IFI30 expression was accompanied by greater infiltration of M2 macrophages, which may be a research direction (25). Also, the five genes used in our prognostic model were found to be enriched mainly in the base excision repair, cell cycle, and glycolysis pathways, which have a clear relationship with tumor occurrence and development (26-29).

From another viewpoint, our model could be used to predict the prognosis of tumors and potential new treatment targets. Abnormal DNA methylation is well known to occur in all kinds of cancers and, in previous studies (30,31), we found that the expression of CALM1, FBXO16, IFI30, and MORC3 was negatively correlated with their methylation. Normal cells can transform to tumor cells through the occurrence of driver mutations and then the process of methylation; from this, it can be inferred that these genes play a role in tumorigenesis.

Liu et al. found that functional Treg cells facilitate the development of an M2 macrophage phenotype by repressing the CD8+ T cell-interferon γ axis and promoting the mitochondrial integrity of M2 macrophages via CD8+ T cells (32). Based on their study, it can be inferred that there is a positive correlation between Treg cells and M2 macrophages; however, in the present study, the opposite was observed, which may indicate that, in HCC, tumorigenesis and progression may occur as a result of changes in these pathways.

We confirmed that the prognostic model we developed can divide patients with HCC into low- and high-risk groups, and then showed that the oncologic outcomes of the low-risk group were better than those of the high-risk group. A nomogram was then built to predict the 1-, 2-, and 3-year survival rates. The calibration plot demonstrated that the resulting curve was very close to the real (observed data) curve, which means well predictive value. Thus, we further analyzed differences in infiltrating immune cells and immune-related gene expression between the low- and high-risk groups. In the present study, high-risk patients were more likely than low-risk patients to have higher levels of immune checkpoint molecules and a poor prognosis.

Many studies have focused on the role of exosomes in intercellular communication in HCC. It is widely accepted that HCC cells communicate with normal cells and promote the development and metastasis of HCC through exosomes. Some exosomes regulate the signal transduction pathway between cells, whereas others can be used as medicines due to the protective effects of the outer membrane (33,34). However, there is a lack of research on exosome-related genes for predicting survival and the effect of immunotherapy in HCC. In recent years, there have been breakthroughs in the treatment of liver cancer, especially in immunotherapy (35,36); the present study adds to the bank of exploratory research in this area and will hopefully contribute to future studies.

However, this study has some limitations. First, our study was retrospective and its findings need to be further verified in prospective studies, especially randomized control trials. Second, the LASSO analysis may have resulted in some important factors that contribute to the prognosis of HCC being overlooked. Third, the expression and the prognostic effects of the five genes used to develop the model require further investigation at the protein level. Fourth, functional experiments are needed to clarify it’s the potential mechanism. Finally, the detailed mechanism between tumor cells and immune cells requires further clarification.

Conclusions

This study has identified a novel signature based on exosome-related genes that has potential as a biomarker for the prognosis of HCC. These results may provide an immunological perspective for the development of precision treatment for HCC.

The article’s supplementary files as