Neuroblastoma (NB) is the most common extracranial solid tumor in children, with the highest incidence and mortality in infancy, accounting for approximately 8–10% of pediatric malignancies (Li et al., 2008; Marshall et al., 2014). It is a highly heterogeneous disease characterized by diverse clinical manifestations ranging from spontaneous regression to progression with therapy resistance. Based on clinical and molecular patterns, NB patients can be classified into low-risk, intermediate-risk, and high-risk groups (Cohn et al., 2009). Although the prognosis of NB has significantly improved over the past decade with the development of multiple treatment methods and the application of immunotherapy targeting NB-specific antigens, such as GD2 (Casey and Cheung, 2020), less than 40% of high-risk NB patients achieve long-term survival (Pinto et al., 2015). Given the significantly divergent outcomes among these patients, an effective prediction model is needed so that effective management and therapy can be tailored to diverse subgroups of NB patients.

There is increasing evidence that several genomic alterations, such as gain of chromosome 17q, hemizygous deletions of 1p and 11q, and MYCN gene amplification, are powerful prognostic markers that are strongly associated with clinical outcomes (Janoueix-Lerosey et al., 2009; Maris, 2010; Schleiermacher et al., 2012). Among these factors, MYCN amplification is one of the earliest biomarkers discovered in NB and is still considered one of the most reliable predictors of aggressive clinical behavior (Park et al., 2013; Campbell et al., 2017). Augmented MYCN is detected in 20–30% of all NB cases, and the overall survival rate for this group remains less than 50% (Huang and Weiss, 2013). Although MYCN plays a pivotal role in NB development, its application as a therapeutic target remains challenging. Although much research has been done on targeting MYCN-related factors, new target genes are still needed to effectively treat this aggressive and heterogeneous disease (Puissant et al., 2013; Gustafson et al., 2014; Huang et al., 2019). Due to the rapid development of high-throughput sequencing, multigene signatures have been shown in many studies to be more effective at predicting prognosis than conventional biomarkers. The genes discovered by bioinformatic approaches can also provide insights into cancer progression and identify pathway profiles involved in tumorigenesis. Several NB prognostic models have been made based on sequencing data (Vermeulen et al., 2009; De Preter et al., 2010; He et al., 2020; Wang et al., 2020). A deeper dive into public data sets may reveal other genes associated with prognosis and new potential therapeutic targets.

Herein, we aimed to build a prognostic signature that could predict the prognosis of NB patients more accurately than traditional clinical factors. We systematically analyzed multiple public NB data sets to identify differentially expressed genes (DEGs) based on MYCN amplification. Then, we successfully established a prognostic five-gene risk prediction model that can accurately predict the outcome of NB. Simultaneously, we performed preliminary validation of the pivotal DEGs using clinical samples. A prognostic nomogram was finally designed based on the gene signature and clinical factors, and this nomogram will help physicians develop more customized treatment plans for different subgroups of NB patients.

In the current study, we developed a novel multigene signature that showed better predictive power than conventional biomarkers. At present, there is still an urgent need for new effective stratification systems and treatment targets for NB. Although it has been well established that MYCN amplification affects the prognosis of NB (Brodeur et al., 1984; Seeger et al., 1985), it remains challenging to target this oncogene in clinical treatment. The development of the sequencing technology allows us to deeply explore the molecular mechanism by which MYCN amplification affects prognosis at the whole-genome level. Relevant NB data were obtained from multiple data sets, and patients were grouped based on MYCN status. We identified 319 DEGs, including 51 upregulated and 268 downregulated genes. GO and KEGG enrichment analyses of these DEGs were performed to further investigate the molecular mechanism of NB. The DEGs were shown to be significantly enriched in biological processes linked to the occurrence and development of NB, such as glial cell differentiation, the catecholamine biosynthetic process, and sympathetic nervous system development. The disialoganglioside GD2, a unique carbohydrate antigen, is one of the most attractive immune surface targets for NB. GD2 is uniformly expressed by NB cells, which makes GD2 an ideal therapeutic target for NB (Suzuki and Cheung, 2015; Sait and Modak, 2017). According to GO analysis, the DEGs in our study were significantly enriched in molecular functions involved in carbohydrate binding. Among these DEGs, there may be potential regulatory genes associated with GD2 immunotherapy. KEGG pathway analysis showed that the DEGs were mainly enriched in cancer-related pathways, such as the amino acid metabolism and cell adhesion processes. Cancer cells require an abundant supply of amino acids to sustain their growth (Vettore et al., 2020), and an increased amino acid metabolism has been observed in chemoresistant NB cell lines (Gunda et al., 2020). The cancer metastatic cascade is dependent on the loss of adhesion between cells, resulting in the dissociation of the cell from the primary tumor (Martin and Jiang, 2009). The subsequent PPI network results further confirmed that these DEGs participated in the NB cell differentiation and metabolism.

We further screened 35 DEGs that significantly affected survival. Afterward, a risk prediction model consisting of five genes (GDPD5, NXPH1, AHI1, CPLX3, and SPAG6) was established. The upregulated genes CPLX3 and SPAG6 appeared to be associated with poor survival, whereas the downregulated genes GDPD5, NXPH1, and AHI1 were identified as tumor suppressors. The outcome was significantly different between NB patients with high and low risk scores. Additionally, the correlation between NB clinical features and the five-gene signature was investigated. Many known risk factors for NB that affect prognosis were shown to be related to a high risk score. Furthermore, the risk score for patients in each INSS stage increased from stage 1 to stage 4, except for stage 4S, which often has a favorable outcome (Kawano et al., 2021). Collectively, these findings suggest that patients with high risk scores exhibit more aggressive features than those with low risk scores. Subsequent external validation further confirmed that the five-gene signature could accurately predict survival in NB patients. To assess its clinical efficacy, we compared the gene signature with INSS staging and MYCN status, which are commonly used in predicting outcome in the clinic. The AUC and DCA curves also confirmed that the five-gene signature had a better prognosis prediction than the INSS stage and MYCN status, both in the short term and in the long term. In the MYCN nonamplification group, the gene model was still able to provide a good indication of the prognosis, which may aid physicians in evaluating prognosis and optimizing treatment options for this cohort of patients.

Due to the complex biological and clinical characteristics of NB, conventional staging systems may be difficult to utilize for calculating tumor risk in an accurate manner. To determine the treatment plan and predict the overall survival rate, doctors and researchers need a comprehensive prognostic evaluation system. Nomograms are being increasingly assessed in medical studies and are often used in the prediction of prognosis (Balachandran et al., 2015; Semenkovich et al., 2021). We integrated our five-gene model and other well-known clinicopathological parameters to create a nomogram. The prognosis for the children diagnosed with NB before 18 months is generally better than those diagnosed later (Vo et al., 2014). In the nomogram, age remained a crucial factor and ranked second after the five-gene model. Patients were divided into low-, medium-, and high-risk groups based on their total nomogram scores. The KM curves showed that the nomogram had excellent discriminative ability. According to the calibration curves, the predictions made by the nomogram were consistent with the observed data. The nomogram offers a visualized scoring system to facilitate medical decision making, which gives high-risk patients a chance to obtain more accurate treatments. Using it, physicians could predict the overall survival for each patient and make an individualized treatment plan. Patients in the high-risk group should be given special attention and more intensive treatment, while excessive treatment could be avoided for patients in the low–medium-risk groups.

Immunotherapy is currently a key hotspot in the field of cancer research. However, due to the difference in immunogenicity from adult tumors, targeted drugs, including immune checkpoint inhibitors, have limited applications in pediatric tumors (Wienke et al., 2021). In recent years, although new targeted therapies such as anti-GD2 therapy have exhibited great promise in treating NB, the limited efficacy and considerable toxicity make it still urgent to expand therapeutic targets and identify groups with potential benefits for NB immunotherapy (Sait and Modak, 2017; Park and Cheung, 2020). In our study, tumor immunological characteristics were assessed by the ESTIMATE algorithm. Higher immune scores and ESTIMATE scores were observed in patients with low risk scores. In light of the difference in overall survival between these two groups, more immune components in the tumor microenvironment may contribute to the better outcome of low-risk patients. In addition, the results calculated by the CIBERSORT algorithm indicated a significantly higher presence of T follicular helper cells (TFh) and activated NK cells in low-risk patients, both of which have been demonstrated to have protective roles in cancer. TFh plays an important role in humoral immunity (Eivazi et al., 2016). A previous study in NB showed that γδ TFh cells might facilitate the maturation of B cells and the production of antibodies (Mou et al., 2017). In studies of colorectal and breast cancer, infiltrating TFh cells have also been found to be positively correlated with prognosis (Bindea et al., 2013; Gu-Trantien et al., 2013). Similarly, a reduced degree of NK cell infiltration in tumor tissues has been found to be associated with a worse outcome in some patient cohorts (Albertsson et al., 2003; Mandal et al., 2016). NK cells play a central role in cancer immunotherapy (Shimasaki et al., 2020). As mentioned above, the application of anti-GD2 immunotherapy in recent years has significantly improved the survival of high-risk NB patients. NK cell-mediated antibody-dependent cellular cytotoxicity (ADCC) is a potent mechanism of anti-GD2 immunotherapy for NB. Anti-GD2 monoclonal antibodies coadministered with NK cell infusion have achieved promising results in several clinical trials. These studies demonstrated that NK cell infusion enhances ADCC directed by the antibody (Federico et al., 2017; Modak et al., 2018). Dinutuximab-mediated autologous killing of NB cells could be enhanced by plasmacytoid dendritic cell-induced NK cell activation (Belounis et al., 2020). Recent findings have shown that exosomes derived from human NK cells contain tumor suppressor miRNAs that are cytotoxic to NB cells and inhibit tumor cell immune escape, suggesting that NK cell-derived exosomes could be used as a complement to NK cell-based immunotherapy (Neviani et al., 2019). Except for these two types of cells, no differences were found between the two groups for other important immune cells in the training set, such as tumor-associated macrophages and CD8⁺ T cells. Interestingly, we obtained different results in the other two validation sets, suggesting that the different data sets selected might have been confounding. Unlike adult cancers, most pediatric solid tumors are considered immunologically cold, characterized by low tumor mutation burden, limited T-cell infiltration, and a lack of responsiveness to immune checkpoint inhibitors (Grobner et al., 2018; Wienke et al., 2021). In light of these differences between childhood tumors and adult tumors, a better understanding of NB immunobiology is needed, which will improve future immunotherapy in NB.

Among the five hub genes, AHI1, NXPH1, and CPLX3 have not been reported in any previous studies on NB. The oncogenic role of AHI1 has been identified in human leukemia and Sezary syndrome (Ringrose et al., 2006; Chen et al., 2013). Mutations in the AHI1 gene could lead to Joubert syndrome, which is a rare genetic disorder characterized by the underdevelopment of cerebellar vermis (Alvarez Retuerto et al., 2008). Moreover, it is also involved in autism (Alvarez Retuerto et al., 2008) and was considered to be a susceptibility gene for schizophrenia (Torri et al., 2010). It seems that AHI1, which has been associated with these brain developmental disorders, has a crucial role in brain growth during early life. As a specific neuronal cell surface protein, NXPH1 is essential for trans-synaptic activation (Missler et al., 2003). By using a genome-wide methylation screen, Faryna et al. found that the methylation level of NXPH1 in normal tissues was significantly higher than that in breast cancer samples, indicating it as a potential useful biomarker for tumor diagnosis (Faryna et al., 2012). In a study of intraductal papillary mucinous neoplasms (IPMNs), NXPH1 was more likely to be methylated in low-grade dysplasia than in high-grade dysplasia (Hong et al., 2012). NXPH1 was incorporated in a 10-gene signature that predicted biochemical recurrence for prostate cancer, and its expression level was upregulated in patients with a Gleason score ≥7 (Wu et al., 2020). In contrast, there was a negative correlation between lymph node metastasis and NXPH1 in pancreatic cancer (Jin and Tsai, 2016). GDPD5 is located on chromosome 11q13, and the presence of an 11q deletion is usually associated with an advanced stage and worse outcome in NB (Jiang et al., 2011). It belongs to the glycerophosphodiester phosphodiesterase family, which is critical for the metabolism of glycerol. In our study, we discovered that overexpression of GDPD5 is associated with better survival, which is in accordance with the result of a previous report in which high GDPD5 expression was found to be strongly correlated with a favorable outcome in neuroblastoma, while low GDPD5 expression was associated with a poor outcome. In addition, the study further revealed that GDPD5 induced the differentiation and inhibited the motility of NB cells via multiple mechanisms (Matas-Rico et al., 2016). GDPD5 was found to be overexpressed in colorectal cancer (CRC) and to promote metastasis and chemoresistance. Additionally, silencing GDPD5 in 5-Fu-resistant CRC cells decreased epithelial-to-mesenchymal transition (EMT) and cell invasion, both of which are essential for CRC metastasis (Feng et al., 2018). Cao et al. silenced GDPD5 in breast cancer cell lines and found a significant decrease in tumor cell viability, migration, and invasion (Cao et al., 2016). SPAG6 is identified as a cancer-testis antigen that regulates multiple functions in various cell types (Zheng et al., 2019). A variety of human cancers have been shown to be associated with SPAG6. A study on the differential analysis of genome-wide methylation of NB cell lines reported that SPAG6 was a possible target of the CpG island methylator phenotype (CIMP), and CIMP was demonstrated to be negatively related to the survival of NB patients (Abe et al., 2005; Abe et al., 2008). However, the role of SPAG6 in NB progression requires further study. CPLX3 is specifically localized to retinal ribbon synapses. It is a key regulator of transmitter release at retinal ribbon synapses (McMahon et al., 1995; Mortensen et al., 2016). At present, there is no research on CPLX3 and its role in tumors. Our results suggest that it may promote the development of NB and affect prognosis. Therefore, the specific role of CPLX3 in NB needs to be further studied. Our research has discovered five genes closely related to MYCN. Some of them have been identified to be involved in development of NB, and some have been studied in other malignant tumors. Next, we will further investigate the relationship between these five pivotal genes and MYCN as well as their roles in the development of NB.

However, there were several imitations to our study. First, since the TARGET-NBL gene expression data set used as the training set mainly contained NB cases in stage 1 and stage 4, we did not include INSS staging in the nomogram. This may limit the predictive power of the nomogram. Second, the bioinformatic analysis of this study was based on public data sets, and the tissues used for verification were retrospectively collected. Thus, an inherent bias in case selection may have influenced the findings.

Overall, we developed a novel five-gene signature and successfully built a prognostic nomogram that accurately predicts the survival of NB patients. It will simplify the stratification of NB patients and provide some guidance for individualized treatment. The five hub genes obtained in this study were based on integrated analyses of multiple data sets, which provides a very high level of reliability. Prior to this study, three of the identified prognostic genes had not been associated with NB. The mechanisms of these genes and their roles in targeted therapy for NB remain to be investigated.