TFAP2D drives neuroblastoma progression: a disulfidptosis-fatty acid metabolism-based molecular subtyping and prognostic model

Introduction

Neuroblastoma (NB) is a malignant tumor arising from neural crest precursor cells of the sympathetic nervous system. It accounts for about 13–15% of all cancer-related deaths among children and adolescents (1-3). This disease is highly heterogeneous in its clinical presentation and prognosis varies greatly among different patients (4). Even with intensive and multimodal therapies (surgery, chemotherapy, radiotherapy, autologous stem cell transplantation), the 5-year survival rate of high-risk patients remains below 50% (5-8), and this is a major obstacle to improve overall cure rate of NB.

The clinical response of NB to immunotherapy remains suboptimal, in contrast to its proven effectiveness against a range of other tumors (9). Currently, although targeted monoclonal antibodies against ganglioside disialoganglioside 2 (GD2), chimeric antigen receptor T-cell therapy (CAR-T) cell therapy, and immune checkpoint inhibitors such as programmed cell death protein 1 (PD-1) and programmed cell death ligand 1 (PD-L1) antibodies have improved the survival of NB patients to some extent, their effects are unstable and much lower compared to the treatment of hematologic tumors and melanoma (10). This limitation stems largely from the highly dynamic and immunosuppressive tumor microenvironment (TME) of NB, where sustained interplay among malignant, immune, and stromal cells critically drives disease progression and metastasis (11,12). Firstly, NB cells show limited capacity to effectively initiate antigen presentation and activate tumor-specific T cells because of their weak antigenicity. Multi-genomic sequencing shows that NB has a low mutation burden, low immunogenicity, and reduced levels of human leukocyte antigen (HLA) on cell membrane. Secondly, abnormal vessels within the tumor impair T cell infiltration into tumors, while Treg cells and other myeloid cells derived from the bone marrow increase in both number and activity within the tumor. Thus, NB is usually described as a “cold tumor” for these immune abnormalities, where most anti-tumor immune cells are unable to exert cytotoxic effects. The key to enhancing the effectiveness of immunotherapy and improving survival in NB lies in transforming the immune microenvironment into a “hot tumor” state through molecular targeting.

Fatty acid metabolism (FAM) is central to antitumor immunity (13,14). Specifically, deficiency in acetyl-CoA carboxylase 1 (ACC1), the rate-limiting enzyme of fatty acid synthesis, compromises T cell proliferation and survival. Conversely, fatty acid oxidation is partly restrained in CD8⁺ T cells, which helps maintain immune function and long-term survival and enables rapid immune responses to antigens (15). Inhibition of cholesterol metabolism pathway can lead to an elevation of intracellular cholesterol in tumor-infiltrating lymphocytes (TILs), leading to facilitated T cell activation and generation of superior immune responses. Moreover, fatty acid accumulation enhances the immunosuppressive function of myeloid-derived suppressor cells and induces resistance to immune checkpoint blockade (16). Fatty acid accumulation also directly regulates the function of dendritic cells, reducing their ability to present antigens and activate T cells within the TME (17). Moreover, FAM is pivotal in this regulatory network for which is also important in NB cell survival and sensitivity to chemotherapy (13,18-20). Currently, there is limited research on the association between FAM and the immune microenvironment in NB. Development and validation of more fatty-acid-metabolism-related prognostic biomarkers require further investigation (21,22).

Unlike classical programmed cell death, such as apoptosis or ferroptosis, disulfidptosis is a novel type of cell death caused by disulfide stress, especially in SLC7A11 high cells starved of glucose (23). When combined with FAM, these processes create a unique metabolic dynamic and imply a more specific vulnerability of TME activation that affects the recruitment of immune cells and immune fitness. A shortage of glucose depletes the cellular nicotinamide adenine dinucleotide phosphate (NADPH) reservoir, setting the stage for the progressive and dysregulated accumulation of molecules linked by disulfide bonds, including disulfide-bonded associations in the actin cytoskeleton. This ultimately results in the collapse of the action network and disulfide-mediated cell death. These data reveal the need for cancer cells to maintain equilibrium between cystine acquisition and glycolytic metabolism. In preclinical models, disulfidptosis in cancer cells can be triggered by glucose inhibition, specifically in those exhibiting elevated SLC7A11 levels, effectively inhibiting tumor growth, and showing no significant toxicity to normal tissues. As for now, signatures for death modes in NB, including autophagy and ferroptosis, have been constructed, but disulfide death has not yet been explored (24-26).

Previous work has provided a comprehensive understanding of the molecular landscape underlying tumor-immune crosstalk. Many studies have identified clinical and molecular prognostic indicators from gene expression or DNA methylation profiles for NB, providing substantial evidence for their high prognostic value compared to traditional clinical features (27-31). However, there has been limited systemic research exploring the molecular characteristics of NB-immune interactions through the use of machine learning algorithms. There is still a lack of effective and widely accepted molecular classification systems to predict prognosis and guide clinical applications (32-35).

This study identifies potential mechanisms underlying NB progression by investigating variations in the expression of FAM genes and disulfidptosis, as well as their relationship with the TME. Furthermore, we use the least absolute shrinkage and selection operator (LASSO) and weighted gene co-expression network analysis (WGCNA) to identify six significantly prognostic genes and construct the disulfidptosis-related fatty acid metabolism (DFAM) score. A novel molecular classification system that shows significant correlations to immune cell infiltration in NB, showed strong prognostic capabilities for NB patients. Additionally, using the weights derived for the screened genes during risk model construction, a critical gene, TFAP2D, is discovered. The malignant transformation potential of TFAP2D is validated in vitro. We present this article in accordance with the MDAR reporting checklist (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-700/rc).

Methods

Study design

NB datasets from three independent sources were obtained. Within these samples, molecular subtypes defined by disulfidptosis-related fatty acids were determined and categorized. We then employed bioinformatic approaches to characterize TME features distinguishing these subtypes. Subsequent systematic analysis evaluated the association between subtype-specific molecular profiles and critical clinical-pathological characteristics, patient outcomes, and established immune checkpoint markers. Finally, DrugBank screening was leveraged to identify potential therapeutic compounds targeting distinct DFAM subtypes. This integrated methodology delineates DFAM heterogeneity, clinical implications, TME interplay, and therapeutic susceptibilities in NB, as summarized in Figure S1.

Data sources

NB phenotypic data and gene expression data were sourced from the NCBI Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/), the ArrayExpress database (https://www.ebi.ac.uk/biostudies/arrayexpress), and The Cancer Genome Atlas (TCGA) (https://www.cancer.gov/ccg/research/genome-sequencing/tcga). After obtaining the raw expression data, we applied background adjustment to the data, followed by quantile normalization. The “Combat” technique was used to merge 3 datasets and remove batch effects. The subsequent analyses comprised a cohort of 971 NB patients, since we removed data from patients who did not have recurrence-free survival (RFS) information. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Consensus clustering analysis of DFAM genes

Gene sets were downloaded from the Molecular Signatures Database (MSigDB) database (http://www.broad.mit.edu/gsea/msigdb/) and earlier publications provided 48 DFAM genes for retrieval (Table S1). Table S2 displays all of these genes’ information. NB patient stratification into molecular subgroups was achieved via consensus clustering of DFAM gene expression, using the “ConsensusClusterPlus” R package based on the following criteria: (I) the optimal number of clusters was determined by evaluating the cumulative distribution function (CDF) curve; (II) no cluster exhibited an excessively small sample size; and (III) post-clustering, intra-cluster correlation increased while inter-cluster correlation decreased. To interrogate the activity of DFAM-related biological processes, we performed gene set variation analysis (GSVA) using the “c2.cp.kegg.v7.2” hallmark gene set from the MSigDB.

Relationship between molecular subtypes and the clinical features and prognosis

We examined the associations between genetic subtypes, clinicopathological features, and prognosis. To determine the clinical relevance of the two consensus clusters, we first characterized patient demographics (including age and gender) and clinicopathological features [International Neuroblastoma Staging System (INSS) stage, Children’s Oncology Group (COG) risk, and MYCN amplification]. We then compared survival rate between the subgroups using Kaplan-Meier analysis implemented with the R packages.

Association between molecular subtypes and TME or immune checkpoints (PD-1/PD-L1) in NB

To comprehensively explore the NB immune landscape, we performed the ESTIMATE algorithm for stromal/immune scoring, combined CIBERSORT to resolve the proportions of 22 immune cell types, and employed single-sample gene set enrichment analysis (ssGSEA) to evaluate the overall immune infiltration level within the TME.

Identification and functional annotation of differential gene expression

Differently expressed genes (DEGs) were identified between subtypes using the package limma (adjusted P value <0.05, |log₂fold-change| >1.5) and WGCNA. Functional enrichment analysis on the overlapping gene set (in order to find associated biological functions) was done using the package clusterProfiler.

The WGCNA R package was used to find prognostic-related modules. Soft-thresholding power of 6 to achieve scale-free networks. The adjacency matrix was converted to a topological overlap matrix (TOM). We used dynamic tree cut (minModuleSize =30, deepSplit =2) to find modules. Modules with a correlation greater than 0.75 were merged. To find the key prognostic module, Pearson correlation was performed. The MEturquoise module showed the strongest correlation and was selected for analysis, based on the screening criteria of gene significance (GS) >0.5 and |module membership (MM)| >0.8.

Construction of the DFAM-related prognostic score

To measure the different DFAM types, a score was computed. DEGs associated with NB RFS were used to perform the univariate Cox regression analysis. Second, an unsupervised clustering method based on the expression of prognostic DFAM genes was used to classify the patients into three subtype groups (subtypes A, B, and C) for further investigation. Ultimately, a 1:1 randomization process was employed to divide all NB patients into training (n=485) and testing (n=485) sets. The training set was then utilized to create the FAM-related prognostic DFAM score. In summary, the “glmnet” R package was employed to perform LASSO-regularized Cox regression, thereby reducing the risk of overfitting in the model derived from prognostic genes associated with DFAM. Following the assessment of each independent variable’s trajectory, a 10-fold cross-validation model was built. Based on this, multivariate Cox analysis was utilized to screen candidate genes and formulate a prognostic DFAM score in the training cohort. A total of 556 patients were trained in the training cohort. The patients was divided into low- and high-risk groups using the median DFAM score as the cutoff. Kaplan-Meier analysis was performed to compare survival between groups. Principal components analysis (PCA) was performed using “ggplot2” R package. Following the same stratification principle, all other cohorts and tests were similarly divided. For each group, we obtained receiver operating characteristic (ROC) curve analysis and Kaplan-Meier survival analysis.

Association between genetic mutations and drug response

Somatic mutations in high and low risk NB groups were analysed using mutation annotation format (MAF) data from TCGA processed using R package Maftools and tumor mutation burden (TMB) of each patient in both groups. The “pRRophetic” package estimates half-maximal inhibitory concentration (IC50) values of standard chemotherapeutic agents used in NB treatment to study differences between treatment effects on patients in both groups.

Establishment and validation of a nomogram scoring system

Each variable is assigned a score that is proportional to its weight, with the total score for each variable being the sum of its individual scores. To evaluate the discrimination of the nomogram, time-dependent ROC curves were generated for 3, 5, and 10 years. Additionally, the predicted survival probabilities were calibrated by comparing them to the observed outcomes at corresponding time points.

Cell culture and gene transfection

NB cell lines SK-N-BE2 and SK-N-AS (Procell, Wuhan, China) were cultured in Dulbecco’s Modified Eagle’s Medium (DMEM)/F12 with 10% fetal bovine serum, ampicillin (100 U/mL) and kanamycin (100 mg/mL) at 37 ℃ in a humidified incubator (5% CO₂) at 37 ℃. The TFAP2D siRNAs were transfected with JetPRIME reagent (Polyplus Transfection, France). Sequences of siRNAs are shown in Table S3.

The quantitative polymerase chain reaction (qPCR) analysis

Total RNA was isolated with TRIzol (Invitrogen, Carlsbad, CA, USA) and cDNA was synthesized using PrimeScript RT Master Mix (Takara, Kusatsu, Japan). qPCR was performed with TB Green^® Premix Ex Taq II (Takara) on a 7500 Real Time PCR System (Applied Biosystems, Foster City, CA, USA). Primer sequences are listed in Table S4. All experiments were performed in triplicate.

Western blot

The proteins of whole cell lysis were transferred onto a polyvinylidene fluoride (PVDF) membrane subsequent to sodium dodecyl sulfate (SDS) gel electrophoresis (Merck Millipore, Billerica, MA, USA). Membranes were first blocked with 5% skim milk solution and incubated with to primary antibodies overnight at 4 ℃, washed and then treated with secondary antibodies for 1 hour at room temperature. Protein signals were visualized using an enhanced chemiluminescence (ECL) kit. Antibodies used in this study are displayed in Table S5.

Cell invasion assay

After 48 hours of transfection/treatment, NB cells 2×10⁴ cells/well) were placed in Matrigel-coated upper chambers. Serum-free medium was used in the upper chamber, whereas 20% FBS medium was used in the lower chamber for invasion. After 24 h, non-invasive cells were removed. The invasive cells in the lower membrane were fixed (4% paraformaldehyde 15 min), stained with 0.1% crystal violet, and visualized under a microscope. Three biological repeat experiments were carried out.

Wound healing assay

NB cells were plated in 6-well plates at a density of 4×10⁵ cells per well and cultured for 48 hours until nearly 100% confluency was reached. A uniform wound was then introduced into the monolayer using a 200 µL pipette tip. After washing with phosphate-buffered saline (PBS) to remove detached cells, serum-free medium was added, and the cells were further incubated for 24 hours. Wound healing was assessed by acquiring images immediately after scratching and at the end of the incubation period. Three biological repeat experiments were carried out.

Statistical analyses

All statistical analyses were conducted with R software (V 4.5.1). P value of less than 0.05 was defined as statistically significant.

Results

Genetic and transcriptional alterations of DFAM genes in NB

The mutation data included 210 NB tumor tissue samples was downloaded from TCGA. 103 (49.05%) samples were found to have multiple mutations with different frequencies in DFAM genes in NB cohort. Among them, genes with the highest frequency of mutation were MUC16 followed by FLG, PTPN11, ABCA13, TTN and FAT2 (Figure 1A). Only 14 (6.67%) samples were detected with mutations in 16/84 prognostic genes (Figure 1B). Next, we mapped the locations for these genes on their respective chromosomes (Figure 1C).

Figure 1 Genetic alterations of DFAM genes in NB. (A) Mutation frequencies of the top 30 DFAM genes in 210 patients from the TCGA_NBL cohort. (B) Mutation frequencies of the DFAM genes with prognostic value in 210 patients from the TCGA_NBL cohort. (C) Locations of hub genes on the 23 chromosomes. DFAM, disulfidptosis-related fatty acid metabolism; NB, neuroblastoma; NBL, neuroblastoma; TCGA, The Cancer Genome Atlas; TMB, tumor mutational burden.

Identification of novel DFAM-related subtypes for NB patients

Figure 2A illustrates the interactive network encompassing various DFAM and regulatory connections, as well as their correlation (table available at https://cdn.amegroups.cn/static/public/tp-2025-aw-700-1.xlsx). To further characterize the expression patterns of DFAM genes in NB, we performed consensus clustering analysis. The analysis indicated that k=2 was the appropriate choice, dividing the cohort into two subtypes (Cluster A and Cluster B) (Figure 2B, Figure S2). PCA demonstrated a clear separation between these two subtypes based on their DFAM gene expression profiles (Figure 2C). The survival analysis revealed a significant in overall survival (OS) between the two clusters with NB patients who were clustered into subtype B having a better prognosis (P<0.001) (Figure 2D). Cluster A preferentially showed a higher probability of MYCN amplification (P<0.05), an advanced International Neuroblastoma Staging System (INSS) stage (P<0.05), and an association with older Age at diagnosis (P<0.05) compared to cluster B (Figure 2E).

Figure 2 Illustration of DFAM subtypes and their association with clinicopathological characteristics. (A) Network representation of DFAM interactions in NB patients. Each connection between DFAMs signifies their interaction, with line thickness denoting the strength of association. Blue and pink lines denote negative and positive correlations, respectively. (B) Heatmap of the consensus matrix displaying the segregation of DFAM samples into two distinct clusters (k=2). (C) PCA results highlighting pronounced differences in transcriptomes between the two identified DFAM subtypes. (D) The survival analysis of the DFAM related subtypes in NB. (E) Comparative analysis distinctions in clinicopathologic features and DFAM related genes expression levels between the two discernible DFAM subtypes in NB patients. DFAM, disulfidptosis-related fatty acid metabolism; EMTAB, ArrayExpress Archive of Functional Genomics Data; NB, neuroblastoma; INSS, International Neuroblastoma Staging System; PCA, principal components analysis; TCGA, The Cancer Genome Atlas.

Characteristics of the TME in two novel NB subtypes

To explore the potential pathway between the two-novel clusters in ArrayExpress Archive of Functional Genomics Data (EMTAB), GSE49710 and TCGA datasets. The GSVA analysis was conducted, and the result showed that cluster A was significantly enriched in cell cycle, RNA degradation, purine metabolism and pyrimidine metabolism pathways (Figure 3A, table available at https://cdn.amegroups.cn/static/public/tp-2025-aw-700-2.xlsx). We further identified DEGs between Cluster A and Cluster B. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses revealed that Cluster A (poor prognosis) was predominantly enriched in DNA replication, cellular senescence and cell cycle-related pathways, which correlates with its association with high proliferation (Figure S3A). We also compared the expression of immune checkpoints between the two clusters. The differential expression analysis of 47 immune checkpoint genes revealed that Cluster B is more immune-active than Cluster A. Tey regulatory genes CD274 (PD-L1), PDCD1 (PD-1), CTLA4, LAG3 and TIGIT upregulated strongly in Cluster B (Figure S3B). Then, we compare the immune landscape of the two subtypes and find Cluster A contains more immune cells [CD8⁺ T cells, natural killer (NK) cells, B cells, monocytic cells, fibroblasts] (Figure 3B). Moreover, we analyzed the relationship between DFAM genes and the TME using immune infiltration estimation methods such as CIBERSORT. We found statistically significant differences in the infiltration density of all 10 subtypes of immune cells, with similarity between NK cells, T cells, and CD8⁺ T cells (Figure 3C). It is notable that cluster B had a higher immuneScore and ESTIMATES score and a higher level of antitumor immune cell infiltration (Figure 3D,3E). These results indicate Cluster B is “Immune-hot”, patients in this group may benefit more from implantable cardioverter-defibrillator (ICD) therapy.

Figure 3 The tumor microenvironments in the two DFAM related subtypes. (A) GSVA depicting the comparative activation (in red) and inhibition (in blue) of distinct biological pathways between the two distinct NB subtypes. (B,C) Evaluation of the abundance of infiltrating immune cell types within the two different subtypes. (D) The TME score including StromalScore, ImmuneScore, ESTIMATEScore for the different DFAM related subtypes. **, P<0.01; ***, P<0.001. DFAM, disulfidptosis-related fatty acid metabolism; EMTAB, ArrayExpress Archive of Functional Genomics Data; GSVA, gene set variation analysis; NB, neuroblastoma; TCGA, The Cancer Genome Atlas; TME, tumor microenvironment.

Novel gene subtypes for NB patients

To investigate the genetic basis of the observed DFAM clusters, we overlapped the DEGs between two different clusters and the most related gene clusters based on the WGCNA for further exploration. First, in order to find the hub node from many ordinary nodes, a scale-free network with soft threshold power 6 was constructed (Figure 4A,4B). The MEturquoise module was the highest correlation with DFAM genes (correlation index =0.77, Figure 4C). GS 0.5 and |MM| 0.8 gave 454 core genes from MEturquoise module (Figure 4D). Then, by intersecting cluster differentially expressed genes with those genes from WGCNA, a total of 19 intersecting genes were obtained (Figure 4E). Based on the 19 differential genes obtained above, we performed Gene ontology and KEGG pathway enrichment analysis (Figure 4F,4G, table available at https://cdn.amegroups.cn/static/public/tp-2025-aw-700-3.xlsx). Then, cluster analysis identified k=3 as the optimal number, leading to the division of the entire cohort into three distinct groups (A, B, and C) (Figure 4H, Figure S4), and the clusters were stratified into three distinct prognostic groups based on survival analysis: Cluster B was associated with the best prognosis, Cluster A with an intermediate prognosis, and Cluster C with a significantly inferior prognosis (Figure 4I). Further characterization and comparison between clusters found that cluster B, according to DEGs, was also associated with a lower probability of MYCN amplification, lower grade of INSS, and lower age (Figure 4J). While, the specific genetic differences among the three groups were shown in Figure 4K with different degrees of inter-group differences.

Figure 4 Construction of the gene subtypes based on the hub genes. (A,B) Establishment of a scale-free network with a soft threshold value of 6. (C) Identification of 6 modules associated with DFAM related genes in NB, which were subsequently merged into 6 distinct modules. (D) Scatterplot analysis of the turquoise module, pinpointing key genes that met the criteria of GS >0.5 and MM >0.8 within the upper right region. (E) The Venn diagram depicting the overlap between the 10 hub genes identified within the turquoise module and cluster differentially expressed genes. (F,G) The Gene Oncology and KEGG pathway analysis for the hub genes. (H) Heatmap of the consensus matrix displaying the segregation-based genes into two distinct clusters (k=3). (I) Kaplan-Meier survival curves depicting the survival analysis for the 3 distinct gene subtypes, assessed using log-rank tests, revealing statistically significant differences (P<0.001). (J) Exploration of associations between clinicopathologic characteristics and the aforementioned gene subtypes. (K) Evaluation of disparities in the expression of hub DFAM related genes within the context of the 3 gene subtypes. *, P<0.05; **, P<0.01; ***, P<0.001. DFAM, disulfidptosis-related fatty acid metabolism; EMTAB, ArrayExpress Archive of Functional Genomics Data; GS, gene significance; INSS, International Neuroblastoma Staging System; KEGG, Kyoto Encyclopedia of Genes and Genomes; ME, module eigengene; MM, module membership; NB, neuroblastoma; TCGA, The Cancer Genome Atlas; WGCNA, weighted gene co-expression network analysis.

Construction and validation of the DFAM prognostic score for NB

We further included 110 prognostic-related genes by LASSO regression analysis. The LASSO regression model, including 6 genes (TFAP2D, BMP7, RET, FAM19A5, DGKB and MAL), with coefficients assigned to six genes are presented in Figure 5A-5C. We calculated the DFAM risk score using the following formula: = (0.142*TFAP2D + 0.115*BMP7 + 0.105*RET +0.096*FAM19A − 0.139*DGKB − 0.139*MAL). Subsequent analysis of prognostic characteristics of the constructed risk score signature showed that patients with low risk could achieve distinct extended survival (Figure 5D). Based on the median of risk score, patients were classified into high- and low-risk groups. High-risk patients had poorer survival (Figure 5E). Survival distributions of each patient are shown in a scatter plot, which showed a negative correlation between risk score and survival rate, showing that a higher score was associated with a worse clinical outcome (Figure 5F, Figure S5A-S5F). The expression of the 6 genes is shown in the heatmap (Figure 5G, Figure S5G,S5H). Area under the curve (AUC) shows that this signature can properly predict the 1-, 3-, and 5-year (Figure 5H, Figure S5I,S5J). Multivariate Cox regression analysis showed that DFAM risk score was an independent prognostic factor of NB (Figure 5I, Figure S5K). Patient distribution based on survival status among the different DFAM subtypes, gene subtypes, and DFAM risk score categories (Figure 5J). Significant differences were observed across DFAM gene subtypes; subtype B exhibited minimal PRG_score values, whereas subtype C displayed maximal levels (Figure 5K). Notably, subtype A demonstrated significantly higher DFAM_scores than subtype B (Figure 5L).

Figure 5 Development of the DFAM score within NB dataset. (A,B) The LASSO regression analysis and partial likelihood deviance assessment of prognostic genes. (C) The coefficient for the genes explored for risk score model construction. (D) Kaplan-Meier survival analysis for the DFAM related risk score model. (E,F) Visualization of the ranked dot and scatter plots displaying the distribution of DFAM scores based on the corresponding patient survival status. (G) The heatmap for the expression of DFAM related gene in risk score model. (H) ROC curves employed for the prediction of sensitivity and specificity regarding 1-, 3-, and 5-year survival based on the DFAM score. (I) The multiCox regression analysis for the DFAM score model. (J) An alluvial diagram illustrating the distribution of molecular subtypes across various DFAM scores and their associated survival outcomes. (K) Assessment of the disparities in DFAM scores among distinct gene subtypes. (L) Evaluation of the discrepancies in DFAM scores across different subtypes. AUC, area under the curve; DFAM, disulfidptosis-related fatty acid metabolism; INSS, International Neuroblastoma Staging System; NB, neuroblastoma; ROC, receiver operating characteristic.

Clinicopathological analysis based on the DFAM score

To elucidate the association between the DFAM score and clinical parameters in the TCGA_NBL dataset, including age, gender, differentiation grade, MYCN status, mitosis karyorrhexis index (MKI) and COG risk. Notably, we found an elevated DFAM score was significantly associated with patients classified within undifferentiated group (P<0.001), age >18 months (P<0.001), stage 4 group (P<0.001) (Figure 6A-6D). The DFAM score was significantly elevated in patients presenting with high MKI, unfavorable histology, MYCN amplification, and classification within the high-risk COG group (Figure 6E-6H). We further evaluated the prognostic value of the DFAM score by examining its relationship with clinical outcomes in the GSE49710 dataset and E-MTAB-8248 also showed the same conclusion (Figure S6A-S6H and Figure S7A,S7B), patients with higher DFAM risk score were positive related with the tumor progression (P<0.001), and the rearrangement TERT status (P=0.001) (Figure S7C,S7D).

Figure 6 Clinicopathological characteristic analysis of the DFAM score for NB patients. (A-H) The influence of the DFAM score on clinical outcomes in the TCGA_NBL dataset to assess the association between the DFAM score and various clinical parameters, including age at diagnosis, gender, differentiation grade, INSS stage, histology, MYCN status, MKI and COG risk states. COG, Children’s Oncology Group; DFAM, disulfidptosis-related fatty acid metabolism; INSS, International Neuroblastoma Staging System; MKI, mitosis karyorrhexis index; NBL, neuroblastoma; TCGA, The Cancer Genome Atlas.

Immunological microenvironment and immunotherapy for NB based on risk score

To assess the biological significance of the risk stratification defined by the risk score model, gene set enrichment analysis (GSEA) and KEGG pathway enrichment analysis were performed to show that high risk related genes mainly enriched the cell cycle, DNA replication, and the glycine, serine and threonine metabolism pathway (Figure 7A), while the low risk group mainly enriched in axon guidance, cytokine-cytokine receptor interaction and T cell receptor signaling pathway (Figure 7B), these mainly contribute to the development of immune-activated microenvironment. To assess the relationship between the risk score and immune cell infiltration, we employed the CIBERSORT algorithm to quantify the proportions of different immune cell types (table available at https://cdn.amegroups.cn/static/public/tp-2025-aw-700-4.xlsx). An inverse association was noted between risk score and elevated immune scores, increased stromal scores and higher ESTIMATEScores (Figure 7C). To investigate the interrelation between the immune cells, the xCell, Microenvironment Cell Populations-counter (MCPcounter), Estimation of Proportions of Immune and Cancer cells (EPIC), and CIBERSORT algorithm were performed to reveal significant correlations for most immune cell types with these genes, including CD4⁺ T cells, CD8⁺ T cells, NK cells, NK T cells, M1 macrophages, and dendritic cells (Figure 7D). Besides, the scatter plot analyses revealed a direct positive correlation of Fatty Acid Risk score with specific infiltrating immune cells, including regulatory T cells, cytotoxic CD8⁺ T cells, CD4⁺ memory T cells, resting NK cells, and antibody-secreting plasma cells (Figure 7E-7I).

Figure 7 Immunological microenvironment and immunotherapy for NB based on risk score. (A,B) The GSEA analysis for Gene Ontology and KEGG pathway between high and low DFAM score. (C) The associations between DFAM score and tumor immune scores. (D) The associations between the prevalence of immune cell populations and the DFMA risk score. (E-I) The scatter plot analyses revealed a direct positive correlation of fatty acid risk score with specific immune cells including Tregs, CD8 T cells, CD4 memory cells, NK resting cells, and plasma cells. DFAM, disulfidptosis-related fatty acid metabolism; EPIC, Estimation of Proportions of Immune and Cancer cells; GSEA, gene set enrichment analysis; KEGG, Kyoto Encyclopedia of Genes and Genomes; MCPcounter, Microenvironment Cell Populations-counter; NB, neuroblastoma; NK, natural killer; Tregs, regular T cells; TME, tumor microenvironment.

To further explore the potential of immunotherapy, the Tumor Immune Dysfunction and Exclusion (TIDE) analysis was performed and revealed that patients in the low-risk group had a significantly lower score than those in the high-risk group (P<0.001; Figure 8A-8C). This observation intimates a potential enhanced responsiveness to immunotherapeutic interventions in low-risk patients, attributable to a decreased propensity for immune escape mechanisms. The TIDE analysis showed that almost 75% high-risk patients showed a poor response to immunotherapy. However, only 54% low risk NB patients were non-responder (P=0.001, Figure 8D), and the AUC for the prediction was 0.656 [95% confidence interval (CI): 0.588–0.722, Figure 8E]. Finally, to explore the TME within the risk group, the expression profiles of activating immune checkpoints, including ICOS, TNFRSF9, CD80, CD86, CD28, CD40, TNFRSF4, TNFSF14, CD27, TNFRSF25, CD40LG (CD154), TNFSF18, TNFRSF18, TNFSF4 (OX40L), and CD200 (Figure 8F). We found that almost all of these immune checkpoints showed low expression in high DFAM score group, indicating a notable distinction in their immunological landscape. This result showed that the FAM score model could predict whether the NB patients were “cold” or “hot” immune microenvironment, and the immunotherapy response for NB patients.

Figure 8 The potential of immunotherapy response for NB patients based on the DFAM score. (A-C) The TIDE analysis based on TCGA_NBL, GSE49710 and E-MTAB-8248 dataset was performed. (D) The distribution ratio of immunotherapy response between high and low risk groups. (E) ROC analysis of the sensitivity and specificity of immune therapy responses. (F) The differential expression of immune checkpoint marker between high- and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001. AUC, area under the curve; CI, confidence interval; DFAM, disulfidptosis-related fatty acid metabolism; EMTAB, ArrayExpress Archive of Functional Genomics Data; NB, neuroblastoma; NBL, neuroblastoma; TIDE, Tumor Immune Dysfunction and Exclusion; TCGA, The Cancer Genome Atlas; Tregs, regular T cells.

Construction of a prognostic nomogram

To further improve the utility of clinical decision-making of DFAM score, we have constructed a prognostic nomogram. The nomogram integrates the DFAM score with various clinicopathological factors (the MYCN amplification status, the INSS stage, and the age of diagnosis) to estimate 1-, 3-, and 5-years survival rates (Figure 9A). Then, the robustness of the nomogram was further affirmed by calibration plots, which indicated a congruence in performance with an ideal model across both the whole NB dataset (Figure 9B). The nomogram demonstrated superior accuracy for 5-year survival rate predictions (AUC =0.796, Figure 9C). After assessing the predictive accuracy and statistical independence of our model, our primary focus shifted towards investigating the potential clinical utility when integrating this model with conventional clinical parameters routinely employed in the management of NB patients. To assess the clinical utility of our model, decision curve analysis (DCA) was performed (Figure 9D). Our findings demonstrate that the nomogram and DFAM score model, in conjunction with MYCN status, age, and INSS stage, enhances prognostic predictions for NB patients. This augmentation leads to a more favorable clinical outcome for these individuals. What’s more, in order to further assess the accuracy and effectiveness of our model predictions, we conducted a comparative analysis with other 3 models [Gupta’s (36), Tian’s (37) and Yang’s (38)]. The results revealed that our model exhibits superior specificity and sensitivity (Figure 10A-10H). The C-index for DFAM score was 0.823, higher than the other predictive signatures (Figure 10I). Subsequently, we compared the treatment responsiveness between the low- and high-risk patient groups. to explore the common chemotherapeutic agents. It found that individuals in the high DFAM risk group exhibited a decreased IC50 value for vinblastine, vorinostat, mitomycin C, cisplatin, doxorubicin, gemcitabine, axitinib, AKT inhibitor VIII, and etoposide (Figure 11A-11I). Patients with low DFAM scores showed increased sensitivity to lapatinib, gefitinib, and imatinib, as evidenced by their significantly lower IC50 values. In concert, these findings underscore the association between DFAM score and drug sensitivity (Figure 11J-11L).

Figure 9 The prognostic nomogram for NB patient outcomes. (A) Nomogram depicting the prognostic model for estimating the probabilities of 1-, 3-, and 5-year survival in NB patients within all cohorts. (B) The calibration plots depicting the performance of the nomogram in predicting 1-, 3-, and 5-year overall survival rates. (C) ROC curve demonstrated superior accuracy of the nomogram in predicting survival rate based on the DFAM score and the NB patient’s INSS stage, age at diagnosis, and MYCN status. (D) DCA to quantify the DFAM score benefit achievable by NB patients. AUC, area under the curve; DCA, decision curve analysis; DFAM, disulfidptosis-related fatty acid metabolism; INSS, International Neuroblastoma Staging System; NB, neuroblastoma; OS, overall survival.

Figure 10 The accuracy and effectiveness of score models for NB patients. (A-C,G) The ROC curves evaluated the accuracy of the nomogram in predicting 1-, 3-, and 5-year overall survival rates. (D-F,H) The Kaplan-Meier survival analysis for all the risk score model mentioned above. (I) The C-index for DFAM score and other risk score models for NB patients. AUC, area under the curve; DFAM, disulfidptosis-related fatty acid metabolism; NB, neuroblastoma; ROC, receiver operating characteristic.

Figure 11 The sensitivity of chemotherapeutic agents in different risk score group. (A-I) A decreased IC₅₀ for vinblastine, vorinostat, mitomycin C, cisplatin, doxorubicin, gemcitabine, axitinib, AKT inhibitor VIII, and etoposide in high DFAM risk score group. (J-L) An elevated IC₅₀ values for lapatinib, gefitinib, and imatinib with high DFAM score. DFAM, disulfidptosis-related fatty acid metabolism; IC₅₀, half-maximal inhibitory concentration; NB, neuroblastoma.

Knockdown of TFAP2D impairs the malignant biological functions in vitro

To explore the potential biological functions of TFAP2D, the Ctrl siRNA and TFAP2D siRNA#1 and #2 were transfected into SK-N-AS and SK-N-BE2 cells for 48 h, and RT-qPCR and Western blot analyses were performed to confirm that the expression of TFAP2D was downregulated by the siRNA segments at both mRNA and protein level (Figure 12A,12B). Then, the Cell Counting Kit-8 (CCK-8) was performed, and the results demonstrated that transfection with TFAP2D siRNAs (#1 and #2) markedly inhibited cell survival relative to the control (P<0.001, Figure 12C). And the transfection of cells with TFAP2D siRNAs, a significant inhibition in the scratch wound closure was observed, indicating a wound healing process. Quantitative analysis revealed that the rate of wound healing in the TFAP2D siRNA#1 and #2 was significantly lower compared to the Ctrl siRNA (P<0.001, Figure 12D). Finally, the transwell assay was performed, and the results showed that the invasion rate in TFAP2D siRNA#1 and #2 was significantly lower than that in Ctrl group (P<0.01, Figure 12E). These results suggested that down regulation of TFAP2D could inhibit cell proliferation, cell migration and invasion.

Figure 12 TFAP2D regulates the malignant phenotype of NB cells in vitro. (A) TFAP2D mRNA level in BE2 and AS cells transfected with TFAP2D-ctrl, TFAP2D siRNA#1 or TFAP2D siRNA#2. GAPDH acted as the internal control of qPCR. (B) TFAP2D protein level in BE2 and AS cells transfected with TFAP2D-ctrl siRNA, TFAP2D siRNA#1 or TFAP2D siRNA#2. Actin acted as the internal control of Western blot. (C) The proliferation ability was significantly inhibited in BE2 and AS cells transfected with TFAP2D siRNA#1 or TFAP2D siRNA#2, compared with the same cells transfected with control siRNAs. (D,E) The motility and invasion ability of BE2 and AS cells transfected with TFAP2D siRNA#1 or TFAP2D siRNA#2 were significantly inhibited. At least three fields were randomly collected. Scale bars, 100 µm. **, P<0.01; ***, P<0.001. GAPDH, glyceraldehyde-3-phosphate dehydrogenase; OD, optical density; qPCR, quantitative polymerase chain reaction.

Discussion

NB exhibits typical “cold tumor” (39) characteristic due to low mutation burden, loss of antigen presentation molecule expression, and blockade of T cell infiltration into the tumor. Therefore, improvement of immunotherapy in NB has been challenged for years. Previous studies demonstrated that small molecule drugs can reprogram the immune microenvironment by promoting the co-localization of Th1 cells infiltrating the tumor with tumor-associated endothelial cells (TECs), transforming NB from an immune-suppressed “cold tumor” to an immune-activated “hot tumor”.

In our previous study, the combination therapy of anlotinib and PD-1 blockade has shown synergistic effects in the treatment of NB. In terms of mechanism, anlotinib promoted the normalization of the vasulation and recruiting more active immune cells to transform the “Cold” immune microenvironment to a “Hot” status (40). The unique feature of the disulfidptosis-FAM axis in our model is that it acts as a metabolic sensor for TME remodeling. Our study revealed a global alteration in the transcription of disulfidptosis and FAM in NB, and identified two distinct molecular subtypes. We show high risk NB patients exhibit metabolic profiles where excessive fat turnover and disulfide stress suppress T-cell growth and promote immunosuppressive states. Unlike broad metabolic checkpoints, targeting this specific link of disulfidptosis and lipid remodeling can overcome resistance to therapy and increase immunogenicity of the NB environment. Based on the expression of six signature genes, we identified two distinct molecular subtypes in NB, which exhibited significantly different TME characteristics. Compared to the immunoactive “hot tumor” B subtype, the A subtype exhibits an immunosuppressive “cold tumor” state, a lower DFAM score, and is also associated with more advanced clinical and pathological features, as well as poorer OS. Furthermore, we identified six gene subtypes based on the DEGs and WGCNA and construct DFAM score, which is able to predict survival risk of NB patients. There are significant differences observed between low- and high-risk patients in terms of clinical and pathological features, prognosis, TME, and drug sensitivity. Finally, by combining DFAM score with tumor stage and other factors, we have developed a nomogram that can be used for prognostic stratification of NB patients. This model is highly practical, more predictively accurate then other similar models, and further enhances the performance of DFAM score.

Effector T cells are key mediators of antitumor immunity. Among them, γδ T cells are capable of directly recognizing and lysing cancer cells. while indirectly inhibiting tumor progression through various mechanisms. γδT in TME can dissolve cancer cells through the perforin-granzyme pathway, and also exert indirect anti-tumor effects through ligand and antibody-dependent cell-mediated cytotoxicity, secretion of cytokines, including interferon-gamma (IFN-γ) and tumor necrosis factor-alpha (TNF-α), or interaction with other immune cells. Tumor-infiltrating T cell density was significantly elevated in cancerous tissues compared to normal tissues, which is associated with better prognosis. Our analysis established a correlation between the DFAM score and immune infiltration, with the low-risk group exhibiting increased infiltration of activated memory CD4⁺ and CD8⁺ T cells, along with γδ T cells. Consistently, differential gene analysis in this group revealed an enrichment of pathways related to T cell receptor and cytokine-receptor interactions.

It has been observed that B cell infiltration in NB tumors is limited. However, more B cells may exist outside the tumor (41,42). Increased infiltration of immune cells has been associated with clinical outcomes in NB patients (43,44). A study (45) demonstrated the correlation of several immune cell types with RFS and OS, highlighting the previously overlooked role of B cell abundance in NB, which is highly correlated with survival and a hot TME. The presence of organized lymphoid structures rich in B cells, possibly originating from tumor-associated tertiary lymphoid structures at the tumor margin, is closely correlated with both prolonged survival and an effective anti-tumor immune response (46,47). Antigen-specific interactions between B cells and T cells promote the cytotoxicity of CD8⁺ T cells in TME. Cytokines secreted by B cells, including IFN-γ and interleukin (IL)-12, further activate anti-tumor CD8⁺ T and NK cells. We observed significantly lower infiltration of naïve B cells in low-risk group and high DFAM score group. Increased infiltration of B cells promotes the reshaping of the NB immune microenvironment, promotes the recruitment of additional immune cells and activates immune-related pathways, thereby enhancing the effectiveness of NB immunotherapy and improving prognosis.

Tumor-associated macrophages (TAMs), including M1 and M2 (48). The low-risk group, characterized by a higher density of M1 macrophages (which are antitumorigenic and produce type I pro-inflammatory cytokines), may therefore be more likely to benefit from immunotherapy (49). On the other hand, M2 macrophages exhibit immunosuppressive effects and promote tumor growth by facilitating matrix remodeling. High levels of M2 macrophages and a restricted gene signature were found to correlate with a worse prognosis (50). In high-risk group, NB cells promoting macrophages predominate, showing high expression of M2 markers (51). We observed an increased infiltration of M1 macrophages in low-risk group, which correlated with a favorable prognosis. Conversely, increased infiltration of M2 macrophages in high-risk group was associated with a poorer prognosis.

The high-risk group showed a poor predicted response rate of immunotherapy of 25% compared to 46% in low-risk group, but showed higher response to traditional chemotherapy drugs such as vorinostat, vinblastine, cisplatin, and doxorubicin. It is suggested a positive association with high PD-L1 expression and a favorable prognosis in high-risk NB patients. PD-1⁺/PD-L1⁺ macrophages in high-risk tumors may comprise two subgroups, with those displaying a PD-1/PD-L1 expression ratio >1 showing a favorable response to PD-1 blockade therapy (52). In this study, we observed elevated co-expression of PD-1 and PD-L1 in the low-risk group, suggesting a potentially enhanced response to immunotherapy.

TFAP2D, a key molecule that holds the highest weight in DFAM score, is one of TFAP2 family, which is mainly associated with developmental and reproductive processes. As for immune system TFAP2D particularly impacts the immune cells and TAM polarization (53-55). TFAP2D mutations are most likely to occur in lung and gastric tumors, which can interact with cancer and drug targets. Overexpression of TFAP2D is associated with prostate cancer progression and contributes to genomic instability (56). TFAP2D not only influences genome organization and chromatin structure through its downstream targets but also emerges as the most promising therapeutic candidate within the AP-2 family for cancer treatment (57,58). Our study confirms that TFAP2D can promote cell proliferation and enhance metastasis, which demonstrates its malignant transformation potential.

Indeed, further analysis and exploration of specific targets related to immunotherapy are needed in order to develop optimal treatment strategies for more clinical benefit based on risk stratification. Furthermore, all analyses in this study were based on data from public databases, and all samples were retrospectively obtained. Although TFAP2D is strongly associated with metastatic progression, we did not conduct in vivo animal experiments to observe tumor growth and metastasis in living cells. Future studies using NB xenograft model and lung metastases may further understand why TFAP2D promotes NB aggressiveness. More prospective research and additional experiments are needed to confirm our findings.

Conclusions

This study constructed a DFAM-related score and identified new molecular subtypes, which emphasized the crucial impact of DFAM in exploring the TME and regulating treatment response in NB. The model provides a novel perspective for prognostic assessment and significantly contributes to the evolution of immunotherapeutic strategies.

Acknowledgments

We appreciate the NCBI GEO, ArrayExpress and TCGA databases for providing expression data of NB samples, and the clinical-pathologic characteristics.

Footnote

Reporting Checklist: The authors have completed the MDAR reporting checklist. Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-700/rc

Data Sharing Statement: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-700/dss

Peer Review File: Available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-700/prf

Funding: This work was supported by the Scientific Research Project of Tianjin Municipal Education Commission (No. 2024KJ175), Tianjin Health Research Project (No. TJWJ2025QN028), the Tianjin Metrology Science and Technology Project (No. 2025TJMT009), and Tianjin Key Medical Discipline Construction Project (No. TJYXZDXK-3-003A).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tp.amegroups.com/article/view/10.21037/tp-2025-aw-700/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Lee AC. Neuroblastoma: the challenge remains. Singapore Med J 2012;53:1-2.
2. Katta SS, Nagati V, Paturi ASV, et al. Neuroblastoma: Emerging trends in pathogenesis, diagnosis, and therapeutic targets. J Control Release 2023;357:444-459. [Crossref] [PubMed]
3. Xie W, Zhang Y, Xu J, et al. Characteristics, treatments, and outcomes of adolescents and adults with neuroblastoma: a retrospective study in China. Ther Adv Med Oncol 2025;17:17588359251337494. [Crossref] [PubMed]
4. Izbicki T, Mazur J, Izbicka E. Epidemiology and etiology of neuroblastoma: an overview. Anticancer Res 2003;23:755-60.
5. Seidinger AL, Fortes FP, Mastellaro MJ, et al. Occurrence of Neuroblastoma among TP53 p.R337H Carriers. PLoS One 2015;10:e0140356. [Crossref] [PubMed]
6. Chung C, Boterberg T, Lucas J, et al. Neuroblastoma. Pediatr Blood Cancer 2021;68:e28473. [Crossref] [PubMed]
7. Maris JM. Recent advances in neuroblastoma. N Engl J Med 2010;362:2202-11. [Crossref] [PubMed]
8. Maggi E, Landolina N, Munari E, et al. T cells in the microenvironment of solid pediatric tumors: the case of neuroblastoma. Front Immunol 2025;16:1544137. [Crossref] [PubMed]
9. Anderson J, Majzner RG, Sondel PM. Immunotherapy of Neuroblastoma: Facts and Hopes. Clin Cancer Res 2022;28:3196-3206. [Crossref] [PubMed]
10. Neviani P, Wise PM, Murtadha M, et al. Natural Killer-Derived Exosomal miR-186 Inhibits Neuroblastoma Growth and Immune Escape Mechanisms. Cancer Res 2019;79:1151-64. [Crossref] [PubMed]
11. Joshi S. Targeting the Tumor Microenvironment in Neuroblastoma: Recent Advances and Future Directions. Cancers (Basel) 2020;12:2057. [Crossref] [PubMed]
12. Pinto NR, Applebaum MA, Volchenboum SL, et al. Advances in Risk Classification and Treatment Strategies for Neuroblastoma. J Clin Oncol 2015;33:3008-17. [Crossref] [PubMed]
13. Tao L, Mohammad MA, Milazzo G, et al. MYCN-driven fatty acid uptake is a metabolic vulnerability in neuroblastoma. Nat Commun 2022;13:3728. [Crossref] [PubMed]
14. Feng X, Zhang J, Yu B, et al. Lipid metabolic reprogramming in tumor-associated macrophages: A key driver of functional polarization and tumor immunomodulation. Crit Rev Oncol Hematol 2025;215:104881. [Crossref] [PubMed]
15. Xu S, Chaudhary O, Rodríguez-Morales P, et al. Uptake of oxidized lipids by the scavenger receptor CD36 promotes lipid peroxidation and dysfunction in CD8(+) T cells in tumors. Immunity 2021;54:1561-1577.e7. [Crossref] [PubMed]
16. Adeshakin AO, Liu W, Adeshakin FO, et al. Regulation of ROS in myeloid-derived suppressor cells through targeting fatty acid transport protein 2 enhanced anti-PD-L1 tumor immunotherapy. Cell Immunol 2021;362:104286. [Crossref] [PubMed]
17. Lian X, Yang K, Li R, et al. Immunometabolic rewiring in tumorigenesis and anti-tumor immunotherapy. Mol Cancer 2022;21:27. [Crossref] [PubMed]
18. Kazantzis M, Stahl A. Fatty acid transport proteins, implications in physiology and disease. Biochim Biophys Acta 2012;1821:852-7. [Crossref] [PubMed]
19. Talapatra J, Reddy MM. Lipid Metabolic Reprogramming in Embryonal Neoplasms with MYCN Amplification. Cancers (Basel) 2023;15:2144. [Crossref] [PubMed]
20. Agostini M, Melino G, Habeb B, et al. Targeting lipid metabolism in cancer: neuroblastoma. Cancer Metastasis Rev 2022;41:255-60. [Crossref] [PubMed]
21. Ding Y, Yang J, Ma Y, et al. MYCN and PRC1 cooperatively repress docosahexaenoic acid synthesis in neuroblastoma via ELOVL2. J Exp Clin Cancer Res 2019;38:498. [Crossref] [PubMed]
22. Farrell EK, Chen Y, Barazanji M, et al. Primary fatty acid amide metabolism: conversion of fatty acids and an ethanolamine in N18TG2 and SCP cells. J Lipid Res 2012;53:247-56. [Crossref] [PubMed]
23. Liu X, Nie L, Zhang Y, et al. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis. Nat Cell Biol 2023;25:404-14. [Crossref] [PubMed]
24. Hu J, Song F, Kang W, et al. Integrative analysis of multi-omics data for discovery of ferroptosis-related gene signature predicting immune activity in neuroblastoma. Front Pharmacol 2023;14:1162563. [Crossref] [PubMed]
25. Kang W, Hu J, Zhao Q, et al. Identification of an Autophagy-Related Risk Signature Correlates With Immunophenotype and Predicts Immune Checkpoint Blockade Efficacy of Neuroblastoma. Front Cell Dev Biol 2021;9:731380. [Crossref] [PubMed]
26. Liu X, Zhuang L, Gan B. Disulfidptosis: disulfide stress-induced cell death. Trends Cell Biol 2024;34:327-37. [Crossref] [PubMed]
27. Bao X, Shi R, Zhang K, et al. Immune Landscape of Invasive Ductal Carcinoma Tumor Microenvironment Identifies a Prognostic and Immunotherapeutically Relevant Gene Signature. Front Oncol 2019;9:903. [Crossref] [PubMed]
28. Deng X, Lin D, Chen B, et al. Development and Validation of an IDH1-Associated Immune Prognostic Signature for Diffuse Lower-Grade Glioma. Front Oncol 2019;9:1310. [Crossref] [PubMed]
29. Li B, Cui Y, Diehn M, et al. Development and Validation of an Individualized Immune Prognostic Signature in Early-Stage Nonsquamous Non-Small Cell Lung Cancer. JAMA Oncol 2017;3:1529-37. [Crossref] [PubMed]
30. Long J, Wang A, Bai Y, et al. Development and validation of a TP53-associated immune prognostic model for hepatocellular carcinoma. EBioMedicine 2019;42:363-74. [Crossref] [PubMed]
31. Jin W, Zhang Y, Liu Z, et al. Exploration of the molecular characteristics of the tumor-immune interaction and the development of an individualized immune prognostic signature for neuroblastoma. J Cell Physiol 2021;236:294-308. [Crossref] [PubMed]
32. Zafar A, Wang W, Liu G, et al. Molecular targeting therapies for neuroblastoma: Progress and challenges. Med Res Rev 2021;41:961-1021. [Crossref] [PubMed]
33. Hu JJ, Wang HM. Coexistence of orbital rhabdomyosarcoma and adrenal neuroblastic tumor in a child. Pediatr Int 2021;63:859-61. [Crossref] [PubMed]
34. Hu JJ, Wang HM. Primary renal ganglioneuroblastoma mimicking Wilms' tumor in a 3-year-old girl. Pediatr Int 2021;63:106-7. [Crossref] [PubMed]
35. Qiu B, Matthay KK. Advancing therapy for neuroblastoma. Nat Rev Clin Oncol 2022;19:515-33. [Crossref] [PubMed]
36. Gupta M, Kannappan S, Jain M, et al. Development and validation of a 21-gene prognostic signature in neuroblastoma. Sci Rep 2023;13:12526. [Crossref] [PubMed]
37. Tian X, Cao F, Li X, et al. Tumor dormancy is closely related to prognosis prediction and tumor immunity in neuroblastoma. Transl Pediatr 2023;12:445-61. [Crossref] [PubMed]
38. Yang J, Han L, Sha Y, et al. A novel ganglioside-related risk signature can reveal the distinct immune landscape of neuroblastoma and predict the immunotherapeutic response. Front Immunol 2022;13:1061814. [Crossref] [PubMed]
39. Rohila D, Park IH, Pham TV, et al. Targeting macrophage Syk enhances responses to immune checkpoint blockade and radiotherapy in high-risk neuroblastoma. Front Immunol 2023;14:1148317. [Crossref] [PubMed]
40. Su Y, Luo B, Lu Y, et al. Anlotinib Induces a T Cell-Inflamed Tumor Microenvironment by Facilitating Vessel Normalization and Enhances the Efficacy of PD-1 Checkpoint Blockade in Neuroblastoma. Clin Cancer Res 2022;28:793-809. [Crossref] [PubMed]
41. Slyper M, Porter CBM, Ashenberg O, et al. A single-cell and single-nucleus RNA-Seq toolbox for fresh and frozen human tumors. Nat Med 2020;26:792-802. [Crossref] [PubMed]
42. Coughlin CM, Fleming MD, Carroll RG, et al. Immunosurveillance and survivin-specific T-cell immunity in children with high-risk neuroblastoma. J Clin Oncol 2006;24:5725-34. [Crossref] [PubMed]
43. Helmink BA, Reddy SM, Gao J, et al. B cells and tertiary lymphoid structures promote immunotherapy response. Nature 2020;577:549-55. [Crossref] [PubMed]
44. Cabrita R, Lauss M, Sanna A, et al. Tertiary lymphoid structures improve immunotherapy and survival in melanoma. Nature 2020;577:561-5.
45. Schaafsma E, Jiang C, Cheng C. B cell infiltration is highly associated with prognosis and an immune-infiltrated tumor microenvironment in neuroblastoma. J Cancer Metastasis Treat 2021;7: [Crossref] [PubMed]
46. Sharonov GV, Serebrovskaya EO, Yuzhakova DV, et al. B cells, plasma cells and antibody repertoires in the tumour microenvironment. Nat Rev Immunol 2020;20:294-307. [Crossref] [PubMed]
47. Coronella JA, Spier C, Welch M, et al. Antigen-driven oligoclonal expansion of tumor-infiltrating B cells in infiltrating ductal carcinoma of the breast. J Immunol 2002;169:1829-36. [Crossref] [PubMed]
48. Orecchioni M, Ghosheh Y, Pramod AB, et al. Macrophage Polarization: Different Gene Signatures in M1(LPS+) vs. Classically and M2(LPS-) vs. Alternatively Activated Macrophages. Front Immunol 2019;10:1084. [Crossref] [PubMed]
49. Kovaleva OV, Rashidova MA, Sinyov VV, et al. M1 macrophages - unexpected contribution to tumor progression. Front Immunol 2025;16:1638102. [Crossref] [PubMed]
50. Asgharzadeh S, Pique-Regi R, Sposto R, et al. Prognostic significance of gene expression profiles of metastatic neuroblastomas lacking MYCN gene amplification. J Natl Cancer Inst 2006;98:1193-203. [Crossref] [PubMed]
51. Larsson K, Kock A, Idborg H, et al. COX/mPGES-1/PGE2 pathway depicts an inflammatory-dependent high-risk neuroblastoma subset. Proc Natl Acad Sci U S A 2015;112:8070-5. [Crossref] [PubMed]
52. Tang XX, Shimada H, Ikegaki N. Macrophage-mediated anti-tumor immunity against high-risk neuroblastoma. Genes Immun 2022;23:129-40. [Crossref] [PubMed]
53. Yuan F, Sun Y, Dai GC, et al. Comprehensive Analysis of Prognostic Value and Immune Infiltration of TFAP2 Family Members in Bladder Cancer from Database and FFPE Sample. J Cancer 2023;14:3050-65. [Crossref] [PubMed]
54. Li X, Monckton EA, Godbout R. Ectopic expression of transcription factor AP-2δ in developing retina: effect on PSA-NCAM and axon routing. J Neurochem 2014;129:72-84. [Crossref] [PubMed]
55. Van Otterloo E, Li W, Garnett A, et al. Novel Tfap2-mediated control of soxE expression facilitated the evolutionary emergence of the neural crest. Development 2012;139:720-30. [Crossref] [PubMed]
56. Fraune C, Harms L, Büscheck F, et al. Upregulation of the transcription factor TFAP2D is associated with aggressive tumor phenotype in prostate cancer lacking the TMPRSS2:ERG fusion. Mol Med 2020;26:24. [Crossref] [PubMed]
57. Kołat D, Zhao LY, Kciuk M, et al. AP-2δ Is the Most Relevant Target of AP-2 Family-Focused Cancer Therapy and Affects Genome Organization. Cells 2022;11:4124. [Crossref] [PubMed]
58. Jin C, Luo Y, Liang Z, et al. Crucial role of the transcription factors family activator protein 2 in cancer: current clue and views. J Transl Med 2023;21:371. [Crossref] [PubMed]