Combinatorial targeting of MTHFD2 and PAICS in purine synthesis as a novel therapeutic strategy

MYCN-amplified (MNA) neuroblastoma is an aggressive neural crest-derived pediatric cancer. However, MYCN is indispensable for development and transcriptionally regulates extensive network of genes. Integrating anti-MYCN ChIP-seq and gene expression profiles of neuroblastoma patients revealed the metabolic enzymes, MTHFD2 and PAICS, required for one-carbon metabolism and purine biosynthesis were concomitantly upregulated, which were more susceptible to metastatic neuroblastoma. Moreover, we found that MYCN mediated the folate cycle via MTHFD2, which contributed one-carbon unit to enhance purine synthesis, and further regulated nucleotide production by PAICS in response to cancer progression. Dual knockdown of the MYCN-targeted gene pair, MTHFD2 and PAICS, in MNA neuroblastoma cells synergically reduced cell proliferation, colony formation, migration ability, and DNA synthesis. By systematically screening the compound perturbagens, the gene expression levels of MTHFD2 and PAICS were specifically suppressed by anisomycin and apicidin across cell lines, and our co-treatment results also displayed synergistic inhibition of MNA neuroblastoma cell proliferation. Collectively, targeting a combination of MYCN-targeted genes that interrupts the interconnection of metabolic pathways may overcome drug toxicity and improve the efficacy of current therapeutic agents in MNA neuroblastoma.


Introduction
Neuroblastoma is an embryonal tumor derived from the sympathetic nervous system and accounts for 15% of pediatric cancer mortality 1,2 . Although the outcome of high-risk patients has improved on the use of therapy, the majority eventually relapses and develops resistance with a 5-year event-free survival <50% 3,4 . High-risk (HR) neuroblastoma is characterized by V-myc avian myelocytomatosis viral oncogene neuroblastoma derived homolog (MYCN) amplification and segmental DNA copy number aberrations which are strongly correlated with poor prognosis, metastasis, and even treatment failure [5][6][7][8][9] . MYCN is a transcription factor that governs downstream signaling pathways for a broad range of fundamental processes during embryogenesis and development 10,11 . Amplification of MYCN is associated with tumor initiation and progression, and occurs in~20% of all neuroblastoma cases, with a higher incidence in INSS stages 3 and 4 compared to stages 1, 2, and 4S, suggesting a contribution in regulating oncogenic pathways 3,12,13 .
Alteration of cellular metabolism to maintain energy status for rapid cell progression has been considered as a feature of cancer cells 14,15 . MYCN has been closely tied to the regulation of neuroblastoma cell growth, and confers the serine-glycine-one-carbon pathway to promote metabolic reprogramming in HR neuroblastoma 16,17 . MYCN-amplified (MNA) neuroblastoma increases cell growth by glutamine transport and metabolism enhancement, thus depletion of glutamine causes the limitation of TCA cycle intermediates and eventually cell death 18 . Moreover, MYCN has been shown to stimulate mitochondrial biogenesis of β-oxidation and glucose metabolism in cancer progression 19,20 .
Cellular metabolism relies on highly regulated sequential enzyme-catalyzed reactions to acquire the essential components such as purines to generate DNA and RNA molecules 21 . Other than nucleotide formation, purines are the constituents of important biomolecules as energy and signal transduction cofactors to drive cellular processes 22 . Interestingly, a co-localization of purinosomes with mitochondria was found under purine-depleted environment suggesting a functional interaction between mitochondria and purine biosynthesis 23 . The mitochondrial folate enzyme methylenetetrahydrofolate dehydrogenase 2 (MTHFD2) is the most consistently overexpressed metabolic gene in cancer which plays a pivotal role in regulating the release of formate into cytoplasm 24 . In addition, serine is a precursor for the synthesis of glycine in one carbon pool by folate pathway at the mitochondria compartment which is further used up for purine biosynthesis 25 .
Collectively, we aimed to explore the relationship of MYCN-mediated pathways by integrating ChIP-seq and gene expression profiles to comprehensively dissect potential MYCN-targeted genes. Here, we elucidated that MNA neuroblastoma favors one carbon metabolism and purine biosynthesis by the up-regulation of MYCN-targeted metabolic genes MTHFD2 and PAICS. Our results also indicated that dual knockdown of MTHFD2/PAICS significantly reduced cell proliferation, colony formation, and migration which might diminish the aggressiveness and tumor progression ability in MNA neuroblastoma. Furthermore, combining the compound perturbagens to MYCN target gene pair displayed synergistic effects in MNA neuroblastoma cells which provides therapeutic opportunities for children with neuroblastoma.

Differential expression analysis
Microarray datasets of neuroblastoma samples were acquired from NCBI Gene Expression Omnibus (GEO) with accession number GSE45547, GSE3446, GSE19274, GSE16254, GSE12460, and GSE16237, including 1065 neuroblastoma. For multiple probes corresponding to a single gene, median expression value of those probes was taken. Total 11,939 common genes of six datasets were considered for downstream analysis. The expression values were log transformed and normalized by quantile normalization across all samples. HR neuroblastoma (age at diagnosis >18 months and stage 4) were classified into HR-MNA (n = 67) and HR-non-MNA (n = 156) groups based on the MYCN status. Significance Analysis of Microarrays (SAM) was used to identify differentially expressed genes between HR-MNA and HR-non-MNA with false discovery rate (FDR) <0.001 26 .
Public data sources and bioinformatics analysis MYCN-bound genes were obtained from our previous work 27 which ChIP-seq was used for genome-wide identification of MYCN regulatory networks. Two independent neuroblastma cohorts (SEQC and TAR-GET) were used for survival and correlation analyses. SEQC cohort was download from GEO with accession number GSE47792 and TARGET cohort was queried via GDC data portal (https://portal.gdc.cancer.gov/). The H3K4me3 and H3K27ac epigenetic profiles were obtained from ENCODE project. KEGG enrichment analysis was performed using the R/Bioconductor package clusterProfiler 28 .

RNA isolation of neuroblastoma cells and patient tissues, and cDNA synthesis
Neuroblastoma cells and tissue samples were homogenized by TRIzol reagent (Invitrogen). Total RNA of cells and tissue samples were prepared by Direct-zol™ RNA MiniPrep kit (ZYMO RESEARCH). RNA concentration and quality were determined by NanoDrop ND-1000 (NanoDrop Technologies). In all, 100 ng of total RNA of each tissue sample or 500 ng of total RNA template of each neuroblastoma cell line sample was reverse transcript to cDNA by using RevertAidTM H Minus First Strand cDNA Synthesis Kit (Thermo Scientific) with oligo (dT) 18 primer. The cDNA samples were stored at −30°C until use. The 21 neuroblastoma tissues were kindly provided by Dr. Wen-Ming Hsu (National Taiwan University Hospital, Taipei, Taiwan; IRB number 201407043RINC).

DNA manipulation and plasmid construction
The MYCN gene was amplified from synthesized cDNA as described previously (Thermo Fisher Scientific). PCR was performed to generate pCMV6-XL4 plasmids (Invitrogen) with a full-length sequence of MYCN. The resulting plasmid (pCMV6-XL4-MYCN) was transformed into Escherichia coli (E. coli) strain DH5α and selected by ampicillin. The plasmid DNA (pCMV6-XL4-MYCN) was prepared and sequenced at the DNA Sequencing Facility (Genomics BioSci. & Tech., Taipei, Taiwan). In all, 2.5 × 10 5 SK-N-AS cells were seeded 24 h before transfection in a 6-well plate using Lipofectamine 3000 (Invitrogen). Transiently transfected cells were harvested at 48 h posttransfection.

qRT-PCR analysis
The cDNA sample was amplified and applied by using CFX Connect™ Real-Time PCR Detection System (Bio-Rad Laboratories). The mRNA expression values were measured by ΔΔCt and normalized to GAPDH. The qRT-PCR primer sequences were listed in Supplementary  Table S1.

Luciferase reporter assay
Promoter regions of the MTHFD2 and PAICS genes were amplified using PCR and cloned into the pGL4.18 vector (Promega) flanked with NheI and HindIII sites.
The sequences of the promoter region primers are listed in Supplementary Table S2. SK-N-AS cells were seeded at 2.5 × 10 5 per 6-well plate for 24 h. Then SK-N-AS cells were co-transfected with either 500 ng of MTHFD2, PAICS promoter luciferase reporters or pGL4.18 empty vector along with 10 ng of pGL4.74 Renilla luciferase plasmid DNA together with 500 ng of MYCN expression plasmid (pCMV6-XL4-MYCN) or control vector (pCMV6-XL4). At 5 h post-transfection, cells were recovered in completed DMEM for 1 h and then cells were maintained in completed DMEM containing 1 μl/ml 70% ethanol or 1 μg/ml tetracycline and incubated for 48 h. At 48 h post-transfection, cells were lysed with passive lysis buffer for 15 min at room temperature and the firefly and Renilla luciferase activities were measured with the Dual-Luciferase Reporter assay system (Promega) according to the manufacturer's instructions.

Cell harvest and extraction for targeted metabolomics assay
Cells were grown in 15-cm culture dish, during which the medium was replaced every day (DMEM supplemented with 10% fetal bovine serum and 2 μg/ml puromycin for stable clones) at 37°C with 5% CO 2 before extraction. Cells were collected at 80% confluence and rapidly rinsed with warm 0.9% NaCl isotonic saline three times before quenching. Then, 1 ml of ice cold water was added and flash frozen in liquid nitrogen and detached using a cell scraper. Cell suspension were lysed by freeze-thaw cycles twice, followed by sonication. Protein concentration were measured by BCA Protein Assay Kit for normalization. Metabolism was quenched and metabolites were extracted by 75:25 methanol:water extraction buffer. In brief, 900 μl of −20°C methanol was add to 300 μl cell lysate and vigorously vortexed and then centrifuged at 14,000 × g for 20 min at 4°C. The supernatants were collected and dried using a centrifugal evaporator and stored at −80°C prior to analysis 29 .
An electrospray ionization interface was used to direct column eluent to the mass spectrometer. Mass spectra in the range m/z 50-500 were obtained by electrospray ionization in negative-ion mode. Initial instrument optimization was performed by infusing a mixture of L-serine

Cell cycle analysis
Stable shRNA expression cells were collected, fixed in 70% ethanol, and stored at −20°C overnight. The cells were washed with PBS and resuspended in PBS containing 100 μg/ml RNase A and 0.1% Triton X-100, then incubated at 37°C for 1 h. Cells were stained with 5 μg/ml propidium iodide (Santa Cruz). The DNA content of the cells was analyzed using a FACSCanto instrument (BD Biosciences Immunocytometry Systems). Ten thousand cells were collected for each measurement. The percentage of cells in different phases of the cell cycle was analyzed using ModFit LT (Verity Software House).

Cell proliferation, colony formation, and migration assays
For MTS analysis, stable shRNA expression cells (shLacZ, shMTHFD2#50, shMTHFD2#53, shPAICS#74, shPAICS#75, and shMTHFD2/PAICS) were seeded in 96well plate at 4000 cells per well with DMEM containing 2 μg/ml puromycin. At each 24 h, 48 h, 72 h time point, 20 μl of MTS/PMS (Promega) were added and incubated for 2 h at 37°C with 5% CO 2 . The cell viability was determined by OD 490 nm using an ELISA reader (BioRad). For colony formation assay, the stable shRNA expression cells were seeded in 6-well plates (600 cells/ well) with DMEM containing 2 μg/ml puromycin and incubated for 2 weeks. After 2 weeks, the colonies were fixed with methanol overnight and then stained by 4% GIEMSA stain for 1 h and analyzed by ImageJ. For cell migration, 3 × 10 4 stable shRNA expression cells with 1% FBS medium were loaded into the inserts, and medium containing 10% FBS was loaded into the lower compartments of an 8-μM pore size Transwell plate (Corning). The cells were incubated at 37°C with 5% CO 2 for 6 h, then fixed for 30 min with 100% methanol and stained with 1% GIEMSA for 30 min. Cotton swabs were used to remove cells from the upper side of the inserts. Images of five different microscope fields of each insert were captured and counted.

Identification of potential compounds suppressing MTHFD2 or PAICS expression
Perturbation datasets corresponding to the treatment of various drugs or compounds for ten cell lines, including A375, A549, HA1E, HCC515, HEPG2, HT29, MCF7, NPC, PC3, and VCAP, were obtained from the Library of Integrated Network-based Cellular Signatures (LINCS) project 30 . To identify the compounds that significantly suppress the gene expressions of MTHFD2 and PAICS based on the perturbation datasets, computational approach was performed as detailed in our previous methods 31,32 . Briefly, the expression changes of all measured transcripts (n = 12,494) were ordered from the most downregulated (negative) to the most upregulated (positive) cross ten cell lines and assessed if the gene of interest was significantly enriched on the down-regulation edge cross all cell types under the perturbation of the given compound. Finally, the expressions of MTHDF2 and PAICS were found to be suppressed by anisomycin and apicidin, respectively. To test the synergistic effect of anisomycin and apicidin combination, the corresponding combination solutions were prepared proportionally at 1:1 molar ratio. The control groups received an equivalent amount of DMSO as corresponding treatment groups. At the end of 24 and 48 h of treatments, 20 μl of MTS/PMS (Promega) were added and incubated for 2 h at 37°C with 5% CO 2 . The cell viability was determined by OD 490 nm using an ELISA reader (BioRad) and normalized with the corresponding control groups. Drug synergy was evaluated using CompuSyn (http://www.combosyn.com) to calculate the combination index (CI) values: synergism (CI < 0.9), additive effect (CI = 0.9-1.1), and antagonism (CI > 1.1) 33 .

Cell apoptosis assessment
In all, 1 × 10 6 SK-N-DZ or SK-N-BE(2)-C cells were seeded in a 10-cm culture dish for 24 h before drug treatments. For apoptosis detection, 0.2 µM anisomycin, 0.2 µM apicidin, and the combination of solutions were prepared proportionally at 1:1 molar ratio. At the end of 24 h of treatment, cells underwent PBS washes before treated with the FITC Annexin V Apoptosis Detection Kit I (BD Pharmingen, USA) according to the manufacturer's instructions. Ten thousand stained cells were acquired on a FACSCanto instrument (BD Biosciences Immunocytometry Systems) for each measurement.

Statistical analysis
Results were expressed as the mean ± standard deviation (SD). The Student's t-test was used to analyze differences between two groups. The differences between groups were considered to be statistically significant when P < 0.05. All experiments were performed in triplicate.

Results
Integrated analysis of multi-layer omics data reveals higher cellular metabolism in MNA neuroblastoma To identify highly confident MYCN-regulated genes, we combined genes that differentially expressed between MYCN status with MYCN-bound genes (Fig. 1a). Here, we compiled a neuroblastoma gene expression cohort from six independent studies and identified 4275 differentially expressed genes between HR-MNA and HR-non-MNA samples. By integrating with 874 MYCN-bound genes that identified in our previous anti-MYCN ChIPseq experiment 27 , of which 427 common genes were found to be regulated by MYCN (Fig. 1b and Supplementary Table S4), and we termed these genes as potential MYCN-regulated genes. A few of these genes have been verified as direct transcriptional targets of MYCN, such as TWIST1, DLL3, and DDX1 (Fig. 1c) [34][35][36] . The functional enrichment analysis of the potential MYCN-regulated genes revealed cancer hallmarkassociated biological processes such as cell cycle and metabolic pathways ( Fig. 1d and Supplementary  Table S5).
Previous studies have revealed that amplification of MYCN increases neuroblastoma cell growth through glutathione biosynthesis and glutamine uptake, thus suggesting MYCN reprograms neuroblastoma metabolism 11,37,38 . In our analysis, one carbon pool by folate, pyrimidine metabolism and purine metabolism were significantly altered by MYCN (Fig. 1e). One carbon pool by folate metabolism regulates the conversion of serine to glycine in mitochondria and releases the production of folate to cytosol for purine synthesis. Interestingly, a colocalization of purinosome and mitochondria was also reported under purine-depleted condition 25 . To investigate whether MNA neuroblastoma favors one carbon metabolism and purine metabolism, we first assessed the intrinsic differences of MNA and non-MNA neuroblastoma cells for serine consumption and nucleoside production. Targeted metabolic assay was performed to quantify the relative levels of serine, nucleoside intermediate (AICAR), and nucleoside monophosphate (IMP and GMP) between MNA SK-N-DZ and non-MNA SK-N-AS cell lines ( Fig. 1f and Supplementary Fig. 1). The levels of AICAR, IMP, and GMP were significantly higher in MNA neuroblastoma compared to non-MNA neuroblastoma cells, suggesting MNA neuroblastoma may require higher demand in purine biosynthesis for rapid cell proliferation (Fig. 1g). In addition, we observed a lower level of serine in MNA neuroblastoma cells which could be converted into glycine in one carbon pool by folate pathway (Fig. 1g). Collectively, these data suggest that MYCN may play a role in regulating metabolism in neuroblastoma through potential MYCN-regulated genes.
MTHFD2 and PAICS expressions are strongly correlated to MYCN status in neuroblastoma We speculate that high levels of cellular metabolism are linked to neuroblastoma progression, thus two highly expressed metabolic genes MTHFD2 and PAICS among one carbon pool by folate and purine metabolism pathways in MNA neuroblastoma prompted us to further examine their relationship with MYCN status (Fig. 1a, c, e). We first evaluated the mRNA expression levels of MTHFD2, PAICS, and MYCN in 21 tissue samples from neuroblastoma patients (MNA: n = 9; non-MNA: n = 12). The mRNA expressions of MTHFD2 and PAICS were positively correlated to MYCN by qRT-PCR with a Pearson correlation coefficient (r) of 0.873 and 0.850 respectively (Fig. 2a). In consistent with our finding, the SEQC cohort results showed that mRNA expressions of MTHFD2 and PAISC were positively correlated with MYCN, and highly expressed in MNA patients ( Fig. 2b and Supplementary Fig. 2). Next, we assessed the protein expression levels in six neuroblastoma cell lines, where  Fig. 2c and Supplementary Fig. 3).
To further confirm the regulation of MYCN with MTHFD2 and PAICS, we assessed the expression response upon MYCN perturbation in MNA and non-MNA cell lines. Overexpression of MYCN was introduced by pCMV6-XL4-MYCN in non-MNA SK-N-AS cells, both PAICS and MTHFD2 were elevated along with MYCN (Fig. 2d). Contrastingly, silencing MYCN in MNA SK-N-DZ and SK-N-BE(2)-C cells using small interfering RNA (siRNA) revealed that PAICS and MTHFD2 expressions were concomitantly downregulated with MYCN (Fig. 2e). Moreover, the tetracycline inducible system was also applied to evaluate the effects on neuroblastoma. The Tet-21N cells express MYCN in the absence of tetracycline, whereas the addition of tetracycline inhibits the tTA to bind to the promoter of MYCN thus gene expression is turned off (Fig. 2f). After 48 h of tetracycline treatment, both MTHFD2 and PAICS were significantly downregulated following MYCN depletion at mRNA and protein levels indicating that MYCN might f The Tet-21N cells express MYCN in the absence of tetracycline, where the addition of tetracycline inhibits tTA to bind to the promoter of MYCN thus gene expression is turned off. g, h After 48 h of tetracycline treatment, both MTHFD2 and PAICS were significantly downregulated following MYCN depletion at mRNA (g) and protein levels (h). The relative mRNA expressions were calculated by qRT-PCR analysis, and normalized to GAPDH gene. ***P < 0.001 influence the transcription efficiency of MTHFD2 and PAICS (Fig. 2g, h).

MYCN transcriptionally targets MTHFD2 and PAICS in neuroblastoma
Our previous anti-MYCN ChIP-seq data revealed that MTHFD2 and PAICS had significant MYCN binding signals comparing with IgG control (Supplementary Fig.  4). To examine whether MYCN directly drives MTHFD2 or PAICS, we applied a dual luciferase reporter assay to determine the MYCN binding sites on MTHFD2 and PAICS promoters (Fig. 1a). According to the ChIP-seq results and examination of consensus E-box motifs, four MTHFD2 (M1, M2, M3, and M4) and three PAICS (P1, P2, and P3) promoter sequences were constructed using the pGL4.18 luciferase vector and co-transfected with pCMV6-XL4-MYCN in MNA SK-N-AS cell lines (Fig. 3a,  b). The luciferase activity of M1, M3, and M4 were significantly enhanced, of which M3 that consisted the MYCN ChIP binding peaks resulted the highest activity with a 42-fold change relative to the empty pGL4.18 vector (Fig. 3c). On the other hand, P1, P2, and P3 constructs containing the ChIP-seq binding peaks were increased by the presence of MYCN vector compared to empty pGL4.18 (Fig. 3d). Moreover, we compared the MYCN overexpressed SK-N-AS with SK-N-AS control (Fig. 3e), and also observed that the binding sequences of MTHFD2 and PAICS constructs were significantly upregulated by MYCN expression (Fig. 3f, g).

MTHFD2 and PAICS regulate purine biosynthesis in MNA neuroblastoma
To explore the impact of the two metabolic enzymes on MNA neuroblastoma, we first constructed stable knockdown of LacZ (control), MTHFD2, PAICS, and MTHFD2/PAICS (dual-knockdown) by delivering shRNA hairpins into MNA SK-N-DZ cells followed by puromycin selection (Supplementary Fig. 5). In one carbon pool by folate pathway, serine is catabolized in mitochondria to produce formaldehyde, and then transforms into essential molecule, 10-formyl-THF, for purine biosynthesis by MTHFD2 39 . To investigate the relationship of one carbon pool by folate and purine biosynthesis in neuroblastoma, we analyzed the relative abundance of metabolites using liquid chromatography mass spectrometry ( Fig. 1a and Supplementary Fig. 6). We found shMTHFD2, shPAICS, and shMTHFD2/PAICS SK-N-DZ cells decreased the consumption of serine, and this observation implies that knockdown of either MTHFD2 or PAICS suppresses the metabolic activity in one carbon pool by folate pathway (Fig. 4a, b).
On the other hand, 5-aminoimidazole-4-carboxamide ribonucleotide (AICAR), the final intermediate for assembling inosine monophosphate (IMP), and guanosine monophosphate (GMP) 40 were markedly reduced in shMTHFD2/PAICS relative to single knockdown. Thus, we suggest that the depletion of AICAR level limits the production of IMP and GMP when both MTHFD2 and PAICS are suppressed (Fig. 4a, b). Of note, the single knockdown of either MTHFD2 or PAICS elevated the levels of AICAR compared to the control. The results suggest that the synergistic effect of MTHFD2 and PAICS are related to the biosynthesis of nucleotides, as well as their interconnection in one carbon by folate and purine metabolism in MNA neuroblastoma cells.
Since the targeted metabolomics analysis indicated an imbalance of purine nucleotides in dual knockdown cells, we hypothesized that the cell cycle perturbations might be induced. Flow cytometry for DNA content was conducted and the distribution of cells in different phases of the cell cycle was analyzed. The results showed that the percentage of cells in the S phase of shMTHFD2 and shPAICS cells increased, by 4.29% and 4.34%, respectively, and that of dual knockdown cells increased by 7.70% relative to shLacZ control of MNA SK-N-DZ cells, indicating an interruption of the transition from S to G2/M phases (Fig.  4c, d). The percentage of cells decreased in both G1 and G2/M phases, along with increased cell population of S phase, especially the dual knockdown cells which decreased by 5.33% and 2.36% of population in the G1 and G2/M phases, respectively (Fig. 4c, d). Thus, purine intermediates depletion and S phase arrest may be correlated with dysregulation of DNA synthesis.

Dual knockdown of MTHFD2 and PAICS suppresses neuroblastoma cell growth
To elucidate the prognostic implications of MTHFD2 and PAICS in neuroblastoma progression, we analyzed the association of MTHFD2 and PAICS expressions with patient survival. Poor overall and event-free survivals were found in highly expressed MTHFD2 or PAICS neuroblastoma patients (Fig. 5a, b). According to the impact of MTHFD2 and PACIS on neuroblastoma patient survival, we ascertain their stimulatory effects on cell proliferation by cell viability and colony formation assays in response to MTHFD2 and PAICS expression manipulation (Figs. 1a, 5c-h and Supplementary Fig. 7). Cell viability was measured at 24, 48, and 72 h by MTS assays, and the results indicated a noticeable decline in cell proliferation in response to MTHFD2 and PAICS. Cell proliferation was strongly affected by MTHFD2 and PAICS; the viability of shMTHFD2 cells declined by 20. 4%, 20.8%, and 29.7%, at 24, 48, and 72 h, and the viability of shPAICS cells also decreased by 19.2%, 17.2%, and 24.0% at 24, 48, and 72 h respectively (Fig. 5d, e). Either shRNA-mediated knockdown of MTHFD2 or PAICS reduced cell viability compared with shLacZ control (Supplementary Fig. 8). However, the most significant reduction in cell viability occurred in dual knockdown of MTHFD2 and PAICS (shMTHFD2/PAICS) where the cell viability declined 48.4%, 50.2%, and 60.3% at 24, 48, and 72 h respectively (Fig. 5d, e).
Colony formation was monitored to investigate the longterm effects of MTHFD2 and PAICS on neuroblastoma cell proliferation (Fig. 5f-h). Single knockdown of MTHFD2 or PAICS alone was associated with~32.1% and 29.6% declension of colony formation ability respectively, whereas dual knockdown of both MTHFD2 and PAICS resulted in a decrease of colony number to 46.5% compared with shLacZ cells (Fig. 5h). These observations suggest that simultaneous knockdown of MTHFD2 and PAICS further impair the MNA cell proliferation, and they may serve as a regulatory purpose in cell growth together.

Dual knockdown of MTHFD2 and PAICS suppresses neuroblastoma cell migration
To further investigate the impacts of MTHFD2 and PAICS on tumor aggressiveness, we analyzed the relationship of MTHFD2 and PAICS expressions with tumor stages (INSS). A higher MTHFD2 or PAICS expression was found in stage 4 compared with stages 1-3 and 4S, which had a similar expression pattern as MYCN in INSS stages (Fig. 6a). The above results suggest that MTHFD2 and PAICS are correlated with metastatic stage and tumor progression. We also observed a differentiated cell morphology in dual knockdown shMTHFD2/PAICS cells in comparison with single knockdown or shLacZ control cells which is linked to a less aggressive neuronlike phenotype in neuroblastoma (Fig. 6b) 4,11 . The shLacZ control cells exhibiting typical rounder morphology compared with single or dual knockdown MNA neuroblastoma cells, whereas the dual knockdown shMTHFD2/PAICS cells exhibited a noticeable outgrowth of neurite lengths and numbers of branch points (Fig. 6b, c). These morphological changes, together with the correlation of MTHFD2 or PAICS expression to INSS, led us to probe whether MTHFD2 and PAICS play roles in malignant transformation and tumor aggressiveness.
Cell migration is a crucial process during metastasis, thus we performed transwell assay to determine the migration ability on MTHFD2 and PAICS (Fig. 6d-f). Dual knockdown of MTHFD2 and PAICS had a significantly poorer migration ability compared with single knockdown of MTHFD2 and PAICS alone where shMTHFD2/PAICS was dramatically decreased by 69.7% compared to shLacZ control (Fig. 6f). As in the single knockdown results, shMTHFD2 cells decreased by 49.3%, while that of shPAICS cells decreased by 38.4% compared to shLacZ control, suggesting silencing of MTHFD2 and PAICS had synergistic effect on migration ability in MNA neuroblastoma (Fig. 6f).

Compound perturbagens suppressing MTHFD2 and PAICS show synergistic effects on MNA neuroblastoma
Anti-metabolic agents have been used in cancer treatment for decades, however, a relatively low number of metabolic inhibitors developed for cancer therapy due to the adverse effects on patients by blocking the synthesis of vital cellular constituents from metabolic pathways 41 . Our strategy for identifying anti-cancer agents is to screen for compounds that particularly target MYCN-regulated metabolic genes across cell lines by using transcriptomic profiles 31,32 . By matching gene signatures with drug signatures, we systematically identified the gene expressions of MTHFD2 and PAICS were significantly suppressed by . h Single knockdown of MTHFD2 or PAICS alone was associated with approximately 32.1% and 29.6% declension, respectively, whereas dual knockdown of both MTHFD2 and PAICS resulted in a decrease of colony number to 46.5% compared with shLacZ cells anisomycin and apicidin across ten cell lines ( Fig. 7a and Supplementary Fig. 9). Anisomycin is a bacterial antibiotic which inhibits protein biosynthesis by blocking peptide bond formation in ribosomes, and activates kinase cascades including JNK, p38 MAPK signaling pathways [42][43][44] . Apicidin is a fungal metabolites shown to inhibit histone deacetylase, and induce p21(WAD1/Cip1), caspase-3 and caspase 9, and DNA fragmentation 45,46 .
To investigate the effects of anisomycin combined with apicidin on MNA neuroblastoma, drug combination assays were performed using SK-N-DZ and SK-N-BE(2)-C cells with various concentrations for 24 and 48 h (Fig.  7b-e and Supplementary Fig.10). The drug combination assay showed synergistic effects on the inhibition of MNA neuroblastoma cell proliferation, as compared with either anisomycin or apicidin treatment alone (Fig. 7b, d). Furthermore, the mean CI values were 0.71 and 0.74 in SK-N-DZ and SK-N-BE(2)-C, respectively, which indicated synergistic drug combinations (Fig. 7c, e). In addition, the non-MNA SK-N-AS and SK-N-SH neuroblastoma cells displayed a less sensitive effects to the anisomycin and apicidin treatments ( Supplementary Fig. 11). Anisomycin and apicidin have been reported to induce apoptosis in certain cell types 46,47 , however, the effects of drug combination on apoptosis in neuroblastoma remains unidentified. Therefore, we performed apoptosis assays and found that both anisomycin and apicidin treatments induced apoptosis in MNA neuroblastoma (Fig. 8a-d).
Notably, the co-treatment results exhibited the most significant cell apoptosis (Fig. 8b, d).
Next, we examined the MTHFD2/PAICS gene combinations with synergistic effects on neuroblastoma patient survival. First, the median expression value of the given gene was used to define the high and low expressions, and then we split neuroblastoma patients into three groups: high-high expression of MTHDF2 and PAICS  (Fig. 8e, f). Altogether, the patient survival yielded consistent results as the drug combination assays that simultaneously suppressing MTHFD2 and PAICS expressions showed a synergism in MNA neuroblastoma patients.

Discussion
MYCN serves as an oncogenic regulator that mediates aberrant signaling pathways and leads to HR neuroblastoma. Recent studies achieved some clinical success to indirectly target MYCN, since developing direct MYCNtarget agent encounters difficulties and has not yet been established 38,48,49 . Here, we addressed highly confident Cancer metabolism is a complex enzyme-catalyzed reaction that involves multiple factors to overcome any single disruption 50 . Thus, we aim to investigate the MYCN-regulated interconnection of metabolic pathways and biological processes. The potential MYCN-targeted genes are found to be involved in purine biosynthesis which consists ten steps and six enzymes from phosphoribosyl pyrophosphate (PRPP) to inosine monophosphate (IMP), guanosine monophosphate (GMP), and adenosine monophosphate (AMP) 51 . Purinosome was coined to describe the reversible six-enzyme complex for purine biosynthesis that assembled and co-localized on mitochondria while the purine levels of cytosol depleted 23,52 . Besides, one carbon pool by folate was also identified to be associated with MNA neuroblastoma, which generated formates, the precursor to 10formyltetrahydrofolate in the mitochondria and exported into cytosol for the need of transformylase enzymes, GRAT and ATIC in purine biosynthesis 21 . On the other hand, loss of MTHFD2 would cause folate trapping within the mitochondria, therefore, folate-dependent pathways may be altered 24,53,54 .
Previous studies indicate that amino acid metabolism, nucleic acid metabolism, and mitochondrial metabolism play important roles in cancer progression 55,56 . Biosynthetic pathways are essential for cancer metabolism since they generate energy and macromolecules that are required for replicative cell division and tumor growth 56 . Our targeted metabolomics analysis indicated that the levels of IMP and GMP were significantly higher in MNA neuroblastoma compared to non-MNA neuroblastoma cells, suggesting MNA neuroblastoma may require higher demand in purine biosynthesis for rapid cell proliferation. The MYCN-targeted gene pair results implied that MYCN could regulate the interconnection of metabolic pathways by MTHFD2 and PAICS in two directions. MYCN transcriptionally regulates the mitochondrial metabolic enzyme MTHFD2 for formate exportation and the bifunctional enzyme PAICS of the purinosome to generate SAICAR for purine biosynthesis.
Treatment that combines two or more specific therapeutic targets might reduce the dosage and side-effect for aggressive therapy in children 57 . Our dual knockdown and co-treatment assays that target MTHFD2 and PAICS had synergistic effects on MNA neuroblastoma which might diminish the aggressiveness and tumor progression ability. In conclusion, establishing MYCN downstream signaling pathways in neuroblastoma may improve our capacity in identifying molecular interactions, disease pathways, and even targeted drug discovery.