Predictive Models for HCC Prognosis, Recurrence Risk, and Immune Infiltration Based on Two Exosomal Genes: MYL6B and THOC2

Introduction

Liver cancer is the sixth most commonly diagnosed cancer worldwide and the fourth leading cause of cancer death with high morbidity and mortality rates;¹ there were 905,677 new cases of liver cancer and 830,180 related deaths in 2020, according to the WHO statistics (https://gco.iarc.fr/today/fact-sheets-cancers). Hepatocellular carcinoma (HCC) accounts for most cases of liver cancer. HCC is a highly malignant tumor, its onset is insidious, and it progresses rapidly. Patients are often in advanced stages (BCLC stage B and stage C) at the time of diagnosis and cannot undergo surgery or transplantation, thereby leading to the average survival time of HCC is only 6 months.^2,3 The efficacy of radiofrequency ablation (RFA) and transcatheter arterial chemoembolization (TACE), which are proposed treatments for advanced HCC, is also unsatisfactory.⁴ HCC is characterized by molecular heterogeneity, and its molecular pathogenesis is very complex.⁵ Although many new treatments for HCC have been proposed in recent years, less progress has been made in significantly improving the outcome of HCC.^6,7 Therefore, there is an urgent need for a better understanding of the pathogenesis and molecular mechanism of HCC and the identification of novel key HCC markers to help determine more effective treatment strategies for patients.

Exosomes are a type of vesicle (30–100 nm) derived from endosomes and are shed by cells during cellular housekeeping and communication.⁸ Exosomes are important mediators of cell interactions, and they transmit information between cells through their membrane-enclosed structure and cytoplasmic content.⁹ Several studies have shown that the release rate and content of exosomes derived from cancer cells are significantly different from those originating from normal cells, and multiple types of cancer cell lines produce more exosomes than normal cells.¹⁰ During the progression of HCC, the secretion of cancer cell-derived exosomes has been shown to promote tumorigenesis, enhance angiogenesis and epithelial to mesenchymal transition, induce drug resistance and interfere with antitumor immune mechanisms, thereby promoting tumor progression and metastasis.^11–13 Exosomes play an important role in the carcinogenic process of HCC, leading to poor patient outcomes; however, studies to date have mainly focused on exploring the role of a single factor of exosomes in the progression of HCC.^14,15 No combination analyses have been conducted of exosome-derived molecules to assess the association between these molecules and the development and prognosis of HCC.

HCC is an immunogenic tumor that develops in a continuous immunomodulatory environment caused by chronic viral inflammation, such as hepatitis B (HBV) and hepatitis C (HCV) infections,¹⁶ and once HCC is recognized by the immune system, tumor cells are vulnerable to attacks by immune cells.^17,18 Immunotherapy has received widespread attention in the antitumor field as a novel treatment option for HCC. However, tumor cells can enhance immunosuppression through immune checkpoint expression changes and immune cell aggregation, which can lead to immune escape and promote HCC progression.^19–21 The immune activity of exosomes can affect the immunoregulatory mechanism of tumor cells, including the regulation of antigen presentation, immune activation, immune suppression, immune surveillance and immune escape.^22,23 Further exploration of the marker molecules of exosomes and their relationship with the immune microenvironment of HCC may provide important insights for immune recognition and therapeutic intervention.

The emergence of high-throughput sequencing technology has increased the convenience of analyzing and identifying marker molecules that have an important impact on the immune microenvironment of HCC and disease prognosis.²⁴ In this study, we aimed to identify the optimal gene set derived from exosomes in HCC with substantial predictive value through comprehensive analyses to construct prognostic, recurrence risk and diagnostic models and identify candidates suitable for immune checkpoint inhibitor immunotherapy and chemotherapy, thereby guiding individualized patient treatment options. Our results will help improve the early detection rate of HCC, develop individualized therapeutic plans and improve disease outcomes.

Materials and Methods

Identification of Differentially Expressed Genes (DEGs) between HCC and Adjacent Nontumor Tissues

The messenger RNA (mRNA) sequences of HCC patients with available clinical information were obtained from the TCGA (http://www.tcga.org/), ICGC (https://dcc.icgc.org/) and GEO databases (https://www.ncbi.nlm.nih.gov/geo/). The patients’ baseline character in the three cohorts were listed in Table S1. Exosome-derived mRNA sequences were downloaded from the exoRbase database (http://www.exorbase.org/), which is a repository containing information such as circular RNA (circRNA), long noncoding RNA (lncRNA) and mRNA sequences derived from RNA-seq data analysis of human blood exosomes, and the data sources include both normal individuals and patients with different diseases. Using the limma R package, DEGs with a threshold of absolute log2-fold change (FC) > 1 and an adjusted P value < 0.05 were chosen.

Weighted Correlation Network Analysis (WGCNA)

WGCNA²⁵ was used to identify the HCC exosome-specific module in the exoRbase database (http://www.exorbase.org/). First, the R package “WGCNA” was applied to build an expression matrix of exosome-derived mRNAs in the exoRbase database. The best “soft threshold capability” was determined graphically as β=9 to ensure a scale-free distribution network, and the expression matrix and β value were adopted to construct a coexpression matrix. The DynamicTreeCut algorithm was employed to classify genes with similar expression levels into the same gene modules to generate coexpression modules. The eigengene function in the “WGCNA” R package was used to calculate the correlation between the module eigengenes (ME) of the module and HCC, and heat maps were used to visualize the correlation between each module and HCC. The “ME greenyellow” module with the highest correlation with HCC was identified as the HCC exosome-specific module, and the Pearson coefficient was calculated to evaluate the correlation between the genes in the module and HCC.

Construction and Validation of a Predictive Prognostic Model Based on the Identified Key Exosome-Derived Genes

Univariate Cox regression analysis was performed to identify exosome-derived genes related to the prognosis of patients with HCC. Genes with P<0.05 were considered statistically significant. The identified genes and other genes that the univariate Cox regression analysis considered not to be statistically significant but that have been clinically identified to have statistical significance in the prognosis of HCC patients were incorporated in subsequent multivariate Cox regression analysis. Through forward selection and backward elimination selection, the optimal gene set was obtained to construct a prognostic model. The prognostic index (PI) of the prognostic model was calculated as follows: (PI) =Σⁿ_i=1(β_{i *} e_i). where β_i represents the correlation coefficient between the gene (i) and HCC prognosis, and e_i represents the expression level of gene. The expression level of genes can be obtained in HTSeq-FPKM in the TCGA database. HCC patients with known survival status were divided into high- and low-risk groups based on the optimal cutoff value determined by X-tile software. Kaplan-Meier (K-M) curve and time-dependent receiver operating characteristic (ROC) curve analyses were performed to evaluate the predictive performance of the prognostic model.

Establishment and Evaluation of a Predictive Nomogram

Univariate and multivariate Cox regression analyses were performed to assess whether the prognostic model can independently predict the prognosis of HCC patients compared to traditional clinical variables (including age, AFP, weight, vascular tumor cell, sex, pathological grade and TNM stage). The hazard ratio (HR) and 95% confidence interval (CI) for each variable were calculated. Variables with P<0.05 were considered statistically significant independent prognostic factors. Rms R software was used to integrate these independent prognostic factors to establish a predictive nomogram and corresponding calibration maps. Calibration and identification were used for the validation of calibration maps. The closer the calibration curve was to the reference line (at 45°), which represents the best prediction value, the better the prediction performance of the nomogram. The consistency index (C Index), calculated using a bootstrap method with 1000 resamples, was also adopted to validate the consistency between the predicted results of the nomogram and the actual results. In addition, ROC curve analysis was performed to evaluate the predictive value of the nomogram compared to a single independent prognostic factor, and decision curve analysis (DCA) was adapted to measure the difference in clinical benefit that could be achieved using the nomogram for prognosis prediction compared to using a single prognostic factor.²⁶

Internal and External Validation of the Expression Profiles of Exosome-Derived Genes

Using Wilcoxon signed-rank tests in Prism 7.0 (GraphPad, San Diego, CA, USA), we compared the expression profile of a single gene in HCC tissues and normal tissues in the TCGA HCC cohort and confirmed gene expression differences in HCC tissues and nodule samples in the GSE6764 HCC cohort. P<0.05 was considered statistically significant. The K-M curve was used to evaluate the influence of the expression of a single gene on the prognosis, and ROC curve analysis was performed to validate the predicted influence of the gene. For external validation, we further confirmed the effect of single gene expression on prognosis and recurrence risk in the Kaplan-Meier Plotter database.²⁷ The effect of single gene expression on immune cell infiltration in tumors was also evaluated.

Estimating the Fractions of Different Types of Immune Cells

CIBERSORT analysis was adopted to quantitatively analyze the absolute abundance of 22 immune cell types in heterogeneous tissue to evaluate the infiltration of immune cells.²⁸ The R package “CIBERSORT” was used to convert the mRNA data from nontumor cells into the tumor microenvironment infiltration level.

Cell Culture and Transfection

SK-HEP1 and HCCLM3 cell lines were purchased from the American Type Culture Collection (ATCC) (Manassas, VA, USA). The cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 100 U/mL penicillin-streptomycin in an incubator with an atmosphere of 5% CO₂ and 95% relative humidity at 37°C.

The transfection process is briefly described as follows: MYL6B siRNAs and THOC2 siRNAs were synthesized by GeneChem (Shanghai, China), and Lipofectamine 2000 (Invitrogen, California, USA) was used to achieve cell transfection with these siRNAs. SK-HEP1 and HCCLM3 cell lines were resuspended in DMEM supplemented with 10% FBS and cultured in a 6-well plate at a density of approximately 40% for 12 hours. One hundred pmol of siRNA and 5 µL of Lipofectamine 2000 were dissolved in 250 µL of Opti-mem serum-free medium, mixed and allowed to stand for 20 minutes. The cells were washed with PBS 3 times, and 1.5 mL Opti-men and the mixture containing the corresponding siRNA and Lipofectamine 2000 were added to the cells. After 6 hours, the mixture was removed and the cells were then cultured in medium supplemented with serum. The sequences of the knockdown and control siRNAs are listed in the supporting information.

Western Blot Assay

RIPA lysis buffer (Invitrogen, Carlsbad, USA) supplemented with 1% protease and phosphatase inhibitor (Thermo, Shanghai, China) was used to obtain protein from SK-HEP1 and HCCLM3 cells. The separation of protein samples with different molecular weights was achieved by 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE). The separated protein was electrotransferred to a polyvinylidene fluoride (PVDF) membrane (Invitrogen), which was then blocked with 5% skimmed milk at room temperature for 90 minutes, incubated with the primary antibody at 4 °C overnight, and subsequently incubated with the HRP-conjugated secondary antibody for 2 hours at room temperature. The iBright FL1500 intelligent imaging system (Invitrogen) was employed to detect the signal and image the blot. The antibody information is listed in the supporting information.

Cell Proliferation Assay

Cell counting kit-8 (CCK-8) assays (Dojindo, Tabaru, Mashikimachi, Japan) and 5-ethynyl-2ʹ-deoxyuridine (EdU) assays (Ruibo, Guangzhou, China) were employed to measure cell proliferation. The CCK-8 assay was performed as follows: SK-HEP1 and HCCLM3 cells were seeded in 96-well plates with 3000 cells and 100 μL DMEM per well. After the cells were incubated for 0, 24, 48 and 72 hours, the medium was changed to fresh medium supplemented with 10% CCK-8 reagent and the cells were incubated for another 2 hours at 37°C. The absorbance of the cells at 450 nm was detected using a microplate reader (Bio-Rad, Berkeley, CA, USA). Each group included three replicate wells. For the EdU assay, we seeded SK-HEP1 and HCCLM3 cells in a 96-well plate at 3000 cells per well. After 24 hours of incubation, 100 μL 50 μM EdU reagent was added to each well, and the incubation was continued for another 2 hours. Then, the cells were fixed with 4% polyoxymethylene for 30 minutes at room temperature and incubated with Apollo staining solution and Hoechst 33342 for 30 minutes. A fluorescence microscope was used to detect cell proliferation.

Cell Migration Assay

A Corning Costar Transwell chamber was used to conduct the cell migration assay. Serum-free DMEM (Gibco, NY, USA) was added to the upper chamber, and DMEM supplemented with 10% fetal bovine serum was added to the lower chamber. The transfected cells were trypsinized, centrifuged and resuspended in the medium in the upper chambers at a density of 1×10⁵ cells per well. After incubation at 37°C for 48 hours, non-migrated cells from the upper side of the upper chamber were removed by a cotton swab, and cells that migrated into the lower side of the upper chamber were stained with 0.1% crystal violet and counted under an optical microscope. Each experiment was repeated three times.

Statistical Analysis

GraphPad Prism 7 software (version 8.0) and R language (version 3.6.3) were applied to statistically analyze the experimental data and generate figures. The quantitative data are presented as the mean ± standard error of the mean (SEM). The comparison of the mean of two independent samples was performed by Student’s t-test or Wilcoxon test. P values less than 0.05 were considered statistically significant.

Results

Identification of Exosome-Derived Genes Strongly Associated with HCC Prognosis

First, we drew a flowchart to vividly depict our experimental process (Figure S1). Through WGCNA, we divided the exosome-derived mRNAs in the exoRbase database into trait modules to collect genes with similar traits (Figure S2), thereby determining HCC exosome-specific modules that contained 168 hub genes (Figure 1A and B). Then, we identified 5775 DEGs in the TCGA database and 150 DEGs in the exoRbase database. Finally, a total of 19 HCC-related exosome-derived hub genes were incorporated in the subsequent analysis through integrated analysis (Figure 1C). The correlation between genes is shown in Figure 1D. GO analysis performed on these genes indicated that the genes were mainly enriched in “cytosolic ribosome”, “polysomal ribosome” and “ribosomal subunit” (Figure 1E). Through the univariate Cox regression analysis, we found that seven genes were closely related to the OS rate of HCC patients (P<0.05), and the subsequent multivariate Cox regression analysis revealed an optimal gene set consisting of two genes (MYL6B and THOC2) for predicting HCC prognosis (Table S2).

Figure 1 Identification of key exosome-derived genes in HCC. (A) Modules with different traits determined via WGCNA. (B) Identification of genes in the exosome-specific module in HCC. (C) Determination of key exosome-derived genes in HCC. (D) Correlation between exosome-derived genes. (E) GO analysis of the exosome-derived genes. WGCNA, weighted gene co-Expression network analysis. Abbreviations: GO, Gene ontology; WGCNA, weighted gene co-Expression network analysis.

Establishment of a Prognostic Model and Validation of Its Predictive Performance

Based on the two exosome-derived genes, we constructed a prognostic model. The prognostic index was calculated as PI = (0.0273 * expression level of MYL6B) + (0.1931 * expression level of THOC2). The optimal cutoff value determined by X-tile software was used as the threshold to divide patients into a high-risk group and a low-risk group. The optimal cutoff was determined to be 1.331 in the HCC cohort in the TCGA database, 3.104 in the ICGC database and 1.863 in the GSE14520 dataset. The HCC cohort in the TCGA database was adopted as the training cohort, and the HCC cohorts in the ICGC database and the GSE14520 dataset were used as validation cohorts to evaluate the predictive performance of the prognostic model. Patients in the high-risk group had a worse prognosis than patients in the low-risk group in both the training cohort (P<0.001) (Figure 2A) and the validation cohorts (P<0.01) (Figure 2D and G). Figure 2B, E and H show the expression levels of the two exosome-derived genes in the patients and their corresponding risk scores. The area under the curve (AUC) at 0.5, 1, 3, and 5 years was 0.66, 0.68, 0.64, and 0.59, respectively, in the training cohort, indicating good predictive performance of the prognostic model (Figure 2C). The AUCs of the two validation cohorts at 0.5, 1, 3 and 4 years reached 0.72, 0.72, 0.68 and 0.63 (Figure 2F) and 0.66, 0.61, 0.62 and 0.59 (Figure 2I), respectively, confirming the good predictive value of the prognostic model.

Figure 2 K-M curve analysis, risk score distribution and time-dependent ROC curve analysis of the prognostic model in the HCC cohorts from TCGA (A–C), ICGC (D–F) and GSE14520 (G–I). (A, D and G) Shorter survival times were presented in the high-risk group than in the low-risk group (P<0.01). (B, E and H) Distribution of prognostic risk scores based on the expression levels of the exosome-derived genes. (C, F and I) Time-dependent ROC curves of the prognostic model. Abbreviation: ROC curve, Receiver operating characteristic curve.

Subsequently, the correlation between the risk scores obtained with the prognostic model and clinical variables was explored. Figure S3A shows the distribution of prognostic risk scores of patients stratified by traditional clinical variables and the expression profiles of the two exosome-derived genes. Patients with pathological grade G3-G4 HCC had a higher prognostic risk index than patients with G1-G2 HCC (Figure S3B) (P<0.05). The prognosis of patients with vascular tumor invasion was also significantly worse than the prognosis of those without vascular tumor invasion (P<0.05) (Figure S3C), consistent with clinically accepted diagnostic criteria. However, there was no significant difference in the prognostic risk among HCC patients with different TNM stages (Figure S3D).

Chemotherapy Response of Patients with Different Prognostic Risk Scores

Chemotherapy is one of the common treatment options for HCC patients. It is a suitable treatment option for HCC patients who cannot undergo surgery to reduce tumor volume and prolong survival time. It is also commonly used in postoperative adjuvant chemotherapy to suppress tumor recurrence and progression.²⁹ The resistance of tumor cells to chemotherapeutics is a key factor that affects the efficacy of anticancer drugs. We adopted the prognostic model to evaluate the response of HCC patients to 266 traditional chemotherapeutic drugs and molecular targeted drugs using data available from the Genomics of Drug Sensitivity in Cancer (GDSC) website (https://www.cancerrxgene.org/). The half-maximum inhibitory concentration (IC50) was used as the criterion for tumor sensitivity to drugs. Figure 3A–L shows the difference between the sensitivity of high-risk patients and low-risk patients to drugs such as methotrexate, erlotinib and vinorelbine. High-risk patients were more resistant to chemotherapy than low-risk patients (P<0.05), suggesting that low-risk patients may respond better to chemotherapy. Through gene set enrichment analysis (GSEA), we explored the signaling pathways that the exosome-derived prognostic model mainly interferes with. Figure 3M presents the top 5 positively regulated pathways, and Figure 3N shows the top 5 negatively regulated pathways. We found that the prognostic model was involved in the regulation of several signaling pathways, including “DNA replication”,³⁰ “homologous recombination”,³¹ and “mismatch repair”,³² which were determined to be associated with cancer cell chemoresistance in previous studies. The expression difference in these signaling pathways may partly explain the difference in chemoresistance among patients with different prognostic risk.

Figure 3 Response of high-risk patients and low-risk patients to chemotherapy. (A–L) Sensitivity of high-risk patients and low-risk patients to traditional chemotherapy drugs and molecular targeted drugs. (M) Top5 signaling pathways positively regulated by the exosome-derived prognostic signature. (N) Top5 signaling pathways negatively regulated by the exosome-derived prognostic signature.

Characteristics of the Immune Microenvironment Between High- and Low-Risk HCC Patient Groups

The blockade of immune checkpoints or pathways is an important form of immunotherapy and a current focus of antitumor therapy research.³³ The relative proportion of immune cell infiltration and corresponding prognostic risk scores of 22 immune cell types from the TCGA HCC cohort are shown in Figure 4A. Compared with that in the low-risk group, a higher proportion of memory B cells, M0 macrophages and follicular helper T cells (Figure 4B–D) and a lower proportion of NK resting cells (Figure 4E) were found in the high-risk group (P<0.05). Figure 4F shows the correlation between the prognostic model and the immune checkpoint expression levels. In the high-risk group, the expression levels of the immune checkpoints PD1, B7H3, B7H5, CTLA4 and TIM3 were significantly higher than those in the low-risk group (P<0.05) (Figure 4G–K). Immune checkpoints, including PD1 and CTLA4, play a vital role in regulating the activation of T cells and exerting immune function and have been found to be overexpressed in multiple cancers, such as HCC, leading to the inhibition of T cells and the promotion of tumor progression.^33,34 Targeted suppression of these abnormally expressed checkpoints effectively improves the immune response mediated by T cells, prompting the immune system to recognize, mobilize against, and eliminate cancer cells. The above findings indicated that patients with higher risk have high immune checkpoint expression levels, and the immune response is significantly improved by immunotherapy, indicating that better clinical benefit can be achieved with immunotherapy for these patients.

Figure 4 Correlation of immune cell infiltration and immune checkpoint expression with the prognostic risk in HCC patients. (A) The infiltration of immune cells in HCC patients and the corresponding prognostic risk score and survival status. (B–E) Differences in the infiltration of memory B cells (B), M0 macrophages (C), follicular helper T cells (D) and resting NK cells (E) in the high- and low-risk HCC patient groups. (P<0.05) (Wilcoxon test). (F) Interaction between prognostic risk and immune checkpoints including PDL1, PD1, CTLA4, B7H3 and TIM3. (G–K) The expression levels of PD1 (G), B7H3 (H), B7H5 (I), CTLA4 (J) and TIM3 (K) in the high-risk group were higher than those in the low-risk group. (P<0.05) (Wilcoxon test).

Establishment and Evaluation of a Corresponding Nomogram in the TCGA HCC Cohort

Patient in the TCGA HCC cohort with available clinical information were included in univariate and multivariate Cox regression analyses to assess the independent predictive performance of the prognostic model compared to that of clinical variables. The results of Cox regression analysis determined that age (HR = 1.687, P<0.05), TNM stage (HR = 1.947, P<0.05) and the prognostic model (HR = 1.792, P<0.05) could independently predict the prognosis of HCC patients (Figure 5A). A predictive nomogram integrating these independent prognostic factors was then constructed to quantitatively predict the prognostic risk of HCC patients (Figure 5B). The calibration curves for 1, 3 and 5 years are shown in Figure 5C–E. The C index of the nomogram reached 0.63, indicating a better prediction consistency with the actual results compared to the use of the single independent prognostic factor age (0.57), TNM stage (0.53) or the prognostic model (0.59). The AUC of the ROC curves at 1, 2, 3, and 5 years further confirmed that the predictive value of the nomogram was superior to any single independent prognostic factor (Figure 5F–H). The results of the DCA suggested that, compared with the use of a single independent prognostic factor, the best clinical benefit was obtained by using the nomogram to predict the prognosis of HCC patients (Figure 5I–K). The above results determined that the nomogram is suitable for clinically predicting the prognosis of HCC patients.

Figure 5 Construction and validation of a nomogram. (A) The independence of the prognostic model and clinical variables for prognosis prediction. (B) A nomogram was constructed to predict 1, 3, and 5 years of OS in HCC patients. (C–E) Calibration curves of the nomogram. (F–H) ROC curves validating the prediction specificity of the nomogram. (I–K) DCA indicated that the nomogram achieved the optimal clinical benefit. OS, overall survival. ROC curve, receiver operating characteristic curve. Abbreviations: DCA, Decision curve analysis; OS, overall survival; ROC curve, receiver operating characteristic curve.

Construction and Validation of a Recurrence Risk Model Based on the Two Exosome-Derived Genes

A recurrence risk model integrating the same two exosome-derived genes was constructed. The relapse index was defined as RI = (0.0236 * MYL6B expression level) + (0.0899 * THOC2 expression level). Patients from the TCGA HCC cohort with recurrence information available were divided into a high-risk group and a low-risk group using the optimal cutoff value, which was determined to be 1.272 using X-tile software. The recurrence probability of patients in the high-risk group was significantly higher than that of patients in the low-risk group (HR = 2.44, P<0.001) (Figure 6A). Figure 6B shows the expression levels of the two exosome-derived genes in the patients and their predicted recurrence risk. ROC curves indicated great specificity and sensitivity of the recurrence risk model (Figure 6C). Patients from the GSE14520 HCC cohort with recurrence information available were used as a validation cohort for the predictive performance of the recurrence risk model. The optimal cutoff value was 1.403. A significantly higher probability of recurrence was found in the high-risk group than in the low-risk group (HR = 1.54, P<0.05) (Figure 6D). The expression profiles of the two exosome-derived genes in HCC patients and the prediction of recurrence risk are shown in Figure 6E. The AUCs of ROC at 0.5, 1, 3, and 5 years were 0.68, 0.64, 0.58 and 0.57, respectively (Figure 6F).

Figure 6 K-M curve analysis, risk score distribution and time-dependent ROC curve analysis of the recurrence model in the HCC cohorts from TCGA (A–C) and GSE14520 (D–F). (A and D) The high-risk group was more likely to experience recurrence than the low-risk group (P<0.05). (B and E) Expression of MYL6B and THOC2 and the corresponding prognostic risk scores of patients. (C and F) ROC curves validating the predictive value of the recurrence model. Abbreviation: ROC curve, Receiver operating characteristic curve.

The characteristics of clinical variables and expression of the two exosome-derived genes in the TCGA HCC patients and the corresponding risk scores predicted by the recurrence risk model are shown in Figure S4A. Patients with pathological grade G3-G4 HCC were more likely to experience recurrence than patients with G1-G2 HCC (P<0.05) (Figure S4B). The probability of recurrence of HCC in patients with vascular tumor invasion was also significantly higher than that in HCC patients without vascular tumor invasion. (P<0.05) (Figure S4C), while Figure S4D shows that there were no significant differences in recurrence probability between patients with TNM stage III–IV and those with TNM stage I–II.

A Nomogram Integrating Independent Recurrence Factors Was Constructed and Validated

An independence analysis through Cox regression determined that age (HR = 2.253, P<0.05), TNM stage (HR = 2.433, P<0.05) and the recurrence model (HR =4.907, P<0.05) were independent predictive factors of recurrence in HCC patients (Figure 7A). Based on these independent predictors, a corresponding nomogram was generated (Figure 7B). The calibration curves at 1 year, 3 years, and 5 years indicated that the predicted results of the nomogram were consistent with the actual results (Figure 7C–E). The C index of the nomogram (0.70)> age (0.61)> recurrence model (0.60)> TNM stage (0.56) indicated that the clinical value of the nomogram to predict the recurrence probability of HCC was superior to that of a single independent predictive factor. The results of the ROC curves were also consistent with the above conclusions (Figure 7F–H). The DCA confirmed that the clinical benefit of the nomogram for recurrence prediction was better than that of any single independent predictor (Figure 7I–K).

Figure 7 Establishment and evaluation of a nomogram. (A) Independent analysis of the recurrence model compared to clinical characteristics. (B) A nomogram was constructed to quantitatively assess the probability of recurrence. (C–E) Calibration curves of the nomogram. (F–H) ROC curves determined that the prediction performance of the nomogram was better than that of a single independent predictor. (I–K) DCA determined that the prediction of the nomogram could achieve the maximum clinical benefit. Abbreviation: DCA, Decision curve analysis.

Diagnostic Models Were Built and Demonstrated to Have Superior Predictive Power

A diagnostic model based on the two exosome-derived genes was constructed to accurately distinguish between patients with early-stage HCC and normal individuals. The diagnostic equation determined by the stepwise logistic regression method was defined as follows: logit (P = HCC) = −23.0370 + (2.2566 × expression level of MYL6B) + (3.8247 × expression level of THOC2). Applying the diagnostic model to the paired HCC samples from the TCGA cohort (consisting of 50 normal samples and 50 HCC samples) resulted in a sensitivity of 76.00% and a specificity of 82.00% (Figure 8A). Figure 8B shows the consistency between the disease outcome predicted by the diagnostic model based on the expression levels of MYL6B and THOC2 and the actual outcome. The AUC of the diagnostic model was 0.8792, confirming the superior performance of the diagnostic model for identifying HCC (Figure 8C). The ICGC HCC cohort (consisting of 202 normal samples and 243 HCC samples) was adopted to validate the diagnostic performance. The sensitivity of the diagnostic results was 76.95%, the specificity was 85.64% (Figure 8D), and the AUC was 0.9077 (Figure 8F), indicating that the diagnostic model can accurately distinguish HCC patients from normal individuals. The consistency between the predicted results based on the expression of MYL6B and THOC2 and the actual results is shown in Figure 8E.

Figure 8 Diagnostic models constructed to distinguish HCC patients from normal individuals (A–F) and HCC from dysplastic nodules (G–L). (A and D) Confusion matrix of binary results indicated that the diagnostic model accurately distinguished HCC patients and normal individuals. (B and E) The expression profiles of MYL6B and THOC2 in HCC patients and the consistency between the predicted disease status and actual disease status (C and F) ROC curves determined the superior predictive performance of the diagnostic model. (G and J) The diagnostic model distinguished HCC and dysplastic nodules with high specificity and sensitivity. (H and K) Expression characteristics of MYL6B and THOC2 in HCC patients, distribution of the predicted disease status and actual disease status. (I and L) ROC curves validated the predictive value of the diagnostic model. Abbreviation: ROC curve, Receiver operating characteristic curve.

The early diagnosis of HCC is mainly based on imaging and pathological examination, but small nodules <2 cm often lead to missed diagnoses because of difficulties in accurate characterization.³⁵ We also constructed a diagnostic model based on MYL6B and THOC2 to distinguish HCC and dysplastic nodules. The diagnostic formula was determined as follows: logit (P-HCC) = −45.5786 + (2.6004 × MYL6B expression level) + (4.0578 × THOC2 expression level). We adopted patient information from the GSE6764 dataset as the training set and the GSE89377 dataset as the validation set to evaluate the predictive value of the diagnostic model. The diagnostic model in the training set reached 91.43% sensitivity and 76.47% specificity (Figure 8G), and the AUC was 0.9328 (Figure 8I). The sensitivity and specificity in the validation set were 82.50% and 50.00% (Figure 8J), respectively, and the AUC was 0.7864 (Figure 8L), indicating that the diagnostic model had good predictive performance. The gene expression profile of HCC patients and the consistency between the predicted disease state and the actual disease state are shown in Figure 8H and K.

The Expression Characteristics of Exosome-Derived Genes Were Validated to Accurately Predict the OS Outcomes of Patients with HCC

To evaluate whether MYL6B and THOC2 were suitable for constructing prognostic, recurrence risk and diagnostic models, we examined the expression profiles of these two genes. MYL6B and THOC2 were obviously more highly expressed in the TCGA HCC cohort samples than in the paired normal samples (Figure 9A and B). The expression levels of MYL6B and THOC2 in HCC patients from the GSE6764 dataset were significantly higher than those in dysplastic nodule samples (Figure 9C and D). The K-M curves indicated that patients with high gene expression had a worse prognosis than patients with low gene expression (Figure 9E and F), and ROC curve analysis confirmed the specificity of MYL6B and THOC2 for predicting patient prognosis (Figure 9G and H). We further explored the association between gene expression, patient prognosis and recurrence risk using the Kaplan-Meier Plotter database. The results confirmed that high gene expression was closely related to shorter overall survival times and progression-free survival times (Figure 9I - L). Figure 9M and N shows the correlation between the expression profiles of MYL6B and THOC2 and immune cell infiltration. In addition, the correlation of the expression of MYL6B and THOC2 and immune checkpoint genes was also explored (Figure S5A and S5F), and compared to patients with low expression of MYL6B and THOC, patients with high expression of MYL6B and THOC tended to have higher expression patterns of immune checkpoint genes, including PD-1, B7H3, CTLA4 and TIM3 (Figure S5B–E and S5G–J), indicating that patients with high expression of MYL6B and THOC are more likely to have an increased immune response and inhibition of tumor progression in response to targeted immune checkpoint immunotherapy. The above results determined the effective predictive value of the two exosome-derived genes MYL6B and THOC2.

Figure 9 The expression characteristics of MYL6B and THOC2 and their correlation with patient prognosis, recurrence, and immune infiltration. (A, B) MYL6B (A) and THOC2 (B) expression was higher in the HCC tissues than in the paired normal tissues in TCGA. (P<0.05) (Wilcoxon test). (C, D) Higher expression of MYL6B (C) and THOC2 (D) was found in HCC tissues than in dysplastic nodules in GSE6764. (P<0.05) (Wilcoxon test). (E, F) Survival analysis of MYL6B (E) and THOC2 (F). (G, H) ROC curves validating the prognostic value of MYL6B (G) and THOC2 (H). (I–L) Shorter survival times (I and J) and progression-free survival times (K and L) were found in patients with high expression of MYL6B and THOC2 than in patients with low expression of MYL6B and THOC2 in the Kaplan-Meier Plotter database. (M, N) Correlation between the expression of MYL6B (M) and THOC2 (N) and immune infiltration.

Figure 10 MYL6B and THOC2 suppress the proliferation and migration of HCC cells. (A–D) Western blot analysis confirmed that the expression of MYL6B and THOC2 was inhibited in the HCC cell lines with transfection of MYL6B and THOC2 siRNA. (E–H) CCK8 assay showed that MYL6B and THOC2 inhibition significantly inhibited the proliferation of HCC cell lines. ***p < 0.001(Student’s t-test). (I and J) The inhibitory effect of MYL6B and THOC2 inhibition on the proliferation of HCC cell lines was also validated in EdU assay. (K and L) The migration of SK-HEP1 cells and LM3 cells was significantly inhibited under the inhibition of MYL6B and THOC2. (M–P) Quantitative statistical results of the function of MYL6B and THOC2 inhibition on the migration of HCC cell lines. **p < 0.01(Student’s t-test). Abbreviation: si-CTR, Control treatment.

Determination of the Oncogenic Effect of MYL6B and THOC2 in HCC Cells

To further validate the functions of MYL6B and THOC2 in HCC cell proliferation and migration, we knocked down the expression of MYL6B and THOC2 in HCC cell lines (including SK-HEP1 and HCCLM3 cells). The results of Western blotting confirmed that the application of MYL6B siRNA and THOC2 siRNA significantly inhibited the expression of MYL6B and THOC2 in HCC cells (Figure 10A–D). The CCK assay (Figure 10E–H) and EdU assay (Figure 10I and J) both revealed that MYL6B inhibition and THOC2 inhibition had a significant inhibitory effect on the proliferation of SK-HEP1 cells and HCCLM3 cells. In addition, the Transwell assay showed that MYL6B inhibition (Figure 10K) and THOC2 inhibition (Figure 10L) obviously suppressed the migration of SK-HEP1 cells and HCCLM3 cells. Figure 10M–P shows the quantitative statistical analysis of the Transwell assay data. The above findings confirm that the inhibition of MYL6B and THOC2 suppressed the proliferation and migration of HCC cells, suggesting that the expression of MYL6B and THOC2 is closely related to the progression of HCC and that MYL6B and THOC2 are valuable for predicting the prognosis of HCC.

Discussion

As one of the most common malignant tumors, HCC represents a significant health and economic burden.³⁶ HCC involves a series of complex genetic and epigenetic alterations.³⁷ Although a variety of staging systems have been proposed and used in clinical decision-making to predict patient prognosis, these systems are mainly based on clinicopathological characteristics without considering the key role of complex molecular pathogenic mechanisms in the carcinogenesis and progression of HCC,^38,39 and patient outcomes have not been significantly improved.⁴⁰ Exosomes are the main type of extracellular vesicle (EV) that can mediate the transfer of protein, DNA and various forms of RNA between cells and participate in cell-to-cell communication.⁴¹ The regulatory dysfunction of exosomes along with multiple complex molecular mechanisms promotes the development of HCC.⁴² The increased level of exosomes released by tumor cells enhances the communication between HCC cells and non-HCC cells. Identifying key markers derived from exosomes and performing a combination analysis may help accurately assess the progression of HCC and predict the prognosis of HCC patients and disease recurrence risk. In this study, we used high-throughput sequencing technology to explore the potential mechanisms of exosome-related DEGs. Then, based on two exosomal genes (MYL6B and THOC2), we constructed an exosomal gene signature with high predictive performance. Then, we found that the overall survival and disease-free survival times were obviously shorter in HCC patients in the high-risk group than in patients in the low-risk group. The risk score of the exosomal gene signature served as an independent prognostic predictive factor. In addition, combined with other clinical features of HCC, a nomogram was built with superior predictive value in HCC patients. These results suggested that the exosomal gene signature may represent new possibilities for developing individualized treatments for HCC patients.

MYL6B is an essential light chain for non-muscle myosin II (NMII) and is involved in cell adhesion, migration, material transport and endocytosis control.^43,44 MYL6B binds to the p53 protein to accelerate p53 degradation and promote the development of HCC.⁴⁵ THOC2 is involved in mediating mRNA export from the nucleus to the cytoplasm, chromosomal arrangement, mitotic process, and genome stability.⁴⁶ The expression of THOC2 has been found to promote the proliferation and invasion of melanoma.⁴⁷ Our subsequent cell experiments also determined that MYL6B and THOC2 inhibition significantly suppressed the proliferation and migration of HCC cells, revealing that the expression of these two genes is positively associated with the progression of HCC. Based on MYL6B and THOC2 expression, the diagnostic model constructed was also determined to be able to distinguish HCC patients from normal individuals and HCC from dysplastic nodules with high specificity and high sensitivity. Our findings lay a foundation for adopting MYL6B and THOC2 as biomarkers for the early diagnostic evaluation of HCC and indicate that the exosomal gene signature may serve as an effective method for diagnosing HCC, especially in the early stages.

Exosomes are also involved in the regulation of the tumor immune response.⁴⁸ Tumor cell exosomes can promote the production of prostaglandin E2, IL-6 and TGF-β by myeloid-derived suppressor cells (MDSCs), thus forming a powerful immune suppression environment in tumor lesions.⁴⁹ Exosomes are also considered to be important mediators of antitumor immune responses and immune surveillance evasion by tumor cells. Some studies have shown that HCC-derived exosomes can reduce the cytotoxicity of T cells and NK cells and promote the aggregation of immunosuppressive M2 macrophages and N2 neutrophils.^50,51 Characterizing the expression of exosome markers may help predict patient response to immunotherapy and identify more effective treatments. In this study, the prognostic model and the recurrence risk model were adopted to evaluate the immune infiltration score of HCC patients with different risk scores. Higher infiltration levels of B memory cells and M0 macrophages and lower infiltration levels of NK cells were found in high prognostic risk patients, and high recurrence risk patients also presented higher infiltration fractions of B memory cells, M0 macrophages and T follicular helper cells and lower infiltration fractions of monocytes, which may influence their response to immunotherapy. Moreover, we found that patients with high prognostic risk expressed higher levels of immune checkpoint genes, such as PD1, B7H3, and CTLA4, in tumor tissues than patients with low prognostic risk. The PD1–PD1 ligand 1 (PDL1) receptor–ligand pair is a crucial immune checkpoint that regulates the tumor microenvironment. This checkpoint controls immune homeostasis by suppressing T cell activation in peripheral cells at the later phase of the cell cycle and is often overexpressed by cancer cells to avoid immune attack.⁵² CTLA4 mainly acts in the early phase of the immune response, participates in the initiation and activation of T cells in lymph nodes, and enhances the immunosuppressive activity of regulatory T cells (Treg cells).⁵³ Studies have shown that anti-CTLA-4 antibodies increase the frequency of tumor-associated antigen-specific cytotoxic T cells in HCC patients and simultaneously exert significant antitumor effects.⁵⁴ B7H3 is mostly expressed on immune cells (such as antigen-presenting cells or macrophages) and tumor cells and has an inhibitory effect on T cells. Its abnormal overexpression helps tumor cells achieve immune escape.³⁴ Our results revealed that high-risk patients tend to have higher expression of immune checkpoints and may be suitable candidates for immunotherapy to improve the immune response and may obtain greater clinical benefit such therapy than low-risk patients.

Considering that chemotherapy is one of the main treatments for patients with advanced HCC and that drug resistance is one of the important factors affecting chemotherapy efficacy,⁵⁵ we also evaluated the sensitivity of patients to common chemotherapy drugs and molecular targeted drugs through the prognostic model. High-risk patients showed increased resistance to chemotherapy drugs and molecular targeted drugs such as methotrexate, erlotinib and vinorelbine compared to low-risk patients, suggesting that chemotherapy efficacy may be more pronounced in low-risk patients than in high-risk patients, and low-risk patients are more suitable for chemotherapy than high-risk patients. In exploring the reasons for differences in the chemotherapy sensitivity of patients, the results of subsequent GSEA showed that the patient’s high prognostic risk is closely related to the signaling pathways involved in the regulation of chemoresistance, such as “DNA replication”,³⁰ “homologous recombination”,⁵⁶ and “mismatch repair”,⁵⁷ which may partly explain the difference in chemoresistance among patients with different prognostic risk.

Inevitably, our research does have some limitations. The difference of cutoff values determined by X-tile between the training set and validation set in the prognosis model may be a limitation that relatively limited the clinical practicability of the prognostic model. In subsequent studies we planned to collect enough clinical samples to further validate the predictive value of this predictive model. The mRNA sequences and corresponding clinical information in this study were mainly obtained from the TCGA and ICGC databases, and the clinical feasibility of the prognostic model in predicting the chemotherapy response and immunotherapy response needs to be prospectively tested in subsequent clinical trials. We will also further explore the role of MYL6B and THOC2 in the progression of HCC in the future.

Conclusion

Prognostic and recurrence models based on two exosome-derived genes (MYL6B and THOC2) were determined to have predictive performance with high specificity and high sensitivity, to generate risk scores that are closely related to traditional clinical variables, and to help identify candidates who are suitable for immunotherapy and chemotherapy. Patients with different risk scores had significant differences in immune infiltration levels, and patients with higher risk scores were more likely to respond to immune checkpoint inhibitor therapy. In addition, low-risk patients may obtain greater clinical benefit from chemotherapy. The diagnostic model constructed with the same two genes could accurately distinguish HCC patients from normal individuals and HCC from hyperplastic nodules. The integration and comprehensive analysis of exosome-derived genes and genomic data of HCC patients provides the possibility to improve the early detection rate of HCC patients and improve their outcomes and provides novel insights into the development of individualized treatment strategies.