Endoplasmic Reticulum Stress-related Classication for Prognosis Prediction in Hepatocellular Carcinoma

Background: Cancer cells under ER stress are common in hepatocellular carcinoma (HCC) and ER stress is strongly associated with poor prognosis. The aim of this study was to discover credible biomarkers for predicting prognosis of HCC based on ER stress-related genes (ERSRGs). Methods: Univariate Cox regression was performed to calculate the association between ERSRGs and survival outcomes of HCC patients in TCGA. Then LASSO-Cox regression strategy and stepwise Cox regression examination were performed to investigate the quality and establish the prognostic characteristics associated with prognosis. Finally, the model was subsequently validated in two additional independent HCC cohorts. Results: A novel seven-gene prognostic risk model based on ERSRGs was constructed and exhibited superior accuracy in forecasting the survival outcomes and 1-, 2-, 3- year survival rate of HCC patients. qRT-PCR was performed to validate the prognostic risk model in an independent clinical cohort containing 59 HCC patients and the results revealed that this signature had a good prognostic performance. Moreover, we found ER stress could affect the immune microenvironment in HCC and immune checkpoint inhibitors (ICIs) treatment was more effective for patients in high-risk subgroup. In addition, we identied 103 tumor-sensitive drugs in the CellMiner database that may be available for the treatment of HCC patients targeting ER stress and constructed a nomogram combining ER stress-related feature, TNM stage, age and gender. Conclusions: Our seven genetic risk model associated with ER stress can accurately predict survival outcome in HCC patients and facilitate the selection of the best individualized treatment targeting ER stress.


Introduction
The endoplasmic reticulum (ER) is a multifunctional organelle consisting of branching tubes and attened vesicles that is the main site of protein synthesis and transport, lipid biosynthesis and calcium storage [1,2]. However, many factors, including inhibition of protein glycosylation, oxidative stress, nutritional de ciencies, imbalance of calcium homeostasis and hypoxia, could reduce the e ciency of ER in processing protein folding and nally lead to ER stress and unresponsive protein response (UPR) [3,4]. UPR plays a crucial role in regulating cellular adaptation to ER stress by increasing ER content, improving ER protein folding capacity and downgrading misfolded proteins [5,6]. Three transmembrane ER sensors, including activating transcription factor 6 (ATF6), protein kinase R (PKR)-like endoplasmic reticulum kinase (PERK) and inositol-requiring enzyme 1 (IRE1α), have been found to determine the triggering of ER stress and subsequent activation of the UPR [7]. With the increasing recognition of ER stress mechanisms, ER stress dysregulation has been found to play an essential role in various human diseases, including cardiometabolic diseases [8][9][10], diabetes [11,12], chronic kidney disease [13,14], Alzheimer's disease [15,16] and cancers [17][18][19][20]. High levels of ER in hepatocytes can rhythmically activate the hepatic UPR to regulate lipid and glucose metabolism, while ER stress and chronic activation of the UPR due to hepatic viral infection or obesity could cause liver dysfunction and metabolic disorders [21]. ER stress plays a crucial role in the pathogenesis of nonalcoholic fatty liver disease (NAFLD) [22] and is strongly associated with survival and death in HCC patients [23]. Recent studies suggested that ER stress and UPR have emerged as new signaling targets for therapeutic interventions in NAFLD and HCC [3,[24][25][26]]. In the current study, a novel prognostic risk model based on ER stress-related genes (ERSRGs) was constructed, which could be effectively used for prognostic classi cation of HCC patients and utilized as a potential target for individualized immunotherapy.

Public datasets and generation of ERSRGs
This study included mRNA expression data and clinical features of HCC patients from two publicly available datasets including TCGA-LIHC and ICGC (LIRI-JP). ERSRGs were available from previous research [4].

Prognostic risk score model construction and Functional Analysis
The univariable cox relapses were to begin with performed to calculate the a liation between ERSRGs and survival results of HCC patients in two cohorts. LASSO-Cox relapse strategy and stepwise Cox relapse examination were at that point performed to survey the over covering prognosis-related qualities and set up prognostic characteristics. Risk score was at last set up based on the premise of directly combining the equation underneath with the mRNA expression level duplicated the multivariate Cox relapse coe cient (β) demonstrate. Risk score = (βmRNA1 × mRNA1) + (βmRNA2 × mRNA2) +…+ (βmRNAn × mRNAn). We strati ed patients in TCGA dataset into two subgroups due to ideal hazard score edge. The prescient control and autonomy of the prognostic signature in TCGA were evaluated by ROC examination, Kaplan-Meier survival examination and cox relative risks relapse investigation. Gene set enrichment analysis (GSEA) between the two subgroups was performed to distinguish the altogether cautioned GO items with FDR < 0.05.

Clinical specimens and Quantitative real-time PCR (qRT-PCR) analysis
Fresh frozen tumor tissues from previously collected HCC patients were selected as an independent validation cohort [27]. qRT-PCR was used to detect the mRNA levels of genes in the model [28]. After the relative mRNAs expression levels were normalized to β-ACTIN and log 2 transformed, patients were strati ed into two subgroups according to the above formula. Primer sequences are showed in Table S1.

Immune status calculation and immune in ltrates analysis
The immune status of each sample was assessed by applying the ESTIMATE algorithm to the TCGA cohort and calculating immune and stromal scores. Association between risk scores and immune, stromal scores were analyzed by Pearson correlation analysis. To explore impacts of the prognostic model on immunotherapies, we calculated the relationship among risk score and 15 potentially available targeted immune checkpoint genes [29]. Furthermore, in order to assess the potential association between prognostic signature and tumor-in ltrating immune cells (TIICs) in HCC microenvironment, TCGA database was used to measure the abundance ratios of 22 types of TIICs through CIBERSORT [30] (http://cibersort.stanford.edu/). Finally, the predictive ability of signi cantly changed TIICs was assessed by Kaplan-Meier survival analysis.

Genetic alterations and TMB analysis
The mutation and CNA data of 350 HCC patients were downloaded from TCGA to analyze the difference of genetic alterations between the high-and low-risk score subgroups with R package "maftools", and the tumor mutation burden (TMB) of each patient was subsequently assessed.

Drug susceptibility analysis
The association between anticancer drug sensitivity and mRNA molecules in our risk model was directly explored in the CellMiner database [31]. 574 in advanced clinical trials and 216 Food and Drug Administration (FDA)-approved drugs were used for follow-up analyses. Drugs with adjusted P value <0.001 and Pearson correlation coe cient >0.3 as cut-off criteria were considered tumor-sensitive drugs.

Statistical analysis
Categorical data were compared with Pearson chi-square test or Fisher exact test whenever appropriate, and quantitative variables were analyzed using independent-samples t test. ROC curve analysis and Kaplan-Meier survival analysis were performed to assess the prediction performance of survival outcomes with R software (Version 4.0.3). Cox proportional model was performed to analyze relationship between prognostic signature and survival outcomes, together with other clinical features. Clinical characteristics of HCC patients in TCGA, ICGC and clinical validation cohorts were showed in Table S2. Results were considered statistically signi cant when P value <0.05.

Identi cation of overlapped ERSRGs in TCGA and ICGC cohorts
As calculated by univariable cox regressions with an adjust P value <0.05, 330 ERSRGs in TCGA and 225 ERSRGs in ICGC had signi cant prognostic relevance, respectively. Then 111 genes were identi ed as overlapped ERSRGs for further analysis ( Figure 1A). The expression levels of the 111 overlapping ERSRGs in normal and tumor tissues in TCGA-LIHC dataset were shown in Figure S1A. Results of Go and KEGG analysis indicated that these overlapped genes were majorly associated with ER stress-related pathway ( Figure S1B).
3.2 Establishment of an ER stress-related signature in TCGA Overlapped ERSRGs were selected by performing the LASSO-Cox regression model based on the minimum value of λ and 19 genes were screened as shown in Figure 1B. These 19 genes were then placed into a stepwise Cox proportional model and nally a prognostic seven-gene signature was identi ed. Risk score = (0.17243378×MAPT) + (0.09459892×NQO1) -(0. 12755413×PON1) -(0. 14029120×PPARGC1A) + (0. 65200737×NGLY1) + (0. 37755837×CDK5) + (0. 19417209×RNF186). Risk scores for HCC patients were calculated with the above formula, and patients were strati ed into high-or low-risk subgroups with an optimal risk score threshold ( Figure 2C). The association between risk score and clinical characteristics including age, gender, grade, stage, vascular invasion, value of AFP, cirrhosis, HBV infection status and tumor status were evaluated. The results revealed that higher risk scores were linked to advanced TNM stage, later grade, later T stage, HBV infection and recurrence (Supplementary Figure S2). Kaplan-Meier survival analysis revealed that patients with higher risk score were signi cantly relevant to poorer survival outcomes ( Figure 2D). In addition, further strati ed survival analysis was applied for different clinical characteristics, and the results demonstrated that this prognostic model could further differentiated patients with different clinical characteristics including age, vascular invasion, grade, recurrence, TNM stage, gender, HBV infection status and AFP value (Supplementary Figure S3). Finally, ROC analysis revealed that this signature had a good prognostic performance with AUCs at 1-, 2-, 3-year of 0.800, 0.732, 0.745 ( Figure 2E). To explore whether the seven-gene signature could be acted as an independent prognostic model for HCC patients, univariable and multivariate Cox analysis were performed. Univariable Cox regression analysis revealed that this seven -gene signature was statistically associated with survival outcomes for HCC patients (HR=3.088, 95%CI 2.200-4.335, P < 0.001, Figure 2F Left). Then statistically signi cant variables obtained above were entered into multivariate Cox regression analysis, which revealed that this signature could be served as an independent prognostic factor for HCC patients (HR=2.203, 95%CI 1.313-3.694, P < 0.001, Figure 2F Right) after adjusting for other clinical features.

Veri cation of the signature in ICGC cohort
To validate the signature, ICGC dataset was applied as a validation cohort. Risk scores of patients were calculated with the same formula, and patients were strati ed into high-or low-risk subgroups in ICGC cohort ( Figure 3A). Kaplan-Meier survival analysis revealed that patients with higher risk score were prominently relevant to poorer OS rate in ICGC cohort ( Figure 3B). ROC analysis revealed that this signature had a good prognostic performance with AUCs at 1-, 2-, 3-year of 0.723, 0.740, 0.736 in ICGC cohort ( Figure 3C). To explore whether the seven-gene signature could be acted as an independent prognostic model for HCC patients, univariable and multivariate Cox analysis were performed. Univariable Cox regression analysis revealed that this seven -gene signature was statistically associated with survival outcomes for HCC patients (HR=3.054, 95%CI 1.882-4.958, P < 0.001, Figure 3D Left). Then statistically signi cant variables obtained above were entered into multivariate Cox regression analysis, which revealed that this signature could be served as an independent prognostic factor for HCC patients (HR=2.463, 95%CI 1.487-4.078, P < 0.001, Figure 3D Right) after adjusting for other clinical features.

Functional analysis and Immune in ltrates analysis
Five Go items with FDR < 0.05 were enriched in this signature, including adaptive immune response, defense response to virus, immune effector process, response to virus and T cell differentiation ( Figure   3A). According to the results of ESTIMATE algorithm, risk scores was signi cantly associated with stromal scores, but not immune scores ( Figure 3B), and patients in low-risk subgroup had higher stromal scores when compared with patients in high-risk subgroup ( Figure 3C), indicating that this signature was closely related to tumor immune status. Furthermore, based on CIBERSORT algorithm, the differences of 22 types of TIICs in two subgroups in TCGA were assessed by Wilcoxon signed-rank test analysis, and the results demonstrated that HCC patients in low-risk score subgroup had modestly increased ratios of resting mast cells, while patients in high-risk score subgroup had signi cantly elevated ratios of Macrophages.M0 cells ( Figure 3D). Furthermore, Kaplan-Meier survival analysis revealed that only Macrophages.M0 cells were prominently relevant to poor survival outcomes in HCC patients ( Figure 3E).

Genetic alterations and TMB analysis
The results of genetic alterations analysis indicated that the mutation rates of the top 10 most signi cantly mutated genes were remarkably different in the two subgroups ( Figure 4A&4B). Subsequently, the TMB of each patient was assessed. We found that the risk score was closely related with TMB and patients in high-risk scores subgroup had signi cantly increased TMB ( Figure 4C).

Immune checkpoint genes analysis
In the following, the expression levels of 15 potentially targetable immune checkpoint genes were compared between the two subgroups in TCGA database, and results showed that patients in high-risk subgroup had signi cantly increased PD1, PD-L1, CCL2, LAG3, CD276, CTLA4, CXCR4, IL1A, PD-L2, TGFB1, OX40 and CD137 ( Figure 5), indicating that immune checkpoint inhibitors (ICIs) treatment was more effective for patients in high-risk subgroup.

Establishment of a nomogram model in TCGA
To investigate the coe cient prediction e ciency of this signature, a nomogram model was established in TCGA dataset, and the result revealed that the nomogram with a C-index of 0.723 could help us provide a quantitative method for predicting the 1-, 2-, 3-year survival rate accurately ( Figure 6A). The overlap between the forecasted and actual probabilities of 1-, 2-, 3-year survival rate in the calibration curves indicated good agreement ( Figure 6B).

Drug susceptibility analysis
Among the 574 in advanced clinical trials and 216 Food and Drug Administration (FDA)-approved drugs, 103 were considered tumor-sensitive drugs (Table S3) and the top 16 most signi cant tumor-sensitive drugs were shown in Figure 7.

Expression levels of genes in risk model
The mRNA expression of the seven genes in HCC samples were explored in the independent clinical cohort by qRT-PCR ( Figure 8A). Risk scores of patients were calculated with the same formula, and patients were strati ed into high-or low-risk subgroups ( Figure 8A). Kaplan-Meier survival analysis revealed that patients with higher risk score were prominently relevant to poorer OS rate ( Figure 8B). ROC analysis revealed that this signature had a good prognostic performance with AUCs at 1-, 2-, 3-year of 0.846, 0.824, 0.686 ( Figure 8C).

Discussion
There is growing evidence that ER stress-mediated cell proliferation, metabolic conversion and genomic destabilization are important in the development of many chronic liver diseases, including alcoholic liver disease (ALD), hepatic brosis, NAFLD, HBV or HCV hepatitis and HCC [32][33][34]. What's more, some synergistic effects between virus infection, alcohol abuse, NAFLD and ER stress were found in the carcinogenesis of HCC [23]. In addition, ER stress can enhance cancer cell immune evasion and promote recurrence and metastasis by affecting the tumor microenvironment (TME) [35,36]. ER stress-mediated UPR induces autophagy via IRE1α, PERK and ATF6 signaling channels and stimulates vascular endothelial growth factor(VEGF)secretion by macrophages, thereby promoting vasculogenesis in TME [37,38]. Studies to date have shown that ER stress plays a substantial role in regulating tumor cell fate through altered metabolic status and has emerged as a novel signaling target for the treatment of HCC. Inhibition of IRE1α, XBP1s and PERK expression could trigger tumor cell death under ER stressed conditions [24][25][26]. Proteasome inhibitor MLN2238 exacerbates ER stress and promotes cycle stagnation and apoptosis [39]. Sorafenib induces increased ER stress and activates cellular autophagy in HerpG2 cells [40]. Therefore, we hypothesized that aberrant expression of ER stress-related genes may have prognostic value for HCC patients.
In the current study, a novel seven-gene prognostic risk model based on ERSRGs was constructed and exhibited superior accuracy in forecasting the survival outcomes and 1-, 2-, 3-year survival rate of HCC patients in TCGA, ICGC and another cohort of individual clinical samples. More importantly, this feature was an independent risk factor for HCC patients when other clinical factors in the three cohorts were taken into account. In addition, signi cant effects of this feature on the immune microenvironment of HCC and the response to immune checkpoint inhibitors were investigated. Patients in high-risk score subgroup had signi cantly elevated ratios of Macrophages.M0 cells and increased TMB, PD1, PD-L1, CCL2, LAG3, CD276, CTLA4, CXCR4, IL1A, PD-L2, TGFB1, OX40 and CD137, indicating that ER stress could affect the immune microenvironment in HCC and immune checkpoint inhibitors (ICIs) treatment was more effective for patients in high-risk subgroup. In addition, we explored the association between expression of genes in risk model and anticancer drug sensitivity in the CellMiner database, and identi ed 103 tumor-sensitive drugs that may be available for the treatment of HCC patients. Finally, to exploit the full potential of this risk model, a nomogram combining ER stress feature, TNM stage, age and gender was constructed and exhibited superior predictive performance. Therefore, our seven genetic risk models associated with ER stress can accurately predict survival outcome in HCC patients and facilitate the selection of the best individualized treatment.
Of the seven genes in HCC patients, MAPT, NQO1, NGLY1, CDK5 and RNF186 expression were positively correlated with poor prognosis, while PON1 and PPARGC1A expression were negatively correlated with poor prognosis. Pathological aggregation of microtubule-associated protein tau (MAPT), a gene promotes microtubule protein assembly and microtubule stabilization, is the underlying mechanism that triggers endoplasmic reticulum ER stress [41]. Abnormal expression of MAPT could forecast the prognosis of patients in prostate cancer [42], glioma [43] and clear cell renal cell carcinoma [44]. NAD(P)H quinone oxidoreductase 1 (NQO1) plays an important role in cellular defense mechanisms against oxidative stress and ER stress and enhances the invasiveness and apoptosis evasion of hepatocellular carcinoma cells [45,46]. Aberrant expression of NGLY1 affects melanoma survival, tumor growth, and response to therapeutic agents in in vivo and in vitro and targeting NGLY1 is a novel anti-melanoma therapeutic modality [47]. Highly expressed Cyclin-Dependent Kinase 5 (CDK5) can promote the proliferation and oncogenicity of hepatocellular carcinoma cells by regulating the phosphorylation and stability of TPX2 [48], while inhibition of CDK5 enhances the sorafenib response in the treatment of HCC [49]. RING nger protein 186 (RNF186) inhibits the advancement of ER stress-mediated colorectal cancer by inhibiting NF-κB activity [50]. Serum PON1 level is a powerful prognostic factor and can be used to evaluate microvascular invasion in hepatocellular carcinoma [51]. PPARGC1A polymorphism favors the prognosis of patients by affecting the susceptibility to HCC [52]. However, relevant studies on MAPT, RNF186 and NQO1 in HCC are scarce and follow-up experiments are still needed to verify it.
In prior studies, a number of features have been identi ed for predicting survival outcomes in HCC patients [53][54][55][56][57][58]. Compared to these features, this seven-gene risk model has some new features: First, the signature is built based on ERSRG and ER stress is prevalent in HCC, which may provide new ideas for the treatment of HCC by targeting ER stress. Furthermore, the signature has been validated in a small clinical cohort by qRT-PCR analysis to guarantee its clinical relevance. Admittedly, our study has some limitations. The diversity and individual variability of HCC patients may reduce the performance of this feature. In addition, the small sample size limits the validation of the model, and future multicenter randomized controlled studies are needed to evaluate this signature. Finally, the expression and prognostic predictive role of these seven genes at the protein level and their speci c mechanisms in HCC need to be further evaluated in the future by additional in vivo and in vitro experiments.

Conclusions
In summary, we constructed an ERSRG-based risk model in this study, which can effectively classify HCC patients for prognostic prediction and individualized immunotherapy targeting ER stress.

DECLARATION OF COMPETING INTEREST
The authors declare that there is no con ict of interests.