Identification of MYCN non-amplified neuroblastoma subgroups points towards molecular signatures for precision prognosis and therapy stratification

Background Despite the extensive study of MYCN-amplified neuroblastomas, there is a significant unmet clinical need in MYCN non-amplified cases. In particular, the extent of heterogeneity within the MYCN non-amplified population is unknown. Methods A total of 1566 samples from 16 datasets were identified in Gene Expression Omnibus (GEO) and ArrayExpress. Characterisation of the subtypes was analysed by ConsensusClusterPlus. Independent predictors for subgrouping were constructed from the single sample predictor based on the multiclassPairs package. Findings were verified using immunohistochemistry and CIBERSORTx analysis. Results We demonstrate that MYCN non-amplified neuroblastomas are heterogeneous and can be classified into 3 subgroups based on their transcriptional signatures. Within these groups, subgroup_2 has the worst prognosis and this group shows a ‘MYCN’ signature that is potentially induced by the overexpression of Aurora Kinase A (AURKA); whilst subgroup_3 is characterised by an ‘inflamed’ gene signature. The clinical implications of this subtype classification are significant, as each subtype demonstrates a unique prognosis and vulnerability to investigational therapies. A total of 420 genes were identified as independent subgroup predictors with average balanced accuracy of 0.93 and 0.84 for train and test datasets, respectively. Conclusion We propose that transcriptional subtyping may enhance precision prognosis and therapy stratification for patients with MYCN non-amplified neuroblastomas.


INTRODUCTION
Neuroblastoma is the most common extra-cranial solid tumour in children, representing 6-10% of all childhood cancers [1].It is an embryonic tumour arising from precursor cells in the sympathetic nervous system and adrenal medulla [2], with a median age of diagnosis of 18 months [3].It can also be present in the neck, chest, abdomen, or pelvis [4].Neuroblastoma is a highly heterogeneous disease, with clinical behaviour ranging from spontaneous regression to drug resistance and metastasis ultimately resulting in death [5].The prognosis of the disease is poor with a 5-year overall survival of approximately 20%, despite more aggressive therapies [6].As a result, risk stratification and personalised treatment approaches in neuroblastomas are urgently needed.
The International Neuroblastoma Risk Group Staging System (INRGSS) defines the high-risk group to include patients with MYCN-amplified tumours and patients > 18 months old with metastatic tumours [7].N-MYC is a key regulator of transcription, which activates genes that affect cancer development.It is widely involved in various pathological processes of neuroblastoma including cell growth [8], apoptosis [9], differentiation [10], angiogenesis [11], tumour invasion, and metastasis [12].
MYCN amplification was identified as the first independent prognostic factor indicating adverse clinical outcomes in neuroblastomas [13,14], which is observed in approximately 20% of cases [15] and accounts for about 40% of high-risk neuroblastomas [16].Despite the extensive study of MYCN-amplified neuroblastomas, there is a significant unmet clinical need in MYCN non-amplified cases.In particular, the extent of heterogeneity within the MYCN non-amplified population is unknown.
Here, we investigated whether transcriptional subtyping of MYCN non-amplified neuroblastomas can identify molecular subtypes with discrete prognosis and therapeutic vulnerabilities.Our analysis suggested that MYCN non-amplified neuroblastomas were heterogeneous and could be classified into 3 subgroups based on their transcriptional profiling.Within them, subgroup 2 had the worst prognosis and this group had a 'MYCN' signature that was potentially induced by the overexpression of Aurora Kinase A (AURKA); whilst subgroup 3 was accompanied by an 'inflamed' gene signature.We propose that transcriptional subtyping may enhance precision prognosis and therapy stratification for patients with MYCN non-amplified neuroblastomas.

RESULTS
Characterisation of molecular subtypes in MYCN nonamplified neuroblastomas Following quality control and eliminating duplicates (Supplementary Figs. 1 and 2; details provided in the Supplementary Methods), a total of 1566 samples from 16 datasets were identified in GEO (Gene Expression Omnibus) and ArrayExpress, in which 313 cases are with MYCN gene amplification (MYCN-AMP) and 1253 cases MYCN non-amplified (MYCN-normal) (Fig. 1a; Supplementary Table 1).Following the removal of batch effects (Supplementary Fig. 3a), 2 clear clusters corresponding to MYCN-AMP and MYCNnormal neuroblastomas, respectively, were visualised using principal component analysis (PCA) (Supplementary Fig. 3b).Samples in the MYCN-normal group (n = 1253) were further randomly divided into a train and a test group with a 7:3 ratio, containing 878 and 375 cases, respectively (Fig. 1a).
Further clinical characterisation of these subtypes identified key distinguishing features.Patients within subgroup 2 were frequently observed in the advanced neuroblastomas according to the International Neuroblastoma Staging System (INSS) and in those defined as 'high risk' [7] (Fig. 2a, b; Supplementary Fig. 4b).We then analysed their overall survival together with MYCN-AMP cases.Patients with MYCN amplification had the worst prognosis (Fig. 2c, d; Supplementary Fig. 4c).Importantly, there was a high degree of variability for overall survival among MYCN non-amplified cases, in which subgroup 2 was associated with a poor prognosis, followed by subgroup 3; while patients within subgroup 1 had the most favourable outcomes.These observations were consistent in both train and test cohorts.In addition, the molecular subtype classification was a strong independent predictor of mortality including in multivariate analysis with the risk classification that uses commonly measured clinical variables to predict mortality in neuroblastomas [7].Using subgroup 1 as a reference, the hazard ratio (HR) and 95% confidence interval (CI) for subgroups 2 and 3 were 20.2 (4.8−85) and 9.2 (2.1−40), respectively (Fig. 2e).Similar results were obtained using univariate or multivariate cox regression analysis with age and INSS stages in MYCN non-amplified neuroblastomas (Supplementary Table 3).A comprehensive multivariate analysis also revealed our subgroups to be independent of genomic features such as 1p, 11q, and 17q (Supplementary Fig. 4d-f).Impressively, the molecular subtype classification alone outperformed INSS stages (Fig. 2f) and shows a comparable prediction accuracy as the risk classification (Supplementary Fig. 4g).
Overall, subgroup 2 and subgroup 3 (to a lesser extent) were associated with poor survival in MYCN non-amplified neuroblastomas, suggesting fundamentally different mechanisms leading to an advanced disease.
Subgroup 2 shows a 'MYCN' signature, potentially induced by Aurora Kinase A (AURKA) overexpression Our above analysis suggested that mechanisms other than gene amplification induce N-MYC activity in subgroup 2. Indeed, the mRNA level of MYCN in subgroup 2 was significantly lower than cases within the MYCN-AMP group (Fig. 4a; Supplementary Table 4; P < 0.0001).To evaluate N-MYC activity in neuroblastoma samples, a total of 87 genes upregulated by N-MYC were selected to classify its activity [19].The MYCN score for each sample was calculated using single-sample gene set enrichment analysis (ssGSEA) based on this 87-gene expression signature.MYCN scores in subgroup 2 were significantly higher than those in subgroups 1 and 3, and were comparable to those in the MYCN-AMP group, although slightly lower (Fig. 4b).Moreover, the MYCN score was an independent predictor of mortality including in multivariate analysis with the risk classification (Fig. 4c; HR: 3.3; P < 0.001).
To investigate the potential mechanism that leads to higher MYCN scores in subgroup 2, correlation analysis coupled with protein−protein interactions network construction was performed (Fig. 4d; Supplementary Table 7).Among the candidate genes, AURKA (Aurora kinase A) was identified to interact with MYCN.AURKA, a serine/threonine kinase regulating the process of mitosis [20], was previously demonstrated to regulate N-MYC protein stability [21].AURKA was expressed at significantly higher levels in subgroup 2 when compared to the other 2 subgroups and its levels were even slightly higher than those in the MYCN-AMP group (Fig. 4e).Classifying MYCN non-amplified neuroblastomas into high and low groups, we demonstrated that the AURKA mRNA levels alone could predict the overall survival (Fig. 4f; HR 4.8; P < 0.0001).In addition, the high level of AURKA was an independent predictor (HR 3, P < 0.001) of mortality including in multivariate analysis with the risk classification (Supplementary Fig. 6).
These findings were further investigated by immunohistochemistry (IHC) staining of N-MYC or AURKA in a custom neuroblastoma        tissue microarray, which contains 94 MYCN non-amplified neuroblastomas.Within them, 22 samples were positive for N-MYC (Fig. 5a), and they had worse survival compared to those with N-MYC negative staining (n = 72) (Fig. 5b; P = 0.03; Supplementary Table 8).In parallel, patients with higher levels of AURKA had unfavourable survival outcomes (Fig. 5c, d; P = 0.00014).Moreover, a higher percentage of patients with high AURKA staining was observed in the N-MYC-positive group compared to the N-MYC-negative group (Fig. 5e; 64% vs. 39%; P = 0.041).Taken together, these results suggested that a "MYCN" signature in subgroup 2 is potentially induced by AURKA overexpression in MYCN non-amplified neuroblastomas.Subgroup 3 is accompanied by an 'inflamed' gene signature Considering immune-related pathways were enriched in subgroup 3, the activity of immune cells and pathways were further systematically explored.ssGSEA was performed to calculate enrichment scores of 46 immune gene sets summarised from two previous studies [22,23], and subgroup 3 showed significantly higher activity of immune cells and pathways compared to the other 2 subtypes as well as MYCN-AMP group (Fig. 6a; Supplementary Table 9).Consistently, cytolytic activity (CYT) or MHC-1 (major histocompatibility complex-1) scores were highest in subgroup 3 (Fig. 6b, c).This was also true when using the ESTIMATE algorithm to evaluate the immune scores, stromal scores, and tumour purity scores in neuroblastomas [24], showing the highest immune and stromal scores, and lowest tumour purity scores in subgroup 3 (Fig. 6d; Supplementary Fig. 7a, b).

S t a g e 1 S t a g e 2 S t a g e 3 S t a g e 4 S t a g e 4 S L o w r i s k H i g h r i s k S t a g e 1 S t a g e 2 S t a g e 3 S t a g e 4 S t a g e 4 S L o w r i s k H i g h r i s k
For a comprehensive assessment of immune cell infiltration, we used CIBERSORTx deconvolution [25] to quantify various immune populations based on a single cell RNA sequencing (scRNA-seq) dataset in MYCN non-amplified neuroblastoma [26] (Supplementary Fig. 7c).While similar immune cell types were present in each subtype, the absolute number of several immune cell populations were markedly increased in subgroup 3, including B cells, myeloid cells, T cells and pDC (plasmacytoid dendritic cells) (Fig. 6e).Finally, to investigate whether subgroup 3 would benefit more from immunotherapy than the other subgroups, we compared the expression matrix of 3 subgroups with published melanoma datasets including response information after treating with immunotherapies [27,28].The SubMap analysis highlighted that patients within subgroup 3 are predicted to respond to anti-PD1 immunotherapy (Fig. 6f; Supplementary Fig. 7d).In addition, Su et al. observed that anlotinib treatment in neuroblastoma mice reprogrammed the immunosuppressive tumour microenvironment (TME) into an immune-stimulatory TME, leading to an extension in the duration of vascular normalisation, and dynamic changes in the expression levels of PD-1 and PD-L1.In addition, it is noteworthy that the combination of anlotinib with PD-1 checkpoint inhibitors counteracted the immune suppression induced by PD-L1 upregulation after monotherapy, ultimately inducing the regression of neuroblastoma [29].Therefore, we reanalysed the RNA-seq data of neuroblastoma syngeneic mouse models treated with vehicle/anlotinib for 9 days.Then, we compared the molecular features of each condition to our subgroups.Interestingly, SubMap analysis revealed that subgroup 3 exhibited a significant similarity in expression profile to mouse models after anlotinib treatment (Fig. 6g; P = 0.032).
Taken together, these results demonstrated that subgroup 3 is accompanied by an 'inflamed' gene signature, and is more likely to benefit from anti-PD1 therapies.

Identification of independent predictors to subgroup patients within MYCN non-amplified neuroblastomas
To identify independent predictors for subgrouping, we applied a multi-cohort analysis pipeline via multiclassPairs [30] (see Supplementary Methods).In total, a random forest model, trained using 928 rules derived from a set of 432 genes (Supplementary Table 10) displayed the ability to predict different subgroups accurately in both train and test sets with an F1 score > 0.74 (Supplementary Table 11).The prediction model and example files can be downloaded from https://zenodo.org/records/10258748.Furthermore, the random forest model successfully stratified patients with MYCN non-amplified neuroblastoma into distinct subgroups 1, 2, and 3 with significant differences in survival across five independent validation sets (GSE49711 [31], TARGET Microarray [32], TARGET RNA-seq, Westermann ALK cohort [33] and Stefan Hüttelmaier cohort [34] respectively) (Fig. 7a, b; Supplementary Table 11).These independent predictors worked consistently between microarray and RNA-seq within GSE47792 (Fig. 7c).

Evaluation of different patient stratification strategies
Finally, we evaluated our subgrouping method (named Hu_Subgroups) together with other reports.van Groningen and colleagues reported that neuroblastoma is composed of 2 superenhancer-associated differentiation states: an 'ADRN' subgroup showing up-regulated genes involved in adrenergic differentiation and an 'MES' subgroup with higher expressions of mesenchymal markers [35].To quantify these characteristics, we calculated the 'ADRN' or 'MES' scores of our subgroups based on a 369-gene 'ADRN' signature or a 485-gene MES signature, respectively.We observed that subgroup 3 showed the highest 'MES' scores and the lowest 'ADRN' scores, consistent with our above findings; while subgroups 1 and 2 had the highest 'ADRN' scores with the lowest 'MES' scores in subgroup 2 (Supplementary Fig. 8a, b).
We also compared Hu_Subgroups with the Valentijn classification [19].All subgroup 1 samples (n = 33) and two-thirds of subgroup 3 (n = 8) belonged to Valentijn's NEG group, while 13 out of 23 subgroup 2 samples were part of Valentijn's POS group (Fig. 8a; Supplementary Table 12).In addition, the multivariate analysis indicated that our subgroup 3 could be an independent variable after being adjusted by Valentijn's classifier (Fig. 8b).
Since 2006, Oberthuer and colleagues have been dedicated to constructing a molecular classification system capable of accurately categorising patients into favourable and unfavourable groups, continually iterating over the following decade [36][37][38][39].The most recent molecular predictors NB-th24 and NB-th44 were introduced in 2017 [40].A comparative analysis between our model and their two models reveals a strong consistency in the favourable and unfavourable outcomes of the respective groupings (Fig. 8c).Specifically, 218 out of 230 subgroup 1 samples and 77 out of 124 subgroup 3 samples were labelled as the favourable group based on SVM_th24.Conversely, more than half of the subgroup 2 samples were categorised as unfavourable (Supplementary Table 12).Similar results were identified in the SVM_th44 comparison (Supplementary Fig. 8c).Additionally, multivariate analysis to determine subgroup 2 could serve as an independent variable after adjusting for Oberthuer's classifier (Fig. 8d; Supplementary Fig. 8d).

Value
Together with other reports, our findings emphasised the extent of inner heterogeneity within the MYCN non-amplified population and the importance of patient stratification.

DISCUSSION
Neuroblastoma remains a challenge in the era of personalised therapy, largely due to inter-and intra-tumoral heterogeneity.Gene amplification in MYCN is the first genetic marker that indicates a highly invasive, advanced neuroblastoma, which has been observed in about 20% of primary and about 40% of highrisk neuroblastoma cases [44].Despite the extensive study of MYCN-amplified neuroblastomas, there is a significant unmet clinical need in MYCN non-amplified neuroblastomas.
In this study, using tumour expression data and ConsensusClus-terPlus, we demonstrate that MYCN non-amplified neuroblastomas are heterogeneous and can be further classified into 3 subgroups based on their transcriptional profiling.Within them, subgroup 2 has the worst prognosis and this group exhibits a 'MYCN' signature that is potentially induced by the overexpression of AURKA.AURKA interacts with both N-MYC and SCF (Fbxw7) ubiquitin ligase, which ubiquitinates N-MYC for degradation.Consequently, overexpression of AURKA counteracts the degradation of N-MYC, leading to the growth of neuroblastoma cells [21,45].Subgroup 3 is accompanied by EMT and an 'inflamed' phenotype, with high expression of genes related to IL2_STAT5 signalling, IL6 JAK STAT3 signalling, interferon-α activation, interferon-γ activation, and inflammatory response, consistent with the association between EMT and immune-related gene expression [46,47].The findings were further confirmed by using CIBERSORTx deconvolution [25] to quantify various immune populations based upon a MYCN non-amplified neuroblastoma scRNA-seq data [26], showing increased percentages of fibroblasts, B cells, myeloid cells, T cells, and pDC (plasmacytoid dendritic cells).
The clinical implications of this subtype classification are significant, as each subtype demonstrates a unique prognosis and vulnerability to investigational therapies.For example, patients in subgroup 1 show the most favourable prognosis with a long-term survival rate above 85%, despite some of them being clinically classified as INSS stage IV or high risk.It might suggest that we should take a more careful and precise evaluation of some patients in reality after the consideration of all clinical information such as age, stage, risk status, or our stratification rather than making a decision based on a single parameter.With regard to therapy stratification, evidence showing significantly high MHC-I and CYT scores in subgroup 3 suggests that patients within this group may benefit from immunotherapy.Our analysis suggests that subgroup 3 is predicted to respond to anti-PD1 immunotherapy.The application of immunotherapy in neuroblastoma has started with treatments such as GD2 monoclonal antibody (dinutuximab) and Chimeric antigen receptor T cells (CAR-T) therapy [48,49].Further studies, including in vitro, in vivo, and clinical validations, are required to investigate if patients within subgroup 3 can benefit from immunotherapy.
In addition, our study suggests that patients within subgroup 2 may benefit from AURKA inhibitors that can disrupt the interaction between AURKA and N-MYC.Indeed, AURKA inhibitors, MLN8054 and MLN8237 (Alisertib), are able to disrupt this interaction, leading to N-MYC degradation and subsequently cell death and differentiation in neuroblastoma cells [45,50].MLN8237 (Alisertib) is currently under phase 2 clinical evaluation in neuroblastoma (NCT01154816).
With the establishment of independent predictors, MYCN nonamplified neuroblastomas were easily classified into one of the 3 subtypes, permitting a realistic scenario in which prospective subtyping is performed in a cohort, wherein patients are assigned to different therapeutics (e.g., subgroup 3 to immunotherapy, subgroup 2 to AURKA inhibitors) based on their subtype.If any one of these predictions demonstrated significant benefit, it would represent the first standard-of-care molecular biomarker selection for MYCN non-amplified neuroblastomas and a foundational step toward personalised therapy for this devastating disease.

Subtype identification
The study design is provided in Fig. 1a with a summary of datasets in the Supplementary Table S1.A detailed description of the approach and further characterisation of the subtypes by PCA, ConsensusClusterPlus, single-sample Gene Set Enrichment Analysis (ssGSEA), weighted gene coexpression network analysis (WGCNA), and CIBERSORTx analysis is provided in the Supplementary Methods.

Analysis of hazard ratio and overall survival
The univariate and multivariate Cox proportional hazards model assessed the hazard ratio of each parameter through the survminer (v0.4.9).We performed a log-rank test to compare Kaplan-Meier survival curves between each subgroup by survival (v3.2-10).Prediction error curves of each prognostic model were generated from pec (v2019.11.03) [51].

Analysis of clinically actionable genes and drug response
To investigate subgroup-specific druggable targets, we performed an integrative analysis to assess the associations between molecular features and the response to anticancer drugs in MYCN non-amplified neuroblastomas.A detailed description of the approach is provided in Supplementary Methods.

Identification of independent predictors
To identify independent predictors for subgrouping, we applied a multicohort analysis pipeline via MetaIntegrator [30] and validated with the machine learning classifier, support vector machine (SVM) (see details in Supplementary Methods).

Tissue microarray (TMA) preparation and immunohistochemistry (IHC)
Separate a small part of the tissue specimen and shape it in a customised mold for chip production and fix it overnight in 4% paraformaldehyde (PFA).Tissue blocks were embedded in paraffin in a prepared array.Then the sample was sliced (5 μm) and adhered to a poly-L-lysine coated glass slide for immunohistochemical staining, which was performed as Fig. 8 A systematic comparison of the subgroup classifier with previously published gene expression classifiers.a Prediction differences in GSE16476 using the subgrouping method from this report (named Hu) or Valentijn and colleagues (Valentijn).b Multivariate analysis of subgroup classification with Valentijn classification in MYCN non-amplified neuroblastomas.HR (hazard ratio), 95% CI (confidence interval), patient number (n), and P values are shown.c Prediction differences in E-MTAB-1781 using the subgrouping method from this report (named Hu) or Oberthuer and colleagues (Oberthuer's svm_th24).d Multivariate analysis of subgroup classification with Oberthuer's svm_th24 classification in MYCN non-amplified neuroblastomas.HR (hazard ratio), 95% CI (confidence interval), patient number (n), and P values are shown.e Prediction differences in GSE49711 using the subgrouping method from this report (named Hu) or Westermann and colleagues (Westermann).Kaplan−Meier plots showing the overall survival in Westermann_MNA-HR (f) or Westermann_MNA-LR (g) patients using the subgrouping method from this report.Numbers below are n (%).P values are indicated.h Prediction differences in GSE49711 using subgrouping method from this report (named Hu) or George and colleagues (George).Kaplan−Meier plots showing the overall survival in George_Immunogenic (i) or George_Neuronal (j) patients using the subgrouping method from this report.The numbers below are n (%).P values are indicated.k Prediction differences in GSE85047 using the subgrouping method from this report (named Hu) or Califano and colleagues (Califano).l Kaplan−Meier plots showing the overall survival in Califano_11q-LOH & MYCNA patients using the subgrouping method from this report.The numbers below are n (%).P values are indicated.
previously described [52,53], using specific antibody against N-MYC (1:600 dilutions; Cell Signalling Technology 51705) and Aurora kinase A (1:200 dilutions; Abcam ab52973).Blindly, with no knowledge of the clinicopathological characteristics of the tumour, the immunoreactivity in tissue sections was observed under three microscopes at random and then evaluated by 3 pathologists.Differences in scoring were discussed until a consensus was reached.The tissue sections were then scored under an optical microscope according to the degree of staining (0−3 points were negative staining, light yellow, light brown, dark brown) and the positive range (1−4 points were 0−25%, 26−50%, 51−75%, 76−100%).Finally, samples were divided into a high-expression group and a low-expression group based on the median of the final staining score.All procedures adhered to the ethical standards set by the Clinical Committee of Xinhua Hospital, Shanghai Jiao Tong University School of Medicine (Approval No: XHEC-D-2016-037).

Fig. 1
Fig. 1 Characterisation of molecular subtypes in the MYCN non-amplified neuroblastomas.a Workflow showing the study design (details provided in the Supplementary Methods).b Consensus clustering of top 50% variable genes of train cohort.c Principal component analysis (PCA) showing neuroblastoma patients with subgroup annotations.

Fig. 2
Fig. 2 Clinical characterisation of subtypes within MYCN non-amplified neuroblastomas identifies key distinguishing features.Graphs showing the frequency (%) of each molecular subtype in different International Neuroblastoma Staging System (INSS) stages or risk status in either train (a) or test (b) cohort.P values are indicated.Kaplan−Meier plots showing the overall survival in each molecular subtype or MYCNamplification (MYCN-AMP) in either train (c) or test (d) cohort.The numbers below are n (%).P values are indicated.e Multivariate analysis of subgroup classification with risk status in MYCN non-amplified neuroblastomas.HR (hazard ratio), 95% CI (confidence interval), patient number (n), and P values are shown.f Prediction error curves (indicating a mean squared error in predicting survival status) are calculated for the subgroup (red) and INSS stage (green).

Fig. 3
Fig. 3 Defining molecular features of 3 subtypes in MYCN non-amplified neuroblastomas.a Heatmap showing differential expression of selected genes.Red indicates up-regulation and blue for down-regulation.Colour bars show subgroup information.b Gene set enrichment analysis (GSEA) in 3 subtypes.*FDR (false discovery rate) <0.25; **FDR < 0.05; ***FDR < 0.01.c Weighted gene co-expression network analysis (WGCNA) showing 3 molecular modules.Nodes are colour-coded according to the WGCNA modules.Representative enriched pathway terms are indicated.d Overlay of the median-cantered log 2 fold change values per subgroup on the network.

Fig. 4
Fig. 4 Subgroup 2 shows a "MYCN" signature, potentially induced by Aurora Kinase A overexpression.Violin plots showing MYCN mRNA levels (a) or MYCN scores (b) in neuroblastomas.P values are indicated.c Multivariate analysis of MYCN score and risk status in MYCN nonamplified neuroblastomas.HR (hazard ratio), 95% CI (confidence interval), patient number (n), and P values are shown.d Protein-protein interaction (PPI) network showing an interaction between AURKA and MYCN.e Violin plot showing AURKA mRNA levels in neuroblastomas.P values are indicated.f Kaplan−Meier plot showing the overall survival in samples with low vs. high AURKA expression.The numbers below are n (%).HR (hazard ratio), 95% CI (confidence interval), patient number (n), and P values are shown.

Fig. 5 N
Fig. 5 N-MYC expression correlates with Aurora kinase A status in MYCN non-amplified neuroblastomas and is indicative of patient survival.a Representative N-MYC staining pattern (negative or positive N-MYC) in MYCN non-amplified neuroblastoma tissue microarray cores.Scale bar: 1 mm (the left column) and 50 μm (the right column).b Kaplan−Meier plot showing the overall survival in samples with negative vs. positive N-MYC expression.The numbers below are n (%).HR (hazard ratio), 95% CI (confidence interval), patient number (n), and P values are shown.c Adjacent tumour sections from representative cases showing N-MYC and Aurora Kinase A expression in MYCN non-amplified neuroblastoma.Scale bars: 50 μm.d Kaplan−Meier plot showing the overall survival in samples with low vs. high Aurora kinase A expression.The numbers below are n (%).HR (hazard ratio), 95% CI (confidence interval), patient number (n), and P values are shown.e Graph showing percentage (%) and numbers of samples with low or high Aurora kinase A in the negative or positive N-MYC group.P = 0.041.

Fig. 6
Fig. 6 Subgroup 3 is accompanied by an "inflamed" gene signature.a Heatmap showing neuroblastoma-associated immune pathways and immune cell signatures in subgroups and MYCN-AMP.Graphs showing the cumulative distribution of CYT (b) or MHC-1 (c) scores in different subgroups and MYCN-AMP.d Violin plots showing immune scores in different subgroups and MYCN-AMP in train, test, or train plus test cohort.e Graph showing cell compositions of each subgroup using CIBERSORTx analysis.f Graph showing differential putative immunotherapeutic response in different subgroups.Bonferroni adjusted P values indicated.g Subclass association (SA) matrix for the comparison between different subgroups and vehicle/anlotinib treated mouse.Bonferroni adjusted P values indicated.