Proteomic analysis reveals chromatin remodeling as a potential therapeutical target in neuroblastoma

Abstract

Background

Neuroblastoma (NB) is the most common solid tumor in children, characterized by high recurrence rates, drug resistance, and significant mortality.

Methods

In this study, we analyzed the proteomic profiles of NB tissue samples alongside other pathological categories, including ganglioneuroma (GN) and ganglioneuroblastoma (GNB). Using weighted gene co-expression network analysis (WGCNA), the core prognostic gene models associated with histopathology of NB were identified. Furthermore, by mapping our core prognostic gene models onto drug-perturbed transcriptome profiles from the L1000FWD and CMap databases, repurposing drug candidates were screened and validated for NB.

Results

Our proteomic analysis reveals that pathways associated with the cell cycle and DNA replication are significantly upregulated in NB, while oxidative phosphorylation, pyruvate metabolism, and the TCA cycle are notably downregulated compared to GNB and GN. By applying WGCNA, we identified a core prognostic gene model strongly associated with the unfavorable subtype and high MKI of NB and primarily related to chromatin binding and mRNA metabolic process. Protein–protein interaction network analysis identified 15 hub genes in this core prognostic module: SMARCA4, SMARCA5, SMARCC2, SMARCC1, PBRM1, BRD3, ARID1A, BRD2, ARID1B, KDM1A, TP53BP1, ALYREF, CBX1, SF3B1, and ADNP, which mainly related to chromatin remodeling. Notably, SMARCA4 and ALYREF are also high-risk genes of mortality and validated as potential prognostic biomarkers for NB. Through repurposing drugs screening, mocetinostat and clofarabine were validated as effective treatments in two NB cell lines.

Conclusion

Mocetinostat and clofarabine offer valuable insights for the development of novel targeted therapies in neuroblastoma.

Introduction

Neuroblastoma (NB), a cancer that arises from the developing sympathetic nervous system, is the most common extracranial solid tumor in children, accounting for approximately 15% of childhood cancer-related deaths [1, 2]. The genetic, morphological, and clinical heterogeneity of NB limits the effectiveness of current treatment modalities [3]. Prior to 2009, evidence-based treatment for high-risk NB involved a multimodal approach, including surgery, local radiotherapy, and combination chemotherapy regimens, often supplemented by consolidation protocols that utilized autologous stem cell transplantation [4]. The introduction of immunotherapy using anti-disialoganglioside (GD2) monoclonal antibodies (mAbs) in combination with chemotherapy has significantly improved the 5-year survival rate for patients with metastatic NB, increasing it from less than 20% to over 50% [2].

The recently revised Children’s Oncology Group (COG) classifier incorporates factors such as age, INRG stage, and tumor histologic and genetic features, including MYCN amplification status and chromosomal aberrations [5]. Although many potential biomarkers identified through next-generation sequencing techniques have been studied, very few have demonstrated independent prognostic significance in multivariable analyses compared to established classifiers like MYCN and stage [6]. To date, large-scale analyses, including those from the Cancer Cell Line Encyclopedia (CCLE), have primarily focused on genetic information, while a thorough exploration of the proteomic profile remains largely unexplored [7]. The protein expression profile of clinical tissues from NB patients is still unknown. Recent studies utilizing advanced high-throughput “omics” technologies have uncovered numerous genetic and genomic alterations, as well as dysregulated pathways that drive the onset, progression, and treatment resistance of NB [8, 9]. Investigating the molecular changes within these pathways may help identify new therapeutic targets and approaches for NB.

Histological features are particularly impactful in predicting outcomes for patients with locoregional disease and for children with metastatic disease [6]. In NB, the degree of differentiation and the mitosis-karyorrhexis index have independent prognostic significance and, alongside age, inform the favorable or unfavorable International Neuroblastoma Pathology Classification (INPC) histology [10, 11]. In this study, we aim to identify key gene models and hub genes associated with pathological histology in NB, based on a proteomic analysis of tumor tissues from NB, GNB, GN, and AG. Furthermore, we plan to conduct computational screening for drug repositioning for NB and validate these potential drugs in NB cell models.

Materials and methods

Study design and clinical characteristics of the patient groups

All subjects were enrolled at the Department of Pediatric Surgery in Hunan Childen’s Hospital and clinically diagnosed as ganglioneuroma (GN), ganglioneuroblastoma (GNB) and NB patients. All patients underwent laboratory and imaging examinations at the hospital. Non-renal NB were excluded. Clinical samples consist of 9 NB samples, 6 GN samples, 4 GNB samples and 5 adrenal gland (AG) samples. Proteins of the tumor tissues were enriched and quantitatively analyzed by data-independent acquisition (DIA) LC–MS/MS. This study was approved by the Medical Ethics Committee of Hunan Children’s Hospital (No: HCHLL-2021-110). The clinical information and protein expression profiles of the patients were collected for prognosis analysis of NB study. The clinical characteristics and histopathological phenotypes are described in the Table S1 of Additional file 1. In this study there is another NB cohort for verification experiment via western blot, including 6 NB, 2 GNB and 5 GN. Their clinical characteristics and histopathological classification are described in the Table S2 of Additional file 1.

Sample preparation and LC–MS/MS analysis

The sample was grinded with liquid nitrogen into cell powder and then transferred to four volumes of lysis buffer (8 M urea, 1% protease inhibitor cocktail), followed by sonication three times on ice. The remaining debris was removed by centrifugation at 12,000g at 4 °C for 10 min. For digestion, the protein solution was reduced with 5 mM dithiothreitol for 30 min at 56 °C and alkylated with 11 mM iodoacetamide for 15 min at room temperature in darkness. The protein sample was then diluted by adding 100 mM TEAB to urea concentration less than 2 M. Protein samples were digested by trypsin in a trypsin-to-protein mass ratio of 1:50 for the first digestion overnight and 1:100 trypsin-to-protein mass ratio for a second 4 h-digestion. The digested peptides were desalted by Strata X C18 SPE column (Phenomenex) and vacuum-dried.

Peptides are fractionated using high pH reverse-phase HPLC. The chromatography column used is an Agilent 300 Extend C18 (5 μm particle size, 4.6 mm inner diameter, 250 mm length). The gradient for peptide fractionation is set at 8%–32% acetonitrile, pH 9, over 60 min. The peptides are then combined into 12 fractions. Following the instructions in the iRT reagent manual, an appropriate amount of iRT reagent is added to each combined fraction, and the samples are vacuum freeze-dried for subsequent processing.

The peptides are dissolved in mobile phase A of the liquid chromatography and then separated using a NanoElute ultra-high-performance liquid chromatography system. Mobile phase A consists of 0.1% formic acid and 2% acetonitrile in water, and mobile phase B is acetonitrile with 0.1% formic acid. The LC gradient is set to 90 min with a flow rate maintained at 450 nL/min. After separation by the UHPLC system, peptides are ionized in the capillary ion source and then injected into the tims-TOF Pro mass spectrometer for data acquisition.

The ion source voltage is set to 1.7 kV, and both the precursor ions and their secondary fragments are detected and analyzed using TOF. Data acquisition is performed in data-independent parallel accumulation and serial fragmentation (dia-PASEF) mode. The first-stage mass spectrum scan range is set from 100–1700 m/z, with 16 PASEF acquisitions following each primary mass spectrum. The second-stage mass spectrum scan range is set from 400–1200 m/z, with a window of 25 m/z per scan.

DDA data is searched using MSFragger (v 2.3). The database used is SwissProt_Human (20,380 sequences), with a decoy database added to calculate the false discovery rate (FDR) due to random matching. The enzyme setting is Trypsin/P, allowing up to 2 missed cleavage sites. The precursor ion mass tolerance is set at 20 ppm, and the fragment ion mass tolerance is 0.02 Da. Cysteine alkylation is set as a fixed modification, with variable modifications set for methionine oxidation and protein N-terminal acetylation. The FDR for protein and PSM identification is set at 1%.

DIA data is processed using Skyline (v 20.1.0) software, importing the corresponding spectral library and adding the iRT parameters under Prediction. Transition precursor ion charges are set to 2, 3, 4, and 5, and fragment ion charges to 1 and 2. The top six ions with the highest intensity in the spectral library are extracted for peptide quantification. After generating the decoy library, DIA data is imported and FDR filtering is applied using the mProphet algorithm. The MSstats R package is used to obtain relative quantification results for proteins.

Proteomic data analysis

Enrichment of biological function of DEGs was performed via online platform metascape (https://metascape.org/gp/index.html). In the analysis the functional set (GO molecular functions) and pathways (GO Biological Processes, Reactome Gene Set, KEGG Pathway, WikiPathways, Canonical Pathways, Hallmark Gene Sets and BioCarta Gene Sets) were involved and applied for enrichment of biological processes.

WGCNA for histopathology-related prognostic gene modules of NB

To find out core clinicopathological prognostic gene modules, weighted gene co-expression network analysis (WGCNA) is used to cluster the genes with similar co-expression mode into a gene module. Then the correlation between gene modules and clinicopathological phenotypes was analyzed for the core prognostic gene modules of NB. From a pathological perspective, there is a Shimada classification method that divides NB into two types: Favorable Histology Group (FH) and Unfavorable Histology Group (UH) [12]. The correlation between gene expression modules and clinical prognosis in the clinical datasets of GN, GNB and NB was analyzed using Spearman correlation analysis. The criteria for selecting significant DEGs for WGCNA were described as following. Firstly, the DEGs(|FC|≥ 2, P < 0.05) of each group NB/AG, GNB/AG, GN/AG, NB + GNB/AG, NB + GNB + GN/AG were integrated together. Then the DEGs must be identified and quantified in at least 21 samples among total 24 proteomic samples. Finally, it leads to 3313 DEGs for further WGCNA analysis. Co-expression gene modules are constructed from the DEGs and WGCNA was performed using Sangerbox (http://sangerbox.com/home.html) online platform. Heatmaps depicting the correlation of each module with clinical prognosis were generated. We used the "plot()" function to generate scatter plots between gene significance (GS) and module membership (MM) within each module in order to understand the importance of highly connected genes within the module. The differentially expressed proteins in NB were subjected to protein–protein interaction (PPI) network analysis using the STRING database (version 11.5, http://string-db.org/) and Cytoscape software (version 3.8.2, https://cytoscape.org/) to identify hub genes. Basic settings of STRING were described as following. Active interaction sources include textmining, experiments, databases, co-expression, neighborhood, gene fusion and co-occurrence. Minimum required interaction score was set as medium confidence (0.400). Max number of interactors to show in the 1st shell was set as “none/query proteins only” and 2nd shell was “none”.

The plugin, cytoHubba (version 0.1), was used to identify hub genes, and the top 15 hubba nodes’ scores were calculated by using the Maximum Clique Centrality (MCC) method and ranked. GO pathway enrichment analysis was conducted on hub genes. P < 0.05 in pathway enrichment analysis was considered statistically significant.

Transcriptome proximity analysis in L1000FWD and CMap

The reversed expression profile of core prognostic modules was used to query L1000FWD online database (http://amp.pharm.mssm.edu/L1000FWD). The similarity score of predicted drugs were given to quantify the similarity between the reversed protein expression of the core prognostic regulated genes and their drug-perturbed gene expression profiles in L1000FWD. Based on the similarity score, p-value, comprehensive score and other factors, the relevant candidate drug molecules were selected, so fa their similarity score absolute value greater than 0.1.

The CMap (https://clue.io) is another online database for large-scale drug query. As a complementary query in L1000FWD, the core prognostic gene modules were used to query CMap. The similarity of the query to each CMap signature was computed and yielded a rank-ordered list of the signatures. A tau ≥ 90 was considered as convinced strong score.

Bioassays

CCK8 assay

The human neuroblastoma cell line SH-N-AS was cultured in EMEM (with NEAA) containing 10% fetal bovine serum. Cells were seeded into 96-well plates (1000 per well), treated for 24 h with different concentration of PI-828, mocetinostat, clofarabine, ethacrynic acid and mafenide (10^–8–10^–5 mol/L, all purchased from MCE). The cells were incubated at 37℃ for 0.5 h with CCK-8 Kit (Abiowell, China). The experiment is performed according to the manual of the kit.

RNA extraction and quantitative real-time PCR (qRT-PCR)

Total RNA from clinical tumor tissues were extracted with trizol. Reverse transcription was performed using the mRNA reverse transcription kit (CW2569, CWBIO, China). qRT-PCR was carried out on the CFX Connect Real-Time PCR system (1855201, Bio-Rad Laboratories, Hercules, CA, USA) using UltraSYBR Mixture (CW2601, CWBIO, China). The expression levels of all target genes were normalized to β-actin. The quantification of mRNA was carried out using the 2^−ΔΔCt method. The primers used are shown in Table 1.

Table 1 Primer sequences for target genes

Western blot (WB)

The proteins from the samples were separated by SDS-PAGE and then transferred to a polyvinylidene difluoride (PVDF) membrane. The blocked membrane was incubated with the indicated primary antibodies. The antibodies used were ALYREF (1:1000, 16690-1-AP, Proteintech, USA), SMARCA4/BRG1 (1:1000, ab110641, Abcam) and β-actin (1:5000, AWA80002, Abiowell, China). Next, the membrane was incubated with HRP-conjugated goat anti-rabbit IgG (SA00001-2, Proteintech, USA). The immunoreactive bands were visualized using a SuperECL Plus ultrasensitive luminescent liquid (awb0005, Abiowell, China).

Data analysis

The survival analysis and immune filtration analysis of the hub genes were performed by using Sangerbox online platform. The data was presented as mean ± standard deviation. All data were collected from at least three independent experiments. Student's t-test or Wilcoxon test was used to analyze the differences between the two groups. One-way analysis of variance (ANOVA) was used for comparisons among multiple groups. Survival curves were depicted using Kaplan–Meier plots, and the log-rank test was used for comparisons. The software used were R and GraphPad Prism (v8.0.1, GraphPad Software, USA). P < 0.05 was considered statistically significant.

Results

Identification of differentially regulated proteins based on proteomic analysis of NB, GN, GNB, and AG tumor tissues

To identify neuroblastoma (NB)-specific biological pathways distinct from other pathological categories, we compared tumor tissues from NB, GNB, and GN with those from AG using proteomic analysis. Principal Component Analysis (PCA) of protein levels in patient samples revealed clear distinctions between the groups and similar protein profiles within each group (Fig. 1A). A total of 8525 proteins were identified and quantified (Fig. 1B). Differentially expressed genes (DEGs) were defined as |Fold Change|≥ 2 and P < 0.05, leading to the identification of 872 upregulated and 1168 downregulated proteins in NB/AG, 2945 upregulated and 124 downregulated proteins in GN/AG, and 1258 upregulated and 688 downregulated proteins in GNB/AG (Fig. 1C). Notably, cell cycle and DNA replication pathways were highly upregulated in NB/AG, while oxidative phosphorylation, pyruvate metabolism, and the TCA cycle were significantly downregulated compared to GNB/AG and GN/AG (Fig. 1D).

Fig. 1

Proteomic analysis of NB study. A PCA analysis of proteome samples in NB study; B the statistical data of identified peptides and proteins in NB proteome study; C the number of regulated proteins; D enrichment of KEGG pathway of proteomic DEGs

Identification of clinical prognosis-related gene modules based on histopathology of NB

NB is classified into FH, and UH phenotypes based on its pathological histology, which are closely associated with NB outcomes. Compared to survival analysis, pathological histology evaluation is more precise and can reflect the severity of NB in real-time. To construct a more precise prognostic gene model for NB, we compared the proteome profiles of NB, GN, GNB, NB + GNB, NB + GNB + GN with adrenal glands (AG) respectively. A total of 3,313 differentially regulated proteins were identified. Using WGCNA, these differentially expressed proteins were clustered into 6 gene modules (Fig. 2A, B, E). Correlation analysis with histopathology classification revealed that high expression of the yellow module was positively associated with the UH subtype and high MKI, while the black module was negatively associated with the UH subtype and high MKI, and the red module was negatively associated with the high MKI phenotype (Fig. 2D, P < 0.05). Therefore, the yellow module is a key upregulated gene expression module associated with poor NB outcomes, enriched in functional pathways such as chromatin binding and mRNA metabolic processes (Fig. 2F). PPI network analysis identified 15 hub genes in the yellow module: SMARCA4, SMARCA5, SMARCC2, SMARCC1, PBRM1, BRD3, ARID1A, BRD2, ARID1B, KDM1A, TP53BP1, ALYREF, CBX1, SF3B1, and ADNP, which is mainly associated with chromatin remodeling (Fig. 2C). Verification in the TARGET NB database, which includes transcriptomic data and clinical survival information from 159 NB patients, confirmed that high expression of ALYREF and SMARCA4 was associated with high mortality risk, consistent with our proteomic results (Fig. 3 A-B). Additionally, experimental validation in another NB cohort showed that RNA expression of TP53BP1, ALYREF, SMARCC2, and SMARCA4 was also upregulated in NB tumor tissues compared to GN tissues (Fig. 3C) and protein expression of ALYREF and SMARCA4 were much higher in NB tumor tissues compared to GN tissues (Fig. 3D, E, Figs. S1, S2 in the Additional file 2).. Furthermore, it is also demonstrated that the protein expression level of ALYREF and SMARCA4 in tumor samples are correlated with their histopathological phenotypes (Table S3 in the Additional file 1), e.g. INPC (International Neuroblastoma Pathology Classification) and INRG (International Neuroblastoma Risk Group) stage. Especially, SMARCA4 is strongly associated with INPC phenotype (Spearman’s ρ correlation coefficient is 0.803, p = 0.001). Therefore, ALYREF and SMARCA4 are potential prognostic biomarkers for NB patients with poor outcomes. Furthermore, high expression of these genes is associated with low immune infiltration (Fig. 3F).

Fig. 2

Construction of a clinicopathology-related prognostic gene modules. A The gene number of each co-expression modules; B correlation membership between gene co-expression modules; C PPI analysis of genes in yellow module and identification of the top 15 hub genes. D Correlation between gene expression modules and clinical prognosis was analyzed through weighted gene Co-expression Network analysis (WGCNA). The x-axis of the scatter plot represents the correlation between module Eigen genes and gene expression profiles, while the y-axis represents the absolute correlation between genes and prognosis (UH). E Correlation between gene co-expression modules and clinical samples; F enrichment of biological process in yellow and black modules

Fig. 3

Validation of ALYREF and SMARCA4 as prognostic genes in NB. A, B Survival analysis of ALYREF and SMARCA4 based on the TARGET-NB dataset; C detection of RNA expression of TP53BP1, ALYREF, SMARCC2, and SMARCA4 in clinical GN and NB samples by qRT-PCR. D, E Detection of ALYREF and SMARCA4 protein expression levels in the clinical GN and NB tumor samples by western blot; F immune filtration analysis of ALYREF and SMARCA4 in NB

Computational screening of drug repositioning candidates for NB

To discover effective and safe therapeutic drugs against NB, we conducted a drug repurposing search based on the proximity of our core prognostic gene module profile to data in L1000FWD and Cmap (Figs. 2F, 4 and Fig. S4). We queried the opposed expression profile of the core prognostic gene set in these databases. In L1000FWD, we identified 50 drugs with a similarity score > 0.1, targeting dopamine receptors, cyclooxygenase, topoisomerase, calcium channels, and histamine receptors. In CMap, we identified 310 potential anti-NB drugs with tau > 90, targeting HDAC, mTOR, mitochondrial oxidative phosphorylation, PI3K, and other pathways. There were 31 overlapping potential drugs in both databases. From these, we selected the top 5 drugs for further analysis: PI-828, mocetinostat, clofarabine, ethacrynic acid, and mafenide. These drugs target HDAC, PI3K, ribonucleotide reductase V, GST, and bacterial pathways.

Fig. 4

Detailed chart for therapeutic targets discovery and drug discovery

Experimental validation of repurposing drugs in neuroblastoma cells

To verify the anti-tumor activities of the predicted drugs, we tested the top 5 overlapping drugs from L1000FWD and CMap in neuroblastoma cell lines SH-N-AS and SH-SY-5Y.. Among them, only mocetinostat and clofarabine significantly increased cell death of both NB cell lines (Fig. 5A–D). Meanwhile, mocetinostat and clofarabine have been demonstrated to induce apoptosis of NB cell line (SH-SY-5Y) significantly, which indicate that these two drugs may induce cell death of NB cells by the way of apoptosis (Fig. 5E, F). Furthermore, treatment with mocetinostat and clofarabine suppressed the expression of SMARCA4, suggesting that these two drugs may target our core prognostic gene SMARCA4 and affect yellow module expression (Fig. 5G, Fig. S3 in Additional file 2). These two drugs had not been previously reported as anti-NB agents, warranting further validation.

Fig. 5

A–D The cellular viability test under the treatment of mocetinostat andclofarabine in human neuroblastoma cell lines SH-SY5Y and SH-N-AS. E Detection of cell apoptosis in cell line SH-SY-5Y under the treatment of mocetinostat and clofarabine by flow cytometry. F Statistical analysis of the sum of Q2 + Q4 for the cell early apoptosis and late apoptosis events under the treatment of mocetinostat and clofarabine. Unpaired t test was utilized to demonstrate statistical difference between DMSO control group and treatment groups. * P < 0.05, ** P < 0.01. G Detection of SMARCA4 protein expression level under the treatment of mocetinostat and clofarabine in cell line SH-N-AS by WB

Discussion

The treatment of neuroblastoma (NB) remains a significant challenge, with a lack of highly effective prognostic biomarkers and therapeutic targets. To our knowledge, this study is the first to analyze NB tissue proteomes and correlate them with their histopathological phenotypes. Through Weighted Gene Co-expression Network Analysis (WGCNA), we identified a yellow gene module enriched in pathways related to chromatin binding and mRNA metabolism, which correlates strongly with unfavorable histology (UH) and high mitosis-karyorrhexis index (MKI) phenotypes. Among the hub genes within this module, ALYREF and SMARCA4 are associated with a high mortality risk in NB, as shown by survival analysis using the TARGET-NB dataset. RNA and protein expression levels of these genes were validated in NB and ganglioneuroma (GN) tissues, supporting ALYREF and SMARCA4 as potential prognostic biomarkers for NB.

ALYREF, a nuclear protein with known roles in tumor growth and progression, is frequently expressed due to chromosome 17q21-ter gain in NB, where it stabilizes the MYCN protein [13, 14]. SMARCA4 encodes the transcription activator BRG1, a frequently altered ATP-dependent catalytic subunit of SWI/SNF chromatin-remodeling complexes involved in transcriptional regulation across cancer types [15]. SMARCA4 expression is consistently upregulated in advanced NB stages, with high levels correlating with reduced event-free and overall survival [16]. Our findings confirm that these two proteins are linked not only to a high mortality risk but also to UH phenotype in NB.

Given the high heterogeneity of NB, high-risk patients often require intensive chemotherapy, yet only less than half experience a cure. Chemoresistance in NB is a pressing issue, highlighting the urgent need for effective and safe therapeutic agents. Our computational drug screening identified mocetinostat and clofarabine, two compounds previously unreported in NB, as promising therapeutic candidates based on their proximity to our core prognostic gene module profile. Mocetinostat, an isoform-selective histone deacetylase (HDAC) inhibitor, has demonstrated anticancer potential. By inhibiting HDAC, mocetinostat promotes histone acetylation, activating gene expression and impacting cancer cell proliferation, apoptosis, angiogenesis, and inflammation [17, 18]. Currently in clinical trials for various cancers, such as follicular lymphoma, Hodgkin’s lymphoma and acute myelogenous leukemia, mocetinostat’s diverse mechanisms make it a promising cancer therapy candidate [19].

Clofarabine, a nucleotide analog, disrupts DNA replication and repair by inhibiting DNA polymerase and ribonucleotide reductase [20]. It interferes with nucleotide metabolism, blocks DNA synthesis, and triggers apoptosis through mitochondrial pathway disruption [21]. Clofarabine has shown high therapeutic efficacy in pediatric patients with refractory or relapsed acute lymphoblastic leukemia (ALL) [22]. Although our study demonstrates clofarabine’s inhibitory effect on NB cell viability and proliferation, further exploration is required to determine its role in high-risk NB therapy. Notably, both mocetinostat and clofarabine involve chromatin remodeling and epigenetic regulation, suggesting that targeting chromatin remodeling could be an effective strategy for inhibiting NB cell proliferation and survival.

In NB studies, genetic alterations in chromatin remodeling related genes have been documented [9]. Our proteomic analysis of differentially expressed proteins in NB tissues, combined with public database insights, points to chromatin remodeling as critical therapeutic target for NB. Moving forward, we aim to investigate the specific regulatory mechanisms of these identified hub genes in NB, with the objective of identifying more direct therapeutic targets.

Conclusion

This study identified a core prognostic gene module in neuroblastoma (NB), primarily associated with RNA processing and chromatin remodeling. SMARCA4 and ALYREF were validated as key prognostic biomarkers for NB patients. Through drug repurposing, we identified five top candidates for NB treatment, with mocetinostat and clofarabine showing efficacy in two NB cell lines. Further validation of these drugs in a murine NB model is strongly recommended, offering valuable insights for developing novel molecular targeted therapies for NB.