Sequencing of bulk RNA data and corresponding clinical information of patients with NB were obtained from Therapeutically Applicable Research to Generate Effective Treatments (TARGET, n = 159), Gene Expression Omnibus (GSE4971, n = 498), and Array Express (AE) E-MTAB-8248 (AE, n = 223).

Based on the results of single gene analysis in the GSE49710 dataset, we performed WGCNA to identify a set of genes associated with poor prognosis in the pathway. Initially, we computed the median absolute deviation (MAD) for each gene using the gene expression profile and excluded the bottom 50% of the genes with the smallest MAD. We also eliminated the sample and gene outliers to construct a scale-free co-expression network. Next, we merged modules with distances below 0.25 and generated scatter plots of module gene significance (GS) versus module membership (MM) for genes in the modules strongly linked to poor prognosis. We applied strict filters of GS $>$0.1 and MM $>$0.8 to identify core genes and applied a weight threshold of 0.1. We then intersected the cell cycle-related and core genes and defined these genes as CCGs (Cell Cycle Genes).

To further investigate the impact of CCG expression patterns on patients with NB, we conducted consensus clustering using the TARGET and GSE49710 databases. The k-means algorithm was used to identify distinct cell cycle-related expression patterns. The “ConsensusClusterPlus” R package was utilized to construct consistent clusters quantitatively, and we performed 1000 iterations to enhance the stability of these groups. Subsequently, we performed gene set variation analysis (GSVA) using the Kyoto Encyclopaedia of Genes and Genome (KEGG) gene set (c2.cp.kegg. v7.4) to examine differences in biological functions across the different subgroups. Furthermore, we examined the association between genomic expression and clinical variables across clusters.

To develop an accurate prognostic model for predicting survival outcome in NB based on the cell cycle, differential analysis was conducted using the Lima 2.13.0 (Pacific Biosciences, Menlo Park, California, USA) on the two subgroups of the TARGET database. Differential genes were identified based on $|$Log Foldchange$|$ $>$2 and p 0.01. From these genes, a final set was selected using univariate Cox regression analysis with a significance level of p 0.01. The output gene sets obtained from the three ML methods (RF, SVM, and XGboost) were merged, and LASSO dimensionality reduction was applied to improve the accuracy that is defined as (3+x). Based on the gene feature selection results of the above four ML algorithms (RF, SVM, XGboost, and 3+x), we then matched with 20 mainstream ML algorithms (Linear Regression, Ridge Regression, RidgeCV, Linear Lasso, Lasso, Elastic Net, BayesianRidge, Logistic Regression, Stochastic Gradient Descent, SVM, k-Nearest Neighbours, Naive Bayes, Decision Tree, Bagging, RF, Extra Tree, AdaBoost, Gradient Boosting, Voting, and Artificial Neural Network) using TARGET as the training set and GSE49710 and AE as the validation sets to construct 80 algorithm combinations to screen the final genes.

We excluded non-coding RNAs and subjected the final gene set to Multivariate Cox regression analysis. This approach allowed us to develop a prognostic model to calculate risk scores. However, due to the small sample size of the TARGET database, we validated our model using two external databases (AE and GSE49710). Then we used two R packages (TimeROC 0.4 (Paul Blanche, Copenhagen, Denmark) and SurvivalROC 1.0.3.1 (Patrick J. Healing, Seattle, Washington, USA)) to construct receiver operating characteristic (ROC) curves and calculate the AUC.

To facilitate the implementation of risk scores in clinical practice, we integrated these scores with four relevant clinical variables from TARGET: age at diagnosis, International NB staging system (INSS), Children’s Oncology Group (COG) risk classification system, and MYCN gene amplification status. Then we plotted a nomogram to predict the overall survival of patients with NB. We used the “rms” R package to generate line plots displaying the predicted survival probability at 1, 3, and 5 years. Additionally, we conducted a calibration curve and decision curve analysis to assess the performance, clinical usefulness, and external validity of the nomogram. Using clinical information from TARGET, GSE49710, and AE, we compared the value of risk scores to other clinical variables.

We built a web calculator based on the risk model using Shiny Server (“Shiny” package in R version 1.8.0). The risk score was calculated based on the expression values of the six signature genes and compared to the cutoff to assign a query sample to either a high- or low-risk group. Based on all the sample information from the above three public databases, we plotted survival curves on a web page for reference.

The proximity algorithm is based on the target gene database of the extracted drug and the hub gene of NB to perform omics correlation analysis on these two types of genes. In the human genome interaction network, the average shortest path of each drug target and disease-related gene recorded in the database is calculated. If the target gene of the drug is closer to the disease-related gene, it means that the drug plays a certain role in the treatment of the disease. The distance between the drug and the disease-related genes is called proximity, which is converted to a z-score. After 1000 random perturbations, the significant p value of each drug can be obtained to screen out the targeted drugs [18, 19, 20, 21].

Transcriptome sequencing analysis was performed on 12 samples using the following workflow: total RNA extraction, mRNA enrichment, double-stranded cDNA synthesis, end repair, A-addition, adapter ligation, fragment selection, PCR amplification (MiniAmp, Thermo Fisher, Waltham, Massachusetts, USA), library detection, and Illumina sequencing. This process generated a total of 132.86 Gb of clean data, with each sample contributing at least 12 Gb of clean data and a Q30 base percentage of 95% or higher.

SK-N-SH, SK-N-BE2, and SH-SY5Y cells were purchased from Wuhan Pricella Biotechnology Co., Ltd. (Wuhan, China). They were supplemented with 10% (v/v) fetal bovine serum (FBS; Hy Clone, Logan, UT, USA) and stored in a 37 °C humidified environment containing 5% CO₂. The c-Myc (#10828-1-AP, Dilution ratio 1:1000), GAPDH (#10494-1-AP, Dilution ratio 1:1000), cyclin-dependent kinase 2 (CDK2) (#10122-1-AP, Dilution ratio 1:2000), cyclin E1 (CCNE1; #11554-1-AP, Dilution ratio 1:500), and Tubulin (#14555-1-AP, Dilution ratio 1:1000) antibodies were purchased from Protein tech (Wuhan, China). CDK4 (#sc-56277, Dilution ratio 1:1000) and CDK6 (#sc-53638, Dilution ratio 1:1000) antibodies were purchased from Santa Cruz (Dallas, TX, USA). The CCND1 (#A19038, Dilution ratio 1:2000) antibodies were purchased from ABclone (Wuhan, China). The cycle detection kit (#C1052), Edu detection kit (#C0078S), and gentian violet (#C0121) were purchased from Beyotime Biotechnology (Shanghai, China). Dimethyl sulfoxide (DMSO; #D2650) was purchased from Sigma-Aldrich (St. Louis, MO, USA). The apoptosis detection kit (#40302ES60) was purchased from Yeasen (Shanghai, China), and eugenol (#HY-N0337) was purchased from MedChemExpress (Shanghai, China). Eugenol was dissolved in DMSO to generate a 10 mM stock solution. The samples were stored at –80 °C for in vitro studies. Mycoplasma testing using Wuhan Procell’s Mycoplasma Detection Kit yielded negative results for all cell lines (no mycoplasma contamination). This study adhered to the principles of the Helsinki Declaration and was approved by the Ethics Committee of Children’s Hospital Affiliated to Chongqing Medical University (No. 2024312). All cell lines underwent short tandem repeat (STR) analysis for identification shortly before use, with the analysis performed by the analytical testing center of Wuhan Pricella Biotechnology Co., Ltd.

The proliferation rate of NB cells was detected using Cell Counting Kit-8 (CCK-8) assays. To assess the cell viability, we cultured 5000 cells/well in a 96-well plate for 48 h. The viability of SK-N-BE2 cells—which grew as both semi-adherent and semi-suspended cultures—was measured every 12 h using a CCK-8 kit (#SB-CCK8-100 mL; Share-Bio, Shanghai, China) every 12 h.

The cells and tissues were lysed with RIPA and PMSF (Product Code: P0013 and ST507, Beyotime Biotechnology, shanghai, China). The proteins were separated by 8–12% sodium dodecyl sulfate polyacrylamide gel electrophoresis and transferred to PVDF membranes. After being blocked in 5% BSA (Product Code: ST023, Beyotime Biotechnology, shanghai, China), the membranes were incubated with primary antibodies overnight at 4 °C, followed by incubation with horseradish peroxidase-conjugated secondary antibodies for 1 h at room temperature. The VIBER FUSION FX6 EDGE imaging system (Product Code: FX6 EDGE, Paris, France) was used to expose the membranes.

The DNA synthesis ability of NB cells was evaluated using the BeyoClick™ Edu Cell Proliferation Kit with Alexa Fluor 488 (#C0071S; Beyotime), according to the manufacturer’s instructions. Notably, we digested and collected all SK-N-BE2 cells before conducting the Edu assays due to the suspended growth of some SK-N-BE2 cells. After Edu staining, cells were visualized via a fluorescence microscope (Ti2-E PFS, Nikon, Japan), and the percentage of Edu-positive cells was calculated as the number of Edu-positive cells out of the total number of cells (100$\times$).

Trans well chambers were used to detect cell migration. Cells (4 $\times$ 10⁴) in 1% FBS MEM were seeded in the upper chamber, while 10% FBS MEM was added to the lower chamber as a chemoattractant. After being cultured for 8 h, the cells were fixed in 4% paraformaldehyde for 15 min and stained with gentian violet for 10 min. At least five areas were selected for imaging under a microscope.

For cell cycle detection, the indicated cells were harvested and fixed in 75% ethanol at 4 °C for 24 h. The cell lysate was washed with phosphate-buffered saline, stained with RNase A and propidium iodide (#C1052; Beyotime) at 4 °C for 30 min, and then analyzed using the CytoFLEX flow cytometer (Beckman Coulter, Brea, CA, USA). For the detection of cell apoptosis, the indicated cells were harvested and resuspended in binding buffer, followed by Annexin-V-FITC and PI staining with an Annexin-V-FITC/PI Apoptosis Detection Kit (#40302ES60; Yeasen).

Statistical analyses were conducted for discrete variables using SPSS software (version 25.0; IBM Co., Armonk, NY, USA) and R software (version 4.1.2) and related packages for data processing, analysis, presentation, and p-value calculations. Statistical significance was defined as p 0.05. Between-group comparisons of continuous variables were performed using the Mann-Whitney Wilcoxon test on results from quantitative PCR analysis.

Based on the public dataset GSE49710 and its association with adverse clinical outcomes (INSS stage 4, COG high-risk group, patient mortality, and MYCN amplification), we performed WGCNA. This identified two gene modules significantly correlated with poor prognosis in patients with NB (Fig. 1A,B). Given our group’s prior focus on cell cycle-related genes (Fig. 1C), we intersected genes from these four WGCNA modules with a predefined cell cycle gene set, yielding hub genes (Fig. 1D). Expression analysis revealed that these hub genes were consistently upregulated in deceased patients (Fig. 1E), suggesting their strong association with adverse clinical outcomes.

Fig. 1.

Groups. (A) After the division of WGCNA into different modules, the similarity between individual modules was visually represented by color intensity: a darker color indicates greater similarity between the two modules. (B) The midnight blue and black modules contain genes corresponding to adverse clinical prognostic outcomes. There are eight groups of genes, for a total of 360. (C) The set of 30 genes identified earlier, which we call “correlated genes”, has been mapped into a protein–protein interaction network. In this network visualization, the BC value determines the color intensity: a darker color indicates a higher BC value, placing the corresponding protein closer to the network’s core. (D) A Venn diagram analysis was used to identify 21 hub genes, defined as those genes present in the intersection of all 8 gene groups and their correlated counterparts. (E) All samples in GSE49710 and TARGET, based on two final clinical outcomes (survival and death) were divided into two groups, and hub genes showed differential expression in both groups. (F) All samples in GSE49710 and TARGET were divided into two non-overlapping categories based on hub genes, C1 and C2 groups. (G) KM survival analysis comparing the C1 and C2 groups within the GSE49710 and TARGET datasets, respectively. (H) The GSVA results of the C1 and C2 groups. WGCNA, weighted gene co-expression network analysis; BC, Betweenness Centrality; TARGET, Therapeutically Applicable Research to Generate Effective Treatments; KM, Kaplan-Meier; GSVA, gene set variation analysis.

Using consensus clustering based on hub genes, we stratified all patients from the GSE49710 and TARGET datasets into two distinct subtypes (Fig. 1F). Kaplan-Meier (KM) survival curves demonstrated significantly divergent survival probabilities between these subtypes (p 0.001; Fig. 1G). GSVA further confirmed differential cell cycle activity between subtypes (Fig. 1H).

To construct a more precise prognostic model, we applied univariate Cox regression analysis to identify 131 differential genes (Supplementary Table 1). These genes were further analyzed using linear and tree models, ultimately resulting in the selection of 24 genes (Supplementary Table 2). Subsequent LASSO regression refined these to 14 genes (AC015818.5, AC090505.2, AC123912.4, C21orf58, CNR1, fatty acid hydroxylase domain containing 2, H2A histone family member X, LAMA4, LEPR, meiotic double-stranded break formation protein 4, retrotransposon-like protein 1, SLC25A10, TERT, Z82196.1), with the “3+X” method proving optimal for feature selection (Fig. 2A). The final coding RNA-based risk score formula was as follows (Fig. 2B,C):

Fig. 2.

Construction of a prognostic model. (A) Eighty ML algorithms were combined to screen for characteristic genes. (B,C) Multicox forest plot and calculated coefficients for the characteristic genes. (D–F) KM curves of the training cohort and validation cohort and the AUC values. (G,H) Performance of the risk score was compared to other clinical and molecular variables in three cohorts. ML, Machine learning; AUC, area under the curve.

Risk score =

(0.035108026) $\times$ [*C21orf58*] + (–0.095135295) $\times$ [*CNR1*] + (–0.025481523) $\times$ [*LAMA4*]

+ (–0.056961113) $\times$ [*LEPR*] + (0.040014248) $\times$ [*SLC25A10*] + (0.027990096) $\times$ [*TERT*]

In both the training (TARGET) and validation cohorts (GSE49710, AE), KM curves showed significantly longer survival in low-risk versus high-risk groups (p 0.0001, Fig. 2D–F). The model exhibited robust predictive accuracy, with AUCs of 0.76 (TARGET), 0.81 (GSE49710), and 0.83 (AE). Risk scores outperformed other clinical variables in outcome prediction (Fig. 2G,H).

Using a threshold of –0.81, we stratified patients into low-risk and high-risk groups. An online calculator was developed for clinical application (http://biowinford.site:3838/Cellcycle_RiskScore_of_NB/).

For example, inputting expression values (0.035, 0.095, 0.025, 0.056, 0.040, 0.028) yielded a risk score of 0.0164, classifying the patient as high-risk with an associated survival curve (Fig. 3A).

Fig. 3.

New drug. (A) Screenshot of online prognostic assessment platform (http://biowinford.site:3838/Cellcycle_RiskScore_of_NB/). (B) Eugenol was tested against three drug targets: ESR1, ESR2, and the AR. Both MYC and CCND1 achieved the highest BC values, with a darker colour indicating a more central position within the network. (C) The 12 groups of samples were sequenced by bulk RNA and then calculated using the online platform and divided into high- and low-risk groups. ESR1, estrogen receptor 1; AR, androgen receptor; MYC-CCND1, MYC-cyclin D1.

To identify targeted therapies for high-risk patients, proximity-based drug screening identified eugenol, which targets estrogen receptor 1 (ESR1), ESR2, and the androgen receptor (AR). Betweenness Centrality (BC) analysis revealed MYC and CCND1 as core nodes in its network (Supplementary Table 3), suggesting their key roles in eugenol’s mechanism (Fig. 3B).

We validated the risk model in 12 NB cell lines using bulk RNA-seq (Fig. 3C), focusing on how risk scores correlated with MYCN and c-MYC expression levels. SH-SY5Y cells​​ (MYCN non-amplified with normal c-MYC expression) scored below –0.81 (indicating low risk) across all four replicates. SK-N-BE2 cells (MYCN-amplified with high c-MYC expression) scored above –0.81 (indicating high risk) across all four replicates. SK-N-SH cells (MYCN non-amplified with inducible c-MYC expression) showed mixed results, with two replicates scoring above –0.81 (high risk), and the other two scoring below –0.81 (low risk) (Fig. 4A,a). These results are in line with our knowledge and expectations for the three cell lines.

Fig. 4.

Mechanism. (A) (a) Grouping of cells following analysis of bulk RNA seq data. (b–e) Visualization of cell visibility or viability across a range of concentration gradients: 100, 80, 60, 50, 25, 12, 6, 3, and 1.5 µM. SH-SY5Y/low-risk SK-N-SH: 66.0 $\pm$ 2.0 µM and 54.4 $\pm$ 1.2 µM (mean $\pm$ SD). SK-N-BE2/high-risk SK-N-SH: 33.5 $\pm$ 1.2 µM and 47.4 $\pm$ 3.5 µM. (f) Comparison of IC₅₀ results between different groups of cells. (B) Western blotting revealed stepwise upregulation of c-MYC and CCND1 protein levels across cell lines with increasing risk scores. The c-MYC Western blot band was excised at 60–80 kDa, the tubulin Western blot band was excised at 45–55 kDa, and the CCND1 Western blot band was excised at 30–40 kDa. (C) After treatment of high-risk cells (BE2 and SH-H) with the CDK4/6 inhibitor abemaciclib mesylate, CDK4/6 expression was reduced, and G1-phase arrest was increased. (D) (a,b) The GAPDH Western blot band was excised at 35–40 kDa, the CDK2 Western blot band was excised at 30–35 kDa, and the CDK4 Western blot band was excised at 30–35 kDa. (E) (a–c) Visualization of cell visibility or viability across a range of concentration gradients: 50, 40, 30, 20, 15, 10, 7.5, 5, and 2.5 µM. (d–f) Visualization of cell visibility or viability across a range of vincristine concentration gradients: 10, 5, 2.5, 1, 0.5, 0.3, and 0.1 µM. (*p 0.05, **p 0.01, ***p 0.001, ****p 0.0001, ^nsp $\mathrm{>}$ 0.05). BE2 stands for SK-N-BE2 cells; SH-H represents a high-risk group of SK-N-SH cells; SH-L represents SK-N-SH cells in the low-risk group; 5Y represents SH-SY5Y cells; NC, Normal control cells; Abe, Abemaciclib mesylate; Cis, Cisplatin; Vin, Vincristine.

Since SK-N-BE2 cells and some SK-N-SH cells were classified as high-risk groups, and SH-SY5Y cells and the remaining SK-N-SH cells were classified as low-risk groups, IC₅₀ values were measured for each of these four cells. The IC₅₀ detection of SH-SY5Y cells and SK-N-SH cells in the low-risk group was repeated three times; and the results were 66, 70, and 69 µM, respectively (for SH-SY5Y cells) and 55, 55.23, and 53 µM, respectively (for SK-N-SH cells) (Fig. 4A,b,c).

The IC₅₀ assay of SK-N-BE2 cells and SK-N-SH cells in the high-risk group was repeated three times, and the results were 33, 32.61, and 35 µM, respectively (SK-N-BE2 cells); 43, 49.31, and 50 µM, respectively (SK-N-SH cells) (Fig. 4A,d,e). The statistical analyses showed that SK-N-BE2 cells were the most sensitive to eugenol, followed by SK-N-SH cells in the high-risk group (Fig. 4A,f). These results reveal the killing effect of eugenol on tumor cells of NB, which is worth further study.

Western blotting indicated progressive increases in c-MYC and CCND1 protein expression with increasing risk scores. Specifically, their expression levels increased gradually across different cell lines: lowest in SH-SY5Y cells, intermediate in SK-N-SH in the low-risk group, higher in SK-N-SH cells in the high-risk group, and highest in SK-N-BE2 cells. c-MYC protein expression was only highly expressed in SK-N-BE2 cells, whereas CCND1 protein was highly expressed in both SK-N-BE2 cells and SK-N-SH cells in the high-risk group. This differential expression pattern suggests that these proteins may be the primary targets through which eugenol exerts its effects (Fig. 4B).

The c-MYC protein increased the levels of the CCND1 protein by both promoting its gene expression directly and by stabilizing it indirectly through non-coding RNA. Because CCND1 forms a complex with the cell cycle proteins CDK4 or CDK6, the CDK4/6 inhibitor abemaciclib mesylate was used. Treatment with this inhibitor resulted in a detected decrease in the expression of CDK4 and CDK6 in high-risk SK-N-SH and SK-N-BE2 cell lines (Fig. 4C). The observed increase in the proportion of cells in the G1 phase of the cell cycle (Fig. 4D) confirmed target engagement. This finding lays the foundation for further research into whether CCND1 is critical to the onset of eugenol’s effects.

To compare eugenol’s efficacy against the conventional chemotherapy drugs cisplatin and vincristine, we measured and compared their respective 48-hour IC₅₀ values in SK-N-BE2, SK-N-SH, and SH-SY5Y cell lines. Interestingly, we found that the IC₅₀ of the two drugs was higher in SK-N-BE2 cells than in the other two cell lines (Fig. 4E), which indicates that high-risk NB tumors are resistant to commonly used medications. Thus, it is necessary to explore drugs that are more sensitive to high-risk NB. Edu incorporation assays demonstrated that eugenol could inhibit the proliferation of NB tumors, even in the low-risk group of SK-N-SH cells (Fig. 5A,a). In addition, the effect of eugenol on cell proliferation in SK-N-SH cells and SK-N-BE2 cells in the high-risk group was better than that of the conventional chemotherapy drugs cisplatin and vincristine; however, there was no such effect in SK-N-SH cells in the low-risk group (Fig. 5A,b). It is reasonable to infer that the c-MYC and CCND1 axes play a central role in the onset of eugenol, and after inhibiting the action of CCND1 with abemaciclib mesylate, it was found that eugenol could not continue to inhibit the proliferative state of cells (Fig. 5A,c).

Fig. 5.

Phenotype. (A) (a) Effect of Eug/Cis/Vin on the proliferation of BE2/SH-H/SH-L cells. (b) Effect of Eug on the cell proliferation ability of BE2/SH-H/SH-L cells compared with that of Cis/Vin. (c) Effect of eugenol on the proliferative ability of cells in high-risk groups after blocking CCND1 with Abe. The scale bar: 100 μm. (B) (a) Effect of Eug/Cis/Vin on the cell migration capabilities of BE2/SH-H/SH-L cells. (b) Effect of Eug on the cell migration capabilities of BE2/SH-H/SH-L cells was compared with that of Cis/Vin. (c) Effect of eugenol on the migration capabilities of cells in high-risk groups after blocking CCND1 with Abe. The scale bar: 100 μm. (C) (a,b) After blocking CCND1, BE2/SH-H cell could not induce apoptosis. (*p 0.05, **p 0.01, ***p 0.001, ****p 0.0001, ^nsp $\mathrm{>}$ 0.05). BE2 stands for SK-N-BE2 cells; SH-H represents high-risk group SK-N-SH cells; SH-L represents SK-N-SH cells in the low-risk group; 5Y represents SH-SY5Y cells. NC, Normal control cells; Abe, Abemaciclib mesylate; Eug, Eugenol; Cis, Cisplatin; Vin, Vincristine.

To further confirm that the c-MYC and CCND1 axes play a central role in the onset of eugenol, we tested the migration capacity of cells by a Trans well assay. Similarly, we found that eugenol had a better effect on cell migration in SK-N-SH cells and SK-N-BE2 cells in the high-risk group than with the conventional chemotherapy drugs cisplatin and vincristine (Fig. 5B,a,b). After blocking CCND1, eugenol could not affect the migration ability of tumor cells (Fig. 5B,c).

Furthermore, after inhibiting the effect of CCND1, eugenol lost the effect of inducing apoptosis in the high-risk group (SK-N-SH cells and SK-N-BE2 cells in the high-risk group) (Fig. 5C), which further confirmed our hypothesis that the c-MYC and CCND1 axes played a central role in the onset of eugenol.