Transcriptome profiling of caspase-2 deficient EμMyc and Th-MYCN mouse tumors identifies distinct putative roles for caspase-2 in neuronal differentiation and immune signaling

Caspase-2 is a highly conserved cysteine protease with roles in apoptosis and tumor suppression. Our recent findings have also demonstrated that the tumor suppression function of caspase-2 is context specific. In particular, while caspase-2 deficiency augments lymphoma development in the EμMyc mouse model, it leads to delayed neuroblastoma development in Th-MYCN mice. However, it is unclear how caspase-2 mediates these differential outcomes. Here we utilized RNA sequencing to define the transcriptomic changes caused by caspase-2 (Casp2−/−) deficiency in tumors from EμMyc and Th-MYCN mice. We describe key changes in both lymphoma and neuroblastoma-associated genes and identified differential expression of the EGF-like domain-containing gene, Megf6, in the two tumor types that may contribute to tumor outcome following loss of Casp2. We identified a panel of genes with altered expression in Th-MYCN/Casp2−/− tumors that are strongly associated with neuroblastoma outcome, with roles in melanogenesis, Wnt and Hippo pathway signaling, that also contribute to neuronal differentiation. In contrast, we found that key changes in gene expression in the EμMyc/Casp2−/− tumors, are associated with increased immune signaling and T-cell infiltration previously associated with more aggressive lymphoma progression. In addition, Rap1 signaling pathway was uniquely enriched in Casp2 deficient EμMyc tumors. Our findings suggest that Casp2 deficiency augments immune signaling pathways that may be in turn, enhance lymphomagenesis. Overall, our study has identified new genes and pathways that contribute to the caspase-2 tumor suppressor function and highlight distinct roles for caspase-2 in different tissues.


Introduction
The role of caspases in tumor suppression has been largely ascribed to their established function in cell death 1 . Recent experimental evidence also suggests that nonapoptotic roles for caspases, including cell proliferation, inflammation, migration, or invasion, can contribute to tumorigenesis 2,3 . These features are well-established hallmarks of cancer, both independently and in co-operation with cell death evasion 4 , and may determine caspase contribution to both tumor suppression and tumor promotion in certain contexts 3 . While the mechanisms that regulate caspase functions in tumorigenesis are still unclear, there is evidence to suggest that caspases mediate their differential functions in a developmental and tissue-specific manner 1,2 .
Correlative evidence supports a tumor suppressor function for caspase-2. The human CASP2 gene on Ch7q34, is part of a region frequently deleted in hematological malignancies 24 . In addition, reduced CASP2 expression has been reported in lymphoma and leukemia and correlates with poor prognosis in AML and ALL 25,26 . Somatic mutations in CASP2, although rare, are found in cases of high-grade colon and gastric cancer, lung, skin, and breast cancer 19,27,28 . Furthermore, in colorectal cancer reduced CASP2 expression is caused by loss of its transcriptional regulator BCL9L, and is a key cause of aneuploidy tolerance, tumor progression, and resistance 29 . In contrast, Casp2 deficiency in mice does not affect tumor onset or progression following 3methylcholanthrene (3-MC)-induced fibrosarcoma or irradiation-driven lymphoma 20 . Our previous studies have also demonstrated that Casp2 loss delays neuroblastoma development in the Th-MYCN mouse model, and that low CASP2 expression correlates with increased survival in human neuroblastoma 30 . These findings suggest that caspase-2 has a context/tissue-specific function in tumor suppression and that both tissue and genomic variability can cooperate with caspase-2 to determine tumor outcome. To understand how caspase-2 mediates its differential functions in tumor suppression, it is important to determine if disruption of additional genes co-operate with caspase-2 in tissue-specific contexts.
In this study, we comparatively analyzed the transcriptomes of Casp2 −/− tumors from EμMyc and Th-MYCN mice to identify changes in transcriptional tumor networks that are influenced by caspase-2 deficiency. While we identified several unique genes that were aberrantly expressed in the Casp2 −/− tumors, we found that Megf6/EGFL3 was the only gene that was differentially expressed in the two tumor types. Our study also identified several enriched pathways and gene signatures specific to EμMyc/Casp2 −/− or Th-MYCN/Casp2 −/− tumors that may be associated with enhanced EμMycinduced lymphomas and/or delayed Th-MYCN-mediated neuroblastoma. The cross-tumor-specific transcriptional aberrations are highly indicative of distinct roles for caspase-2 during neuronal and B-cell development that perhaps influence its tumor suppressor function. Comparisons of the differentially expressed genes (DEGs) were performed between (i) WT and Casp2 −/− tumors to identify genes associated with tumorigenesis following Casp2 loss and (ii) between the EμMyc/Casp2 −/− and Th-MYCN/Casp2 −/− tumor samples, to narrow down candidate genes affected by Casp2 deficiency and/or that enhance lymphomagenesis in EμMyc/Casp2 −/− mice. These DEGs are summarized as volcano plots (Fig. 1b,  Supplementary Tables S2a-f). From the heat-maps, it is clear that Casp2 deficiency changes the expression of many genes, with more downregulated genes in Th-MYCN/ Casp2 −/− compared to Th-MYCN tumors (Fig. 1c) and more upregulated genes in EμMyc/Casp2 −/− versus EμMyc tumors (Fig. 1d, Supplementary Tables S3a-b).

RNA-seq analysis of
In total, 13,714 genes were included for analysis in each comparison group (Supplementary Tables S2a-f  Th-MYCN/Casp2 −/− comparison identified four unique downregulated genes (Ncor2, Diras2, Ednrb, Shb) and one upregulated (Eya2) gene, with roles in transcription (Ncor2), signal transduction (Ednrb, Shb), melanogenesis (Ednrb), and DNA damage repair (Eya2). Two of the unique EμMyc/Casp2 −/− genes have associated roles in transcription regulation (Six5) and calmodulin binding (Nrgn) and have been previously associated with lung squamous cell carcinoma, breast cancer (Six5 31,32 ), and T-cell lymphoma (Nrgn 33 ). Tdrd5 was excluded as it was only increased in one out of three biological replicates. Significant downregulation of Casp2 was also verified in the EμMyc/Casp2 −/− comparison and Casp2 gene deletion was confirmed by genotyping and immunoblotting in all tumor samples.
To identify common genes that are significantly altered in the Casp2 −/− samples, the DEGs were subdivided into upregulated and downregulated lists, illustrated in a Venn diagram (Fig. 2a). This identified 86 uniquely down-  Table S4b). There was also a single overlapping gene; Megf6 (Multiple EGF-Like Domains 6) differentially altered in EμMyc/Casp2 −/− and Th-MYCN/Casp2 −/− tumors, and this was validated by quantitative PCR in additional tumor samples (Fig. 2b).
We used R2: Genomics Analysis and Visualization Platform, to examine the correlation between Megf6 transcript expression and clinical outcome in various human neuroblastoma and B-cell lymphoma publicly available expression array data sets. Our previous data showed that Casp2 levels correlated with clinical outcome in a subset of MYCN non-amplified human neuroblastomas 30 . Here, we also found that higher Megf6 transcript levels are associated with poor outcome in this neuroblastoma subset (Fig. 2c), and this is consistent with our RNA-seq data and the delayed tumorigenesis observed in the Th-MYCN/Casp2 −/− mouse model 30 . We did not find any significant association in MYCN-amplified neuroblastoma. In contrast, we found there was a trend for lower Megf6 expression being predictive of poorer B-cell lymphoma outcome (e.g., Xiao -420, fRMA-u133p2 dataset; P = 0.061) (Fig. 2d). Notably, high Megf6 expression in B-cell lymphomas also showed a somewhat poor survival outcome, indicating that Megf6 levels are probably not predictive of lymphoma progression.  Table S5a). The top BP was cell differentiation (GO:0030154) (Fig. 3a, Supplementary Table S5a), with many of these genes also being associated with the BPs of transcriptional regulation from RNA polymerase II promoter (GO:0000122) and multicellular organismal development (GO:0007275). This indicates that Casp2 loss can affect neuronal differentiation and gene transcription in neuroblastoma. Of note, response to hypoxia (GO:0001666) and ischemia (GO:0002931), processes previously associated with caspase-2 function 34-37 were also enriched.
Consistent with the enriched BPs, KEGG pathway analysis identified T-cell receptor signaling and primary immunodeficiency as the most significantly altered pathways followed by cell adhesion molecules and Rap1 signaling (Fig. 4b, Supplementary Table S6b). All pathways were significantly enriched when analyzed using GSEA (Fig. 4a, c), suggesting that caspase-2 may have a role in regulating T-cell signaling and/or Rap1 signaling to influence lymphoma development in EμMyc mice.

Differential expression of genes in caspase-2 deficient tumors associated with increased survival in human neuroblastoma
To identify genes associated with delayed neuroblastoma development in Th-MYCN/Casp2 −/− mice, we examined various molecular biomarkers of neuroblastoma development and prognosis 39,40 . The expression  Table S7a), we did note that several genes adjacent to Th on chromosome 7 (Ascl2, Igf2, and H19) were also somewhat decreased (Supplementary Table S7b). This suggests this gene region may be affected in some cells.
Gene ontology analysis of the significantly altered neuroblastoma genes indicated their involvement in synaptic transmission, monoamine transport and cell growth regulation (Fig. 5c, Supplementary Table S7b). Interestingly, many of these genes are associated with melanogenesis, dopaminergic synapse, Wnt, and Hippo signaling (Fig. 5d, Supplementary Tables S7a and S7c). Downregulation of melanogenesis is associated with reduced expression of the transcription factor Mitf, along with Wnt signaling pathway components (Fzd1, Lef1, Camk2a). We noted a trend of altered expression of other Wnt signaling components (Wnt5a, Dkk2, Bmp6; FDR > 0.1) and their regulators, including Gata2 and its target homeobox transcription factors Six3 and Hesx1 (Supplementary Tables S2a, S7a), which also have roles in neural determination 42,43 . Increased expression of Wnt5a and Fzd1 and decreased expression of Lef1 and Camk2a were validated in different tumor samples by quantitative PCR (Fig. 6a, b). Together, these data identify altered expression levels of several Wnt-associated genes and suggest a fine-tuning of Wnt signaling components in  Table S8a). Many oncogenes have previously been associated with EμMyc-mediated tumor progression and/or B-cell lymphoma [47][48][49][50] . Several of these genes are clustered on nearby chromosome regions, including Ch.4 (Megf6-Prdm1, Arhgef16), Ch.6 (Cd4-Lag3, Cd8b-Cd8a), Ch.9 (Cd3d, Cd3g, Cd3e), Ch.10 (Dgka-Pme1), and Ch.11 (Sgsh, Slc26a11, Card14) (Supplementary Table S8b) indicative of possible common regulatory elements or amplification of these regions in the EμMyc/Casp2 −/− tumors. These findings highlight changes in the expression of multiple cancer-associated genes and/or gene regions that likely contribute to enhanced EμMyc lymphomagenesis onset.
Analysis of the 147-gene signature in EμMyc/Casp2 −/− tumors, identified axonogenesis, T-cell signaling, protein kinase B signaling, and signal transduction as the top enriched BPs and MFs (Fig. 7b, Supplementary  Table S8c). Consistent with this, T-cell receptor signaling, primary immunodeficiency and hematopoietic cell lineage were the top enriched KEGG pathways (Fig. 7c, Supplementary Table S8d), with Rap1 signaling and cell adhesion molecules also identified from our total DEG list (Fig. 4e). This indicates the DEGs in EμMyc/Casp2 −/− tumors are largely established cancer-associated genes.
Interestingly, the pathways melanogenesis and Wnt signaling were identified in the DEGs from both  (Fig. 5d) tumor comparisons. Increased expression of Tcf7, and its co-operating transcription factor Lef1, were further confirmed by quantitative PCR analysis in different EμMyc/Casp2 −/− tumor samples (Fig. 6a).  Table S5d). Analysis of the downregulated genes (n = 374) identified 13 unique BPs, with neuron maturation being the most significantly decreased, perhaps consistent with its primary role in Th-MYCN/Casp2 −/− tumors (Benjamini corrected P < 0.05). Other notable downregulated pathways (e.g., response to fatty acid, ubiquitination, negative regulation of epithelial cell proliferation, and cell cycle arrest) highlight multiple uniquely deregulated processes, previously been linked to Casp2 loss, and suggest their role in contributing to enhanced lymphomagenesis in EμMyc/Casp2 −/− mice. KEGG pathway analysis identified increased Rap1 signaling and cytokine-cytokine receptor interaction as the most significantly enriched and unique pathways (P<0.001, Fisher's exact test) (Fig. 8b, Supplementary  Tables S6c-d). In addition to T-cell receptor signaling identified from the above analysis (Fig. 4b), enrichment of NFκB signaling and cytosolic-DNA sensing were also confirmed by GSEA (Fig. 8c) (Supplementary  Table S6d). Interestingly, endocytosis was the top pathway associated with downregulated genes (Benjamini  Supplementary Table S6d. Dotted line indicates significance cut off (-log10(P-value) Fisher's exact test). c GSEA of the upregulated and downregulated pathways from (b). Black bars represent the position of members of the category in the ranked list together with the running enrichment score (plotted in green). NES = normalized enrichment score, FDR = false discovery rate corrected P = 0.049), a process that has emerged as a key hallmark of cancer, by affecting cell surface expression of critical molecules and signaling processes 53 . Calcium signaling, and Hippo signaling pathways were also uniquely enriched (Fig. 8b, Supplementary Table S6d). Consistent with this, reduced expression of Camk2a and Camk2b was validated in different EμMyc/Casp2 −/− tumor samples (Fig. 6b). These findings highlight several new pathways that are deregulated in EμMyc/Casp2 −/− tumors that can contribute to tumorigenesis.

Discussion
We used transcriptome profiling to characterize gene expression changes and pathways affected by Casp2 deficiency, in tumors from Th-MYCN neuroblastoma 30 and EμMyc lymphoma transgenic mouse models 7 . Our study has identified specific caspase-2 expression signatures in the different tumor types that likely contribute to the distinct tumor outcomes 7,30 . The diversity in gene expression highlights the (1) heterogeneity in tumor cell populations and/or clonal co-operation between subclones within the same tumor (2) tumor type-specificity (3) aberrations in components of several different signaling pathways, and/or (4) subtle differences that make it difficult to detect single causative pathways.
While the changes in gene expression between the Th-MYCN/Casp2 −/− and EμMyc/Casp2 −/− comparison groups were distinct, it was interesting that Megf6 was the only gene differentially expressed in the Casp2 −/− tumors. Megf6 is implicated in neural system disorders, such as ataxia 54 and recently, the epithelial-to-mesenchyme transition (EMT) to promote metastasis in colorectal cancer via TGFβ signaling 55 . Consistent with a tumor progression function, we found higher expression of Megf6 in EμMyc/Casp2 −/− and lower expression in Th-MYCN/Casp2 −/− tumors. Lower levels of Megf6 expression also correlated with better survival outcome in human neuroblastoma, suggesting it may contribute to delayed neuroblastoma in the Th-MYCN/Casp2 −/− mouse 30 . The lack of clinical correlation with Megf6 expression in B-cell lymphoma may be indicative of a tissue-specific function for Megf6 and/or that Casp2 loss can cooperate with high Megf6 expression to augment lymphomagenesis in the EμMyc/Casp2 −/− model. Nevertheless, these findings provide premise to further investigate Megf6 role and potential co-operation with caspase-2 function in tumorigenesis.
The delayed neuroblastoma onset in Th-MYCN/Casp2 −/− mice, previously highlighted the tissue/context specific role for caspase-2 in tumor suppression 30 . Our study has now identified cell differentiation as the most significantly altered biological process in Th-MYCN/Casp2 −/− tumors, associated with downregulation of genes that regulate neuronal differentiation and cell fate determination (Ascl2, Gldn, Mitf, Notch3, Shb, Syt17) and many being transcription factors (Ascl2, Batf3, Cbfa2t3, Elf4, Gata2, Mitf). A role for caspase-2 in neuronal differentiation, remodeling and protection against ischemic injury, has been previously implicated 34,56,57 . Silencing Casp2 can increase expression of neuronal differentiation markers 57 , and increased differentiation in Casp2 −/− neurons, would be consistent with favorable neuroblastoma outcome 40 and delayed tumor onset in Th-MYCN/Casp2 −/− mice 30 . Increased neuronal differentiation is also associated with high ROS levels 58 , a feature of Casp2 −/− cells 14,15 , that can further contribute to neuronal differentiation 58 . Identification of caspase-2 substrate(s) and/or interacting partners in differentiating NCCs will be key to defining its role in neurons and elucidate its function in neuroblastoma.
Our analysis identified a 41-gene signature in Th-MYCN/Casp2 −/− tumors, associated with neuroblastoma outcome, and involved in altered melanogenesis, Wnt and Hippo pathway signaling. Melanogenesis is known to be inhibited by excessive ROS 59 and may therefore be affected by the increased oxidative stress and ROS in Casp2 −/− cells 14,15 . Melanogenesis is also associated with aggressive neuroblastoma, characterized by higher tyrosinase activity and associated increase in DOPA synthesis [60][61][62] . The key regulator of tyrosinase is Mitf, was significantly reduced in Th-MYCN/Casp2 −/− compared to Th-MYCN tumors, and may suggest a defect in DOPA synthesis in Th-MYCN/Casp2 −/− tumors, consistent with the reduction in Dopaminergic synapse signaling. In addition, Wnt signaling regulates melanogenesis by stabilization and signaling through Mitf 59 . While a role for Wnt signaling in neuroblast differentiation has been described 43,63 , it has diverse functions in cell proliferation, polarity and migration, so its role in neuroblastoma is likely complex 63 . Nevertheless, the finding here, show that loss of Casp2 affects several components regulating the Wnt signaling pathway in neuroblastoma, that could impact neuronal differentiation and possibly neuroblastoma onset in Th-MYCN/Casp2 −/− mice.
Our previous studies demonstrated that loss of Casp2 accelerates EμMyc-driven lymphomagenesis 7 . This study has now identified aberrant expression of 147 cancerassociated genes including increased expression of 72 oncogenes and reduced expression of 9 TSG in the EμMyc/Casp2 −/− tumor samples. While many of these genes may play co-operating roles with Casp2 loss to enhance EμMyc-mediated lymphomagenesis, aberrant expression of some genes are likely a consequence of enhanced tumor growth caused by Casp2 loss. Interestingly, the key pathways associated with EμMyc/Casp2 −/− tumors, was increased immune response, in particular Tcell activation and signaling pathways. Higher expression of CD4, CD25/IL2ra and Lag3 in EμMyc/Casp2 −/− tumors, are indicative of activated regulatory T-cells (Tregs) at the tumor site, which are associated with both favorable and unfavorable prognosis in B-cell lymphomas 64 . They can facilitate tumor immune evasion by directly suppressing B-lymphoma cells, and also function to inactivate tumor specific CD4 + /CD8 + T-cells, leading to an immunosuppressive network associated with increased T-cell tolerance and reduced T-cell mediated killing 64 . This is a key factor in uncontrolled tumor growth and poor clinical outcome in many solid tumors 64 . Several clinical studies have also described the presence of T-cell markers (including CD3, CD4, CD8) in diffuse large B-cell lymphoma (DLBCL), but the prognostic significance of this is not clear 65 . Our data are highly suggestive that the active immune response plays a key role in augmenting tumorigenesis in EμMyc/Casp2 −/− mice. A role for caspase-2 in regulating inflammasome signaling and the innate immune response has been previously reported 6 . However, it is also likely that the potential role of caspase-2 in hematopoietic cell differentiation, (e.g., myeloid cell differentiation) 66 may contribute to immunosurveillance mechanisms that determine tumor outcomes.
The comparison of gene expression profiles in EμMyc/ Casp2 −/− tumors to Th-MYCN/Casp2 −/− tumors identified a significant enrichment of genes associated with Rap1 signaling. Rap1 regulates cell adhesion; a process also enriched in the EμMyc/Casp2 −/− tumors, and perhaps indicates a link between these two processes. While Rap1 can promote tumor cell migration, invasion, and metastasis 67 , these features have not been observed in EμMyc/Casp2 −/− mice, indicating possible alternative roles for Rap1 signaling in EμMyc/Casp2 −/− tumors. This may include Rap1 roles in regulating cytoskeletal dynamics and/or cell proliferation, previously associated with loss of caspase-2 function 6,12 . Alternatively, Rap1 signaling plays diverse roles in hematopoietic cells, including lymphocyte activation, migration/trafficking, immunological synapse 67 , B-cell development, and selftolerance 68 . While the specific and diverse functions of Rap1 are context dependent, a role for caspase-2 in regulating Rap1 signaling is highly relevant to lymphoma onset and progression and will be important to further define the contributing roles of Rap1 and caspase-2 in Bcell lymphomagenesis.
In summary, this is the first description of tissue specific and differential gene expression signatures associated with caspase-2 deficiency that can influence tumorigenesis. In particular, our findings indicating that loss of Casp2 impacts neuronal differentiation, is likely a key contributor to delaying neuroblastoma onset and reduced adrenal tumor development in Th-MYCN/Casp2 −/− mice. In contrast, increased immune processes, including hematopoietic cell lineage determination, are potential effectors in EμMyc/ Casp2 −/− lymphomagenesis. In addition, we have identified aberrant expression of several oncogenes and tumor suppressor genes that potentially co-operate with Casp2 loss, to determine tumor outcome. These findings reinforce the distinct contribution of signaling pathways associated with tumor onset and progression in both Th-MYCN/Casp2 −/− and in EμMyc/Casp2 −/− mice.

Tumor samples
Lymphomas from EuMyc and EuMyc/Casp2 −/− mice were taken from cervical lymph nodes 7 . Tumor samples from Th-MYCN and Th-MYCN/Casp2 −/− mice were selected from paraspinal thoracic tumors based on absence or reduced blood vascularization 30 . Four individual tumors per genotype were used for mRNA-seq analysis.

RNA extraction
Total RNA was extracted from frozen tumor tissue using Trizol® reagent (Life Technologies) according to the manufacturer's protocol. RNA was resuspended in diethylpyrocarbonate-(DEPC) treated water and quantified using a NanoDrop1000™ spectrophotometer (Thermo Scientific). RNA samples were quality tested using a 2100 Bioanalyser (Agilent Technologies) to confirm RNA integrity number (RIN) > 7.
mRNA sequencing and bioinformatics analysis RNA-seq and bioinformatics analysis was carried out at the ACRF Cancer Genomics Facility (Centre for Cancer Biology, UniSA). Briefly, polyA + mRNA was enriched using oligo-dT beads and samples were barcoded for pooled sequencing on an Illumina HiSeq 2000 (Agilent Technologies). Short, single-end reads were carried out (1 × 50 bp flow cells) with four samples per lane. This yielded~20-30 million raw reads per sample.
Sequence quality was analyzed using FastQC and the resulting data sets were aligned using the STAR alignment algorithm (https://www.ncbi.nlm.nih.gov/pubmed/ 23104886) to the UCSC mm10/GRCm38 version of the Mus musculus reference genome. Gene counts were obtained from the second strand of the STAR output and the resulting files were transferred to the R statistical environment. The edgeR Bioconductor package 69 was used to fit a linear model and samples were log transformed for normalization and ranking and then expressed as log2-fold change of expression (Log2FC). Data are available in Gene Expression Omnibus (GEO) under accession GSE124051.
Multidimensional scaling plots were generated to explore the relationships in the data before and after processing. (Supplementary Figure S1a and Supplementary Table S1c). All biological replicates were similar to each other except for EμMyc/Casp2 −/− sample #212, which fell outside the EμMyc cluster. Normalization of data was able to correct for this, however a scatterplot matrix analysis of the EμMyc samples indicated that sample #212 was a strong outlier (Supplementary Figure S1a) and did not correlate well with other biological replicates (Supplementary Figure S1b). For this reason, EμMyc/Casp2 −/− sample #212 was removed from further differential expression analysis.

Differential expression analysis of individual genes
Differential expression analysis of individual genes was carried out using edgeR. We divided samples into the two tumor types EµMyc and MYCN groups, and also into two genotype groups Casp2 +/+ and Casp2 −/− . Four-way comparisons in gene expression levels were made, including: Raw counts were extracted for each of these samples and edgeR was used to calculate the DEGs between the two phenotypic groups and expressed as log2-fold change (Log2FC) of expression between conditions being tested. Significant changes in gene expression levels were determined using adjusted Pvalues according to the Benjamini-Hochberg method to control for the false discovery rate (FDR) (Supplementary Tables S1c and S2a-f). For heat-maps, samples were clustered based on gene expression levels and false discovery rate (FDR < 0.05), and the most variable genes plotted (Fig. 1c, d, Supplementary Figures S2d-f, Supplementary Tables S3a-e).

Enrichment analysis of biological processes and pathways
Analysis of Gene Ontology and enriched biological pathways was carried out on DEGs with log2FC > 1 and ≤ −1 (i.e., 2^(1); equivalent to an actual fold change of 2) . (FDR < 0.05), using the Functional Annotation Tool in the Database for Annotation, Visualization and Integrated Discovery (DAVID) Bioinformatics Resources 6.8 (https:// david.ncifcrf.gov/home.jsp) 70,71 . Enriched pathways were analyzed from the KEGG (Kyoto Encyclopedia of Genes and Genomes) and REACTOME database in DAVID. The threshold value of gene counts was set at 3, and the EASE score was set at 0.1 and pathways with Fisher's exact P < 0.05, were considered to be enriched. Where indicated, pathways were also determined using PANTHER 72 (http://www.pantherdb.org/about.jsp) and EnRichR (http://amp.pharm.mssm.edu/Enrichr/) 73,74 databases. Gene Set Enrichment Analysis (GSEA) of pathways and genes was performed using Bioconductor fast pre-ranked GSEA (fgsea) package in R.

Real-time quantitative PCR (qPCR)
Total RNA extracted from frozen tumor tissue, was reverse transcribed using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems) with MultiScribe™ Reverse Transcriptase and oligo-dT primers. Gene expression was normalized to the housekeeping gene β-actin and then expressed as fold change compared to EµMyc or Th-MYCN/Casp2 +/+ tumor samples, using the 2 −ΔΔCt method. See Supplementary  Table S9 for primer sequences.

Statistical analyses
Statistical analysis was carried out using GraphPad Prism software (v 6.0, San Diego, CA, USA). Data are expressed as mean ± SEM and two-tailed unpaired t-test with Welch's correction was used for pair-wise comparisons, unless otherwise indicated. Heat-maps were generated using R software, from natural log2 transformed data. Venn diagrams were generated using online software jvenn (http://jvenn.toulouse.inra.fr/app/example.html) 75 .