Aberrant DNA methylation distorts developmental trajectories in atypical teratoid/rhabdoid tumors

Abstract

Atypical teratoid/rhabdoid tumors (AT/RTs) are pediatric brain tumors known for their aggressiveness and aberrant but still unresolved epigenetic regulation. To better understand their malignancy, we investigated how AT/RT-specific DNA hypermethylation was associated with gene expression and altered transcription factor binding and how it is linked to upstream regulation. Medulloblastomas, choroid plexus tumors, pluripotent stem cells, and fetal brain were used as references. A part of the genomic regions, which were hypermethylated in AT/RTs similarly as in pluripotent stem cells and demethylated in the fetal brain, were targeted by neural transcriptional regulators. AT/RT-unique DNA hypermethylation was associated with polycomb repressive complex 2 and linked to suppressed genes with a role in neural development and tumorigenesis. Activity of the several NEUROG/NEUROD pioneer factors, which are unable to bind to methylated DNA, was compromised via the suppressed expression or DNA hypermethylation of their target sites, which was also experimentally validated for NEUROD1 in medulloblastomas and AT/RT samples. These results highlight and characterize the role of DNA hypermethylation in AT/RT malignancy and halted neural cell differentiation.

Figure 1.. Characterization of DNA methylation differences among AT/RTs, MBs, and PLEXs reveals AT/RT hypermethylation across all AT/RT subclasses and the genomic regions affected by large-scale DNA methylation changes.

(A) Illustration of our data analysis and integration approach. The number of samples in each cohort is shown on the left. Data were used to call differentially methylated regions (DMRs) and differentially expressed (DE) genes. Data from the Gene Transcription Regulation Database provided transcription factor DNA binding information. DNA methylation data from pluripotent stem cells and normal samples were used as references and to study normal neural cell differentiation. (B) Tumor types are separated into tumor subgroups (omitted from Capper et al [2018]) based on DNA methylation, when the 10,000 most variable regions measured in both i450k and RRBS data (see the Materials and Methods section) were used for the tSNE visualization. RRBS samples are positioned adjacent to i450k samples representing the tumor subgroups that matched their clinical diagnosis. (C) Venn diagrams showing the number of DMRs in each comparison. Tumor type–specific DMRs are marked into the intersecting areas. A higher number of DMRs were detected in the RRBS than in the i450k data. For i450k results, DMRs were filtered using DNA methylation data from normal brain samples. (D) For AT/RTs, larger numbers of hypermethylated than hypomethylated regions were detected in all the comparisons in both i450k and RRBS data. The numbers of DMRs and direction of DNA methylation change for each comparison in both datasets. (E) AT/RT subgroups showed the highest DNA methylation among the pooled DMRs when compared to other tumor types and normal control samples (CONTR). Average DNA methylation of the probes hitting each i450k DMR is visualized as tumor subgroup-wise violin plots. (F) k-Means clustering analysis revealed DMR clusters that are specifically hypermethylated in AT/RTs. The median DNA methylation of the DMRs in each cluster was used to summarize the DNA methylation patterns. None of the DMR clusters showed AT/RT subtype–specific DNA methylation patterns, but there were DNA methylation differences between tumor subtypes within both MBs and PLEXs. (G) Several topologically associating domains were influenced by large-scale DNA methylation differences, especially in the AT/RT-MB comparison. Karyoplot visualizes the topologically associating domains that harbor large-scale DNA methylation differences, that is, several DMRs that were predominantly either hyper- or hypomethylated in the comparison (see the Materials and Methods section). Color indicates a comparison in which a difference was observed.

Figure S1.. Characterization of DMRs.

(A) Venn diagram of i450K DMRs obtained without filtering based on tumor location or normal tissue samples. (B) DMRs in different comparisons show similar patterns, but a higher proportion of DMRs were covering interCGI and intergenic regions in RRBS data than in i450K data. Normalized proportions of DMRs hitting to annotatr built-in CpG island annotations (in orange), and genomic annotations and enhancers (in green) are visualized. The i450K and RRBS results were made comparable by annotating each bp from DMRs separately and normalizing the annotation counts with the total length of DMRs in each comparison. InterCGI means areas not hitting to CpG islands, shelves, or shores, that is, the “open sea.” 1–5 kb contains regions 1–5 kb upstream of the TTS. (C) Distribution of DMRs into autosomal chromosomes. Most of the DMRs are hypermethylated in AT/RTs and hypomethylated in MBs when compared to other tumor types. (D) Methylation distribution in pooled DMRs from RRBS data. (E) Same as in (D) but with subgroup information.

Figure S2.. Heatmaps for k-means clustering.

(A, B, C, D, E, F, G, H, I, J) Corresponding cluster heatmaps for k-means clustering. The methylation class shows the tumor type and subtype.

Figure 2.. Neural cell differentiation–related TFs and key epigenetic regulators are enriched in tumor type–specific DMRs.

(A) Most of the AT/RT-specific DMRs were hypermethylated (99% and 79% in i450K and RRBS data, respectively), whereas hypomethylated DMRs were more commonly observed in MBs (85% and 88% in i450K and RRBS data, respectively) and PLEXs (98% and 71% in i450K and RRBS data, respectively). Four-field plots for each tumor type show the number of hypermethylated and hypomethylated regions with respect to two other tumor types. Tumor-specific DMRs have the DNA methylation change in the same direction when compared to two other tumor types (e.g., hypermethylated in AT/RTs when compared to MBs and to PLEXs). (B) Most of the transcription factors and other transcriptional regulators (jointly referred to as TFs) are specifically enriched in AT/RT-hyper, MB-hyper, or PLEX-hypo DMRs and largely linked to neural differentiation, SWI/SNF, and PRC2. Upper part: upset plot showing the number of enriched TFs for regions that are hypermethylated or hypomethylated in an AT/RT-, MB-, and PLEX-specific manner. Some transcriptional regulators were enriched in several tumor-specific DMR groups. Lower part: the number of TFs in manually annotated function-related theme groups is shown for each upset plot column. The color of the heatmap shows the fraction of TFs in each theme (row). (C) Binding sites of neural TFs measured in brain tumors and other neural samples were enriched in AT/RT-hyper DMRs, whereas those measured in pluripotent stem cells were enriched in MB-hyper DMRs. GTRD TF binding data were categorized based on the measured sample type into the listed subsets (at the bottom of the plot), and the enrichment of TF binding sites in all the DMRs with AT/RT- or MB-specific DNA methylation was calculated for each of the GTRD subsets separately. Category “All*” means all the reported TF binding sites, so the full GTRD data. Results for the most relevant TFs are shown after organizing them into the theme groups listed in Fig 2B. The dot is not marked when a given TF or other regulator is not measured in a given GTRD subset. (D) PRC2 subunits rarely co-localize with neural TFs and other regulators in AT/RT-hypermethylated sites. Heatmap visualization of the enrichment P-value (one-sided Fisher’s exact test) for co-localization. All the adjusted P-values of 0.01 or higher are marked in white. Themes for each TF are annotated on the right-hand side of the heatmap. (E, F) In our CUT&RUN sequencing analysis, the DNA binding sites of NEUROD1 in the MB cell line overlapped regions that are hypermethylated in AT/RTs and hypomethylated in MBs (AT/RT versus MB-hyper), whereas no NEUROD1 binding sites in AT/RT samples overlapped with DMRs. Heatmaps showing the NEUROD1 binding sites located in different types of DMRs in i450K (E) and RRBS (F) data across the analyzed cell lines.

Figure S3.. Annotations for cancer specific DMRs and average methylation in relevant TFBSs.

(A) Genomic annotations for cancer specific DMRs in all tumors for RRBS and i450k data. A higher proportion of MB-hypermethylated DMRs (7.6% and 32% in RRBS and i450k data, respectively) were located in CpG islands, when compared to MB-hypomethylated, AT/RT-hypermethylated, or PLEX-hypomethylated DMRs (1.9–5.9% and 11–12% in RRBS and i450k data, respectively). (B) Binding sites of neural differentiation factors NEUROG2, ASCL1, and PAX7 are more methylated in AT/RTs irrespective of the tumor subtype. Of the TFs linked to histone lysine methylation, EZH2 (involved in histone H3 lysine 27 trimethylation) behaved similarly, but a distinct DNA methylation pattern was observed for MIER1 and EHMT2 (involved in histone H3 lysine 9 methylation) binding sites. Violin plots visualizing the distribution of DNA methylation with selected TF binding sites in all i450k DMRs. The tumor subgroups are presented separately.

Figure S4.. Expression heatmaps of enriched TFs.

(A) Heatmaps visualizing the expression of TFs associated with tumor type–specific hypomethylated DMRs in RNA-seq and microarray data. (B) Same as in (A) but for hypermethylated DMRs.

Figure S5.. Expression boxplots for most relevant enriched TFs.

(A) Neural differentiation–related TFs enriched in DMRs hypermethylated in AT/RTs (AT/RT-hyper) showed differential expression between tumor types. The expression of selected TFs enriched in tumor-specific DMRs is visualized in boxplots (*, **, and *** refer to P < 0.05, P < 0.01, and P < 0.001, respectively). The upper part is for RNA-seq data and the lower for GSE42658 array data. (B) Pluripotency-related TF expression and the expression of selected TFs enriched in tumor-specific DMRs are visualized in boxplots (*, **, and *** refer to P < 0.05, P < 0.01, and P < 0.001, respectively). the left part is for RNA-seq data and the right for GSE42658 array data.

Figure S6.. Enrichment dotplots for all the TFs in categorized GTRD data.

(A) Enrichment of SWI/SNF, H3K27me, pluripotency, and H3K27ac TF binding sites when GTRD data were split based on the sample of origin in the listed categories (bottom). The dot is not marked when the given TF is not measured in the given GTRD category. (B) Same as in (A) but with neural TF binding sites.

Figure S7.. Co-localization heatmaps for different regions.

(A) Co-localization heatmap of TFs with binding sites enriched in regions hypermethylated in AT/RTs. All the adjusted P-values of 0.001 or higher are marked as white. (B) Same as in (A) but for regions hypomethylated in PLEXs. (C) Same as in (A) but for regions hypomethylated in MBs.

Figure S8.. Co-localization heatmap of TFs with binding sites enriched in regions hypermethylated in MBs.

All the adjusted P-values of 0.001 or higher are marked as white.

Figure S9.. CUT&RUN shows NEUROD1 binding mostly in MB cells.

(A) Venn diagram presenting the overlap between NEUROD1 peaks that share the binding site with NEUROD1 GTRD binding site or have NEUROD1 motif present in the peak. (B) Venn diagram presenting overlapping NEUROD1 peaks from different cell lines. (C, D) Heatmaps showing the NEUROD1 binding sites located in different types of DMRs in i450K (C) and RRBS (D) data across the analyzed cell lines. (E) Violin plot showing the average methylation of cell lines in NEUROD1 binding sites from an MB-3021 cell line. (F) Western blot with a NEUROD1 antibody used in CUT&RUN. Beta-tubulin was used as a loading control.

Figure 3.. AT/RTs harbor pluripotent stem cell–like and AT/RT-unique DMRs, which are associated with the DNA binding of relevant transcriptional regulators.

(A) Pluripotent stem cells (PSCs), primary adult brain, and primary fetal brain (FB) are separated from tumor samples based on DNA methylation in tSNE visualization, when the 10,000 most variable regions in i450k data were used for visualization. (B) When using the same set of pooled DMRs as in Fig 1E, the median DNA methylation level of PSCs is most similar to AT/RTs. (C) AT/RT-hyper DMRs were mostly AT/RT-unique or PSC-like, whereas MB DMRs were MB-unique or FB-like. Very few MB-hypermethylated DMRs were associated with large-scale differences in DNA methylation. Tumor type–specific DMRs (Fig 2A) were categorized based on DNA methylation levels in PSC and FB samples. The bar plot on the left shows the number of DMRs in different categories. Annotations show whether DMRs are PSC-like (P), FB-like (F), or unique (different from PSCs and FB) and whether DNA methylation changes during cell differentiation from PSC to FB. The proportion of DMRs in large-scale DNA methylation differences within annotated DMR categories is shown in blue on the right. The number of DMRs is marked in the figure. (D) DMR category–related DNA binding patterns revealed transcriptional regulators (TFs) involved in tumor-unique, normal cell–like, and differentiation-related regulation of DNA methylation. PRC2 subunits were enriched in the AT/RT-unique DMRs, whereas neural TFs were enriched in both AT/RT-unique and PSC-like DMRs with varying enrichment patterns. TF binding site enrichment was calculated separately for each normal cell differentiation–related DMR category (bottom). TFs were organized into the themes listed in Fig 2B. The dot is not marked when a given TF is not measured in a given GTRD category.

Figure S10.. All the categories from PSC and FB comparisons.

(A) Tumor type–specific DMRs (Fig 2A) were categorized based on DNA methylation levels in PSC and FB samples. The bar plot on the left shows the number of DMRs in different categories. Annotations show whether DMRs are PSC-like (P), FB-like (F), or unique (different from PSCs and FB) and whether DNA methylation changes during cell differentiation from PSCs. The proportion of DMRs in large-scale methylation differences within annotated DMR categories is shown in blue on the right. The number of DMRs is marked in the figure. This figure has all the possible groups (compared with Fig 3C). (B) TF binding site enrichment was calculated separately for each DMR category (bottom). TFs were organized into the themes listed in Fig 2B. The dot is not marked when a given TF is not measured in a given GTRD category. This has all the groups shown in (A).

Figure 4.. DNA methylation associated with the differential expression of genes relevant for neural differentiation and oncogenesis.

(A, B, C, D) Differential DNA methylation (DM) was associated with differential gene expression (DE). Gene expression and DNA methylation patterns were studied in four contexts: differential gene expression alone (A) and DE coupled with DM in the genomic neighborhood (±200 kb from the transcription start site [TSS] within the same topologically associating domain) (B), DE coupled with DM in gene-linked enhancer (C), and DE coupled with DM in the gene promoter (2 kb upstream and 500 bp downstream from the TSS) (D). Venn diagrams show the numbers of genes behaving similarly in both sequencing and array data. Differentially expressed genes associated with differential DNA methylation (B, C, D) are called DM-DE genes. Only cases where the sign of DM change was opposite to DE were included in the figure. (E) DM-DE genes in AT/RT comparisons show generally high DNA methylation among AT/RTs. Sample-wise heatmaps show the levels of DNA methylation (average methylation of variable sites) and gene expression. The rightmost heatmap summarizes in which comparison the DM-DE gene was detected, what was the direction of DNA methylation change (hyper/hypo), and the genomic location of the DMR. (F) Expression patterns of selected DM-DE genes. *P < 0.05, **P < 0.01, ***P < 0.001. (G) Hypermethylated DMRs in relevant genes, which are hypermethylated and underexpressed in AT/RTs. CXXC5 and TCEA3 are AT/RT-specifically suppressed DM-DE genes, and NEUROG1, EBF3, and NEUROD2 are DM-DE genes in the AT/RT-MB comparison. Distal DMRs are connected to the TSS via an arch. Oncoprint indicates which relevant TFs have binding sites in these regions in selected GTRD categories. The color of the DMR indicates whether the DMR is PSC-like and whether it is demethylated during neural cell differentiation. The number in front of the DMR indicates the k-means cluster which DMR belongs to (see Fig 1F). Gray DMRs were not included in TF binding and DMR cluster analysis as they were not AT/RT-specific.

Figure S11.. Oncoprint showing the DM-DE genes, which are differentially expressed in the tumor comparison, and the direction of differential gene expression.

For each comparison, different possible DMR locations (Gene neighbourhood within TAD, Promoter [2 kb], Enhancer [Fantom 5], Enhancer [GeneHancer]) are visualized separately.

Figure S12.. Exression and methylation heatmaps for DM-DE genes.

(A) Expression and methylation heatmaps for DM-DE genes using public microarray data. Expression on the left and methylation on the right. (B) Expression and methylation heatmaps for DM-DE genes using sequencing data (RNA-seq, RRBS). Expression on the left and methylation on the right.

Figure S13.. Expression and methylation heatmaps for tumor-specific DM-DE genes.

(A) Expression and methylation heatmaps for tumor-specific DM-DE genes using microarray data. Expression on the left and methylation on the right. (B) Same as in A but with sequencing data (RNA-seq and RRBS). Expression on the left and methylation on the right.

Figure S14.. Expression boxplots for selected DM-DE and NEUROG/NEUROD target genes.

(A) Expression of genes presented in Fig 4F from microarray data (GEO accession GSE42658). (B) Expression of AT/RT-unique DMRs with EZH2 binding site target genes from RNA-seq data. (C) NEUROG/NEUROD target genes from microarray data (GEO accession GSE42658). (D) NEUROG/NEUROD target genes from RNA-seq data.

Figure S15.. Methylation–expression correlation plots for selected genes.