Novel endoplasmic reticulum stress-related gene signature unveils CDKN3 as a prognosticator in neuroblastoma

Introduction

Neuroblastoma (NB) is the predominant extracranial solid tumor in children, constituting 8% to 10% of all childhood malignancies and accounting for approximately 15% of childhood cancer deaths (1). NB exhibits significant heterogeneity, with biological and clinical manifestations ranging from spontaneous regression to refractory disease (2). Moreover, genetic factors, such as polymorphisms in ALKBH5 (3), TET2 (4), have been demonstrated to significantly affect both the risk and clinical prognosis of NB. While patients with low to intermediate-risk NB generally have a favorable prognosis, those diagnosed with high-risk NB frequently present with metastases at diagnosis, missing the critical window for optimal treatment (5). Despite recent advancements enhancing treatment outcomes for this high-risk group, the five-year survival rate remains below 50% (6). Currently, there is a lack of global consensus on the staging and typing of NB, with standards varying significantly across different regions and among diverse NB collaborative groups. This highlights the importance of identifying reliable biomarkers for prognosis prediction and exploring targeted approaches for precision therapy.

For NB patients, MYCN amplification is an important clinical biomarker associated with poor prognosis (7). This genetic alteration is widely acknowledged driver of high-risk NB (8). The MYCN protein promotes cellular proliferation and growth by regulating the expression of downstream genes while concurrently inhibiting cell differentiation and apoptosis, thereby facilitating tumor progression (9). Recent studies have shown that MYCN influences and regulates endoplasmic reticulum (ER) stress, which impacts multiple facets of cancer (10-12).

ER plays a critical role in protein synthesis, folding, and transport in eukaryotes, accounting for at least one-third of these activities and regulating protein homeostasis (13). To ensure the correct folding and assembly of proteins, multiple regulatory mechanisms are in place within the ER (14). Factors such as inflammation, hypoxia, and malnutrition impair ER homeostasis, resulting in the accumulation of improperly folded or misfolded proteins in the ER and inducing ER stress (15). ER stress is primarily managed through the activation of three unfolded protein response (UPR) signaling pathways: IRE1α, PERK, and ATF6 (16,17). Activation of the UPR serves as a pivotal cellular mechanism that mitigates the misfolding and aggregation of proteins. It can precipitate a cascade of events culminating in apoptosis or the emergence of other detrimental pathological alterations, allowing cells to adapt to changes in their internal and external environments. Once, the ER stress is too severe or persists for too long, it can lead to apoptosis or other pathological changes (18,19). Many studies have demonstrated that continuous ER stress plays a role in the onset and progression of various tumors (20). Additionally, genes related to ER stress have been identified as central to cancer malignancy (21). However, the role of ER stress in the occurrence and progression of NB is yet to be understood.

In this study, we identified differentially expressed ER stress-related genes linked to MYCN amplification. Based on these findings, we developed a prognostic model comprising five genes that accurately predicts NB outcomes. We conducted a comprehensive analysis assessing prognosis, clinical characteristics, tumor microenvironment (TME), and drug response across various patient subgroups. Additionally, we performed an initial validation of a critical gene in the risk model using clinical samples. Ultimately, we constructed a prognostic nomogram incorporating this gene signature and relevant clinical factors. This nomogram is designed to assist physicians in crafting tailored treatment strategies for distinct NB patient subgroups. We present this article in accordance with the TRIPOD reporting checklist (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-142/rc).

Methods

Data acquisition

Gene expression datasets and corresponding clinical information for NB were retrieved from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) under accession number GSE49710 (n=498) and from ArrayExpress (https://www.ebi.ac.uk/biostudies/arrayexpress) under accession number E-MTAB-8248 (n=223). A total of 721 NB samples from these datasets were randomly divided into a training set and a test set, following a 7:3 ratio (505 samples for training, 216 samples for testing).

Acquisition of differentially expressed ER stress-related genes

ER stress-related genes were extracted from the GeneCards database (https://www.genecards.org/) with a relevance score of ≥7, yielding a total of 785 genes for subsequent analysis (available online: https://cdn.amegroups.cn/static/public/tp-2025-142-1.xlsx). To identify ER stress-related differentially expressed genes (DEGs) in MYCN amplified and non-amplified patients, all ER stress-related genes were explored between two groups using the “limma” R package (22), applying screening criteria of P<0.05 and |log₂fold change (FC)| >1. The results were visualized using a volcano plot created with the “ggplot2” R package. Subsequently, we performed clustering analysis on ER stress-related DEGs, categorizing them based on MYCN amplification and non-amplification, and the outcomes were visualized using a heat map plotted with the “heatmap” R package.

Construction and validation of ER stress-related signature

The training set comprising 505 samples was utilized for univariate Cox regression analysis based on ER stress-related DEGs. Next, the candidate model genes were incorporated into the least absolute shrinkage and selection operator (LASSO) regression, aimed at minimizing the potential overfitting impact of the signature. Multivariate analysis was then used to identify hub prognostic genes and develop an ER stress-related gene signature. The formula was as follows: Riskscore=∑i=1nCoefi×ERSRGi. In this formula, Coef represents the coefficient derived from Cox analyses. The samples were divided into high- and low-risk groups according to the median risk score. The predictive performance of the signature was evaluated using Kaplan-Meier (KM) analysis in both the training and testing sets. Receiver operating characteristic (ROC) curves were used to assess the accuracy of the signature.

Enrichment of differentially expressed genes by gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and gene set enrichment analysis (GSEA)

Functional enrichment analysis of ER stress-related DEGs, including GO, KEGG, and GSEA was conducted using the “clusterProfiler” R package to explore the biological functions and pathways associated with these genes (23).

Immune infiltration analysis

We employed bioinformatics tools, including Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) (24), xCell (25), Microenvironment Cell Population counter (MCPcounter) (26), and Estimating the Proportions of Immune and Cancer cells (EPIC) (27), to analyze immune infiltration in NB patients. The Wilcoxon rank-sum test was used to assess differences in immune infiltration between high- and low-risk groups, while Spearman correlation analysis was conducted to determine the correlation coefficients of immune infiltration. Additionally, ssGSEA was utilized to calculate scores for infiltrating immune cells and to evaluate the activity of immune-related pathways. The Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm was applied to estimate the abundance of stromal and immune cells in the TME by analyzing gene expression in NB tissues (28).

Potential therapeutic drugs

Data on anticancer drugs were obtained from the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) (29). Using the ‘pRRophetic’ package in R software (30), sensitivity scores for various drug groups were calculated and differences in drug sensitivity between high-risk and low-risk groups were analyzed.

Immunohistochemistry (IHC)

This study received ethics approval from the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (approval No. E20210027) and was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Written informed consent was obtained from the legal guardians or next of kin of all participants. Forty-six paraffin-embedded NB samples were obtained from patients at Tianjin Medical University Cancer Institute and Hospital, with related clinicopathological and follow-up information also collected. Tissue slides were prepared, heated, and dehydrated. After antigen retrieval, endogenous peroxidases were quenched and the slides were blocked with 5% goat serum for one hour before overnight incubation with the primary antibody CDKN3 (Bioss, Beijing, bs-5743R). The following day, the slides were incubated with the secondary antibody and subsequently stained with diaminobenzidine (DAB). Hematoxylin was employed for nuclear counterstaining, and images were captured using microscopy.

Prognostic analysis and construction of a predictive nomogram

Univariate and multivariate Cox regression analyses were conducted to identify independent prognostic factors. A nomogram based on risk score and all independent prognostic factors was constructed to predict overall survival for NB patients across the entire cohort. ROC curves, decision curve analysis (DCA), and calibration curves were used to assess the accuracy of the nomogram. Stratified analysis was employed to assess the accuracy of prognostic predictions in relation to various clinicopathological features.

Cell culture

SK-N-SH and SK-N-BE(2) cells were cultured in DMEM (Gibco) and DMEM/F12 (Gibco) supplemented with 10% FBS (Procell) and 1% penicillin-streptomycin (Solarbio). All cell lines were incubated at 37 ℃ with 5% CO₂. NB cells used in this study were harvested during logarithmic growth phase at 70–80% confluence.

Reverse transcription quantitative polymerase chain reaction (RT-qPCR)

Total RNA was isolated from cells or tissues using Trizol reagent (Invitrogen). Following isolation, RNA was reverse-transcribed into cDNA using StarScript III RT Master Mix (GenStar, A233). RT-qPCR amplification was performed on an ABI QuantStudio 5 (Q5) instrument using Real Star Power SYBR qPCR Mix (GenStar, A311) according to the manufacturer’s protocol. Primer sequences as follows: CDKN3: Forward: 5'-TCCGGGGCAATACAGACCAT-3'; Reverse: 5'-GCAGCTAATTTGTCCCGAAACTC-3'. GAPDH: Forward: 5'-GAAGGTGAAGGTCGGAGTC-3'; Reverse: 5'-GAAGATGGTGATGGGATTTC-3'. The GAPDH served as the endogenous control. Relative gene expression was calculated using the 2^−ΔΔCt method with triplicate technical replicates per sample.

Cell proliferation assay

NB cells were transfected with CDKN3 siRNA#2 and control siRNA. Then, cells were plated at in 96-well plates (SK-N-SH: 3,000 cells/well and SK-N-BE2: 8,000 cells/well) in growth media and incubated for 48 h. Using the Cellaview System (AF-100), cell proliferation was continuously measured through automated phase imaging (12-hour intervals) and quantified as confluence percentage by the integrated software.

Statistical analysis

Statistical analyses and data visualization were conducted using R software (version 4.2.1; The R Foundation for Statistical Computing, Vienna, Austria). Two-tailed unpaired Student’s t-tests and Wilcoxon tests were used to assess differences between groups. Pearson’s Chi-squared test or Fisher’s exact tests were utilized to evaluate differences in proportions. P values of less than 0.05 were deemed statistically significant.

Results

Selection of ER stress-related DEGs and construction of a risk model

A total of 19 ER stress-related DEGs were identified based on MYCN amplification status (Figure 1A). Among these, 12 genes were downregulated and 7 were upregulated. Clustering analysis revealed that 7 genes clustered in the MYCN-amplified group, while the remaining 12 genes clustered in the non-MYCN-amplified group (Figure 1B). Univariate Cox regression analysis showed that these 19 genes were closely related to the prognosis of NB patients, with the 7 genes having hazard ratios (HR) greater than 1, indicating high-risk for NB patients (Figure 1C). Subsequent LASSO regression analysis (Figure 1D,1E) identified 5 hub ER stress-related DEGs for the risk score model (PINK1, IL7, CDKN3, C1S, MMP9). Notably, only CDKN3 expression was positively associated with MYCN amplification and acted a risk factor for NB progression.

Figure 1 Identification and clustering of ER stress-related DEGs in neuroblastoma. (A) Volcano plot displaying the distribution of ER stress-related DEGs between MYCN amplified and non-amplified NB patients. Red dots represent upregulated genes, and green dots represent downregulated genes. (B) Heatmap showing the clustering of ER stress-related DEGs in MYCN amplified and non-amplified NB patients. Each column represents a patient sample, and each row represents a gene. (C) Forest plot from univariate Cox regression analysis indicating the HR of the 19 identified ER stress-related DEGs. (D) LASSO regression analysis to select the hub ER stress-related DEGs. (E) Cross-validation plot for tuning parameter selection in the LASSO model. CI, confidence interval; DEGs, differentially expressed genes; ER, endoplasmic reticulum; FC, fold change; HR, hazard ratio; LASSO, least absolute shrinkage and selection operator; NB, neuroblastoma.

Development and Validation of an ER stress-related five-gene risk signature

The risk score calculation formula is as follows: risk score = (−3.3605 × PINK1 expression) + (−0.8256 × IL7 expression) + (6.4808 × CDKN3 expression) + (2.381 × C1S expression) + (−1.8428 × MMP9 expression). All NB cases in the training set were categorized into high- and low-risk groups based on the median risk score (Figure 2A). The heatmap highlights the differential expression of these five genes across different risk categories, while the scatter plot demonstrates that an increase in risk score correlates with higher mortality rates and decreased overall survival times (Figure 2A). Furthermore, KM survival curves (Figure 2B) show that the gene model effectively predicts prognosis, with patients in the high-risk group experiencing significantly worse survival outcomes compared to those in the low-risk group (P<0.001). Additionally, ROC curves (Figure 2C) indicate area under the curve (AUC) values for 1-, 3-, and 5-year overall survival (OS) predictions of 0.947, 0.829, and 0.812, respectively, demonstrating the robust predictive performance of the risk model. The testing set was used to evaluate the prognostic efficacy of the risk score. Similar observations to the training set were made, with the high-risk group exhibiting shorter OS and notable survival disparities between the groups (Figure S1A), thereby confirming the accuracy of the five-gene signature in classifying patient risk. Finally, we successfully differentiated between high- and low-risk groups in the entire set using the gene risk model (Figure S1B).

Figure 2 Development and validation of the er stress-related five-gene risk signature. (A) Heatmap of the expression levels of the five genes (PINK1, IL7, CDKN3, C1S, MMP9) in high- and low-risk groups, along with the corresponding risk score distribution and survival status. (B) KM survival curves comparing overall survival between high- and low-risk groups in the training set. (C) ROC curves showing the predictive accuracy of the five-gene risk signature for 1-, 3-, and 5-year OS in the training set. KM, Kaplan-Meier; OS, overall survival; ROC, receiver operating characteristic.

Correlation between risk scores and clinicopathological features

The correlation analysis between the risk model and clinical characteristics highlighted significant variations in risk scores across different age groups, International Neuroblastoma Staging System (INSS) stages, histology, MYCN amplification status, tumor progression, NB risk group (Figure 3A-3G). High risk scores were positively correlated with several common clinical and pathological factors that indicate poor prognosis, such as older age, advanced clinical stages, unfavorable histology and MYCN amplification. However, no notable differences in risk scores were observed between genders (Figure 3H).

Figure 3 Correlation between risk scores and clinicopathological features. (A) Boxplot showing the distribution of risk scores across different age groups. (B) Boxplot displaying risk scores according to INSS stages. (C) Boxplot indicating the relationship between risk scores and histological classifications. (D) Risk score distribution in relation to MYCN amplification status. (E) Boxplot correlating risk scores with tumor progression. (F) Distribution of risk scores across different neuroblastoma risk groups. (G) Boxplot showing the correlation between risk scores and disease-specific mortality. (H) Boxplot comparing risk scores between genders. INSS, International Neuroblastoma Staging System; FH; UFH.

Prognostic analysis of risk model and clinicopathological features

NB patients were divided into high- and low-risk categories based on their median risk score. Subgroup survival analysis revealed that, irrespective of age (Figure 4A,4B) or INSS stages (Figure 4C,4D), patients in the high-risk category consistently exhibited a poorer prognosis than those in the low-risk category. Furthermore, within the subset of patients clinically diagnosed as high-risk NB (Figure 4E) and those experiencing tumor progression (Figure 4F), cases with a high-risk score still demonstrated a significant reduction in OS.

Figure 4 Subgroup survival analysis based on the five-gene risk signature. (A) KM survival curves for patients age <18 months comparing high- and low-risk groups. (B) KM survival curves for patients age ≥18 months comparing high- and low-risk groups. (C) KM survival curves for patients with INSS stage 1, 2, 3 and 4S comparing high- and low-risk groups. (D) KM survival curves for patients with INSS stage 4 comparing high- and low-risk groups. (E) KM Survival curves for clinically considered high-risk NB patients comparing high- and low-risk groups. (F) KM survival curves for patients experiencing tumor progression comparing high- and low-risk groups. INSS, International Neuroblastoma Staging System; KM, Kaplan-Meier; NB, neuroblastoma.

GSEA analysis

The analysis indicated that the high-risk group showed enrichment in pathways related to chromosome segregation, nuclear chromosome segregation, regulation of chromosome segregation, regulation of mitotic nuclear division, sister chromatid segregation, cell cycle, DNA replication, homologous recombination, ribosome, and spliceosome (Figure S2A,S2B). Conversely, the low-risk group exhibited enrichment in pathways associated with cell junction assembly, neurotransmitter transport, regulation of trans-synaptic signaling, vesicle-mediated transport in synapses, presynaptic processes, axon guidance, calcium signaling, cell adhesion molecules (CAMs), hematopoietic cell lineage, and lysosome (Figure S2C,S2D).

Tumor immunity relevance of the ER stress‑related signature

To investigate the correlation between ER stress-related risk scores and the TME in NB, we deployed a series of analytical algorithms. We found that T cell infiltration had a negative correlation with the risk score. In the high-risk group, the infiltration of CD8⁺ and CD4⁺ T cells was significantly lower (Figure 5A). Additionally, the high-risk group presented a distinct immune profile characterized by a diminished ESTIMATE Score, indicating a lower level of stromal and immune cell infiltration in their TME (Figure 5B). Moreover, NB patients with high levels of immune cell infiltration, including CD8⁺ Tem cells, CD4⁺ Tem cells, activated dendritic cells (aDCs), natural killer T (NKT) cells, CD4 memory resting T cells, myeloid dendritic cells, and natural killer (NK) cells, had a better prognosis (Figure S3).

Figure 5 Immune infiltration and drug sensitivity analysis based on ER risk groups. (A) Scatterplot comparing the infiltration levels of immunocyte between high- and low-risk groups. (B) Boxplot showing the ESTIMATE algorithm results for stromal and immune cell infiltration in the tumor microenvironment of high- and low-risk groups. (C) Boxplot indicating the expression levels of crucial immune checkpoints in high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001. (D) Boxplot comparing the predicted IC₅₀ values for various anticancer drugs between high-risk and low-risk groups. ER, endoplasmic reticulum; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; IC₅₀, half maximal inhibitory concentration.

To further explore the immune status of the two groups, we estimated the enrichment value of different immune cell infiltrations and immune functions. 16 types of immunocytes were predominantly enriched in the low-risk group, whereas mast cells were notably upregulated in the high-risk group (Figure S4A). Additionally, the ssGSEA scores for immune functions revealed that most of immune functions exhibited higher enrichment scores in the low-risk group compared to the high-risk group (Figure S4B). An analysis of crucial immune checkpoint expressions uncovered significant upregulation of immune suppressive molecules such as CD276 (B7-H3), CD70, CD86, lymphocyte activation gene 3 (LAG3), tumor necrosis factor receptor superfamily member 8 (TNFRSF8), and V-set domain containing T cell activation inhibitor 1 (VTCN1) in the high-risk group (Figure 5C). This suggests that patients in the high-risk group may exhibit a compromised anti-tumor immune response.

Prediction of sensitive drugs

Utilizing the GDSC database, we predicted the drug responsiveness of NB patients by comparing half maximal inhibitory concentration (IC₅₀) values between high- and low-risk groups. Notably, the IC₅₀ values for ibrutinib, gefitinib, crizotinib, cisplatin, lapatinib, and irinotecan were significantly lower in the high-risk group, indicating enhanced sensitivity to these drugs among these patients (Figure 5D).

Independent prognostic analysis and construction of a nomogram

Univariate and multifactorial analyses demonstrated that a high ER stress-related risk score is associated with poor prognosis in NB patients and serves as an independent prognostic indicator (Figure 6A,6B). To enhance prognostic accuracy, we integrated the ER stress risk model with other independent prognostic factors, including MYCN status, age, and INSS stage, to develop a comprehensive nomogram. This nomogram assigns specific values to each case, where lower values indicate more favorable patient outcomes (Figure 6C). The reliability of the nomogram was confirmed through ROC curves and calibration curves (Figure 6D,6E). The nomogram exhibited superior predictive performance compared to using the ER stress risk model or clinical indicators alone. The AUC values for 1-, 3-, and 5-year OS predictions were 0.905, 0.831, and 0.805, respectively (Figure S4C). DCA further highlighted the predictive nomogram’s advantage, demonstrating the greatest clinical benefit (Figure 6F).

Figure 6 Independent prognostic analysis and nomogram construction. (A) Forest plot from univariate Cox regression analysis showing the prognostic significance of the five-gene risk score and other clinical variables. (B) Multivariate Cox regression analysis indicating the independent prognostic value of the five-gene risk score. (C) Nomogram integrating the five-gene risk score and clinical variables to predict overall survival for neuroblastoma patients. (D) ROC curves evaluating the predictive accuracy of the nomogram for overall survival. (E) Calibration curves demonstrating the agreement between predicted and observed survival outcomes using the nomogram. (F) DCA illustrating the clinical benefit of the nomogram. DCA, decision curve analysis; ROC, receiver operating characteristic.

Elevated CDKN3 expression in NB correlates with poor prognosis

To further validate our findings, we assessed the expression of CDKN3 in NB tissues from 46 patients using IHC (Figure 7A). In this cohort, patients with high CDKN3 expression demonstrated significantly shorter progression-free survival (PFS) and OS (Figure 7B,7C), affirming that high CDKN3 expression is an important indicator of poor prognosis in NB patients, consistent with its identified risk level in our model. As anticipated, CDKN3 expression positively correlated with MYCN amplification (Figure 7D). Additionally, its expression was positively associated with higher INSS stages, metastasis at diagnosis, and elevated levels of neuron-specific enolase (NSE)—all markers indicative of poor prognosis (Figure 7E-7G).

Figure 7 Immunohistochemical validation of CDKN3 expression in neuroblastoma tissues. (A) Representative IHC images showing CDKN3 expression in neuroblastoma tissues. (B) KM survival curves comparing DFS between patients with high and low CDKN3 expression. (C) KM survival curves comparing OS between patients with high and low CDKN3 expression. (D) Boxplot showing the correlation between CDKN3 expression and MYCN amplification status. (E) Boxplot indicating the relationship between CDKN3 expression and INSS stages. (F) Boxplot displaying the association between CDKN3 expression and metastasis at diagnosis. (G) Boxplot showing the correlation between CDKN3 expression and NSE levels. IHC, immunohistochemistry; KM, Kaplan-Meier; DFS, disease-free survival; OS, overall survival; INSS, International Neuroblastoma Staging System; NSE, Neuron-Specific Enolase.

Downregulation of CDNK3 expression inhibits the proliferation of NB cells

To further investigate the function of CDKN3 expression on NB proliferation, we transfected NB cells with CDKN3 siRNAs. We found that only CDKN3 siRNA#2 effectively suppressed CDKN3 expression (Figure 8A,8B). Consequently, CDKN3 siRNA#2 was selected for further experiments. SK-N-SH and SK-N-BE2 cells were transfected with CDKN3 siRNA#2 or control siRNA. then cells were seeded into 96-well plates and incubated in a CellView system for live-cell imaging over 48 hours. The result showed that the cell density in the control was significantly higher than that in the CDKN3-knockdown group (Figure 8C,8D). Furthermore, compared to the control group (black curve), the confluence rate curve for the siCDKN3#2 group (red curve) was markedly reduced (Figure 8E,8F). Collectively, these results suggested that down-regulation of CDKN3 significantly inhibits the proliferation of NB cells.

Figure 8 Down-regulation of CDKN3 inhibit the cell proliferation of NB cells. (A,B) RT-qPCR analysis of CDKN3 knockdown efficiency in SK-N-SH and SK-N-BE2 cells. (C,D) Images of siCDKN3 and control group cells at 0 and 48 h. (E,F) Cell confluence curves of NB cells transfected with siCDKN3 and control from 0 to 48 h. ns, no significant; *, P<0.05; ***, P<0.001; ****, P<0.0001. NB, neuroblastoma; RT-qPCR, reverse transcription quantitative polymerase chain reaction.

Discussion

Treating high-risk NB remains a formidable challenge, emphasizing the need for precise stratification and accurate prognosis assessment of NB patients. There is growing evidence that ER stress is linked to the malignant behavior of NB and treatment resistance (31,32). For instance, París-Coderch et al. demonstrated that ABTL0812 induces NB cell death through the induction of ER stress, activation of the UPR, autophagy, and apoptosis (33). Consequently, in this study, we aimed to explore the relationship between ER stress and NB prognosis.

Considering that MYCN amplification is the most well-defined genetic marker of high risk in NB and is closely associated with ER stress, we used MYCN amplification status to identify ER stress-related DEGs. We then developed a risk model based on five key ER stress-related DEGs—PINK1, IL7, CDKN3, C1S, and MMP9—to evaluate survival outcomes in NB. This ER stress-related signature effectively categorizes the risk levels in NB patients. Our survival analysis indicated that patients in the high-risk group had worse prognoses regardless of their age and INSS stage, and the signature was identified as an independent prognostic factor for NB. Furthermore, we constructed a nomogram that integrates the risk score with several clinical variables to enhance the predictive accuracy of patient outcomes. The effectiveness of our nomogram was substantiated through calibration plots, confirming its utility in clinical decision-making for NB patients.

The five genes identified in our study were linked to the onset and progression of various malignancies. PTEN-induced kinase 1 (PINK1) is a serine/threonine kinase closely associated with mitochondrial function and cellular autophagy (34). Yao et al. suggested that PINK1-mediated autophagy, also known as mitophagy, plays a crucial role in maintaining mitochondrial integrity by eliminating damaged mitochondria, thereby supporting cellular health and function (35). PINK1 has been reported to inhibits glioblastoma growth (36), while in lung cancer, increased PINK1 expression promotes proliferation and chemoresistance (37). In NB, PINK1 has been identified as potentially harboring pathogenic germline variants that are significantly enriched among NB cases, indicating a role in genetic susceptibility to the disease (38). Additionally, inhibiting the degradation of PINK1 can enhance its neuroprotective effects in human NB cells (39). Interleukin 7 (IL-7) is a cytokine essential for the adaptive immune system. IL-7 can directly or indirectly influence tumor cells and exert anti-tumor effects by enhancing tumor eradication or boosting adaptive immunity (40). It inhibits melanoma growth and affects the invasion and growth of bladder, prostate, and lung cancer cells (41-43). In NB, IL-7 enhances the antitumor activity of chimeric antigen receptor T-cell (CAR T) cells targeting ganglioside D2 (GD2)-positive cells (44). CDKN3 has been proven to regulate cell survival and proliferation across various types of cancer (45). Recently, research revealed that CDKN3 could induce differentiation in NB cells and that knocking down CDKN3 could inhibit their proliferation (44). In this study, Vernaza et al. discovered that MYCN directly regulates CDKN3 transcription in NB cells, with CDKN3 expression levels showing a positive correlation with MYCN expression in NB specimens (46). These findings align with our own, wherein CDKN3 was positively correlated with MYCN and identified as an indicator of poor prognosis. Therefore, CDKN3 not only serves as a valuable prognostic indicator for NB but also represents a potential therapeutic target. C1S, a component of the C1 complex initiator of the classical complement pathway, is linked to various cancer development (47-49). A previous study indicates that C1S is implicated in the secretion of complement components by human NB cell lines, suggesting a role in immune response modulation within the TME of NB (50). However, further studies on the relationship between C1S and NB are lacking. Matrix metalloproteinase 9 (MMP9), a key member of the MMP family, is pivotal in degrading the extracellular matrix and enhancing tumor malignancy (51). It is upregulated in advanced stages of NB, where it is associated with increased angiogenesis and tumor progression (52,53). By facilitating the breakdown of the extracellular matrix and supporting angiogenesis, MMP-9 critically contributes to the invasive and metastatic characteristics of NB, thereby potentially aiding the dissemination of tumor cells. In conclusion, our research has identified five genes closely associated with MYCN and ER stress. Several of these genes have been implicated in the development of NB, while others have been examined in different malignant tumors. Further studies are necessary to clarify the functions of these ER stress-related genes in NB.

The advent of immunotherapy has introduced new hope for patients with NB. However, the effectiveness of treatments such as immune checkpoint inhibitors is often limited due to NB’s characteristic low immune cell infiltration (54). While anti-GD2 monoclonal antibodies have shown efficacy in high-risk NB cases (55), their effectiveness can vary, with poor outcomes in some patients due to insufficient CD8 T-cell infiltration (56). Our study indicates that high-risk NB patients, characterized by low CD4⁺ and CD8⁺ T-cell infiltration and poor prognosis, may benefit less from immunotherapy. We also discovered that elevated expression of immune checkpoints—CD276 (B7-H3), CD70, CD86, LAG3, TNFRSF8, and VTCN1—in the ER high-risk group contributes to tumor immune evasion, suggesting a more complex and immunosuppressive TME in these patients. Notably, CD276 (B7-H3), which is overexpressed in various human cancers including NB and often associated with poor outcomes, protects NB cells from NK cell cytotoxic activity (57). Due to its role in tumor immune evasion, B7-H3 has become a novel target for immunotherapy, with B7-H3-targeted CAR-T cells recently approved for NB treatment. Ultimately, our study has demonstrated that the five-biomarker risk classifier can effectively differentiate levels of tumor immune infiltration and predict NB response to immunotherapy. Further exploration of the mechanisms linking ER stress with immune regulation in NB cells may enhance our understanding and treatment of this immunologically ‘cold’ tumor.

Current targeted therapies are not universally effective for all NB patients, prompting our assessment of drug efficacy based on risk-group stratification. Our research discovered that an ER stress-related signature strongly correlates with sensitivity to specific targeted drugs—crizotinib, ibrutinib, gefitinib, and lapatinib, which respectively target anaplastic lymphoma kinase (ALK), Bruton’s tyrosine kinase (BTK), epidermal growth factor receptor (EGFR), and EGFR/human epidermal growth factor receptor 2 (HER2). Notably, these drugs are more effective in ER high-risk patients. Given ALK’s role in enhancing cell survival, migration, and proliferation in NB, substantial efforts are being made to develop ALK inhibitors like crizotinib, which is currently undergoing phase III clinical trials for NB (4). Meanwhile, research has shown that BTK amplifies ALK signaling, and a combined treatment of ibrutinib with crizotinib significantly enhances their inhibitory effects (58). EGFR is expressed in NB cell lines and tumors, and its inhibition by gefitinib induces apoptosis, but gefitinib alone has limited efficacy against NB xenografts in vivo (59,60). Tee et al. highlighted that a combination of the EGFR/HER2 inhibitor lapatinib with YM155 forms an effective therapeutic approach for NB, regardless of ATP binding cassette subfamily B member 1 (ABCB1) status (61). Our results indicate that the ER stress-related signature may act as a predictive biomarker for NB patients treated with these targeted drugs, potentially aiding in the selection of treatment strategies, though further research is required to confirm these implications.

Our study faces certain limitations. Our analysis utilized retrospective data from a public database, which lacked comprehensive details for all NB cases. The NB tissue samples were also retrospectively collected. Therefore, it is crucial to validate these findings with expanded research that includes larger sample sizes and comprehensive prospective studies. Furthermore, the regulation of ER by these five ER-related genes in NB has not been confirmed, necessitating additional in vivo and in vitro experiments for verification.

Conclusions

We have developed and validated an ER stress-related gene signature that accurately predicts the prognosis of NB patients. A nomogram based on this signature, combined with clinical indicators, effectively forecasts patient outcomes. Further investigations indicate that this signature may also predict responses to immunotherapy and sensitivity to anti-tumor drugs, thus enhancing clinical decision-making and treatment planning for NB patients. Additionally, we have conducted preliminary validation of CDKN3 in NB clinical samples, suggesting its potential as a new target in ER stress-related research for NB. Our study highlights the importance of ER stress-related genes as prognostic markers and affirms the critical role of ER stress in NB progression.

Acknowledgments

We extend our gratitude to the GEO and Array Express for providing essential gene expression data.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-142/rc

Data Sharing Statement: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-142/dss

Peer Review File: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-142/prf

Funding: This research was supported by Scientific Research Project of Tianjin Municipal Education Commission (No. 2024KJ175), the Tianjin Binhai New Area Health Science and Technology Project (grant No. 2023BWKQ016), the Introduction of Talented Personnel and Doctoral Startup Funding (No. B2304), and Tianjin Key Medical Discipline (Specialty) Construction Project (Nos. TJYXZDXK-009A and TJYXZDXK-012A).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-142/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study received ethics approval from the Ethics Committee of Tianjin Medical University Cancer Institute and Hospital (approval No. E20210027) and was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Written informed consent was obtained from the legal guardians or next of kin of all participants.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Park JR, Eggert A, Caron H. Neuroblastoma: biology, prognosis, and treatment. Hematol Oncol Clin North Am 2010;24:65-86. [Crossref] [PubMed]
2. Chung C, Boterberg T, Lucas J, et al. Neuroblastoma. Pediatr Blood Cancer 2021;68:e28473. [Crossref] [PubMed]
3. Guan Q, Lin H, Hua W, et al. Variant rs8400 enhances ALKBH5 expression through disrupting miR-186 binding and promotes neuroblastoma progression. Chin J Cancer Res 2023;35:140-62. [Crossref] [PubMed]
4. Lin L, Wang B, Zhang X, et al. Functional TET2 gene polymorphisms increase the risk of neuroblastoma in Chinese children. IUBMB Life 2024;76:200-11. [Crossref] [PubMed]
5. Maris JM, Hogarty MD, Bagatell R, et al. Neuroblastoma. Lancet 2007;369:2106-20. [Crossref] [PubMed]
6. Zafar A, Wang W, Liu G, et al. Molecular targeting therapies for neuroblastoma: Progress and challenges. Med Res Rev 2021;41:961-1021. [Crossref] [PubMed]
7. Huang M, Weiss WA. Neuroblastoma and MYCN. Cold Spring Harb Perspect Med 2013;3:a014415. [Crossref] [PubMed]
8. Campbell K, Gastier-Foster JM, Mann M, et al. Association of MYCN copy number with clinical features, tumor biology, and outcomes in neuroblastoma: A report from the Children's Oncology Group. Cancer 2017;123:4224-35. [Crossref] [PubMed]
9. Ruiz-Pérez MV, Henley AB, Arsenian-Henriksson M. The MYCN Protein in Health and Disease. Genes (Basel) 2017;8:113. [Crossref] [PubMed]
10. Zhang T, Li N, Sun C, et al. MYC and the unfolded protein response in cancer: synthetic lethal partners in crime? EMBO Mol Med 2020;12:e11845. [Crossref] [PubMed]
11. Schmidt S, Gay D, Uthe FW, et al. A MYC-GCN2-eIF2α negative feedback loop limits protein synthesis to prevent MYC-dependent apoptosis in colorectal cancer. Nat Cell Biol 2019;21:1413-24. [Crossref] [PubMed]
12. Yaari-Stark S, Shaked M, Nevo-Caspi Y, et al. Ras inhibits endoplasmic reticulum stress in human cancer cells with amplified Myc. Int J Cancer 2010;126:2268-81. [Crossref] [PubMed]
13. Wu J, Qiao S, Xiang Y, et al. Endoplasmic reticulum stress: Multiple regulatory roles in hepatocellular carcinoma. Biomed Pharmacother 2021;142:112005. [Crossref] [PubMed]
14. Wang M, Kaufman RJ. Protein misfolding in the endoplasmic reticulum as a conduit to human disease. Nature 2016;529:326-35. [Crossref] [PubMed]
15. Urra H, Dufey E, Avril T, et al. Endoplasmic Reticulum Stress and the Hallmarks of Cancer. Trends Cancer 2016;2:252-62. [Crossref] [PubMed]
16. Di Prisco GV, Huang W, Buffington SA, et al. Translational control of mGluR-dependent long-term depression and object-place learning by eIF2α. Nat Neurosci 2014;17:1073-82. [Crossref] [PubMed]
17. Hollien J, Weissman JS. Decay of endoplasmic reticulum-localized mRNAs during the unfolded protein response. Science 2006;313:104-7. [Crossref] [PubMed]
18. Lin JH, Li H, Yasumura D, et al. IRE1 signaling affects cell fate during the unfolded protein response. Science 2007;318:944-9. [Crossref] [PubMed]
19. Rutkowski DT, Arnold SM, Miller CN, et al. Adaptation to ER stress is mediated by differential stabilities of pro-survival and pro-apoptotic mRNAs and proteins. PLoS Biol 2006;4:e374. [Crossref] [PubMed]
20. Chen X, Cubillos-Ruiz JR. Endoplasmic reticulum stress signals in the tumour and its microenvironment. Nat Rev Cancer 2021;21:71-88. [Crossref] [PubMed]
21. Oakes SA. Endoplasmic Reticulum Stress Signaling in Cancer Cells. Am J Pathol 2020;190:934-46. [Crossref] [PubMed]
22. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015;43:e47. [Crossref] [PubMed]
23. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
24. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
25. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 2017;18:220. [Crossref] [PubMed]
26. Kimes PK, Liu Y, Neil Hayes D, et al. Statistical significance for hierarchical clustering. Biometrics 2017;73:811-21. [Crossref] [PubMed]
27. Metsalu T, Vilo J. ClustVis: a web tool for visualizing clustering of multivariate data using Principal Component Analysis and heatmap. Nucleic Acids Res 2015;43:W566-70. [Crossref] [PubMed]
28. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
29. Yang W, Soares J, Greninger P, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res 2013;41:D955-61. [Crossref] [PubMed]
30. Geeleher P, Cox NJ, Huang RS. Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol 2014;15:R47. [Crossref] [PubMed]
31. Zhang D, Wang F, Pang Y, et al. Down-regulation of CHERP inhibits neuroblastoma cell proliferation and induces apoptosis through ER stress induction. Oncotarget 2017;8:80956-70. [Crossref] [PubMed]
32. Chen YY, Chen G, Fan Z, et al. GSK3beta and endoplasmic reticulum stress mediate rotenone-induced death of SK-N-MC neuroblastoma cells. Biochem Pharmacol 2008;76:128-38. [Crossref] [PubMed]
33. París-Coderch L, Soriano A, Jiménez C, et al. The antitumour drug ABTL0812 impairs neuroblastoma growth through endoplasmic reticulum stress-mediated autophagy and apoptosis. Cell Death Dis 2020;11:773. [Crossref] [PubMed]
34. Li Y, Chen H, Xie X, et al. PINK1-Mediated Mitophagy Promotes Oxidative Phosphorylation and Redox Homeostasis to Induce Drug-Tolerant Persister Cancer Cells. Cancer Res 2023;83:398-413. [Crossref] [PubMed]
35. Yao J, Wang J, Xu Y, et al. CDK9 inhibition blocks the initiation of PINK1-PRKN-mediated mitophagy by regulating the SIRT1-FOXO3-BNIP3 axis and enhances the therapeutic effects involving mitochondrial dysfunction in hepatocellular carcinoma. Autophagy 2022;18:1879-97. [Crossref] [PubMed]
36. Agnihotri S, Golbourn B, Huang X, et al. PINK1 Is a Negative Regulator of Growth and the Warburg Effect in Glioblastoma. Cancer Res 2016;76:4708-19. [Crossref] [PubMed]
37. Liu L, Zuo Z, Lu S, et al. Silencing of PINK1 represses cell growth, migration and induces apoptosis of lung cancer cells. Biomed Pharmacother 2018;106:333-41. [Crossref] [PubMed]
38. Pugh TJ, Morozova O, Attiyeh EF, et al. The genetic landscape of high-risk neuroblastoma. Nat Genet 2013;45:279-84. [Crossref] [PubMed]
39. Liu Y, Lear TB, Verma M, et al. Chemical inhibition of FBXO7 reduces inflammation and confers neuroprotection by stabilizing the mitochondrial kinase PINK1. JCI Insight 2020;5:e131834. [Crossref] [PubMed]
40. Wang C, Kong L, Kim S, et al. The Role of IL-7 and IL-7R in Cancer Pathophysiology and Immunotherapy. Int J Mol Sci 2022;23:10412. [Crossref] [PubMed]
41. Park SL, Lee EJ, Kim WJ, et al. p27KIP1 is involved in ERK1/2-mediated MMP-9 expression via the activation of NF-κB binding in the IL-7-induced migration and invasion of 5637 cells. Int J Oncol 2014;44:1349-56. [Crossref] [PubMed]
42. Qu H, Zou Z, Pan Z, et al. IL-7/IL-7 receptor axis stimulates prostate cancer cell invasion and migration via AKT/NF-κB pathway. Int Immunopharmacol 2016;40:203-10. [Crossref] [PubMed]
43. Ming J, Jiang G, Zhang Q, et al. Interleukin-7 up-regulates cyclin D1 via activator protein-1 to promote proliferation of cell in lung cancer. Cancer Immunol Immunother 2012;61:79-88. [Crossref] [PubMed]
44. Li G, Zhang Q, Han Z, et al. IL-7 and CCR2b Co-Expression-Mediated Enhanced CAR-T Survival and Infiltration in Solid Tumors. Front Oncol 2021;11:734593. [Crossref] [PubMed]
45. Campbell GJ, Hands EL, Van de Pette M. The Role of CDKs and CDKIs in Murine Development. Int J Mol Sci 2020;21:5343. [Crossref] [PubMed]
46. Vernaza A, Cardus DF, Smith JL, et al. Identification of CDKN3 as a Key Gene that Regulates Neuroblastoma Cell Differentiation. J Cancer 2024;15:1153-68. [Crossref] [PubMed]
47. Daugan MV, Revel M, Russick J, et al. Complement C1s and C4d as Prognostic Biomarkers in Renal Cancer: Emergence of Noncanonical Functions of C1s. Cancer Immunol Res 2021;9:891-908. [Crossref] [PubMed]
48. Ye J, Yang P, Yang Y, et al. Complement C1s as a diagnostic marker and therapeutic target: Progress and propective. Front Immunol 2022;13:1015128. [Crossref] [PubMed]
49. Yu K, Yang H, Lv QL, et al. Construction of a competitive endogenous RNA network and analysis of potential regulatory axis targets in glioblastoma. Cancer Cell Int 2021;21:102. [Crossref] [PubMed]
50. Veerhuis R, Janssen I, De Groot CJ, et al. Cytokines associated with amyloid plaques in Alzheimer's disease brain stimulate human glial and neuronal cell cultures to secrete early complement proteins, but not C1-inhibitor. Exp Neurol 1999;160:289-99. [Crossref] [PubMed]
51. Lee CJ, Jang TY, Jeon SE, et al. The dysadherin/MMP9 axis modifies the extracellular matrix to accelerate colorectal cancer progression. Nat Commun 2024;15:10422. [Crossref] [PubMed]
52. Somasundaram DB, Aravindan S, Major R, et al. MMP-9 reinforces radiation-induced delayed invasion and metastasis of neuroblastoma cells through second-signaling positive feedback with NFκB via both ERK and IKK activation. Cell Biol Toxicol 2023;39:1053-76. [Crossref] [PubMed]
53. Fietta A, Fusco P, Germano G, et al. Neuroblastoma-derived hypoxic extracellular vesicles promote metastatic dissemination in a zebrafish model. PLoS One 2024;19:e0316103. [Crossref] [PubMed]
54. Wienke J, Dierselhuis MP, Tytgat GAM, et al. The immune landscape of neuroblastoma: Challenges and opportunities for novel therapeutic strategies in pediatric oncology. Eur J Cancer 2021;144:123-50. [Crossref] [PubMed]
55. Cheung IY, Cheung NV, Modak S, et al. Survival Impact of Anti-GD2 Antibody Response in a Phase II Ganglioside Vaccine Trial Among Patients With High-Risk Neuroblastoma With Prior Disease Progression. J Clin Oncol 2021;39:215-26. [Crossref] [PubMed]
56. Sait S, Modak S. Anti-GD2 immunotherapy for neuroblastoma. Expert Rev Anticancer Ther 2017;17:889-904. [Crossref] [PubMed]
57. Bottino C, Vitale C, Dondero A, et al. B7-H3 in Pediatric Tumors: Far beyond Neuroblastoma. Cancers (Basel) 2023;15:3279. [Crossref] [PubMed]
58. Li T, Deng Y, Shi Y, et al. Bruton's tyrosine kinase potentiates ALK signaling and serves as a potential therapeutic target of neuroblastoma. Oncogene 2018;37:6180-94. [Crossref] [PubMed]
59. Bondarev N, Ivanenko K, Khabusheva E, et al. MGL S3 Chimeric Enzyme Drives Apoptotic Death of EGFR-Dependent Cancer Cells through ERK Downregulation. Int J Mol Sci 2022;23:12807. [Crossref] [PubMed]
60. Ma G, Tan C, Shan Y, et al. An insulin growth factor-I/II-neutralizing monoclonal antibody in combination with epidermal growth factor receptor inhibitors potently inhibits tumor cell growth. J Cancer 2022;13:1830-6. [Crossref] [PubMed]
61. Tee AE, Ciampa OC, Wong M, et al. Combination therapy with the CDK7 inhibitor and the tyrosine kinase inhibitor exerts synergistic anticancer effects against MYCN-amplified neuroblastoma. Int J Cancer 2020;147:1928-38. [Crossref] [PubMed]