Endoplasmic reticulum stress regulators exhibit different prognostic, therapeutic and immune landscapes in pancreatic adenocarcinoma

Abstract Endoplasmic reticulum stress (ERS) and unfolded protein response are the critical processes of tumour biology. However, the roles of ERS regulatory genes in pancreatic adenocarcinoma (PAAD) remain elusive. A novel ERS‐related risk signature was constructed using the Lasso regression analysis. Its prognostic value, immune effect, metabolic influence, mutational feature and therapeutic correlation were comprehensively analysed through multiple bioinformatic approaches. The biofunctions of KDELR3 and YWHAZ in pancreatic cancer (PC) cells were also investigated through colony formation, Transwell assays, flow cytometric detection and a xenograft model. The upstream miRNA regulatory mechanism of KDELR3 was predicted and validated. ERS risk score was identified as an independent prognostic factor and could improve traditional prognostic model. Meanwhile, it was closely associated with metabolic reprogramming and tumour immune. High ERS risk enhanced glycolysis process and nucleotide metabolism, but was unfavourable for anti‐tumour immune response. Moreover, ERS risk score could act as a potential biomarker for predicting the efficacy of ICBs. Overexpression of KDELR3 and YWHAZ stimulated the proliferation, migration and invasion of SW1990 and BxPC‐3 cells. Silencing KDELR3 suppressed tumour growth in a xenograft model. miR‐137 could weaken the malignant potentials of PC cells through inhibiting KDELR3 (5′‐AGCAAUAA‐3′). ERS risk score greatly contributed to PAAD clinical assessment. KDELR3 and YWHAZ possessed cancer‐promoting capacities, showing promise as a novel treatment target.


| INTRODUC TI ON
Pancreatic adenocarcinoma (PAAD), a common abdominal cancer, is associated with high malignancy and mortality.According to the American Cancer Society, PAAD accounts for >400,000 deaths in the United States each year, contributing to 4.7% of all cancerrelated deaths. 1 The 5-year overall survival rate (OSR) of patients with PAAD is <35%. 2 Approximately 10% of patients with PAAD are diagnosed with localized carcinoma and are eligible for surgical excision. 3The incidence of pancreatic cancer (PC) has been increasing. 3Recent advances in PAAD treatment have not improved the therapeutic outcomes.FOLFIRINOX (fluorouracil, leucovorin, irinotecan and oxaliplatin) and gemcitabine, which are the mainstream chemotherapy regimens for PAAD, have not increased the median disease-free survival to >25 months. 4Additionally, programmed cell death protein 1/ligand 1 (PD-1/L1) inhibitors such as pembrolizumab exert beneficial effects in <14.2% of patients. 5Thus, there is an urgent need to develop new therapeutic strategies and modify the currently used prognostic assessment system.Endoplasmic reticulum (ER) stress (ERS), a hallmark of cancer, is involved in tumour onset and progression.
Extrinsic and intrinsic factors, such as hypoxia, increased mutation rate, elevated demand for protein synthesis, and oncogene activation disrupt the homeostasis of protein synthesis in ER, which is called ER stress (ERS). 6To maintain proteostasis, cells activate a unique mechanism, termed the unfolded protein response (UPR).UPR can promote protein folding and enhance the clearance of misfolded proteins that have accumulated due to ERS. 7 Additionally, UPR is regulated by three stress transducers, namely ATF6, EIF2AK3 and ERN1. 8ERS and UPR are reported to regulate tumour cell proliferation, progression, immune regulation, immunotherapy response and drug resistance. 9For example, the ERS-related gene SERP1 is associated with the prognosis and immune response in skin cutaneous melanoma. 10ERS-induced exosomal miR-27a-3p promotes immune escape in breast cancer by regulating CD274 expression in macrophages. 11These findings indicate that targeting ERS can be a potential therapeutic strategy for cancer.
KDELR3 belongs to the KDEL (Lys-Asp-Glu-Leu) receptor (KDELR) family, whose members regulate protein transport between Golgi and ER. 12 Mechanistically, KDELR members can bind to the heterotrimeric G protein in the Golgi complex, promoting initial protein transport. 13Previous studies have reported that KDELR3 functions as an adaptive gene in response to unfolded or misfolded protein accumulation in the ER.Thus, KDELR3 is a critical regulator of ERS and UPR and is associated with the pathogenesis of various human diseases, such as cancers.For example, a microarray study demonstrated that KDELR3 is a biomarker of atherosclerosis (AS) onset. 14In uveal melanoma (UM), KDELR3 is associated with prognosis, immune infiltration and chemoresistance. 15The roles of KDELR3 in PAAD have not been elucidated.
Several studies have elucidated the roles of ERS regulatory genes in cancer clinical assessment using high-throughput approaches.
Zhang et al. constructed an ERS signature to evaluate immune features in glioma. 16The prognostic value of ERS signature in various cancers, such as glioma, 17 hepatocellular carcinoma (HCC), 18 lung adenocarcinoma (LUAD) 19 and renal cell carcinoma, 20 has been previously evaluated.However, the clinical significance of ERS-related genes in PAAD has not been well evaluated.This study comprehensively examined the roles of ERS-related genes in the prognosis, tumour immune microenvironment (TIM), and mutational and metabolomic features of PAAD.To screen the candidates who can benefit from immune checkpoint blockade (ICB) therapy, the correlation of ERS risk score with the efficacy of ICB therapy was examined.
Furthermore, the pro-oncogenic roles of two core ERS signature genes, KDELR3 and YWHAZ, were verified using in vitro and in vivo experiments.A novel regulatory axis miR-137/KDELR3 was also be proved to play critical functions in PAAD progression.The findings of this study can contribute to improve the clinical assessment and treatment of PAAD.

| Establishment of ERS-related gene set
This study established a comprehensive ERS-related gene set using the Molecular Signatures Database (MSigDB) (https:// www.gseamsigdb.org/ gsea/ msigdb/ ).Some reviews 8,9,21 have suggested that UPR is the major mechanism activated in response to ERS, which is mediated by ERN1, ATF6 and EIF2AK3.Hence, the following keywords were queried in the MSigDB: 'Endoplasmic reticulum stress', 'Unfolded protein response', 'IRE1', 'ATF6', and 'PERK'.In total, 11 core gene sets were obtained (Table S4).After removing the overlapping genes, the final gene set comprised 379 ERS-related genes (Table S5).To further confirm the reliability of the ERS gene set, biological functional analyses were performed using the Metascape online tool (http:// metas cape.org/ ).

| Construction of ERS-related risk signature
The ERS-related risk signature was constructed using a four-step process.The ERS-related DEGs were identified using the 'Limma' package in R software (Ver 4.1.2)based on the following criteria: adjusted p < 0.05 and |Log 2 fold-change (FC) ≥ 1|.Next, the ERS regulatory genes with prognostic value were identified using Cox univariate regression analysis.The DEGs were then intersected with prognostic genes using a Venn diagram.To ensure the accuracy of modelling, the genes whose expression patterns were inconsistent with their prognostic values were excluded.Finally, the genes obtained after these three screening rounds were subjected to least absolute shrinkage and selection operator (LASSO) regression analysis to construct a novel ERS risk signature for PAAD.

| Gene Set Enrichment Analysis (GSEA)
GSEA was employed to unravel three issues.First, the differential UPRs between non-cancerous and tumour samples were evaluated.Next, the relationships between ERS risk score and ERS intensity were explored.Third, the effects of ERS risk score on PAAD metabolic pathways, including glycolysis, amino acid (AA) metabolism, nucleotide metabolism and lipid metabolism, were examined.The detailed description of the candidate gene sets is provided in Table S6.GSEA was performed with 1000 permutations and the gene symbols parameter set to 'No_collapse'.The phenotype labels were set as high-risk score samples versus lowrisk score samples.

| Survival analyses
Survival analyses were performed using the Kaplan-Meier method.
The data of samples with a short follow-up duration (<30 days) were excluded.The cohort was divided into high-risk and low-risk groups based on the optimal cut-off value of the ERS risk score, which was determined using the Cutoff Finder online tool (http:// molpa th.chari te.de/ cutoff). 23Cox univariate and multivariate analyses were performed to identify the independent prognostic factors.The predictive accuracy of the ERS risk signature was evaluated using the receiver operating characteristic curve.Decision curve analysis (DCA) was performed to assess whether the ERS score can improve the traditional prognostic model based on TNM staging or clinical stage system.The corresponding modelling process relied on the multivariate logistic regression method.The clinical subgroup analyses were performed to map the applicable range of the ERS prognostic model.The data of patients with M1-stage tumours were not included in the subgroup analyses owing to the limited number of samples (n = 5).A nomogram comprising TNM staging and ERS risk score was constructed to predict the 1-year, 2-year and 3-year OSRs.The prognostic accuracy of the ERS signature was assessed using the calibration curves.The ICGC-PACA-CA, ICGC-PACA-AU, GSE62452, GSE28735 and GSE57495 datasets served as the validation cohorts.

| Immune analyses
The infiltration levels of 22 immune cells in each PAAD sample were calculated using the cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm. 24The differential abundances of immune cells between the high-risk and low-risk groups were compared using the 'Limma' R package.The activities of 13 immune-related pathways were quantified using single-sample GSEA.ESTIMATE algorithm was used to measure the tumour purity and the stromal and immune cell admixture in each PAAD sample. 25,26TIMER, a web server, was used for comprehensively analysing the tumour-infiltrating immune cells (https:// cistr ome.shiny apps.io/ timer/ ). 27The correlation between somatic copy number alterations (SCNAs) and the abundance of immune cells was determined.

TA B L E 1
The purpose and description of the datasets used in this study.

| Mutational analyses
The somatic mutation data of 149 PAAD samples were obtained from TCGA database.'Data Type' and 'Workflow Type' were set to 'Masked somatic mutation' and 'VarScan', respectively.The tumour mutational burden (TMB) value of each PAAD sample was calculated as follows: TMB = total mutation frequency/38. 28The somatic mutation data were visualized and curated using the 'maftools' R package. 29The cBioPortal database (http:// cbiop ortal.org) provided the mutational frequency and pattern of seven hub ERS-related genes across four PC projects (n = 532 samples). 30

| Therapeutic correlation analyses
The correlation between the efficacy of ICB therapy and ERS risk score was analysed.TMB is a useful indicator for the selection of therapeutics for ICB therapy. 31A high TMB is correlated with favourable responses to ICB therapy.Hence, the differential TMBs between the high-risk and low-risk groups were determined.The tumour immune dysfunction and exclusion (TIDE) score predicts the response of patients to anti-PD-1/L1 and anti-CTLA4 treatments based on the estimation of T-cell dysfunction and tumour immune evasion. 32The TIDE score of each risk group was calculated using an online tool (http:// tide.dfci.harva rd.edu/ login/ ).Patients exhibiting upregulated expression levels of immune checkpoints (ICs) exhibit favourable responses to ICB therapy. 33,34Hence, the correlation between the expression levels of six critical ICs and the ERS risk score was analysed.T-cell functions and cytotoxic effects are closely associated with the efficacy of immunotherapy. 35,36Thus, the differential T-cell functions and cytotoxic effects between the high-risk and low-risk groups were examined.ImmuCellAI can be used to predict the response to ICB therapy based on the abundance of 24 immune cells in samples (http:// bioin fo.life.hust.edu.cn/ ). 37IMvigor210, a real clinical cohort, was used to examine the differential ERS risk scores in patients exhibiting differential outcomes to treatment with the PD-1 inhibitor atezolizumab. 38

| Human Protein Atlas (HPA) database analysis
The HPA database provides the proteome profiles of 32 different cancers in the form of immunohistochemical images (https:// www.prote inatl as.org/ ). 39The histological expression levels of seven hub ERS regulators in non-cancerous and PC tissues were analysed.

| Cell culture and transfection
In vitro experiments were performed with three PC cell lines (SW1990, BxPC-3 and PANC-1) and one healthy pancreatic cell line (HPDE6-C7 targeting two ERS signature genes (sh-KDELR3 and sh-YWHAZ) and their overexpression (OE-KDELR3 and OE-YWHAZ) lentivirus constructs were purchased from HanHeng Biotechnology (Shanghai, China).The specific sequences of these constructs are shown in Table S7.The cells were transfected with lentiviruses (HanHeng Biotechnology, Shanghai, China).

| Quantitative real-time polymerase chain reaction (qRT-PCR)
Total RNA was extracted using chloroform and TRIzol reagent (TaKaRa, Japan).The optical density (OD) at 260 and 280 nm of the sample was measured using an ultraviolet spectrophotometer (Nanodrop 2000 spectrophotometer).An OD 260 to OD 280 ratio in the range of 1.8-2.0indicates a high RNA purity.These RNA samples were subjected to further reverse transcription using the PrimeScript RT reagent kit (TaKaRa, Japan).qRT-PCR analysis was performed using SYBR-Green PCR Reagent (Takara, Japan) with the ABI Prism 7900 system.The PCR conditions were as follows: 95°C for 30 s (pre-denaturation), followed by 40 cycles at 95°C for 5 s, 95°C for 15 s, 60°C for 30 s and 95°C for 15 s.GAPDH was used as the internal reference.The relative gene expression was calculated using the 2 −ΔΔC T method.The primer sequences used in qRT-PCR analysis are shown in Table S8.

| Immunohistochemical staining (IHC)
IHC was performed on four pairs of clinical samples from Second Affiliated Hospital of Xi'an Jiaotong University to validate expressions of KDELR3 and YWHAZ in PAAD.All patients signed informed consents before IHC detection.IHC process was analogous to a previous study. 40Briefly, paraffin sections was prepared through dehydration, embedding, dewaxing and rehydration.Subsequently, tissue sections were treated with 1.5% H 2 O 2 solution to inactivate endogenous peroxidase.Citric acid buffer was used for antigen repair.After blocking by primary and secondary antibodies, tissue sections were stained by ABC method.

| Flow cytometric analysis
Transfected cells were fixed with precooled 70% anhydrous ethanol at −20°C.Next, the cells were resuspended in phosphate-buffered saline (PBS) and incubated with 100 μL RNase A (50 μL/mL) at 37°C for 30 min.The cells were then stained with propidium iodide (50 μg/ mL, BD Biosciences) in the dark for 30 min.Flow cytometric analysis was performed using the Becton Dickinson FACScan system.

| Colony formation assay
At 48 h post-transfection, cells in the logarithmic growth phase were seeded into six-well plates at a density of 1 × 10 3 /well.The colonies were fixed and stained with Giemsa stain in methanol.The number of colonies in five random visual fields was counted under the microscope.

| Transwell migration and invasion assays
Transfected cells (1 × 10 4 per well) were seeded into the 24-well Transwell chambers (Corning, NY, USA).The medium supplemented with 0.1% FBS was added to the upper chambers, while the medium supplemented with 10% FBS was added to the lower chambers.After culturing for 24 h, the cells on the upper surface of the membrane were removed using PBS.The cells that invaded the membrane were fixed with paraformaldehyde and stained with 0.1% crystal violet.
The number of cells in five random visual fields was counted under a high-magnification microscope (100-fold).The upper chambers were precoated with Matrigel to perform the invasion assays.

| Dual-luciferase assay
Dual-luciferase assay was utilized to verify the miRNA binding sites with KDELR3.The recombinant vectors of wild type (WT) and mutation type (MUT) of KDELR3 3′UTR region were designed and synthesized by Genechem (Genechem Incorporation, Shanghai, China).
The detailed protocol was similar to the previous description. 41iefly, miRNA vectors and KDELR3 recombinant vectors were cotransfected 293T cells.Then, luciferase reagent was added in lysed cells.The relative fluorescence intensity of each experimental group was measured using a microplate reader.

| Xenograft assay
Female BALB/c nude mice aged 6 weeks were used to perform the tumour xenograft experiments.BxPC-3 cells that were stably transfected with sh-KDELR3 were subcutaneously injected into the flanks of each mouse.The tumour volume was calculated as follows: tumour volume = 0.5 × (tumour length) × (tumour width) 2 .Tumour length and width were measured using a vernier calliper once every 5 days.
After 2 weeks, all mice were euthanized, and the xenograft tumours were excised.This study was approved by the Ethics Committee of the Second Affiliated Hospital of Xi'an Jiaotong University.

| Statistical analysis
All statistical analyses were performed using the R software (version 4.1.2) and GraphPad Prism (version 8.0.1).The continuous variables between the groups were compared using the t-test or Kruskal-Wallis test.Meanwhile, the categorical variables were compared using the Wilcoxon rank sum test.Correlation analyses were performed using the Spearman method.Prognostic meta-analysis was performed using Review Manager 5.2 software (The Cochrane Collaboration, Oxford, UK) with the Mantel-Haenszel method.The odds ratio value was applied as the prognostic evaluation index.The I 2 value was used to assess the statistical heterogeneity of meta-analysis.The fixed effect model was used if the I 2 value was <50%, while the random effect model if the I 2 value was ≥50%.The overall effects were tested using Z-test.The in vitro experiments were repeated three times independently.Differences were considered significant at p < 0.05.

| Establishment of a comprehensive and reliable ERS-related gene set
The workflow of this study is shown in Figure 1A The ERN1 pathway promotes the expression of transcription factors (TFs), especially XBP1. 42Additionally, ERN1 can directly target mRNAs and miRNAs through regulated IRE1-dependent decay (RIDD), regulating the translation process. 43The EIF2AK3 pathway downregulates protein translation by phosphorylating EIF2A, which alleviates the burden of protein folding in ER. 21In the ATF6 pathway, ATF6 is transported to the Golgi apparatus and cleaved by the proteases SP1 and SP2.Subsequently, cleaved ATF6 releases its functional fragment ATF6f, which functions as a TF 44 and promotes the degradation of misfolded proteins.
Based on the mechanisms of ERS and UPR, 11 gene sets were obtained from the MSigDB to establish a comprehensive ERSrelated gene set (Figure 2A).The protein-protein interaction (PPI) network of 379 ERS regulatory genes is shown in Figure S1A.Critical ERS-related genes, such as ERN1, EIF2AK3 and ATF6, were located in the core module (Figure 2B).Biological function analyses revealed that selected ERS-related genes were enriched in ERS and UPR processes.These findings indicated that the ERS gene set established in this study was reliable (Figure 2C).

| UPR is upregulated in PC samples
The differential activities of UPR between non-cancerous and PC samples were analysed using GSEA.As shown in Figure 2D-I, the UPR pathways in PC samples were significantly upregulated when compared with those in non-cancerous samples.In particular, only the ATF6 signalling pathway was markedly upregulated in PC samples (Figure 2G), whereas the ERN1 and EIF2AK3 pathways were unaffected (Figure 2H,I).Based on these findings, we hypothesized that the upregulation of UPR, especially the ATF6 pathway, is a biological hallmark of PAAD.

| Construction of a novel ERS-related risk signature for PAAD
Several ERS regulators (24.8%; 94/379) were differentially expressed between non-cancerous and tumour samples (Figure 3A,B).Most of these DEGs were upregulated in PC samples.Approximately 30% of ERS-related genes (28.8%; 109/379) were associated with the survival outcomes of patients with PAAD (Figure S1B).In total, 32 ERSrelated genes exhibited differential expression and had prognostic values (Figure 3C).The expression patterns of some intersection genes were not consistent with their prognostic values.For example, EEF2 was upregulated in tumour samples (Log2FC = 1.356) but predicted favourable prognosis (hazard ratio (HR) = 0.451).These 'contradictory genes' (n = 11) were not included in LASSO regression analysis to prevent their interference with modelling.As shown in Figure 3D, when partial-likelihood deviance (PLD) was the minimum, the best value of Log(λ) was slightly larger than −3.At this time, the model fitting degree of ERS risk signature was optimal, and it was consisted of seven variables.Similarly, when Log(λ) took the optimum value (Around −3), there were only seven genes whose coefficients did not decay to 0 (Figure 3E).  to predict the 1-year, 2-year and 3-year survival rates of patients with PAAD (Figure 4S).For example, the 2-year OSR of a 50-year-old patient who was diagnosed as T2N0M0 was predicted to be approximately 80% when the patient was classified into the low-risk group.
In contrast, the 2-year OSR of this patient was predicted to be <60% when the patient was classified into the high-risk group.The calibration plots indicated that the predicted survival rates were consistent with the actual survival rates (Figure S2).

| ERS risk signature can be adapted to multiple external cohorts
To enable the application of the ERS risk signature, the prognostic

| ERS risk score indicates metabolic changes, especially glycolysis and nucleotide metabolism
Metabolic changes are the core hallmarks of cancers.As shown in Table 2, aerobic glycolysis, nucleotide metabolism, AA metabolism and lipid metabolism are involved in tumour progression processes, such as cell proliferation and therapy resistance.This study investigated the effects of ERS risk levels on various metabolic pathways.The metabolic reprogramming toward glycolysis, also called the 'Warburg effect', indicates the preference of tumour cells for inefficient glucose utilization via glycolysis rather than oxidative phosphorylation. 45EA revealed that four glycolysis-related genes were significantly enriched in the high-risk group (Figure 6A).This study also analysed the correlation between the ERS risk score and four rate-limiting enzymes of glycolysis.The expression levels of PKM, HK1 and HK2 in the high-risk group were significantly higher than those in the low-risk group (Figure S3A).However, the PFKFB3 expression levels were not significantly different between the high-risk and low-risk groups.The ERS risk score was positively correlated with the expression levels of PKM, HK1 and HK2 (Figure S3B-D) but not with those of PFKFB3 (Figure S3E).Thus, a high ERS risk score may be correlated with active glycolysis and the upregulation of relevant rate-limiting enzymes.
Active nucleotide metabolism, which contributes to the occurrence and progression of tumours, 46 was upregulated in the high-risk group (Figure 6B).The enrichment of AA metabolism and lipid metabolism was not markedly different between the high-risk and low-risk groups (Figure 6C

| Differential landscapes of TIM between the high-risk and low-risk groups
The abundances of 22 immune cells in each PAAD sample are shown in Figure S4A.Immune cell enrichments significantly varied between the high-risk and low-risk groups.The levels of infiltrating B cells, CD8+ T cells, activated memory CD4+ T cells, activated dendritic cells and neutrophils were downregulated, whereas those of macrophages were upregulated in the high-risk group (Figure 7A).As shown in Table 3, these changes in the abundance of immune cells suppressed anti-tumour immune responses and promoted immune tolerance.For example, CD8+ T cells function as the most powerful effectors in the anti-cancer immune response and are the targets for cancer immunotherapy. 47Hence, the downregulation of CD+ T cells inhibits the anti-tumour immune process.
Among the immune-related signalling pathways, the activities of type-II interferon (IFN) response, T-cell functions, cytolytic activity and checkpoint were suppressed in the high-risk group (Figure 7B).
Type-II IFN is a pleiotropic molecule with anti-proliferative and pro-apoptotic properties. 48The suppression of type-II IFN leads to unfavourable anti-tumour immune responses.Similar findings were observed in 'ESTIMATE' analyses.The stromal, immune and ESTIMATE scores in the high-risk group were markedly lower than those in the low-risk group (Figure 7C).The tumour purity in the high-risk group was significantly higher than that in the low-risk group, which can be attributed to the suppression of anti-tumour immune responses in the high-risk group (Figure 7D).The SCNAs than those in the low-risk group (Figure 8F,G).As shown in Table 4, the upregulation of these genetic mutations promotes the malignant behaviours of PAAD, which may be the intrinsic reason for the poor prognosis of the high-risk group.The mutation profile of ERS signature genes was conserved across multiple PC projects (Figure 8H).
The mutation frequencies of these genes were less than 10%.

| ERS risk score is a potential biomarker for predicting ICB therapy response
ICB therapy offers a novel paradigm for cancer treatment.

Bioinformatic analyses revealed the potential correlation between
ERS risk score and the efficacy of ICB therapy.The correlation patterns varied depending on the analysis (Figure 9).A high ERS risk score indicated a good response to ICB therapy.As shown in Figure 9A,B, the TMB in the high-risk group was significantly higher than that in the low-risk group and was positively correlated with the ERS risk score.A high TMB results in the production of tumour neoantigens and consequently activates the immune system to recognize tumour cells.Thus, patients with high TMB benefit from ICB therapy. 49,50The TIDE score in the high-risk group was markedly lower than that in the low-risk group (Figure 9C).This suggests that patients in the high-risk group are less likely to exhibit T-cell dysfunction and immune evasion (Figure 9C).
A high ERS risk score was associated with unfavourable ICB therapy response (Figure 9D-J).The expression levels of BTLA, TIGIT, LAG3 and HACR2, but not CD274 and CTLA4, in the highrisk group were downregulated when compared with those in the low-risk group (Figure 9D).Additionally, the expression levels of BTLA, TIGIT, LAG3 and HACR2 were negatively correlated with the ERS risk score (Figure 9E-J).The upregulation of ICs is considered to be a marker of favourable ICB therapy response. 33,51us, patients with a high-risk score may not be sensitive to ICB therapy.Moreover, the functions of immunotherapy-related signalling pathways, including 'checkpoint', 'T-cell functions' and 'cytolytic activity', were downregulated in the high-risk group, resulting in unfavourable ICB therapy response (Figure 9K).The predictive results of the online tool 'ImmuCellAI' were consistent with these findings.Patients in the high-risk group were predicted to exhibit unfavourable ICB therapy responses.In contrast, the predicted response rate in the low-risk group was 4.2% (4/95) (Figure 9L).
A real clinical cohort (IMvigor 210) was used to further investigate the correlation between ERS risk score and the efficacy of ICB (atezolizumab) therapy (Figure 9M,N).The ERS risk score was not significantly different between the responder and non-responder groups (Figure 9M).The clinical response rates were similar between the high-risk (24.2%) and low-risk (21.5%) groups (Figure 9N).

| ERS signature genes are differentially expressed in PC tissues
The histological expression levels of seven ERS signature genes were confirmed using the HPA database (Figure S5A).The BCL2L1, CAV1, KDELR3 and YWHAZ expression levels were upregulated in the tumour tissues and downregulated in the healthy pancreatic tissues.
However, NCCRP1 was upregulated in tumour tissues and not detected in healthy tissues.This was not consistent with the mRNA expression patterns of these proteins.Furthermore, USP13 was downregulated in both cancer and healthy tissues.

| KDELR3 promotes malignant behaviours and modulates the cell cycle of PC cells
Several studies have examined the roles of five ERS signature genes in PAAD (Table 5).However, limited studies have examined the roles of KDELR3, USP13 and YWHAZ.KDELR3 and YWHAZ were selected for further analysis as they exhibited the high coefficient values in the ERS risk model (0.115 and 0.227).

TA B L E 3
The immune microenvironment profiles of the high-risk group.

TA B L E 4
The clinical implications of characteristic mutations in multiple cancers.
detections on 2 pairs of clinical samples also supported above findings (Figure 10A).sh-KDELR3 and OE-KDELR3 effectively altered the mRNA and protein levels of KDELR3 in PC cell lines (Figure 10B,C).The number of colonies in the OE-KDELR3transfected group was higher than that in the other groups, indicating that KDELR3 overexpression upregulated the proliferation of PC cells (Figure 10D,E).In contrast, cell proliferation was downregulated in the sh-KDELR3-transfected group (Figure 10D,E).Moreover, KDELR3 overexpression and knockdown exerted contrasting effects on the cell cycle of PC cells.
Flow cytometric analysis revealed that KDELR3 overexpression significantly promoted the transition of cell cycle from the G0/G1 phase to the S phase, whereas KDELR3 knockdown arrested the cell cycle at the G0/G1 phase (Figure 10F-H).

| YWHAZ also has promotive effects on the malignant behaviours of pancreatic carcinoma cells
As for YWHAZ, IHC detection confirmed that YWHAZ was upregulated in tumour samples compared to normal ones (Figure 12A).
Specific overexpression vector and shRNA could effectively alter the mRNA expressions of YWHAZ in BxPC3 and SW1990 cells (Figure 12B).The effectiveness of above genetic tools was also validated through western blot assays (Figure 12C).Clonogenic assays revealed that overexpression of YWHAZ promoted the proliferation of PC cells, whereas silencing YWHAZ inhibited this process (Figure 12D,E).Moreover, PC cells with YWHAZ overexpression exhibited more aggressive migratory and invasive capacities through Transwell assays (Figure 12F-H).Quantitative data analyses also confirmed that migrative and invasive cells in overexpression groups were significantly higher than that in control groups, but stained cells in silencing groups were obviously decreased (Figure 12G-I).
Clearly, YWHAZ enhanced the malignant potentials of PC cells.

| KDELR3 knockdown suppresses tumour growth in the xenograft model
The tumour burden in nude mice injected with sh-KDELR3transfected cells was lower than that in nude mice injected with negative control-transfected cells (Figure 13A).Mice were sacrificed to obtain the xenograft tumour.KDELR3 knockdown markedly suppressed xenograft tumour growth (Figure 13B).Western blot test manifested that the protein expressions of KDELR3 in xenograft tumours from the KDELR3 deletion group (sh-KDELR3, S1-S6) were significantly lower than that in the negative control group (sh-vector, C1-C6) (Figure 13C), which ensured the accuracy of the experiments.Tumour weight and volume in mice injected with sh-KDELR3transfected cells were significantly lower than that in mice injected with sh-vector-transfected cells (Figure 13D,E).These findings indicate that KDELR3 knockdown suppresses PAAD growth.

| Potential upstream regulatory mechanisms of KDELR3
This study examined the potential upstream regulatory mechanisms of KDELR3.The upstream miRNAs of KDELR3 were predicted.Six candidate miRNAs were identified from three miRNA databases (Figure 14A), The regulatory network of these miRNAs is shown in Figure 14B.miRNAs negatively regulate downstream target genes by binding to their 3′-untranslated region (UTR).Thus, the candidate miRNAs were also hypothesized to negatively regulate KDELR3.Correlation analyses revealed that miR-137 negatively regulated KDELR3 (Figure 14C-H, Cor = −0.225,p = 0.003).Meanwhile, the binding site of miR-137 in KDELR3 was predicted using the TargetScanHuman database (Figure 14I).The findings of this study indicated that miR-137 suppresses KDELR3 expression by targeting its 3′-UTR (5′-AGCAAUAA-3′).
The upstream TFs of KDELR3 were obtained from the hTFtarget and Genecards databases.In total, 13 TFs were predicted as the potential regulators of KDELR3.The regulatory network of these TFs was mapped (Figure 14J,K).As TFs activate the transcription of downstream genes, the correlation of their expression with KDELR3 expression was examined.The expression levels of HDAC1, BHLHE40 and FOXM1 were positively correlated with those of KDELR3.FOXM1 exhibited the highest correlation coefficient (Cor = 0.342, Figure 14L).
Moreover, GSEA analysis was utilized to explore the potential signalling pathways of KDELR3 in PAAD progression.The results revealed that Wnt and TGFβ signalling pathways significantly enriched in patients with high KDELR3 expression, whereas the enrichments of other signalling pathways were not subjected to KDELR3 expression levels (Figure 14M).

| DISCUSS ION
PC is a major health and economic burden for society.The median survival time of patients with PC is <28 months owing to the high malignancy and metastasis of PC.Currently used treatments are not effective against PC.In response to ERS, the UPR process and various biological processes, such as cell proliferation, oncogene stimulation and mutation induction, are activated in the cells. 52ERS and UPR exhibit pleiotropic functions in cancer progression. 21The roles of ERS and UPR in the immune process, prognosis and progression of PAAD have not been elucidated.This study aimed to evaluate the roles of ERS and UPR in PAAD pathogenesis.Prognostic assessment is crucial for managing patients with cancer.The currently used prognostic evaluation system is associated with some limitations.The American Joint Committee on Cancer eighth edition staging system cannot predict the prognosis of patients with T stage/lymph node-negative PAAD. 53Although the eighth edition of the TNM staging system addresses the issues of extrapancreatic extension, its prognostic accuracy is not superior to that of the seventh edition of the TNM staging system. 54This study constructed a novel ERS risk signature, which increased the decision benefit of the traditional TNM staging system (Figure 4H).
In addition to serving as an independent prognostic factor of PAAD, the ERS risk score can distinguish the differential prognoses of patients with no metastasized lymph nodes (N0 stage) (Figure 4P).
These characteristics addressed the limitations of the TNM system.
Therefore, the ERS risk score can accurately predict the prognosis of patients with PAAD.The increased infiltration levels of macrophages in the high-risk group promoted the onset of immunological tolerance (Figure 7A, Table 3).Thus, the ERS risk score can predict the immune status of patients with PAAD.
ERS can serve as a predictive biomarker for cancer immunotherapy response.Cell and animal experiments have demonstrated that the inhibition of UPR mediators exerts anti-cancer effects. 56This study investigated the correlation between ERS risk score and the efficacy of ICB therapy from multiple perspectives.Although the analytical results were not consistent, the findings of this study suggested that a high ERS risk score adversely affects the efficacy of ICB therapy.This can be attributed to two main reasons.First, TMB may be the consequence of ICB therapy response and not the cause. 58Chen et al. 59 reported that camrelizumab response was not correlated with TMB.Second, the efficacy of immunotherapy is dependent on the immune status of patients.The anti-cancer immune processes were suppressed in the high-risk group (Figure 7).Previous studies have reported that ERS adversely affects the efficacy of immunotherapy by modulating the TME. 6Therefore, a high ERS risk score is an unfavourable biomarker for ICB efficacy.
Metabolic reprogramming is a core hallmark of cancer.This study demonstrated that the metabolic pathways in PC, especially glycolysis, varied between the high-risk and low-risk groups.In tumours, the metabolism is reprogrammed toward glycolysis ('Warburg effect'). 602][63][64] As glycolysis provides the biological precursors required for cell proliferation and confers therapy resistance to tumours, 63 active glycolysis can drive cancer development.
The glycolysis process was significantly enriched in the high-risk group (Figure 6A), which may be the metabolic driver contributing to poor prognosis in high-risk patients.
This study demonstrated that oncogenic mutations were upregulated in the high-risk group.The mutation frequencies of KRAS, TP53 and CDKN2A in the high-risk group were higher than those in the low-risk group (Figure 8F,G).KRAS mutation is a critical genetic alteration with diagnostic value in various tumours, such as PC, LUAD and colorectal cancer. 65Cheng et al. 66 reported that KRAS mutation was closely associated with the upregulation of circulating regulatory T cells and indicated poor prognosis and advanced clinical stage in PAAD.Patients with CDKN2A mutation are at a high risk of developing melanomas and PC. 67The risk of cancerrelated death in patients with melanoma harbouring CDKN2A mutation is three times higher than that in patients with melanoma not harbouring CDKN2A mutation. 680][71] These findings indicate that the characteristic mutations in the high-risk group are the critical epigenetic factors involved in PC progression.
KDELR3 encodes a member of the KDEL ER protein retention receptor family.The ER is a key site for the synthesis of biomacromolecules, such as proteins and lipids.Hence, KDELR3 is involved in the pathogenesis of various human diseases.For example, KDELR3 may serve as a diagnostic and therapeutic biomarker for AS.
. An ERS-related gene set was established based on some crucial reviews.The core mechanisms of UPR are shown in Figure 1B.Intrinsic (oncogene activation, high demand for cell proliferation, genomic instability and increased mutation frequency) and extrinsic (such as hypoxia, shortage of nutrients and acidosis) factors promote the accumulation of misfolded and unfolded proteins in ER.Subsequently, the folding of proteins and the clearance of unfolded proteins are upregulated through UPR, which is mediated by three major signalling pathways.

F I G U R E 1
ERS is associated with cancer pathogenesis.(A) The flow chart of this study.(B) The core mechanism of ERS.ERS, endoplasmic reticulum stress.Finally, a novel ERS-related risk signature was constructed based on 22 core ERS-related genes (Figure 3F).The risk score was calculated as follows (Figure 3H): ERS risk score = [(0.191× BCL2L1 relative expression) + (0.092 × CAV1 relative expression) + (0.088 × DDIT4 relative expression) + (0.115 × KDELR3 relative expression) + (0.073 × NCCRP1 relative expression) + (−0.085 × USP13 relative expression) + (0.227 × YWHAZ relative expression).Principal component analysis (PCA) revealed that the ERS risk score can explain up to 57.5% of the prognostic variance, suggesting a good fit of the ERS model (Figure 3G).Meanwhile, the associations of ERS risk score with the intensity of intracellular ER stress were investigated using GSEA method.The enrichments of ERS and UPR were significantly increased in PAAD samples with high ERS score, indicating that high ERS score was concomitant with high ERS intensity (Figure 3H-K).

3. 4 |F I G U R E 4 F I G U R E 5 | 11 of 26 LIU
Figure 4D).Additionally, the ERS risk score exhibited the best performance for predicting 3-year OSR (AUC = 0.807, Figure4E).Univariate Cox analysis indicated that age, histological grade, value of the ERS risk signature was validated in five external cohorts.The ERS risk score could effectively stratify the prognosis of patients in the ICGC-PACA-AU (HR = 2.16, p = 0.01), GSE62452 (HR = 1.81, p = 0.044) and ICGC-PACA-CA (HR = 1.53, p = 0.023) cohorts.A high-risk score indicated unfavourable survival outcomes (Figure 5A-C) and decreased PFS (Figure 5J).The predictive accuracy of the ERS model was moderate in these cohorts (AUC = 0.641-0.656)(Figure 5E-G).The prognostic values of the ERS model were not observed in the GSE28735 and GSE57495 cohorts (Figure 5D,H,I,K).Survival meta-analysis was performed to comprehensively evaluate the prognostic value of the ERS risk score in five validation cohorts.As shown in Figure 5L, the overall mortality risk in the high-risk group was markedly higher than that in the low-risk group (Z = 2.29, p = 0.02).Although the ERS risk signature did not exhibit prognostic values in two validation cohorts, the overall effects of the ERS risk signature were closely associated with PAAD prognosis.The adjusted funnel plot suggested no publication bias (Figure 5M).

F I G U R E 9 | 17 of 26 LIU
The potential correlation between ERS risk score and the efficacy of ICB therapy.(A) The differential TMB between the high-risk and low-risk groups.(B) The correlation between TMB and ERS risk score was analysed using the Spearman method.(C) The differential TIDE scores between the high-risk and low-risk groups.(D) The differential expression levels of ICs between the high-risk and low-risk groups.(E-J) The correlation of IC expression levels with ERS risk score.(K) The differential activities of four IC-related immune pathways.(L) The ICB therapy response rate in different risk groups was predicted using the ImmuCellAI online tool.(M) The differential risk scores between the responder and non-responder groups in the IMvigor 210 cohort.(N) The response rates in different risk groups.ERS, endoplasmic reticulum stress; ICs, immune checkpoints; ICB, immune checkpoint blockade; TMB, tumour mutation burden; TIDE, tumour immune dysfunction and exclusion.*p < 0.05, **p < 0.01, ***p < 0.001.et al.KDELR3 upregulated the migration and invasion of PC cells.The numbers of migrated cells in the OE-KDELR3-transfected group were significantly higher than those in other groups.In contrast, the number of migrated cells was the least in the sh-KDELR3-transfected group (Figure11A,B).The results of the Transwell invasion assays revealed that KDELR3 overexpression markedly increased the invasive ability of BxPC-3 and SW1990 cells (Figure11C,D).In contrast, KDELR3 knockdown downregulated cell invasion (Figure11C,D).These findings suggest that KDELR3 promotes PC progression.

3. 15 || 19 of 26 LIU
miR-137 mediates the PC progression through inhibiting KDELR3It was speculated that miR-137 could pre-transcriptional inhibit KDELR3 through a series of bioinformatic analyses.Subsequently, we conducted multiple experiments in vitro to demonstrate above hypothesis.PCR tests confirmed that the mimics and inhibitors of miR-137 could effectively manipulate its expressions in BxPC3 and SW1990 cells (Figure15A).Meanwhile, miR-137 mimics significantly downregulated KDELR3 expressions, whereas its inhibitors increased KDELR3 expressions, indicating miR-137 could negatively regulate KDELR3 expression (Figure15B).The dualluciferase analysis manifested that the luciferase intensity was significantly reduced when cells were co-transfected with the F I G U R E 1 0 KDELR3 promotes the proliferation and regulates the cell cycle of PC cells.(A) Histological expression of KDELR3 in tumour and adjacent normal tissues through IHC detection (B, C) The transfection efficiency of sh-KDELR3 and OE-KDELR3 constructs.(D, E) The results of the colony formation assays revealed the effects of sh-KDELR3 and OE-KDELR3 on cell proliferation.(F-H) Flow cytometric analysis revealed that KDELR3 overexpression increased the proportion of cells at the S phase and decreased the proportion of cells at the G0/G1 phase.IHC, immunohistochemical staining; OE, overexpression; sh, short hairpin RNA.*p < 0.05, **p < 0.01, ***p < 0.001.et al.WT recombinant vector of KDELR3 and miR-137 mimics, confirming that miR-137 could bind to the 3′-UTR region of KDELR3 (5′-AGCAAUAA-3′, Figure 15C).Clonogenic assays revealed that miR-137 mimics significantly weakened the proliferative abilities of PC cells, but overexpression of KDELR3 could partially reverse above inhibitory effects (Figure 15D,E).Similarly, the cells transfected with miR-137 mimics also exhibited attenuated migrative and invasive performance, while overexpression of KDELR3 partially recovered the malignant capacities of PC cells (Figure 15F-I).All these findings depicted an upstream regulatory mode of KDELR3, namely miR-137 retarded PC progression through targeting KDELR3.

F
I G U R E 11 KDELR3 enhances the migration and invasion of PC cells.(A) Transwell migration assays were performed to evaluate cell migration in different groups.(B) The numbers of migrated cells in different groups.(C) Transwell invasion assays were performed to evaluate cell invasion in different groups.(D) The counts of invasive cells in different groups.OE, overexpression; PC, pancreatic cancer; sh, short hairpin RNA.*p < 0.05, **p < 0.01, ***p < 0.001.

F I G U R E 1 2 | 21 of 26 LIU
YWHAZ enhances the malignant potentials of PC cells.(A) Histological expression of YWHAZ in tumour and adjacent normal tissues through IHC detection (B, C) The transfection efficiency of sh-YWHAZ and OE-YWHAZ constructs.(D, E) The results of the colony formation assays revealed the effects of sh-YWHAZ and OE-YWHAZ on cell proliferation.(F) Transwell migration assays were performed to evaluate cell migration in different groups.(G) The numbers of migrated cells in different groups.(H) Transwell invasion assays were performed to evaluate cell invasion in different groups.(I) The counts of invasive cells in different groups.IHC, immunohistochemical staining; OE, overexpression; sh, short hairpin RNA.*p < 0.05, **p < 0.01, ***p < 0.001.et al.Several studies have reported that cancer cells in which ERS is induced can modulate immune cell function although the underlying mechanisms have not been elucidated.55ERS and UPR are associated with the TME.For example, ERS promotes the development of immune tolerance by activating immunosuppressive molecules and stimulating the accumulation and immunosuppressive functions of myeloid-derived suppressor cells.56,57In this study, the enrichment of CD8 T cells and the cytotoxic effects of T cells were significantly downregulated in the high-risk group (Figure7A,B).

F I G U R E 1 3
The effect of KDELR3 on xenograft tumour behaviour.(A, B) KDELR3 knockdown suppresses PC growth in the xenograft model.(C) Expression levels of KDELR3 in mice xenografts.C, negative control group, sh-vector; S, KDELR3 deletion group, sh-KDELR3.(D) The differential tumour weight between mice injected with shvector-transfected cells and those injected with sh-KDELR3-transfected cells.(E) The differential tumour volume between mice injected with sh-vector-transfected cells and those injected with sh-KDELR3transfected cells.PC, pancreatic cancer; sh, short hairpin.*p < 0.05, **p < 0.01, ***p < 0.001.F I G U R E 1 4 Potential upstream regulatory mechanisms of KDELR3.(A) Six candidate miRNAs were identified from three public databases.(B) Regulatory network of six predicted miRNAs.(C-H) Correlation of the expression of six predicted miRNAs with that of KDELR3.(I) Predicted binding site between miR-137 and KDELR3.(J) Thirteen candidate TFs were identified from two public databases.(K) Regulatory network of 13 predicted TFs.(L) Correlation of the expression of 13 predicted TFs with that of KDELR3.(M) Potential effects of KDELR3 on the activities of eight classic tumour-related signalling pathways based on GSEA analyses.miRNA, microRNA; TF, transcription factor.*p < 0.05; **p < 0.01. 22 Correlation between ERS risk score and metabolic reprogramming.
,D).Thus, nucleotide metabolism, but not AA metabolism and lipid metabolism, is upregulated in the high-risk group.TA B L E 2Abbreviations: ERS, endoplasmic reticulum stress; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MHC, major histocompatibility complex; NA, not available; NK, natural killer; PAAD, Pancreatic adenocarcinoma.
Functions of ERS signature genes in PAAD.
72Limited studies have examined the functions of KDELR3 in cancer.Marie et al. demonstrated that KDELR3 deletion impairs the metastasis of malignant melanoma.72This stdy, for the first time, demonstrated the cancer-promoting properties of KDELR3 in PAAD, providing novel insights into the mechanisms underlying PAAD progression.
73i et al. reported that KDELR3 was significantly upregulated in HCC cells,73which was consistent with its expression pattern in PC.This suggested that KDELR3 exerts oncogenic effects in various cancers.