Characterization of an endoplasmic reticulum stress‐related signature to evaluate immune features and predict prognosis in glioma

Abstract Endoplasmic reticulum (ER) stress has considerable impact on cell growth, proliferation, metastasis, invasion, angiogenesis and chemoradiotherapy resistance in various cancers. However, the effect of ER stress on the outcomes of glioma patients remains unclear. In this study, we established an ER stress risk model based on The Cancer Genome Atlas (TCGA) glioma data set to reflect immune characteristics and predict the prognosis of glioma patients. Survival analysis indicated that there were significant differences in the overall survival (OS) of glioma patients with different ER stress‐related risk scores. Moreover, the ER stress‐related risk signature, which was markedly associated with the clinicopathological properties of glioma patients, could serve as an independent prognostic indicator. Functional enrichment analysis revealed that the risk model correlated with immune and inflammation responses, as well as biosynthesis and degradation. In addition, the ER stress‐related risk model indicated an immunosuppressive microenvironment. In conclusion, we present an ER stress risk model that is an independent prognostic factor and indicates the general immune characteristics in the glioma microenvironment.


| INTRODUC TI ON
Glioma is the most prevalent primary brain tumour in adults. As well as being highly invasive, it is characterized by diffuse infiltration, fuzzy boundaries and aggressive proliferation. 1 With the development of molecular biological techniques, our understanding of the pathogenesis of glioma has greatly improved, and clinically, important genetic changes have been identified. In the revised classification of central nervous system (CNS) tumours proposed by the World Health Organization (WHO) in 2016, gliomas were classified based on a combination of histological findings and molecular findings, namely isocitrate dehydrogenase (IDH) mutations, 1p19q codeletion, H3 Lys27Met and RELA-fusion. 2 However, although the diagnosis and treatment of gliomas have been greatly improved, the outcome of glioma patients is still unfavourable, with a mortality rate of nearly 80% during the first year of diagnosis. 3 For glioblastomas (GBMs), median survival is only approximately 15 months. 4 Thus, there is an urgent need to improve the diagnosis and treatment of glioma.
As the largest organelle in eukaryotic cells, ER is a membrane structure composed of branched tubules and flat sacs, and is a major site of protein synthesis, processing and transport. 5 However, the protein-processing capacity of the ER is finite. When the protein folding capacity of the ER is exceeded, the cell is considered to be in a state of ER stress state. 6 Many factors can reduce the efficiency of protein folding and lead to ER stress, including oxidative stress, nutrient deprivation, proteotoxicity, hypoxia and metabolic stress, as well as impaired calcium balance. 7,8 It is generally believed that ER stress is triggered by three branches of transmembrane ER sensors: IRE1α, PERK and ATF6. Misfolded proteins are continuously monitored by these three receptors, and when the concentration of misfolded proteins reach a certain level, the sensors trigger the ER stress response. 9 Several studies have confirmed that chronic ER stress is a typical feature of many diseases, including tumours. 10 In tumours, the high metabolism and proliferation of tumour cells lead to ischaemia and hypoxia within the tumour, causing tumour cells to enter a state of continuous ER stress. Signal transduction and regulation induced by the ER stress state promote tumour growth, angiogenesis, immune escape and chemoradiotherapy resistance. 11,12 ER stress can stimulate tumour cells to secrete a heatresistant factor, which can act on leucocytes around tumour cells, thereby changing the local immune characteristics of the tumour and promoting the growth of tumour cells. 13 In addition, macrophages secrete VEGF in response to treatment with the conditioned medium of tumour cells under ER stress, thus enhancing angiogenesis within the tumour microenvironment. 14 In gliomas, ER stress can change the metabolic state of tumour cells to promote tumorigenesis and therapy resistance. 15,16 Therefore, ER stress has become a new target for glioma therapy. 17 General ER stress activators 18 and the selective ATF6, 19 PERK 20 or GRP78 21 activators have been confirmed to regulate multiplication and to facilitate apoptosis of glioma cells. However, some small molecules that interfere with protein processing and folding can inhibit the growth of tumour cells. 22 Therefore, aberrant expression of ER stress-related genes may have prognostic value for glioma patients and can be exploited as potential therapeutic targets.
In this study, we developed an ER stress-related risk model, which can not only accurately predict the outcomes of glioma patients, but also distinguish the immune characteristics of glioma. In addition, we established a nomogram that integrates the prognosis model with clinicopathological factors (age, gender, grade, IDH mutation status, 1p19q codeletion status and MGMT promoter methylation status) and found its performance in estimating 1-, 3-and 5-year survival rates of glioma patients is excellent.

| Data sets and data collection
GeneCards (https://www.genec ards.org/) is a searchable, integrative database that provides comprehensive, user-friendly information on all annotated and predicted human genes. ER stress-related genes were extracted from GeneCards, and genes with a relevance score ≥7 were selected. TCGA RNA-seq transcriptome data and clinical information were obtained from TCGA database (http://cance rgemo me.nih.gov/). The Chinese Glioma Genome Atlas (CGGA) mRNA expression data (mRNAseq_325 and mRNA-array) and corresponding clinicopathological features were collected from the CGGA database (http://www.cgga.org.cn). The GSE16011 data set was procured from the Gene Expression Omnibus (GEO) database (http://www. ncbi.nlm.nih.gov/geo/). Raw data were processed via the R package, affy, and robust multi-array analysis (RMA) was used for background correction and normalization. Biospecimens and clinical data from related research were used as supplements. 23 Patient characteristics, including detailed data for each research object, are shown in Table S1.

| Samples and quantitative real-time PCR analysis
From June 2018 to September 2019, we gathered 12 glioma tissues (4 cases each of grade II, III and IV glioma) from patients at the First Hospital of China Medical University. The clinical characteristics of 12 patients were shown in Table S2. Total RNA extraction and quantitative real-time PCR (qRT-PCR) was performed as previously described. 24 After extraction using TRIzol reagent (Invitrogen/Thermo Fisher Scientific), total RNA was transcribed to first-strand cDNA and then underwent qRT-PCR (SYBR Green Master Mix). Each sample was assayed in triplicate.
After being normalized to GAPDH expression levels, the expression values were log 2 transformed. Primer sequences of target genes are shown in Table S3. This research was approved by the Medical Ethics Committee of the First Affiliated Hospital of China Medical University. All participants provided written informed consent.

| Western blot assay
Western blot assay was performed as previously described. 24 Briefly, the total protein of glioma tissues was extracted using a protein extraction kit (Beyotime). Equivalent amounts of protein (25-50 μg) were then electrophoresed and transferred to polyvinylidene fluoride membranes (0.45 μm; Millipore). After being blocked, membranes were incubated with antibodies against ATF6

| ER stress-related risk signature construction and validation
We first conducted univariate Cox regression and Kaplan-Meier (KM) analysis using the survival R package to identify ER stressrelated genes associated with patients' overall survival (OS) time from TCGA, CGGA (CGGA refers to mRNAseq_325 data unless otherwise specified), CGGA (mRNA-array) and GSE16011 data sets.
Only when the p values of both analysis methods were ≤0.05 were the genes included in the next step. The intersection genes related to OS from the above four data sets were analysed by least absolute shrinkage and selection operator (LASSO) regression, using the glmnet R package in TCGA database, to narrow the range of prognosisrelated genes. Subsequently, the Akaike information criterion (AIC) method of multivariate Cox regression analysis was performed using the survival package to establish an optimal ER stress-related risk signature based on linear integration of the regression coefficient obtained from the multivariate Cox regression analysis and expression level of the selected ER related genes. The risk score was computed as follows: where Exp i is the expression value of the ER stress-related genes and Coef i is the corresponding regression coefficient calculated by multivariate Cox regression analysis. TCGA data were used as the training cohort, and CGGA and GSE16011 data were used for the validation cohorts.

| Functional enrichment analysis
The GSVA package in R was applied to estimate the Gene Ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways that correlated with the risk signature.
The GSVA package scored the GO biological processes and KEGG pathways in each sample, and then by comparing the score differences in different risk groups, we identified the different biological processes enriched in the high-and low-risk groups. The limma package in R was adopted to identify differentially expressed genes and gene sets in different groups. To further verify the GO processes and KEGG pathways related to the signature, GO and KEGG analyses were performed for the differentially expressed genes using the clusterProfiler package in R. All heat maps were drawn using the Pheatmap package in R.

| Evaluation of immune cell fractions and immune subtypes
Myeloid-derived suppressor cells (MDSCs), 25  pheochromocytoma and paraganglioma (PCPG), hepatocellular carcinoma (LIHC) and gliomas, while low-grade gliomas (LGG) consist mostly of subtype C5. 30 The glioma samples were classified into different immune subtypes using the ImmuneSubtypeClassifier package in R. 30

| Independent prognostic role of the risk signature
To identify whether the ER stress-related risk signature depended on other clinicopathological factors (including age, gender, tumour grade, IDH mutation status, 1p19q codeletion status and MGMT promoter methylation status) in predicting patients' OS, univariate and multivariate Cox regression analyses were performed using the survival R package. The results of independent prognostic factor analysis were displayed in the form of forest plots using the forest-

| Development and assessment of the nomogram
The nomogram shows the probability of clinical events through simple graphs of statistical prediction models to form a personalized prediction model. Age, gender, tumour grade, IDH mutation status, 1p19q codeletion status, MGMT promoter methylation status and ER stress-related risk scores were combined to develop a nomogram using the R packages survival and rms. Calibration curves and decision curve analysis (DCA) were adopted to evaluate the accuracy of the nomogram in predicting one-, three-and five-year survival rates of glioma patients. 31,32 For the calibration curve, the higher was the degree of agreement between the predicted curve and the actual curve, the stronger was the predictive ability of the nomogram. For

| Statistical analyses
R software (version 3.6.3) and GraphPad Prism v7.00 (GraphPad Software Inc.) were used as statistical analysis tools in this study. Quantitative data are presented as the mean ± standard error of the mean (SEM) or standard deviation (SD). The Wilcoxon test was applied to compare the statistical differences between the two groups, and the Kruskal-Wallis H test was employed to compare multiple groups. Statistical significance was defined as P < .05. Venn diagrams were drawn using jvenn. 33 The other plots were constructed using R software or GraphPad Prism.

| Sixteen ER stress-related genes were identified to develop risk model
A total of 787 ER stress-related genes with a relevance score of ≥7 were extracted from the GeneCards database to generate prognostic gene signatures (Table S4). These genes were subjected to univariate Cox regression and KM analyses. In the TCGA, CGGA, CGGA (mRNA-array) and GSE16011 data sets, 593, 504, 422 and 330 genes, respectively, were significantly related to OS of glioma patients ( Figure 1A, Tables S5-8). The overlapping genes (190 genes, Table S9) were included in the LASSO regression analysis to avoid overfitting problems in the risk signature ( Figure 1B,C, Table S10).
The AIC method of multivariate Cox regression analysis was applied to the genes returned from the LASSO regression analysis

| The ER stress-related signature is associated with clinicopathological features
We investigated whether the ER stress-related risk signature was correlated with the clinicopathological features of glioma patients.
First, we examined the prognostic effect of the signature in LGG and GBM patients using KM analysis. As displayed in Figure 3A-F, the high-risk group showed reduced OS in terms of both LGG and GBM in the TCGA, CGGA and GSE16011 data sets. Additionally, the risk scores were notably different among stratified patients, with high-risk scores in high-grade glioma, wild-type IDH, 1p19q noncodeletion, mesenchymal subtype and MGMT promoter unmethylated patients ( Figure 3G-N, Figure S4A-D). These results suggest that the ER stress-related risk signature can accurately distinguish different clinicopathologic features of glioma patients.

| Functional annotation of the risk signature
Many studies on ER stress have used the intracellular expression levels of related proteins, such as ATF6, HSPA5, XBP1 and ATF4, as indicators to examine the intensity of ER stress in cells or tissues. [34][35][36] To explore ER stress status in different glioma groups, we measured the expression levels of these markers in the TCGA, CGGA and GSE16011 data sets. Most of them had higher expression levels in the high-risk group, demonstrating that the ER stress in this group was markedly more intense than that in the low-risk group ( Figure 4A-C). To further verify the relationship between the risk score and the ER stress intensity of gliomas, we used qRT-PCR to measure the expression levels of the 16 risk genes in the 12 glioma samples. The samples were scored based on the expression levels of these genes and grouped according to the score.
We then used Western blotting to measure the expression levels of ATF6, EIF2α, p-EIF2α and p-IRE1α in the samples from the high-and low-score groups. The Western blotting results were consistent with those of the above databases analysis (Table S11, Figure S5). These results further confirmed the correlation between the signature and ER stress activation. Next, we used the anti-tumour immunity and promotes tumour cell escape from immunosurveillance. 13,37,38 Here, the enrichment in immune and inflammatory response-related pathways and processes in the high-risk group further demonstrates that the immune environment of the high-risk group is complex.

| High ER stress-related risk score demonstrates an immunosuppressive feature
The process of tumour eradication by the immune system involves the cancer-immunity cycle. 39 We explored the expression characteristics of genes that inhibited this cycle in the TCGA, CGGA and GSE16011 data sets. These genes were acquired from the Tracking Tumor Immunophenotype website (http://biocc.hrbmu.edu.cn/ TIP/index.jsp). 40 In the three cohorts, most genes were highly expressed in the high-risk group ( Figure 5A-C). TGFB1, VEGFA, ARG1, FGL2 and IL10 are secreted immunosuppressive factors in glioma, [41][42][43][44][45][46] while CD95L and CD70 are glioma cell-surface immunosuppressive factors. 47,48 As shown in Figure 5D-O and Figure S8A-I, they were all obviously overexpressed in the highrisk groups. Furthermore, some immune cells, such as MDSCs, Tregs, NK cells, and neutrophils, infiltrate into the tumour microenvironment and exhibit immunosuppressive effects to promote tumour occurrence, progression, and therapy resistance. [25][26][27][28] These immunosuppressive cells were also found to be enriched in the ER stress-related high-risk groups ( Figure 5P-R). To further explore the characteristics of the immune microenvironment in glioma of the high-and low-risk groups, the ImmuneSubtypeClassifier package was used to divide the samples into different immune subtypes in the three cohorts. We found that in both the high-and the lowrisk groups, the main subtypes were C4 and C5, but the high-risk group had considerably more C4 subtypes than the low-risk group ( Figure 5S). The prognosis of the C4 immune subtype in tumours is worse than that of the C5 immune subtype, 30 which was consistent with the prognosis of high-and low-risk glioma patients. This further confirmed the accuracy of the ER stress-related risk signature in predicting immune subtypes and prognosis in gliomas.
In addition to the cancer-immunity cycle inhibitors mentioned above, immune checkpoints can also suppress the immune system's ability to clear tumours. 49 In recent years, immune checkpoints have become potential therapeutic targets for many malignant tumours and play an important role in tumour immunotherapy. 50 Comparison of the expression of immune checkpoints in the high-and low-risk groups showed that most immune checkpoints were up-regulated in the high-risk groups in the TCGA, CGGA and GSE16011 cohorts ( Figure 6A-C). Considering the critical roles played by PD1 and PD-L1 in tumour immunosuppression and immunotherapy, we separately investigated the relationship between their expression levels and the ER stress-related risk score. We found that the expression levels of PD1 and PD-L1 were obviously positively correlated with the risk score ( Figure 6D-G), and their expression levels in the high-risk group were markedly higher than those in the low-risk group ( Figure 6H-K). These data suggest that the ER stress-related risk signature can accurately predict the immune characteristics of glioma.

| Construction and validation of the nomogram
In addition to the ER stress risk score, there are numerous known prognostic factors for glioma, such as age, gender, WHO grade, IDH mutation status, 1p19q codeletion status and MGMT promoter methylation status. Therefore, it was necessary to examine whether the ER stress risk signature could independently predict prognosis. In the TCGA training cohort, univariate Cox analysis revealed that ER stress-related risk score was negatively correlated with the OS of glioma patients. Moreover, age, grade, IDH mutation status, 1p19q codeletion status and MGMT promoter status were also significantly related to OS ( Figure 7A). Subsequent multivariate Cox regression analysis indicated that the ER stress-related risk signature, grade and IDH mutation status were significantly correlated with OS ( Figure 7B). These results were further confirmed in the CGGA and GSE16011 data sets ( Figure S9A-D). These findings indicate that the ER stress-related risk signature constructed using the TCGA data set was an independent prognostic factor for glioma patients.
We established a nomogram to predict 1-, 3-and 5-year OS in the TCGA data set, by integrating ER stress risk signature, age, gender, WHO grade, IDH mutation status, 1p19q codeletion status and MGMT promoter methylation status. In the nomogram, each signature was assigned points according to its risk contribution to OS

| D ISCUSS I ON
Substantial data suggest that a specific intensity of ER stress promotes multiple mechanisms of cancer progression, including cancer cell survival and metastasis, therapeutic resistance, and angiogenesis. 11,12 ER stress can help cancer cells evade immunity and facilitate metastases. 38 Moreover, ER stress in tumour cells can also affect tumour immunity. Cancer cells under ER stress release unknown factors to induce ER stress in macrophages in the microenvironment, inducing them to release pro-inflammatory factors. Simultaneously, these unknown factors also promote bone marrow-derived dendritic cells to release immunosuppressive factors, such as Arginase, to inhibit the antigen presentation ability of CD8 + T cells. 13,37 In breast cancer, ER-stressed tumour cells can up-regulate miR-27a-3p content in exosomes and promote PD-L1 expression in macrophages. 51 Although ER stress in tumour cells has the afore-mentioned effects on immune cells in the tumour microenvironment, the specific mechanism has not yet been fully elucidated. One possible explanation is that ER-stressed tumour cells can induce ER stress in immune cells in the tumour microenvironment. Due to ER stress signalling, such as IRE1α, and PERK, the function of ER-stressed immune cells is inhibited. 13,52,53 Collectively, these findings indicate that a complete understanding of the specific mechanisms of ER stress in the regulation of anti-tumour immunity may greatly improve the efficiency of tumour immunotherapy.
In gliomas, ER stress has an important influence on tumorigenesis, progression and therapy response. 54  However, our study also had some limitations. The expression and prognostic predictive effect of the 16 genes at the protein level require additional assessment. In addition, further research is needed to confirm the specific functional mechanisms of the ER stress-related risk signature in glioma. Although it showed outstanding performance in distinguishing glioma survival differences and immune characteristics, the accuracy of the risk model in discriminating between normal brain tissue and glioma tissue remains to be investigated. In summary, our research provides important resources for elucidation of the specific role of ER stress in glioma. China for technical assistance.

CO N FLI C T O F I NTE R E S T
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.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets used in this study can be downloaded from http://www.