Identification of the role of endoplasmic reticulum stress genes in endometrial cancer and their association with tumor immunity
BMC Medical Genomics volume 16, Article number: 261 (2023)
Endometrial cancer (EC) is one of the worldwide gynecological malignancies. Endoplasmic reticulum (ER) stress is the cellular homeostasis disturbance that participates in cancer progression. However, the mechanisms of ER Stress on EC have not been fully elucidated.
The ER Stress-related genes were obtained from Gene Set Enrichment Analysis (GSEA) and GeneCards, and the RNA-seq and clinical data were downloaded from The Cancer Genome Atlas (TCGA). The risk signature was constructed by the Cox regression and the least absolute shrinkage and selection operator (LASSO) analysis. The significance of the risk signature and clinical factors were tested by time-dependent receiver operating characteristic (ROC) curves, and the selected were to build a nomogram. The immunity correlation was particularly analyzed, including the related immune cells, pathways, and immune checkpoints. Functional enrichment, potential chemotherapies, and in vitro validation were also conducted.
An ER Stress-based risk signature, consisting of TRIB3, CREB3L3, XBP1, and PPP1R15A was established. Patients were randomly divided into training and testing groups with 1:1 ratio for subsequent calculation and validation. Based on risk scores, high- and low-risk subgroups were classified, and low-risk subgroup demonstrated better prognosis. The Area Under Curve (AUC) demonstrated a reliable predictive capability of the risk signature. The majority of significantly different immune cells and pathways were enriched more in low-risk subgroup. Similarly, several typical immune checkpoints, expressed higher in low-risk subgroup. Patients of the two subgroups responded differently to chemotherapies.
We established an ER Stress-based risk signature that could effectively predict EC patients’ prognosis and their immune correlation.
EC is a worldwide gynecological malignancy, ranking as the fourth most common cancer and the sixth top cause of cancerous death in females. Although the median diagnostic age is already more than 61 years, EC still tends to occur in the younger population in recent decades . The most significant risk factor for EC is exposure to estrogen, whether endogenous or exogenous, or even obesity. Other correlated factors include delayed menopause, early-onset menarche, nulliparity, diabetes, Lynch Syndrome, and Cowden Syndrome . Despite many obvious achievements in EC treatment, such as surgery, radiotherapy, cytotoxic chemotherapy, and hormonal therapy , the morbidity and mortality rates of EC patients remain quite high globally. So, the action of identifying molecular mechanisms and then searching for novel treatment targets is still urgent for efficiently improving or even remarkably advancing the EC treatment.
As demonstrated, the tumoral abnormalities will cause the microenvironment mess, such as ischemia, hypoxia, oxidative stress, nutrient imbalance, and DNA damage, so as to disturb the homeostasis of endoplasmic reticulum, following misfolded or unfolded protein accumulation, thus leading to an intracellular state named ER Stress . In turn, the continuous activation of the ER Stress will alter the tumorigenesis through both the transcriptional and translational pathways . Furthermore, ER Stress can reprogramme the function of immune cells and thus apparently temper anti-cancer immunity . ER Stress induces tumor cells to release factors that change the leucocytes around tumor cells, thus promoting tumor growth . ER-stressed tumor cells secrete PD-L1 that drives M2 macrophage polarization, facilitating the escape from immune surveillance . On the other hand, ER-stressed M1-like macrophages secrete DAMPs that lead to tumor cell death .
So far, continuously accumulated studies have presented that ER Stress should have played a critical role in endometrial carcinogenesis . Protein and mRNA levels of several ER Stress-induced apoptosis indicators, including GRP78 and ATF6, are apparently increased in endometrial cancer immunohistochemical samples . And GRP78 plasma membrane localization is elevated in endometrial cancer . Therefore, cancer therapy focused on ER Stress is promising for deeper exploration.
In this paper, we constructed a bioinformatic method-based risk signature to explore the relevance and effect of ER Stress on EC progression. First, by combining ER Stress-related gene expression with clinical data, significantly differentially expressed genes between normal endometrial tissues and EC tissues were selected. After Cox analyses and LASSO regression analysis, 4 prognostic genes were filtered out. Based on these 4 genes that were identified as screening biomarkers of EC, a risk signature (risk score) was established by secondary clustering, which helped to discover molecular characteristics and immune infiltration specific to ER Stress of EC patients.
According to the survival prediction, EC patients with higher risk scores represented worse prognoses which indicated the clinical implication of our risk signature. For further analyses, results showed patients with higher risk scores were with less immune infiltration to evade the immune surveillance, as well as a lower rate of immune checkpoints, demonstrating that risk signature was capable to forecast the EC patients’ immunotherapies response. Furthermore, special potential chemotherapies were identified based on risk scores. 4 risk signature genes were then experimentally verified in EC cell lines and EC tissues. One challenge our risk signature face is that despite the strong correlation between immune and our risk signature, our risk signature was not significantly correlated with Tumor Mutation Burden (TMB), which is a potential biomarker in tumor immunotherapies selection [13, 14]. Overall, our study provided an effective prognostic signature based on ER Stress-related genes and helped personalized therapy for EC patients.
Materials and methods
Normal and tumor datasets
The transcriptome RNA-seq data of the EC samples (554 tumor datasets) and the normal samples (23 normal datasets) were obtained from TCGA database (https://portal.gdc.cancer.gov/repository) on 1st January 2023. The clinical features of EC patients were also obtained from the TCGA.
Identification of ER stress-related differentially expressed genes
304 ER Stress-related genes (ERGs) were obtained from GSEA (http://www.gsea-msigdb.org/gsea/index.jsp), while 818 ERGs from GeneCards (http://www.genecards.org/) and 178 ERGS were finally selected after finishing intersection. The differentially expressed genes (DEGs) were selected based on the mRNA expression difference between the EC patients and the control via “limma” R package , and then following |log2fold change (FC)|>1 and false discovery rate (FDR) < 0.05. The Protein-Protein Interaction (PPI) networking was conducted by utilizing search tools for the retrieval of interacting genes (STRING) (http://string-db.org/) .
Establishment and validation of the risk signature based on differentially expressed genes
To improve the predictive value of the DEGs in the EC patients selected, we already combined the DEGs with survival data from TCGA. The univariate Cox (uni-Cox) regression analysis was performed to filter prognostic ERGs from the DEGs with P < 0.05. 6 ERGs were identified for subsequent analysis. LASSO and the multivariate Cox (multi-Cox) regression analyses  were applied to establish a prognostic signature with “glmnet” R package . 4 ERGs were eventually selected to calculate the risk score. The formula for performing risk score [19, 20] is presented below:
where the expr is an expression of the ERGs, and the coef is the corresponding regression coefficient calculated by multi-variant Cox regression analysis.
Next, we summarized 523 samples and then randomly divided them into a training group and a testing group with a 1:1 ratio, and the extra portion was placed in the training group. The training group was mainly used to establish a risk signature, while the testing one was then utilized to validate the risk signature. The samples from both the training and testing groups were respectively separated into a high-risk subgroup and a low-risk subgroup according to the median risk score. The Kaplan-Meier analysis was performed on both the functioned groups above and ROC curves of the overall survival (OS) were employed to assess both the sensitivity and effectiveness of the risk signature. The patients with missing OS value or its value < 30 days were eliminated to reduce the statistical bias during comparison.
Independent prognostic factor analysis
We tested whether the clinical characteristics (age, stage, grade) and the risk score could be used as independent prognostic factors through uni-Cox and multi-Cox regression analyses. By using “timeROC” and “survminer” R packages, both the time-dependent ROC curves and AUC were respectively constructed for comparing the efficiency of different predictive factors and the time-dependent survival rate.
Function prediction analysis
To identify the possible functions of the ERGs chosen, we conducted both the gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) [21,22,23] enrichment analyses based on the DEGs between the high- and low-risk subgroups by applying “clusterProfiler” R package . The KEGG analysis was employed to find out the significantly enriched pathways, and then GO analysis consists of three enrichment parts: cellular component (CC), molecular function (MF), and biological process (BP), respectively.
Nomogram and calibration analysis
The essential issues including the nomogram, the integrating risk score, age, stage, and grade along with a consistency index, were created to illustrate the prediction efficiency of 1-, 2-, and 3-year OS using the “rms” R package. The calibration curves were used to visualize the results and then evaluate the prediction consistency of the nomogram. The diagonal line (45°) can be recognized according to the best prediction value.
Tumor immune analysis
We integrated the risk scores and the immune factors for predicting the relationship between the tumor immune and the ERGs. The immune cells and the immune pathway infiltrating scores were calculated according to the TCGA data via “gsva” R package . The expression and survival differences of the high- and low-score patients in each immune cell/pathway type were further explored. A contrast was conducted according to the median immune cell or immune pathway enrichment level. The immune cell composition in the subgroups was analyzed carefully through the Cell type Identification By Estimating Relative Subsets of RNA Transcripts (CIBSORT) algorithm (http://cibersort.stanford.edu) with “corrplot” R package. The relationship between immune cell infiltration and ERGs expression was explored based on TIMER database (Tumor Immune Estimation Resource) , as well as the gene copy number variation (CNVs) of ERGs. The analyzed immune cells include B cells, CD4 + T cells, CD8 + T cells, macrophages, neutrophils, and dendritic cells.
Potential treatment compounds
To find out a suitable potential clinical treatment that could be used for EC patients, we used the Genomics of Drug Sensitivity in Cancer (GDSC) (https://www.cancerrxgene.org/) database to calculate the half-maximal inhibitory concentration (IC50) of the compounds for predicting the sensitivity via “pRRophetic” R package .
Verification of target ERGs in databases
The protein expression in EC tissues and normal tissues was verified according to the HPA (The Human Protein Atlas) . And the mRNA expression was examined from the CCLE (Cancer Cell Line Encyclopedia).
Quantitative real-time polymerase chain reaction
A quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) was used to verify the expression of the target ERGs in EC cell lines. The total RNA was extracted from the EC cells with TRIZOL reagent (Takara, Otsu, Japan) according to the protocol, and then a 20 µl corresponding cDNA was reversely transcribed with Hiscript@ QRT 157 SuperMix (Vazyme, Nanjing, China). The PCR reaction steps were then followed as different configurations: 95 °C for 30 s, 95 °C for 5 s, and 60 °C for 1 min with 40 cycles. The qRT-PCR process was performed through the CFX Connect Real-Time PCR Detection System (Bio-Rad, 161 Hercules, CA, USA) with SYBR green supermix (Vazyme, Nanjing, China). β-Actin acted as an endogenous control to normalize the expression of each target ERG, and the relative expression levels were calculated with the 2−ΔΔCt method (ΔCt = ΔCt target – ΔCt β-actin). The primer sequences were already listed in Table 1.
The statistical analyses were performed with R in version 4.1.0. The classification of variables in the training and testing groups was decided by the Person chi-square test. The overall survival was analyzed via Kaplan–Meier method along with log-rank test. An independent prognostic analysis was conducted by uni- and multi-Cox regressions. The analyses of several typical characteristics such as the clinicopathological factors, the risk score, and the immune infiltration levels, were assessed using the Wilcoxon test. In general, P < 0.05 was considered as being significantly different.
Identifying 4 prognostic-related differentially expressed ER stress genes
The process of this study was exhibited in the flow chart (Fig. 1). We already obtained 23 normal samples and 554 tumor samples with mRNA expression, and the clinical data from TCGA in total. 178 ER Stress-related gene expression was compared practically between the normal and tumor groups, and then 41 ERGs were selected as differentially expressed genes based on |log2FC|>1 and FDR < 0.05, among them, 12 samples (CAV1, ITPR1, LRRK2, THBS1, ATF3, BCL2, SERP2, MAP3K5, CLU, PRKN, CREB3L2, PPP1R15A) were downregulated, while 29 samples (LONP1, GRINA, EIF2AK1, HYOU1, AUP1, BAX, EDEM2, AIFM1, HM13, PDIA3, ERMP1, DNAJC10, XBP1, CXCL8, ATP2A1, MANF, DERL3, CALR, BAK1, ERO1A, PPP1CA, PDIA4, CREB3L3, P4HB, SDF2L1, TRIB3, RNF183, PDIA2, ERN2), were upregulated. The expression level of the 41 DEGs was presented in the heatmap (Fig. 2A). We further conducted the PPI manipulation (Fig. 2B), which presented the direct interaction between the proteins encoded by DEGs to better understand the potential roles of ER Stress in manipulating EC progression. 132 ERGs were identified as hub genes with a restriction of the interaction score as the highest confidence of 0.9, among them 32 were DEGs (Table 2).
Construction and validation of the 4 ERGs-based risk signature for prognostic effectiveness
To identify the candidates for constructing a risk signature, we first filtered 6 genes related to the prognosis of EC according to the uni-Cox regression analysis with P < 0.05 (BAK1, CREB3L3, DNAJC10, PPP1R15A, TRIB3, XBP1) (Fig. 3A). Then, the LASSO regression analysis confirmed that these 6 genes are related to the survival of the EC patients (Fig. 3B, C). To avoid overfitting the prognostic signature, we applied the minimum likelihood of deviance of the first-rank value of Log(λ) during LASSO regression analysis. We also performed a multi-Cox regression analysis, and then 4 genes were finally identified as significant risk signature genes (CREB3L3, PPP1R15A, TRIB3, XBP1). Among them, PPP1R15A and XBP1 are related to the protective effect with hazard ratio (HR) < 1, whereas CREB3L3 and TRIB3 are the increased risk factors with HR > 1. The calculation formula of the risk scoring is as below:
After screening out, we first sorted out 523 clinical samples from TCGA. Next, we randomly divided the patients into the training and testing groups with 1:1 ratio (of 262 samples in the training group and 261 samples in the testing group). The risk scores of the training and testing groups were then calculated separately. Based on the median score, all samples from both the training and testing groups were classified into high- and low-risk subgroups. As demonstrated, patients from the low-risk subgroup already exhibit a better prognosis with statistically lower risk and fewer deaths than the high-risk subgroup (Fig. 4A, B). Furthermore, the OS between the high- and low-risk subgroups were significantly different with P < 0.05 (Fig. 4C). The AUCs of the training and testing cohorts were 0.711 or 0.678 at 1 year, 0.696 or 0.625 at 2 years, 0.690 or 0.610 at 3 years, respectively (Fig. 4D). The above data imply that the risk score can be used as a valid factor for effectively predicting the prognosis of EC patients. Both the PCA and tSNE analyses indicate that the risk signature constructed by us can be utilized to distinguish high-risk patients from low-risk patients (Fig. 4E, F). The heatmaps constructed by us can also be used to visualize the variance of the prognostic ERGs expression between the high- and low-risk subgroups (Fig. 4G). It should be noted that the risk score can be used as an independent predictor from uni-Cox and multi-Cox regression analyses that the HR of the risk score and 95% confidence interval (95% CI) was 1.422, and 1.244–1.627 (P < 0.05) in uni-cox regression while 1.717, and 1.009–1.360 (P < 0.05) in multi-cox regression (Fig. 5A).
Nomograph construction and the clinical correlation
According to 4 independent prognostic factors such as risk score, age, stage, and grade (all p < 0.05 in uni-Cox and multi-Cox), we built a clinical nomogram to predict the prognosis of the EC patients through a comprehensive score (Fig. 5B). We also utilized the 1- and 2- and 3-year calibration plots to present or prove that the nomogram had a reliable prediction efficiency, which can be used to exhibit excellent accordance between the predictive and actual survival in both the training and testing groups (Fig. 5C, D). As shown, the diagonal line presents the ideal prediction. Furthermore, the performed ROC analyses show that the AUC of the 1-and 2- and 3-year predictive survival according to the nomograph are 0.691and 0.689and 0.676 in the training group, and then 0.677 and 0.623 and 0.608 in the testing group, which means that the nomograph can be effectively used to predict (Fig. 5E, F).
Enriched functions based on the risk signature
To further illustrate the biological functions and the enriched pathways of the risk signature, we performed both the KEGG pathway analysis and the GO enrichment analysis between the high- and low-risk subgroups. The KEGG demonstrates that the pathways enriched in the low-risk subgroup are: protein export, alpha-linolenic acid metabolism, fatty acid metabolism, ether lipid metabolism, and peroxisome, respectively, whereas the pathways in the high-risk subgroup are: DNA replication, RNA polymerase, type II diabetes mellitus, glycosaminoglycan biosynthesis chondroitin sulfate and heparan sulfate (Fig. 6A). Moreover, the GO showed that genes are mainly enriched in the response to the factors such as endoplasmic reticulum stress, endoplasmic reticulum lumen, and endoplasmic reticulum protein-containing complex, which indicates that the genes chosen by us are closely related to the ER Stress to some extent (Fig. 6B).
Relationship between 4 ERGs-based risk signature and clinical characteristics
The correlation between the survival probability and different clinicopathological characteristics was further interrogated. The data demonstrate that the patients from the low-risk subgroup present much better survival overall, and the risk scores are significantly related to the issues such as age, grade, and stage features (Table 3). The patients with a higher stage and grade level are more likely to present significant differences based on the risk score (Fig. 6C). Instructively, the risk signature constructed by us is convictive and reliable for diagnosing and treating EC patients.
Immune infiltration differences between two subgroups
Generally, immune abnormalities will play a critical role in oncogenesis, thus we continuously analyzed the enrichment of the immune cells and then immune pathways based on the risk score through a single sample gene set enrichment analysis (ssGSEA) method. According to the box plot, cells such as B cells naïve, plasma cells, T cells CD4 memory resting, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, mast cells activated, and neutrophils, are all sensitive to the risk score with p < 0.05 criteria. Among them, only macrophages M1 and macrophages M2 are higher in the high-risk subgroup (Fig. 7A). Continuously, we analyzed the survival probability of the samples including immune cells, dendritic cells activated, NK cells activated, plasma cells, T cells CD8, and T cells regulatory (Tregs), which all exhibited a significant difference between the high- and low-risk subgroups (Fig. S1A). The dendritic cells were activated and the plasma cells showed a better prognosis in the low-risk subgroup, whereas the NK cells were activated, and the T cells CD8 and T cells regulatory (Tregs) were in the high-risk subgroup. It is evident that the risk signature was closely related to immune cell infiltration, thus we performed the CIBERSORT to further explore the tumor microenvironment (TME). The compositions of 22 immune cells in the samples and the relative percentage filtered in the high- and low-risk subgroups were analyzed further (Fig. 7B, C). The high-risk subgroup had a higher proportion of anti-inflammatory macrophage M2 cells, while the low-risk subgroup had an increased proportion of resting memory CD4 T cells, activated dendritic cells and mast cells, indicating the potential anti-tumor effects of the low-risk subgroup. Meanwhile, we further investigate the correlation between the abundance of several typical immune cells and the mRNA expression of 4 ERGs through TIMER database (Fig. 8). Macrophages are significantly related to all 4 ERGs. TRIB3 expression level is only correlated with the abundance of another immune cell, neutrophil. CREB3L3 is strongly correlated with B cells and CD4+ T cells with similarly positive tendencies. Except for CD4 + T cells, XBP1 shows significantly negative correlations with CD8 + T cells, neutrophils, and dendritic cells. PPP1R15A is positively related to the abundance of neutrophils and dendritic cells. Compared with the diploid/normal group, we further figure out that immune cell infiltration is also related to CNVs of 4 ERGs (Fig. 9). Arm-level deletion of TRIB3, XBP1, and PPP1R15A decreases the infiltration of immune cells. Compared with the arm-level deletion, the arm-level gain of CREB3L3 is more likely to decrease the infiltration of CD8+ T cells, neutrophils, and dendritic cells. High amplification of PPP1R15A mainly affects B cells, CD8+ T cells, macrophages, and dendritic cells. TRIB3, CREB3L3, and XBP1 are lack of deep deletion.
Besides the immune cells, we have also explored the relationship between the immune pathways and the risk score. The box plot also shows the typical characteristics of the microstructures including B_cells, CCR, checkpoint, iDCs, aDCs, HLA, neutrophils, NK_cells, T_helper_cells, TIL, Treg, T_cell_co-inhibition, T_cell_co-stimulation, Type_I_IFN_Response, and Type_II_IFN_Response, are significantly different according to the risk score (Fig. 7D). The survival analysis of the immune pathways was further conducted, and the results indicate that the immune pathways with prominent survival differences, all have better survival in the high-risk subgroups than in the low-risk subgroup (Fig. S1B). Because the immune checkpoints are significantly different between the two subgroups, the expression distinction was then analyzed, where only a few are expressed higher in the high-risk subgroup, like CD 40. As demonstrated, the majority of the immune checkpoints are higher in the low-risk subgroup, especially corresponding to CTLA4 and CD28. Unfortunately, there is no significant difference in PD-1 or PD-L1 (Fig. 7E).
The above outcomes suggest that the risk score will exhibit a considerable impact on the prognosis of various immune cells and pathways, which may helpfully predict new individualized immunotherapy for EC patients.
Consensus clusters selection based on the expression level
The consensus clustering analysis of all 523 EC samples in the TCGA cohort for recognizing the relationship between the expression levels of 4 risk-signature ERGs and subtypes, was performed. We already selected the clustering variable (k) = 3 from 2 to 9 based on a higher intragroup correlation and a lower intergroup correlation (Fig. 10A), with 225 cases in cluster 1, 142 cases in cluster 2, and 156 cases in cluster 3. Survival analysis indicates that there exists an evident difference in survival among the three subgroups above. As compared to Cluster 2 and Cluster 3, Cluster 1 has a much better prognosis (Fig. 10B).
The issues including the information on mRNA expression, clusters, and clinical characteristics concerning grade and stage, and age, were integrated and displayed in the heatmap. As shown, the parameter grade and stage were significantly different among the three clusters (Fig. 10C). Furthermore, we further explored the expression differences of 30 immune checkpoints in three clusters, among them, the samples of CD244, TNFSF9, Lag3, CD200, TNFSF15, HHLA2, ICOSLG, TNFRSF25, CD40, TNFRSF14, TNFSF14, CD44, and BTNL2, are more statistically significant than others (p < 0.001) (Fig. 10D).
Diverse potential chemotherapies of two subgroups
Because chemotherapy is one of the effective and widely-accepted treatments for EC patients, we further did a prediction of the sensitivity of potential chemotherapies based on IC50 for each sample in both the low- and high-risk subsgroups. Generally, the patients with a lower risk score are more sensitive to medicines such as AKT inhibitor VIII, Bicalutamide, Docetaxel, and Temsirolimus. While the higher-risk score patients should be more sensitive to Pyrimethamine, Shikonin, Rucaparib, and Veliparib (Fig. 11). These findings might provide new options for future clinical treatment.
Dual verification of the 4 ERGs-based risk signature
The protein level of 4 risk-signature ERGs (CREB3L3, PPP1R15A, TRIB3, XBP1) was verified according to the HPA database (Fig. 12A). Unfortunately, we couldn’t obtain successfully a protein expression of XBP1 in the HPA database. Besides this gene, we find that both the CREB3L3 and TRIB3 are upregulated, but PPP1R15A is downregulated on the protein level. The mRNA levels of 4 ERGs are different in 40 EC cell lines (Fig. 12B).
We then performed qRT-PCR to verify the mRNA expression of 4 risk-signature ERGs (CREB3L3, PPP1R15A, TRIB3, XBP1) in EC and normal tissues, separately. The expression level of these genes differs significantly between the tumor and normal tissues, but CREB3L3 and TRIB3, and XBP1 are overexpressed, while PPP1R15A is expressed conversely (Fig. 12C). The primer sequences are provided in Table 1. Overall, the experimental results further confirm the reliability of the risk signature (*p < 0.05, **p < 0.01, and ***p < 0.001).
EC is one of the highly heterogeneous malignancies developed in the female genital tract. A great number of studies have reported that endoplasmic reticulum stress is closely related to the progression and outcome of different types of cancer, including the EC [11, 29], but the mechanisms have not yet been fully discovered. Furthermore, most studies investigated the impact of ER Stress on cancer biological activities, just a few focused on the prognostic significance of ER Stress-related genes, particularly in the case of EC. Therefore, it is necessary to develop a predictive risk signature based on the relationship between both EC and ER Stress. Here, we have built a stable prognostic signature to predict the prognostic ER Stress-related genes and then divided them into different subgroups, to systemically learn the relationship among genes, tumor microenvironment, and immune activities. It should be noted that we further predict the potential chemotherapies that may be effectively utilized in future clinical treatment.
First of all, we successfully constructed a reliable 4-gene ER Stress-related prognostic signature for EC with TCGA gene expression and clinical files. 4 prognostic ERGs (CREB3L3, PPP1R15A, TRIB3, and XBP1) were screened from 178 ER Stress-related genes via Cox and LASSO regression analyses. Two subgroups were classified based on the risk score. A nomograph was employed to prove the reliable predictive value of the risk signature. OS results suggested a strong prognostic correlation between EC and ER Stress. The prediction ability of the risk signature was obviously affected by the clinical features that patients with lower risk scores presented better prognoses.
The functional enrichment highlighted the interactive mechanisms between ER Stress and endometrial carcinogenesis. GO enrichment suggested that ERGs are mainly involved in apoptotic regulation besides ER Stress related to signaling pathways. The mutual effect between ER Stress and apoptosis in driving cancers has been well explained . KEGG enrichment presenting pathways enriched in the high-risk subgroup show more carcinogenesis correlation, like DNA replication. The metabolism-related pathways are highly enriched in the low-risk subgroup, such as fatty acid, protein, and RNA metabolism, to confirm that the patients from the high-risk subgroup exhibit unsatisfactory prognoses and immune-suppressive status.
Considering the previous studies have demonstrated the critical role of immune-related activities in endometrial cancer progression, we further analyze the differences in immune infiltration between the subgroups. From ssGSEA results, B cells naïve, plasma cells, and resting memory CD4+ T cells are evidently downregulated in the high-risk subgroup. Both naïve B cells and plasma cells are critical subpopulations of B cells, naïve B cells are the precursor of functional B cells and plasma cells are recognized as the basis of humoral immunity [31, 32]. B cells and plasma cells are associated with better survival, and particularly B-cell markers can prolong survival specifically in high-grade tumors . Memory T cells play important roles in T cell persistence and tumor immunotherapy efficacy . Thus, this observation may demonstrate the higher chance of the high-risk subgroup to be a suppressive immune status than the lower-risk subgroup. By contrast, M1 and M2 macrophages and mast cells are upregulated in the high-risk subgroup, indicating their relationship with the unfavorable prognosis, which has been proved by previous studies .
CIBERSORT analysis demonstrates the top one cell enriched in the low-risk subgroup with a significant difference is resting memory CD4 T cells, indicating the protective role in the EC progression. Similar to the ssGSEA result, the high-risk subgroup had a higher proportion of anti-inflammatory macrophage M2 cells. So far, the immune pathway analysis has paid more attention to CCR and DCs, and IFN responses. According to the differences in checkpoints, we further find that several crucial immune checkpoints, such as CTLA 4, and CD 28, are differentially expressed. The drug targeting immune checkpoints could be used alone or combined with other appropriate treatments, and thus present an advantage of lower toxicity and much higher effectiveness to some extent . In general, ERGs will influence immune infiltration and prognosis. It should be noted that understanding the mechanisms is not enough for the EC patients’ treatment decision-making. Thus, the potential chemotherapies are continuously discussed. Patients in the high-risk subgroup are more sensitive to the components such as AKT inhibitor VIII, Bicalutamide, Docetaxel, and Temsirolimus, whereas the patients in the low-risk subgroup will benefit from the medications such as Pyrimethamine, Shikonin, Rucaparib, and Veliparib.
CREB3L3 (cAMP-responsive element-binding protein 3 like 3) is a membrane-bound transcription factor located in the ER and structurally similar to ATF6 that could be proteolytically activated by ER Stress . CREB knockdown in macrophages could downregulate M2 marker genes, then improve insulin resistance, suggesting that CREB is important in maintaining insulin sensitivity in white adipose tissue via its initiation of the innate immune system .
PPP1R15A (protein phosphatase 1 regulatory subunit 15 A), also known as GADD34, is a stress-inducible eIF2α phosphatase, which can accelerate cell death by boosting protein synthesis and activating death-related pathways . PPP1R15A has been proven to be a hypoxia/autophagy-related gene in breast cancer respectively [40, 41]. However, there are just a few studies focusing on the action of PPP1R15A in endometrial cancer progression to date, which requires more studies.
TRIB3 is one of the Tribbles Pseudokinase homologs with a weak ATP affinity. Despite the absence of kinase activity, TRIB3 is an adaptor/scaffold protein for numerous functional proteins . TRIB3 promotes endometrial malignant actions via elevating CTNNB1 transcription . And it’s found that TRIB3 expression change is part of the anti-endometrial cancer activities in one ongoing clinical trial .
XBP1 (X-box binding protein 1) belongs to the CREB/ATF basic leucine zipper (bZIP) protein family. XBP1 can be alternatively spliced into two isoforms, the unspliced isoform (XBP1-u) is spliced into the spliced isoform (XBP1-s) under ER Stress . XBP1 is highly enriched in endometrioid endometrial tumors . One recent study indicates that XBP1 potentially distinguishes the polymerase epsilon exonuclease (POLE) from the copy number (CN)-low endometrial cancer subtype . Therefore, XBP1 may be a strong promoter of endometrial carcinogenesis. Otherwise, some studies suggest the dual roles of XBP1. XBP1-u inhibits autophagy  and protects cells from oxidative stress . XBP1-s promotes cell death by altering calcium levels under strong ER Stress, while under non-fatal ER Stress, XBP1-s protects cell survival by transcriptionally regulating UPR target genes . Together, it may explain that XBP1 is identified as the protective factor in our risk signature, but expressed higher in EC tissues than normal tissues.
Despite already presenting an obvious value of the risk signature developed by us, some limitations and deficiencies still exist. First off, the RNA-sequence and the sample data used by us are only acquired from TCGA, which might cause bias in the following calculations. It may be a better approach by combining the data from other resources to derive comprehensive results. Second, the specific mechanisms based on ERGs should be dug out through more experiments. Last but not least, as the patients from the high-risk subgroup present a relatively worse prognosis, closer follow-ups should be required to further prove the risk signature and nomograph viability.
In this paper, an effective ER Stress-related risk signature based on 4 prognostic genes is successfully established by us. The risk signature can accurately predict the prognosis and help treatment-decision making for EC patients. The fundamental relationship between immunity, TME, and EC-specific ER Stress, is discussed carefully through the risk signature. We also verified the reliability of the risk signature with basic experiments.
Availability of data and materials
- ER Stress:
endoplasmic reticulum stress
The cancer genome atlas
Gene set enrichment analysis
The human protein atlas
Cancer Cell Line Encyclopedia
The genomics of drug sensitivity in cancer
ER Stress-related genes
Differentially expressed genes
Receiver operating characteristic
Area under the curve
the Kyoto Encyclopedia of Genes and Genomes
the Gene Ontology
Copy Number Variation
Half-maximal inhibitory concentration
Real-time quantitative RT-PCR
cAMP-responsive element-binding protein 3 like 3
Protein phosphatase 1 regulatory subunit 15 A
Tribbles pseudokinase 3
X-box binding protein 1
Crosbie EJ, Kitson SJ, McAlpine JN, Mukhopadhyay A, Powell ME, Singh N. Endometrial cancer. Lancet. 2022;399(10333):1412–28.
Passarello K, Kuria S, Villanueva V. Endometrial Cancer: an overview of Pathophysiology, Management, and Care. Semin Oncol Nurs. 2019;35(2):157–65.
Tung HJ, Huang HJ, Lai CH. Adjuvant and post-surgical treatment in endometrial cancer. Best Pract Res Clin Obstet Gynaecol. 2022;78:52–63.
Oakes SA, Papa FR. The role of endoplasmic reticulum stress in Human Pathology. Annu Rev Pathol. 2015;10:173–94.
Wang G, Yang ZQ, Zhang K. Endoplasmic reticulum stress response in cancer: molecular mechanism and therapeutic potential. Am J Transl Res. 2010;2(1):65–74.
Cubillos-Ruiz JR, Bettigole SE, Glimcher LH. Tumorigenic and immunosuppressive effects of endoplasmic reticulum stress in Cancer. Cell. 2017;168(4):692–706.
Mahadevan NR, Rodvold J, Sepulveda H, Rossi S, Drew AF, Zanetti M. Transmission of endoplasmic reticulum stress and pro-inflammation from tumor cells to myeloid cells. Proc Nati Acad Sci USA. 2011;108(16):6561–6.
Yi Y, Jiao P, Wang Z, Chen M, Du H, Xu L, et al. Endoplasmic reticulum stress promotes the release of exosomal PD-L1 from head and neck cancer cells and facilitates M2 macrophage polarization. Cell Commun Signal. 2022;20(1):12.
Wei W, Zhang Y, Song Q, Zhang Q, Zhang X, Liu X, et al. Transmissible ER stress between macrophages and tumor cells configures tumor microenvironment. Cell Mol Life Sci. 2022;79(8):403.
Ulianich L, Insabato L. Endoplasmic reticulum stress in endometrial cancer. Front Med (Lausanne). 2014;1:55.
Bifulco G, Miele C, Jeso BD, Beguinot F, Nappi C, Carlo CD, et al. Endoplasmic reticulum stress is activated in endometrial adenocarcinoma. Gynecol Oncol. 2012;125(1):220–5.
Cali G, Insabato L, Conza D, Bifulco G, Parrillo L, Mirra P, et al. GRP78 mediates cell growth and invasiveness in endometrial cancer. J Cell Physiol. 2014;229:1417–26.
Yixuan W, Xin L, Jiayin W, et al. A joint model considering measurement errors for optimally identifying tumor mutation burden threshold. Front Genet. 2022;13:915839.
Yixuan W, Jiayin W, Wenfeng F, et al. TMBserval: a statistical explainable learning model reveals weighted tumor mutation burden better categorizing therapeutic benefits. Front Immunol. 2023;14:1151755.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021;49(D1):D605–12.
Francesco B, Terrie W, Richard B, et al. The relationship between different settings of medical service and incident frailty. Exp Gerontol. 2018;108:209–14.
Friedman J, Hastie T, Tibshirani R. Regularization Paths for generalized Linear Models via Coordinate Descent. J Stat Softw. 2010;33(1):1–22.
Jinghua W, Xinqi Q, Jiayu H, et al. Development and validation of a novel mitophagy-related gene prognostic signature for glioblastoma multiforme. BMC Cancer. 2022;22(1):644.
Liqiang G, Qiong W, Zhen M, et al. Identification of immune-related genes that predict prognosis and risk of bladder cancer: bioinformatics analysis of TCGA database. G (Albany NY). 2021;13(15):19352–74.
Kanehisa M, Goto S. KEGG: Kyoto Encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.
Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28:1947–51.
Kanehisa M, Furumichi M, Sato Y, Kawashima M, Ishiguro-Watanabe M. KEGG for taxonomy-based analysis of pathways and genomes. Nucleic Acids Res. 2023;51:D587–92.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
Li T, Fu J, Zeng Z, Cohen D, Li J, Chen Q, et al. TIMER2.0 for analysis of tumor-infiltrating immune cells. Nucleic Acids Res. 2020;48(W1):W509–14.
Geeleher P, Cox N, Huang RS. pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS ONE. 2014;9(9):e107468.
Uhlén M, Fagerberg L, Hallström BM, Lindskog C, Oksvold P, Mardinoglu A, et al. Proteomics. Tissue-based map of the human proteome. Science. 2015;347(6220):1260419.
Tierney KE, Ji L, Dralla SS, Yoo E, Yessaian A, Pham HQ, et al. Endoplasmic reticulum stress in complex atypical hyperplasia as a possible predictor of occult carcinoma and progestin response. Gynecol Oncol. 2016;143(3):650–4.
Verfaillie T, Garg A, Agostinis P. Targeting ER stress induced apoptosis and inflammation in cancer. Cancer Lett. 2013;332(2):249–64.
Michael PC, Mary MT. Memory B cells and plasma cells: the differentiative continuum of humoral immunity. Immunol Rev. 2021;303(1):72–82.
Sarah AO, Jennifer AW, Margaret LH, et al. Plasma cell development and survival. Immunol Rev. 2010;237(1):140–59.
Mandal G, Biswas S, Anadon C, Yu X, Gatenbee C, Prabhakaran S, et al. IgA-Dominated Humoral Immune responses govern patients’ outcome in Endometrial Cancer. Cancer Res. 2022;82(5):859–71.
Makoto A, Minako I, Tanakorn S, et al. Memory T cell, exhaustion, and tumor immunity. Immunol Med. 2020;43(1):1–9.
Ribatti D, Finato N, Crivellato E, Marzullo A, Mangieri D, Nico B, et al. Neovascularization and mast cells with tryptase activity increase simultaneously with pathologic progression in human endometrial cancer. Am J Obstet Gynecol. 2005;193(6):1961–5.
Rowshanravan B, Halliday N, Sansom D. CTLA-4: a moving target in immunotherapy. Blood. 2018;131(1):58–67.
Zhang K, Shen X, Wu J, Sakaki K, Saunders T, Rutkowski D, et al. Endoplasmic reticulum stress activates cleavage of CREBH to induce a systemic inflammatory response. Cell. 2006;124(3):587–99.
Luan B, Yoon YS, Lay JL, Kaestner KH, Hedrick S, Montminy M. CREB pathway links PGE2 signaling with macrophage polarization. Proc Natl Acad Sci U S A. 2015;112(51):15642–7.
Liu L, Ito S, Nishio N, Sun Y, Chen N, Tanaka Y, Isobe KI. GADD34 facilitates cell death resulting from Proteasome Inhibition. Anticancer Res. 2015;35(10):5317–24.
Wang J, Wang Y, Xing P, Liu Q, Zhang C, Sui Y, Wu C. Development and validation of a hypoxia-related prognostic signature for breast cancer. Oncol Lett. 2020;20(2):1906–14.
Zhong S, Chen H, Yang S, Feng J, Zhou S. Identification and validation of prognostic signature for breast cancer based on genes potentially involved in autophagy. PeerJ. 2020;8:e9621.
Stefanovska B, André F, Fromigué O. Tribbles Pseudokinase 3 regulation and contribution to Cancer. Cancers (Basel). 2021;13(8):1822.
Wang W, Hong G, Chien P, Huang Y, Lee H, Wang P, et al. Tribbles Pseudokinase 3 contributes to Cancer Stemness of Endometrial Cancer cells by regulating β-Catenin expression. Cancers (Basel). 2020;12(12):3785.
Muñoz-Guardiola P, Casas J, Megías-Roda E, Solé S, Perez-Montoyo H, Yeste-Velasco M, et al. The anti-cancer drug ABTL0812 induces ER stress-mediated cytotoxic autophagy by increasing dihydroceramide levels in cancer cells. Autophagy. 2021;17(6):1349–66.
Luo X, Alfason L, Wei M, Wu S, Kasim V. Spliced or Unspliced, that is the question: the Biological Roles of XBP1 Isoforms in Pathophysiology. Int J Mol Sci. 2022;23(5):2746.
Berger A, Korkut A, Kanchi R, Hegde A, Lenoir W, Liu W, et al. A Comprehensive Pan-Cancer Molecular Study of gynecologic and breast cancers. Cancer Cell. 2018;33(4):690–705e9.
Kim K, Hwangbo S, Kim H, Kim Y, No J, Suh D, et al. Clinicopathologic and protein markers distinguishing the polymerase epsilon exonuclease from the copy number low subtype of endometrial cancer. J Gynecol Oncol. 2022;33(3):e27.
Zhao Y, Li X, Cai M, Ma K, Yang J, Zhou J, et al. XBP-1u suppresses autophagy by promoting the degradation of FoxO1 in cancer cells. Cell Res. 2013;23:491–507.
Martin D, Li Y, Yang J, Wang G, Margariti A, Jiang Z, et al. Unspliced x-box-binding protein 1 (XBP1) protects endothelial cells from oxidative stress through Interaction with histone deacetylase 3. J Biol Chem. 2014;289:30625–34.
This work was supported by the International Science and technology cooperation project of Hubei Province (GJHZ2021000024), the National Natural Science Foundation of China (81974409), the Fundamental Research Funds for the Central Universities (2019kfyXKJC072), the Young Scientists Fund of the National Natural Science Foundation of China (82103157).
Ethics approval and consent to participate
Consent for publication
The authors declare no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Further survival analysis linked with immune infiltration. A Kaplan-meier analysis of significantly different immune cells. B Kaplan-meieranalysis of significantly different immune pathways.
About this article
Cite this article
Zhang, T.a., Zhang, Q., Zhang, J. et al. Identification of the role of endoplasmic reticulum stress genes in endometrial cancer and their association with tumor immunity. BMC Med Genomics 16, 261 (2023). https://doi.org/10.1186/s12920-023-01679-5