A mitochondria-related genes associated neuroblastoma signature - based on bulk and single-cell transcriptome sequencing data analysis, and experimental validation

Background: Neuroblastoma (NB), characterized by its marked heterogeneity, is the most common extracranial solid tumor in children. The status and functionality of mitochondria are crucial in regulating NB cell behavior. While the significance of mitochondria-related genes (MRGs) in NB is still missing in key knowledge.

Materials and methods: This study leverages consensus clustering and machine learning algorithms to construct and validate an MRGs-related signature in NB. Single-cell data analysis and experimental validation were employed to characterize the pivotal role of FEN1 within NB cells.

Results: MRGs facilitated the classification of NB patients into 2 distinct clusters with considerable differences. The constructed MRGs-related signature and its quantitative indicators, mtScore and mtRisk, effectively characterize the MRGs-related patient clusters. Notably, the MRGs-related signature outperformed MYCN in predicting NB patient prognosis and was adept at representing the tumor microenvironment (TME), tumor cell stemness, and sensitivity to the chemotherapeutic agents Cisplatin, Topotecan, and Irinotecan. FEN1, identified as the most contributory gene within the MRGs-related signature, was found to play a crucial role in the communication between NB cells and the TME, and in the developmental trajectory of NB cells. Experimental validations confirmed FEN1’s significant influence on NB cell proliferation, apoptosis, cell cycle, and invasiveness.

Conclusion: The MRGs-related signature developed in this study offers a novel predictive tool for assessing NB patient prognosis, immune infiltration, stemness, and chemotherapeutic sensitivity. Our findings unveil the critical function of FEN1 in NB, suggesting its potential as a therapeutic target.

1 Introduction

Neuroblastoma (NB), the most common extracranial solid tumor in children, originates from embryonic neural crest cells and accounting for 15% of all childhood cancer deaths (1, 2). Characterized by its marked heterogeneity, NB presents a varied clinical spectrum. Patients with low to intermediate-risk NB exhibit a survival rate exceeding 95%, with some cases even showing spontaneous regression without the need for therapeutic intervention; while the long-term survival rate for individuals with high-risk NB remains dismal, falling below 50% (3). There is also significant intratumor heterogeneity between cells within the same individual NB patient, and a hallmark feature of high-risk NB is the presence of multiple cell subsets (4). Therefore, individualized precision treatment is particularly important in NB.

NB is traditionally classified as an immunosuppressive “cold” tumor, characterized by low immunogenicity and a poor response to immunotherapeutic interventions (5, 6). While recent advancements in immunotherapy have significantly improved survival rates for several highly immunogenic adult solid tumors, the treatment efficacy for NB remains substantially challenged by its immunosuppressive microenvironment, with the majority of pediatric patients deriving minimal benefit from current immunotherapeutic approaches (7, 8). Consequently, identifying strategies to transform the immunosuppressive “cold” tumor into an immunostimulatory “hot” tumor, conducive to tumor immune microenvironment (TIME) activation, represents a critical and urgent task for enhancing the efficacy of immunotherapy in the clinical management of NB.

Mitochondria are increasingly recognized for their critical roles in the etiology and advancement of malignant tumors, acting through a plethora of mechanisms (9, 10). Their status and functionality are crucial in regulating tumor cell apoptosis, cell cycle progression, metabolic pathways, and so on (11, 12). The interaction between tumor cells and the tumor microenvironment (TME) is also modulated by mitochondrial dynamics, which extends to affecting the efficacy of immune cells within the TME, facilitating immune evasion, and contributing to the development of resistance to treatments (13, 14). The significance of mitochondria-related genes (MRGs) in malignancies, including but not limited to NB, is evident through their substantial impact on patient prognosis (15–17). Research into MRGs-related prognostic signatures in cancers such as lung adenocarcinoma, stomach adenocarcinoma, and breast cancer has shown promising results (18–20). However, the exploration of such prognostic models in NB is still absent, underscoring a critical gap in current knowledge and presenting a clear opportunity for groundbreaking contributions to personalized cancer therapy.

This study embarks on constructing a prognostic model for NB using MRGs through a series of bioinformatics methods and machine learning algorithms, aiming to categorize patients for more targeted clinical management and therapeutic strategies. This study delves into the application of the MRGs-related signature to delineate the TIME of NB patients, assessing tumor cell stemness, and evaluating chemotherapy drug sensitivity. This comprehensive approach seeks to enhance precision in patient classification, thereby facilitating clinical benefits. Moreover, through single-cell transcriptomic analysis and experimental validation, this research explores the significant role of FEN1, the most critical molecule within the MRGs-related signature, in NB, suggesting FEN1 as a potential therapeutic target and offering new avenues for treatment strategies.

2 Materials and methods

2.1 Data sources

The data of bulk RNA sequencing in GSE49710 (21) and single cell RNA sequencing in GSE137804 (22) were acquired from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). The microarray data E-MTAB-8248 was obtained from ArrayExpress database (https://www.ebi.ac.uk/biostudies/arrayexpress) (23). The genomic data of NBL (neuroblastoma) project in TARGET (Therapeutically Applicable Research to Generate Effective Treatments) database was downloaded from the Genomic Data Commons (GDC, https://portal.gdc.cancer.gov) (24). The Supplementary Table 1 presents the clinical baseline characteristics of the 4 datasets included in this study. The list of 2,030 MRGs was derived from the study by J. Chang, et al. (19).

2.2 Screening of differentially expressed genes

The “limma” package was used to screen differentially expressed genes (DEGs) in this study (25). Linear models were fitted using the lmFit function of the “limma” package and subsequently empirical Bayesian methods were applied using the eBayes function to stabilize the variance estimates. The P value < 0.05 and | log2 Fold change (FC) | > 1 were defined as the threshold for DEGs in this study.

2.3 Unsupervised clustering

The consensus clustering method and “ConsensusClusterPlus” package was performed to discover stable and consistent cluster structures in this study (26). One to nine clustering iterations were implemented in the datasets. In each iteration, the data are randomly split into subsets and then the K-means clustering algorithm is applied. Using the consistency matrix, the optimal number of clusters was determined by evaluating the consistency and stability of clusters under different cluster numbers (27).

2.4 Survival analysis

This study employed the “survival” package to conduct Kaplan-Meier (K-M) survival analysis (28, 29), a non-parametric method used to estimate the survival function from time-to-event data. And the “survminer” package was used to visualize survival estimates and generate the survival curves in this study (30).

2.5 Construction of prognostic signature

Prior to constructing a prognostic signature, this study initially employed K-M survival analysis and univariate Cox proportional hazards regression, both with OS as the endpoint, to screen for genes significantly associated with prognosis. In univariate COX regression models, the criteria were: P value < 0.05, and the 95% confidence interval (CI) of hazard ratio (HR) was consistently distributed ipsilateral to 1. Following the preliminary survival analysis, the machine learning algorithm Least Absolute Shrinkage and Selection Operator (LASSO) regression was utilized in GSE49710 dataset to determine the optimal number of genes and their respective coefficients for the prognostic model. The LASSO method was conducted using the “glmnet” package in R (31, 32). A score was obtained by linearly combining the mRNA expression levels of selected genes, each weighted by their respective coefficients derived from the LASSO regression analysis, which is termed the mitochondria-related risk score (mtScore).

2.6 Validation of prognostic signature

In order to evaluate the signature’s predictive capability across diverse patient populations, the same statistical method was applied to the GSE49710 dataset for internal validation and to two independent datasets E-MTAB-8248 and TARGET-NBL for external validation. K-M survival analysis, receiver operating characteristic (ROC) curve analysis and correlation analysis of key clinical features based on mtScore were all validated in the above different datasets. The area under the curve (AUC) value of ROC curves exceeded 0.70 was considered to be efficient prediction (33).

2.7 Analysis of immune infiltration

For the comprehensive evaluation of the immune cell infiltration within the TME, 4 prominent computational methods were employed: ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumour tissues using Expression data) (34), EPIC (Estimating the Proportions of Immune and Cancer cells) (35), MCPcounter (Microenvironment Cell Populations-counter) (36), and CIBERSORT (Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts) (37). In addition, we also evaluated the infiltration of 28 kinds of immune cells provided by the study from Q. Jia, et al. (38). The infiltration abundance in TME of 28 different types of immune cells was calculated using the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm as described by D. A. Barbie, et al. (39).

2.8 Calculation of mRNA expression-based stemness index

To quantify the degree of NB cellular dedifferentiation, which is indicative of stemness characteristics within tumor samples, we employed the mRNA expression-based stemness index (mRNAsi). This index was calculated following the methodology developed by T. M. Malta, et al., leveraging a machine learning model predicated on the one-class logistic regression (OCLR) algorithm (40). The gene expression profile of GSE49710 used in this study was mapped against the stemness signature to calculate the mRNAsi score for each sample. The mRNAsi scores range from 0 to 1, with higher values indicating a closer resemblance to the pluripotent state, thereby suggesting higher tumor cell stemness.

2.9 Assessment of chemotherapeutic response

To evaluate the predictive value of mtRisk for drug sensitivity, the Genomics of Drug Sensitivity in Cancer (GDSC, https://www.cancerrxgene.org/) database was utilized to analyze the drug response of patients with varying mtRisk levels to 3 commonly used clinical drugs for NB patients (Cisplatin, Topotecan, and Irinotecan) (41). The half-maximal inhibitory concentration (IC50) serving as a gauge for drug potency computed by the he “DrugResponse” package (42).

2.10 Single-cell data pre-processing and single-cell communication analysis

In this study, single-cell RNA sequencing data from 16 samples, comprising 160,847 cells, were utilized for single-cell analysis. Initial data processing included quality control measures, notably the exclusion of cells characterized by an exceptionally low number of detected genes or elevated mitochondrial gene expression, followed by normalization and mitigation of batch effects. Subsequent data analysis involved dimensionality reduction using the Uniform Manifold Approximation and Projection (UMAP) algorithm, as implemented in the “Seurat” package in R (43). Cell types were annotated based on the cell markers recommended within the GSE137804 dataset.

To investigate the expression patterns of the FEN1 gene within tumor cells, all tumor cells were categorized into two groups based on the median expression level of FEN1: FEN1 high expression group (FEN1-High) and FEN1 low expression group (FEN1-Low). Inter-cellular communication was analyzed separately for FEN1-High and FEN1-Low tumor cells in relation to other cells within the TME. This analysis was facilitated using CellPhoneDB software (version 2.0; Wellcome Sanger Institute, Hinxton, Cambridge, UK) (44).

2.11 Single-cell pseudotime trajectory analysis

This study further explored the dynamic expression of FEN1 within the developmental trajectory of NB tumor cells. This study employed pseudotime trajectory analysis to simulate the continuum of cell differentiation states and to chart the progression of tumor cells from their origin to mature states. For this analysis, the “Monocle” package in R was utilized (45). Cells were ordered in a pseudotime sequence, an inferred temporal continuum that represents the maturation or progression of cells through a developmental pathway, based on their gene expression profiles. FEN1 expression levels were then quantitatively assessed across the pseudotime to elucidate the gene’s dynamic expression patterns during the development and differentiation of tumor cells. This differential expression analysis across pseudotime states aimed to identify significant changes in FEN1 expression, employing methods incorporated within the Monocle framework.

2.12 Overexpression and knockdown of FEN1 in NB cell

Human NB cell SH-SY5Y purchased from Meisen Chinese Tissue Culture Collections were cultured in Minimum Essential Medium/Ham’s F12 (MEM/F12) medium supplemented with 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin (PS), under a humidified atmosphere containing 5% CO₂ at 37°C to ensure healthy cell growth.

The overexpression (OE) and knockdown (KD) of FEN1 were achieved through infection procedure. Five distinct groups were structured in this study: Vector, FEN1, Scramble, sh-FEN1#1, and sh-FEN1#2. The pCDH-puro-FEN1 plasmid encoding human FEN1 was defined as OE of FEN1 and named “FEN1”. The empty pCDH-CMV-MCS-EF1-puro plasmid was defined as a negative control (NC) for OE (named “Vector”). Two short hairpin RNAs (shRNAs) targeting FEN1 was synthesized and inserted into the plasmid to generate the FEN1 KD vectors defines plsi-puro-FEN1–1 and plsi-puro-FEN1–2 (named “sh-FEN1#1” and “sh-FEN1#2”). The shRNAs sequences are shown in Supplementary Table 2. The NC for KD was constructed from a scrambled shRNA inserted into plsi-ctrl-puro plasmid named “Scramble”. The cloned plasmids and packaging plasmids (psPAX2 and pMD2-VSVG) were transfected into 293 T cells to synthesize the lentiviral particles. The NB cells were infected with the collected lentiviral particles. The total RNA and total protein of the infected NB cells were collected, and the efficiency of FEN1 OE and KD was verified by quantitative real-time PCR (qRT-PCR) and Western blot analysis.

2.13 Cell counting kit-8 assay

The NB cells were seeded in 96-well plates and cultured at 5% CO₂, 37°C atmospheres. The incubation was continued for 3 hours after adding 10 μL of cell counting kit-8 (CCK-8) solution (CA1210, Solarbio, China) to each well. The absorbance at 450 nm was measured at four distinct time points: 0-, 24-, 48-, and 72-hours post-adherence, using a microplate reader. The relative cell proliferation activity was calculated according to the following formula:

$$\text{Relative proliferation}=\left({\mathit{OD}}_{450mm}^{OEorKD}-{\mathit{OD}}_{450mm}^{blank}\right)/\left({\mathit{OD}}_{450mm}^{NC}-{\mathit{OD}}_{450mm}^{blank}\right)$$

2.14 Plate cloning assay

The single-cell suspensions were prepared using 0.25% trypsin-EDTA solution. The NB cells were seeded in 6-well plates at a density of 1,000 cells per well. After one week of culture, colonies were fixed with 4% paraformaldehyde for 15 minutes at room temperature and subsequently stained with 0.1% crystal violet for 15 minutes. Excess stain was removed by washing the plates with distilled water, and the plates were allowed to air dry. Colonies consisting of more than 50 cells were counted manually under a light microscope.

2.15 Mitochondrial membrane potential ΔΨm assay with JC-1

The mitochondrial membrane potential (MMP) is a critical parameter in the regulation of cell apoptosis, serving as a key indicator of cell health (46). The decline of MMP is a hallmark event in the early stage of apoptosis. The JC-1, an ideal fluorescent probe widely used to detect the MMP ΔΨm (47), was applied as an indicator of apoptosis in this study (C2003, Beyotime, China). In detail, each of the 5 groups of NB cells were trypsinized, collected, and washed twice with cold phosphate-buffered saline (PBS). Cells were resuspended in 500 µL of PBS and subsequently mixed with 500 µL of JC-1 staining solution. The mixture was then incubated for 20 minutes at 37°C in the dark to allow for staining. After incubation, cells were washed twice with dye buffer and immediately analyzed by flow cytometry. For the detection of JC-1 monomers, the analysis conducted through the FITC channel. Conversely, the assessment of JC-1 aggregates, was performed with the PE channel for detection. A minimum of 20,000 events were recorded for each sample.

2.16 Cell cycle analysis

The NB cells in 5 groups were harvested and washed twice with cold PBS. Cells were then fixed in 70% ethanol at 4°C overnight. After fixation, cells were washed with PBS and then resuspended in 100 µL of RNase A solution (CA1510, Solarbio, China) and incubated at 37°C for 30 minutes. The cells were stained with propidium iodide (PI) (CA1510, Solarbio, China) for 30 minutes at 4°C in the dark. For the analysis of DNA content, the emitted fluorescence of PI-stained cells was detected in the PE channel using a flow cytometer. Data acquisition was performed for at least 20,000 cells per sample to ensure statistical relevance. The Dean-Jett-Fox model, a built-in algorithm within FlowJo, was employed to fit the DNA content histogram and quantitatively assess the proportions of cells in G0/G1, S, and G2/M phases of the cell cycle (48).

2.17 Transwell invasion assay

The invasive potential of the NB cells was assessed using Transwell permeable supports with 8.0 µm pore polystyrene membrane inserts. The upper surface of the insert was coated with 50 µL of Matrigel at a concentration of 2 mg/mL and allowed to solidify at 37°C incubator for 1 hour to form a thin layer of matrix barrier mimicking the extracellular matrix. Then, 1x10⁵ cells in 200 µL of serum-free medium were placed into the upper chamber, and 600 µL of medium containing 10% FBS was added to the lower chamber as a chemoattractant. After 48 hours of incubation at 37°C, non-invading cells on the upper surface of the membrane were gently removed with a cotton swab. Cells that had invaded through the Matrigel and reached the lower surface of the membrane were fixed with 4% paraformaldehyde, stained with 0.1% crystal violet, and counted under a light microscope in five randomly selected fields per well.

2.18 Statistical analysis

Continuous variables were expressed as mean ± standard deviation (SD). For comparisons between two groups, the Student’s t-test was employed. Categorical variables were presented as numbers (percentages) and analyzed using the Chi-square test. Correlations between continuous variables were evaluated using Pearson’s correlation coefficient. The P value < 0.05 was considered statistically significant for all tests. All experiment were performed in triplicate.

The bioinformatics analysis was carried out using R software (version 4.3.3; R Foundation for Statistical Computing, Vienna, Austria). Post-acquisition flow cytometry data processing was conducted using FlowJo (version 10.8.1; BD Biosciences, San Jose, California USA). Part of the statistical analysis and the generation of corresponding figures were performed with GraphPad Prism (version 9.0; GraphPad Software, San Diego, California USA). Image processing and assembly tasks were accomplished using Adobe Photoshop 2023 and Adobe Illustrator 2023 (Adobe Systems Incorporated, San Jose, California USA).

3 Results

3.1 MRGs-based clustering of NB patients into 2 distinct clusters with unique differences

In this study, we identified 1,694 DEGs between NB patients with and without MYCN amplification in GSE49710 dataset, employing a threshold of P value < 0.05 and |log2FC| > 1, with 730 up-regulated genes and 964 down-regulated genes (Figure 1A). Intersection of these 1,694 MYCN status DEGs with 2,030 MRGs yielded 105 MRGs specifically relevant to NB (Figure 1B). Unsupervised consensus clustering based on the expression profiles of 105 MRGs in the GSE49710 dataset stratified 498 NB patients. The consensus cumulative distribution function (CDF) plot suggested that the optimal k value was 2 (Figure 1C). Consequently, 498 NB patients were categorized into 2 clusters: Cluster A with 361 patients and Cluster B with 137 patients (Figure 1D).

Figure 1

Figure 1 Consensus clustering in GSE49710 based on MRGs associated with NB. (A) Volcano plot displaying DEGs between MYCN-amplified and non-amplified NB patients in GSE49710 (Genes with P value < 0.05 and |log2FC| > 1 are highlighted). (B) Venn diagram illustrating the intersection of MYCN status DEGs with MRGs, identifying 105 MRGs specifically associated with neuroblastoma. (C, D) Consensus clustering of NB patients into clusters A and B based on the expression of 105 MRGs, with k=2 as the optimal cluster number. MRGs, mitochondria-related genes; NB, neuroblastoma; DEGs, differentially expressed genes; FC, fold change; CDF, cumulative distribution function.

Significant disparities were observed between patients in Clusters A and B in terms of survival, clinical characteristics, and immune cell infiltration. Principal component analysis (PCA) distinctly separated the 2 clusters, validating the classification robustness (Figure 2A). The expression heatmap of the 105 MRGs in Clusters A and B was showed in Supplementary Figure 1. The K-M survival analysis indicated that NB patients in Cluster B had significantly worse overall survival (OS) compared to those in Cluster A (P < 0.001) (Figure 2B). Further analysis of clinical features showed substantial statistical differences between the 2 clusters in key clinical indicators (Figure 2C). Patients for progression and INSS stage 4 (an independent risk factor for NB) (49) were predominantly found in Cluster B (P < 0.0001 for both). Similarly, patients with clinical risk factors, MYCN amplification, and age below 18 months (a factor associated with poorer prognosis) were significantly concentrated in Cluster B (P < 0.0001 for all comparisons). Detailed relations between each clinical characteristic and distribution across Clusters A and B are depicted using Sankey diagrams in Supplementary Figure 2. Additionally, a marked difference in the TME between Clusters A and B was uncovered through 4 distinct immune infiltration analysis algorithms. The ESTIMATE algorithm suggested that patients in Cluster A had higher scores overall in terms of ESTIMATE, immune, and stromal scores compared to Cluster B (Figure 2D). Results from the EPIC algorithm revealed significant statistical differences in cell proportions of all 7 cell types between the two clusters (Figure 2E). In parallel, the analysis using MCPcounter indicated that the cell abundance of T cells, Cytotoxic lymphocytes, B lineage, NK cells, Monocytic lineage, and Myeloid dendritic cells in Cluster A was statistically higher compared to Cluster B (Figure 2F). Finally, the CIBERSORT analysis also reflected that Cluster A patients exhibited higher cell proportions for various immune cells compared to Cluster B (Figure 2G).

Figure 2

Figure 2 Significant disparities between Cluster A and B in GSE49710. (A) The PCA scatter plot demonstrating the segregation of NB patients into 2 clusters, Cluster A (blue) and B (red), based on MRG expression profiles. (B) The Kaplan-Meier survival curves depicting the OS probability for patients in clusters A and B. (C) The heatmap displaying the distribution of progression status, INSS stage, clinical risk, MYCN amplification status, and age in patients within Clusters A and B. (D) Box plots representing the ESTIMATE scores in Clusters A and B. (E) Box plots illustrating the proportion of various immune cells as analyzed by the EPIC algorithm in Clusters A and B. (F) Box plots depicting the cell abundance of different immune cell types as analyzed by the MCPcounter algorithm. (G) Box plots detailing the cell proportion of various immune cells as analyzed by the CIBERSORT algorithm in Cluster A and B. PCA, principal component analysis; NB, neuroblastoma; MRG, mitochondria-related genes; OS, overall survival; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumour tissues using Expression data; EPIC, Estimating the Proportions of Immune and Cancer cells; MCPcounter, Microenvironment Cell Populations-counter; CIBERSORT, Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts. (*P<0.05, **P<0.01, ***P<0.001, ****P<0.0001).

To further validate the broad applicability of the clustering based on MRGs, the same cluster analysis was conducted on the E-MTAB-8248 dataset. The consensus CDF plot suggested an optimal number of clusters k = 2 (Figure 3A), subdividing the 223 NB children in the dataset into 138 in Cluster A and 85 in Cluster B (Figure 3B). The heatmap of the expression of 105 MRGs for Clusters A and B within the E-MTAB-8248 dataset is presented in Supplementary Figure 3. The PCA confirmed robust separation between the clusters (Figure 3C). Echoing the results from the GSE49710 dataset, K-M survival analysis within E-MTAB-8248 also demonstrated poorer OS for patients in Cluster B (P<0.001) (Figure 3D). Substantial differences were also observed between clusters A and B in important clinical characteristics within E-MTAB-8248 (Figure 3E). Chromosome 1p aberrations, recognized as markers of poor prognosis in NB patients (49), were more frequently observed in patients of Cluster B (P<0.0001). Similarly, patients associated with poor prognostic factors such as INSS stage 4, MYCN amplification, and age under 18 months were predominantly found in Cluster B (all comparisons P<0.0001). Supplementary Figure 4 employs Sankey diagrams to detail the distribution of these significant clinical features between clusters A and B in the E-MTAB-8248 dataset.

Figure 3

Figure 3 Consensus clustering of E-MTAB-8248 dataset into 2 distinct clusters. (A) The consensus CDF plot identifying the optimal cluster number (k=2). (B) Consensus matrix heatmap at k=2, displaying the robust bifurcation of the dataset into Clusters A and B. (C) The PCA demonstrating a clear separation between the two clusters, validating the clustering approach. (D) Kaplan-Meier survival analysis revealing a significant survival disadvantage for Cluster B compared to Cluster A (p<0.001). (E) The clinical characteristic heatmap exhibiting distinct profiles between Cluster A and B, with chromosome 1p status, INSS stage, MYCN amplification status, and age at diagnosis, underscoring the clinical relevance of the clustering. Chr: chromosome; CDF, cumulative distribution function; PCA, principal component analysis. (***P<0.001, ****P<0.0001).

3.2 Construction and internal validation of the MRGs-related signature

In an effort to delineate the 2 clusters formed by MRG expression, this study constructed an MRGs-related signature to quantify the distinction through a scoring mechanism. Initially, 1,497 DEGs (618 up-regulated genes and 879 down-regulated genes) between Clusters A and B in the GSE49710 dataset (Figure 4A), and 830 DEGs (369 up-regulated genes and 461 down-regulated genes) between Clusters A and B in the E-MTAB-8248 dataset (Figure 4B) were identified (using a threshold of P < 0.05 and |log2FC| > 1). A Venn diagram depicts the 33 intersecting genes found between the 1,497 DEGs in GSE49710, the 830 DEGs in E-MTAB-8248, and the 105 MRGs specific to NB (Figure 4C). Further, the prognostic implications of the 33 intersecting genes were validated in both the GSE49710 and E-MTAB-8248 datasets, using K-M analysis with OS as the endpoint, categorized by the gene median expression. In GSE49710, K-M analysis revealed statistical survival differences between patients with high and low expression of each of the 33 genes (Supplementary Figure 5). Similarly, in the E-MTAB-8248 dataset, 31 of the 33 genes showed statistically significant survival correlations in K-M analysis (Supplementary Figure 6). Subsequent analysis involved univariate Cox regression to further screen for genes significantly associated with prognosis. Each of the 31 genes underwent univariate Cox regression analysis with OS as the endpoint in both GSE49710 and E-MTAB-8248 (Figure 4D), with only one gene (marked in red in Figure 4D) not showing statistical significance in E-MTAB-8248. Consequently, the remaining 30 genes proceeded to the next phase of analysis. The LASSO regression analysis was applied in GSE49710 to refine the selection to 10 genes, assigning coefficients to each (Figures 4E, F). The mtScore was defined as the linear combination of the mRNA expression levels of the 10 genes, weighted by their respective coefficients provided by LASSO analysis within GSE49710. The formula is as follows: mtScore = (-0.243142218 × the expression of DNM3) + (-0.173523906 × the expression of AGBL4) + (-0.156591746 × the expression of CROT) + (-0.155479922 × the expression of SLC22A4) + (-0.074329063 × the expression of TP63) + (-0.001579582 × the expression of PID1) + (0.029877545 × the expression of HK2) + (0.072465418 × the expression of DLGAP5) + (0.196643607 × the expression of TERT) + (0.570707127 × the expression of FEN1).

Figure 4

Figure 4 Construction of the MRGs-related signature. (A) Volcano plot displaying DEGs between Cluster A and B in GSE49710 (Genes with P value < 0.05 and |log2FC| > 1 are highlighted). (B) Volcano plot illustrating DEGs between Cluster A and B in E-MTAB-8248 (Genes with P value < 0.05 and |log2FC| > 1 are highlighted). (C) Venn diagram demonstrating the 33 intersecting genes found between the 1,497 DEGs in GSE49710, the 830 DEGs in E-MTAB-8248, and the 105 MRGs specific to NB. (D) Forest plots of HR for the 33 intersecting genes from GSE49710 and E-MTAB-8248 datasets, indicating their association with OS. P value < 0.05, and the 95% CI of HR was consistently distributed ipsilateral to 1 was considered statistically significant. (E) LASSO coefficient profiles of the 30 candidate genes. (F) Cross-validation for tuning parameter selection in the LASSO model used in (E). MRGs, mitochondria-related genes; DEGs, differentially expressed genes; FC, fold change; NB, neuroblastoma; HR, hazard ratio; OS, overall survival; CI, confidence interval; LASSO, Least Absolute Shrinkage and Selection Operator.

The mtScore for all 498 patient samples in the GSE49710 dataset was calculated in this study. Following the computation of mtScores, patients were dichotomized into two risk categories, low and high mtRisk, based on the median mtScore value. Internal validation of the predictive value of mtScore and mtRisk was conducted in the GSE49710 dataset. Initially, Figures 5A and B demonstrated that mtScore and mtRisk could effectively discriminate between Clusters A and B, with Cluster B patients exhibiting higher mtScores. The PCA indicated a clear distinction between high and low mtRisk groups (Figure 5C). A bipartite plot of mtScore distribution revealed a concentration of dead patients within the high mtScore group (Figure 5D). The K-M analysis with endpoints of OS (Figure 5E) and Event-Free Survival (EFS) (Figure 5F) showed that high mtRisk NB patients faring worse in both OS and EFS compared to their low mtRisk counterparts. The heatmap in Figure 5G displayed the expression patterns of the 10 genes used in mtScore calculation, with TERT, HK2, DLGAP5, and FEN1 being overexpressed in high mtRisk patients, while PID1, TP63, DNM3, AGBL4, CROT, and SLC22A4 showed lower expression in the high mtRisk group. The ROC curve illustrated the prognostic prediction capability and accuracy of mtScore. The AUCs of the ROC curve for OS at 3, 5, and 10 years were 0.910, 0.911, and 0.907, respectively (Figure 5H), while for EFS, they were 0.824, 0.819, and 0.843 (Figure 5I). In contrast, the AUCs of the MYCN, a well-established biological indicator of poor prognosis in NB, prediction for OS at 3, 5, and 10 years were only 0.769, 0.692, and 0.672 (Figure 5J). Furthermore, the high and low mtRisk patient groups correlated well with key clinical features (Figure 5K), with significantly more patients with progression, INSS stage 4, clinical risk, MYCN amplification, and age under 18 months in the high mtRisk group (P < 0.0001 for all 5 comparisons). The detailed distribution of mtScores across different clinical features is presented as violin plots in Supplementary Figure 7.

Figure 5

Figure 5 Internal validation of the MRGs-related signature in GSE49710. (A) Combined violin and box plot with overlaid scatter plot illustrating the distribution of mtScores across clusters A and B. (B) Sankey diagram demonstrating the effective discrimination of Clusters A and B using mtRisk. (C) PCA plot showing the clear separation between the patients of high and low mtRisk groups. (D) Bipartite plot of mtScore against OS time, highlighting the survival status (Alive or Dead) of patients. (E) Kaplan-Meier survival curve for OS, between the high and low mtRisk groups. (F) Kaplan-Meier survival curve for EFS, between the high and low mtRisk groups. (G) Heatmap showing the gene expression profiles of the 10 genes constituting the mtScore across the mtRisk stratified patient groups. (H) ROC curve of mtScore for OS. (I) ROC curve of mtScore for EFS. (J) ROC curve of MYCN for OS. (K) Heatmap displaying the association of clinical features with the mtRisk groups. MRGs, mitochondria-related genes; PCA, principal component analysis; OS, overall survival; EFS, event-free survival; ROC, receiver operating characteristic; AUC, area under the curve. (****P<0.0001).

3.3 External validation of the MRGs-related signature

To ascertain the general applicability of the MRGs-related signature and its quantitative indices, mtScore and mtRisk, further external validation was undertaken in the E-MTAB-8248 and TARGET-NBL datasets. Within the E-MTAB-8248 dataset, mtScore corresponded well with the identified Cluster A and B. Both the Figures 6A and B effectively demonstrated that mtScore and mtRisk could distinguish between Clusters A and B, with Cluster B associated with higher mtScores and mtRisk. the PCA underscored that patients in the E-MTAB-8248 dataset could be well-separated into high and low mtRisk groups (Figure 6C). The bipartite distribution plot revealed that dead patients predominantly occupied the high mtScore sector (Figure 6D). K-M analysis for OS (Figure 6E) and EFS (Figure 6F) were performed in E-MTAB-8248, demonstrating that high mtRisk patients had poorer prognoses compared to their low mtRisk counterparts (P < 0.01 for both). The heatmap in Figure 6G illustrates the expression of the 10 genes comprising the mtScore in high and low mtRisk patients within the E-MTAB-8248 dataset. The ROC curves displayed the predictive capacity of mtScore, with AUC values for OS at 3 years (0.837), 5 years (0.857), and 10 years (0.864) (Figure 6H), and for EFS at 3 years (0.765), 5 years (0.775), and 10 years (0.783) (Figure 6I). In comparison, the AUCs of the MYCN prediction for OS were only 0.705, 0.674, and 0.566 at the same time points (Figure 6J). The mtRisk also showed significant statistical correlation with key clinical features in the E-MTAB-8248 dataset (Figure 6K), where more malign clinical features such as chromosome 1p aberration, INSS stage 4, MYCN amplification, and age <18 months were significantly more prevalent in the high mtRisk group (P<0.0001 for all comparisons). The distribution of mtScores across different clinical feature groups is presented as violin plots in Supplementary Figure 8.

Figure 6

Figure 6 External validation of MRGs-related signature in the E-MTAB-8248 and TARGET-NBL datasets. (A) Combined violin, box, and scatter plot demonstrating mtScore distribution across Clusters A and B in the E-MTAB-8248 dataset. (B) Sankey diagram demonstrating the correlation between mtRisk stratification and Cluster A and B designation in E-MTAB-8248. (C) PCA plot distinctly separating high and low mtRisk patient groups in the E-MTAB-8248 dataset. (D) Scatter plot of mtScore against OS time in E-MTAB-8248, distinguishing patients by survival status. (E) Kaplan-Meier survival curve for OS in E-MTAB-8248 stratified by mtRisk. (F) Kaplan-Meier survival curve for EFS in E-MTAB-8248stratified by mtRisk. (G) Heatmap of the expression of 10 genes constituting the mtScore in E-MTAB-8248, differentiated by mtRisk groups. (H) ROC curve of mtScore predictive capacity for OS in E-MTAB-8248. (I) ROC curve of mtScore predictive capacity for EFS in E-MTAB-8248. (J) ROC curve of MYCN predictive capacity for OS in E-MTAB-8248. (K) Heatmap correlating mtRisk with key clinical features in E-MTAB-8248. (L) PCA plot distinctly separating high and low mtRisk patient groups in the TARGET-NBL dataset. (M) Kaplan-Meier survival curve for OS in TARGET-NBL categorized by mtRisk. (N) Kaplan-Meier survival curve for EFS in TARGET-NBL categorized by mtRisk. MRGs, mitochondria-related genes; PCA, principal component analysis; OS, overall survival; EFS, event-free survival; ROC, receiver operating characteristic; AUC, area under the curve; Chr, chromosome. (****P<0.0001).

External validation of mtScore and mtRisk in the TARGET-NBL dataset reiterated the robust predictive power of the MRGs-related signature. PCA delineated a clear distinction between high and low mtRisk groups in the TARGET-NBL dataset (Figure 6L). K-M survival analysis confirmed the consistent predictive power of mtRisk for prognosis, with high mtRisk patients exhibiting worse outcomes in both OS (Figure 6M) and EFS (Figure 6N) endpoints (P<0.001 and P<0.01, respectively). Beyond prognosis, mtRisk was also correlated with important clinical features in the TARGET-NBL dataset (Supplementary Figure 9A). High-risk COG risk group, unfavorable histology, high MKI, MYCN amplification, INSS stage 4, and age <18 months were statistically more frequent in high mtRisk patients. Supplementary Figures 9B–G display the mtScore distributions for different clinical feature groups as violin plots.

3.4 Predictive efficacy of MRGs-related signature for TIME, stemness, and chemosensitivity

The MRGs-related signature and its quantitative markers, mtScore and mtRisk, developed in this study, not only predict the prognosis of NB patients but also show significant relevance to immune infiltration in the TME. Various algorithms were applied to assess the indication of mtScore and mtRisk towards immune infiltration within the GSE49710 dataset. The ESTIMATE algorithm indicated that patients with low mtRisk had higher scores, with higher total ESTIMATE score, immune score, and stromal score compared to patients with high mtRisk (Figure 7A). According to the EPIC algorithm, the proportion of immune cells such as CD4⁺ T cells, CD8⁺ T cells, and macrophages was statistically higher in the TIME of patients with low mtRisk compared to those with high mtRisk (Figure 7B). The MCPcounter algorithm suggested that the cell abundance of T cells, CD8⁺ T cells, cytotoxic lymphocytes, NK cells, monocytic lineage, myeloid dendritic cells, endothelial cells, and fibroblasts was higher in the low mtRisk patient group than in the high mtRisk group (Figure 7C). CIBERSORT analysis revealed statistically significant differences in the cell proportions of naive B cells, memory B cells, plasma cells, resting memory CD4⁺ T cells, follicular helper T cells, resting NK cells, activated NK cells, monocytes, M0 macrophages, M2 macrophages, resting mast cells, and activated mast cells between the high and low mtRisk groups (Figure 7D). Additionally, a heatmap displayed the correlation between mtScore, the 10 genes constituting mtScore, and 28 types of immune cells (Figure 7E).

Figure 7

Figure 7 Immune profile, stemness, and drug Sensitivity analysis related to MRGs-related signature in GSE49710. (A) Box plots representing the ESTIMATE scores stratified by mtRisk groups. (B) Box plots illustrating the proportion of various immune cells as analyzed by the EPIC algorithm, stratified by mtRisk groups. (C) Box plots depicting the cell abundance of different immune cell types as analyzed by the MCPcounter algorithm. (D) Box plots detailing the cell proportion of various immune cells as analyzed by the CIBERSORT algorithm stratified by mtRisk groups. (E) Heatmap displaying the correlation between the mtScore, the ten genes comprising mtScore, and 28 immune cell types. (F) Scatter plot demonstrating a positive association between mtScore and mRNAsi. (G–I) Violin plots illustrating the estimated IC50 values for Cisplatin (G), Topotecan (H), and Irinotecan (I), comparing high and low mtRisk groups. MRGs, mitochondria-related genes; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumour tissues using Expression data; EPIC, Estimating the Proportions of Immune and Cancer cells; MCPcounter, Microenvironment Cell Populations-counter; CIBERSORT, Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts; IC50, half-maximal inhibitory concentration. (*P<0.05, **P<0.01, ***P<0.001).

Furthermore, results demonstrated that mtScore is significantly positively correlated with tumor cell stemness. A scatter plot in Figure 7F exhibited the relationship between mtScore and mRNAsi values in the GSE49710 dataset. The scatter plot revealed a strong positive correlation, with an R value of 0.77 (P < 0.0001). A trend line drawn through the data points further emphasized this positive linear relationship.

Moreover, mtRisk may also predict the sensitivity to 3 commonly used drugs in NB patients. Results showed that patients in the low mtRisk group had statistically significantly higher IC50 values for Cisplatin (Figure 7G), Topotecan (Figure 7H), and Irinotecan (Figure 7I) compared to the high mtRisk group, suggesting increased sensitivity to these 3 drugs in the high mtRisk patients (P < 0.001 for 3 comparison).

3.5 Analysis of FEN1’s essential role through single-cell transcriptome sequencing data

This study further validated the significant function and role of FEN1 in NB, the gene with the highest contribution in the MRGs-related signature. Quality control, normalization, and batch effect removal were initially applied to single-cell transcriptomic sequencing data from 160,847 cells of 16 NB patients sourced from GSE137804. Subsequent dimensionality reduction via UMAP clustered the cells into 30 clusters (Figure 8A), and cell type annotation using known markers identified 8 cell types similar to the original GSE137804 study (22), including tumor cells, T cells, B cells, endothelium, plasmacytoid dendritic cell (pDC), Schwann cells, fibroblasts, and myeloid cells (Figure 8B). The markers used for cell annotation are presented in Supplementary Figure 10A.

Figure 8

Figure 8 Single-cell transcriptomic data analysis revealing the role of FEN1 in NB tumor cells. (A) UMAP dimensionality reduction and clustering of 160,847 cells from 16 neuroblastoma patient samples. (B) Annotated UMAP Clustering of Cell Types. (C) Cell-cell communication counts network of FEN1-high tumor cells with the surrounding microenvironment. (D) Cell-cell communication counts network of FEN1-low tumor cells with the surrounding microenvironment. (E) Heatmap illustrating significant ligand-receptor pairs between FEN1-high tumor cells and other cell types. (F) Heatmap illustrating significant ligand-receptor pairs between FEN1-low tumor cells and other cell types. (G) Pseudotime trajectory analysis of NB tumor cells, with color gradient indicating progression through pseudotime. (H) Cells are categorized into 3 developmental states along the pseudotime trajectory. (I) Trend of FEN1 gene expression across pseudotime, demonstrating its dynamic changes during tumor cell development. NB, neuroblastoma; UMAP, uniform manifold approximation and projection; pDC, plasmacytoid dendritic cell; Commun Prob, communication probability.

Tumor cells were segregated into FEN1-high and FEN1-low groups based on the median expression of FEN1. Cell-cell communication analysis further investigated differences between FEN1-high and FEN1-low tumor cells in their interactions with the TME. Compared to FEN1-high tumor cells, FEN1-low tumor cells exhibited increased and stronger communication with surrounding cells. Figures 8C, D respectively show the communication counts of tumor cells (FEN1-high) and tumor cells (FEN1-low) with adjacent cells. Supplementary Figures 10B, C display the communication weights of tumor cells (FEN1-high) and tumor cells (FEN1-low), respectively. Additionally, the exploration of interacting pairs in cell communication revealed that tumor cells primarily interact with other cells via the MIF-(CD74+CXCR4) axis; conversely, other cells predominantly communicate with tumor cells through the PTN-NCL axis (Figures 8E, F).

Pseudotime trajectory analysis illustrated the temporal expression changes of FEN1 in the development of NB tumor cells. Figure 8G depicts the pseudotime trajectory of NB tumor cells, with the cell developmental trajectory divided into 3 states in Figure 8H. The gene expression trend of FEN1 over pseudotime is shown in Figure 8I, indicating that FEN1 expression is higher in the early stages of NB tumor cell development than in later stages.

3.6 Functional validation of FEN1 in NB cell lines

The current study further investigated the important role of FEN1 in NB cell by modulating its expression through lentiviral-mediated OE or KD in human NB cells. Successful modulation of FEN1 at both RNA and protein levels was confirmed by qRT-PCR (Figure 9A) and Western blot (Figure 9B), achieving the desired OE and KD effects.

The impact of FEN1 on NB cell proliferation was assessed using the CCK-8 assay (Figure 9C). Cells with FEN1 OE showed significantly higher proliferation rates compared to the vector group (P<0.0001). Cells with FEN1 KD (sh-FEN1#1 and sh-FEN1#2), demonstrated significantly reduced proliferation compared to the scramble group (P<0.0001 for both comparisons). Plate clonogenic assays corroborated these findings (Figure 9D), where FEN1 OE increased colony formation in NB cells compared to NC (P<0.01), and KD led to a significant decrease in colony numbers (P<0.01 for both sh#1 and sh#2).

Figure 9

Figure 9 Functional validation of FEN1 OE and KD in NB cell lines. (A) qRT-PCR Validation of FEN1 OE and KD Efficiency at RNA Level in NB Cells. (B) Western Blot Confirmation of FEN1 OE and KD at Protein Level in NB Cells. (C) Line graph depicting the relative proliferation of NB cells, measured by CCK-8 assay. (D) Representative images and quantification of colony formation assay results. (E) Flow cytometry analysis of JC-1 staining to assess MMP (ΔΨm), with the Aggregates/Monomers ratio indicating changes of apoptosis. OE, overexpression; KD, knockdown; NB, neuroblastoma; MMP, mitochondrial membrane potential. (**P<0.01, ***P<0.001, ****P<0.0001).

To explore changes in apoptosis following FEN1 modulation, we quantified ΔΨm using the JC-1 dye. An increase in red/green fluorescence ratio, indicating higher J-aggregates formation, reflects higher MMP and thus, a lower level of apoptosis. Figure 9E shows that cells with FEN1 OE had a significantly higher Aggregates/Monomer ratio compared to NC (P<0.001), whereas FEN1 KD cells displayed a significant reduction, indicating increased apoptosis (P<0.01 for both sh#1 and sh#2).

Cell cycle analysis also yielded positive findings (Figure 10A), with FEN1 KD cells showing an increased proportion in the G2/M phase, suggesting G2/M arrest compared to the NC (P<0.0001 for sh#1, P<0.001 for sh#2). Cells OE FEN1 exhibited a significant increase in the S phase proportion compared to NC (P<0.0001), indicative of heightened DNA synthesis and replication activity.

Figure 10

Figure 10 Impact of FEN1 modulation on NB cell cycle progression and invasiveness. (A) Flow cytometry cell cycle analysis presenting the situation of NB cells in 5 group. (B) Transwell invasion assays with representative images and quantitative analysis displaying the invasive capacity of NB cells in 5 groups. NB, neuroblastoma; ns, not significant. (**P<0.01, ***P<0.001, ****P<0.0001).

Furthermore, FEN1 may influenced the invasive capacity of NB cells. Transwell assays revealed changes in NB cell invasion following FEN1 OE and KD (Figure 10B). Compared to NC, FEN1 OE enhanced the invasiveness of NB cells, with increased cell numbers traversing the membrane (P<0.01). Conversely, FEN1 KD reduced NB cell invasion, with fewer cells penetrating the membrane (P<0.01 for both sh#1 and sh#2).

4 Discussion

Our study has successfully constructed a prognostic signature based on MRGs and developed quantitative indices, namely mtScore and mtRisk. Through extensive internal and external validation, the MRGs-related signature exhibited superior prognostic predictive effect and value over the traditional molecular marker MYCN in NB patients. Furthermore, it demonstrated predictive capability for immune infiltration in the TME, tumor cell stemness, and sensitivity to specific chemotherapeutic agents. Furthermore, through single-cell transcriptomic analysis, we underscored the pivotal role played by FEN1, the most contributive molecule in the MRGs-related signature, in the molecular crosstalk and developmental trajectory of NB cells. The experimental validation results underscored that FEN1 expression significantly affects processes such as cell proliferation, apoptosis, cell cycle progression, and invasiveness in NB cells.

In this study, MRGs effectively stratified NB patients into 2 distinct groups with significant differences in prognosis and immune infiltration characteristics. An MRGs-related signature and its quantitative metrics, mtScore and mtRisk, were developed to characterize these 2 patient groups, revealing that lower mtScore and mtRisk are associated with significantly better prognoses. Similar works have been reported in colorectal cancer, hepatocellular carcinoma, and other cancers (50, 51), but this is the first time an MRGs-related signature has been identified and explored in NB, to our knowledge. It’s noteworthy that the predictive efficiency of this MRGs-related signature for NB patients surpasses that of MYCN, traditionally recognized as the best genetic marker for forecasting outcomes in NB (52, 53). This highlights the potential of MRGs-related signature to provide a more comprehensive understanding of NB prognosis beyond conventional markers.

In line with the above results, MRGs-related signature not only predicts prognosis, but also identifies key clinical features to a certain extent. MYCN amplification, which is closely related to the poor prognosis of high-risk NB patients (54), was found to be significantly correlated with high mtRisk in this study. INSS4 stage, as an independent risk factor for NB (49), was found to be significantly associated with high mtRisk in this study. Similarly, high mtRisk was profoundly associated with a spectrum of adverse clinical features, including clinical risk, progression, chromosomal 1p aberration, high MKI, unfavorable histology, and high COG risk. The congruence of the MRGs-related signature’s predictive capacity for both prognosis and key clinical features in NB patients underscores its utility beyond mere prognostic estimation.

The Further analysis of immune infiltration implies that patients with lower mtScore and mtRisk tend to have a more active immune environment within their tumors. This heightened immune activity is hypothesized to be a key factor contributing to their better prognosis, suggesting a link between mitochondrial function, immune engagement, and cancer outcome. Exploring the TIME and transitioning from “cold” to “hot” tumors could significantly enhance therapeutic efficacy in NB, a tumor traditionally marked by immune suppression and modest responses to immunotherapy (5). The establishment of the MRGs-related signature provides a promising avenue for identifying and modulating the TME. By targeting specific MRGs to activate the immune landscape within NB, it may be possible to convert these traditionally “cold” tumors into “hot” tumors, potentially making them more amenable to immunotherapeutic interventions. This approach, aligning with findings in other cancers, underscores the critical interplay between mitochondrial dynamics and immune responsiveness in determining cancer prognosis and treatment outcomes (14, 55). The results of this study may provide a new target for immunotherapy of NB and other tumors (56–59).

In addition, the significant linear correlation between mtScore and mRNAsi underscores the potential of our developed MRGs-related signature to effectively indicate the stemness of tumor cells in patients. Our team’s previous work has developed an mRNAsi-based risk score for NB, which demonstrated excellent performance in predicting patient prognosis, immune infiltration, and treatment response (60). This study extends the prognostic utility of mtScore as a marker of stemness, potentially impacting the clinical management of NB significantly. Patients with elevated mtScores might be identified as harboring a higher burden of tumor stem cells, likely to undergo aggressive disease progression and exhibit poor responses to standard therapies. This observation aligns with the outcomes of poor prognosis and reduced sensitivity to certain chemotherapies among patients with high mtScores. Such insights potentially could facilitate the stratification of patients into more personalized treatment regimens.

In the realm of clinical pharmacotherapy, our findings highlight a significant role of mtRisk in shaping therapeutic responses. At first glance, patients categorized under high mtRisk appear to present a formidable challenge in treatment management due to their poor prognosis, immune-suppressive TME, and pronounced tumor cell stemness. However, a pivotal discovery of our study is the heightened sensitivity of high mtRisk patients to 3 clinically prevalent drugs for NB: Cisplatin, Topotecan, and Irinotecan. All 3 drugs are internationally recognized for the treatment of NB patients (61–63). This enhanced drug responsiveness, surpassing that of low mtRisk patients, uncovers a nuanced aspect of mtRisk’s clinical implications. This holds some promise for the therapeutic management of patients with high mtRisk, and also re-emphasizes the importance of precise individualized treatment in highly heterogeneous NB patients.

The MRGs-related signature in NB is composed of ten genes: FEN1, TERT, DLGAP5, HK2, PID1, TP63, SLC22A4, CROT, AGBL4, DNM3. FEN1 (Flap endonuclease 1) is crucial in DNA replication and repair, and its overexpression is linked to poor prognosis in various cancers, indicating its potential as a target for cancer therapy (64). TERT (Telomerase reverse transcriptase), the catalytic subunit of telomerase, affects telomere length by affecting telomerase activity, and considered to be a useful marker in diagnosis and prognosis of various cancers and a new therapy approach (65). DLGAP5 (Discs large homolog associated protein 5) is involved in mitotic spindle assembly, and its overexpression is associated with tumor progression and adverse outcomes in cancer patients, highlighting its role in cell division and potential as a therapeutic target (66, 67). HK2 (Hexokinase 2) catalyzes the first step of glycolysis and its upregulation in tumors is linked to enhanced glycolytic metabolism typical of cancer cells, suggesting its involvement in the Warburg effect and as a target for metabolic therapy (68). The role of PID1 (Phosphotyrosine interaction domain-containing protein 1) in cancer involves modulating lipid metabolism and mitochondrial function, indicating its potential impact on tumor metabolic reprogramming and its association with cancer progression (69, 70). A member of the p53 family, TP63 is implicated in the development and progression of several cancers, where it can influence cell cycle regulation, apoptosis, and the immune response in the TME (71, 72). SLC22A4, a solute carrier protein, has been associated with drug disposition and response in cancer therapy, reflecting its role in modulating chemotherapeutic efficacy and resistance mechanisms in tumors (73). CROT is involved in fatty acid metabolism, and alterations in its expression are linked to changes in cancer cell metabolism and potential effects on tumor growth and patient prognosis (74). AGBL4 (ATP/GTP binding protein-like 4), as a neuronal differentiation marker, participates in neuronal differentiation by promoting mitochondrial axonal growth and axonal transport (75). DNM3 (Dynamin 3), involved in exosomes, endocytosis, and tumor metastasis, is considered as a tumor suppressor gene in a variety of cancers such as non-small-cell lung cancer (NSCLC), hepatocellular carcinoma, papillary thyroid carcinoma, and colon cancer (76–79). Therefore, the MRGs-related signature developed in this study encompasses a range of factors related to patient prognosis, the TIME, apoptosis, cell cycle progression, mitochondrial function, neuronal differentiation, and resistance mechanisms in NB. This comprehensive signature holds potential for guiding therapeutic strategies and prognosis assessment in NB patients.

The single-cell transcriptomics data analysis further explored the role of the most important gene in the MRGs-related signature, FEN1, providing profound insights into this gene’s multifaceted functions in NB. The enhanced communication between FEN1-low tumor cells and the surrounding cells, as revealed through our analysis, suggests a more dynamic interplay within the TME. Taken together with the results of bulk sequencing data analysis (FEN1 is an oncogene associated with poor prognosis in NB patients and positively contributes to mtScore), we may be able to make an important hypothesis. High FEN1 expression could lead to a reduced need for external support from the TME, suggesting that these cells might have developed autonomous signaling pathways that promote proliferation, resist apoptosis, and enhance invasion capabilities without the extensive need for stromal or immune cell interactions. This autonomy could be a factor in their increased progression and aggressiveness. The differential communication patterns, particularly through the MIF-(CD74+CXCR4) and PTN-NCL axes, may underscore specific pathways amenable to therapeutic intervention. Furthermore, pseudotime trajectory analysis elucidating the temporal changes in FEN1 expression across the developmental trajectory of NB tumor cells underscores its significance. The observation that FEN1 expression is higher in the early stages of tumor cell development suggests a role in the initial phases of tumorigenesis or in maintaining a stem-like state of the tumor cells. This finding is supported by research from Z. Peng, et al., which confirmed FEN1’s capability to promote stemness in tumor cells (80).

FEN1 is an enzyme characterized by its multifunctional enzymatic activities, playing a pivotal role in the processes of DNA replication and repair (81). The exonuclease activity of FEN1 is critical for the maturation of lagging strand Okazaki fragments during DNA replication. This activity facilitates the removal of RNA primers at the termini of these fragments, as well as the trimming of damaged ends during DNA repair mechanisms (82, 83). As an endonuclease, FEN1 recognizes and cleaves flap structures, which are intermediates formed during DNA replication. The precise excision of these flap structures by FEN1 ensures the continuity and fidelity of DNA synthesis. This endonuclease activity is particularly crucial for genomic stability, as it aids in the accurate resolution of structural anomalies encountered during DNA replication and repair processes (84). FEN1 has been implicated in the malignancy progression of various cancers, including gastric cancer, NSCLC, and cholangiocarcinoma (85–87). This study extends the understanding of FEN1’s oncogenic role by exploring its functional significance in NB cell lines through OE and KD experiments. Our findings illustrate the multifaceted role of FEN1 in modulating NB cell behaviors, emphasizing its critical involvement in cell proliferation, apoptosis, cell cycle progression, and invasiveness. The pivotal role of FEN1 in supporting NB cell proliferation and colony formation underscores its important contribution to the proliferation and growth of NB. This aligns with FEN1’s established roles in DNA replication and repair mechanisms (88). Furthermore, the regulation of FEN1 significantly impacts NB cell apoptosis and cell cycle dynamics. Increased apoptosis in FEN1 KD cells highlights a potential vulnerability that could be therapeutically exploited to induce cell death in NB cells. This observation aligns with existing research indicating FEN1’s capacity to inhibit tumor cell apoptosis (89). Cell cycle analysis revealed G2/M phase arrest in FEN1 KD cells and an increase in the proportion of cells in the S phase following FEN1 OE, further evidencing FEN1’s regulatory role in cell cycle progression. Such findings suggest that FEN1 is not only indispensable for the replication process but also crucial for the proper progression of the cell cycle, likely by ensuring the fidelity of DNA replication and repair prior to mitotic entry. Alterations in the invasive capacity of NB cells post-FEN1 modulation underscore the gene’s role in tumor metastasis. Enhanced invasiveness in FEN1 OE cells could reflect an increased ability to facilitate tumor spread. Coupled with resistance to apoptosis and increased proliferation, these capabilities position FEN1 as a key driver of NB aggressiveness and metastatic potential. In summary, our findings advocate for FEN1’s important role in NB cell proliferation, survival, and invasiveness, and emphasize that FEN1 may serve as a potential therapeutic target for NB.

This study presents a pioneering effort in constructing an MRGs-related signature for NB, revealing potential avenues for further research and clinical application. However, there are limitations to consider. Of the ten genes comprising the constructed MRGs-related signature, only one was confirmed in subsequent single-cell data analysis and experimental validation. It is imperative to highlight that the validated gene, FEN1, is not only a part of the signature but its most crucial and significantly contributory element, underscoring its pivotal role within the MRGs-related signature for NB. The roles of the remaining nine genes in NB require further exploration through more extensive data analysis and experimental investigation. Additionally, the validation of the signature and its associated genes would benefit from a broader array of NB cell lines, animal models, and clinical patient samples to confirm their function and impact more definitively. This approach underscores the importance of comprehensive validation to strengthen the findings and potential clinical applications of genomic signatures in cancer research.

In conclusion, this study shows for the first time that MRGs can divide NB patients into two clusters that differ significantly in terms of survival prognosis, clinical features, and TIME. On this basis, this study developed MRGs-related signature and its quantitative indicators mtScore and mtRisk to characterize the above two clusters. The MRGs-related signature constructed in this study can successfully distinguish heterogeneous NB patients in different clusters, which is of great significance for the targeted and precise treatment of NB patients with different characteristics. Notably, the MRGs-related signature can predict the prognosis of NB patients, and the predictive performance is better than that of MYCN. The MRGs-related signature is significantly associated with malignant clinical features including MYCN amplification status. Besides, the MRGs-related signature may indicate the immune infiltration in TIME of NB patients to a certain extent, which may be of great significance for distinguishing “hot” tumors from “cold” tumors and predicting the response of immunotherapy. The MRGs-related signature was also adept at representing tumor cell stemness, and sensitivity to the chemotherapeutic agents Cisplatin, Topotecan, and Irinotecan. Furthermore, the important role of FEN1, the most important gene in MRGs-related signature, in NB was demonstrated by single-cell data analysis and experimental validation in this study. The development of mtScore and mtRisk provides a new perspective and evidence for the precise treatment, prognosis prediction, the conversion of “cold” and “hot” tumors, and the activation of TIME of NB patients. The important role of FEN1 demonstrated in this study also provides a potential new target for the treatment of NB.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

CW: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. JT: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. YJ: Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation. ZL: Writing – review & editing, Investigation, Data curation. JY: Writing – review & editing, Visualization, Methodology. YBJ: Writing – review & editing, Visualization, Methodology. YX: Writing – review & editing, Investigation, Data curation. BG: Writing – review & editing, Supervision, Resources, Funding acquisition. QD: Writing – review & editing, Supervision, Resources, Funding acquisition. QZ: Writing – review & editing, Supervision, Resources, Funding acquisition.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study was funded by National Key Research and Development Program of China (2018YFC1313000) and Tianjin Key Medical Discipline (Specialty) Construction Project (TJYXZDXK-009A).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1415736/full#supplementary-material

References

1. Cheung NK, Dyer MA. Neuroblastoma: developmental biology, cancer genomics and immunotherapy. Nat Rev Cancer. (2013) 13:397–411. doi: 10.1038/nrc3526

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Ngan ES. Heterogeneity of neuroblastoma. Oncoscience. (2015) 2:837–8. doi: 10.18632/oncoscience.216

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lundberg KI, Treis D, Johnsen JI. Neuroblastoma heterogeneity, plasticity, and emerging therapies. Curr Oncol Rep. (2022) 24:1053–62. doi: 10.1007/s11912–022-01270–8

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Gomez RL, Ibragimova S, Ramachandran R, Philpott A, Ali FR. Tumoral heterogeneity in neuroblastoma. Biochim Biophys Acta Rev Cancer. (2022) 1877:188805. doi: 10.1016/j.bbcan.2022.188805

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kennedy PT, Zannoupa D, Son MH, Dahal LN, Woolley JF. Neuroblastoma: an ongoing cold front for cancer immunotherapy. J Immunother Cancer. (2023) 11. doi: 10.1136/jitc-2023–007798

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Wienke J, Dierselhuis MP, Tytgat GAM, Künkele A, Nierkens S, Molenaar JJ. The immune landscape of neuroblastoma: Challenges and opportunities for novel therapeutic strategies in pediatric oncology. Eur J Cancer. (2021) 144:123–50. doi: 10.1016/j.ejca.2020.11.014

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Galluzzi L, Chan TA, Kroemer G, Wolchok JD, López-Soto A. The hallmarks of successful anticancer immunotherapy. Sci Transl Med. (2018) 10. doi: 10.1126/scitranslmed.aat7807

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Anderson J, Majzner RG, Sondel PM. Immunotherapy of neuroblastoma: facts and hopes. Clin Cancer Res. (2022) 28:3196–206. doi: 10.1158/1078–0432.Ccr-21–1356

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Wallace DC. Mitochondria and cancer. Nat Rev Cancer. (2012) 12:685–98. doi: 10.1038/nrc3365

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Zong WX, Rabinowitz JD, White E. Mitochondria and cancer. Mol Cell. (2016) 61:667–76. doi: 10.1016/j.molcel.2016.02.011

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Badrinath N, Yoo SY. Mitochondria in cancer: in the aspects of tumorigenesis and targeted therapy. Carcinogenesis. (2018) 39:1419–30. doi: 10.1093/carcin/bgy148

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Pangou E, Sumara I. The multifaceted regulation of mitochondrial dynamics during mitosis. Front Cell Dev Biol. (2021) 9:767221. doi: 10.3389/fcell.2021.767221

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Bai R, Cui J. Mitochondrial immune regulation and anti-tumor immunotherapy strategies targeting mitochondria. Cancer Lett. (2023) 564:216223. doi: 10.1016/j.canlet.2023.216223

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Klein K, He K, Younes AI, Barsoumian HB, Chen D, Ozgen T, et al. Role of mitochondria in cancer immune evasion and potential therapeutic approaches. Front Immunol. (2020) 11:573326. doi: 10.3389/fimmu.2020.573326

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Chai C, Chen Y, Luo Y, Zhang H, Ye Z, He X, et al. Mitochondria-associated gene expression perturbation predicts clinical outcomes and shows potential for targeted therapy in neuroblastoma. Front Pediatr. (2023) 11:1094926. doi: 10.3389/fped.2023.1094926

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Wang Y, Song F, Zhang X, Yang C. Mitochondrial-related transcriptome feature correlates with prognosis, vascular invasion, tumor microenvironment, and treatment response in hepatocellular carcinoma. Oxid Med Cell Longev. (2022) 2022:1592905. doi: 10.1155/2022/1592905

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Liu W, Du Q, Mei T, Wang J, Huang D, Qin T. Comprehensive analysis the prognostic and immune characteristics of mitochondrial transport-related gene SFXN1 in lung adenocarcinoma. BMC Cancer. (2024) 24:94. doi: 10.1186/s12885-023-11646-z

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Li YP, Liu GX, Wu ZL, Tu PH, Wei G, Yuan M, et al. A novel mitochondrial-related gene signature for the tumor immune microenvironment evaluation and prognosis prediction in lung adenocarcinoma. J Immunol Res. (2022) 2022:5366185. doi: 10.1155/2022/5366185

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Chang J, Wu H, Wu J, Liu M, Zhang W, Hu Y, et al. Constructing a novel mitochondrial-related gene signature for evaluating the tumor immune microenvironment and predicting survival in stomach adenocarcinoma. J Transl Med. (2023) 21:191. doi: 10.1186/s12967–023-04033–6

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wang Y, Wang DY, Bu KN, Gao JD, Zhang BL. Prognosis prediction and risk stratification of breast cancer patients based on a mitochondria-related gene signature. Sci Rep. (2024) 14:2859. doi: 10.1038/s41598-024-52981-w

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Su Z, Fang H, Hong H, Shi L, Zhang W, Zhang W, et al. An investigation of biomarkers derived from legacy microarray data for their utility in the RNA-seq era. Genome Biol. (2014) 15:523. doi: 10.1186/s13059-014-0523-y

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Dong R, Yang R, Zhan Y, Lai HD, Ye CJ, Yao XY, et al. Single-cell characterization of Malignant phenotypes and developmental trajectories of adrenal neuroblastoma. Cancer Cell. (2020) 38:716–33.e6. doi: 10.1016/j.ccell.2020.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Roderwieser A, Sand F, Walter E, Fischer J, Gecht J, Bartenhagen C, et al. Telomerase is a prognostic marker of poor outcome and a therapeutic target in neuroblastoma. JCO Precis Oncol. (2019) 3:1–20. doi: 10.1200/po.19.00072

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Pugh TJ, Morozova O, Attiyeh EF, Asgharzadeh S, Wei JS, Auclair D, et al. The genetic landscape of high-risk neuroblastoma. Nat Genet. (2013) 45:279–84. doi: 10.1038/ng.2529

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Monti S, Tamayo P, Mesirov J, Golub T. Consensus clustering: a resampling-based method for class discovery and visualization of gene expression microarray data. Mach learning. (2003) 52:91–118. doi: 10.1023/A:1023949509487

CrossRef Full Text | Google Scholar

28. Therneau TM. A package For Survival Analysis in R. (2024).

Google Scholar

29. Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. New York: Springer (2000).

Google Scholar

30. Kassambara A, Kosinski M, Biecek P. Survminer: Drawing Survival Curves using 'ggplot2'. (2021).

Google Scholar

31. Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. J Stat Software. (2010) 33:1–22. doi: 10.18637/jss.v033.i01

CrossRef Full Text | Google Scholar

32. Tibshirani R. Regression shrinkage and selection via the lasso. J R Stat Soc Ser B: Stat Methodology. (1996) 58:267–88. doi: 10.1111/j.2517-6161.1996.tb02080.x

CrossRef Full Text | Google Scholar

33. Heagerty PJ, Lumley T, Pepe MS. Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics. (2000) 56:337–44. doi: 10.1111/j.0006-341X.2000.00337.x

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Racle J, Gfeller D. EPIC: A tool to estimate the proportions of different cell types from bulk gene expression data. Methods Mol Biol. (2020) 2120:233–48. doi: 10.1007/978–1-0716–0327-7_17

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. (2016) 17:218. doi: 10.1186/s13059–016-1070–5

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Jia Q, Wu W, Wang Y, Alexander PB, Sun C, Gong Z, et al. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer. Nat Commun. (2018) 9:5361. doi: 10.1038/s41467-018-07767-w

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature. (2009) 462:108–12. doi: 10.1038/nature08460

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Malta TM, Sokolov A, Gentles AJ, Burzykowski T, Poisson L, Weinstein JN, et al. Machine learning identifies stemness features associated with oncogenic dedifferentiation. Cell. (2018) 173:338–54.e15. doi: 10.1016/j.cell.2018.03.034

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, et al. Genomics of Drug Sensitivity in Cancer (GDSC): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res. (2013) 41:D955–61. doi: 10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Geeleher P, Cox NJ, Huang RS. Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol. (2014) 15:R47. doi: 10.1186/gb-2014–15-3-r47

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, et al. Comprehensive integration of single-cell data. Cell. (2019) 177:1888–902.e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc. (2020) 15:1484–506. doi: 10.1038/s41596-020-0292-x

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. (2017) 14:979–82. doi: 10.1038/nmeth.4402

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Zamzami N, Kroemer G. The mitochondrion in apoptosis: how Pandora's box opens. Nat Rev Mol Cell Biol. (2001) 2:67–71. doi: 10.1038/35048073

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Salvioli S, Ardizzoni A, Franceschi C, Cossarizza A. JC-1, but not DiOC6(3) or rhodamine 123, is a reliable fluorescent probe to assess delta psi changes in intact cells: implications for studies on mitochondrial functionality during apoptosis. FEBS Lett. (1997) 411:77–82. doi: 10.1016/s0014–5793(97)00669–8

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Dean PN, Jett JH. Mathematical analysis of DNA distributions derived from flow microfluorometry. J Cell Biol. (1974) 60:523–7. doi: 10.1083/jcb.60.2.523

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Irwin MS, Naranjo A, Zhang FF, Cohn SL, London WB, Gastier-Foster JM, et al. Revised neuroblastoma risk classification system: A report from the children's oncology group. J Clin Oncol. (2021) 39:3229–41. doi: 10.1200/jco.21.00278

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Xie YH, Jiang HZ. The exploration of mitochondrial-related features helps to reveal the prognosis and immunotherapy methods of colorectal cancer. Cancer Rep (Hoboken). (2024) 7:e1914. doi: 10.1002/cnr2.1914

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Shi Y, Huang G, Jiang F, Zhu J, Xu Q, Fang H, et al. Deciphering a mitochondria-related signature to supervise prognosis and immunotherapy in hepatocellular carcinoma. Front Immunol. (2022) 13:1070593. doi: 10.3389/fimmu.2022.1070593

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Brodeur GM, Seeger RC, Schwab M, Varmus HE, Bishop JM. Amplification of N-myc in untreated human neuroblastomas correlates with advanced disease stage. Science. (1984) 224:1121–4. doi: 10.1126/science.6719137

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Huang M, Weiss WA. Neuroblastoma and MYCN. Cold Spring Harb Perspect Med. (2013) 3:a014415. doi: 10.1101/cshperspect.a014415

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Otte J, Dyberg C, Pepich A, Johnsen JI. MYCN function in neuroblastoma development. Front Oncol. (2020) 10:624079. doi: 10.3389/fonc.2020.624079

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Wang SF, Tseng LM, Lee HC. Role of mitochondrial alterations in human cancer progression and cancer immunity. J BioMed Sci. (2023) 30:61. doi: 10.1186/s12929-023-00956-w

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Zhao Y, Liao P, Huang S, Deng T, Tan J, Huang Y, et al. Increased TOX expression associates with exhausted T cells in patients with multiple myeloma. Exp Hematol Oncol. (2022) 11:12. doi: 10.1186/s40164–022-00267–0

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Huang S, Liang C, Zhao Y, Deng T, Tan J, Zha X, et al. Increased TOX expression concurrent with PD-1, Tim-3, and CD244 expression in T cells from patients with acute myeloid leukemia. Cytometry B Clin Cytom. (2022) 102:143–52. doi: 10.1002/cyto.b.22049

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Huang S, Liang C, Zhao Y, Deng T, Tan J, Lu Y, et al. Increased TOX expression concurrent with PD-1, Tim-3, and CD244 in T cells from patients with non-Hodgkin lymphoma. Asia Pac J Clin Oncol. (2022) 18:143–9. doi: 10.1111/ajco.13545

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Chan GC, Chan CM. Anti-GD2 directed immunotherapy for high-risk and metastatic neuroblastoma. Biomolecules. (2022) 12. doi: 10.3390/biom12030358

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Xia Y, Wang C, Li X, Gao M, Hogg HDJ, Tunthanathip T, et al. Development and validation of a novel stemness-related prognostic model for neuroblastoma using integrated machine learning and bioinformatics analyses. Transl Pediatr. (2024) 13:91–109. doi: 10.21037/tp-23–582

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Krystal J, Foster JH. Treatment of high-risk neuroblastoma. Children (Basel). (2023) 10. doi: 10.3390/children10081302

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Zhang D, Kaweme NM, Duan P, Dong Y, Yuan X. Upfront treatment of pediatric high-risk neuroblastoma with chemotherapy, surgery, and radiotherapy combination: the CCCG-NB-2014 protocol. Front Oncol. (2021) 11:745794. doi: 10.3389/fonc.2021.745794

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Park JR, Scott JR, Stewart CF, London WB, Naranjo A, Santana VM, et al. Pilot induction regimen incorporating pharmacokinetically guided topotecan for treatment of newly diagnosed high-risk neuroblastoma: a Children's Oncology Group study. J Clin Oncol. (2011) 29:4351–7. doi: 10.1200/jco.2010.34.3293

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Zheng L, Jia J, Finger LD, Guo Z, Zer C, Shen B. Functional regulation of FEN1 nuclease and its link to cancer. Nucleic Acids Res. (2011) 39:781–94. doi: 10.1093/nar/gkq884

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Dratwa M, Wysoczańska B, Łacina P, Kubik T, Bogunia-Kubik K. TERT-regulation and roles in cancer formation. Front Immunol. (2020) 11:589929. doi: 10.3389/fimmu.2020.589929

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Shi M, Guo H, Bai Y, Niu J, Niu X, Sun K, et al. Upregulated mitosis-associated genes CENPE, CENPF, and DLGAP5 predict poor prognosis and chemotherapy resistance of Acute Myeloid Leukemia. Cancer biomark. (2022) 35:11–25. doi: 10.3233/cbm-203170

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Tang N, Dou X, You X, Shi Q, Ke M, Liu G. Pan-cancer analysis of the oncogenic role of discs large homolog associated protein 5 (DLGAP5) in human tumors. Cancer Cell Int. (2021) 21:457. doi: 10.1186/s12935–021-02155–9

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Garcia SN, Guedes RC, Marques MM. Unlocking the potential of HK2 in cancer metabolism and therapeutics. Curr Med Chem. (2019) 26:7285–322. doi: 10.2174/0929867326666181213092652

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Erdreich-Epstein A, Robison N, Ren X, Zhou H, Xu J, Davidson TB, et al. PID1 (NYGGF4), a new growth-inhibitory gene in embryonal brain tumors and gliomas. Clin Cancer Res. (2014) 20:827–36. doi: 10.1158/1078–0432.Ccr-13–2053

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Yang J, Li S, He J, Xu Q, Xie M, Yang C, et al. Dual role of PID1 in regulating apoptosis induced by distinct anticancer-agents through AKT/Raf-1-dependent pathway in hepatocellular carcinoma. Cell Death Discovery. (2023) 9:139. doi: 10.1038/s41420–023-01405–1

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Muller PA, Vousden KH. p53 mutations in cancer. Nat Cell Biol. (2013) 15:2–8. doi: 10.1038/ncb2641

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Gonfloni S, Caputo V, Iannizzotto V. P63 in health and cancer. Int J Dev Biol. (2015) 59:87–93. doi: 10.1387/ijdb.150045sg

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Jong NN, McKeage MJ. Emerging roles of metal solute carriers in cancer mechanisms and treatment. Biopharm Drug Dispos. (2014) 35:450–62. doi: 10.1002/bdd.1903

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Fernández-Hernando C, Suárez Y, Rayner KJ, Moore KJ. MicroRNAs in lipid metabolism. Curr Opin Lipidol. (2011) 22:86–92. doi: 10.1097/MOL.0b013e3283428d9d

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Fontanella RA, Ghosh P, Pesapane A, Taktaz F, Puocci A, Franzese M, et al. Tirzepatide prevents neurodegeneration through multiple molecular pathways. J Transl Med. (2024) 22:114. doi: 10.1186/s12967-024-04927-z

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Lu Q, Ni Y, Wang W, Wang L, Jiang T, Shang L. Dynamin 3 inhibits the proliferation of non-small-cell lung cancer cells by suppressing c-MET-GBR2-STAT3 complex formation. Front Cell Dev Biol. (2021) 9:641403. doi: 10.3389/fcell.2021.641403

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Inokawa Y, Nomoto S, Hishida M, Hayashi M, Kanda M, Nishikawa Y, et al. Dynamin 3: a new candidate tumor suppressor gene in hepatocellular carcinoma detected by triple combination array analysis. Onco Targets Ther. (2013) 6:1417–24. doi: 10.2147/ott.S51913

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Lin S, Tan L, Luo D, Peng X, Zhu Y, Li H. Linc01278 inhibits the development of papillary thyroid carcinoma by regulating miR-376c-3p/DNM3 axis. Cancer Manag Res. (2019) 11:8557–69. doi: 10.2147/cmar.S217886

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Ma Y, Guan L, Han Y, Zhou Y, Li X, Liu Y, et al. siPRDX2-elevated DNM3 inhibits the proliferation and metastasis of colon cancer cells via AKT signaling pathway. Cancer Manag Res. (2019) 11:5799–811. doi: 10.2147/cmar.S193805

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Peng Z, Wang S, Wen D, Mei Z, Zhang H, Liao S, et al. FEN1 upregulation mediated by SUMO2 via antagonizing proteasomal degradation promotes hepatocellular carcinoma stemness. Transl Oncol. (2024) 44:101916. doi: 10.1016/j.tranon.2024.101916

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Hosfield DJ, Mol CD, Shen B, Tainer JA. Structure of the DNA repair and replication endonuclease and exonuclease FEN-1: coupling DNA and PCNA binding to FEN-1 activity. Cell. (1998) 95:135–46. doi: 10.1016/s0092–8674(00)81789–4

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Ishimi Y, Claude A, Bullock P, Hurwitz J. Complete enzymatic synthesis of DNA containing the SV40 origin of replication. J Biol Chem. (1988) 263:19723–33. doi: 10.1016/S0021-9258(19)77695-9

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Goulian M, Richards SH, Heard CJ, Bigsby BM. Discontinuous DNA synthesis by purified mammalian proteins. J Biol Chem. (1990) 265:18461–71. doi: 10.1016/S0021-9258(17)44775-2

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Storici F, Henneke G, Ferrari E, Gordenin DA, Hübscher U, Resnick MA. The flexible loop of human FEN1 endonuclease is required for flap cleavage during DNA replication and repair. EMBO J. (2002) 21:5930–42. doi: 10.1093/emboj/cdf587

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Wang K, Xie C, Chen D. Flap endonuclease 1 is a promising candidate biomarker in gastric cancer and is involved in cell proliferation and apoptosis. Int J Mol Med. (2014) 33:1268–74. doi: 10.3892/ijmm.2014.1682

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Zhang K, Keymeulen S, Nelson R, Tong TR, Yuan YC, Yun X, et al. Overexpression of flap endonuclease 1 correlates with enhanced proliferation and poor prognosis of non-small-cell lung cancer. Am J Pathol. (2018) 188:242–51. doi: 10.1016/j.ajpath.2017.09.011

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Yuwei X, Bingzi D, Zhaowei S, Yujie F, Wei Z, Kun L, et al. FEN1 promotes cancer progression of cholangiocarcinoma by regulating the Wnt/β-catenin signaling pathway. Dig Liver Dis. (2024) 56:695–704. doi: 10.1016/j.dld.2023.08.050

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Xiao Y, Yin M, Wang Y, Shi R, Mao S, Hua Y, et al. Post-translational modifications by acyl groups regulate FEN1's activities and play essential roles in cell proliferation and DNA repair. Protein Pept Lett. (2023) 30:597–607. doi: 10.2174/0929866530666230529152209

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Wu H, Yan Y, Yuan J, Luo M, Wang Y. miR-4324 inhibits ovarian cancer progression by targeting FEN1. J Ovarian Res. (2022) 15:32. doi: 10.1186/s13048–022-00959–5

PubMed Abstract | CrossRef Full Text | Google Scholar