Identification of the role of endoplasmic reticulum stress genes in endometrial cancer and their association with tumor immunity

doi:10.21203/rs.3.rs-2827667/v1

Download PDF

Research Article

Identification of the role of endoplasmic reticulum stress genes in endometrial cancer and their association with tumor immunity

https://doi.org/10.21203/rs.3.rs-2827667/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 25 Oct, 2023

Read the published version in BMC Medical Genomics →

You are reading this latest preprint version

Background: 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.

Method: 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.

Result: 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.

Conclusion: We established an ER Stress-based risk signature that could effectively predict EC patients' prognosis and their immune correlation.

Endometrial cancer

Endoplasmic reticulum stress

Bioinformatic analysis

Risk signature

Prognosis

Immune infiltration

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 [1]. 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 [2]. Despite many obvious achievements in EC treatment, such as surgery, radiotherapy, cytotoxic chemotherapy, and hormonal therapy [3], 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 [4]. In turn, the continuous activation of the ER Stress will alter the tumorigenesis through both the transcriptional and translational pathways [5]. Furthermore, ER Stress can reprogramme the function of immune cells and thus apparently temper anti-cancer immunity [6]. ER Stress induces tumor cells to release factors that change the leucocytes around tumor cells, thus promoting tumor growth [7]. ER-stressed tumor cells secrete PD-L1 that drives M2 macrophage polarization, facilitating the escape from immune surveillance [8]. On the other hand, ER-stressed M1-like macrophages secrete DAMPs that lead to tumor cell death [9].

So far, continuously accumulated studies have presented that ER Stress should have played a critical role in endometrial carcinogenesis [10]. Protein and mRNA levels of several ER Stress-induced apoptosis indicators, including GRP78 and ATF6, are apparently increased in endometrial cancer immunohistochemical samples [11]. And GRP78 plasma membrane localization is elevated in endometrial cancer [12]. Therefore, cancer therapy focused on ER Stress is promising for deeper exploration.

In this paper, a bioinformatic method for constructing an ER Stress-based risk signature and then verifying the expression of the target prognostic genes by in vitro experiments are proposed by us for exhibiting a new insight into the correlation between ER Stress and EC.

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 [13], and then following |log₂fold change (FC)|>1 and false discovery rate (FDR) < 0.05. The protein-protein intersection (PPI) networking was conducted by utilizing search tools for the retrieval of interacting genes (STRING) (http://string-db.org/) [14].

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 [15]. 4 ERGs were eventually selected to calculate the risk score. The formula for performing risk score is presented below

$${\text{r}\text{i}\text{s}\text{k} \text{s}\text{c}\text{o}\text{r}\text{e}}^{}={\sum }_{k}^{n}{\text{e}\text{x}\text{p}\text{r} (\text{E}\text{R}\text{G}\text{s}）}^{k}{\times \text{c}\text{o}\text{e}\text{f} (\text{E}\text{R}\text{G}\text{s}）}^{k}$$

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) enrichment analyses based on the DEGs between the high- and low-risk subgroups by applying “clusterProfiler” R package [16]. 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 [17]. 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) [18], 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 [19].

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) [20]. 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 seconds, 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 3.

Statistical Analysis

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.

dentifying 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 |log₂FC|>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 protein-protein interaction (PPI) manipulation (Fig. 2B), which already presents a direct interaction between the ERGs. 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 1).

Table 1

List of 32 DEGs with 0.9 interaction score between EC tissue and normal endometrial tissue in TCGA.
AIFM1	CLU	ERN2	PDIA3
ATF3	CREB3L2	HM13	PDIA4
AUP1	CREB3L3	HYOU1	PPP1CA
BAK1	CXCL8	ITPR1	PPP1R15A
BAX	DERL3	LRRK2	SDF2L1
BCL2	DNAJC10	MANF	THBS1
CALR	EDEM2	MAP3K5	TRIB3
CAV1	EIF2AK1	P4HB	XBP1

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:

Risk Score = (0.3965)×(CREB3L3 exp.)+(-0.4152)×(PPP1R15A exp.)+(0.2223)×(TRIB3 exp.)+(-0.3508)×(XBP1 exp.).

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 effectively utilized to well 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 2). 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.

Table 2

Endometrial cancer patient clinical characteristics from TCGA
Variables	Number of cases	Percentage(%)
Age(year) < 65/≥65	288/247	53.8/46.2
Grade G1/G2/G3	99/1212/316	18.4/22.7/58.8
Stage I-II/III-IV	390/147	72.6/27.4

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, the samples such as B cells naive, 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 are activated and the plasma cells show a better prognosis in the low-risk subgroup, whereas the NK cells are activated, and the T cells CD8 and T cells regulatory (Tregs) are in the high-risk subgroup. It is evident that the ERGs are closely related to immune cell infiltration, thus we already 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). 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 3. Overall, the experimental results further confirm the reliability of the risk signature. (*p < 0.05, **p < 0.01, and ***p < 0.001)

Table 3

The primers sequence of 4 prognostic-related ERGs
Gene	Primer Sequences
TRIB3-F	TGCGTGATCTCAAGCTGTGT
TRIB3-R	GCTTGTCCCACAGGGAATCA
CREB3L3-F	ATCTCCTGTTTGACCGGCAG
CREB3L3-R	GTCGTCAGAGTCGGGGTTTG
XBP1-F	AGGAGTTAAGACAGCGCTTGGGGATGGAT
XBP1-R	CTGAATCTGAAGAGTCAATACCGCCAGAAT
PPP1R15A-F	CTGGCTGGTGGAAGCAGTAA
PPP1R15A-R	TATGGGGGATTGCCAGAGGA

As indicated, the 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, 21], 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.

In this study, we have summarized 178 genes related to ER stress from GSEA and GeneCards. The mRNA expression files were downloaded and the clinical files were screened out from the TCGA. At first, we identified 4 prognostic ERGs (CREB3L3, PPP1R15A, TRIB3, and XBP1) and then constructed a risk signature based on those 4 prognostic ERGs via Cox and LASSO regression analyses. We then constructed a nomograph that already proves the reliable predictive value of the risk signature. In combination with survival data, the OS suggests a strong correlation between EC and ER Stress. The prediction of the risk signature is obviously affected by the clinical characteristics of the patients with a low-risk score has presented a better survival than a high-risk score. Based on the risk signature, we can easily divide the EC patients into high- and low-risk subgroups. In addition, from GO and KEGG pathway enrichment analyses, we find that the ERGs are closely related to processes including apoptosis, metabolism, and immunological activities. Furthermore, the ERGs are closely correlated with immune infiltration.

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 [22]. Luan B et al. found that a precise knockdown of CREB or CRTC2/3 in macrophages could downregulate M2 marker gene expression, and then improve insulin resistance, suggesting that CREB will play an important role in maintaining insulin sensitivity in white adipose tissue via its initiation of the innate immune system [23].

PPP1R15A (protein phosphatase 1 regulatory subunit 15A), also known as GADD34, is a stress-inducible eIF2α phosphatase, which can accelerate cell death by boosting protein synthesis and activating death-related pathways [24]. Wang and Zhong et al. firmly established that PPP1R15A is a hypoxia-related and autophagy-related gene in breast cancer respectively [25, 26]. To date, there are just a few studies focusing on the action of PPP1R15A in endometrial cancer progression, which warrant more study.

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 other functional proteins [27]. Wang et al. demonstrated that the oncogenic mechanism of TRIB3 in endometrial cancer could elevate CTNNB1 transcription [28]. Guardiola et al. reported that anti-endometrial cancer molecular undergoing the phase 2 clinical trial could induce TRIB3 expression, so as to result in a macroautophagy/autophagy-mediated cancer cell death. And the change in TRIB3 level could be detected in the blood of patients enrolled in the trial [29].

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 its spliced isoform (XBP1-s) under endoplasmic reticulum stress, and both isoforms are r involved in various biological pathways [30]. Berger et al. screened the TCGA database and then summarized that XBP1 can be highly enriched in endometrioid endometrial tumors [31]. A more recent study indicates that XBP1 is a potential marker to distinguish the polymerase epsilon exonuclease (POLE) from the copy number (CN)-low endometrial cancer subtype [32]. Zhao et al. reported that the depletion of XBP1 inhibits a release of the pro-tumor cytokines of re-educating tumor-associated macrophages (TAMs), but enhances the macrophage phagocytosis, therefore driving CRC progression [33]. As a result, XBP1 may be a strong candidate for the cause of endometrial cancer treatment.

So, functional enrichment between the subgroups was conducted to figure out the possible mechanisms between ER Stress and endometrial cancer progression. 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 [34]. 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 already demonstrating a critical role of immune-related activities in endometrial cancer progression, we further analyze the differences in immune infiltration and pathways between the subgroups. The cell level such as B cells, plasma cells, and resting memory CD4 T cells, are predominantly downregulated in the high-risk subgroup, meaning the suppressive immune status. Mandal et al. found that B cells and plasma cells could predict better survival, and particularly B-cell markers can prolong survival specifically in high-grade tumors [35]. By contrast, M1 and M2 macrophages and mast cells are upregulated in the high-risk subgroup than in the low-risk, and thus indicate their relationship with the unfavorable prognosis, which is already proved by previous findings [36].

CIBERSORT analysis demonstrates a top one enriched in the low-risk subgroup with a significant difference is a kind of resting memory CD4 T cells, indicating a protective role in the EC progression. 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 [37]. In general, ERGs will influence immune infiltration and prognosis. It should be noted that knowing the mechanisms is not enough for the EC patient clinical treatments. So, the potential chemotherapies are continuously discussed. As demonstrated, patients in the high-risk subgroup are more sensitive to the issues 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.

Notably, XBP1 is mainly recognized as the promoter of carcinogenesis. However, recent studies have suggested the dual roles of XBP1-u and XBP1-s. XBP1-u could degrade FoxO1 and inhibit autophagy [38]. XBP1-u also participates in the HDAC3/mTOR/pAkt1/Nrf2/HO-1 pathway and protects cells from oxidative stress [39]. When ER stress varies, XBP1-s play different roles. XBP1-s promotes cell death by altering calcium levels under strong ER Stress, while under non-fatal ER stress, XBP1-s reduces ER stress and promotes cell survival by transcriptionally regulating UPR target genes [30]. Together, it may explain the reason that XBP1 is identified as the protective factor in our risk signature, but expressed higher in tumor samples.

In summary, ER Stress should be a critical factor in the occurrence and development of endometrial cancer according to our works. 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 risk signaturebased on 4 ER Stress-related genes for predicting the prognosis of EC patients is successfully established by us. A fundamental relationship between immunity, TME, and the potential EC treatments, is analyzed carefully according to some key predictions and detailed data analyses, and finally verified through basic experiments. It should be noticed that further experiments must be conducted to sufficiently verify our findings and the explode deeper mechanisms. It can be expected that the method will highlight the continuous development of EC.

EC: endometrial cancer

ER stress: endoplasmic reticulum stress

TCGA: The cancer genome atlas

GSEA: Gene set enrichment analysis

HPA: The human protein atlas

CCLE: Cancer Cell Line Encyclopedia

GDSC: The genomics of drug sensitivity in cancer

ERGs: ER Stress-related genes

DEGs: Differentially expressed genes

OS: Overall survival

ROC: Receiver operating characteristic

AUC: Area under the curve

PPI: Protein-protein interaction

KEGG: the Kyoto Encyclopedia of Genes and Genomes

GO: the Gene Ontology

CNVs: Copy Number Variation

IC50: Half-maximal inhibitory concentration

qRT-PCR: Real-time quantitative RT-PCR

CREB3L3: cAMP-responsive element-binding protein 3 like 3

PPP1R15A: protein phosphatase 1 regulatory subunit 15A

TRIB3: Tribbles pseudokinase 3

XBP1: X-box binding protein 1.

Acknowledgments

Not applicable.

Authors' contributions

JZ and QZ designed the bioinformatic study. RS and RZ obtained data. TASZ and STW performed analyses. TASZ wrote and SGL revised the manuscript. HBW and QZ supervised the study. All authors read and approved the final manuscript.

Funding

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).

Data availability

All the data involved in this study were collected from open sources, TCGA (http://portal.gdc.cancer.gov), GDSC (http://www.cancerrxgene.org/), TIMER (https://cistrome.shinyapps.io/timer/), CCLE (http://www.broadinstitute.org/ccle), and HPA (http://www.proteinatlas.org/) databases.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no conflicts of interest.

Crosbie EJ, Kitson SJ, McAlpine JN, Mukhopadhyay A, Powell ME, Singh N. Endometrial cancer. Lancet. 2022;399(10333):1412–1428.
Passarello K, Kuria S, Villanueva V. Endometrial Cancer: An Overview of Pathophysiology, Management, and Care. Semin Oncol Nurs. 2019;35(2):157-165.
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 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 Transmission of endoplasmic reticulum stress and pro-inflammation from tumor cells to myeloid cells. Proc Nati Acad Sci USA. 2011;108(16):6561-6566.
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, InsabatoEndoplasmic 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–1426.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, ShiW, 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, PyysaloS, 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-D612.
Friedman J, Hastie T, Tibshirani Regularization Paths for Generalized Linear Models via Coordinate Descent. J Stat Softw. 2010;33(1):1-22.
Yu G, Wang LG, Han Y, He clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284-287.
Hänzelmann S, Castelo R, Guinney GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013;14:7.
Li T, FuJ, 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-W514.
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, YessaianA, 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-654.
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-599.
LuanB, Yoon YS, Lay JL, Kaestner KH, Hedrick S, Montminy CREB pathway links PGE2 signaling with macrophage polarization. Proc Natl Acad Sci U S A. 2015;112(51):15642-7.
LiuL, 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.
WangJ, 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-1914.
ZhongS, Chen H, Yang S, Feng J, Zhou Identification and validation of prognostic signature for breast cancer based on genes potentially involved in autophagy. PeerJ. 2020;8:e9621.
StefanovskaB, André F, Fromigué Tribbles Pseudokinase 3 Regulation and Contribution to Cancer. Cancers (Basel). 2021;13(8):1822
WangW, 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-GuardiolaP, 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-1366.
Luo X, AlfasonL, Wei M, Wu S, Kasim Spliced or Unspliced, That Is the Question: The Biological Roles of XBP1 Isoforms in Pathophysiology. Int J Mol Sci. 2022;23(5):2746.
BergerA, 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-705.e9
KimK, 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
ZhaoY, Zhang W, Huo M, Wang P, Liu X, Wang Y, et al. XBP1 regulates the protumoral function of tumor-associated macrophages in human colorectal cancer. Signal Transduct Target Ther. 2021;6(1):357
VerfaillieT, Garg A, Agostinis Targeting ER stress induced apoptosis and inflammation in cancer. Cancer Lett. 2013;332(2):249-64.
MandalG, 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-871.
RibattiD, 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.
RowshanravanB, Halliday N, Sansom CTLA-4: a moving target in immunotherapy. Blood. 2018;131(1):58-67
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–30634.

No competing interests reported.

supplementary.docx

Download PDF

Journal Publication

published 25 Oct, 2023

Read the published version in BMC Medical Genomics →

Editorial decision: Major revision
28 Jul, 2023
Reviews received at journal
20 Jul, 2023
Reviewers agreed at journal
10 Jul, 2023
Reviewers invited by journal
09 Jul, 2023
Editor assigned by journal
09 Jul, 2023
Editor invited by journal
06 Jun, 2023
Submission checks completed at journal
20 Apr, 2023
First submitted to journal
17 Apr, 2023

You are reading this latest preprint version

Identification of the role of endoplasmic reticulum stress genes in endometrial cancer and their association with tumor immunity

Status:

Journal Publication

Version 1

Abstract

Figures

Introduction

Materials And Methods

Normal and Tumor Datasets

Identification Of Er Stress-related Differentially Expressed Genes

Establishment And Validation Of The Risk Signature Based On Differentially Expressed Genes

Independent Prognostic Factor Analysis

Function Prediction Analysis

Nomogram And Calibration Analysis

Tumor Immune Analysis

Potential Treatment Compounds

Verification Of Target Ergs In Databases

Quantitative Real-time Polymerase Chain Reaction

Statistical Analysis

Results

dentifying 4 prognostic-related differentially expressed ER Stress genes

Construction And Validation Of The 4 Ergs-based Risk Signature For Prognostic Effectiveness

Nomograph Construction And The Clinical Correlation

Enriched Functions Based On The Risk Signature

Relationship Between 4 Ergs-based Risk Signature And Clinical Characteristics

Immune Infiltration Differences Between Two Subgroups

Consensus Clusters Selection Based On The Expression Level

Diverse Potential Chemotherapies Of Two Subgroups

Dual Verification Of The 4 Ergs-based Risk Signature

Discussion

Conclusion

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1