
==== Front
Discov Oncol
Discov Oncol
Discover Oncology
2730-6011
Springer US New York

39292372
1318
10.1007/s12672-024-01318-0
Analysis
Identification and validation of a novel five-gene signature in high-risk MYCN-not-amplified neuroblastoma
http://orcid.org/0000-0003-2772-0393
Wang Jin-Xia 12
Zhang Hong-Yang 12
Yan Zi-Jun 12
Cao Zi-Yang 12
Shao Jing-Bo shaojb@shchildren.com.cn

3
http://orcid.org/0000-0002-6265-6801
Zou Lin zoulin74@126.com
zoulin@shchildren.com.cn

12
1 grid.16821.3c 0000 0004 0368 8293 Clinical Research Unit, Shanghai Children’s Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai, 200062 China
2 grid.16821.3c 0000 0004 0368 8293 Institute of Pediatric Infection, Immunity, and Critical Care Medicine, Shanghai Children’s Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai, 200062 China
3 grid.16821.3c 0000 0004 0368 8293 Department of Hematology and Cancer, Shanghai Children’s Hospital, Shanghai Jiao Tong University School of Medicine, Shanghai, 200062 China
18 9 2024
18 9 2024
12 2024
15 45622 11 2023
5 9 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. 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-nc-nd/4.0/.
Objective

High-risk neuroblastoma patients often have poor outcomes despite multi-treatment options. The risk stratification of high-risk MYCN-not-amplified (HR-MYCN-NA) patients remains difficult. This study aims to identify a gene set signature that can help further stratify HR-MYCN-NA patients for a potential personalized therapeutic strategy.

Methods

Three microarrays and one single-cell RNA sequence dataset were acquired and analyzed. Firstly, the prognostic-related genes (PRGs) in HR-MYCN-NA tumor cells were identified using TARGET-NB and GSE137804 datasets. Then, the prognostic model was established by LASSO-Cox regression, and verified in external cohort (GSE49710, GSE45547). Moreover, a time-dependent receiver operating characteristic curve (ROC) and area under the ROC (AUC) was used to assess survival prediction. A nomogram was established to predict the 1-, 3- and 5-year overall survival (OS) of HR-MYCN-NA patients.

Results

In the training set, a five-PRGs signature, which include GAL, GFRA3, MARCKS, PSMD13, and ZNHIT3 genes, was identified and successfully stratified HR-MYCN-NA patients into ultra-high risk (UHR) and high-risk (HR) subtypes (HR = 4.29, P < 0.001). ROC curve analysis confirmed its predictive power (AUC = 0.74–0.82), suggesting a good predictive efficacy. Consistently, high-risk scores also predicted worse OS (HR = 2, P = 0.033) in the external validation dataset (AUC = 0.67–0.71). Moreover, the overall C-index of the nomogram was 0.75 (P < 0.001), which indicated good agreement between the observed and predicted survival rates. Further integrating the five PRGs signature with clinical factors, these 5 gene signature (HR = 4.45, P < 0.001) and tumor grade (HR = 4.15, P = 0.02) were found to be independent prognostic factors for HR-MYCN-NA patients.

Conclusion

The novel five PRGs signature could well predict the survival of HR-MYCN-NA patients, which may provide constructive information for these subsets.

Supplementary Information

The online version contains supplementary material available at 10.1007/s12672-024-01318-0.

Keywords

High-risk neuroblastoma
MYCN-not-amplified
Prognostic-related genes
Signature
Shanghai Sailing Program22YF1437200 Wang Jin-Xia http://dx.doi.org/10.13039/501100001809 National Natural Science Foundation of China 82270160 Zou Lin issue-copyright-statement© Springer Science+Business Media, LLC 2024
==== Body
pmcIntroduction

Neuroblastoma (NB) is one of the most common extracranial solid tumors in children, which accounts for 8–10% of pediatric malignant tumors and 15% of tumor-related deaths [1, 2]. Due to the biological and clinical heterogeneity, different treatment strategies are adopted according to risk stratification [3]. The 5-year overall survival (OS) of NB patients with low- or intermediate-risk is about 80–90%, while for high-risk neuroblastoma (HR-NB), the 5-year OS is less than 50% despite intensive multimodal therapy [4]. Therefore, risk stratification is essential to determine the appropriate treatment for HR-NB children.

Currently, there are two commonly used risk stratification systems. The International Neuroblastoma Risk Group (INRG) classification system, with 7 potential prognostic factors, including age, tumor stage, MYCN status, histology, DNA ploidy, and the 11q aberration, stratified patients into four risk groups: very low, low, intermediate, and high risk [3]. By the INRG criteria, patients with 5-year EFS less than 50% were classified as high risk and will treated with intensive, multimodality strategies. While the Children’s Oncology Group (COG) used age, INSS stage, histology, MYCN status, and DNA ploidy to stratify patient [5]. Recently, in 2021, COG reported the revised risk classification system (version 2), which incorporated segmental chromosome aberrations (1p and 11q) as an additional genomic biomarker [6]. These criteria are very important for treatment, and a proper stratification of the high-risk NB patients will help physicians to develop more customized treatment strategy for different subgroups, which may reduce overtreatment or under treatment of those individuals.

MYCN amplification (MYCN-A) is one of the strongest predictors of poor prognosis of high-risk neuroblastoma, such children are prone to invasive progression, and the OS is less than 40% [7]. However, MYCN gene amplification is detected in approximately 30–50% of HR-NB cases [7–9], indicating that about 50–70% of HR-NB are MYCN-not-amplified (MYCN-NA) tumors. Moreover, NB patients without MYCN-A are not always favorable [10, 11]. Shahab reported that the 5-year progression-free survival (PFS) of children with MYCN-NA tumors with a favorable prognosis was 79%, while the PFS of children with an unfavorable prognosis was only 16% [12]. ALK amplification is also a prognostic feature of NB. Studies showed that 1–2% of HR-NB are associated with ALK amplification, and ALK amplification tends to occur simultaneously with MYCN amplification [13, 14]. Moreover, TERT rearrangements or overexpression was found to define a subgroup of approximately 20% HR-NB, these rearrangements occurred only in high-risk MYCN-NA patients [15–17]. These indicates the genomic profile may play a critical role in HR-NB risk stratification.

Recent years, with high-throughput technologies have been used to analyze molecular characters of NB, it is generally believed that gene signatures play a critical role in the prognosis of NB. A proper stratification of HR-NB patients is important to make individualized treatment plans for different subgroups. However, there is still a lack of survival stratification for the high-risk MYCN-NA neuroblastoma. Here, we established a prognostic related gene signature for HR-MYCN-NA tumors both from the TARGET-NB and GSE137804 datasets, and validated in the external GSE49710, GSE45547 cohort. We further constructed a nomogram that combined with 5 gene signatures to predict the 1-, 3- and 5-year OS, and used Cox regression to identify the independent prognostic factors and clinical relationship of the signature in HR-MYCN-NA patients.

Methods

Data collection

Four public datasets, include TARGET-NB, GSE137804, GSE49710, and GSE45547, were used. The TARGET-NB was downloaded from the Therapeutically Applicable Research to Generate Effective Treatments (https://ocg.cancer.gov/programs/ target). The clinical data included age, gender, race, MYCN, stage, ploidy, risk status, and OS, and 143 HR-MYCN-NA tumors were enrolled in our analysis. The GSE49710, GSE45547 were downloaded from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). The clinical data, including age, gender, MYCN status, stage, risk status, OS, and 83 HR- MYCN-NA patients were applied in our study in GSE49710. In addition, a single-cell RNA (scRNA) sequencing data (GSE137804) [18] was downloaded from GEO too. This database involved 16 NB patients and 160,910 cells, which include 8 HR-MYCN-NA tumors, as applied in our research.

Identify prognostic-related genes (PRGs) of HR-MYCN-NA patients

The TARGET-NB and scRNA dataset (GSE137804) were used to identify PRGS in neuroblastoma tumor cells. Firstly, for TARGET-NB dataset, univariate Cox regression analysis was applied to screen out marker genes related to OS in HR-MYCN-NA patients, and genes with P < 0.05 were selected for subsequent analysis. Then, we used GSE137804 dataset to select malignant tumor cells putatively by the author and identified marker genes The “Seurat” package was employed to analyze scRNA-seq data. The “FindAllMarkers” function was used to identify the marker genes of malignant tumor cells. Adjusted P < 0.05 and log2(fold change) > 0.25 were used as the cutoff threshold values to identify marker genes. Furthermore, we overlapped the two above gene lists to determine the prognostic marker genes in HR-MYCN-NA malignant tumor cells (Details of the genes lists please see supplementary S1-S2).

Construction of PRGs signatures for HR-MYCN-NA patients

We used TARGET-NB database as the training set to construct the prognosis model. To minimize the prognostic gene signature overfitting, we applied the least absolute shrinkage and selection operator (Lasso) combined with Cox proportional hazards regression model, and the penalty parameter was estimated by tenfold cross-validation [19]. Then genes related prognostic model was constructed by multivariate Cox regression. The prognostic related gene signatures were calculated according to the following formula: risk score = Ʃ (Coe × Exp), where Coe, the coefficients, represented the weight of the respective signature and Exp represented the expression value. Then, patients were divided into ultra-high risk and high-risk group with the median risk score as the cut-off value. The Kaplan–Meier (K–M) was used to evaluate the differences in survival rate between the two groups. Heat map was generated in Tree View, with the normalized z-score displayed in each row. Furthermore, a time-dependent ROC curve (tROC), and the area under the ROC curve (AUC) were used to assess the sensitivity and specificity of the PRGs signature [20].

Validation of five gene signatures

To verify the predictive value of PRGs signature, the GSE49710 was used as the external verification cohort. The risk scores were calculated by the same formula. Survival and tROC curve were performed as described above. Moreover, the GSE45547 database was used to validate the expression of five genes in NB patient tissue, and Gene Ontology (GO) was used to analysis genes functional enrichment by the “clusterProfiler” R package.

Association of the PRGs signatures and clinical features

A nomogram to predict the 1-year, 3-year and 5-year OS was constructed according to the results of multivariate Cox regression analysis. We used the calibration plot to evaluate the prognostic accuracy of the nomogram. In addition, univariate and multivariate Cox analysis were used to assess the impact of risk score on OS and the clinical characteristics (age, gender, stage and grade). tROC curve was performed to compare the accuracy of the prediction between the clinical features and risk score for HR-MYCN-NA patients.

Statistical analysis

Statistical analysis was performed by R 3.6.1. The Kaplan–Meier survival curves were performed by the survminer package. The tROC curves and AUC analysis were performed by the ‘time ROC’. The glmnet R package was used for LASSO regression. The gene expression, survival status and survival time were evaluated by ggrisk package. The predictive value of risk score and clinical characteristics was assessed by Cox regression analysis. P-values and hazard ratio (HR) with 95% confidence interval (CI) were generated by log-rank tests. A P < 0.05 was considered statistically significant.

Results

Identification of PRGs in HR-MYCN-NA tumor cells

The flow chart of this study is shown in the Fig. 1. There are 143 clinical samples of the HR-MYCN-NA tumors in the TARGET-NB dataset. The result of univariate Cox regression model showed that 1743 candidate genes were significantly associated with the OS of HR-MYCN-NA tumors, which include 646 risk gene factors and 1097 protective gene factors (supplementary S1). For the scRNA-seq dataset (GSE137804), there are 8 clinical samples of HR-MYCN-NA tumors. A total of 3254 marker genes expressed in tumor cells was obtained (supplementary S2). Then, we overlapped these two marker gene list to obtain a total of 435 genes. To further reduce the number of candidate genes, we selected the top 200 genes, and used LASSO regression analysis to identify 16 genes. Furthermore, multivariate Cox regression was conducted and 5 genes was applied to establish a risk model. Kaplan–Meier analysis indicated that higher expression of GAL, PSMD13 and lower expression of MARCKS, GFRA3, ZNHIT3 are significantly associated with poorer OS (P < 0.05, Fig. 2).Fig. 1 The overall workflow of this study

Fig. 2 Association of overall survival and gene expression of GAL, PSMD13, MARCKS, GFRA3, and ZNHIT3 in the TARGET-NB cohort

Construction and evaluation PRGs signatures for HR-MYCN-NA patients

In the TARGET-NB cohort, the PRGs‐based signatures were constructed by multivariate Cox regression. The 5-gene model formula was as follows: Risk score = (0.3156) × GAL + (− 0.2791) × GFRA3 + (− 1.4725) × MARCKS + (1.2956) × PSMD13 + (− 0.8855) × ZNHIT3. Then, according to this formula, we calculated risk scores of each patient. Based on the median score, 143 HR-MYCN-NA neuroblastoma cases were classified into an ultra-high risk group (UHR, n = 71) and high-risk group (HR, n = 72). The risk scores, survival status and gene expression heatmap of these prognostic signatures were shown in Fig. 3A–C. Kaplan–Meier curve indicated that patients in the UHR group showed significantly poorer OS than those in the HR group (HR = 4.29, Fig. 3D). The median OS was 2.2 years for the ultra-high risk patients. The predictive ability and accuracy of the gene signature were measured according to ROC curve analysis, and the AUC of the signature predicting the 1-, 3- and 5-year OS rates were 0.74, 0.82 and 0.81, respectively, indicating that this prognostic model exhibited a good sensitivity and specificity (Fig. 3E).Fig. 3 Construction of the PRGs signature in the TARGET-NB cohort. A The distribution of risk score in HR-MYCN-NA patients. B Survival status and survival time of HR-MYCN-NA patients. C Heatmap of the prognostic-related genes expression; D Survival curves for the ultra-high risk and high-risk groups; E tROC analysis of risk scores to predict the OS in the TARGET-NB cohort

Validation of the prognostic signature for HR-MYCN-NA patients

To validate the robustness of the PRGs signature, firstly, GSE49710 including 83 high-risk MYCN-NA tumor samples were used. According to the median risk score, we divided cases into UHR (n = 41) and HR groups (n = 42). Consistent with the results derived from the TARGET-NB database, the Kaplan–Meier curve presented that patients in the UHR group exhibited significantly poorer OS than those in the HR group (HR = 2, 95%CI [1.1, 3.8], P = 0.033, Fig. 4D). The median OS was 3.7 years for the ultra-high risk patients. The risk scores, survival status and gene expression heatmap of these prognostic signatures were shown in Fig. 4A–C. The AUCs for 1-, 3- and 5-year OS were 0.67, 0.68 and 0.71, respectively (Fig. 4E).Fig. 4 Validation of the prognostic gene signature in the GSE49710cohort. A The distribution of risk score in HR- MYCN-NA patients. B Survival status and survival time of HR-MYCN-NA patients. C Heatmap of the prognostic-related genes expression. D Survival curves for the ultra-high risk and high-risk groups. E tROC analysis of risk scores to predict the OS in the GSE49710 cohort

Secondly, the GSE45547 database was used to validate the expression of five genes. A significantly higher expression of PSMD13, and lower expression of MARCKS and GFRA3 (P < 0.0001) were observed in HR-NB (stage IV or IVs) than low- or intermediate-risk NB (stage I–III) patients (P < 0.0001), which are consistent with the results derived from the TARGET-NB database. But there was no difference of gene GAL and ZNHIT3 expression (P > 0.05) (Fig. 5). GO analyses indicated the five genes were enriched in signaling receptor binding or activity (supplementary S3).Fig. 5 Validation of the expression distribution of GAL, PSMD13, MARCKS, GFRA3, and ZNHIT3 gene in different stage of NB patients in the GSE45547 cohort

Clinical values of prognostic signature for HR-MYCN-NA patients

Univariate and multivariate Cox regression analysis were applied to evaluate the independent prediction ability of prognostic signatures between the gene signature and other common clinical factors, including age, race, gender, histology type, ploidy, MKI, grade, and primary site of NB patients. Univariate Cox analysis indicated that ploidy (HR = 1.78, P = 0.03), tumor grade (HR = 2.52, P = 0.03), and PRGs signature-related risk group (HR = 4.28, P < 0.001), were markedly associated with OS. While the multivariate analysis, gene signature-related risk group (HR = 4.45, P < 0.001) and tumor grade (HR = 4.15, P = 0.02) were independent prognostic factors for HR-MYCN-NA patients (Fig. 6).Fig. 6 Univariate and multivariate Cox regression to identify independent risk factors for OS in HR-MYCN-NA patients

Establishment of a nomogram for HR-MYCN-NA patients

To broaden the clinical application and usability of the gene signatures, we integrated survival time, survival status and PRGs signature, which had most significant values in the multivariate Cox analysis. A nomogram to predict 1-, 3- and 5-year OS of HR-MYCN-NA patients was constructed, which can easily be used to calculate the expected survival rate of individual patients (Fig. 7A). For example, assuming tumor grade was differentiation, and the expression value of 5 genes in a tumor are GAL = 7.5, GFRA3 = 6.0, MARCKS = 10.0, PSMD13 = 7.9, ZNHIT3 = 7.0. According to the nomogram, the corresponding risk scores of these five genes are 16.08, 18.05, 35.50, 25.51, and 52.98 respectively. Then the total points were 148.11, and the probability of 1-, 3- and 5-year OS was 92%, 50% and 30%. In addition, the calibration curve was used to show the predictive ability and accuracy of the nomogram. The overall C-index of the nomogram model was 0.75(95%CI [0.69–0.80], P < 0.001), which demonstrated good agreement between the predicted and observed survival rates (Fig. 7B).Fig. 7 Construction of a nomogram based on the PRGs signature and tumor grade. A Nomogram for HR-MYCN-NA patients estimating the 1-, 3-, and 5-year OS. B Calibration curves of nomogram

Discussion

High-risk neuroblastoma is a heterogeneous disease; however, these patients have been treated in a similar strategy without further risk stratification [21, 22]. Therefore, accurately stratifying HR-NB patients into subtypes is crucial for individualized treatment. This study focus on high risk NB patients without MYCN amplification (HR-MYCN-NA), and successfully stratified these patients into two subgroups through combining both microarray and scRNA-sequence data, which provide information for new potential personalized treating neuroblastoma. To our knowledge, this is the first study aimed at screen for prognostic gene sets that related with survival of HR-MYCN-NA cases.

In this study, we integrated transcriptome datasets both from bulk microarray and single cell RNA-Seq, and finally selected 5 PRGs to construct a prognostic model, which successfully stratified high-risk MYCN-NA patients into ultra-high risk and high-risk groups (P < 0.0001). Moreover, a nomogram to predict the 1-, 3- and 5- year OS was established, and generated a high C-index 0.75. Previous studies also evaluated the significance of risk stratification models in predicting the prognosis of HR-NB patients. Nomogram only included clinical factors, such as MYCN status, LDH, presence of bone marrow metastases, has got an AUC = 0.63 [23]. Another study included INSS stage, MKI, or MKI, MYCN status to predict EFS, OS respectively, and the C-index of the nomogram was 0.62 and 0.65 respectively [24]. Moreover, with high-throughput technologies development, molecular characters has played a critical role in the prognosis of HR-NB. Multi-omics data integration revealed two prognostic subtypes in HR-NB, and generated high C-index (0.74 for EFS and 0.71 for OS) [25]. A nomogram combined with gene signature and clinical factors (age, tumor histology, and MYCN status) got a C-index 0.717 [26].Our study was consistent with these previous studies. The slightly differences of predictive performance may be due to the different parameters, and predictive models established by combining genomic feathers with traditional clinical indicators seems to have higher predictive performance.

The classification of high-risk patients is remained difficult. Especially the HR-MYCN-NA tumors has remain been one of the challenging subgroups for risk stratification [27]. Ohira’ research found that LDH ≥ 1400 U/L, age ≥ 18 months, and 1p/11q losses or 17q gain were strongly correlation with a poor survival in HR- MYCN-NA cases in Japanese patients [28]. While Moreno et al. found that MYCN status, LDH, and metastasis to bone marrow, were prognostic biomarkers in high-risk neuroblastoma; But age ≥ 5 years was not associated with prognosis in INRG project [23]. However, in our study, we found that age ≥ 2 years (P = 0.98) was not associated with OS, while gene signature (HR = 4.45, P < 0.001) and tumor grade (HR = 4.15, P = 0.02) were independent prognostic factors of HR-MYCN-NA patients. These may be difference in sample sizes and patients with different genetic backgrounds. The common feature of these studies is that they contain genome or transcriptome information, which indicate that combination of the tumor genomic/transcriptome features with clinical factors may be a promising predictor for HR-MYCN-NA patients, this was in accordance with previous findings from the INRG and COG projects [3, 6]. However, we couldn’t get the detailed genomic information from the public databases in the current study. And we will conduct the further study based on our own NB patients’ cohort in the future.

Two (GAL, PSMD13) of the 5 prognostic genes were highly expressed in the UHR group, while three (GFRA3, MARCKS, and ZNHIT3) of them are highly expressed in the HR group. Among these 5 genes, only GAL, GFRA3, and MARCKS have been found to play some roles in neuroblastoma prognosis [29, 30], while PSMD13 and ZNHIT3 have not been reported to be related with NB previously. There is some evidence that GAL may plays a critical role in the prognosis of high risk NB. Galanin (GAL) is a 29 to 30-amino acid peptide, and has a wide range of effects, especially in neuromodulatory and endocrine functions. The functions of GAL are mediated by three G protein-coupled receptors GALRs (GAL1R, GAL2R, GAL3R) [31]. Studies indicated that galanin and its receptors are involved in tumorigenesis, the invasion and migration of tumor cells, and related to tumor stage, metastasis, and recurrence in many types of tumors [32]. In neuroblastoma, GAL1R and GAL3R are highly expressed in NB tissues, while the presence of GAL2R is not as common as GAL1R [33]. Moreover, GAL2R mediates apoptosis in SH-SY5Y NB cells, but the GAL anti-proliferative potency was 100-fold higher of overexpressing GAL2R than overexpressing GAL1R in SH-SY5Y NB cells [34]. For the association of GAL expression and patients' prognosis, a Previous study found that GAL upregulated in MYCN-A compared to MYCN-NA tumors, and was a significant OS predictor [30]. More interestingly, we find that expression patterns of GAL in MYCN non-amplified tumors can further stratify patients into clinically significant subtypes. Other studies also found that high expression of GAL is associated with poor prognosis in colorectal cancer [35], endometrial cancer, and lung cancer (https://www.proteinatlas.org/ENSG00000069482-GAL/pathology). The above suggests that GAL and its receptors are important in the development and prognosis of NB tumors. For GFRA3, Rimas et al. found that GFRA3 often overexpressed paired with ALK in MYCN non-amplified stage 4 neuroblastoma [29]. Other studies found that high expression of GFRA3 is unfavorable in carcinoma [36]. For MARCKS, pre-clinical studies reported that MARCKS expression or phosphorylation in neuroblastoma cells may plays a role in neurotoxicity [37]. Clinical study indicated that MARCKS gene was associated with poor outcome in acute myeloid leukemia [38]. The function of the PSMD13 and ZNHIT3 have not been reported in neuroblastoma. The relationship of these genes and their roles in HR-MYCN-NA NB need to further elucidation.

With the development of high-throughput sequencing technology, an effective method to predict disease prognosis and find new therapeutic targets is comprehensively analysis of multi-omics data [25, 39, 40]. However, traditional tumor tissue transcriptome sequencing technology (bulk Seq) is based on multicellular gene sequences. Due to the high heterogeneity in tissues, it cannot identify the internal characteristics of different cell subsets. Recently, NB single-cell transcriptome sequencing has been reported and the identification of NB tumor cells from the single-cell transcriptome level provides new ideas and resources for the development of novel high-risk NB gene set markers [41]. To further stratify HR MYCN-NA patients, microarray and scRNA sequencing data were integrated and analyzed. Prognostic-related genes (PRGs) of HR-MYCN-NA patients in NB tumor cells were identified and successfully stratified these patients into clinical subgroups, which indicated that these genes may be potential targets that interfere with the functions of the prognostic cell types.

There are some limitations to this study. First, limited by the clinical information included in the public data sets, we were not able to establish the nomogram by combining the gene signature and clinical factors. Except for gene signature (HR = 4.45, P < 0.001), we only found one clinical factor (tumor grade, HR = 4.15, P = 0.02) was an independent prognostic factor. Second, we did not conduct experimental studies to further confirm our findings. The specific function of these identified prognostic genes and their potential mechanisms in HR MYCN-NA tumor progression need to be investigated by further experimental studies. We will perform the in vitro and/or in vivo assays to validate the 5 PRG genes and their mechanism for NB in our next study. Third, this study is retrospective, because we focus on HR-MYCN-NA tumors, although we used different cohorts as training and validation datasets, the sample size in this study was still relatively small. Therefore, prospective studies with large sample sizes are still necessary in future research. Despite these drawbacks, the combination of scRNA-Seq and microarray data, and the external validation of the findings by independent datasets provide a high level of confidence.

In the present study, a novel five-gene signature based on prognostic models was constructed to identify UHR subgroup from the HR-MYCN-NA cases. To our knowledge, this is the first study to further identify UHR subsets from the HR-MYCN-NA patients by combining bulk-Seq and scRNA-Seq data. The integrative classification of HR-MYCN-NA neuroblastoma may help clinicians better predict patients’ prognosis, and develop personalized treatment programs.

Supplementary Information

Supplementary file1

Acknowledgements

The authors are very thankful to the TARGET and GEO Program for providing the opening and high-quality resources for researchers. We are also grateful for the bioinformatics support provided by the website of Assistant for Clinical Bioinformatics.

Author contributions

Lin Zou and Jing-Bo Shao contributed to the study’s conception and design. Jin-Xia Wang analyzed the data and wrote the manuscript. Lin Zou and Jin-Xia Wang revised the manuscript. Zi-Yang Cao, Zi-Jun Yan, Hong-Yang Zhang contributed to the collection and assembly of data.

Funding

This study was supported by the Shanghai Sailing Program (22YF1437200), and partially supported by the Natural Science Foundation of China (82070167, 82270160).

Data availability

The data of this study are available in the TARGET and GEO Program.

Declarations

Ethics approval and consent to participate 

Because the TARGET and GEO database is publicly accessible worldwide, therefore, we did not provide the approval of an institutional review board in the current study.

Competing interests

The authors declare that they have no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

1. Almstedt E Elgendy R Hekmati N Integrative discovery of treatments for high-risk neuroblastoma Nature Commun 2020 11 1 71 10.1038/s41467-019-13817-8 31900415
Almstedt E, Elgendy R, Hekmati N, et al. Integrative discovery of treatments for high-risk neuroblastoma. Nature Commun. 2020;11(1):71.31900415
2. Newman EA Abdessalam S Aldrink JH Update on neuroblastoma J Pediatr Surg 2019 54 3 383 389 10.1016/j.jpedsurg.2018.09.004 30305231
Newman EA, Abdessalam S, Aldrink JH, et al. Update on neuroblastoma. J Pediatr Surg. 2019;54(3):383–9.30305231
3. Cohn SL Pearson AD London WB The International Neuroblastoma Risk Group (INRG) classification system: an INRG Task Force report J Clin Oncol 2009 27 2 289 297 10.1200/JCO.2008.16.6785 19047291
Cohn SL, Pearson AD, London WB, et al. The International Neuroblastoma Risk Group (INRG) classification system: an INRG Task Force report. J Clin Oncol. 2009;27(2):289–97.19047291
4. Shohet J Foster J Neuroblastoma BMJ (Clin Res ed) 2017 357 j1863 10.1136/bmj.j1863
Shohet J, Foster J. Neuroblastoma. BMJ (Clin Res ed). 2017;357: j1863.
5. Liang WH Federico SM London WB Tailoring therapy for children with neuroblastoma on the basis of risk group classification: past, present, and future JCO Clin Cancer Inform 2020 4 895 905 10.1200/CCI.20.00074 33058692
Liang WH, Federico SM, London WB, et al. Tailoring therapy for children with neuroblastoma on the basis of risk group classification: past, present, and future. JCO Clin Cancer Inform. 2020;4:895–905.33058692
6. Irwin MS Naranjo A Zhang FF Revised neuroblastoma risk classification system_ a report from the Children's Oncology Group J Clin Oncol 2021 39 29 3229 3241 10.1200/JCO.21.00278 34319759
Irwin MS, Naranjo A, Zhang FF, et al. Revised neuroblastoma risk classification system_ a report from the Children’s Oncology Group. J Clin Oncol. 2021;39(29):3229–41.34319759
7. Otte J, Dyberg C, Pepich A, Johnsen JI. MYCN function in neuroblastoma development. Front Oncol 2021; 10.
8. Wei JS Kuznetsov IB Zhang S Clinically relevant cytotoxic immune cell signatures and clonal expansion of T-cell receptors in high-risk MYCN-not-amplified human neuroblastoma Clin Cancer Res 2018 24 22 5673 5684 10.1158/1078-0432.CCR-18-0599 29784674
Wei JS, Kuznetsov IB, Zhang S, et al. Clinically relevant cytotoxic immune cell signatures and clonal expansion of T-cell receptors in high-risk MYCN-not-amplified human neuroblastoma. Clin Cancer Res. 2018;24(22):5673–84.29784674
9. Schwab M Westermann F Hero B Berthold F Neuroblastoma: biology and molecular and chromosomal pathology Lancet Oncol 2003 4 8 472 480 10.1016/S1470-2045(03)01166-5 12901961
Schwab M, Westermann F, Hero B, Berthold F. Neuroblastoma: biology and molecular and chromosomal pathology. Lancet Oncol. 2003;4(8):472–80.12901961
10. Bénard J Raguénez G Kauffmann A MYCN-non-amplified metastatic neuroblastoma with good prognosis and spontaneous regression: a molecular portrait of stage 4S Mol Oncol 2008 2 3 261 271 10.1016/j.molonc.2008.07.002 19383347
Bénard J, Raguénez G, Kauffmann A, et al. MYCN-non-amplified metastatic neuroblastoma with good prognosis and spontaneous regression: a molecular portrait of stage 4S. Mol Oncol. 2008;2(3):261–71.19383347
11. Schmitt-Hoffner F van Rijn S Toprak UH FOXR2 stabilizes MYCN protein and identifies non-MYCN-amplified neuroblastoma patients with unfavorable outcome J Clin Oncol 2021 39 29 3217 3228 10.1200/JCO.20.02540 34110923
Schmitt-Hoffner F, van Rijn S, Toprak UH, et al. FOXR2 stabilizes MYCN protein and identifies non-MYCN-amplified neuroblastoma patients with unfavorable outcome. J Clin Oncol. 2021;39(29):3217–28.34110923
12. Asgharzadeh S Pique-Regi R Sposto R Prognostic significance of gene expression profiles of metastatic neuroblastomas lacking MYCN gene amplification JNCI-J Natl Cancer I 2006 98 17 1193 1203 10.1093/jnci/djj330
Asgharzadeh S, Pique-Regi R, Sposto R, et al. Prognostic significance of gene expression profiles of metastatic neuroblastomas lacking MYCN gene amplification. JNCI-J Natl Cancer I. 2006;98(17):1193–203.
13. Bellini A Pötschger U Bernard V Frequency and prognostic impact of ALK amplifications and mutations in the European neuroblastoma study group (SIOPEN) high-risk neuroblastoma trial (HR-NBL1) J Clin Oncol 2021 39 30 3377 3390 10.1200/JCO.21.00086 34115544
Bellini A, Pötschger U, Bernard V, et al. Frequency and prognostic impact of ALK amplifications and mutations in the European neuroblastoma study group (SIOPEN) high-risk neuroblastoma trial (HR-NBL1). J Clin Oncol. 2021;39(30):3377–90.34115544
14. Rosswog C Fassunke J Ernst A Genomic ALK alterations in primary and relapsed neuroblastoma Br J Cancer 2023 128 8 1559 1571 10.1038/s41416-023-02208-y 36807339
Rosswog C, Fassunke J, Ernst A, et al. Genomic ALK alterations in primary and relapsed neuroblastoma. Br J Cancer. 2023;128(8):1559–71.36807339
15. Valentijn LJ Koster J Zwijnenburg DA TERT rearrangements are frequent in neuroblastoma and identify aggressive tumors Nat Genet 2015 47 12 1411 1414 10.1038/ng.3438 26523776
Valentijn LJ, Koster J, Zwijnenburg DA, et al. TERT rearrangements are frequent in neuroblastoma and identify aggressive tumors. Nat Genet. 2015;47(12):1411–4.26523776
16. Peifer M Hertwig F Roels F Telomerase activation by genomic rearrangements in high-risk neuroblastoma Nature 2015 526 7575 700 704 10.1038/nature14980 26466568
Peifer M, Hertwig F, Roels F, et al. Telomerase activation by genomic rearrangements in high-risk neuroblastoma. Nature. 2015;526(7575):700–4.26466568
17. Huang M Zeki J Sumarsono N Epigenetic targeting of TERT-associated gene expression signature in human neuroblastoma with TERT overexpression Cancer Res 2020 80 5 1024 1035 10.1158/0008-5472.CAN-19-2560 31900258
Huang M, Zeki J, Sumarsono N, et al. Epigenetic targeting of TERT-associated gene expression signature in human neuroblastoma with TERT overexpression. Cancer Res. 2020;80(5):1024–35.31900258
18. Wei M Ye M Dong K Dong R Mass screening for neuroblastoma in infants Discov Med 2020 30 160 63 70 33382962
Wei M, Ye M, Dong K, Dong R. Mass screening for neuroblastoma in infants. Discov Med. 2020;30(160):63–70.33382962
19. Simon N Friedman JH Hastie T Tibshirani R Regularization paths for Cox's proportional hazards model via coordinate descent J Stat Softw 2011 39 5 1 13 10.18637/jss.v039.i05 27065756
Simon N, Friedman JH, Hastie T, Tibshirani R. Regularization paths for Cox’s proportional hazards model via coordinate descent. J Stat Softw. 2011;39(5):1–13.27065756
20. Heagerty PJ Lumley T Pepe MS Time-dependent ROC curves for censored survival data and a diagnostic marker Biometrics 2000 56 2 337 344 10.1111/j.0006-341X.2000.00337.x 10877287
Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. 2000;56(2):337–44.10877287
21. Qiu B Matthay KK Advancing therapy for neuroblastoma Nat Rev Clin Oncol 2022 19 8 515 533 10.1038/s41571-022-00643-z 35614230
Qiu B, Matthay KK. Advancing therapy for neuroblastoma. Nat Rev Clin Oncol. 2022;19(8):515–33.35614230
22. Tolbert VP Matthay KK Neuroblastoma: clinical and biological approach to risk stratification and treatment Cell Tissue Res 2018 372 2 195 209 10.1007/s00441-018-2821-2 29572647
Tolbert VP, Matthay KK. Neuroblastoma: clinical and biological approach to risk stratification and treatment. Cell Tissue Res. 2018;372(2):195–209.29572647
23. Moreno L Guo D Irwin MS A nomogram of clinical and biologic factors to predict survival in children newly diagnosed with high-risk neuroblastoma: An International Neuroblastoma Risk Group project Pediatr Blood Cancer 2021 68 3 e28794 10.1002/pbc.28794 33205902
Moreno L, Guo D, Irwin MS, et al. A nomogram of clinical and biologic factors to predict survival in children newly diagnosed with high-risk neuroblastoma: An International Neuroblastoma Risk Group project. Pediatr Blood Cancer. 2021;68(3): e28794.33205902
24. Feng L Kan Y Wang W Development and validation of a nomogram for predicting survival in intermediate- and high-risk neuroblastoma of the Children's Oncology Group risk stratification J Cancer Res Clin 2023 149 18 16377 16390 10.1007/s00432-023-05398-3
Feng L, Kan Y, Wang W, et al. Development and validation of a nomogram for predicting survival in intermediate- and high-risk neuroblastoma of the Children’s Oncology Group risk stratification. J Cancer Res Clin. 2023;149(18):16377–90.
25. Zhang L Lv C Jin Y Deep learning-based multi-omics data integration reveals two prognostic subtypes in high-risk neuroblastoma Front Genet 2018 9 477 10.3389/fgene.2018.00477 30405689
Zhang L, Lv C, Jin Y, et al. Deep learning-based multi-omics data integration reveals two prognostic subtypes in high-risk neuroblastoma. Front Genet. 2018;9:477.30405689
26. Xia Y Li X Tian X Zhao Q Identification of a five-gene signature derived from MYCN amplification and establishment of a nomogram for predicting the prognosis of neuroblastoma Front Mol Biosci 2021 8 769661 10.3389/fmolb.2021.769661 34950701
Xia Y, Li X, Tian X, Zhao Q. Identification of a five-gene signature derived from MYCN amplification and establishment of a nomogram for predicting the prognosis of neuroblastoma. Front Mol Biosci. 2021;8: 769661.34950701
27. Matthay KK Maris JM Schleiermacher G Neuroblastoma Nat Rev Dis Prim 2016 2 16078 10.1038/nrdp.2016.78 27830764
Matthay KK, Maris JM, Schleiermacher G, et al. Neuroblastoma. Nat Rev Dis Prim. 2016;2:16078.27830764
28. Ohira M Nakamura Y Takimoto T Retrospective analysis of INRG clinical and genomic factors for 605 neuroblastomas in Japan: A report from the Japan Children's Cancer Group Neuroblastoma Committee (JCCG-JNBSG) Biomolecules 2021 12 1 18 10.3390/biom12010018 35053166
Ohira M, Nakamura Y, Takimoto T, et al. Retrospective analysis of INRG clinical and genomic factors for 605 neuroblastomas in Japan: A report from the Japan Children’s Cancer Group Neuroblastoma Committee (JCCG-JNBSG). Biomolecules. 2021;12(1):18.35053166
29. Orentas RJ Sindiri S Duris C Paired expression analysis of tumor cell surface antigens Front Oncol 2017 7 173 10.3389/fonc.2017.00173 28871274
Orentas RJ, Sindiri S, Duris C, et al. Paired expression analysis of tumor cell surface antigens. Front Oncol. 2017;7:173.28871274
30. Liu Z Grant CN Sun L Miller BA Spiegelman VS Wang HG Expression patterns of immune genes reveal heterogeneous subtypes of high-risk neuroblastoma Cancers (Basel) 2020 12 7 1739 10.3390/cancers12071739 32629858
Liu Z, Grant CN, Sun L, Miller BA, Spiegelman VS, Wang HG. Expression patterns of immune genes reveal heterogeneous subtypes of high-risk neuroblastoma. Cancers (Basel). 2020;12(7):1739.32629858
31. Jurkowski W Yazdi S Elofsson A Ligand binding properties of human galanin receptors Mol Member Biol 2013 30 2 206 216 10.3109/09687688.2012.750384
Jurkowski W, Yazdi S, Elofsson A. Ligand binding properties of human galanin receptors. Mol Member Biol. 2013;30(2):206–16.
32. Sánchez ML Coveñas R The galaninergic system: a target for cancer treatment Cancers 2022 14 15 3755 10.3390/cancers14153755 35954419
Sánchez ML, Coveñas R. The galaninergic system: a target for cancer treatment. Cancers. 2022;14(15):3755.35954419
33. Perel Y Amrein L Dobremez E Galanin and galanin receptor expression in neuroblastic tumours: correlation with their differentiation status Br J Cancer 2002 86 1 117 122 10.1038/sj.bjc.6600019 11857022
Perel Y, Amrein L, Dobremez E, et al. Galanin and galanin receptor expression in neuroblastic tumours: correlation with their differentiation status. Br J Cancer. 2002;86(1):117–22.11857022
34. Berger A Lang R Moritz K Galanin receptor subtype GalR2 mediates apoptosis in SH-SY5Y neuroblastoma cells Endocrinology 2004 145 2 500 507 10.1210/en.2003-0649 14592962
Berger A, Lang R, Moritz K, et al. Galanin receptor subtype GalR2 mediates apoptosis in SH-SY5Y neuroblastoma cells. Endocrinology. 2004;145(2):500–7.14592962
35. Nagayoshi K Ueki T Tashiro K Galanin plays an important role in cancer invasiveness and is associated with poor prognosis in stage II colorectal cancer Oncol Rep 2015 33 2 539 546 10.3892/or.2014.3660 25504183
Nagayoshi K, Ueki T, Tashiro K, et al. Galanin plays an important role in cancer invasiveness and is associated with poor prognosis in stage II colorectal cancer. Oncol Rep. 2015;33(2):539–46.25504183
36. Wu ZS Pandey V Wu WY Prognostic significance of the expression of GFRα1, GFRα3 and syndecan-3, proteins binding ARTEMIN, in mammary carcinoma BMC Cancer 2013 13 34 10.1186/1471-2407-13-34 23351331
Wu ZS, Pandey V, Wu WY, et al. Prognostic significance of the expression of GFRα1, GFRα3 and syndecan-3, proteins binding ARTEMIN, in mammary carcinoma. BMC Cancer. 2013;13:34.23351331
37. Shiraishi M Hangai M Yamamoto M Alteration in MARCKS phosphorylation and expression by methylmercury in SH-SY5Y cells and rat brain Environ Toxicol Pharmacol 2014 37 3 1256 1263 10.1016/j.etap.2014.04.025 24835554
Shiraishi M, Hangai M, Yamamoto M, et al. Alteration in MARCKS phosphorylation and expression by methylmercury in SH-SY5Y cells and rat brain. Environ Toxicol Pharmacol. 2014;37(3):1256–63.24835554
38. Lebedev TD Vagapova ER Popenko VI Two receptors, two isoforms, two cancers: comprehensive analysis of KIT and TrkA expression in neuroblastoma and acute myeloid leukemia Front Oncol 2019 9 1046 10.3389/fonc.2019.01046 31681584
Lebedev TD, Vagapova ER, Popenko VI, et al. Two receptors, two isoforms, two cancers: comprehensive analysis of KIT and TrkA expression in neuroblastoma and acute myeloid leukemia. Front Oncol. 2019;9:1046.31681584
39. Duffy DJ Krstic A Halasz M Integrative omics reveals MYCN as a global suppressor of cellular signalling and enables network-based therapeutic target discovery in neuroblastoma Oncotarget 2015 6 41 43182 43201 10.18632/oncotarget.6568 26673823
Duffy DJ, Krstic A, Halasz M, et al. Integrative omics reveals MYCN as a global suppressor of cellular signalling and enables network-based therapeutic target discovery in neuroblastoma. Oncotarget. 2015;6(41):43182–201.26673823
40. Meng X Feng C Fang E Combined analysis of RNA-sequence and microarray data reveals effective metabolism-based prognostic signature for neuroblastoma J Cell Mol Med 2020 24 18 10367 10381 10.1111/jcmm.15650 32683778
Meng X, Feng C, Fang E, et al. Combined analysis of RNA-sequence and microarray data reveals effective metabolism-based prognostic signature for neuroblastoma. J Cell Mol Med. 2020;24(18):10367–81.32683778
41. Dong R, Yang R, Zhan Y, et al. Single-cell characterization of malignant phenotypes and developmental trajectories of adrenal neuroblastoma. Cancer Cell 2020; 38(5):716–33 e6.
