Endoplasmic Reticulum Stress-Related Signature for Predicting Prognosis and Immune Features in Hepatocellular Carcinoma

Hepatocellular carcinoma (HCC) with cancer cells under endoplasmic reticulum (ER) stress has a poor prognosis. This study is aimed at discovering credible biomarkers for predicting the prognosis of HCC based on ER stress-related genes (ERSRGs). We constructed a novel four-ERSRG prognostic risk model, including PON1, AGR2, SSR2, and TMCC1, through a series of bioinformatic approaches, which can accurately predict survival outcomes in HCC patients. Higher risk scores were linked to later grade, recurrence, advanced TNM stage, later T stage, and HBV infection. In addition, 20 fresh frozen tumors and normal tissues from HCC patients were collected and used to validate the genes expressed in the signature by qRT-PCR and immunohistochemical (IHC) assays. Moreover, we found the ER stress-related signature could reflect the infiltration levels of different immune cells in the tumor microenvironment (TME) and forecast the efficacy of immune checkpoint inhibitor (ICI) treatment. Finally, we created a nomogram incorporating this ER stress-related signature. In conclusion, our constructed four-gene risk model associated with ER stress can accurately predict survival outcomes in HCC patients, and the model's risk score is associated with the poor clinical classification.


Introduction
The endoplasmic reticulum (ER) is a multifunctional organelle consisting of branching tubes and flattened vesicles that are 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 deficiencies, imbalance of calcium homeostasis, and hypoxia, could reduce the efficiency of ER in processing protein folding and finally lead to ER stress and unfolded 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]. Chronic activation of the UPR caused by ER stress, viral infection, or hepatic obesity may lead to liver dysfunction and disturbances in lipid and glucose metabolism [21]. ER stress plays a crucial role in the pathogenesis of the 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   Journal of Immunology Research (ERSRGs) was constructed, which could be effectively used for prognostic classification 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 three publicly available datasets including TCGA-LIHC, ICGC (LIRI-JP), and GSE14520. ERSRGs were available from previous research [4]. The general clinical characteristics of HCC patients are shown in Table S1.

Molecular Subtype Identification by Nonnegative Matrix
Factorization (NMF) Algorithm. The 365 HCC samples in TCGA were grouped using the NMF algorithm with the criteria "brunet" and 50 iterations based on ERSRGs. The number of clusters (K) ranged from 2 to 6, with cophenetic, dispersion, and profile being used to determine the ideal number of clusters. Kaplan-Meier survival analysis was also undertaken to see if there were any changes in survival across the NMF subtypes.   5 Journal of Immunology Research investigation. Gene set enrichment analysis (GSEA) between the two subgroups was performed to distinguish the altogether cautioned GO items with FDR < 0:05.

Quantitative Real-Time PCR (qRT-PCR) and
Immunohistochemical (IHC) Analysis. Twenty fresh-frozen tumors and normal tissues from HCC patients who underwent liver tissue resection were collected, and all patients signed informed consent. All pathological data were evaluated and codiagnosed by two independent pathologists. All methods were performed following relevant guidelines and regulations. qRT-PCR was used to detect the mRNA levels of genes in the model [27]. Primer sequences are shown in Table S2. IHC assays were used to explore the protein expression of SSR2 in normal and HCC tissues. SSR2 antibody was obtained from Proteintech (China). Two pathologists independently assessed the results.

Immune Status Calculation and Immune Infiltrate
Analysis. The immune status of each sample was assessed by applying the ESTIMATE algorithm to the TCGA cohort and calculating immune and stromal scores. The association between risk scores and immune and stromal scores was analyzed by Pearson correlation analysis. To explore the impacts of the prognostic model on immunotherapies, we calculated the relationship between risk score and 15 potentially available targeted immune checkpoint genes [28]. Furthermore, to assess the potential association between prognostic signature and tumor-infiltrating immune cells (TIICs) in the HCC microenvironment, the TCGA database was used to measure the abundance ratios of TIICs through CIBERSORT [29] (http://cibersort.stanford.edu/), Timer [30], Quantiseq [31], and xCell [32].
2.6. 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 [33]. 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 coefficient > 0:3 as cut-off criteria were considered tumor-sensitive drugs.

Statistical Analysis.
Categorical data were compared with the Pearson chi-square test or Fisher exact test whenever appropriate, and quantitative variables were analyzed using the 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 the relationship between prognostic signature and survival outcomes, together with other clinical features. Results were considered statistically significant when the P value < 0.05.

Molecular Subtype Identification Based on ERSRGs.
By using the NMF algorithm, the optimal number of clusters was identified as three based on the cophenetic (Figure 1(a)). Then, 365 HCC samples were separated into       9 Journal of Immunology Research three subcategories based on the ERSRGs (Figure 1(b)), and substantial differences were found between patients in the three groupings (Figure 1(c)). When two-by-two comparisons were made between the three groups, only cluster 1 and cluster 3 were significantly different from each other (Figure 1(d)). Interestingly, the stromalscore, immunescore, and ESTIMATEscore were different among the three clusters (Figure 1(e)).

Identification of Overlapped ERSRGs in the Three
Cohorts. As calculated by univariable Cox regressions with an adjusted P value < 0.05, 329 ERSRGs in TCGA, 225 ERSRGs in ICGC, and 152 ERSRGs in GSE14520 cohort had significant prognostic relevance, respectively (Figure 2(a)). Then, 28 overlapped ERSRGs are mainly enriched in ER stress-related biological processes and are used for further analysis (Figure 2(b)).

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 9 genes were screened as shown in Figure 2(c). These genes were then placed into a stepwise Cox proportional model, and finally, a prognostic fourgene signature was identified. Risk score = ð0:07012933 × AGR2Þ + ð0:42892985 × SSR2Þ -ð0:12903062 × PON1Þ + ð 0:40769605 × TMCC1Þ. Only SSR2 showed differential expression in normal and tumor tissues ( Figure S1(a)), although all four genes had prognostic significance ( Figure S1(b)). The expression of the four genes in HCC cell lines was shown in Figure S1(c). Risk scores for HCC patients were calculated with the above formula, and patients were stratified into high-or low-risk subgroups with an optimal risk score threshold (Figure 2(d)). Kaplan-Meier survival analysis revealed that patients with higher risk scores were significantly relevant to poorer survival outcomes (Figure 2(e)), and ROC analysis revealed that this signature had a good prognostic performance with AUCs at 1-, 2-, and 3-year of 0.769, 0.738, and 0.715 (Figure 2(f)). The association between risk score and clinical characteristics including age, gender, grade, stage, and vascular invasion and the value of AFP, cirrhosis, HBV infection status, and tumor status were evaluated. The results revealed that higher risk scores were linked to   Journal of Immunology Research later grade (Figure 3(a)), recurrence (Figure 3(b)), advanced TNM stage (Figure 3(c)), later T stage (Figure 3(d)), and HBV infection (Figure 3(e)). As for age, gender, and cirrhosis status, no significant differences were found (Figures 3(f)-3(h)).  Figure S2). Finally, the results of GO and KEGG functional analysis of the differential genes in the high-and low-risk score groups are shown in Supplementary Figure S3. To explore whether the fourgene signature could be acted as an independent prognostic model for HCC patients, univariable and multivariate Cox analyses were performed, and results 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) in TCGA.

Verification of the Signature in the ICGC and GSE14520
Cohorts. To validate the signature, ICGC and GSE14520 datasets were applied as validation cohorts. Risk scores of patients were calculated with the same formula, and patients were stratified into high-or low-risk subgroups in the ICGC (Figure 4(a)) and GSE14520 cohort (Figure 4(f)). We found lower scores in patients with HCC who were alive (Figure 4(b)) or at earlier TNM stages (Figure 4(c)) in the ICGC dataset. Kaplan-Meier survival analysis revealed that patients with higher risk scores were prominently relevant to poorer OS rates in the ICGC cohort (Figure 4(d)). ROC analysis revealed that this signature had a good prognostic performance with AUCs at 1-, 2-, and 3-year of 0.752, 0.716, and 0.698 in the ICGC cohort (Figure 4(e)). The same results were found in the GSE14520 dataset (Figures 4(g) and 4(h)). Finally, the results of univariable and multivariate Cox analyses revealed that this signature could be served as an independent prognostic factor for HCC patients (ICGC:

Genetic Alterations and TMB Analysis.
The results of genetic alterations analysis indicated that the mutation rates of the top 10 most significantly mutated genes were remarkably different in the two subgroups ( Figure 5(a)). Subsequently, the TMB of each patient was assessed. We found that the risk score was closely related to TMB, and patients in the high-risk scores subgroup had significantly increased TMB ( Figure 5(b)).
3.6. Immune Infiltrates and Immune Checkpoint Gene Analysis. As shown in Figure 6(a), according to the results of the Timer algorithm, HCC patients in the high-risk score subgroup had modestly increased ratios of B cells, CD4 + T cells, neutrophils, and myeloid dendritic cells. The results of the CIBERSORT demonstrated that HCC patients in the low-risk score subgroup had modestly increased ratios of resting CD4 memory cells and monocyte cells, while patients in the high-risk score subgroup had significantly elevated ratios of T helper cells and Treg cells. Moreover, the results of the Quantiseq demonstrated that HCC patients in the low-risk score subgroup had modestly increased ratios of neutrophil and uncharacterized cells, while patients in the high-risk score subgroup had significantly elevated ratios of macrophages M1 and M2 cells (Figure 4(a)). The results of the xCell demonstrated that HCC patients in the lowrisk score subgroup had modestly increased ratios of CD8 central memory cells, endothelial cells, hematopoietic stem cells, macrophage cells, macrophage M2 cells, and NK cells, while patients in the high-risk score subgroup had significantly elevated ratios of CD4 Th1 and Th2 cells. In the 14 Journal of Immunology Research following, we found that patients in the high-risk subgroup had significantly increased PD1, PD-L1, CD276, CTLA4, CXCR4, OX40, and CD137 (Figures 6(b)-6(h)), indicating that immune checkpoint inhibitor (ICIs) treatment was more effective for patients in high-risk subgroup.

Establishment of a Nomogram Model in TCGA.
To investigate the coefficient prediction efficiency of this signature, a nomogram model was established in the TCGA dataset, and the result revealed that the nomogram with a Cindex of 0.723 could help us provide a quantitative method for predicting the 1-, 2-, and 3-year survival rate accurately (Figure 7(a)). The overlap between the forecasted and actual probabilities of 1-, 2-, and 3-year survival rates in the calibration curves indicated good agreement (Figure 7(b)).

Comparison with the Previous Signatures.
Comparing the prediction potential among several genetic signatures can help researchers learn more about their prognostic significance. When compared with other previous signatures [34][35][36][37], as shown in Figure S4, our four-gene signature provided the best survival prediction capacity with fewer genes.
3.9. Drug Susceptibility Analysis. Among the 574 in advanced clinical trials and 216 Food and Drug Administration-(FDA-) approved drugs, 55 were considered tumorsensitive drugs (Table S3), and the top 16 most significant tumor-sensitive drugs were shown in Supplementary Figure S5.

Expression Levels of Genes in the Risk Model.
In keeping with the results of the GEPIA analysis [38], we found that only SSR2 was differentially expressed in tumor and normal samples by qRT-PCR (Figure 8(a)). Furthermore, after calculating the patients' risk scores using the same formula, we found that the scores could distinguish well between normal and tumor tissue (Figure 8(b)). Considering that only SSR2 was differentially expressed, we analyzed it further. We also found that the expression of SSR2 was associated with poor progression in HCC patients (Figures 8(c) and 8(d)). The protein expression of SSR2 was overexpressed in HCC tissues compared to normal tissues (Figure 8(e)). Finally, the results of univariable and multivariate Cox analyses revealed that SSR2 could be served as an independent prognostic factor for HCC patients (HR = 2:1:759, 95% CI 1.191-2.599, P = 0:004) in TCGA.

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 fibrosis, NAFLD, HBV, HCV hepatitis, and HCC [39][40][41]. What is 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 microenvi-ronment (TME) [42,43]. 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 [44,45]. 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 stress conditions [24][25][26]. Proteasome inhibitor MLN2238 exacerbates ER stress and promotes cycle stagnation and apoptosis [46]. Sorafenib induces increased ER stress and activates cellular autophagy in HerpG2 cells [47]. Therefore, we hypothesized that aberrant expression of ER stress-related genes may have prognostic value for HCC patients.
In the current study, a novel four-gene prognostic risk model based on ERSRGs was constructed and exhibited superior accuracy in forecasting the survival outcomes and 1-, 2-, and 3-year survival rate of HCC patients in TCGA, ICGC, and GSE14520 cohorts. 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, significant effects of this feature on the immune microenvironment of HCC and the response to immune checkpoint inhibitors were investigated. Patients in the high-risk score subgroup had significantly increased TMB values, PD1, PD-L1, CD276, CTLA4, CXCR4, OX40, and CD137, indicating that ER stress could affect the immune microenvironment in HCC, and immune checkpoint inhibitor (ICIs) treatment was more effective for patients in high-risk subgroup. In addition, we explored the association between the expression of genes in the risk model and anticancer drug sensitivity in the CellMiner database and identified 55 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 was constructed and exhibited superior predictive performance. Therefore, our four genetic risk models associated with ER stress can accurately predict survival outcomes in HCC patients and facilitate the selection of the bestindividualized treatment.
Of the four genes in HCC patients, SSR2, TMCC1, and AGR2 expressions were positively correlated with poor prognosis, while PON1 expression was negatively correlated with poor prognosis. As an important member of the paraoxonase (PON) family, PON1 plays a very important role in protecting mammalian organisms from oxidative stress [48,49]. Endothelial dysfunction in the body can be brought on by glycosylated PON1-inducing ER stress [50]. PON1 is also implicated in the apoptosis that ER stress causes in tumor cells [51] and serum PON1 level is a powerful prognostic factor and can be used to evaluate microvascular invasion in HCC [52]. An endoplasmic reticulum-based protooncogene called anterior gradient-2 (AGR2) is in charge of maintaining the ER homeostasis. High AGR2 expression improves cell folding and quality control, which are essential for the response to ER stress [53]. AGR2 is also involved in the development and progression of multiple tumors, including pancreatic [54], prostate [55], colorectal [56], breast [57], and endometrial [58] cancers. AGR2 may serve as a potential drug target to improve drug sensitivity during cancer treatment [59]. AGR2 also plays a very important role in the proliferation and metastasis of HCC cells [60]. Signal sequence receptor 2 (SSR2) is a protein involved in the unfolded protein response of the endoplasmic reticulum. Through the translocation of proteins necessary for URP, SSR2 may be particularly implicated in melanoma cells' resistance to ER stress [61]. Our study found that SSR2 is upregulated and is strongly associated with poor prognosis in HCC patients. Upregulated SSR2 could be involved in hepatocarcinogenesis and metastasis through the regulation of epithelial-mesenchymal transition (EMT) [62]. The transmembrane coiled-coil domain (TMCC1) is widely present in vertebrates and lower organisms and belongs to the representative members of the TMCC family. TMCC1 localizes to the rough ER through its C-terminal transmembrane domain and binds to ribosomal proteins through its cytoplasmic region, participating in the regulation of ER stressassociated proteins [63]. Although TMCC1 has not been reported in the literature to be associated with HCC, in our study, we found that HCC patients with high expression of TMCC1 had a poor prognosis. Of course, this requires later collection of multiple HCC cases with complete survival data to validate this result. Overall, the four genes involved in our prognostic model are all related to ER stress, and a more in-depth study of their mechanism of action in HCC may provide new ideas for immunotherapy targeting ER stress in the future.
Compared with the 7-gene prognostic model constructed in our previous study presented as a preprint [64] and other previous signatures [34][35][36][37], this 4-gene prognostic model predicted better with fewer genes. Of course, our research has certain limitations. This feature's effectiveness may be hampered by the diversity and individual heterogeneity of HCC patients. Furthermore, more in vivo and in vitro research is needed to further examine the expression and prognostic predictive relevance of these four genes at the protein level, as well as their unique processes in HCC.

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.

Data Availability
The data used to support the findings of this study are available from the corresponding authors upon request.

Ethical Approval
This study was supported by the Ethics Committees of Zhengzhou University. All methods were performed following the relevant guidelines and regulations.

Consent
Written informed consent was obtained from all patients.

Disclosure
An earlier version of the manuscript has been presented as a preprint according to the following link: https://www .researchsquare.com/article/rs-818967/v1.

Conflicts of Interest
All authors declare no conflicts of interest.

Authors' Contributions
ZGH designed the study, analyzed the data, and wrote the manuscript and the manuscript. SJP downloaded the data and completed the IHC experiments. The final manuscript has been approved by all authors. ZGH and SJP contributed equally to this work and should be considered co-first authors. Figure S1: expression levels and prognostic values of the four genes in normal and HCC tissues explored in GEPIA (a, b) and CCLE database (c). Figure S2: this prognostic model could further differentiate patients with different clinical characteristics. Figure S3: the results of GO and KEGG functional analysis of the differential genes in the high-and lowrisk score groups. Figure S4: comparison of the predictive power of the four-gene prognostic model with other prognostic models for survival. Figure S5: top 16 most important tumor-sensitive drugs. Table S1: clinical characteristics of HCC patients involved in the study. Table S2: the sequences of the qRT-PCR primers used in this study.