Integrative Genome-Scale Analysis Identifies Epigenetic Mechanisms of Transcriptional Deregulation in Unfavorable Neuroblastomas

Neuroblastoma originates from precursor cells of the sympathetic nervous system. It is the most frequent solid tumor of early childhood with a remarkable variation in clinical and biological behavior ranging from spontaneous regression to rapid progression in spite of intensive multimodal chemotherapy. The molecular basis of neuroblastoma pathogenesis is still poorly understood. Genetic alterations seen in high-risk tumors include amplification of the proto-oncogene MYCN, activation of the ALK gene, heterozygous deletions of 1p or 11q and gain of 17q. However, the molecular etiology of a substantial portion of aggressive neuroblastomas remains largely enigmatic. Recurrent somatic mutations are rare (1–3), suggesting that epigenetic mechanisms may drive neuroblastoma development and progression.

DNA methylation, an epigenetic modification via methylation of cytosin carbon 5, is a major mechanism in cell differentiation and neoplastic transformation. Promoter-associated CpG islands have been identified as frequent targets of transcriptionally relevant methylation events. Accumulating evidence, however, suggests that nonpromoter methylation may be also actively involved in gene regulatory processes and is an abundant phenomenon among somatically acquired methylation changes in human cancer (4, 5). In neuroblastoma, candidate-based approaches revealed methylation of several genes including CASP8 (6), which plays an important role in the TNF-related apoptosis pathway. Global approaches based on epigenetic drug-induced re-expression analyses and/or affinity-based capture methods in neuroblastoma cell lines revealed further potential mechanisms and targets of DNA methylation in neuroblastoma (7, 8). Further reports suggested the presence of a CpG island methylator phenotype (CIMP), defined as widespread simultaneous CpG island methylation, which is significantly associated with poor outcome (9, 10). Genome-wide DNA methylation analyses in primary neuroblastomas using methylated DNA immunoprecipitation (MeDIP) or DNA methylation arrays identified methylation events representing new candidate epigenetic biomarkers (11–14). An integrated approach, combining methylome, transcriptome, and genome data from a large cohort of primary neuroblastomas with chromatin modification analyses has, however, not yet been applied.

Here, we analyzed DNA methylation in 105 neuroblastomas using the Illumina 450k methylation array, which covers promoter and gene body sites of 99% of Refseq-annotated genes. These data were complemented by expression profiles and copy number information from the same tumors and combined with primary neuroblastoma- and cell line–derived global histone modification data. We identified functional programs targeted by epigenetic mechanisms in high-risk neuroblastomas and provide evidence for methylation of intragenic enhancers being a mechanism of high-risk–associated transcriptional dysregulation. Furthermore, our data suggest an active contribution of PRC2 components to the downregulation of tumor-suppressive genes that are hypermethylated in high-risk patients. In line with this, a combination of drugs targeting the repressive effect of both DNA methylation and PRC2 efficiently reinduced programs epigenetically impaired in high-risk disease. Our results provide fundamental new insights into epigenetic deregulation of aggressive neuroblastoma and may open new diagnostic and therapeutic avenues for children with high-risk disease.

All patients were enrolled in the German Neuroblastoma Trial (NB97+NB2004) and diagnosed between 1998 and 2011. Treatment was according to the trial guidelines. Informed consent was obtained from the patients' parents and risk stratification was according to the German trial protocol for risk adapted treatment of children with neuroblastoma; our cohort for methylation and expression analysis included 40 low-risk, 9 intermediate-risk, and 56 high-risk patients. Clinical disease stage was assessed according to the International Neuroblastoma Staging System (INSS): stage 1, n = 10; stage 2, n = 9; stage 3, n = 10; stage 4, n = 56; stage 4S, n = 20. Age at diagnosis ranged from 0 to 24.6 years (median age 1 year) and amplified MYCN was seen in 33 neuroblastomas. Chromatin immunoprecipitation DNA-sequencing (ChIP-seq) was done for two additional stage 4 high-risk neuroblastomas, one MYCN-nonamplified (termed HR neuroblastoma I), one MYCN-amplified (HR neuroblastoma II). Tumor cell content was >60% for all samples used.

All cell lines were kindly provided by Larissa Savelyeva (German Cancer Research Center) in 2012–2015 and authenticated via multiplex-FISH karyotyping and short tandem repeat DNA typing at the DSMZ (German Collection of Microorganisms and Cell Cultures). Culture conditions and treatment with all-trans retinoic acid (ATRA) were as described previously (15). Be(2)-C and IMR5-75 were treated with 5-aza-2′-deoxycytidine (DAC; 1 μmol/L, 72 hours) and EPZ-6438 (1 μmol/L 96 hours) in combination or alone. DAC was resubstituted every 24 hours. Controls were treated with solvent (DMSO). For modulation of MYCN levels, stable neuroblastoma cell models were used where MYCN can be either up- or downregulated upon addition of tetracycline [SH-SY5Y-MYCN and IMR5-75-shMYCN (small hairpin RNA targeting MYCN), respectively; refs. 16, 17]. IMR5-75-shMYCN with and without MYCN knockdown (±Tet) were synchronized via thymidine (2 mmol/L) block for 18 hours before release by wash out.

Expression profiles were generated using customized 4 × 44 K oligonucleotide microarrays (Agilent Technologies) and genome-wide DNA methylation was assessed using Infinium HumanMethylation450 (450k) BeadChips (Illumina). See Supplementary Materials and Methods for details.

Genomic annotation of CpGs was done using the assignGenomeAnnotation program of the HOMER tool suite (http://homer.salk.edu/homer). GpGs with gene context annotations “EXON”, “INTRON”, “3UTR” and “TTS” (transcription termination site) were defined as intragenic (gene body), whereas CpGs with annotations “PROMOTER-TSS” and “5UTR” were defined as 5′-associated. For association of DNA methylation and gene expression data, only CpGs with gene context annotation were considered. MassARRAY analysis (Sequenom) at two selected CpGs within CDKN2D revealed strong correlation with methylation levels assessed by 450k arrays [Pearson correlation coefficients: 0.85 (95% confidence interval (CI), 0.79–0.90) and 0.89 (95% CI, 0.84–0.92; ref. 18].

Copy number alterations were assessed in the 450k Infinium array data using a previously described custom approach integrating both methylated and unmethylated signals (19). Resulting profiles were manually curated and overall genomic patterns used for patient risk stratification according to Janoueix-Lerosey and colleagues (20).

ChIP-seq was performed as described previously (21) with changes as described in Supplementary Materials and Methods. RNA-seq was performed as described in Supplementary Materials and Methods.

Patients were clustered on the basis of DNA methylation and resulting patient clusters were analyzed for MYC(N) activity and survival probability as described in Supplementary Materials and Methods. To identify correlations between DNA methylation and gene expression we used maximally selected Wilcoxon rank-sum statistics estimating CpG methylation cut-off points that separate patient groups with differential expression of the corresponding gene (22). Oligonucleotide probes annotated to the same gene were analyzed separately. To identify CpGs whose expression-associated methylation is associated with patient subgroups of differential disease risk, ORs for high-risk disease were estimated. To identify GpGs whose methylation levels are robustly associated with both expression of the corresponding gene and high-risk disease, the following filters were applied: (i) P value ≤0.05 for maximally selected Wilcoxon rank-sum statistics for methylation cut-off point selection, (ii) minimum 2-fold change in expression between cut-off point separated patient subgroups, (iii) minimum mean methylation difference between cut-off point separated patient subgroups: 0.1 (β-values), and (iv) P value ≤0.05 for testing the null hypothesis that the OR estimating the association of methylation separated patient subgroups with high-risk disease is 1 (Fisher exact test). Gene sets associated with the identified CpGs were investigated as described in Supplementary Materials and Methods. To identify differential DNA methylation in drug-treated cell lines, after filtering out CpGs with a delta β-value <0.1 across all samples, a linear regression model was used and significant transcripts were selected on the basis of the moderated t statistic (23). Benjamini–Hochberg procedure was applied to estimate the false discovery rate (FDR). A low-stringency P value threshold (FDR < 0.2) in combination with ranking by M-value fold changes (corresponding to delta β-values > 0.1) was used to identify differential CpG methylation (24). Statistical analyses were performed using the R/Bioconductor software environment.

ChIP-seq single-end reads were aligned to the hg19 genome using Bowtie and only uniquely aligned reads were kept. BAM-Files of aligned reads were further processed using the deepTools suite. Input files were subtracted from the treatment files using the bamCompare tool, applying the SES method for normalization of signal to noise. Resulting signals were normalized to an average 1× coverage to produce signal (bigWig) files. H3K4me1, H3K27ac, and H3K4me3 peaks were called using the MACS 1.4 tool using default parameters. Enhancers were defined as overlapping H3K4me1 and H3K27ac peaks with a minimal distance of 2 kb to the closest H3K4me3 peak, a criterion imposed to rule out selection of (unannotated) promoters. We tested for enriched colocalization of CpG subsets with enhancer elements or ENCODE-derived ChIP-seq peaks using hypergeometric test and all 450k array-represented CpGs with corresponding gene context annotation as background population. The HOMER tool suite (version 4.6) was used for H3K27me3 bedgraph file generation. Differential H3K27me3 coverage in treated versus nontreated cells was estimated by SICER (25) using a window size of 400 and a gap size of 400. Candidate H3K27me3 islands with an FDR < 0.01 were considered as significantly differentially covered and were annotated to the nearest transcription start site.

mRNA-seq raw data were mapped to hg19 by Tophat (2.0.5). Number of reads per gene was calculated using HTSeq (0.5.3p9). Normalization was done using DESeq2 (R/Bioconductor). Fold changes were calculated from normalized data. A Wald test (DESeq2) was used to test for differential expression in treated versus nontreated cell lines and P values were adjusted by Benjamini–Hochberg procedure. An adjusted P value < 0.05 was considered to indicate differential expression. For time course data in synchronized IMR5-75 shMYCN cells, a likelihood ratio test (DESeq2) was used. RNA-seq expression data from 498 primary neuroblastomas (26) were analyzed using the R2 Genomics Analysis and Visualization Platform (http://r2.amc.nl).

RNA-seq, ChIP-seq, as well as expression and methylation array data were deposited in Gene Expression Omnibus (GEO) under accession numbers GSE73518, GSE80197, GSE80397, GSE79859, GSE80243, and GSE80445. ChIP-seq data from primary tumors were deposited at the DKFZ data management platform and may be accessed upon request. Public expression data submitted under GEO accessions GSE35218 and GSE62564 as well as ArrayExpress accession E-MTAB- 2691 were used. UCSC accession IDs for ENCODE ChIP-seq tracks from the neuroblastoma cell line SK-N-SH were: wgEncodeEH003373 (Mx1), wgEncodeEH003393 (p300, SC-584), wgEncodeEH003227 (TCF12), wgEncodeEH003249 (GATA3), wgEncodeEH003302 (RXRA), wgEncodeEH003237 (NFIC), wgEncodeEH003226 (JunD), wgEncodeEH003297 (FOSL2), wgEncodeEH003243 (FOXM1), wgEncodeEH003298 (MEF2A), wgEncodeEH003286 (TEAD4), wgEncodeEH003270 (p300, SC-585), wgEncodeEH003261 (ZBTB33), wgEncodeEH003242 (ELF1), wgEncodeEH002269 (NRSF), wgEncodeEH003375 (RFX5), wgEncodeEH003300 (Pbx3), wgEncodeEH003376 (Rad21), wgEncodeEH003377 (SMC3), wgEncodeEH003269 (USF-1), wgEncodeEH003371 (CTCF), wgEncodeEH003374 (Nrf1), wgEncodeEH003228 (YY1), wgEncodeEH003248 (GABP), wgEncodeEH002271 (Sin3Ak-20), wgEncodeEH002270 (Pol2), wgEncodeEH002301 (TAF1) and wgEncodeEH003299 (Max). UCSC accession IDs for ENCODE ChIP-seq tracks from the neuroblastoma cell line SH-SY5Y were: wgEncodeEH002031 (GATA3, SC269) and wgEncodeEH001770 (GATA2).

Genome-wide DNA methylation patterns of 105 neuroblastomas were assessed using the Illumina 450k methylation array interrogating > 485,000 methylation sites per sample at single-nucleotide resolution. Hierarchical clustering based on the 1,000 most variable probes identified two distinct patient clusters, termed cluster 1 and cluster 2 (Fig. 1A and Supplementary Fig. S1). DNA methylation patterns suggested subgroups in both clusters, such as 2s, a large group within cluster 2 displaying a remarkably homogeneous methylation pattern. Subgroup 2s was strongly enriched with low-risk patients, as stratified according to the German NB2004 trial criteria. Subgroup 2s patients were diagnosed at <1.5 years of age and had mainly localized (INSS 1-3) or 4s tumors that lacked amplified MYCN. Segmental chromosomal aberrations, as assessed by interpretation of copy number data derived from 450k arrays, were rare in subgroup 2s. When overall genomic alteration pattern was used to stratify into risk groups (20), 32 of the 34 patients with tumors harboring only numerical aberrations (Type A, associated with good prognosis) mapped to subgroup 2s. The remaining patients of cluster 2 had MYCN-nonamplified tumors, similar to patients in subgroup 2s, but were otherwise characterized by high-risk disease and enriched for other variables associated with poor patient outcome, including age at diagnosis >1.5 years, stage 4 disease, 11q deletions, or other segmental chromosomal alterations (poor prognosis types B+D ≙ segmental alterations and MYCN-nonamplified). Of the 40 patients in cluster 1, 37 were high-risk (NB2004), and stage 4 disease and an age at diagnosis of >1.5 years were strongly enriched. All MYCN-amplified tumors mapped to cluster 1, being in line with an impact of MYCN on DNA methylation patterns. MYCN-nonamplified cluster 1 tumors had significantly higher MYC(N) target gene activity (defined in Supplementary Materials and Methods) than cluster 2 tumors (Fig. 1B), which might be partially due to c-MYC activity as indicated by elevated MYC mRNA levels in a subgroup of the nonamplified cluster 1 tumors (Fig. 1A). This may suggest that MYC(N) deregulation affects DNA methylation even in the absence of genomic amplification. Testing the association between DNA methylation-based subgroups and patient survival revealed poor outcome for cluster 1 patients and intermediate outcome for cluster 2 patients (Fig. 1C) with excellent outcome for subgroup 2s patients (Supplementary Fig. S2). Multivariate survival analysis suggested methylation-based clustering as an independent predictor of outcome (P = 0.036), with variance inflation factors for MYCN status and methylation cluster indicating a moderate degree of collinearity (Supplementary Table S1), probably due to a correlation between these two predictors (Spearman ρ = 0.86). Together, global DNA methylation can identify groups of neuroblastoma patients, which are strongly associated with differential outcome and either presence or absence of amplified MYCN.

Previous analyses suggest that CpG island methylation in the gene body of protocadherin beta family (PCDHB) members is closely associated with the presence of a high-risk-associated neuroblastoma CIMP (9, 10). To investigate whether the global methylation patterns identified in our study are synonymous with this PCDHB-defined CIMP, we performed hierarchical clustering of the 105 patients using only 450k data representing gene body methylation of PCDHB family members 2–18 (119 CpGs). Of the two clusters identified, one was defined as CIMP-positive by PCDHB hypermethylation, which was corroborated by hypermethylation of two additional CIMP marker genes, CYP26C1 and HLP (Fig. 2; refs. 9, 10). The CIMP cluster contained 32 predominantly high-risk patients and was strongly associated with variables of unfavorable tumor biology, including amplified MYCN, higher age at diagnosis and advanced tumor stage. While not all patients with poor outcome were detected, the majority of patients that died from disease mapped to the CIMP cluster. Of the 32 CIMP patients, the majority mapped to global methylation-defined patient groups associated with high-risk (Fig. 1), supporting that PCDHB family methylation status may help estimate global methylation patterns associated with high-risk neuroblastoma biology.

To estimate the effect of differential DNA methylation on gene expression in neuroblastomas, we integrated 450k methylome data with mRNA expression profiles derived from customized 44 k oligonucleotide microarrays (Agilent Technologies) of the same 105 tumors. Maximally selected Wilcoxon rank-sum statistics were used to estimate CpG methylation cut-off points separating patient groups with differential expression of coannotated genes. Methylation levels of 4,074 CpGs were significantly associated with expression of the corresponding genes (P < 0.05). Of these, hypermethylation of 1,462 and hypomethylation of 1,283 CpGs were significantly associated with high-risk disease (P ≤ 0.05, Supplementary Table S2). Four CpG categories were identified within this high-risk (HR) disease-associated methylation/expression (Fig. 3A; Supplementary Table S2): (i) 860 CpGs representing 341 genes exhibited hypermethylation with downregulation of the corresponding gene (HyperDownHR), (ii) 396 CpGs representing 167 genes exhibited hypomethylation with upregulation of the corresponding gene (HypoUpHR), (iii) 602 CpGs representing 178 genes exhibited hypermethylation with upregulation of the corresponding gene (HyperUpHR), and (iv) 887 CpGs representing 273 genes exhibited hypomethylation with downregulation of the corresponding gene (HypoDownHR).

CpGs whose methylation was negatively associated with gene expression (HyperDownHR and HypoUpHR) mapped closer to the transcription start site (TSS; median distance: +2.9 kb and +3.4 kb, respectively) compared with CpGs with a positive association (HyperUpHR and HypoDownHR; median distance to TSS: +15.5 kb and +38.8 kb, respectively; Fig. 3B). With respect to gene context annotation, CpGs whose methylation was negatively associated with expression of the corresponding gene had a higher fraction of TSS and 5′-UTR-annotated CpGs (HyperDownHR and HypoUpHR; 27% and 28%, respectively) compared with CpGs whose methylation was positively associated with expression of the corresponding gene (HyperUpHR and HypoDownHR; 13% and 8%, respectively). Conversely, CpG categories with a positive association between methylation and expression had a higher fraction of gene body-annotated CpGs (HyperUpHR and HypoDownHR; 87% and 92%, respectively) compared with CpG categories with a negative association (HyperDownHR and HypoUpHR; 73% and 72%, respectively; P < 0.001, Fig. 3A). These results are in line with the established role of promoter/TSS-associated DNA methylation in transcriptional repression and agree with previous genome-wide epigenomic studies linking gene body methylation to active transcription (5). Together, we identified subgroups of CpGs whose methylation is significantly associated with expression of corresponding genes in neuroblastomas of differential risk groups. The nature of this association is strongly dependent on gene context and relative distance to the TSS.

Despite gene body methylation being largely positively associated with expression, a substantial fraction of negatively associated CpGs (HyperDownHR and HypoUpHR) mapped to the gene body. To investigate whether this negative association is due to DNA methylation targeting specific intragenic regulatory DNA elements, we mined ChIP-seq data from neuroblastoma cell lines deposited in the ENCODE (Encyclopedia of DNA Elements, UCSC) database. These analyses revealed a highly significant enrichment of intragenic HyperDownHR and HypoUpHR CpGs with sites bound by enhancer-associated protein p300 in the SK-N-SH neuroblastoma cell line (Fig. 4A; P < 0.001), indicating that methylation of intragenic enhancers might be relevant for the negative association between gene body methylation and gene expression seen in these CpG subgroups. To investigate the role of intragenic enhancers in additional neuroblastoma models, we used ChIP-seq for defining enhancers (H3K4me1/H3K27ac positive, H3K4me3 negative) in two stage 4 primary high-risk neuroblastomas and three high-risk neuroblastoma-derived cells lines, Be(2)-C, SH-SY5Y and IMR5-75. Intragenic CpGs whose methylation was negatively associated with gene expression (HyperDownHR and HypoUpHR) were significantly enriched with enhancer elements in both cell lines and primary tumors (Fig. 4B; all P < 0.02, except HypoUpHR CpGs with primary HR NB I enhancers P = 0.1). These data suggest that methylation of intragenic enhancer elements may regulate expression of genes relevant for high-risk neuroblastoma development, highlighting a role of nonpromoter DNA methylation in neuroblastoma pathogenesis.

To elucidate the biology of genetic programs regulated by promoter and/or enhancer methylation in neuroblastomas of differential risk, we further investigated the HypoUpHR and HyperDownHR gene groups. Among genes that where highly expressed in high-risk neuroblastomas with concomitant CpG demethylation (HypoUpHR), we identified several genes previously described to be involved in the biology of aggressive neuroblastoma including CXCR4 (27), GAL (28), LRRN1 (29), ODC1 (30), TWIST1 (31), and WHSC1 (32). Notably, the HypoUpHR genes included four members (DDX43, PRAME, TEX14, TMEM108) of the Cancer/Testis Antigens and genes for two proteins previously suggested as neuroblastoma antigens or therapy targets, NEK2 (33) and NPY (34). This suggests that aberrant demethylation in aggressive neuroblastomas is associated with (i) upregulation of genes contributing to neuroblastoma progression and (ii) induction of possible targets for immunotherapeutic intervention.

The 341 genes downregulated in high-risk patients with concomitant CpG hypermethylation (HyperDownHR) represented the largest group of genes among all possible expression/methylation association combinations. Some of these HyperDownHR genes have been previously implicated in neuroblastoma-relevant aberrant methylation, including ABCB1, CACNA1G, CD44, DUSP23, PRDM2, RBP1, SFRP1 (reviewed in ref. 35), CHD5 (36), and NTRK1 (37). Four genes (KRT19, PRPH, CNR1, QPCT) were part of an eight-gene DNA methylation-based prognostic biomarker (8), substantiating that the prognostic value of these methylation markers is tightly linked to downregulation of the corresponding genes. Mapping of HyperDownHR genes to the catalog of nonsilent mutations identified via global neuroblastoma sequencing approaches (1–3) revealed an overlap of 85 genes that are likely to be targeted in neuroblastoma by rare somatic mutations and epigenetic silencing (Supplementary Table S3). These include DLC1, a tumor suppressor gene recurrently mutated in neuroblastomas (2, 3) and involved in Rac/Rho signaling, a pathway essential for neuritogenesis that is frequently impaired in high-risk neuroblastomas (3). HyperDownHR genes also included direct targets of neuroblastoma-specific structural aberrations, such as the neurotransmission-associated ASIC2, whose disruption by a constitutional translocation has been shown (38), and RGS5, reported to be downregulated via MYCN-activated miRNAs and affected by focal homozygous deletion (39). Chromosomal distribution analysis of HyperDownHR genes revealed significant enrichment of two cytogenetic regions, 5q31 (P < 0.006), encompassing the neuroblastoma CIMP-predictive PCDHB gene family cluster, and 1p36 (P = 0.018, Table 1), which is frequently deleted in neuroblastoma and other cancers. To provide further insight into molecular processes mediated by HyperDownHR genes, we performed GO term and pathway enrichment analysis using the Database for Annotation, Visualization and integrated Discovery tool (DAVID, Table 1; Supplementary Table S4). Significantly enriched GO terms included terms like “regulation of apoptosis” and “positive regulation of lymphocyte differentiation”, with the latter likely reflecting differential lymphocytic infiltration among neuroblastoma subtypes. A predominant fraction of significantly enriched GO terms, however, related to neuronal differentiation or function (e.g., “neuron projection”, “synapse”, “axon”, “synaptic vesicle”, “nervous system development”). In line with this, the Reactome pathway most strongly overrepresented among HyperDownHR genes was “Synaptic Transmission” (REACT_13685, P < 0.008). Involvement of HyperDownHR genes in neuronal differentiation is further supported by their time-dependent upregulation in Be(2)-C cells induced to differentiate via retinoic acid treatment (Supplementary Fig. S3). Intriguingly, genes occupied by Polycomb repressive complex 2 (PRC2) components EED and SUZ12 and marked by H3K27me3 in embryonic stem cells (ESC; ref. 40) were significantly enriched in the HyperDownHR group (all P < 1 × 10⁻¹⁰, Table 1), indicating that silencing HyperDownHR genes may contribute to cell stemness. Overall, investigation of HyperDownHR genes revealed a strong representation of established tumor suppressors and genes previously found to be impaired in neuroblastomas by genetic events, together with a significant enrichment of ESC PRC2 targets and genes involved in neuronal differentiation.

Enrichment of HyperDownHR genes with ESC PRC2 targets corresponds to the idea that genes silenced in ESCs by Polycomb group proteins are prone to hypermethylation in cancer (41). To estimate whether the PRC2-mediated H3K27 trimethylation of HyperDownHR genes is a transient developmental feature or whether it contributes to silencing of HyperDownHR genes in established neuroblastoma cells, H3K27me3 ChIP-seq was performed in two stage 4 primary high-risk neuroblastomas and two high-risk neuroblastoma-derived cells lines, Be(2)-C and SH-SY5Y. H3K27me3 was strongly enriched at HyperDownHR genes compared with global gene-associated H3K27me3 occupancy in both cell lines and primary tumors (Fig. 5A–D). Integration of H3K27me3 and DNA methylation data (450k array) in Be(2)-C identified several H3K27me3-covered HyperDownHR genes with high methylation levels at putative regulatory regions. Among these, promoter methylation as exemplified by HENMT1 (Supplementary Fig. S4A) was less frequent than methylation of intragenic enhancers as exemplified by SPOCK2 (Supplementary Fig. S4B) and SLC18A2 (Supplementary Fig. S4C). Our data suggest that PRC2 activity contributes to the ongoing repression of hypermethylated genes in established high-risk neuroblastomas, which needs to be considered for therapeutic approaches aiming at derepression of HyperDownHR genes.

c-MYC has been shown to induce expression of PRC2 components by various mechanisms including direct transcriptional activation (42, 43), and expression of the PRC2 component gene EZH2 is elevated in high-risk neuroblastomas (44). To investigate the influence of amplified MYCN on PRC2 component expression in primary neuroblastomas, we analyzed transcriptomes from 498 tumors. Unsupervised clustering based on expression of the PRC2 core genes EED, EZH2, SUZ12 and RBBP7 (RbAp46) identified a group of tumors characterized by PRC2 hyperactivity and strong prevalence of amplified MYCN (Supplementary Fig. S5A). In line with this, EED, EZH2, and RBBP7 were significantly higher expressed in MYCN-amplified versus -nonamplified neuroblastomas (all P < 0.001, Supplementary Fig. S5B). Analyzing the impact of MYCN on PRC2 component expression in neuroblastoma cells, we found significantly higher expression of EED, EZH2, and RBBP7 in synchronized MYCN-amplified IMR5-75 S-phase cells compared with MYCN-knockdown IMR5-75 cells at equivalent time points (all P < 0.05, Supplementary Fig. S5C). To estimate whether a MYCN-associated PRC2 component induction might contribute to the observed H3K27me3 enrichment at HyperDownHR genes, ChIP-seq was performed in neuroblastoma cells upon MYCN modulation. In SH-SY5Y cells, MYCN induction led to increased H3K27me3 coverage at HyperDownHR genes (Fig. 6A), which was accompanied by a time-dependent repression of HyperDownHR gene expression (Fig. 6B; Supplementary Table S5). The remaining global transcriptome was affected to a substantially lesser extend as demonstrated by strong enrichment of HyperDownHR genes among genes downregulated upon MYCN induction (Fig. 6B, bottom). In line with this, MYCN knockdown in IMR5-75 cells led to decreased H3K27me3 coverage at HyperDownHR genes, which was accompanied by a specific induction of HyperDownHR gene expression (Supplementary Fig. S6 and Fig. 6C; Supplementary Table S5). The activation of PRC2 components by MYCN, together with our observation that global DNA methylation is strongly associated with genomic MYCN status in primary neuroblastomas, indicates that MYCN may prime a subset of genes for epigenetic silencing in high-risk neuroblastomas.

As reinducing genetic programs silenced in high-risk neuroblastoma is likely to be therapeutically beneficial, we investigated the potential of epigenetically active drugs to derepress HyperDownHR genes in the neuroblastoma cell lines Be(2)-C and IMR5-75. In both cell lines, treatment with the 5-aza-2′-deoxycytidine (DAC) DNA-demethylating agent led to a preferential induction of HyperDownHR genes compared with global gene induction (Fig. 6D; Supplementary Table S6). The overall potency of DAC to significantly induce gene transcription, however, was substantially higher in IMR5-75 compared with Be(2)-C (induced gene fraction HyperDownHR 25.2% vs. 5%, global 7.9% vs. 1.4%, respectively). Preferential induction of HyperDownHR genes was also seen upon treatment with EPZ-6438, a specific inhibitor of the H3K27 methyltransferase EZH2. With EPZ-6438, the overall effect on gene transcription was more pronounced in Be(2)-C compared with IMR5-75 (induced gene fraction HyperDownHR 11.7% vs. 6.5%, global 2.7% vs. 0.9%, respectively). Most efficient and preferential HyperDownHR reinduction (P < 0.001 in both cell lines), however, was achieved by combined treatment with DAC and EPZ-6438, supporting the notion that DNA methylation and H3K27 trimethylation cooperate to repress HyperDownHR genes (fraction of HyperDownHR genes induced: 24.6% vs. 5.7% global in Be(2)-C, 34.3% vs. 10.9% global in IMR5-75). To monitor the epigenetic changes associated with this induction, we investigated gene-associated H3K27me3 coverage (ChIP-seq) and DNA methylation (450k arrays) in Be(2)-C and IMR5-75 cells after DAC/EPZ-6438 combination treatment (Fig. 6E; Supplementary Table S7). Of 81 HyperDownHR genes induced in Be(2)-C, 47 (58%) lost either K27me3 methylation, DNA methylation or both. Of 114 HyperDownHR genes induced in IMR5-75, 51 (45%) lost either K27me3 methylation, DNA methylation, or both. Loss of both repressive marks was relatively rare with three HyperDownHR genes in Be(2)-C and two HyperDownHR genes in IMR5-75 (exemplified by DPP6 in Be(2)-C, Supplementary Fig. S7). To further investigate the idea of preferential HyperDownHR induction by combination treatment in a larger neuroblastoma panel, we reanalyzed genome-wide expression data assessed in 17 neuroblastoma cell lines treated with DAC and trichostatin A (TSA; ref. 45), a histone deacetylase (HDAC) inhibitor previously shown to restore expression of PRC2-repressed genes (46). After genome-wide testing for differential expression upon treatment across all 17 cell lines, gene-set enrichment analysis revealed significant enrichment of HyperDownHR genes among DAC/TSA-induced genes (P = 0.039, Supplementary Fig. S8, Supplementary Table S8). Together, these data suggest that genetic programs silenced in high-risk neuroblastoma can be efficiently and preferentially derepressed in a broad range of neuroblastoma subtypes by combining epigenetic drugs that target both PRC2- and DNA methylation-mediated repression.

In an integrative approach, we investigated methylomes, transcriptomes, and copy number aberrations of 105 neuroblastomas together with tumor- and cell line-derived chromatin modification data. Clustering based on DNA methylation profiles identified patient clusters and subgroups strongly associated with key clinical parameters and specific genetic alterations in the tumors, which suggests that deregulated methylation is a fundamental feature of high-risk neuroblastomas. Multivariate survival analysis suggested that 450k methylome data contain prognostic information that could complement current risk stratification. Further studies in an extended patient cohort need to refine (i) which loci are most informative as predictive biomarkers, (ii) to which extend DNA methylation can complement the correlated predictive factor MYCN status, and (iii) which patient subgroups would benefit most from DNA methylation-based stratification.

Clustering based on methylation of PCDHB family members largely reflected global DNA methylation, suggesting that PCDHB methylation could function as a marker for genome-wide DNA methylation patterns associated with high-risk neuroblastoma. Detection platform differences may contribute to the observed limitation of PCDHB methylation-based clustering to detect all patients with poor outcome, as the 450k probe design only incompletely covers the CIMP marker sites previously identified (9, 10). PCDHB family methylation has, until now, been generally thought to be a surrogate marker for synchronized methylation events acting on other genes that contribute to neuroblastoma development, instead of having a direct regulatory role (47). Of note, we identified hypermethylation of several PCDHB family members to be associated with downregulation of their transcripts in high-risk neuroblastomas, which argues for a more direct role of PCDHB impairment in neuroblastoma pathogenesis.

As expected, we frequently observed methylation in the promoter and surrounding regions that was negatively associated with expression. Positive correlation between DNA methylation and expression was predominantly associated with the gene body, a phenomenon detected in various biological contexts, including development, differentiation and cancer (5). A surprisingly large fraction of negatively associated DNA methylation sites, however, also mapped to the gene body. This might be partially due to methylation of promoter downstream correlating regions (pdCR), which were recently shown to extend up to tens of kilobases downstream of the transcriptional start site (4), a gene context not annotated as promoter-associated in our analysis. However, we revealed another likely regulatory mechanism to be transcriptional inhibition via intragenic enhancer methylation. A growing body of evidence suggests that enhancer methylation may be a pivotal element in both physiologic transcriptional regulation and disease-associated dysregulation, and, strikingly, expression variation in a large fraction of genes appears to be exclusively associated with enhancer methylation while promoter methylation remains unaffected (5, 48). Modulation of intragenic enhancer methylation has also been described to be the dominant methylome change during cell differentiation (5). We conclude that the high prevalence of expression-associated intragenic enhancer methylation being differential across neuroblastomas of varying biology suggests that this epigenetic event accounts for a significant proportion of subtype-specific intertumor diversity and contributes to high-risk neuroblastoma development.

Our analysis identified a catalog of genes that are strong candidates for being transcriptionally deregulated by aberrant methylation in high-risk neuroblastomas. The substantial overlap with genes that are targeted by genetic events in neuroblastomas indicates that epigenetic and genetic mechanisms converge partially on the same targets to contribute to neuroblastoma pathogenesis. Given the low incidence of most somatic mutations identified by genome-wide sequence analyses in neuroblastomas (1–3), our data may be used to prioritize candidates for future analyses.

From the subgroup of genes hypomethylated and upregulated in high-risk neuroblastomas (HypoUpHR), several were previously shown to be activated in aggressive neuroblastomas, but little was known about the mechanisms of their deregulation. Our data suggest that these mechanisms include aberrant demethylation during high-risk neuroblastoma development. HypoUpHR genes should be strongly enriched for potential therapeutic targets. Of special interest in this context is the high-risk–specific hypomethylation and overexpression of DDX43, PRAME, TEX14, and TMEM108, of which PRAME has been previously identified as a candidate for therapeutic intervention in neuroblastoma (49). These genes all encode Cancer/Testis Antigens, proteins whose expression is typically restricted to the testis but that are aberrantly expressed in many cancers, constituting attractive targets for biomarkers and immunotherapy. Our results support the idea that DNA hypomethylation is a driver of Cancer/Testis Antigen derepression in neuroblastoma as has been previously proposed for other tumor entities (50).

The most prominent high-risk–associated phenomenon was hypermethylation in combination with transcriptional downregulation of genes (HyperDownHR). A subset of this gene group has been previously reported to mediate tumor-suppressive functions and/or be targeted in neuroblastomas by methylation or genetic events including point mutation and structural aberrations (1–3, 8, 35–39). Positional gene enrichment analysis revealed an overrepresentation of 1p36 genes, which included the well-established dosage-dependent neuroblastoma suppressor candidates CAMTA1, CHD5, and KIF1B (51), suggesting that DNA methylation may cooperate with genomic deletion to shift gene product dosage to levels below those needed for oncosuppression. Strong enrichment of neuron-specific pathways and GO terms, together with the observed induction of HyperDownHR genes in differentiating neuroblastoma cells, indicates that these genes contribute to neuronal development and function. Together, this suggests that hypermethylation and concomitant downregulation of the identified gene group provides a selective advantage during high-risk neuroblastoma development by facilitating circumvention of tumor suppressive and prodifferentiation programs.

A subset of genes targeted by PRC2 in ESCs and undifferentiated early cell compartments is thought to be transferred to a tighter repressive state by methylation during tumor development. This gene subset, previously termed “DNA methylation module”, often comprises developmental regulators that may differ across tumor entities, and whose silencing may contribute to the stem-like state of cancer cells (52). Enrichment of both PRC2 ESC targets and neuronal development genes among HyperDownHR genes indicates that we identified the neuroblastoma-specific DNA hypermethylation module, which defines stem cell qualities such as unlimited self-renewal and evasion from differentiation. We identified overrepresentation of H3K27me3 at HyperDownHR genes in established neuroblastoma cells, suggesting that PRC2 continues to contribute to silencing of these genes, a phenomenon of great therapeutic implication also recently acknowledged in other cancer entities (53, 54). In line with this, combined treatment targeting both DNA methylation and PRC2 (EPZ-6438/DAC) efficiently and preferentially reinduced this gene group, while the potency of EPZ-6438 and DAC alone was lower and diverged in the tested cell lines. This suggests that the efficiency of such combined treatment is robust with respect to the (epi)genetic background of cells, which is further supported by significant reexpression of HyperDownHR genes in a diverse panel of 17 neuroblastoma cell lines by TSA/DAC combination. A large fraction of EPZ-6438/DAC-induced genes lost H3K27me3 coverage or DNA methylation but loss of both marks was not common, which may suggest that on the level of individual genes, either H3K27me3 or DNA methylation is the main repressive mechanism. However, on the level of cancer-silenced gene groups, targeting both modifications seems to be required for efficient induction of the genetic program. We propose that combining drugs targeting PRC2- and DNA methylation–mediated repression could be a rational strategy for high-risk neuroblastoma therapy. This concept is strengthened by recent data showing that a modification by both of these marks is strongly overrepresented in cancer cells compared with their normal counterparts and, thus, may serve as a cancer-specific target (54).

c-MYC stimulates PRC2 components in ESCs to keep bivalent genes silent, maintaining ESCs in an undifferentiated state (43), and it was shown to activate the PRC2 component EZH2 in cancer cells (42). Our data indicate that MYCN plays a similar role in neuroblastoma cells by directly or indirectly activating the expression of PRC2 components, which adds to MYCN's function as recruiter of EZH2 to the promoter of target genes (55). Considering the observed increase in H3K27me3 coverage at HyperDownHR genes in MYCN-overexpressing neuroblastoma cells and the established functional crosstalk between PRC2 activity and DNA methylation, our data point to a model where deregulated MYCN guides DNA methylation via setting PRC2-defined pre-marks for permanent silencing of a set of genes whose expression is incompatible with progression to high-risk neuroblastoma. Together, our integrative analysis has uncovered a wealth of epigenetic alterations that enhance the understanding of neuroblastoma heterogeneity and shed light on potential regulatory mechanisms involved in neuroblastoma development and progression. The strikingly subtype-specific methylation patterns are promising points of departure for prognostic purposes whereas the deregulated genetic programs and particularly the epigenetic mechanisms involved in their deregulation are attractive targets for therapeutic intervention.

No potential conflicts of interest were disclosed.

Conception and design: K.-O. Henrich, S. Bender, E. Bell, S.M. Pfister, F. Westermann

Development of methodology: K.-O. Henrich, D. Dreidax, M. Gartlgruber, C. Plass, A. Benner, F. Westermann

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): K.-O. Henrich, S. Bender, D. Dreidax, M. Gartlgruber, M. Parzonka, L. Wehrmann, M. Fischer, D.J. Duffy, E. Bell, A. Torkov, P. Schmezer, F. Westermann

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): K.-O. Henrich, S. Bender, M. Saadati, M. Gartlgruber, C. Shao, C. Herrmann, M. Wiesenfarth, M. Parzonka, D.J. Duffy, P. Schmezer, C. Plass, T. Höfer, A. Benner, S.M. Pfister, F. Westermann

Writing, review, and/or revision of the manuscript: K.-O. Henrich, S. Bender, M. Saadati, D. Dreidax, M. Gartlgruber, M. Parzonka, M. Fischer, C. Plass, A. Benner, S.M. Pfister, F. Westermann

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Gartlgruber, M. Parzonka, F. Westermann

Study supervision: K.-O. Henrich, S. Bender, F. Westermann

We thank the tumor bank team of the University Hospital of Cologne (Germany) for providing primary tumor samples, Steffen Bannert, Jochen Kreth, Elisa Hess, and Young-Gyu Park for excellent technical assistance and Kathy Astrahantseff for critical reading of the manuscript.

This work was supported by the German Cancer Aid (grant no. 110122 to F. Westermann), the German Ministry of Science and Education (BMBF) as part of the e:Med initiative (grant no. 01ZX1307D to M. Fischer and F. Westermann), the BMBF MYC-NET (grant no. 0316076A to F. Westermann, S.M. Pfister, and S. Bender), the European Union grant no. 259348 (F. Westermann), the German Cancer Research Center (DKFZ) intramural program for interaction projects (F. Westermann and S.M. Pfister), the DKFZ – Heidelberg Center for Personalized Oncology (HIPO) & National Center for Tumor Diseases (NCT) Precision Oncology Program (S.M. Pfister), and intramural funding through NCT 3.0: ENHANCE – Enhancers and non-coding (epi-)mutations (F. Westermann, C. Herrmann, C. Plass, S.M. Pfister).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.