Diagnosis and prognosis models for hepatocellular carcinoma patient's management based on tumor mutation burden.

Introduction: The development and prognosis of HCC involve complex molecular mechanisms, which affect the effectiveness of its treatment strategies. Tumor mutational burden (TMB) is related to the efficacy of immunotherapy, but the prognostic role of TMB-related genes in HCC has not yet been determined clearly.
Objectives: In this study, we identified TMB-specific genes with good prognostic value to build diagnostic and prognostic models and provide guidance for the treatment of HCC patients.
Methods: Weighted gene co-expression network analysis (WGCNA) was applied to identify the TMB-specific genes. And LASSO method and Cox regression were used in establishing the prognostic model.
Results: The prognostic model based on SMG5 and MRPL9 showed patients with higher prognostic risk had a remarkedly poorer survival probability than their counterparts with lower prognostic risk in both a TCGA cohort (P < 0.001, HR = 1.93) and an ICGC cohort (P < 0.001, HR = 3.58). In addition, higher infiltrating fractions of memory B cells, M0 macrophages, neutrophils, activated memory CD4 + T cells, follicular helper T cells and regulatory T cells and higher expression of B7H3, CTLA4, PD1, and TIM3 were present in the high-risk group than in the low-risk group (P < 0.05). Patients with high prognostic risk had higher resistance to some chemotherapy and targeted drugs, such as methotrexate, vinblastine and erlotinib, than those with low prognostic risk (P < 0.05). And a diagnostic model considering two genes was able to accurately distinguish patients with HCC from normal samples and those with dysplastic nodules. In addition, knockdown of SMG5 and MRPL9 was determined to significantly inhibit cell proliferation and migration in HCC.
Conclusion: Our study helps to select patients suitable for chemotherapy, targeted drugs and immunotherapy and provide new ideas for developing treatment strategies to improve disease outcomes in HCC patients.

Introduction

As the most common type of primary liver cancer, hepatocellular carcinoma (HCC) is listed as the sixth most common tumor and the fourth leading cause of cancer-related death worldwide, bringing a heavy health burden to society [1]. Since HCC is usually diagnosed at an advanced stage, it is not suitable for surgical treatment, such as radical liver resection or liver transplantation; as such, the prognosis of patients is usually poor, and the probability of recurrence is high [2]. The molecular mechanisms leading to the evolution of HCC tumors are extremely complex, involving multiple genetic abnormalities, including genomic instability, single nucleotide polymorphisms (SNPs), somatic mutations, and dysregulation of signaling pathways related to HCC occurrence and development [3], [4]. Analysis of the complex relationship between genes, proteins and other molecules aids the understanding of the underlying oncogenic molecular mechanisms of HCC; in addition, specific gene expression characteristics obtained from such analysis may help accurately predict the risk of HCC and be translated into novel diagnostic, prognostic and treatment indicators in clinical interventions, providing new options for treatment [5], [6].

In recent years, immunotherapy, including immune checkpoint inhibitors, has attracted wide attention in the field of anticancer therapy as a new treatment option [7]. Many studies have confirmed that multiple types of tumor cells, including HCC cells, carry new antigens induced by gene mutations on the surface that may be recognized by the immune system, making them targets for activated immune cells [8], [9], [10]. However, there is currently no reliable molecular predictive marker that can clarify which HCC patients respond to immunotherapy. Tumor mutation burden (TMB) has been proven to be related with the efficacy of immunotherapy; the higher the TMB is, the better the tumor remission effect and clinical benefits obtained from immunotherapy [11], [12]. TMB is interpreted as the total number of nonsynonymous mutations occurring in each coding region in the tumor genome [13]. Some studies have shown that some types of cancer (such as non-small-cell lung cancer and colorectal cancer) with high TMB values are likely to have a poor prognosis [14], [15], but the prognostic value of the TMB in HCC has not yet been determined clearly.

The rise of high-throughput sequencing technology creates an opportunity to find specific gene expression characteristics that can effectively predict the risk and prognosis in HCC [16]. In this study, we conducted comprehensive analysis on the genomic sequence in HCC and identified genes bearing important mutations contributing to TMB and closely related with the prognosis in HCC patients to establish diagnostic and prognostic models for HCC. We assessed the capability of the prognostic model in predicting the response to chemotherapy, targeted drugs and immune checkpoint inhibitors in HCC patients. Our findings will help develop individualized therapeutic plans for HCC patients and improve their disease outcomes.

Materials and Methods

Determination of differentially expressed hub genes in HCC

The mRNA sequence and corresponding clinical information of HCC patients were obtained from the TCGA database, including 19,676 annotated mRNA sequences and data from 370 HCC tissue samples and 50 normal tissue samples. The limma R package was performed to determine differentially expressed genes (DEGs) with cut-off values of absolute log 2-fold change (FC) > 1 and adjusted P value < 0.05.

Calculation of TMB and definition of high- and low-TMB groups

Somatic mutation profiles were downloaded from the General Genomic Data (GDC) website. Maftools were adopted to evaluate somatic mutations. The formula used to calculate the TMB was as follows: TMBi = 1.0 * NTMi + 1.5 * TMi [17].

TMBi represents the TMB of sample i. NTMi represents the total number of nontruncated mutations, including missense, in-frame deletion, in-frame insertion and nonstop deletions. TMi represents the total number of truncated mutations, including nonsense, frame-shift deletion, frame-shift insertion, and splice-site deletions. Silent mutations were not integrated into the calculation of TMB because they do not cause any downstream changes. According to Foundation Medicine official reports, TMB levels were divided into three groups: low (1–5 mutations/mb), intermediate (6–19 mutations/mb), and high (≥20 mutations/mb) [18]. In this study, we classified samples with TMB < 6 mutations/mb into the low TMB group and samples with TMB ≥ 20 mutations/mb into the high TMB group.

WGCNA for identifying the TMB-specific module

WGCNA was applied in identifying the HCC TMB-specific module using the R package “WGCNA” [19], [20]. The expression profiles of DEGs between HCC patients and normal samples in the TCGA database were used as input for the WGCNA, and TMB was calculated and defined as the phenotype of the sample. A signed scale-free co-expression gene network was ensured through setting the power of β = 10 and scale-free R2 of 0.90 as the soft threshold parameters. A co-expression matrix was established based on the β value and the input gene expression matrix to classify genes with similar expression patterns into the same gene module, thereby generating co-expression modules. The correlation between the module Eigengenes (ME) and TMB was calculated using the Eigengenes function. A heat map was adopted to visualize the correlation between each module and TMB. And the dark green module depicted the highest correlation with TMB among these modules was identified as the TMB-specific module. Then, the correlation matrix between the TMB and the genes in the TMB-specific module was calculated.

Gene set variation analysis (GSVA) of the genes in the TMB-specific module

GSVA was accomplished with the R package “GSVA” [21] to evaluate differentially enriched metabolic pathways in which genes in the TMB-specific module were involved. Pathways used for GSVA were obtained from the Molecular Signatures Database (MSigDB) (http://software.broadinstitute.org/gsea/msigdb).

Determination of genes in the TMB-specific module closely related to the prognosis in HCC

Cox regression and the LASSO method were adopted to assess the prognostic value of the genes in the HCC TMB-specific module. We first used univariate Cox regression to choose genes associated with the prognosis of HCC. Genes with P < 0.05 represents statistical significance. The LASSO method with an L1 penalty was performed to determine the genes with the greatest impact on the prognosis of HCC. LASSO applies an L1 penalty to reduce the regression coefficient and reduce the number of indicators with a final weight of nonzero, thus determining the indicators with the highest contribution [22]. In this study, LASSO to cut down the number of genes was accomplished by the glmnet package in R using 10-fold cross validations and 1000 iterations. Related parameters were set as follows: max_iter = 1000 and cv = 10. After 1000 iterations of LASSO, the higher the nonzero coefficient, the stronger the capability of the corresponding gene in predicting prognosis was. The selected genes were then incorporated in a multivariate Cox regression model and the gene set with the best prognostic value of HCC was determined through forward selection and backward elimination.

Establishment and validation of a prognostic model

A prognostic model was constructed from the gene set identified by the multivariate Cox regression. The prognostic score formula was organized as follows: prognostic index (PI) = (β1 * expression level of SMG5) + (β2 * expression level of MRPL9). The median was set as a threshold to divide HCC patients with survival data into high-risk and low-risk groups. The predictive performance of the prognostic model was examined by K-M curves and ROC curve analysis.

Construction and assessment of a predictive nomogram

Independence analysis of the prognosis-predicting ability of the prognostic model compared to traditional clinical prognostic variables (containing age, AFP, weight, vascular tumor cell invasion, sex, pathologic grade and TNM stage) was carried out using univariate and multivariate Cox regression analyses. The confirmed independent predictive factors were integrated to construct a nomogram and corresponding calibration map via “rms” R software. The calibration curves and consistency index (C-index), which was calculated via a bootstrap method with 1000 resamples, were applied in measuring the consistency between the predicted results of the nomogram and the actual results. ROC curves were generated to compare the specificity and sensitivity of the prediction ability of the nomogram with single independent predictors. Decision curve analysis (DCA) was employed to analyze the clinical benefits achieved by the predictive results of the nomogram compared to those of single independent predictors.

Validation of the prognostic value of gene expression for HCC

Further validation of the prognostic value of SMG5 and MRPL9 was performed in the HCC cohort from the UALCAN database. The Wilcoxon signed-rank test was employed to analyze the expression differences of SMG5 and MRPL9 between tumor tissues and nontumor tissues and examine the effect of the expression of the two genes on the survival time and recurrence time of HCC patients.

Determination of immune cell infiltration

CIBERSORT analysis was adopted to quantitatively convert transcriptome data into the absolute abundance of immune cell and matrix cell types in heterogeneous tissues to evaluate the proportions of 22 human immune cell subsets [23]. Standard annotation files were used to organize gene expression profiles. mRNA data converting into infiltration levels of nontumor cells in tumor microenvironments was realized using the R package “CIBERSORT”.

Cell culture and siRNA treatment

Human hepatocellular carcinoma cell lines (including SK-HEP1 and LM3) were obtained from the American Type Culture Collection (ATCC) (Manassas, VA, USA). The cells were cultured in DMEM containing 10% fetal bovine serum (FBS), 2 mM L-glutamine and 100 U/ml penicillin–streptomycin solution at 37 °C in 5% CO2 in an incubator.

SMG5 and MRPL9 siRNAs were synthesized by GeneChem (Shanghai, China) and transfected into cells using Lipofectamine 2000 (Invitrogen, California, USA) according to the manufacturer’s instructions. The cells were cultured in basic DMEM media for 6–8 h before the media was for DMEM media containing FBS and penicillin–streptomycin.

Western blot analysis

RIPA lysis buffer (Invitrogen) containing PMSF (Bio-Rad, Shanghai, China) was used to collect proteins from SK-HEP1 and LM3 cells. 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) was applied in separating protein samples and a polyvinylidene fluoride membrane (PVDF) membrane (Invitrogen, Carlsbad, USA) was used to transfer the separated protein. The membrane was blocked in 5% skim milk at room temperature for 2 h in a shaker, and then incubated with the primary antibody at 4 °C overnight and subsequent incubation with the secondary antibody for 2 h. The blots were detected and imaged using the iBright FL1500 intelligent imaging system (Invitrogen, Carlsbad, USA). The details of antibodies in this study are listed in the supporting information.

Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

Total RNA was extracted from cultured SK-HEP1 cells and LM3 cells using TRIzol® reagent (Invitrogen, Carlsbad, USA) following the manufacturer's protocol. Then, the RevertAid First Strand cDNA Synthesis kit (Thermo Fisher Scientific, USA) was used for reverse transcribing total RNA to cDNA, and qRT-PCR was performed using Real-time PCR Master Mix (SYBR Green; TOYOBO, Japan). The relative expression (fold change) of the target genes was determined using the 2-ΔΔCT methodology. β-Actin served as the internal control. Each experiment was repeated at least three times. The primers used for qRT-PCR were purchased from RiboBio Co.Ltd. The primers used for qRT-PCR are listed in the supporting information.

Cell proliferation assay

After transfection with SMG5 and MRPL9 siRNA for 48 h, SK-HEP1 and LM3 cells were cultured in 96-well plates (3000 cells/plate in 200 µl DMEM). The proliferative capacity of treated cells was detected at 0, 24, 48 and 72 h. Cell Counting Kit-8 (CCK8) reagent (Yeasen, Shanghai, China) was added to each plate according to the kit instructions, and the OD450 value was analyzed by a microplate spectrophotometer (Thermo Fisher Scientific, MA, USA). Cell proliferation was also measured using 5-ethynyl-2′-deoxyuridine (EdU) agent (Ruibo, Guangzhou, China), and the EdU assay was conducted following the manufacturer's protocol.

Transwell migration assay

SK-HEP1 and LM3 cells were transfected with SMG5 and MRPL9 siRNA for 48 h and cultured in 24-well culture plates with 8 mm pore-containing membrane inserts to measure cell migration capacity. Serum-free DMEM was added to the upper chamber, and DMEM (Gibco, NY, USA) supplemented with 10% fetal bovine serum was added to the lower chamber. Cells were incubated at 37 °C for three days. Cells in the lower chamber (below the membrane) were stained with 0.4% trypan blue (Invitrogen) and counted under a light microscope. Each experiment was repeated three times.

Statistical analyses

R language (version 3.5.2) and GraphPad Prism 7 software (version 7.0) were used to perform statistical analysis and generate figures. Differences in data with a normal distribution between paired groups were compared using Student’s t-tests, and the results are exhibited as the mean ± standard deviation (SD). Differences in data with a nonnormal distribution between paired groups were compared using the Wilcox rank sum test, and the results are expressed as the median (interquartile range). P < 0.05 was considered to represent statistical significance.

Results

Acquisition of differentially expressed hub genes in TMB in HCC

The TMB characteristics for the TCGA HCC cohort and ICGC HCC cohort are presented in Fig. 1A,B. We identified 6,800 HCC-related DEGs from TCGA (Fig. 1C). Through WGCNA, we grouped these genes into modules to aggregate genes with similar traits (Figure S1). Then, the “dark green” module with the highest correlation with TMB (Cor = 0.15, P < 0.01) was identified as the HCC TMB-specific module (Fig. 1D). Fig. 1E indicates that there is a close correlation between the expression characteristics of genes in the TMB-specific module and TMB level. Then, GSVA was performed on the genes in this module. A total of 50 hallmark pathways obtained from the Molecular Signatures Database (MSigDB) (http://software.broadinstitute.org/gsea/msigdb) were incorporated in the analysis (Table S1). We found that the gene set was mainly enriched in the heme metabolism and apoptosis signaling pathways in HCC patients (Fig. 1F).

Fig. 1

Identification of DEGs affecting TMB in HCC. A and B Characteristics of TMB for the TCGA HCC cohort (A) and ICGC HCC cohort (B). C DEGs in the TCGA HCC cohort. D Modules with different traits identified via WGCNA. E Determination of genes in the HCC TMB-specific module. F GSVA for genes in the HCC TMB-specific module.

Comprehensive analysis to determine the association of genes and the prognosis of HCC

Cox regression and the LASSO method were performed on 370 HCC samples with survival data from TCGA to determine the prognostic relevance of the genes in the HCC TMB-specific module. The univariate Cox regression results indicated that the expression of 41 genes was clearly correlated with the prognosis of HCC (P < 0.05) (Table S2). Then, 5 of these genes were further selected through the LASSO method with 10-fold cross validation in 1000 substitution samplings (Figure S2). The subsequent multivariate Cox regression analysis finally determined that 2 genes (SMG5 and MRPL9) could serve as markers for the prediction of prognosis in HCC (P < 0.05) (Table S3).

Construction of the prognostic model and validation of its predictive performance

A prognostic model based on the 2 genes (SMG5 and MRPL9) was constructed. The prediction score of the model was calculated as follows: PI = (0.369 * expression level of SMG5) + (0.404 * expression level of MRPL9). A total of 370 HCC samples with survival data from TCGA were employed as a training set, and 243 HCC samples with survival information from ICGC were adopted as a validation set to assess the predictive performance of the prognostic model. HCC samples were divided into high-risk and low-risk groups using the median as a threshold. K-M curves indicated that patients with high risk had a significantly lower overall survival (OS) probability than their counterparts in both the training set (P < 0.001, HR = 1.93) (Fig. 2A) and validation set (P < 0.001, HR = 3.58) (Fig. 2D). Fig. 2B and 2E show the expression profiles of SMG5 and MRPL9 and the prognostic risk of HCC patients. The AUC values of ROC curves at 0.5, 1, 3, and 5 years reached 0.72, 0.75, 0.66, and 0.67, respectively, in the training set (Fig. 2C) and 0.62, 0.69, 0.72 and 0.71, respectively, in the validation set (Fig. 2F), confirming the high specificity and sensitivity of the prognostic model for predicting the prognosis of HCC. We then further explored the correlation between the prognostic model and TMB in the HCC cohort from TCGA and confirmed that the prognostic score of HCC patients was positively correlated with TMB (Figure S3).

Fig. 2

Survival analysis results, distribution of risk scores, and predictive performance in the training set (A-C) and validation set (D-F). A and D K-M curves showing the prognosis differences between the high-risk and low-risk groups. B and E Distribution of prognostic risk and expression of SMG5 and MRPL9 in patients with HCC. C and F ROC curves for validating the specificity and sensitivity of the prognosis model.

Construction and validation of a nomogram integrating independent predictive factors

370 HCC samples with clinical information from TCGA was adopted to evaluate the independent performance of the prognostic score in predicting prognosis compared with that of clinical characteristics (including age, AFP, weight, vascular tumor cell invasion, sex, pathologic grade and TNM stage) using univariate and multivariate Cox regression analyses. The results determined that age (P < 0.05, HR = 1.735), TNM stage (P < 0.05, HR = 2.218) and the prognostic score (P < 0.05, HR = 1.905) were independent predictive factors of prognosis in HCC (Fig. 3A) (Table S4). A nomogram based on these independent predictive factors was established to quantify survival probability in HCC patients at 1, 3 and 5 years (Fig. 3B). The calibration curves of the nomogram were close to the 45° line (Fig. 3C–E), and the C-index was 0.65, indicating high consistency between the nomogram's prediction and actual results. ROC curves confirmed that the specificity and sensitivity of the nomogram for predicting prognosis were superior to those of any single independent predictive factor (Fig. 3F–H). DCA also demonstrated that prognosis prediction through the nomogram provided the best clinical benefits, suggesting that the nomogram has operational value in clinical practice (Fig. 3I–K).

Fig. 3

The independence of the prognosis model its correlation with clinical pathological features in prognosis prediction. A Forest map showing independent predictive factors for prognosis in HCC. B Nomogram for predicting the survival probability in HCC patients at 1, 3 and 5 years. C-E Calibration charts for validating the predictive accuracy of the 1-year, 3-year, and 5-year survival probabilities of the nomogram. F-H ROC curves comparing the predicted performance of the nomogram and single independent predictive factors. I-K Evaluation of the clinical benefits that the nomogram can achieve.

Determination of differences in immune cell infiltration and immune checkpoint expression in the high-risk and low-risk groups

We determined the level of immune cell infiltration in the TCGA cohort, and the median was used to divide HCC patients with prognostic data into high-risk and low-risk groups. Fig. 4A shows the expression of multiple subtypes of human leukocyte antigen (HLA) in the high-risk and low-risk groups. The distributions of prognosis risk and the 22 infiltrating immune cell fractions in tumor tissues are shown in Fig. 4B. Comparing with the low-risk group, the high-risk group exhibited a higher proportion of infiltrating memory B cells, M0 macrophages, neutrophils, activated memory CD4 + T cells, follicular helper T cells and regulatory T cells within tumor tissues (Fig. 4C–H). We found that the T cell receptor signaling pathway enrichment score was obviously higher in high-risk patients comparing with low-risk patients (Figure S4).

Fig. 4

The correlation between immune cell infiltration and the expression of immune checkpoints and prognostic model components. A HLA subtype expression in high-risk and low-risk patients. B Distribution of prognostic risk and immune cell infiltration within tumor tissues in patients with HCC. C-H Violin charts revealing the relationship between the fraction of immune cells and prognostic risk score (C memory B cells; D M0 macrophages; E neutrophils; F activated memory CD4 T cells; G follicular helper T cells; H regulatory T cells). I Distribution of prognostic risk score and immune checkpoint expression in patients with HCC. J The association between prognosis risk and immune checkpoints. K-N Column charts showing the expression of immune checkpoints in high-risk and low-risk patients (K PD1; L B7H3; M CTLA4; N TIM3).

As an innovative therapy, immunotherapy based on immune checkpoint inhibitors is gradually being used in the treatment of various types of advanced cancer. Tumors usually upregulate immune checkpoints to avoid detection and killing by the immune system, while immune checkpoint inhibitors can rejuvenate host immune cells and target the adaptive immune system to tumor cells [24]. Fig. 4I shows expression profiles of immune checkpoint molecules and prognostic risk for 370 HCC samples from TCGA. The correlation between prognosis risk score and the expression of immune checkpoints is visualized in Fig. 4J. We found that a positive correlation was exhibited between prognosis risk score and immune checkpoints such as B7H3 and CTLA4 (Table S5). Then, we further revealed that higher relative expression of PD1, B7H3, CTLA4 and TIM3 was presented in patients with higher prognostic risk comparing with patients with lower prognostic risk (Fig. 4K–N), suggesting that patients with high risk may obtain more clinical benefits from immunotherapy.

Prediction of therapeutic response to chemotherapy drugs and targeted drugs according to prognostic risk in patients with HCC

Chemotherapy and targeted drugs are common treatments for patients with HCC, especially suitable for patients with HCC who cannot undergo surgery to reduce tumor volume and prolong their survival as is postoperative adjuvant chemotherapy to inhibit the recurrence and progression of tumor [25], but drug resistance is a critical issue influencing the effectiveness of anticancer drugs. We evaluated the therapeutic response of HCC patients to 266 anticancer drugs on the Genomics of Drug Sensitivity in Cancer (GDSC) website with reference to half the maximum inhibitory concentration (IC50). The responses of HCC patients to traditional chemotherapeutic drugs and molecularly targeted chemotherapeutic drugs such as methotrexate, vinblastine, erlotinib and gefitinib are presented in Fig. 5A–O. We found that more resistance to chemotherapy was exhibited in high-risk patients were than in low-risk patients (P < 0.05). The top 5 positively regulated pathways are presented in Fig. 5P, and the top 5 negatively regulated pathways are presented in Fig. 5Q. Through GSEA, we found that the prognosis model was related to signaling pathways involved in physiological processes including the cell cycle, mismatch repair and spliceosome signaling, and biological processes including the cell cycle [26], mismatch repair [27] and spliceosome signalling [28] were shown to have a regulatory effect on cancer cell chemoresistance in previous studies, and differences in these pathways in patients may lead to different prognosis risks and different resistance to chemotherapy.

Fig. 5

Response of HCC patients to chemotherapy drugs. A-O Differences in response to chemotherapy drugs between high-risk and low-risk patients. P Top 5 signaling pathways positively regulated by the TMB-specific prognostic model. Q Top 5 signaling pathways negatively regulated by the TMB-specific prognostic model.

Construction and validation of diagnostic models to accurately distinguish patients with HCC from normal subjects and those with dysplastic nodules

For improving the probability of early detection in HCC, a stepwise logistic regression method was performed to build a diagnostic model. The diagnostic formula is determined as follows: logit (P-HCC) = -7.819 + (0.160 × SMG5 expression level) + (1.014 × MRPL9 expression level). 50 normal samples and 50 paired HCC samples from TCGA were used as a training set, and 202 normal samples and 210 paired tumor samples from ICGC were adopted as a validation set to determine the performance of the diagnostic model. The results suggested that the diagnostic model achieved 88.00% specificity and 84.00% sensitivity in the training set (Fig. 6A) and 92.57% specificity and 85.71% sensitivity in the validation set (Fig. 6C). The AUCs of the training set and validation set were 0.928 (Fig. 6B) and 0.950 (Fig. 6D), respectively, confirming that the diagnostic model made a accurate distinction between HCC and normal subjects. Fig. 6C and 6G show the distribution of SMG5 and MRPL9 expression in HCC samples, as well as the predicted disease status and actual status. High expression of SMG5 and MRPL9 was exhibited in HCC patients. Next, we explored the relationship between SMG5 and MRPL9 and revealed that there is a positive correlation between the expression levels of the two genes (P < 0.001) (Fig. 6F and H).

Fig. 6

A diagnostic model for differentiating HCC from normal (A-H) and dysplastic nodule (I-P) samples. A and C Confusion matrix of the binary results in the diagnostic model for distinguishing HCC and normal subjects. B and D ROC curves confirming the predictive accuracy of the diagnostic model. E and G Expression levels of SMG5 and MRPL9 in patients with HCC: distribution of the predicted results and actual results. F and H The correlation between the expression of SMG5 and MRPL9. I and K The specificity and sensitivity of the diagnostic model for distinguishing HCC lesions from dysplastic nodules. J and L ROC curves validating the predictive performance of the diagnostic model. M and O Expression characteristics of SMG5 and MRPL9 in patients with HCC: distribution of the predicted results and actual results. N and P Positive correlation of SMG5 and MRPL9 expression levels.

Currently, early diagnosis of HCC mainly relies on radiology examination combined with pathological diagnosis, but small nodules<2 cm have an increased likelihood of being missed because they are difficult to characterize via pathological or radiology examination [29]. We built a diagnostic model integrating SMG5 and MRPL9 in a training set (GSE6467) with 35 HCC samples and 17 dysplastic nodule samples to detect whether HCC and dysplastic nodules could be accurately distinguished. The diagnostic formula was identified as follows: logit (P-HCC) = -51.308 +(2.967 × SMG5 expression level) + (3.966 × MRPL9 expression level). GSE98620, which contains data from 49 HCC samples and 24 dysplastic nodule samples, was adopted as a validation set to further validate the diagnostic capability of the diagnostic model. The results determined that the diagnostic model could correctly distinguish HCC from dysplastic nodules, with specificity of 82.35% and sensitivity of 91.43% in the training set (Fig. 6I) and specificity of 54.17% and sensitivity of 77.55% in the validation set (Fig. 6K). The AUC reached 0.948 in the training set (Fig. 6J) and 0.776 in the validation set (Fig. 6L), indicating high consistency between predicted diagnostic status and actual status. The expression of SMG5 and MRPL9 in patients, predicted disease status and actual status are shown in Fig. 6M and P. Consistent with the above finding, the expression levels of SMG5 and MRPL9 were positively correlated (P < 0.001) (Fig. 6N and 6P).

Determination of the correlation of SMG5 and MRPL9 with prognosis and recurrence of HCC

We explored the expression characteristics of SMG5 and MRPL9 through the online database Oncomine and the data visualization platform GEPIA for external validation. We found that in the Oncomine database, higher expression levels of SMG5 and MRPL9 was exhibited in tumor tissues compared to normal tissues (Fig. 7A,B), and the results of analysis in GEPIA were similar (Fig. 7C,D). A log-rank test and K-M curve analysis were conducted to further assess the predictive value of SMG5 and MRPL9 expression for prognosis and recurrence in HCC. The results determined that patients with HCC and higher expression levels of SMG5 and MRPL9 had remarkedly shorter survival time (Fig. 7E,F) and lower relapse-free survival (RFS) time (Fig. 7G,H) than those with lower expression levels of SMG5 and MRPL9. The above findings confirmed that the prognostic model based on SMG5 and MRPL9 is reasonable and holds value in clinical practice. Besides, we used the data of HCC patients in the TCGA-LIHC cohort to evaluate the correlation between SMG5 and MRPL9 and immune cells in the TIMBER database. Fig. 7I,J shows that the expression of these two genes conducted correlation with the infiltration level of immune cells, such as B cells, CD4 + T cells, macrophages, and neutrophils.

Fig. 7

The value of SMG5 and MRPL9 in predicting the prognosis and recurrence of HCC. A-B The expression of SMG5 (A) and MRPL9 (B) as assessed in Oncomine. C-D The expression of SMG5 (C) and MRPL9 (D) as assessed in GEPIA. E-F Survival analysis for SMG5 (E) and MRPL9 (F). G-H Recurrence analysis for SMG5 (G) and MRPL9 (H). I-J The correlation of the expression of SMG5 (I) and MRPL9 (J) with the infiltration of different immune cells.

Measurement of the oncogenic effect of SMG5 and MRPL9 in HCC cells

To further explore the role of SMG5 and MRPL9 on the proliferation and migration of HCC cells, the expression of SMG5 and MRPL9 in human LO2 hepatocytes and different HCC cell lines was detected using qPCR. The results showed that SMG5 and MRPL9 exhibited the highest expression in SK-HEP1 and LM3 cells (Figure S5A and S5B), so we used SKHEP1 and LM3 to perform follow-up experiments. We knocked down the levels of SMG5 and MRPL9 in HCC cell lines (including SK-HEP1 and LM3 cells). Western blotting showed that the administration of SMG5 siRNA and MRPL9 siRNA significantly inhibited the expression of SMG5 and MRPL9 in cells (Fig. 8A–D). The CCK assay indicated that both SMG5 inhibition (Fig. 8E and G) and MRPL9 inhibition (Fig. 8F and 8H) significantly inhibited the proliferation of HCC cells. The EdU assay also determined that SMG5 and MRPL9 inhibition showed a significant inhibitory effect on the proliferation of SK-HEP1 cells and LM3 cells. (Fig. 8I,J). In addition, the results of transwell analysis revealed that SMG5 and MRPL9 inhibition obviously inhibited the migration of SK-HEP1 cells and LM3 cells (Fig. 8K–L). The quantitative statistics are presented in Fig. 8M–P. The findings above confirmed that the expression of SMG5 and MRPL9 had a close association with the proliferation and migration of HCC cells and further confirmed the value of SMG5 and MRPL9 in predicting the prognosis of HCC.

Fig. 8

The effect of SMG5 and MRPL9 on the progression of HCC. A–D Western blot analysis confirmed that the expression of SMG5 and MRPL9 was inhibited by SMG5 and MRPL9 siRNA administration. E-H CCK8 assay indicated that SMG5 and MRPL9 inhibition significantly suppressed the proliferation of SK-HEP1 cells (E-F) and LM3 cells (G-H). I-J EdU assay revealed that SMG5 and MRPL9 inhibition showed a significant inhibitory effect on the proliferation of SK-HEP1 cells and LM3 cells, respectively. K-L Transwell migration assays confirmed that SMG5 and MRPL9 inhibition obviously inhibited the migration of SK-HEP1 cells and LM3 cells. M−P Quantitative statistical results of the effects of SMG5 and MRPL9 expression on the migration of SK-HEP1 cells (M−N) and LM3 cells (O-P). Data are shown as the mean ± SD of at least three independent experiments. *P < 0.05, **P < 0.01, ***P < 0.001.

Discussion

HCC is a major health problem worldwide with its increasing morbidity and mortality [30]. HCC exhibits obvious molecular heterogeneity, including various somatic genome mutations in addition to epigenetic modifications [2]. However, to date, no proposed molecular markers can accurately predict the prognosis or recurrence of HCC. TMB is a characteristic representing the total number of somatic coding mutations in tumors [31]. There is increasing evidence that TMB has value as a potential biomarker in non-small-cell lung cancer, rectal cancer and other cancer types [18], [32], [33]. Therefore, identifying and utilizing genes that carry mutations that provide important contributions to TMB may help optimize the early diagnosis of HCC and predict the prognostic risk of patients.

The increased ease of use of high-throughput sequencing technology opens up opportunities for us to better understand the pathological process of HCC and identify key genes involved in the carcinogenesis of HCC [34]. In this study, we focused our research on those mutated genes mostly intricately involved in TMB and tried to clarify the prognostic and diagnostic value of these key genes in HCC through a series of comprehensive analyses to accurately stratify the disease risk of patients with HCC and develop individualized treatments. We used integrated analysis, Cox regression, and the LASSO method to determine that SMG5 and MRPL9 have strong prognostic value in HCC. Knockdown of SMG5 and MRPL9 was determined to significantly inhibit cell proliferation and migration in HCC, suggesting that the expression of these two genes has a close association with the progression of HCC. A prognostic model based on the two genes could independently predict the prognosis of HCC patients with high specificity and sensitivity. SMG5 is a regulator of the nonsense-mediated messenger RNA (mRNA) decay (NMD) pathway in cells, which is a cell quality control and post-transcriptional gene regulation mechanism. SMG5 is thought to have important regulatory effects on the viability of most multicellular organisms [35]. MRPL9 is a nuclear gene encoding protein components in the ribosome in mitochondria, which is essential for mitochondrial function [36].

HCC has an insidious onset and progresses rapidly, and it is usually diagnosed in advanced stages (BCLC stage B and C). Chemotherapy is one of the most important treatments for patients with advanced HCC, especially for those who cannot undergo surgical resection, have a poor response to TACE, and have extrahepatic metastases or vascular invasion [37]. However, tumor cells may develop intrinsic or acquired resistance during the process of chemotherapy, resulting in poor chemotherapy effects and affecting the prognosis of patients [38]. Changes in gene expression profiles associated with somatic mutations in HCC play a significant role in regulating the sensitivity of tumor cells to chemotherapy [39]. In addition to drug-resistant regulatory factors and carcinogenic drivers, somatic mutations can also produce uniquely altered proteins, including new antigens to which specific antitumor immunity can be raised. These mutation-derived antigens can be recognized by the immune system and targeted by immunotherapy [40], [41]. The more somatic mutations a tumor has, the more neoantigens it may form. Clinical studies have confirmed a significant correlation between the production of neoantigens and immune-mediated clinical responses [42]. TMB is thus also regarded as a potential marker for immunotherapy and can provide a useful estimate of tumor neoantigen load [43]. Through comprehensive analysis of genes that play an important role in TMB, potential markers may be found to identify patients sensitive to chemotherapy and suitable for immunotherapy.

In our study, we used a TMB-based prognostic model to stratify HCC risk and assess immune cell infiltration and immune checkpoint expression in patients with different risks and test patient response to different chemotherapy drugs; by so doing, we further explored the effect of the expression of SMG5 and MRPL9 on the response of HCC patients to immunosuppressive therapy and chemotherapy. We found that compared with low-risk patients, high-risk patients showed higher levels of immune cell infiltration and higher expression levels of immune checkpoints in tumor tissues, and high-risk patients were more resistant to chemotherapy drugs. Subsequent GSEA suggested that the genes involved in the prognostic model were related to the signaling pathways related to physiological processes such as the cell cycle [26], mismatch repair [27] and the spliceosome [28], which have been found to have a regulatory effect on the chemoresistance of cells in previous studies. It was reported that SMG5 is involved in regulating the expression of a series of genes with a wide range of cellular activity functions, including cell cycle progression [44], which may partly explain the differences in chemoresistance among patients with different prognostic risks. The findings above indicate that in clinical management, high-risk patients may be the beneficiaries of targeted immunotherapy, while low-risk patients can benefit more from chemotherapy.

The early diagnosis of HCC mainly relies on radiological examination and pathological examination. The development of imaging technology has improved the recognition rate of HCC, enabling small HCC tumors (<2 cm) to be diagnosed at an early stage for further histological identification [45], [46]. However, increasing evidence suggests that small nodules are difficult to characterize. There is no clear boundary in morphology and histology between high-grade proliferative nodules and well-differentiated HCC lesions, and it is difficult to accurately distinguish them [47], [48]. Therefore, it is urgent to find molecular markers to objectively and accurately identify HCC. Two diagnostic models including SMG5 and MRPL9 expression were confirmed to accurately distinguish HCC from normal and dysplastic nodule samples (<2 cm) in this study, suggesting that the model may be beneficial to improve the early detection rate of HCC and is conducive to early clinical intervention in patients with HCC, thus improving the prognosis and likelihood of relapse of patients. Our findings lay the foundation for using SMG5 and MRPL9 as biomarkers for the early diagnosis of HCC patients.

Inevitably, there is some limitations in our study. First, when using the LASSO method to filter genes in the early analysis, the limitations of the LASSO method itself led to the potential for missing some genes with similarly important contributions when the regression coefficients were adjusted. The clinical characteristics considered in the independence analysis of the prognostic model and the construction of the nomogram are traditionally recognized as important factors affecting the carcinogenesis of HCC, while some clinical variables with similar contributions, such as geographical factors, dietary patterns, and exposure to endogenous or exogenous chemicals, were not included in the study due to insufficient sample information, which may affect the results. When identifying HCC lesions and nodules, the small sample size also affected the evaluation of the performance of the diagnostic model. Besides, we will further confirm the performance of the prognostic model for evaluating immunotherapy and chemotherapy response in follow-up clinical trials.

Conclusion

The prognostic model based on SMG5 and MRPL9 expression was proven to have accurate prediction capability, and it helped to identify patients suitable for immunotherapy. Patients with higher prognostic risks may obtain better clinical benefits through immunotherapy and achieve greater tumor remission. Prognostic risk was confirmed to have a close correlation with chemotherapeutic response. Low-risk patients had a significant survival advantage compared to high-risk patients receiving chemotherapy. Diagnostic models constructed using SMG5 and MRPL9 can make an accurate distinction between HCC and normal and dysplastic nodule samples. The comprehensive analysis of multidimensional genomic data opens up new ideas for determining patient treatment strategies and improving disease outcomes.

Data and materials availability

The data and materials applied in supporting the findings in this study are available from the corresponding author upon request.

Author Contributions

Jiansong Ji and Risheng Yu conceived and designed the experiments. Bufu Tang and Jinyu Zhu performed the experiments. Zhongwei Zhao, Chenying Lu, Siyu Liu and Shiji Fang analyzed the data. Liyun Zheng, Nannan Zhang, Minjiang Chen and Min Xu contributed analysis tools. Jinyu Zhu and Bufu Tang wrote the paper. Jiansong Ji and Risheng Yu edited the paper. All authors read and approved the final manuscript.

Compliance with Ethics Requirements

This article does not contain any studies with human or animal subjects.

Declaration of Competing Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.