Comprehensive Analysis and Validation Reveal Potential MYCN Regulatory Biomarkers Associated with Prognosis of Neuroblastoma

Neuroblastoma (NB) is an embryonic malignant tumor that occurs in the sympathetic nervous system. The treatment effect of patients in the high-risk group is poor, and relapse and treatment failure can occur even with multiple combination treatments. MYCN is an independent prognostic factor for NB, and the proportion of MYCN amplication in tumor tissues of high-risk patients reaches 40-50%. Hence, exploring new MYCN targeting molecules is a meaningful way to treat high-risk NB patients. The microarray datasets were obtained from Gene Expression Omnibus, and differentially expressed genes were identied. Gene Ontology, Kyoto Encyclopedia of Genes and Genomes, and miRPathDB were used for function and pathway analysis. The STRING and Cytoscape were used to construct a protein-protein interaction network and modular analysis. The miRNet and NetworkAnalyst databases were used to predict and construct target gene-miRNA and target gene-TF networks. The R2 database was used for expression, correlation, and prognosis analysis. The diagnostic effect of the biomarkers was predicted by ROC analysis, and rt-PCR was used to verify the hub genes. Finally, used specic siRNA and plasmid to knockdown and overexpressed MYCN, the mRNA levels were veried. identied preliminarily conrmed be Our results are expected to become novel biomarkers for the treatment of high-risk NB patients.


Introduction
Neuroblastoma (NB) is the most common extracranial solid malignant tumor in children [1]. It predisposes to the adrenal medulla but can also occur in the sympathetic nervous system of the neck, chest, and abdomen [2,3]. NB usually occurs in younger children, with a median age of 17 months at diagnosis [4][5][6].
Its heterogeneity characterizes NB, and tumors can metastasize without treatment and develop resistance to therapy [7]. Its clinical manifestations range from asymptomatic masses to primary tumors caused by local invasion or widespread diseases [8]. According to reports, the incidence of NB in childhood cancer is second only to leukemia and brain tumors, accounting for about 8% of childhood malignancies and 15% of all childhood cancer deaths [4].
Studies have shown that the prognosis of NB patients is related to factors such as age at diagnosis, tumor stage, and oncogene MYCN ampli cation [9,10]. According to the clinical manifestations and the biological characteristics of the tumor, NB patients can be divided into a low-risk group, a medium-risk group, and a high-risk group [6]. Among them, the age of onset of patients in the low and moderate risk groups is generally within 1.5 years of age, and the lesions are mostly localized. Generally, there is less ampli cation of MYCN, so the response to surgery and chemotherapy is better [11]. In some patients with low risk and good biological characteristics, self-healing can often be achieved, and the long-term survival rate can reach more than 90% [12]. The age of onset of patients in the high-risk group is mostly over 1.5 years old, and they generally have manifestations of MYCN expansion [6]. Extensive metastatic lesions often appear at the time of diagnosis, even in the case of intensive chemotherapy combined with surgery, radiotherapy, and autologous bone marrow stem cell transplantation. The long-term survival rate is still less than 40% [10].
MYCN is a member of the MYC proto-oncogene family, and the prognosis of tumor patients is poor when overexpressed [9]. MYCN can promote tumor progression by inducing cell proliferation, dedifferentiation, and mitochondrial metabolism disorders [13]. According to reports, about 20% of the tumor tissues of NB patients have MYCN ampli cation, and the proportion of MYCN ampli cation in tumor tissues of highrisk patients is as high as 40-50%. MYCN is an independent prognostic factor for NB patients [14,15].
Studies have shown that molecular markers were related to the occurrence and development of tumors and can be used for early tumor screening. However, the problem is that the expression of molecular markers in tumors does not have reasonable speci city [16]. Hence, it is necessary to explore new MYCN targeting molecules further to provide new strategies for diagnosing and treating high-risk NB patients.
This study downloaded the GEO database's microarray datasets GSE12460, GSE13136, GSE14880, GSE16237, and GSE73517. The differentially expressed genes (DEGs) were screened by comparing the gene expression between MYCN ampli ed and MYCN non-ampli ed NB cells. STRING and Cytoscape were used to construct protein-protein interaction (PPI) networks and perform modular analysis of DEGs. Then, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and miRPathDB were used to perform function and pathway analysis, miRNet and NetworkAnalyst were used to construct a target gene-miRNA regulatory and target gene-TF regulatory network. The R2 database was used to con rm the relationship between the genes and MYCN Correlation, and the importance of prognosis has determined FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 as hub genes. From the microarray dataset GSE73517 to verify the results were reliable, the diagnostic effect of the hub genes was predicted by Receiver Operator Characteristic Curve (ROC) analysis. Finally, the use of cell transfection technology to knockdown or overexpression MYCN through PCR experiments preliminarily con rmed that MYCN regulates the expression of hub genes FBXO9, RNF19B, and TRIM36. In conclusion, this study attempted to explore promising new targeting molecules for MYCN that predicts the prognosis of NB. Targeting FBXO9, RNF19B, and TRIM36 in NB cells may represent a novel strategy for NB therapeutic intervention.

Identi cation of differentially expressed genes
Based on the Affy package in the R language (http://cran.r-project.org/), the robust multi-array average (RMA) algorithm converts the original microarray dataset into expression values and then performs background correction and quintiles Normalization and probe summary. The paired t-test of the limma package of the R language was used to analyze the DEGs between the MYCN-ampli cation and the MYCN non-ampli cation NB cells. The volcanos were generated using the ggplot2 package of the R language, and Benjamini-Hochberg's method was used to control the False Discovery Rate (FDR). P-value < 0.05 and | log2FC | > 1 were the critical values for screening DEGs.

GO and KEGG enrichment analysis of DEGs
Use the online tool Venny (v 2.1, https://bioinfogp.cnb.csic.es/tools/venny/) to construct a Venn diagram to identify overlapping modules and genes in the DEGs of the four microarray datasets. DAVID (http://david.abcc.ncifcrf.gov/) is an annotation, visualization, and comprehensive discovery database. It is an online tool for gene function classi cation used to analyze and evaluate biological gene functions [22]. miRPathDB (v 2.0, https://mpd.bioinf.uni-sb.de) is a new dictionary on miRNAs and target pathways, and the database contains many miRNAs, different miRNA target sets, and a wide range of functional biochemical categories [23]. Here, we used the DAVID database to perform GO and KEGG enrichment analysis to study the function of DEGs, and miRPathDB was used for pathway veri cation of miRNAs. P < 0.05 was regarded as the critical point with statistical signi cance.
PPI network construction and module analysis STRING (v 11.5, https://string-db.org) is an online tool to search for interacting genes and proteins. DEGs were entered into the database to construct a PPI network that shows physical and functional interactions. In this study, proteins with a score (minimum required interaction score) > 0.9 were selected for PPI network construction. In addition, Cytoscape (v 3.8.2, https://cytoscape.org) was used to visualize the network, calculate the node degree through Network Analyzer, and draw the PPI network with different colors and sizes to show the adjustment and node degree. The MCODE plug-in (degree cutoff = 2, node score cutoff = 0:2, and K − core = 2) was used to identify the hub modules in the network, and the BiNGO plug-in and DAVID were used to analyze the function and pathway enrichment of the genes in the hub modules, which can better explain the reliability of the results. P < 0.05 was used as the critical value for screening. miRNA prediction of target genes miRNet (v 2.0, https://www.mirnet.ca/) database was used to predict the miRNA associated with the target genes, construct a target gene-miRNA interaction network, and Cytoscape was used to visualize it to characterize further the relationship between the target genes and miRNAs Correlation.
Transcription factor prediction of target genes NetworkAnalyst (v 3.0, https://www.networkanalyst.ca/) is a series of network-based tools for statistical meta-analysis, visual data mining, and integration through the rapid generation of biological networks. It supports meta-analysis of gene lists, and data integration is achieved through robust statistical procedures, followed by visual inspection in the PPI network [24]. Used NetworkAnalyst to predict the Transcription Factor (TF) of the target genes to characterize the correlation between the potential target genes and TFs.

Analysis of Hub Genes
The R2 database (http://r2platform.com) is a genomics analysis and visualization platform that provides a biologist-friendly interface for high-throughput data. It was developed in the Cancer Genomics The cells were digested with ordinary trypsin, and the cell suspension containing 1.5 × 10 5 cells/ml was inoculated in a 6-well plate (2 ml/well), and the cell density during transfection was 50%-60%. Washed the cells twice with PBS, added RPMI 1640 medium without double antibodies (penicillin and streptomycin), small interfering RNAs (MYCN siRNA, siRNA control) (Tongyong, Anhui, China), and MYCN overexpression plasmid (NM_001293228, pDsRed2) -N1 vector, XhoI/KpnI digestion, negative control plasmid) (Genechem, Shanghai, China) was used to transfect NB cells. The target sequences of siRNAs were as follows: MYCN siRNA CGGAGATGCTGCTTGAGAA, used jetPRIME transfection reagent for transfection, refer to the instructions for operation. After 24 hours of transfection, the cells were collected to extract RNA, and RT-PCR detected the mRNA level.

Real-time quantitative PCR
The total RNA was extracted with the kit reagent (Invitrogen), the RNA was 5μg, and the GoScript™ reverse transcription system kit (Promega, UK) was used. The SYBR premix Ex Taq™ II (TaKaRa Clontech) was used following the manufacturer's instructions, and the cDNA after reverse transcription was used for qPCR. The 2(−ΔΔCt) method was used for gene expression analysis. The primers were used to detect the target genes (Table 1). All primers were purchased from Sangon Biotech (Shanghai, China). The PCR reaction program was as follows: 95 ℃ 30 s, 95 ℃ 5 s, 60 ℃ 34 s, a total of 40 cycles. Use the ΔΔCt method for analysis, the target genes mRNA = 2-ΔΔCt.

Statistical Analysis
Use GraphPad 9 software for statistical analysis. The data were expressed as mean ± standard deviation (x̄ ± SD). The t-test method was used to compare the data between the two groups. P < 0.05 was considered statistically signi cant.

Microarray datasets identi cation of DEGs
To identify DEGs, we performed background correction and normalization processing on the NB microarray datasets GSE12460, GSE13136, GSE14880, and GSE16237. Filtered DEGs through the limma package in the R language, with P-value < 0.05 and | log2FC | > 1 as the ltering conditions. Among them, the GSE12460 obtained a total of 1479 DEGs, including 402 up-regulated and 1077 down-regulated DEGs, the GSE13136 obtained a total of 2778 DEGs, including 1223 up-regulated and 1555 downregulated DEGs, the GSE14880 obtained a total of 2068 DEGs, including 535 up-regulated and 1533 down-regulated DEGs, GSE16237 obtained a total of 3842 DEGs, including 496 up-regulated and 3346 down-regulated DEGs. The volcanos were used to show the differential expression of genes in the microarray datasets (Fig 1A).

Identi cation of overlapping DEGs
To identify the overlapping DEGs, we used the tool Venny to analyze the four microarray datasets and took the overlapping part of the four microarray datasets. In this study, we identi ed a total of 431 overlapping DEGs, including 90 up-regulated DEGs and 341 down-regulated DEGs (Fig 1B).

GO and KEGG enrichment of overlapping DEGs
To further understand the functions and pathways involved in overlapping DEGs, we used the DAVID database for enrichment analysis. The results showed that in Biological Processes (BP), DEGs were signi cantly enriched in cell adhesion, peripheral nervous system development, extracellular matrix organization, muscle organ development, and cell proliferation, in Cellular Components (CC), DEGs Signi cantly enriched in the neuromuscular junction, neuron projection, cell surface, lopodium, and axon, regarding Molecular Function (MF), DEGs were signi cant in collagen binding, glycosaminoglycan binding, structural constituent of muscle, receptor binding, actin-binding Enrichment. In addition, KEGG analysis showed that DEGs were signi cantly enriched in Proteoglycans in cancer, cell adhesion molecules, MicroRNAs in cancer, Glycine, serine, and threonine metabolism, and Axon guidance ( Fig 1C).

PPI network construction and module analysis
To establish the protein-protein interaction, we used the STRING database to construct a PPI network for overlapping DEGs, involving a total of 194 nodes and 401 edges (Fig 2A). The non-interacting proteins were ltered out by Cytoscape, and the algorithm analysis of MCODE was performed, and the module with the highest MCODE score = 9.0 was selected as the hub module ( Fig 2B). Enrichment analysis through BiNGO plug-in showed that DEGs in the module were mainly enriched in ligase activity, zinc ion binding, acid-amino acid ligase activity, ligase activity, forming carbon-nitrogen bonds, transition metalion binding, male germline stem cell division, asymmetric cell division, leg morphogenesis, and ubiquitin ligase complex. The visualization results are shown in Fig 3A. In addition, KEGG enrichment analysis showed that these DEGs were related to Acute myeloid leukemia, Transcriptional misregulation in cancer, MicroRNAs in cancer, and Pathways in cancer. The visualization results are shown in Fig 3B. miRNA prediction of target genes To predict the miRNA related to the target genes, we used the miRNet database to perform predictive analysis, imported the prediction results into Cytoscape for visualization, and constructed a target gene-miRNA interaction network (Fig 4). The results showed that FBXO9 interacts with 42 miRNAs (for example, hsa-mir-329-3p), FBXO17 interacts with 40 miRNAs (for example, hsa-mir-124-3p), HECW2 interacts with 43 miRNAs (for example, hsa-mir-19a-3p), MIB2 interacts with 17 miRNAs (for example, hsa-mir-335-5p), RNF19B interacts with 117 miRNAs (for example, hsa-mir-17-5p), RNF213 interacts with 140 miRNAs (for example, hsa-mir-101-3p), TRIM36 interacts with 40 miRNAs (for example, hsa-mir-25-3p), TRIM71 interacts with 163 miRNAs (for example, hsa-let-7a-5p), ZBTB16 interacts with 53 miRNAs (for example, hsa-mir-15a-5p). The miRNAs in the network were further ranked by degree method, and the top ve miRNAs with the highest scores are shown in Table 2. In addition, we further validated the pathways in which these miRNAs were involved through the miRPathDB database, with P-value < 0.05 as the ltering conditions. The results showed that these miRNAs were signi cantly annotated to MicroRNAs in cancer, mTOR signaling, Proteoglycans in cancer, Pathways in cancer, and TGF-beta Signaling Pathway. Most of the pathways, including the ve pathways listed earlier, were supported by experimental validation (especially strong experimental validation). The 20 signi cant pathways are shown in Table 3.

Transcription factor prediction of target genes
To predict the TFs related to the target genes, we used the NetworkAnalyst database to conduct predictive analysis and constructed the target gene-TF interaction network (Fig 5A). The results showed that FBXO9 interacts with 38 TFs (for example, KLF16), FBXO17 interacts with 13 TFs (for example, ZEB1), MIB2 interacts with 35 TFs (for example, SIX5), RNF19B interacts with 17 TFs (for example, EGR2), RNF213 interacts with 51 TFs (for example, IRF1), TRIM36 interacts with 10 TFs (for example, CBX8), TRIM71 interacts with 19 TFs (for example, SOX5), ZBTB16 interacts with 4 TFs (for example, WT1). Further GO enrichment analysis of TFs in the network showed that BP was mainly enriched in peptidyl-lysine modi cation, covalent chromatin modi cation, and histone modi cation, CC was signi cantly enriched in transcription regulator complex, RNA polymerase II transcription regulator complex, and transcription repressor complex. MF was mainly involved in DNA-binding transcription activator activity, DNA-binding transcription activator activity, RNA polymerase II-speci c and DNAbinding transcription repressor activity, RNA polymerase II-speci c. The visualization results are shown in Fig 5B. In addition, KEGG enrichment analysis showed that these TFs were related to Transcriptional misregulation in cancers, TGF-beta signaling pathway, Thyroid hormone signaling pathway. The visualization results are shown in Fig 5C.

Survival analysis of hub genes
To study the correlation between gene expression and the prognosis of NB patients, we further analyzed the data with a sample size of 649 through the R2 database. Finally, we screened out FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 among these eight DEGs. The genes closely related to the patient's prognosis were used as hub genes. Compared with NB patients with relatively low hub genes expression in tumor tissues, NB patients with relatively high hub genes expression have a higher overall survival rate (P < 0.001) (Fig 6). We further analyzed the correlation between hub genes expression and MYCN expression in NB tissues. The results showed that the expression of hub genes and MYCN in NB tumor tissues were negatively correlated (P < 0.001) (Fig 7A). In addition, the expression of hub genes in MYCN non-ampli ed tumor tissues was higher than that in MYCN-ampli ed tumor tissues (P < 0.001) (Fig 7B).

Diagnostic Effectiveness of Biomarkers for Neuroblastoma
To make the results more reliable, we used the microarray dataset GSE73517 to verify. Expression levels for MYCN, FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 are shown in the heatmap ( Fig  8A). Consistent with the above results, FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 have signi cantly higher expression in MYCN-nonampli ed NB cells (P < 0.01) (Fig 8B). We analyzed the correlation between seven effective biomarkers and MYCN. Correlation results are shown in Fig 8C. In addition, the microarray dataset GSE73517 was used to verify the diagnostic effectiveness of NB biomarkers through ROC analysis. AUC more than 0.800 was considered to have excellent speci city and sensitivity to diagnose NB. As shown in Fig 8D,   To study whether MYCN regulates the hub genes, we designed MYCN siRNA and overexpression plasmids. MYCN siRNA was transfected into MYCN ampli ed BE2 cells, and MYCN overexpression plasmid was transfected into MYCN non-ampli ed AS cells. Then detected the changes in the expression of hub genes at the mRNA levels. The results showed that after MYCN siRNA was transfected into BE2 cells, the expression of FBXO9, RNF19B, and TRIM36 increased at the mRNA levels, which were consistent with the predicted results ( Fig 10A). After the MYCN overexpression plasmid was transfected into AS cells, the expression of FBXO9, RNF19B, and TRIM36 decreased at the mRNA levels ( Fig 10B). These data indicated that MYCN has a counter-regulatory effect on the expression of FBXO9, RNF19B, and TRIM36 in NB cells.

Discussion
NB is the most common extracranial tumor in childhood, which has different clinical manifestations and courses according to the biological characteristics of the tumor[6]. Studies have shown that maintaining proliferation is a fundamental feature of cancer cells [25,26]. In many types of cancer, this ability can be acquired by the MYC family to promote the growth or abnormal activation of tumor transcription factors [27,28]. MYC family member MYCN ampli cation is a cancer-causing event in high-risk NB development, leading to poor prognosis in NB patients [29]. In addition, abnormal MYCN activation is also associated with the pathogenesis of neuroendocrine cancer in adults, including glioblastoma, small cell lung cancer, and pancreatic cancer [30]. Hence, the MYCN status determined in NB diagnosis is a signi cant adverse prognostic factor. Studying the role and mechanism of new targeting molecules of MYCN is a critical way to explore the treatment of high-risk NB patients.
Bioinformatics analysis technology has been increasingly used to nd new therapeutic targets and diagnostic markers for various cancers in recent years [31][32][33]. Compared with a single cohort study, a multi-cohort study has a lower false-positive rate and false-negative rate [34]. In this study, to reduce the impact of batch effects and increase the credibility of DEGs identi cation, we selected four microarray datasets from the same platform. Here, we have identi ed a total of 431 overlapping DEGs, including 90 up-regulated DEGs and 341 down-regulated DEGs. The enrichment analysis at the GO level shows that these DEGs were mainly involved in cell adhesion, extracellular matrix organization, cell proliferation, and receptor binding. In addition, KEGG analysis showed that DEGs were signi cantly enriched in Proteoglycans in cancer, cell adhesion molecules, and MicroRNAs in cancer. These rich functions and pathways provided insights into the molecular mechanism of NB pathogenesis.
To identify the hub genes in DEGs, we constructed a PPI network, analyzed by the MCODE plug-in in Cytoscape, and selected the highest-rated module in the network as the hub module (including FBXO9, FBXO17, HECW2, MIB2, RNF19B, RNF213, TRIM36, TRIM71, and ZBTB16). Based on the hub module, we have carried out an enrichment analysis of genes. The results suggested that these DEGs were closely related to many signaling pathways that regulate the development of cancers (For example, Pathways in cancer, mTOR signaling, and PI3K-AKT Signaling in Cancer), which also provided potential value for further research on hub modules.
In the above study, MicroRNAs in cancer and Transcription misregulation in cancer were signi cant enrichment results for hub modules. Studies have shown that miRNAs and TFs play an essential regulatory role in developing various cancers, including childhood cancer NB [35][36][37]. In addition, miRNAs and TFs were reported to play an indispensable role in regulating NB proliferation, apoptosis, invasion, and metastasis [38][39][40][41][42]. Here, we further predicted and analyzed the genes in the hub modules of miRNAs and TFs. In the miRNA prediction analysis, we found FBXO9, FBXO17, HECW2, MIB2, RNF19B, RNF213, TRIM36, TRIM71, and ZBTB16 interacted with a total of 615 miRNAs. In addition, we further veri ed through the miRPathDB database that these miRNAs were involved in regulating MicroRNAs in cancer, mTOR signaling, Proteoglycans in cancer, Pathways in cancer, and TGF-beta Signaling Pathway. Using the NetworkAnalyst database for TFs predictive analysis, we found FBXO9, FBXO17, MIB2, RNF19B, RNF213, TRIM36, TRIM71, and ZBTB16 participated with a total of 187 TFs regulation. Enrichment analysis from TFs in the network shows that these TFs were closely related to Transcriptional misregulation in cancers, TGF-beta signaling pathways, and Thyroid hormone signaling pathway. From enrichment analysis of DEGs and hub modules to miRNAs and TFs enrichment veri cation, we found that these biomarkers were closely related to various cancer signaling pathways.
In combination with the R2 database, we further evaluated the prognosis of genes expression in NB patients. The hub genes of FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 were identi ed as closely related to the high prognosis of NB patients. Some related reports have been made in previous studies. For example, Liu JA et al. pointed out that Fbxo9 plays an essential role in mediating Sox10 in destroying Neurog2 protein and directing neural crest progenitor cell lineage to glial cells [43]. Hynes-Smith RW et al. pointed out that the lack of FBXO9 enhanced proteasome activity and promoted the aggressiveness of acute myeloid leukemia [44]. Yang L et al. pointed out that the participation of circRNA HECW2 in the non-autophagy effect of ATG5 is related to endothelial-mesenchymal transition [ [54]. Although FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 have not been reported in NB, based on their role in tumor regulation, this also provides a potential value for their role in NB to a certain extent.
To study the effectiveness of these biomarkers in diagnostic NB, we used the microarray dataset GSE73517 to verify. Consistent with the above analysis results, FBXO9, HECW2 MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 were signi cantly higher in MYCN nonampli ed NB cells. Correlation analysis found that the seven effective biomarkers (FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16) related to MYCN. In addition, the diagnostic effectiveness of the NB biomarkers was veri ed by ROC analysis. The results also clearly indicated that FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 were of excellent speci city and sensitivity to diagnose NB.
To verify the results of bioinformatics analysis, we conducted in vitro experiments. First, we detected the mRNA levels of the hub genes in NB cells. The results showed that FBXO9, HECW2, MIB2, RNF19B, RNF213, TRIM36, and ZBTB16 have higher expression in MYCN non-ampli ed NB cells, which was consistent with our analysis. In addition, MYCN siRNA was transfected into BE2 cells, the expression of FBXO9, RNF19B, and TRIM36 was increased at the mRNA levels, which was consistent with the analysis results. After the MYCN overexpression plasmid was transfected into AS cells, the expression of FBXO9, RNF19B, and TRIM36 was decreased at the mRNA levels, which further con rmed our results. At present, although there is no report of a regulatory relationship between FBXO9, RNF19B, and TRIM36 and MYCN, the prediction based on the databases and the veri cation of mRNA levels essentially explained the interaction between FBXO9, RNF19B, and TRIM36 and MYCN. In addition, FBXO9, RNF19B, and TRIM36 are expected to become the novel target genes of MYCN, but this requires more experimental data to be con rmed.