Introduction

Neuroblastoma (NB) is the most common extracranial malignancy in children and originates from the embryonic neural crest with insidious onset and rapid progression, accounting for 15% of childhood cancer deaths1,2. Clinically, NB can be divided into low-risk (LR), intermediate-risk (IR) and high-risk (HR) types according to the Children's Oncology Group (COG) classification based on the demarcation of age at diagnosis, the International Neuroblastoma Staging System (INSS) stage, the tumor tissue MYCN status, the International Neuroblastoma Pathology Committee (INPC) classification and ploidy. In addition, NB is characterized by a heterogeneous disease spectrum ranging from patients with widespread tumors that spontaneously regress or differentiate without treatment to treatment-resistant tumors with metastatic spread despite intensive multimodal treatment approaches. The 5-year event-free survival rate and 5-year overall survival rate was 91.3% and 97.5% in LR-NB patients, respectively; 85.1% and 96.7% in IR-NB patients, respectively; and 37.7% and 48.9% in HR-NB patients, respectively3,4,5. Therefore, it is urgent to systematically study the occurrence and development of HR-NB, not only for improving the understanding of biological functions, but to also improve therapeutic strategies for HR-NB.

Studies have shown that the lack of HR-NB risk classification diagnostic models and effective therapeutic targets are the main reasons for the significantly lower survival rate than low- and intermediate-risk NB (LIR-NB)6. Alessandra Dondro et al. established a diagnostic model that facilitated the early diagnosis of NB by multiparameter flow cytometry, which can provide new personalized treatments for children with NB and improved survival rate7. By analyzing the expression of nucleolin on the surface of NB cells, Chiara Brignole et al. found that nucleolin was an innovative NB therapeutic cell target and a promising diagnostic model for clinical application was established based on nucleolin8. Through systematic studies, Chiao-Hui Hsieh et al. found that aurora kinase inhibitors interfered with carbohydrate and fatty acid metabolism pathways, resulting in metabolic imbalance. The mitochondrial yellow enzyme acyl-CoA dehydrogenase may be a potential therapeutic target for MYCN-amplified neuroblastoma9. Therefore, systematic research on HR-NB to find its diagnostic biomarkers and abnormal metabolic pathway is expected to improve the survival rate of HR-NB.

It is worth noting that the development of omics has brought new ideas for the diagnosis and treatment of diseases. Metabolomics aims to characterize all small molecules in a sample to accurately reflect the biological metabolic characteristics of disease states, which is beneficial to understanding the pathophysiological processes in disease progression and to find new biomarkers for disease diagnosis and prognosis10. Shivanand Pudakalakatti et al. performed a metabolomics analysis of brain tumor patients and found that platelet-associated lactate, acetate, glutamine, glutamate, succinate, alanine, and pyruvate could be used as biomarkers for brain cancer11. Alessio Imperiale et al. found lower concentrations of glucose, serine, and glycine and increased levels of choline-containing compounds, taurine, lactate, and alanine in small intestinal neuroendocrine tumors by metabolomics analysis, and these metabolites could act as biomarkers, which are of great significance for the development of new targeted therapies in the future12. Transcriptomics uses high-throughput sequencing methods to study all mRNAs transcribed in specific cells, tissues or individuals at a specific time and state from the overall level, which can reveal differences in gene expression and structure in different functional states, and elucidate molecular mechanisms13,14. Transcriptomics analysis of breast tumors by Jun Wang et al. found that in estrogen receptor (+) breast cancer patients, larger body size was associated with upregulation of genes related to the tumor necrosis factor-α/mediated nuclear factor kappa B signaling pathway. In estrogen receptor (−) breast cancer patients, larger body size was additionally associated with downregulation in genes involved in interferon α and interferon γ immune response and phosphatidylinositol 3-kinase/AKT/mammalian target of rapamycin signaling15. Abel Sousa et al. described the differences between sexes in gastric or thyroid normal and tumor tissues in detail by transcriptomics. They found hundreds of sex-biased genes and the peroxisome proliferator-activated receptor signaling pathway was found in normal gastric and thyroid tissues. The abundance of specific sex-biased genes, especially with incidences of overexpression in females, revealed molecular differences and commonalities between sexes, providing new insights into the potential differential risks of these cancers16. In particular, Shancheng Ren et al. discovered metabolic pathway alterations in prostate cancer by combining metabolomics and transcriptomics, and found abnormal expression of cysteine and methionine metabolism, nicotinamide adenine dinucleotide metabolism and hexosamine biosynthesis. In addition, the metabolite sphingosine exhibited high specificity and sensitivity in distinguishing prostate cancer from benign prostatic hyperplasia, promoting the development of new diagnostic biomarkers and therapeutic targets, which will help to distinguish prostate cancer from benign prostatic hyperplasia17. Therefore, the combined analysis of metabolomics and transcriptomics has important application value in HR-NB and with this method, altered metabolic pathways and diagnostic biomarkers in HR-NB can be unearthed in order to establish early diagnosis models and new therapeutic targets for HR-NB.

In this study, a metabolomic analysis of a total of 96 plasma clinical samples and 55 clinical NB tissue samples was conducted, and metabolomics data and transcriptomics data was combined to implement a comprehensive network analysis of NB in order to explore the abnormal pathways of HR-NB and build an early diagnostic model based on candidate biomarkers. The innovation of this work is summarized as follows: (1) By systematic jointing analysis, the transcriptional and metabolic differences between LIR-NB and HR-NB was assessed to uncover the dysregulated network of HR-NB. (2) A noninvasive plasma-based HR-NB risk classification diagnostic model was established using reverse transcription chain reaction (RT-PCR) based approaches. The novel discovery of the HR-NB dysregulation network and NB risk classification diagnostic model is expected to have crucial implications for the robust risk rating of HR-NB and development of a targeted therapy in the future.

Methods and materials

Moral approval

After obtaining approval from the Ethics Committee of Henan Children's Hospital (Approval Number: 2019-H-K11), plasma samples, tissue samples, and relevant clinical data were obtained from patients undergoing surgery at Henan Children's Hospital. Informed consent was obtained from the patients for all samples. All methods were performed in accordance with relevant guidelines and regulations.

Sample collection

A total of 96 plasma (58 cases of HR-NB, 38 cases of LIR-NB) samples and 55 NB tissue samples (32 cases of HR-NB, 23 cases of LIR-NB) were collected and processed at the Henan Children's Hospital from October 2018 to January 2022. Inclusion criteria consisted of: (1) confirmed diagnosis of pathological NB; (2) the risk grade was clinically assessed based on COG classification; (3) informed consent was signed by children or their parents. Exclusion criteria consisted of: (1) complications of other diseases; (2) informed consent was not signed by children or their parents. Plasma samples were obtained from fasting plasma of NB patients on the morning of surgery and immediately frozen in a −80 °C refrigerator for metabolomics analysis. NB tissue samples were obtained from the patient's surgically resected tumor tissue and stored directly in liquid nitrogen for transcriptomic analysis. Tables S1 and S2 showed that there were no significant differences in age and gross tumor volume between HR-NB and LIR-NB whereas the tumor metastasis condition, MYCN amplification, and radiological risk factors were significantly different between HR-NB and LIR-NB samples.

Metabolomics analysis via high performance liquid chromatography-mass spectrometry (HPLC–MS)

The plasma samples were taken from −80 °C storage and immediately thawed in a 4 °C refrigerator. After 10 s of vortexing, 150 μL plasma was added to a 1.5 mL microcentrifuge tube and 450 μL acetonitrile was added at 4 ℃. After 5 min of vortexing at 3000 r/min followed by centrifugation at 13,000 r/min for 15 min (4 ℃), 300 μL supernatant was carefully extracted. Quality control (QC) samples were prepared by mixing equal amounts of supernatant from all samples and were used to evaluate the stability of the overall experimental results. All of the extracts were analyzed using an Agilent 6210 time-of-flight MS system with an Agilent 1100 HPLC, a photodiode array detector, a high-resolution-time-of-flight-MS with an electrospray ionization source and an Agilent workstation. The chromatographic separation was performed on an Agilent Poroshell 120 EC‐C18 (2.7 μm, 3.0 × 100 mm) column. Metabolomics data were collected as follows: mobile phase: A = 0.1% formic acid water, B = 0.1% formic acid acetonitrile, elution conditions: 0–3 min, 5–60% B; 3–25 min, 60–90% B; 25–30 min, 90–100% B; 30–40 min, 100% B. Settings consisted of an injection volume of 10 μL, a column temperature of 30 °C and the flow rate of 0.3 mL/min. Mass spectrometry (MS) negative and positive mode conditions: nitrogen as drying gas, nitrogen temperature 325℃, flow rate 12 L/min, atomization pressure 35 psi; capillary voltage: positive mode 4,000 V, negative mode 3,500 V; fragmentation voltage: positive mode 215 V, negative mode 175 V, separator voltage 60 V; mass acquisition range: all the negative modes were 0.05–1.5 KDa. The samples were analyzed by HPLC–MS to obtain the original data files. Agilent Masshunter HPLC–MS software was used to convert the original data files into a common format. Based on R language platform, XCMS software package was used to conduct retention time (RT) calibration, peak recognition, noise filtering and peak matching for the obtained .mzData format files, and to set the allowable deviation of mass charge ratio and RT (mass/charge ratio tolerance = 0.025DA, RT tolerance = 0.5 min). The metabolites with RT deviation of 0.5 min and mass number deviation of 0.025 Da were considered to be the same metabolites. Finally, the data matrix containing mass/charge ratio, RT, peak area and other information were obtained. Metabolites were identified using both primary and secondary MS techniques. Firstly, the acquired primary MS information was subjected to targeted secondary MS analysis to obtain secondary MS information and provide a reference for subsequent qualitative analysis. Next, according to the precise mass number of excimer ions such as [M+H]+ions and the high-resolution target MS/MS spectrum, combined with the fragmentation laws of various metabolites, and through the online database METLIN (http://metlin.scripps.edu/), HMDB (http://hmdb.ca/) retrieval and literature retrieval methods, analysis and derivation of the possible structure of differential metabolites, get the information of candidate metabolites. The metabolomics analysis, partial least-squares discrimination analysis (PLS-DA), heatmap, volcano map, enrichment analysis, pathway analysis and biomarkers, were conducted by the MetaboAnalyst (https://www.metaboanalyst.ca/MetaboAnalyst/home.xhtml).

Transcriptomics detection through RNA-sequencing (RNA-seq) analysis

Total RNA was extracted using TRIzol reagent according to the manufacturer’s protocol. RNA purity and quantification were evaluated using the NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). RNA integrity was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The libraries were constructed using TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. The transcriptome sequencing and analysis were conducted by OE Biotech Co., Ltd. (Shanghai, China). The libraries were sequenced on an Illumina HiSeq X Ten platform and 150 bp paired-end reads were generated. About 48.349 M raw reads for each sample were generated. Raw data (raw reads) of fastq format were firstly processed using Trimmomatic18 and the low quality reads were removed to obtain clean reads. Around 47.459 M clean reads for each sample were retained for subsequent analyses. The clean reads were mapped to the human genome (GRCh38) using HISAT219. Fragments per kilobase of exon model per million mapped fragments (FPKM)20 of each gene was calculated using Cufflinks21 and the read counts of each gene were obtained by HTSeq-count22. Differential expression analysis was performed using the DESeq (2012) R package23. P value < 0.05 and |log2(fold change)|> 1 were set as the threshold for significantly differential expression. Hierarchical cluster analysis of differentially expressed genes was performed to demonstrate the expression pattern of genes in different groups and samples. Open database sources, including the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG)24,25,26, MetaboAnalyst, Human Metabolome Database and National Center for Biotechnology Information were used to identify metabolic pathways.

Joint analysis of the metabolomics and transcriptomics

Finally, comprehensive transcriptomics and metabolomics analyses were performed using the MetaboAnalyst 5.0. The joint-pathway analysis module was selected to perform topological analysis by entering official gene symbols and compound names with optional fold changes to assess the potential importance of individual molecules (i.e., nodes) according to their position in the network.

Detection the HR-NB plasma potential biomarkers using quantitative PCR based approach

The 24 differential genes were selected as potential candidate biomarkers based on transcriptomics results and literature. Corresponding primers for the differential genes were designed through the website of the National Center for Biotechnology Information (www.ncbi.nlm.nih.gov). Quantitative PCR (qPCR) was performed using the HiScript III All-in-one RT SuperMix kit and AceQ qPCR SYBR Green Master Mix kit (Vazyme, Nanjing) based on kit instructions. The housekeeping gene NAGK was selected as an internal control for mRNA abundance. Fold changes in the levels of target gene mRNA were determined using the formula 2−ΔΔCt.

Develop the risk diagnostic model of NB using qPCR based approach

Logistic regression analysis was used to establish a regression equation for the test set (27 cases of HR-NB, 18 cases of LIR-NB) and validated against the validation set (15 cases of HR-NB, 15 cases of LIR-NB). SPSS software version 21.0 (IBM Corp., Armonk, New York) was used for data processing and GraphPad Prism 8.0 (GraphPad Software Inc., San Diego, USA) was used for mapping.

Results and discussion

The general idea of this research is shown in Fig. 1. Metabolomics analysis was performed on 96 plasma samples (58 cases of HR-NB, 38 cases of LIR-NB) and transcriptomics analysis was performed on 55 NB tissue samples (32 cases of HR-NB, 23 cases of LIR-NB). By integrating metabolomics and transcriptomics data, using PLS-DA, heatmap, enrichment analysis, pathway analysis and other methods for analysis and processing, the abnormal pathway network of HR-NB was comprehensively analyzed and potential clinical therapeutic targets were discovered. Finally, a diagnostic model was established, which is expected to be used for the early diagnosis of HR-NB.

Figure 1
figure 1

Outline of research method. Metabolomics analysis was performed on plasma samples and transcriptomics analysis was performed on NB tissue samples. By integrating metabolomics and transcriptomics data, the abnormal pathway network of HR-NB was analyzed and potential clinical therapeutic targets were discovered. Finally, a diagnostic model was established.

Full size image

The Metabolome differences between HR-NB and LIR-NB

To explore the differences in metabolites among HR-NB, IR-NB and LR-NB, plasma metabolomics analysis was first conducted using a non-targeted metabolomics-based approach. The PCA plot demonstrates the strong clustering of QC samples in both positive and negative modes, providing evidence for the robustness of our findings (Fig. S1). In order to comprehensively understand the metabolomics of NB with different risk levels, the PLS-DA among HR-NB, IR-NB and LR-NB was conducted. As shown in Fig. S2, the IR-NB group and the LR-NB group were clustered together and were situated far away from the HR-NB group, indicating that the LR-NB and IR-NB groups were closely related. However, the plasma metabolites in the IR-NB and LR-NB groups were significantly different from the HR-NB metabolites both in positive mode and negative mode, consistent with the clinical trend of high 5-year survival rates in children with LIR-NB3,4,5. Therefore, in order to better distinguish HR-NB, IR-NB, and LR-NB, IR-NB and LR-NB were merged into one group (LIR-NB). Following, the differences between LIR-NB and HR-NB metabolites were systematically studied.

In order to intuitively express metabolomic differences between LIR-NB and HR-NB, cluster analysis on the plasma metabolites of NB was conducted based on the correlation of compounds and presented in the form of a heatmap (Figure S3), illustrating the differences between the two groups. In order to understand the metabolomics of HR-NB and LIR-NB, the PLS-DA among HR-NB and LIR-NB was conducted first. In positive mode, cumulative R2Y was at 0.722 and Q2 was at 0.575 (Fig. 2A) and in negative mode, cumulative R2Y was at 0.705 and Q2 was at 0.506 (Fig. 2B), indicating that the model had excellent predictive ability. Volcano plots show the metabolites in the HR-NB and LIR-NB groups in positive and negative modes (Fig. 2C,D). Plasma metabolites with fold change > 1.5 and P < 0.05 in the volcano plot were identified as differential metabolites. Therefore, a total of 31 metabolites changed significantly in positive mode, including 21 up-regulated and 10 down-regulated metabolites (Table 1). In negative mode, 13 differential compounds were identified, including 4 up-regulated compounds and 9 down-regulated compounds (Table 2). According to the differential metabolite results, the expression of phospholipids, amino acids and carnitine was significantly altered in HR-NB. Phospholipids play an important role in the composition and determination of biophysical properties of biological membranes, including their mobility, cell signaling, cell–cell interactions in tissues, and molecular transport27,28. Therefore, we believed that phospholipids played an important role in the development of HR-NB due to the marked alteration of phospholipids. Lysophosphatidylcholine (LPC) has been found to be a biomarker for several tumors, including prostate and lung cancer. LPC has also been associated with inflammation, oxidative stress, insulin resistance, apoptosis, lipid remodeling and signal transduction lipogenesis, which is consistent with our metabolomic findings of LPC alterations, thus we suggested that LPC played an important role in HR-NB29,30,31,32,33. Furthermore, in order to further visualize the differential metabolites, heatmaps of differential compounds were drawn according to the correlation of differential compounds (Fig. 2E,F). The heatmap shows the up-regulation of compounds such as LysoPc(0:0/18:0) and the down-regulation of metabolites such as SM(d18:2/14:0), indicating that these compounds have a certain correlation to HR-NB. Hence, a total of 44 differential metabolites were found in metabolomics, demonstrating the significant differences between HR-NB and LIR-NB.

Figure 2
figure 2

Plasma metabolomics between HR-NB and LIR-NB. (A) Comparison of orthographic projections of HR-NB and LIR-NB in PLS-DA 3D plot in positive mode. (B) Comparison of orthographic projections of HR-NB and LIR-NB in PLS-DA 3D plot in negative mode. (C) Volcano plot of metabolite of HR-NB vs. LIR-NB in positive mode. (D) Volcano plot of metabolite of HR-NB vs. LIR-NB in negative mode. (E) The heatmap shows clear distinction of metabolites between HR-NB and LIR-NB patients in positive mode. (F) The heatmap shows clear distinction of metabolites between HR-NB and LIR-NB patients in negative mode. (MetaboAnalyst 5.0, https://www.metaboanalyst.ca/).

Full size image
Table 1 Differential expressed metabolites in HR-NB vs. LIR-NB in positive mode.
Full size table
Table 2 Differential expressed metabolites in HR-NB vs. LIR-NB in negative mode.
Full size table

The altered pathways and biomarkers between HR-NB and LIR-NB

In order to discover abnormal metabolic pathways based on the discovery of abnormal metabolites, enrichment analysis and pathway analysis were conducted. According to the 44 most altered metabolites, the mitochondrial β-oxidation of long chain saturated fatty acids, β-oxidation of very long chain fatty acids, betaine metabolism, carnitine synthesis, oxidation of branched chain fatty acids, plasmalogen synthesis, mitochondrial β-oxidation of short chain saturated fatty acids, methionine metabolism, fatty acid metabolism and glycine and serine metabolism were enriched in the enrichment analysis (Fig. 3A). Recent studies have identified fatty acid oxidation (FAO) as an important source of metabolites that promote cancer growth, showing high activity in various cancers such as breast cancer, glioma, ovarian cancer, which is consistent with the results of our enrichment analysis, FAO is also significantly altered in HR-NB, indicating the important role of FAO in HR-NB34,35,36,37. Interestingly, it has been shown that the activation of FAO can protect cancer cells from effects such as glucose deprivation and hypoxia whereas most cancer cells use glycolysis as the main energy source, therefore we deemed that the glycolysis in HR-NB activated FAO and affected the metabolic pathway of HR-NB. In order to further search for potential abnormal metabolic pathways and visualize the results, pathway analysis was performed. Mitochondrial β-oxidation of long-chain saturated fatty acids, β-oxidation of very long chain fatty acids, plasmalogen synthesis, betaine metabolism, oxidation of branched chain fatty acids and methionine metabolism were significantly changed in the pathway analysis (Fig. 3B). Therefore, 10 significantly altered metabolic pathways were identified by enrichment analysis and 6 significantly altered metabolic pathways were identified by pathway analysis, which are helpful to understand the abnormality of NB pathway network.

Figure 3
figure 3

Altered pathways and biomarkers in metabolomics. (A) Enrichment analysis of differential metabolism reveals various metabolic changes. (B) Pathway analysis reveals significant abnormalities in the pathway. (C) ROC curves and boxplots of metabolic biomarker SM (d16:1/16:0). (D) ROC curves and boxplots of metabolic biomarker SM (d18:2/14:0).

Full size image

To explore the plasma biomarkers of HR-NB in metabolomics and provide a non-invasive strategy for NB risk classification, receiver operating characteristic (ROC) curves of differential metabolites were analyzed. The iconic biomarkers SM (d16:1/16:0) and SM (d18:2/14:0) were found (Fig. 3C,D). Other biomarkers are shown in Fig. S4, Table S3. Area under curve (AUC) of the ROC for all biomarkers were > 0.7, indicating that these metabolites may be potential biomarkers of HR-NB. Both SM (d16:1/16:0) and SM (d18:2/14:0) belong to the sphingolipids, a class of lipids that are expressed in all eukaryotic cells. They are relatively abundant in neuronal tissue cells, enriched in the plasma membrane of cells, and more strongly enriched in certain plasma membrane structural domains, such as the myelin sheath of oligodendrocytes or the parietal membrane of enterocytes38. In addition, they play a signaling role in many cellular functions, such as growth, cell cycle progression, differentiation and apoptosis39,40. Therefore, we inferred that alterations in sphingolipid content were associated with the development of NB and might be a valid biomarker for NB. Hence, 10 altered metabolic pathways were found in the enrichment analysis, 6 altered metabolic pathways were found in the pathway analysis, and 7 biomarkers were found by the biomarker analysis. These findings are helpful to understanding the abnormality of NB pathway network and targeted therapy.

Transcriptomics analysis uncovers the abnormal expression gene between HR-NB and LIR-NB

To further analyze the differences between HR-NB and LIR-NB, transcriptomics differences between 32 HR-NB tissues and 23 LIR-NB tissues were investigated. The detailed results of total RNA concentration, A260/A280, A260/A230, 28S/18S and RNA integrity number of the extracted samples are shown in Table S4. The RNA integrity number values extracted in this study were all greater than 7. The pre-processing results of sequencing data quality showed that the RawBases values of each sample ranged from 6.49 to 7.76 G, the CleanBases values ranged from 6.00 to 7.22 G, and the percentage of Q30 bases in each sample ranged from 92.59 to 95.18%. The GC content of each sample ranged from 44.29 to 50.53% (Table S5). Combining with FPKM (Figure S5) and the total number of mRNAs detected in the samples (Figure S6), the RNA quality of both groups met the standard and could be used for subsequent analysis.

In order to visually express the transcriptomics differences between LIR-NB and HR-NB, the tissue RNA of NB was clustered based on RNA correlation and presented it in the form of a heatmap (Fig. S7), illustrating the obvious differences between the two groups. NB tissue RNAs with P < 0.05 and |log2(fold change)|> 1 in the volcano plot (Fig. 4A) were identified as differential genes. A total of 1,408 differential genes were found, including 1,116 for the number of up-regulated genes and 292 for the number of down-regulated genes (Fig. 4B). The top 20 up-regulated and down-regulated genes in HR-NB and LIR-NB tissues are shown in Tables 3 and 4. It has been suggested that in FABP4-mediated macrophages can promote the proliferation and migratory phenotype of NB cells by inactivating the NF-κB-IL1α pathway through ubiquitination of ATPB, which increases the migration, invasion and tumor growth of NB cells41 and in turn, is consistent with our results where the level of FABP4 increased 14.238-fold. Therefore, we surmised that highly expressed FABP4 may be associated with the risk of NB. In addition, the study by Zhu, Y et al. found that the expression level of UNC5D in the LR-NB group was significantly higher than that in the HR-NB group, suggesting that a high level of UNC5D expression was significantly associated with a good prognosis42, which was consistent with our results where UNC5D expression exhibited a significantly different fold change of 0.281. Therefore, we thought that low expression of UNC5D might be significantly associated with aggressive tumor behavior. The top 100 differential genes were displayed by a cluster analysis heatmap (Fig. 4C), which was used to more intuitively distinguish the differences between HR-NB and LIR-NB, indicating that there were significant differences between the two groups. RT-PCR is a standard method for measuring gene expression due to its high sensitivity and convenience for validating RNA-seq results. Eleven of our selected differential genes were validated by RT-PCR. The RT-PCR results were consistent with transcriptomic results, illustrating the reliability of our transcriptomics results (Fig. 4D). Therefore, 1,408 significantly different genes were found between HR-NB and LIR-NB by transcriptomics, illustrating the differences between the two groups and the results were validated by RT-PCR.

Figure 4
figure 4

Tumor tissue transcriptomics and validation between HR-NB and LIR-NB. (A) The volcano plot shows differentially expressed genes across transcriptomics. Red represents up-regulated mRNAs, green represents down-regulated mRNAs, and gray represents mRNAs with no change in expression. (B) Differentially expressed genes of HR-NB compared to LIR-NB in tumor tissues. (C) The heatmap shows segregation of HR-NB and LIR-NB based on mRNA profile. (oebiotech, https://cloud.oebiotech.com/task/) (D) Expression trends of genes in RT-PCR are consistent with transcriptomics results. 1 to 11 represent: FABP4, CXCL9, MFAP4, CD3, IL-10, CDK2, CIP2A, CHL1, TCF7L2, UNC5D, ERBB3.

Full size image
Table 3 The top 20 genes significantly up-regulated in HR-NB vs. LIR-NB.
Full size table
Table 4 The top 20 genes significantly down-regulated in HR-NB vs. LIR-NB.
Full size table

KEGG and GO analysis between HR-NB and LIR-NB in transcriptomics

To discover abnormal pathways on the basis of differential genes, GO analysis and KEGG analysis were performed. GO annotation analysis of the obtained differential genes was performed to analyze the metabolic pathways in which the genes are located to determine their possible biological functions. The GO enrichment analysis top 30 bar chart shown in Fig. 5A. GO classification bar chart is shown in the Figure S8A. Regarding biological processes, the top three regulated expression were extracellular matrix organization, chemokine-mediated signaling pathway, neutrophil chemotaxis. Regarding cellular component, the top three significantly regulated expressions were extracellular space, extracellular region, extracellular matrix. Regarding molecular function, the top three significantly up-regulated expressions were chemokine activity, extracellular matrix structural constituent, CCR chemokine receptor binding. In addition, KEGG analysis (Figs. 5B, S8B) was performed and cytokine receptor interaction was found to be the most altered pathway, with 69 mRNAs enriched. Cytokines act by binding to specific receptors in the plasma membrane of target cells. The exploration of cytokine receptor interaction is important for understanding the pathogenesis of various human diseases, especially cancers, as well as for identifying potential therapeutic targets. In this study, it was confirmed that cytokine receptor interaction plays an important role in HR-NB. Therefore, more research on cytokine receptor interaction provides a new perspective to understand the molecular mechanism of HR-NB. Relationships between differential genes were visualized by protein–protein interaction circle diagrams (Figure S9), illustrating the close connection between genes such as HBB, HBD, HBA1, HBG2, etc. Therefore, the complex pathways in the transcriptomics between HR-NB and LIR-NB were excavated by GO analysis and KEGG analysis, such as glycan biosynthesis and metabolism and lipid metabolism, which provided experimental support for the subsequent targeted therapy of NB.

Figure 5
figure 5

GO and KEGG analysis between HR-NB and LIR-NB in transcriptomics. (A) GO analysis of biological processes, cellular components and molecular functions of up- and down-regulated genes. (B) Bubble chart of the KEGG pathways (top 20) that differentially expressed genes significantly involved in. The dot color represents the p-value and the dot size represents the number of differential genes.

Full size image

Integrated transcriptomics and metabolomics analyses between HR-NB and LIR-NB

In order to facilitate the systematic study of NB by linking important metabolites and genes through shared metabolic pathways, joint-pathway analysis was used. Nine significantly altered pathways were revealed by integrated analysis of transcriptomics and metabolomics data (Table 5). These altered pathways can be visualized in Fig. 6A, which includes glycerolipid metabolism, retinol metabolism, arginine biosynthesis and linoleic acid metabolism with P-values < 0.05 and impact ≥ 0.5, indicating that dysregulation of these pathways may lead to increased COG grading in NB. The differentially expressed genes associated with these pathways are shown in Table 6.

Table 5 Joint analysis pathways of differential metabolites and genes.
Full size table
Figure 6
figure 6

Integrated transcriptomics and metabolomics analyses of NB metabolic pathways. (A) Joint-pathway analysis of differential genes and differential metabolism. (B) The glycine, serine and threonine metabolism pathway, (C) the glycerolipid metabolism pathway and (D) the glycolysis or gluconeogenesis pathway with altered significantly genes in HR-NB compared to LIR-NB. Significant overexpression in red, significant downexpression in green.

Full size image
Table 6 Related differentially expressed genes by joint-pathway analysis.
Full size table

According to the results of the joint analysis, three pathways were found to be of interest: glycine, serine and threonine metabolism; glycolysis or gluconeogenesis and glycerolipid metabolism, each of which with significant differences in gene expression. In the network diagram of glycine, serine and threonine metabolism (Fig. 6B), the expression of GXT, ALAS2, AOC3, BHMT, CHDH and SDS were up-regulated in HR-NB compared to LIR-NB, which was consistent with metabolomic results (decreased levels of glycine and threonine and marked alterations in amino acid metabolism), but inconsistent in terms of serine (increased levels of serine). Glycine, serine and threonine are important precursors for protein, nucleic acid and lipid synthesis and participate in carbohydrate metabolism pathways. Glycolysis provides 3-phosphoglycerate to promote serine biosynthesis. Serine can be interconverted with glycine through serine hydroxymethyltransferase and excess serine biosynthesis drives tumorigenesis. Therefore we theorized that the activation of HR-NB glycolysis may promote serine synthesis whereas the interconversion of serine and threonine may cause glycine, serine, and threonine results to vary inconsistently43,44,45,46,47. In the glycerolipid metabolism pathway (Fig. 6C), the expressions of AGPAT2, DGAT2, DGKB, GPAM and PNPLA2 were affected. Glycerolipid, an abundant cellular lipid with physiological roles in energy metabolism and membrane structure, is synthesized by modifying the 3-phosphoglycerol backbone through acylation and dephosphorylation reactions48,49,50. GPAT, AGPAT, and DGAT are key genes in this pathway by regulating the glycerol-3-phosphate pathway to control intracellular glycerolipid levels enzymes. Studies have shown that overexpression of GPAT3 leads to increased formation of triacylglycerols, and AGPAT-2 is also overexpressed in various cancers where specific inhibition of AGPAT can induce apoptosis in cancer cells51,52. These studies are consistent with the increased expression of GPAT, AGPAT, and DGAT genes in joint pathway analysis. Therefore, changes in lipid metabolism are closely related to NB and glyceride metabolism may be an important pathway for the development of HR-NB. Finally, the expression of FBP1, HK2, PCK1 was up-regulated and the expression of PKLR was down-regulated in the altered glycolysis or gluconeogenesis pathway (Fig. 6D). In contrast to normal cells relying on mitochondrial oxidative phosphorylation to release energy, most cancer cells use glycolysis as their primary energy source, which produces ATP faster than oxidative phosphorylation to power rapid cell division53. This is consistent with our results, indicating the abnormality of glycolysis in HR-NB. Therefore, these abnormally expressed genes may be used as targets for clinical therapy of NB. Other altered pathways are shown in Figure S10. Nine abnormal pathways were found through combined metabolomics and transcriptomics analysis and three interesting pathways were further mined and speculated to play a role in expanding the scope of targeted therapy.

Classification of NB with selected transcriptome candidate biomarkers

Currently, the risk stratification methods for HR-NB are intricate, clinical presentations are atypical, imaging examinations pose diagnostic challenges, and histological diagnosis remains controversial. HR-NB necessitates comprehensive evaluation involving imaging studies, pathology assessments and other investigations, resulting in the majority of cases being detected at advanced stages. In order to improve the low rate of early diagnosis of HR-NB, a more efficient early diagnostic model was established as a supplement to existing early diagnostic methods. The 24 candidate genes screened according to the results of transcriptomics and literature were analyzed by RT-PCR in order to find biomarkers that could be used for diagnosis (Table S6). The results for individual candidate genes were calculated using Eq. 2−ΔΔCt and the sensitivity and specificity were further calculated. However, the results showed that the areas under the ROC curve of the S100A9, CDK2 and UNC5D were 0.685, 0.564 and 0.809, respectively, indicating that the detection performance of individual biomarkers was limited (Figure S11). Therefore, a diagnostic model combining the three biomarkers was established and the regression equation of the diagnostic model Y = 1.495 − 0.510 X1 (S100A9) − 0.713 X2 (CDK2) + 0.647 X3 (UNC5D) was obtained by logistic regression analysis. In the validation set, the sensitivity and specificity of the diagnostic model were 80% and 80%, respectively, the cutoff value was 0.025 (Fig. 7A) and the AUC of the ROC was 0.836, indicating the effectiveness of the diagnostic model (Fig. 7B). In addition, a close association between these 3 genes and NB was also found. Overexpression of S100A9 enhanced the proliferation, migration and invasion of NB cells54,55,56,57, suggesting that S100A9 may be involved in the development of NB tumors58. In this study, S100A9 increased 7.23-fold in RNA-seq, therefore the high expression of S100A9 may be related to the risk classification of NB. The CDK2 gene encodes a member of the serine/threonine protein kinase family of proteins that can regulate cell cycle progression59,60, which is consistent with our results of altered glycine, serine and threonine metabolism. Zhu et al. found that the expression level of UNC5D mRNA in the LR-NB group was significantly higher than that in the HR-NB group, suggesting that high levels of UNC5D expression were significantly associated with a good prognosis42. Our results are consistent with their study, showing that UNC5D expression was decreased in the HR-NB group and increased in the LIR-NB group. It has also been indicated that low UNC5D expression was significantly associated with aggressive tumor behavior and overexpression of UNC5D significantly inhibited malignant cell behavior, including cell proliferation and migration as well as tumor growth61. Therefore, S100A9, CDK2 and UNC5D are viable biomarkers that can be utilized in combination as an early diagnostic model for predicting the plasma risk classification of NB. This approach holds promise for non-invasive, low-cost detection of NB at an early stage.

Figure 7
figure 7

Establishment of HR-NB early diagnosis model. (A) The prediction accuracies by the S100A9, CDK2 and UNC5D in test phase and validation set are compared between HR-NB and LIR-NB. (B) AUC value for prediction of HR-NB based on 3 plasma genes (S100A9, CDK2, UNC5D), AUC = 0.836.

Full size image

Conclusions

In this study, 96 clinical plasma samples and 55 clinical NB tissue samples were analyzed and 1,408 differential genes and 44 differential metabolites were identified. The metabolomics and transcriptomics characteristics of HR-NB patients were demonstrated and a comprehensive network analysis was carried out. Compared with LIR-NB, HR-NB showed significant differences in glycerolipid metabolism, retinol metabolism, arginine biosynthesis and linoleic acid metabolism. The combination of S100A9, CDK2 and UNC5D was subsequently selected as the risk stratification early diagnostic model for HR-NB. The area under the ROC was 0.837 and the sensitivity and specificity were both 80.0%, indiciated that the diagnostic model could be used for early diagnosis of HR-NB and can be therapeutic targets in the future. In summary, the dysregulated network was examined by joint analysis of metabolomics and transcriptomics and a diagnostic model for HR-NB was developed.