- Research
- Open access
- Published:
Integrated multi-omic analysis and experiment reveals the role of endoplasmic reticulum stress in lung adenocarcinoma
BMC Medical Genomics volume 17, Article number: 12 (2024)
Abstract
Background
Lung cancer is a highly prevalent malignancy worldwide and is associated with high mortality rates. While the involvement of endoplasmic reticulum (ER) stress in the development of lung adenocarcinoma (LUAD) has been established, the underlying mechanism remains unclear.
Methods
In this study, we utilized data from The Cancer Genome Atlas (TCGA) to identify differentially expressed endoplasmic reticulum stress-related genes (ERSRGs) between LUAD and normal tissues. We performed various bioinformatics analyses to investigate the biological functions of these ERSRGs. Using LASSO analysis and multivariate stepwise regression, we constructed a novel prognostic model based on the ERSRGs. We further validated the performance of the model using two independent datasets from the Gene Expression Omnibus (GEO). Additionally, we conducted functional enrichment analysis, immune checkpoint analysis, and immune infiltration analysis and drug sensitivity analysis of LUAD patients to explore the potential biological function of the model. Furthermore, we conducted a battery of experiments to verify the expression of ERSRGs in a real-world cohort.
Results
We identified 106 ERSRGs associated with LUAD, which allowed us to classify LUAD patients into two subtypes based on gene expression differences. Using six prognostic genes (NUPR1, RHBDD2, VCP, BAK1, EIF2AK3, MBTPS2), we constructed a prognostic model that exhibited excellent predictive performance in the training dataset and was successfully validated in two independent external datasets. The risk score derived from this model emerged as an independent prognostic factor for LUAD. Confirmation of the linkage between this risk model and immune infiltration was affirmed through the utilization of Gene Set Enrichment Analysis (GSEA), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. The q-PCR results verified significant differences in the expression of prognostic genes between cancer and paracancer tissues. Notably, the protein expression of NUPR1, as determined by immunohistochemistry (IHC), exhibited an opposite pattern compared to the mRNA expression patterns.
Conclusion
This study establishes a novel prognostic model for LUAD based on six ER stress-related genes, facilitating the prediction of LUAD prognosis. Additionally, NUPR1 was identified as a potential regulator of stress in LUAD.
Introduction
Lung cancer poses a significant worldwide health issue, wherein non-small cell lung cancer (NSCLC) comprises approximately 80-85% of the reported instances [1]. According to global cancer statistics, over 2 million new cases of lung cancer are diagnosed each year [2, 3]. Among the various subtypes of NSCLC, lung adenocarcinoma (LUAD) represented approximately 55-60% of cases [4]. Despite significant advancements in immune checkpoint inhibitors and anti-angiogenesis therapies that have improved survival rates, the 5-year survival rate for patients remains around 20% [5,6,7]. While several conventional clinical models have been utilized to predict the prognosis of LUAD, the inherent heterogeneity of the disease limits their ability to provide accurate results [8]. Consequently, there is a need to develop new prognostic signatures that can enhance the prognosis assessment for LUAD patients.
The endoplasmic reticulum (ER) is a multifunctional organelle responsible for protein folding, lipid biosynthesis, and calcium storage [9, 10]. Notably, it serves as a central hub for protein quality control, enabling adaptation to adverse synthesis, external stimuli, and other detrimental events. ER stress has been implicated in the development and progression of various human malignancies, as it affects multiple cancer hallmarks [11]. External adverse factors can disrupt the integrity of the ER, leading to the accumulation of unfolded or misfolded proteins within its lumen, a condition known as ER stress. This triggers the activation of the unfolded protein response (UPR) [12, 13]. In several cancer types, overexpression of ER stress indicators has been associated with poor prognosis and clinical outcomes [14]. Wei et al. conducted a study confirming that the activation of ER stress signals plays a significant role in the initiation and progression of liver cancer [15]. Furthermore, they discovered that suppressing the ER stress response enhances cellular susceptibility to cisplatin therapy in NSCLC [16]. A recent study demonstrated that ER stress induces oral squamous cell cancer cells to secrete exosome PD-L1, leading to upregulated PD-L1 expression in macrophages and driving the polarization of M2 macrophages [17]. However, a comprehensive understanding of ER stress in LUAD, including the interplay between ER stress regulators and the tumor immune microenvironment (TIME), remains elusive.
To investigate and assess the clinical significance of ER stress in LUAD, a comprehensive analysis of endoplasmic reticulum stress-related genes (ERSRGs) was conducted in this study. Additionally, a predictive model based on ERSRGs was constructed to evaluate its prognostic value in LUAD patients. Functional enrichment analysis revealed a correlation between ERSRGs and immune infiltration. The findings of this study offer insights into the potential molecular mechanisms underlying LUAD and provide valuable prognostic information for clinical management.
Materials and methods
Data collection
The clinical information and RNA sequencing data for the bioinformatics analysis were obtained from publicly available databases, namely The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/). A total of 453 samples diagnosed with LUAD were included in the training set, ensuring the availability of complete clinical information, including survival time, survival status, age, and gender. Additionally, to further validate the findings, datasets consisting of 352 patients diagnosed with LUAD was acquired from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). These datasets comprised of two independent cohorts, namely GSE31210 (246 samples) and GSE37745 (106 samples), from which mRNA expression matrices were extracted. Subsequently, these datasets were integrated and defined as the validation sets.
Clinical sample collection
For quantitative polymerase chain reaction (Q-PCR) experiments, tissue samples were obtained from a cohort of eight patients who underwent pulmonary lobectomy at the Affiliated Cancer Hospital of Nantong University between August 2022 and April 2023. The tissue samples comprised both LUAD and paired collateral cancer specimens. Additionally, six patients diagnosed with LUAD were included for immunohistochemistry experiments, and their samples were sourced from Nantong Tumor Hospital. Prior to their participation in the study, all patients provided written informed consent. The research protocol was approved by the Ethics Committee of the Affiliated Cancer Hospital of Nantong University.
Quantitative polymerase chain reaction
Eight pairs of cancer and adjacent non-cancerous tissues were collected for qPCR. Total RNA extraction was conducted utilizing TRIzol reagent (Thermo Fisher SCIENTIFIC, USA), following the manufacturer’s instructions. Reverse transcription of mRNA was accomplished using the EvoM-MLV reverse transcription kit (Accurate Biology, China) [18]. For reverse transcription of mRNA, the HiScript III RT SuperMix for qPCR (Vazyme, China) was employed. The primers used in this study were procured from Sangon Biotech, and their specific sequences are provided in Table 1.
Western blot
Besa-2b, A549, H1299, H1975, and PC9 cell cultures (MeisenCTCC, China) were collected and underwent lysis using phenylmethylsulfonyl fluoride (1:100, Beyotime, Shanghai, China) in conjunction with cell lysis buffer. The protein was then separated using electrophoresis and transferred to a polyvinylidene fluoride (PVDF) membrane (Invitrogen, USA). Following a 2-hour blocking step with 5% skim milk, the PVDF membranes underwent an incubation period at 4 °C overnight with NUPR1 antibody (1:400, Proteintech, China). After three wash cycles, blots were incubated with horseradish peroxidase-conjugated secondary antibodies (1:1000, Proteintech, China). The determination of relative expression involved dividing the target protein band density by the density of Tubulin.
Immunohistochemistry
The tumor tissue microarray, obtained from the Affiliated Tumor Hospital of Nantong University, was utilized for the validation of the queue. Immunohistochemistry (IHC) was performed following previously established protocols [19]. In brief, the tissue sections were incubated with a primary Anti-NUPR1 antibody (dilution 1:100; catalog number 15056-1-AP, Proteintech, China) and subsequently processed using the appropriate detection system. A scanning microscope (Nikon, Japan) was employed for capturing high-resolution images of the stained sections. The evaluation of NUPR1 staining was conducted by two independent pathologists, who were blinded to the corresponding clinical information. They assessed the staining intensity, distribution, and cellular localization of NUPR1 in a semi-quantitative manner using established scoring criteria. Any discrepancies between the two pathologists were resolved through consensus discussion.
Differentially expressed genes associated with ER stress
Differential gene expression analysis was performed on the TCGA-LUAD dataset using the R-package “limma” to identify genes that were differentially expressed between LUAD and healthy samples. The criteria for differential expression were set as |log2FoldChange| > 1 and adj. p < 0.05. The resulting gene set was then intersected with the ERSRGs, leading to the identification of 106 ER stress-related differentially expressed genes (DEGs).
To further identify ER stress-related prognostic genes, univariate Cox regression analysis was conducted on the TCGA-LUAD dataset. The genes were subjected to least absolute shrinkage and selection operator (LASSO) regression using the R package “glmnet” to identify genes associated with overall survival (OS). Multiple factor stepwise regression was then applied, and the smallest lambda value was considered as the optimal value. Risk score: \( {\sum }_{\text{i}=1}^{\text{n}}=\) Coef (gene) * Expression (gene). A nomogram was constructed using the R packages “rms” and “survival” to predict the survival of LUAD patients. The nomogram included various variables such as age, sex, TN staging, histological grading, and risk score. The accuracy of the nomogram was verified by plotting calibration curves at 1-year, 3-year, and 5-year intervals using the R package “rms”.
Based on the median risk score, patients were divided into high-risk and low-risk subgroups. Kaplan-Meier (K-M) curves, time-dependent receiver operating characteristic (ROC), and a riskscore plot were generated through multivariate Cox regression analysis to illustrate the distribution and survival status of LUAD patients in the two risk groups.
Consensus clustering analysis
Cluster analysis was performed on a cohort of 106 patients using the Pearson correlation distance measure. To ensure robustness, the clustering process was repeated 10 times on 80% of the samples. The optimal number of clusters was determined by analyzing the empirical cumulative distribution function graph.
Functional enrichment analysis
Gene Ontology (GO) analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis, and Gene Set Enrichment Analysis (GSEA) were employed to investigate potential mechanisms and pathways associated with the two clusters and riskscore subgroups [20].
Immune infiltration analysis and single cell sequencing analysis
The degree of immune cell infiltration and the presence of immune checkpoints were compared across different riskscore subgroups. The “ESTIMATE” package was utilized to calculate the stromal score, immune score, and tumor purity for LUAD patients. To ensure the stability of the results, the “TIMER”, “EPIC”, and “MCP-counter” tools were employed.
Additionally, single-cell RNA sequencing data from the GSE127465 dataset was screened to identify relevant information. The TISCH platform, specifically its built-in tSNE algorithm available at http://tisch1.comp-genomics.org/, was employed for dimensionality reduction and visualization of the identified clusters.
Cell culture
Beas-2b cells, representing a human normal lung epithelial cell line, underwent cultivation in RPMI 1640 medium supplemented with 10% fetal calf serum (Gibco, Grand Island, NY, USA) and 1% Penicillin-Streptomycin (NCM Biotech, China). H1299, H1975, and PC9 cells were cultured in RPMI 1640 medium with 10% fetal bovine serum (FBS) at 37 °C. Similarly, A549 cells were nurtured in F-12 K medium (Gibco, Grand Island, NY, USA) supplemented with 10% FBS. Cell lines were meticulously maintained in their respective culture media to facilitate optimal growth and experimental conditions.
Cell counting Kit-8 assay and transwell
H1299 cells were inoculated in 96-well plates at a density of 5*103 cells per well. Subsequently, the cells were exposed to Trifluoperazine dihydrochloride (MedChemExpress, China) at a concentration of 20 μmol /mL and incubated at 37 °C for 24, 48, and 72 h. Following the respective incubation periods, an absorbance reading at 450 nm was obtained using a microplate reader (Thermo Fisher Scientific, Waltham, MA, USA). This measurement was conducted after a 2-hour incubation at 37 °C with 10 μL of Cell Counting Kit-8 (CCK-8; Bimake, Houston, TX, USA) reagent in each well. To standardize the concentration of H1299 cells to 1 × 105/mL, serum-free RPMI-1640 medium was employed. The upper chamber received 200 μL of RPMI-1640 medium containing 10% fetal bovine serum, while the lower chamber was supplemented with 600 μL of RPMI-1640 medium containing 10% fetal bovine serum. The cell culture was maintained at 37℃ with 5% CO2 for 24 h. After chamber removal, the cells in the upper chamber were delicately swabbed with a cotton applicator, fixed with 4% paraformaldehyde for 15 min, stained with crystal violet at room temperature for 25 min, and excess staining solution was rinsed off with PBS. Observation of cells that traversed the membrane was performed under a microscope. For each sample, three randomly selected fields of view were photographed and enumerated, and the average value was computed.
Drug sensitivity analysis
The half-maximum inhibitory concentration (IC50) metrics for chemotherapeutic agents were acquired from the Genomics of Cancer Drug Sensitivity (GDSC) database, accessible at https://www.cancerrxgene.org/. Subsequently, “PRrophytic” R package was performed to calculate the drug susceptibility between different subgroups in R software. The outcomes are visually represented through box plots.
Statistical analysis
Statistical analyses and data visualization were conducted using R software version 4.1.0 and GraphPad Prism version 9.5.1. Both univariate and multivariate Cox regression analyses were employed to examine the impact of various factors on the prognosis of LUAD. Statistical significance was defined as p-values < 0.05.
Results
The screening and characterization of ERSRGs in LUAD
The study flow-chart is shown in Fig. 1. Differential expression analysis between the LUAD and control groups (TCGA cohort) was conducted using the limma package in R, leading to the identification of 20,724 DEGs. Among these DEGs, 14,019 were up-regulated, while 6,705 were down-regulated (Fig. 2A). To investigate the potential involvement of ER stress in LUAD, a Venn analysis was performed to identify the overlap between the LUAD-related DEGs and ER stress genes (Fig. 2B). The interaction network among the 106 genes is visualized in Fig. 2C, and the intersecting genes were functionally annotated using GO and KEGG analyses. As depicted in Fig. 2D, these genes are enriched in various biological processes, such as response to ER stress, response to topologically incorrect protein, ER-associated Degradation (ERAD) pathway, cell components including ER protein-containing complex, integral component of ER membrane, ER ubiquitin ligase complex, and molecular functions such as ubiquitin-like protein ligase binding, ubiquitin ligase binding and ubiquitin protease binding. The critical functions of the ERSRGs include ubiquitin mediated proteolysis, protein processing in ER, B cell receptor signaling pathway and so on (Fig. 2E).
Consensus clustering analysis of ERSRGs in LUAD
Based on the identified set of 106 genes, the cohort of LUAD patients from the training queue (TCGA cohort) was subjected to Cluster Analysis to classify them into two subgroups. Optimal stability was observed when K = 2, resulting in the classification of 232 patients into cluster 1 and 221 patients into cluster 2 (Fig. 3A-C). Principal component analysis demonstrated a distinct separation of samples into two clusters (Fig. 3D). Notably, the OS rate of patients in cluster 2 was significantly higher than that of cluster 1 (P = 0.03; Fig. 3E). These findings supported the subdivision of LUAD patients into two distinct molecular subtypes associated with differing survival outcomes. Furthermore, a volcano plot depicting the logFC and FDR values of 620 upregulated genes and 321 downregulated genes across the two clusters was generated (Fig. 3F). Subsequent GO enrichment analysis revealed that these genes were significantly enriched in specific molecular processes, such as mitotic sister chromatid segregation, humoral immune response, regulation of humoral immune response, condensed chromosome, centromeric region, and others (Fig. 3G; Table 2). Additionally, GSEA revealed that the enriched pathways were predominantly associated with immune infiltration (Fig. 3H). Consequently, it is reasonable to hypothesize that the influence of risk scores may impact the prognosis of LUAD by modulating the immune microenvironment.
Prognostic signature was constructed and validation based on ERSRGs in LUAD
The initial step in developing a prognostic model involved the identification of candidate prognostic ERSRGs through univariate Cox regression analysis. As depicted in Fig. 4A, the OS of LUAD patients exhibited a significant correlation with 18 ERSRGs. Subsequently, LASSO analysis was employed to detect and screen 15 DEGs associated with ER stress (Fig. 4B). Furthermore, a multiple factor stepwise regression analysis was performed, resulting in the selection of 6 genes for constructing the prognostic model (Fig. 4C). The risk scoring model was established using the following formula: Riskscore = (-0.894871321) × EIF2AK3 + (0.268950582) × BAK1 + (-0.127961954) × NUPR1 + (0.255795681) × VCP + (0.393370303) × MBTPS2 + (-0.292131004) × RHBDD2. Based on the median risk score, patients from the TCGA-LUAD cohort were stratified into a high-risk subgroup (n = 220) and a low-risk subgroup (n = 220) to facilitate further investigations. Subsequent K-M analysis and risk survival status plot revealed that the high-risk subgroup exhibited a worse prognosis, whereas the low-risk subgroup demonstrated prolonged survival (Fig. 4D and F). The prognostic models were assessed by calculating the area under the curve (AUC) for 1-year, 3-year, and 5-year survival, yielding values of 0.68, 0.69, and 0.70, respectively (Fig. 4E).
Assessment and external validation for ERSRGs-signature
The risk distribution curve, survival status, and expression heatmap of the external validation sets (GSE37745 and GSE30210) demonstrated that patients with low-risk scores exhibited significantly longer survival times compared to those with high-risk scores, thus validating the findings from the training set (Fig. 5A and B). To further consolidate the prognostic model, the clinical information and genetic characteristics from TCGA were integrated, and a comprehensive multi-factor Cox regression model was developed, resulting in the construction of a nomogram (Fig. 5C). Calibration plots were employed to assess the predictive accuracy of the nomogram, revealing excellent agreement between the predicted and observed OS rates at 1, 3, and 5 years (Fig. 5D). Moreover, the nomogram model was subjected to decision curve analysis (DCA) to evaluate its clinical utility and potential benefits (Fig. 5E-G). Collectively, the risk score, when combined with the ERSRGs-signature, pathological stage, and N-stage, emerged as an independent and robust prognostic indicator, providing enhanced prognostic value for patients with LUAD.
Exploring immune infiltration patterns and single-cell analysis of ERSRGs-signature in LUAD
To unravel the potential functions and pathways associated with prognostic features, we conducted comprehensive enrichment analyses of Gene Set, GO, and KEGG pathways. The results revealed that the genes linked to prognostic features were predominantly enriched in pathways related to immunoinfiltration. Hence, we proceeded to explore the heterogeneity of immune microenvironments among ERSRGs-signature (Fig. 6A-C). Initially, we assessed the correlation between gene expression and immune infiltration in LUAD and observed significant variations in the expression of different genes across immune cells (Fig. 6D). Subsequently, we employed the TIMER and EPIC algorithms to investigate immune infiltration patterns between the low and high-risk subgroups. The low-risk subgroup exhibited significantly elevated expression of B cells, CD4 T cells, CD8 T cells, and macrophage cells compared to the high-risk group (Fig. 6E and F). To validate the stability and robustness of these findings, we utilized additional algorithms, namely MCP-counter and ESTIMATE, which yielded consistent results (Fig. 6G and H). Furthermore, we observed substantial differences in the expression of immune checkpoints between the two subgroups (Fig. 6I).
Subsequently, we conducted single-cell sequencing analysis of the ERSRGs-signature. Cluster analysis was performed, and Fig. 7A depicted the cluster display using t-distributed stochastic neighbor embedding (tSNE), where each color represented a distinct cell type identified within the clusters. Each cell was represented by a scatter plot, and the numbers in the figure corresponded to the cluster numbers. It was evident that there are 25 distinct cell populations. Figure 7B presented the annotation of clusters based on marker analysis, revealing significant differences in gene expression among different immune cells. After applying tSNE dimensionality reduction, the mRNA distribution of BAK1, EIF2AK3, MBTPS2, NUPR1, RHBDD2, and VCP was shown in Fig. 7C-H. Finally, we analyzed the differential expression of ERSRGs in the various immune cell clusters. Among them, BAK1 exhibited the lowest expression in immune cells, while VCP demonstrated the highest expression (Fig. 7I).
Overall, the riskscore demonstrated an inverse correlation with the level of immune infiltration, providing novel insights into the relationship between ERSRGs and the immune status of LUAD.
Validation of the expression levels of ERSRGs in LUAD
To further investigate the association between the prognostic ERSRGs-signature and LUAD, in vitro experiments were conducted using qPCR analysis on peritumoral and tumor tissues. The findings revealed a significant upregulation of BAK1 and EIF2AK3 expression in LUAD tissues, whereas NUPR1, RHBDD2, and VCP exhibited the opposite trend (Fig. 8A-G). Moreover, the Human Protein Atlas (HPA) database analysis showed higher expression levels of BAK1A and EIF2AK3 in LUAD tissues compared to normal tissues (Fig. 8G). However, NUPR1 data was unavailable in the HPA database. Therefore, to explore the protein expression of NUPR1 in LUAD patients, IHC analysis was performed at Nantong Cancer Hospital. Interestingly, the protein expression of NUPR1, as determined by IHC, exhibited an opposite pattern compared to the mRNA expression patterns (Fig. 9A).
Validation of NUPR1 under experiment
In this study, we implemented a comprehensive validation of NUPR1 within authentic laboratory conditions. Initially, IHC analysis was conducted on pathological specimens obtained from 6 LUAD patients. The results revealed a conspicuous aggregation of NUPR1 within cancerous tissue compared to adjacent non-cancerous tissues (Fig. 9A and B). Subsequently, both RNA and protein expression levels of NUPR1 were scrutinized in normal lung epithelial cells and four distinct LUAD cell lines. Surprisingly, NUPR1 RNA exhibited its highest expression in normal cell lines (Fig. 9C), aligning with our bioinformatics analysis outcomes. In contrast, NUPR1 protein displayed heightened expression levels in LUAD cells (Fig. 9D and E, Fig. S2 and S3). We postulated that potential post-translational modifications may underlie this incongruity. To gain deeper insights into the functional role of NUPR1 in LUAD progression, we procured NUPR1 inhibitors and executed cell proliferation and transwell experiments. The results starkly indicated that upon NUPR1 inhibition, both cell proliferation and invasive capacity were markedly attenuated (Fig. 9F and G). This unequivocally underscores the contributory role of NUPR1 protein in the advancement of LUAD.
Correlation between risk score and IC50 values for therapeutic agents
The impact of risk scores on the IC50 values of a set of 30 distinct drug molecules was systematically assessed to discern their therapeutic efficacy. Except for BI-2536 and WIKI4, all other drugs exhibited higher resistance in the high-risk group (Fig. 10 and Fig. S1). This observation underscores the potential utility of our prognostic model in guiding the use of therapeutic agents.
Discussion
LUAD represents the most prevalent subtype of lung cancer, a grave malignancy arising from the accumulation of various genetic mutations. These mutations lead to uncontrolled proliferation of lung cells and the subsequent formation of tumors. Upon recognition by the immune system, these transformed cancer cells elicit an immune response aimed at their elimination [21]. Nonetheless, immune escape not only expedites tumor progression but also impairs the efficacy of cancer immunotherapy [22, 23]. The ER pathway serves as a critical regulator of ER homeostasis. Disruption of ER function triggers a phenomenon referred to as “ER stress” [24]. In the context of tumorigenesis, the rapid proliferation rate of cancer cells necessitates heightened activity of ER protein folding, assembly, and transport, thereby inducing physiological stress within the ER [25]. The ER stress response is believed to confer cellular protection and is implicated in tumor growth and adaptation to challenging environments [26]. Sustained ER stress represents a novel characteristic of cancer, resulting from various metabolic and carcinogenic abnormalities that disrupt protein-folding homeostasis in aggressive immune cells. Constitutive activation of the ER stress response enables malignant cells to adapt to carcinogenesis and environmental stressors by coordinating multiple immune regulatory mechanisms and promoting malignant progression concurrently [27]. Nonetheless, the precise relationship between ER stress and the immune microenvironment remains inadequately investigated.
In our study, we initially screened 106 genes associated with ER stress to identify differential expression patterns between cancer and para-cancer samples. K-Medoids clustering was employed for this purpose. The differential genes in the two resulting clusters were primarily enriched in processes related to the adaptive immune system, humoral immune response, and regulation of humoral immune response. Notably, patients belonging to cluster 1 exhibited a significantly longer survival time compared to those in cluster 2. This discrepancy in prognosis suggests a potential correlation with immune response. Through a series of statistical analyses, including univariate regression, LASSO, and logistic stepwise regression, we identified 6 key ERSRGs. Subsequently, we constructed a novel prognostic risk spectrum based on the expression signature of these six genes (referred to as ERSRGs). This risk spectrum allowed us to classify patients with LUAD into distinct risk subgroups, based on their respective median risk scores. Importantly, a higher risk score was associated with worse prognosis for the patients.
The prognostic features of interest encompass 6 ERSRGs, specifically EIF2AK3, MBTPS2, RHBDD2, VCP, NUPR1, and BAK1. Among these, EIF2AK3, NUPR1, and RHBDD2 demonstrated protective characteristics, while MBTPS2, VCP, and BAK1 were strongly associated with poor prognosis. To assess their expression levels, qPCR analyses were conducted on cancer and para-cancer samples from 8 patients diagnosed with LUAD. The results revealed significant differential expression of EIF2AK3, RHBDD2, VCP, NUPR1, and BAK1, with NUPR1 and RHBDD2 exhibiting the most pronounced differences. EIF2AK3 has been identified as an immune-related prognostic gene in breast cancer, exerting a role in tumor cell apoptosis and facilitating sustained protective antitumor immunity [28]. MBTPS2, a membrane-embedded zinc metalloprotease, activates signaling proteins involved in transcriptional control of sterol and the ER stress response [29], thus promoting the progression of prostate cancer [30] and colorectal cancer [30]. The RHBDD2 (Rhomboid domain containing 2) gene is found to be overexpressed in advanced stages of colorectal cancer (CRC) and potentially modulates the UPR pathway, thereby favoring cell migration, adhesion, and proliferation [31]. VCP (valosin-containing protein) is crucial for maintaining mitochondrial function, and in prostate cancer cells, it employs self-aggregation to inhibit mitochondrial activity, thereby evading cell death during nutrient deprivation and promoting malignancy [32]. In a cohort study, Tao et al. demonstrated that NUPR1 serves as a protective factor in the survival prognosis of LUAD [33], while Li et al. suggested NUPR1 to be a potential risk gene [34]. NUPR1, a nuclear protein, plays a critical role in redox reactions [35], and macrophages have been implicated as the most relevant immune cells associated with NUPR1 expression in bladder cancer [36]. Furthermore, the mechanism through which BAK1 promotes cisplatin resistance in NSCLC is believed to involve the inhibition of cell apoptosis [37]. In summary, all 6 identified genes contribute to tumor development and progression by modulating pathways associated with tumor metabolism, with NUPR1 considered particularly significant.
Nuclear Protein 1 (NUPR1) is a small, highly basic transcriptional regulator involved in the regulation of diverse cellular processes, such as DNA repair, ER stress, and oxidative stress response. The cellular localization of NUPR1 appears to be associated with pathological conditions. Prominent cytoplasmic staining has been observed in large papillary tumors, tumors exhibiting lymph node metastasis, and NSCLC [38]. Our IHC analysis corroborated these findings. However, intriguingly, our real-world cohort study revealed that, in contrast to mRNA expression, NUPR1 accumulates in cancerous tissues, contributing to the malignant progression of cancer, which necessitates further investigation. Garcia Montero et al. reported that under various stress conditions, NUPR1 mRNA expression was rapidly, strongly, and transiently stimulated [39]. Cancer cells endure and adapt to various types of stressful environments over prolonged periods [40], leading us to speculate that NUPR1 mRNA may be consumed more in cancerous tissues compared to adjacent tissues. Additionally, interestingly, the protein expression of NUPR1 has been shown to positively correlate with cell density [41]. Considering that cancer arises from unregulated and excessive cell division and proliferation, resulting in higher cell density [42], we hypothesize that NUPR1 expression is relatively elevated in cancer cells characterized by higher cell density compared to adjacent cells with relatively fewer cells.
To verify the broad applicability of the risk assessment element group, we conducted validation using external datasets GSE31210 and GSE37745. The signature exhibited robust predictive performance not only in the internal dataset but also in the validation sets. Evidence from ROC curves and K-M analysis demonstrated the remarkable predictive effect of the ERSRGs on the prognosis of LUAD patients. Importantly, even after stratifying clinical features, this signature remained significantly prognostic in LUAD patients. Therefore, we propose that ER stress-related features possess excellent predictive performance for OS and could serve as independent prognostic indicators for LUAD. To facilitate clinical application, we constructed a nomogram model and verified its accuracy using calibration diagrams.
Previous research has highlighted the role of ER stress in promoting immune escape and facilitating metastasis [43, 44]. Subsequent GSEA, GO, and KEGG analyses of the two subgroups revealed enrichment in immune-related pathways. Notably, tumor purity has been identified as negatively correlated with immune response, suggesting its potential as an indicator of the immune response level in the tumor microenvironment [45]. To explore this further, we employed four different immune scoring algorithms, and all results consistently indicated that individuals classified as low-risk exhibited higher expression levels of B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and endothelial cells. The density of CD8+ T cells and mature dendritic cells has been closely associated with the survival rate of lung cancers, with higher CD8+ T cell density correlating with better 5-year survival rates [46], consistent with our findings. Additionally, we observed decreased expression of immune checkpoint genes in the high-risk group, which may be attributed to immune cell dysregulation. Therefore, our new prognostic model holds potential to not only assess the survival prognosis of LUAD but also shed light on the immune microenvironment.
Several limitations should be acknowledged in this study. Firstly, the model primarily relies on data from the TCGA database and the Nantong cohort, thus its generalizability to other datasets may be limited. Therefore, a prospective multicenter cohort study is necessary to validate the findings and ensure their applicability to diverse populations. Secondly, in order to comprehensively elucidate the underlying reasons for the discordance between NUPR1 mRNA and protein expression levels, further evidence from additional experiments and investigations is required.
Overall, this study presents a prognostic model based on six genes associated with ER stress. The model exhibits utility in predicting the survival outcomes of patients with LUAD and offers insights into tumor immune infiltration to some extent. Furthermore, the identification of key genes provides novel insights into the molecular mechanisms underlying LUAD.
Data Availability
The datasets generated during and/or analyzed during the current study are available in the TCGA (http://www.genome.ucsc.edu/) and GEO (https://www.ncbi.nlm.nih.gov/geo/) Databases, with specific accession numbers GSE31210 and GSE37745. The original contributions presented in the study are included in the article/Additional files, and further inquiries can be directed to the corresponding authors.
References
Cao W, Chen HD, Yu YW, Li N, Chen WQ. Changing profiles of cancer burden worldwide and in China: a secondary analysis of the global cancer statistics 2020. Chin Med J. 2021;134(7):783–91.
Pine SR. Rethinking gamma-secretase inhibitors for treatment of non-small-cell lung cancer: is notch the target? Clin cancer Research: Official J Am Association Cancer Res. 2018;24(24):6136–41.
Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Cancer J Clin. 2021;71(3):209–49.
Wang N, Chen WQ, Zhu WX, Xing XM, Lu AP, Yang L. [Incidence trends and pathological characteristics of lung cancer in urban Beijing during period of 1998–2007]. Chin J Prev Med. 2011;45(3):249–54.
Aisner DL, Sholl LM, Berry LD, Rossi MR, Chen H, Fujimoto J, Moreira AL, Ramalingam SS, Villaruz LC, Otterson GA, et al. The impact of smoking and TP53 mutations in lung adenocarcinoma patients with targetable mutations-the lung cancer mutation consortium (LCMC2). Clin cancer Research: Official J Am Association Cancer Res. 2018;24(5):1038–47.
Lazzari C, Karachaliou N, Bulotta A, Viganó M, Mirabile A, Brioschi E, Santarpia M, Gianni L, Rosell R, Gregorc V. Combination of immunotherapy with chemotherapy and radiotherapy in lung cancer: is this the beginning of the end for cancer? Therapeutic Adv Med Oncol. 2018;10:1758835918762094.
de Langen AJ, Smit EF. Therapeutic approach to treating patients with BRAF-mutant lung cancer: latest evidence and clinical implications. Therapeutic Adv Med Oncol. 2017;9(1):46–58.
Zuo S, Wei M, Zhang H, Chen A, Wu J, Wei J, Dong J. A robust six-gene prognostic signature for prediction of both disease-free and overall survival in non-small cell lung cancer. J Translational Med. 2019;17(1):152.
Wang M, Kaufman RJ. Protein misfolding in the endoplasmic reticulum as a conduit to human disease. Nature. 2016;529(7586):326–35.
Schwarz DS, Blower MD. The endoplasmic reticulum: structure, function and response to cellular signaling. Cell Mol Life Sci. 2016;73(1):79–94.
Pavlović N, Heindryckx F. Exploring the role of endoplasmic reticulum stress in hepatocellular carcinoma through mining of the human protein atlas. Biology. 2021;10(7).
Wang S, Kaufman RJ. The impact of the unfolded protein response on human disease. J Cell Biol. 2012;197(7):857–67.
Ron D, Walter P. Signal integration in the endoplasmic reticulum unfolded protein response. Nat Rev Mol Cell Biol. 2007;8(7):519–29.
Dalton LE, Clarke HJ, Knight J, Lawson MH, Wason J, Lomas DA, Howat WJ, Rintoul RC, Rassl DM, Marciniak SJ. The endoplasmic reticulum stress marker CHOP predicts survival in malignant mesothelioma. Br J Cancer. 2013;108(6):1340–7.
Lebeaupin C, Vallée D, Hazari Y, Hetz C, Chevet E, Bailly-Maitre B. Endoplasmic reticulum stress signalling and the pathogenesis of non-alcoholic fatty liver disease. J Hepatol. 2018;69(4):927–47.
Lai ST, Wang Y, Peng F. Astragaloside IV sensitizes non-small cell lung cancer cells to cisplatin by suppressing endoplasmic reticulum stress and autophagy. J Thorac Disease. 2020;12(7):3715–24.
Wang L, Liu Y, Zhang X, Ye Y, Xiong X, Zhang S, Gu L, Jian Z, Wang H. Endoplasmic reticulum stress and the unfolded protein response in cerebral ischemia/reperfusion injury. Front Cell Neurosci. 2022;16:864426.
Liu Y, Chen Z, Lin W, Zhou Y, Liu Z, Zhao R, Chen Y, Wu B, Chen A, Lin C. Role of hippocampal circKcnk9 in visceral hypersensitivity and anxiety comorbidity of irritable bowel syndrome. Front Cell Neurosci. 2022;16:1010107.
Liu Y, Lin W, Yang Y, Shao J, Zhao H, Wang G, Shen A. Role of cuproptosis-related gene in lung adenocarcinoma. Front Oncol. 2022;12:1080985.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Spella M, Stathopoulos GT. Immune resistance in lung adenocarcinoma. Cancers. 2021;13(3).
Chong C, Müller M, Pak H, Harnett D, Huber F, Grun D, Leleu M, Auger A, Arnaud M, Stevenson BJ, et al. Integrated proteogenomic deep sequencing and analytics accurately identify non-canonical peptides in tumor immunopeptidomes. Nat Commun. 2020;11(1):1293.
McGranahan N, Rosenthal R, Hiley CT, Rowan AJ, Watkins TBK, Wilson GA, Birkbak NJ, Veeriah S, Van Loo P, Herrero J, et al. Allele-specific HLA loss and immune escape in lung cancer evolution. Cell. 2017;171(6):1259–1271e1211.
Wang WA, Groenendyk J, Michalak M. Endoplasmic reticulum stress associated responses in cancer. Biochim Biophys Acta. 2014;1843(10):2143–9.
Lee AS. GRP78 induction in cancer: therapeutic and prognostic implications. Cancer Res. 2007;67(8):3496–9.
Yadav RK, Chae SW, Kim HR, Chae HJ. Endoplasmic reticulum stress and cancer. J cancer Prev. 2014;19(2):75–88.
Mogilenko DA, Haas JT, L’Homme L, Fleury S, Quemener S, Levavasseur M, Becquart C, Wartelle J, Bogomolova A, Pineau L, et al. Metabolic and innate Immune cues merge into a specific inflammatory response via the UPR. Cell. 2019;177(5):1201–1216e1219.
Wang X, Huang H, Liu X, Li J, Wang L, Li L, Li Y, Han T. Immunogenic cell death-related classifications in breast cancer identify precise immunotherapy biomarkers and enable prognostic stratification. Front Genet. 2022;13:1052720.
Oeffner F, Fischer G, Happle R, König A, Betz RC, Bornholdt D, Neidel U, Boente Mdel C, Redler S, Romero-Gomez J, et al. IFAP syndrome is caused by deficiency in MBTPS2, an intramembrane zinc metalloprotease essential for cholesterol homeostasis and ER stress response. Am J Hum Genet. 2009;84(4):459–67.
Tibbo AJ, Hartley A, Vasan R, Shaw R, Galbraith L, Mui E, Leung HY, Ahmad I. MBTPS2 acts as a regulator of lipogenesis and cholesterol synthesis through SREBP signalling in prostate cancer. Br J Cancer. 2023;128(11):1991–9.
Palma S, Raffa CI, Garcia-Fabiani MB, Ferretti VA, Zwenger A, Perez Verdera PV, Llontop A, Rojas Bilbao E, Cuartero V, Abba MC, et al. RHBDD2 overexpression promotes a chemoresistant and invasive phenotype to rectal cancer tumors via modulating UPR and focal adhesion genes. Biochim et Biophys acta Mol Basis Disease. 2020;1866(8):165810.
Ogor P, Yoshida T, Koike M, Kakizuka A. VCP relocalization limits mitochondrial activity, GSH depletion and ferroptosis during starvation in PC3 prostate cancer cells. Genes to Cells: Devoted to Molecular & Cellular Mechanisms. 2021;26(8):570–82.
Shu L, Liu S, Tao Y. Development and validation of a prognosis prediction model based on 18 endoplasmic reticulum stress-related genes for patients with lung adenocarcinoma. Front Oncol. 2022;12:902353.
Li F, Niu Y, Zhao W, Yan C, Qi Y. Construction and validation of a prognostic model for lung adenocarcinoma based on endoplasmic reticulum stress-related genes. Sci Rep. 2022;12(1):19857.
Zhan Y, Zhang Z, Liu Y, Fang Y, Xie Y, Zheng Y, Li G, Liang L, Ding Y. NUPR1 contributes to radiation resistance by maintaining ROS homeostasis via AhR/CYP signal axis in hepatocellular carcinoma. BMC Med. 2022;20(1):365.
Zhang L, Gao S, Shi X, Chen Y, Wei S, Mi Y, Zuo L, Qi C. NUPR1 imparts oncogenic potential in bladder cancer. Cancer Med. 2023;12(6):7149–63.
Wang H, Huang H, Wang L, Liu Y, Wang M, Zhao S, Lu G, Kang X. Cancer-associated fibroblasts secreted miR-103a-3p suppresses apoptosis and promotes cisplatin resistance in non-small cell lung cancer. Aging. 2021;13(10):14456–68.
Liu S, Costa M. The role of NUPR1 in response to stress and cancer development. Toxicol Appl Pharmcol. 2022;454:116244.
Garcia-Montero A, Vasseur S, Mallo GV, Soubeyran P, Dagorn JC, Iovanna JL. Expression of the stress-induced p8 mRNA is transiently activated after culture medium change. Eur J Cell Biol. 2001;80(11):720–5.
O’Malley J, Kumar R, Inigo J, Yadava N, Chandra D. Mitochondrial stress response and cancer. Trends in cancer. 2020;6(8):688–701.
Huang C, Santofimia-Castaño P, Liu X, Xia Y, Peng L, Gotorbe C, Neira JL, Tang D, Pouyssegur J, Iovanna J. NUPR1 inhibitor ZZW-115 induces ferroptosis in a mitochondria-dependent manner. Cell Death Discovery. 2021;7(1):269.
Patra D, Bhavya K, Ramprasad P, Kalia M, Pal D. Anti-cancer drug molecules targeting cancer cell cycle and proliferation. Adv Protein Chem Struct Biology. 2023;135:343–95.
Tanaka A, Sakaguchi S. Regulatory T cells in cancer immunotherapy. Cell Res. 2017;27(1):109–18.
Alifano M, Mansuet-Lupo A, Lococo F, Roche N, Bobbio A, Canny E, Schussler O, Dermine H, Régnard JF, Burroni B, et al. Systemic inflammation, nutritional status and tumor immune microenvironment determine outcome of resected non-small cell lung cancer. PLoS ONE. 2014;9(9):e106914.
Liu D, Schilling B, Liu D, Sucker A, Livingstone E, Jerby-Arnon L, Zimmer L, Gutzmer R, Satzger I, Loquai C, et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med. 2019;25(12):1916–27.
Marty R, Kaabinejadian S, Rossell D, Slifker MJ, van de Haar J, Engin HB, de Prisco N, Ideker T, Hildebrand WH, Font-Burgada J, et al. MHC-I genotype restricts the oncogenic mutational landscape. Cell. 2017;171(6):1272–1283e1215.
Acknowledgements
All authors sincerely thank the National Biotechnology Information Center (NCBI) for its contribution to the creation and distribution of the GEO database.
Funding
This study was supported by Nantong Municipal Health Commission scientific research project (MS22022118).
Author information
Authors and Affiliations
Contributions
W.G. funding the work. Initial draft was completed by L.Y. and L.Y. conducted all experiments and carried out the data analysis. Q.H., Y.Y., Z.X., W.C., P.X. collected the clinical data. L.W. directed and critically revised the article. All authors approved the final version for submission.
Corresponding authors
Ethics declarations
Ethics approval and consent to participate
The study was conducted in accordance with the Declaration of Helsinki and approved by the Institutional Review Board of Affiliated Tumor Hospital of Nantong University (2022-No.063) for studies involving humans. All patients in the study signed written informed consent forms.
Consent for publication
Not applicable.
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
Supplementary Material 1: Fig. S1
(A-N) Therapeutic drugs showed significant IC50 differences in high- and low-risk groups. Fig. S2 Original, unprocessed versions for WB: NUPR1. Fig. S3 Original, unprocessed versions for WB: Tubulin
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Liu, Y., Lin, W., Qian, H. et al. Integrated multi-omic analysis and experiment reveals the role of endoplasmic reticulum stress in lung adenocarcinoma. BMC Med Genomics 17, 12 (2024). https://doi.org/10.1186/s12920-023-01785-4
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12920-023-01785-4