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

DNA is hypermethylated in AT/RTs when compared to MBs or PLEXs

To study oncogenic epigenetic regulation in AT/RTs, we collected genome-wide DNA methylation Illumina microarray data (i450K) from 497 tumors and unmatched microarray expression data from 110 tumors; 89 normal brain DNA methylation samples were used as controls in i450k-based DNA methylation analysis. In addition, we generated matched data from 10 tumors with reduced representation bisulfite sequencing (RRBS) and RNA sequencing (RNA-seq) (Fig 1A, Table S1). We thus conducted two distinct analyses, one based on array data and the other on sequencing data, and used them to select the method-independent results in various subsequent analyses. tSNE analysis identified and confirmed the tumor types and subtypes of the cohorts (Fig 1B). MB samples were separated into WNT, SHH child/adult, SHH infant, group 3, and group 4 entities (Capper et al, 2018), largely following the latest 2021 WHO classification (Louis et al, 2021), whereas AT/RT and PLEX subtypes were closely positioned, as expected. Our in-house RRBS samples were positioned near the microarray sample clusters representing their clinical diagnosis.

Differential DNA methylation analyses between AT/RTs, MBs, and PLEXs resulted in 4,931 and 64,983 differentially methylated regions (DMRs) in the i450k (methylation difference: beta fold change >= 0.20; see details from the Materials and Methods section) and RRBS (methylation difference: beta fold change >= 0.25) data, respectively (Fig 1C). Higher DMR counts in RRBS data analysis were predictable because i450k DMRs were, on average, longer and RRBS has a wider genomic coverage (Meissner et al, 2005; Pidsley et al, 2013). Tumor location information and normal samples were used to reduce tumor location–related differences in i450k DNA methylation data (see the Materials and Methods section), decreasing the DMR counts, as well (Fig 1C versus Fig S1A). There were a higher proportion of intergenic DMRs in RRBS than in i450k data (Fig S1B, Table S1). DMRs were distributed across all chromosomes (Fig S1C). Most (>70%) of the DMRs were more methylated in AT/RTs and less methylated in MBs than in the other tumor types (Fig 1D). Consistently, the highest DNA methylation levels were detected in AT/RTs and lowest in MBs (Figs 1E and S1D and E).

In i450k data with a higher sample count, DNA methylation in DMRs was consistently higher in the AT/RT subgroups (MYC, SHH, and TYR) than in the subgroups of other tumor types (Wilcoxon rank sum test, P < 0.001, effect size Cliff’s delta |δ| > 0.2 for all comparisons except 0.18 for AT/RT-SHH versus PLEX, PED B and 0.13 for AT/RT-MYC versus PLEX, PED B) (Fig 1E, Table S1). In PLEXs, the highest DNA methylation was detected in the most aggressive subgroup, PLEX, PED B (Fig 1E). To subcategorize the i450k DMRs, we ran k-means clustering, which resulted in 10 clusters (Figs 1F and S2), of which three were ATRT-specific (clusters 1, 2, and 5), one PLEX-specific (cluster 9), and two MB-specific (clusters 7 and 10). Clustering also revealed tumor subgroup–specific DNA methylation patterns, especially for MBs (clusters 4, 6, and 8) and the AD subgroup of PLEXs. However, ATRT subgroups behaved similarly across clusters (Fig 1F) and also predominantly across individual DNA methylation sites (Fig S2A–J).

DNA methylation can be focal or spread to larger genomic regions, which are typically bordered by topologically associating domain (TAD) boundaries (Szabo et al, 2019). These two methylation types have different regulatory outcomes (Buitrago et al, 2021). The DMRs present in both RRBS and i450k data were involved in large-scale DNA methylation differences within specific TADs, especially in chromosomes 1, 11, 13, 17, and 19 in the AT/RT-MB comparison (Fig 1G).

Distinct neural differentiation–related transcriptional regulators are enriched in AT/RT- and MB-hypermethylated sites

To study the function and possible upstream regulation of regions showing differential DNA methylation, we used the chromatin immunoprecipitation sequencing (ChIP-seq) data collected from the Gene Transcription Regulation Database (GTRD) (Yevshin et al, 2019) for identifying the transcription factors and other transcriptional regulators (referred to as TFs from now on) whose binding sites are over-represented in the analyzed DMRs (adjusted P < 0.05, one-sided Fisher’s exact test). To focus our analysis on cancer-specific epigenetic regulation, we included only DMRs, which were hypermethylated or hypomethylated in a cancer-specific manner; that is, they behaved similarly in both comparisons (being, e.g., hypermethylated in AT/RTs when separately compared to MBs and to PLEXs) (Fig 2A). Only TFs enriched with both i450k and RRBS-based analyses were analyzed further.

Of the 183 enriched TFs, 136 (74%) were unique to one tumor type, 30 (16%) associated with two tumor types, and 17 (9.3%) to all three (Fig 2B, Table S2). Altogether, 77 TFs were enriched in AT/RT-specific hypermethylated (AT/RT-hyper) DMRs, 114 TFs in MB-hyper DMRs, and 10 TFs in MB-hypomethylated (MB-hypo) DMRs (Fig 2B). A higher proportion of MB-hyper DMRs were located in CpG islands (Fig S3A, Table S1), which might partly explain a high number of TFs enriched in these regions. No AT/RT subgroup specificity was observed in the DNA methylation patterns of differentially methylated TF-binding sites (TFBSs) (selected examples in Fig S3B). DNA methylation has been shown to reduce the DNA binding of 33/49 (67%) TFs with DNA binding activity and curated information (Table S2) (Yin et al, 2017).

To summarize the regulatory roles of the enriched TFs, we used a literature search to functionally annotate the TFs into different themes (Table S2), which were chosen based on TF functions and gene families that recurrently popped up when manually going through the enriched TFs in the first phase of the analysis (Merino et al, 2015; Swartling et al, 2015; Northcott et al, 2017; Erkek et al, 2019). Notably, 15 neural differentiation and neural cell–related TFs were uniquely enriched in AT/RT-hyper DMRs (Fig 2B). These included NEUROD1, ASCL1, and MYCN. When reanalyzing the TF binding after dividing the DMRs into clusters with k-means clustering (Fig 1F), AT/RT-unique TFs BCHE, MYCN, NFIB, NFIC, and PITX3 were enriched in AT/RT-specific DMR clusters 1 and/or 5 but none in AT/RT-specific cluster 2 (Table S2). Four neural differentiation–related TFs were enriched in both AT/RT-hyper and MB-hypo DMRs (Fig 2B, Table S2). Of these, NEUROG2 and HAND2 were enriched in cluster 2, whereas PHOX2B and PAX7 were enriched in clusters 1 and/or 5 (Fig 1F, Table S2).

We also noticed that 16 neural differentiation–related TFs were uniquely enriched in MB-hyper sites (Fig 2B). However, five of these (REST, SIN3A, ZEB1, TBX1, and ZHX2) inhibit neural differentiation (Table S2). In addition, four of the remaining TFs were active in the later phases of neural differentiation, whereas 7 of 15 neural TFs enriched in ATRT-hyper sites were active in earlier stages (during differentiation from PSCs and from neural stem cells) (Table S2). Thus, neural cell differentiation–related TFs in AT/RTs and MBs represent different functions.

When considering both RNA-seq and array expression data, the expression of 164 enriched TFs (90%) was not significantly associated with any tumor type (Figs S4A and B and S5A and B, Table S3). ERF was the only differentially expressed (DE) TF in AT/RTs compared with other tumors (adjusted P < 0.01, absolute logFC > 1), whereas 13 TFs, including EBF3 and NEUROD1, showed MB-specific differential expression. Neural differentiation TFs enriched in AT/RT-hypermethylated regions were expressed in AT/RTs, and some TFs (e.g., NEUROG2 and HAND2) even showed the highest expression in AT/RTs (Fig S5A).

GTRD data originate from a variety of cell types, and TF binding is also cell type–dependent (Yevshin et al, 2019). To take into account the stage of neural development, brain tumor–derived data, and cancer association, we categorized GTRD data based on the sample type in which the binding data of the TFs have been obtained. This resulted in 11 sample groups, which represented our analyzed tumor types, other brain tumors, neural progenitors, and more differentiated neural sample types, as well as PSCs, other tumor types, and all the other sample types (Table S2) (Yevshin et al, 2019). We then reran the TFBS enrichment analysis for cancer-specific DMRs using these 11 binding site subsets (called as categorized GTRD analysis from now on) and summarized the obtained results as a dot plot (Fig 2C, full version in Fig S6A and B). This approach helped us to investigate whether DNA hypermethylation was mainly covering TFBSs in neural cells, so sites that are bound by the TF in (more) differentiated cells of the neural lineage, whether TF binding has been measured from the analyzed tumor types, and whether TFs are shown to bind to these sites already in PSCs that are expected to represent the earliest phases of (neural) cell differentiation. Altogether, 16 of the analyzed TFs were both neural differentiation– and neural cell–related and measured in brain tumors or neural cells (Fig 2C). Of these 16 TFs, 11 (69%), including NEUROD1, were enriched in sites specifically hypermethylated in AT/RTs (AT/RT-hyper sites), MYCN was enriched in AT/RT-hypo sites, and ZEB1 was enriched in MB-hyper regions (Fig 2C, Table S2). On the contrary, 12 of 16 (75%) neural TFs enriched only in MB-hyper sites were measured or significant only in pluripotent cells or non-neural cell types (Fig S6A and B, Table S2). These included REST and SIN3A, which inhibit neural differentiation and whose binding sites in PSCs, but not those in “other brain and neural samples” or “other brain tumors,” were enriched in MB-hyper regions. Finally, the binding sites of 9 of 85 (11%) TFs uniquely enriched in MB-hyper regions were related to the PSC state and early differentiation (Fig 2B), and they significantly overlapped MB-hyper DMRs when TF binding was measured in PSCs (Fig 2C, Table S2). This enrichment of pluripotency factors in MB-hyper regions is implying that pluripotency programs are not active or are disrupted in MB cells, making it more likely that programs related to more mature cell states have taken over. Together, these findings suggest that the binding sites of neural TFs are hypermethylated in AT/RTs, which might interfere with their DNA binding, whereas TFBSs in PSCs are hypermethylated in MBs, thus representing a more advanced neural differentiation state.

DNA binding sites of SWI/SNF, PRC2, and EP300 harbor unique DNA methylation patterns in AT/RTs

As the SWI/SNF chromatin remodeling complex is incomplete in AT/RTs because of the full inactivation of SMARCB1 or SMARCA4, we separately analyzed its subunits in the categorized GTRD data. Interestingly, the binding sites of SMARCA4- and the PBAF (polybromo-associated BAF, mammalian SWI/SNF complex)-specific TF ARID2 in rhabdoid tumor samples were enriched in AT/RT-hypo and MB-hyper regions (Fig 2C) (Mashtalir et al, 2018), whereas SMARCA4 binding sites in other brain and neural cells were enriched in AT/RT-hyper and MB-hypo regions (Fig 2C, Table S2). Consistently, the acetyltransferase EP300 binding sites in normal neural precursor and differentiated neural cells were hypermethylated in AT/RTs, and EP300 binding sites in pluripotent cells were hypermethylated in MBs (Fig 2C). These results suggest that the neural SWI/SNF-EP300 target sites have been methylated in AT/RTs and unmethylated in MBs. Furthermore, the SWI/SNF complex remains able to promote active chromatin in the regions to which it binds even in the absence of SMARCB1.

PRC2 and EZH2 are antagonizing SWI/SNF and have been reported to populate neural SWI/SNF binding sites in AT/RTs (Erkek et al, 2019). In our analysis, EZH2 was enriched in AT/RT-hyper, MB-hyper, or MB-hypo regions in the uncategorized GTRD analysis, suggesting different modes of regulation (Fig 2C, Table S2). When we analyzed TF binding in AT/RT-hyper regions after dividing them into DMR clusters (Fig 1F), EZH2 was enriched in cluster 1, which shows the highest DNA methylation in AT/RTs and the lowest in MBs (Table S2). When examining the TF enrichment results in the categorized GTRD data, all the EZH2 binding sites, but those originating from other brain tumors, were enriched in AT/RT-hyper regions, including EZH2 binding sites in neural progenitors and more differentiated neural samples (Fig 2C).

PRC2 subunits and neural differentiation–related TFs rarely co-localize

We performed TF binding co-localization analysis for ATRT-hyper, MB-hyper, MB-hypo, and PLEX-hypo DMRs (Figs 2D, S7A–C, and S8). Most TFs related to neural differentiation share their binding sites and were clustered, but significant co-localization (adjusted P < 0.01, one-sided Fisher’s exact test) was rarely observed between PRC2 members and TFs related to neural differentiation (Figs 2D, S7A–C, and S8). In AT/RT-hyper regions, only the neural TFs BCHE and TFAP2A co-localized with PRC2 subunits. TFAP2A is known to regulate neural crest induction and cranial placode specification (Wang et al, 2011). BCHE is expressed early in the nervous system development and plays a role in neural stem cell development (Tiethof et al, 2018). The results suggest that PRC2-related DNA methylation largely affects parts of the genome other than the sites that recruit neural TFs during differentiation, BCHE and TFAP2A representing possible TFs that act at their intersection.

CUT&RUN reveals reduced binding of NEUROD1 in AT/RT DNA–hypermethylated sites

To experimentally validate our hypothesis of inadequate DNA binding of NEUROD1 in AT/RTs, we performed cleavage under targets and release using nuclease (CUT&RUN), sequencing experiment targeting NEUROD1 in one MB and three AT/RT cell lines (see the Materials and Methods section). For the analysis, we selected only the CUT&RUN peaks that overlapped with any of the previously reported NEUROD1 binding sites in the GTRD, or contained NEUROD1 binding motif in the sequence (Grant et al, 2011; Yevshin et al, 2019), or fulfilled both of these criteria (Fig S9A). We detected less binding of NEUROD1 in all of the three AT/RT samples (28, 1, and 12 peaks) compared with the MB sample (275 peaks) (Table S2). None of the peaks detected in AT/RT samples were both reported in the GTRD and to having NEUROD1 binding motif, whereas 215 (78%) peaks detected in the MB sample fulfilled both criteria (Fig S9A). Nearly all (273, 99.3%) of NEUROD1 binding sites detected in MBs were detected only in that sample in our CUT&RUN analysis (Fig S9B). NEUROD1 binding sites detected in AT/RT samples harbored no overlap with any of the DMRs or regions with cancer-specific DNA methylation (Figs 2E and F and S9C and D, Table S2). However, in the MB sample, NEUROD1 was binding to DMRs that were differentially methylated in the AT/RT versus MB comparison or were AT/RT-specific or MB-specific. The vast majority of sites in the AT/RT versus MB comparison (46/48 and 11/11 in RRBS and i450k data, respectively) harbored DNA hypermethylation in AT/RTs (Figs 2E and F and S9C and D, Table S2). This supports our previous claims, as we see NEUROD1 binding at sites hypomethylated in MBs. Similarly, the NEUROD1 binding sites in the MB sample that were overlapping with AT/RT-specific DMRs were predominantly hypermethylated in AT/RTs when compared to MBs and PLEXs (10/11 and 9/9 in RRBS and i450k data, respectively) and those overlapping with MB-specific DMRs were predominantly hypomethylated in MBs when compared to AT/RTs and PLEXs (27/28 and 7/7 in RRBS and i450k data, respectively), further supporting our hypothesis (Fig 2E and F, Table S2). Furthermore, we generated RRBS data from the cell lines used in the CUT&RUN analysis, which also confirmed that most of the NEUROD1 binding sites, which were detected in the MB sample, show high DNA methylation in all the AT/RT samples (beta value >0.5 in all the AT/RT samples in 258 of 275 [94%] sites) (Fig S9E). Although NEUROD1 RNA expression is on average higher in MBs than in ATRT (Fig S5), its protein levels were similar in all the studied cell lines (Fig S9F), so the differential expression of NEUROD1 cannot explain the observed differences in DNA binding.

Together, these findings support the notion that DNA methylation in AT/RTs is disrupting the binding of DNA methylation–sensitive TFs to these sites, thus altering or preventing the activation of neural differentiation–related programs.

AT/RTs harbor both unique and PSC-like DNA methylation patterns

Next, we analyzed DMRs in the context of normal neural cell differentiation using i450k DNA methylation data from PSC and primary FB samples (Nazor et al, 2012; Colunga et al, 2019) (Fig 3A). This allowed us to define whether tumor-specific DNA methylation levels actually reflect either low (PSC-like) or higher (FB-like) cell differentiation state or are unique to each tumor type, and whether they are altered during the differentiation from PSCs to FB. When visualizing the DNA methylation levels in DMRs (Fig 1E), the AT/RT was most similar to PSCs (Fig 3B). We observed MB-unique, AT/RT-unique, PSC-like, and FB-like DMRs, part of which change their DNA methylation levels during normal neural differentiation (Fig 3C, left side, Fig S10A, Table S3). Of the 898 DMRs hypermethylated in AT/RTs, 322 (39%) were similarly methylated in PSCs as in AT/RTs and demethylated during differentiation from PSCs to FB (group 3, Fig 3C). Only 31 (3.5%) of AT/RT-hypermethylated DMRs were FB-like (group 5), highlighting the higher similarity of AT/RTs to PSCs than FB. However, there were also 300 (33%) regions with AT/RT-unique DNA hypermethylation. Of these, 124 (14%) were not related to differentiation (group 2) and 176 (20%) demethylated during neural differentiation (group 1). To summarize, 75% of AT/RT DMRs were demethylated upon neural differentiation, and they include both PSC-like (39%) and AT/RT-unique (61%) regions (Fig 3C).

In contrast, 94% MB-hyper DMRs and 73% MB-hypo DMRs were unique to MBs (Fig 3C, Table S3). Only 17 (15%) unique MB-hyper and as many as 445 (48%) unique MB-hypo DMRs were hypomethylated during differentiation (including 321 [34%] DMRs, which were also PSC-like in the AT/RT comparison). Furthermore, 221 (25%) MB-hypomethylated DMRs were FB-like and hypomethylated during differentiation (group 13). Thus, a high proportion of even tumor type–unique DMRs change their DNA methylation during differentiation and are differentiation-related in a manner that underlines the similarity of MBs to FB and AT/RTs to PSCs.

To determine whether the observed DNA methylation differences were more linked to focal or large-scale regulation of DNA methylation, we calculated the proportion of DMRs involved in large-scale methylation (i.e., either hyper- or hypomethylated DMRs overpopulating certain TAD regions, Fig 1G) in each DMR category (Fig 3C, right side, Fig S10A). Notably, only 9 (8%) of MB-hyper DMRs represented large-scale DNA methylation differences, whereas 430 (50%) and 503 (55%) of AT/RT-hyper and MB-hypo DMRs, respectively, were involved in large-scale methylation (Fig 3C). Thus, MB-hyper DMRs appear to target specific cytosines or focal areas, such as specific promoters and enhancers (consistently with their enrichment in CpG islands), whereas both focal regulation and large-scale regulation of DNA methylation are observed for AT/RT-hyper or MB-hypo regions.

AT/RT-unique DNA methylation is associated with PRC2

We analyzed TF binding separately in the differentiation-related DMR categories presented in Fig 3C to define TFs relevant to different regulatory programs (Figs 3D and S10B). In general, the TFs linked to SWI/SNF and H3K27ac were associated with AT/RT-unique sites (Fig 3D). Furthermore, the binding sites of EZH2 and other PRC2 members (they contribute to histone H3 lysine 27 [H3K27] methylation, so we annotated them as H3K27me TFs) were enriched only in AT/RT-unique DMRs (Fig 3D). Many of the AT/RT-unique DMRs, which are demethylated during differentiation from PSCs to FB and carry EZH2 binding site(s), were linked to genes, such as IGSF9B, DOK7, and CAMK2B, which are involved in dendrite growth and differentiation from neural progenitor cells to mature state (Table S3) (Okada et al, 2006; Woo et al, 2013; Wu et al, 2017). These are the key processes the neuronal BAF complex (nBAF, mammalian neuronal SWI/SNF complex) is activating during neural development (Kadoch & Crabtree, 2015). However, some of the genes harboring AT/RT-unique differentiation–demethylated DMRs with EZH2 binding, such as TNK1, are found to be relevant in cancer development (Hong et al, 2020) and have not been associated with nBAF or late neural development.

The SWI/SNF subunit SMARCA4, whose binding in other brain and neural cells was associated with AT/RT-hyper DMRs (Fig 2C), was enriched in non–differentiation-related AT/RT-unique DMRs (Fig 3D). SMARCA4 enrichment remained significant (P = 0.0001, Fisher’s exact test), when only the binding sites of SMARCA4 in “other brain and neural” samples were included in the analysis. Interestingly, when using these binding sites, SMARCA4 was also significantly enriched in regions, which were similarly hypermethylated in AT/RTs and PSCs and hypomethylated during differentiation (group 3) (P = 1 × 10⁻⁷, Fisher’s exact test). Some members of the SWI/SNF complex, such as SMARCA4 and CHD7, were also enriched in the MB-unique DMRs that are hypermethylated in MBs and non–differentiation-related (group 10). In the case of SMARCA4, these sites might represent SMARCA4 binding sites in rhabdoid tumors (significantly enriched in MB-hyper DMRs, Fig 2C), which have been DNA-methylated in an MB-specific manner, and they can potentially be also linked to the altered functions of the SWI/SNF complex in MBs in the context of BRG1 mutation (Yi & Wu, 2018). However, fewer SWI/SNF TFs were enriched in MB-specific than AT/RT-specific hypermethylated DMRs.

Contrary to PRC2 and SWI/SNF subunits, more neural TFs (n = 12) were enriched in PSC-like AT/RT-hyper DMRs (Fig 3D). Notably, neural TFs BCHE and TFAP2A, the only neural TFs co-localizing with PRC2 subunits in AT/RT-hypermethylated sites (Fig 2D), were enriched only in AT/RT-unique DMRs (groups 1 and 2) (Fig 3D), linking them closely to PRC2-related DNA methylation.

In conclusion, our results suggest that AT/RT-specific DNA methylation is covering the sites that are bound by SWI/SNF and PRC2 complexes in neural progenitors and more differentiated neural cells and that the DNA hypermethylation in most of these sites is AT/RT-unique, not PSC-like. This aberrant DNA methylation has the potential to drive oncogenic epigenetic regulation, which is unique to AT/RTs. In addition, we detected AT/RT-hyper regions that were similarly methylated in PSCs, harbor binding sites for several neural regulators, and were originally methylated via a PRC2-independent mechanism.

DNA methylation contributes to neural cell differentiation–related gene expression

Next, we associated DNA methylation with gene expression based on the results obtained in both microarray- and sequencing-based analysis. A total of 1,305 DE genes were detected (adjusted P < 0.01, logFC >= 1 in both RNA-seq and microarray data) between tumor types (Fig 4A). In the following steps, these DE genes were separately integrated with RRBS and i450k DMRs. Microarray expression datasets represented all the MB subgroups with similar subgroup proportions as in the i450k data (Tables S1 and S4) (Rathi et al, 2020). Altogether, 134 (10%) DE genes harbor an opposite change in DNA methylation (called DM-DE genes) in their genomic neighborhood in at least one tumor comparison (Fig 4B). Furthermore, 42 and 67 DM-DE genes carry DNA methylation differences in the gene promoter and the gene-annotated enhancer region, respectively (Fig 4C and D). Most of the DM-DE genes (142 of 160 genes, 89%) were detected in only one tumor comparison (Fig S11). In addition, a parallel increase/decrease in both expression and DNA methylation was detected in 130 genes, of which 70 were in AT/RT comparisons (62 in AT/RT-MB and 8 in the AT/RT-PLEX comparison, representing different genes) (Table S4).

DM-DE genes were predominantly hypermethylated in AT/RTs and underexpressed when compared to PLEXs, MBs, or both (Figs 4E, S11, and 12A and B). CXXC5 and TCEA3 genes were uniquely hypermethylated and down-regulated in AT/RTs, and seven genes (including neural genes NEUROG1, NBEA, and BARHL1) were uniquely hypomethylated and overexpressed in MBs (Fig 3E and F and S13A and B). The number of AT/RT-specific genes was partly reduced because decreased DNA methylation in both MBs and PLEXs led to increased expression in only one of them.

We identified NEUROG1 and NEUROD2 as DM-DE genes hypermethylated in AT/RTs (Fig 4E), and NEUROG2 and NEUROD1 as TFs whose binding sites were enriched specifically in AT/RT-hypermethylated regions (Figs 2C and S6). Because these TFs are key regulators of neural differentiation (Massari & Murre, 2000; Dixit et al, 2014), we visualized their expression together with the expression of their target genes TBR1, RCOR2, BLHE22, NHLH1, NHLH2, HES6, DAB1, and CDK5R1 (Fig S14A–D). Except for NEUROG2, which had the highest expression in AT/RTs, the highest expression of all these genes was observed in MBs (Fig S14A–D). This result suggests that their induction, which should occur during neural differentiation, has been successful mainly in MBs.

AT/RT-specific DNA methylation in genes suppressed in AT/RTs

To better understand the processes involved in AT/RT-related DNA hypermethylation, we more closely inspected the genomic neighborhoods of CXXC5, TCEA3, NEUROG1, NEUROD2, and EBF3 genes (Fig 4G). Typically, one overlapping DMR was detected in the gene loci with both used measurement techniques (RRBS and i450k) with the exception of the CXXC5 locus, which is located in a large-scale AT/RT-hyper region (in Chr5q31.2) and includes several DMRs (Figs 1G and 4G, S14A–D, and S15, Table S4). An AT/RT-specific DM-DE gene CXXC5 (Fig 4G) also harbors binding sites of several neural TFs, SMARCA4 binding sites in rhabdoid tumor cells, and both PSC-like and AT/RT-unique DMRs. In addition, EP300 binding sites measured from neural cells were also found around CXXC5, representing possible neural SWI/SNF-EP300 target sites described earlier (Fig 4G). TCEA3, also an AT/RT-specific DM-DE gene, carries AT/RT-unique, non–differentiation-related DMRs and neural EZH2 binding sites at the gene promoter.

Interestingly, both NEUROG1 and NEUROD2 were hypermethylated and down-regulated in AT/RTs when compared to MBs in our DM-DE gene analysis and they harbor DMRs, which were hypermethylated in AT/RTs compared with PSCs and bound by BRD9 in rhabdoid tumors (Fig 4G). In addition, NEUROG1 included EZH2 binding sites measured from several different cell types. NEUROG1 and NEUROD2 were also lowly expressed in all PLEX subgroups (Figs 4E and F, S13A and B, and S14A–C) and are likely to be silenced in this differentiation lineage either via alternative ways to suppress gene activity, such as H3K27me3, or by the absence of transcriptional activators in the regulatory regions. In MBs, these two TFs behaved in a complementary manner, as NEUROG1 was induced in all the other MB subgroups except in the SHH subtype, whereas NEUROD2 induction was more clearly detected in SHH and G4 subtypes (Figs 4E and S13A and B).

Finally, we wanted to investigate whether NEUROD1 binding was detected in the DMRs related to the DM-DE genes. Based on our CUT&RUN data, NEUROD1 binds in the MB cell line to the site that was connected to higher EBF3 expression in MBs (Fig 4G, Table S4). This site was DNA-hypomethylated in MBs and is now shown to be bound by NEUROD1 in MBs, but not in AT/RTs. Interestingly, both SMARCA4 and EP300 bind to the same region in “other brain and neural samples” in the GTRD, thus suggesting that it is one of the genomic sites that these regulators are targeting during neural differentiation. Furthermore, this site has been annotated as a distal enhancer by the Encyclopedia of DNA Elements (ENCODE) project. EBF3 has a role in neurogenesis, and it has also been shown to have tumor suppressor activity (Tao et al, 2015; Tiethof et al, 2018). Interestingly, EBF3 has been shown to be regulated by NEUROD1 (Seo et al, 2007), and our analysis positioned NEUROD1 more specifically to its regulatory region.

Overall, our results suggest that AT/RT-related DNA methylation is involved in suppressing these genes, which can contribute to impaired further neural development characterized by AT/RTs.

Sequencing cohort

Tumor samples analyzed with next-generation sequencing were collected for this study from 10 patients operated at Tampere University Hospital. The study was performed in line with the principles of the Declaration of Helsinki. Approval for this study was granted by the Regional Ethics Committee of Tampere University Hospital (decisions: R13050, date 9.4.2013, R14024, date 28.3.2014, and R07042, date 20.9.2017) and Valvira (decision V/78697/2017). Written consent was obtained from patients except for 1 MB and 2 PLEX samples, which were collected as part of the different cohorts. An experienced neuropathologist evaluated the tumor samples and determined the histopathological type and grade according to the criteria presented by the WHO 2016 (Louis et al, 2016).

Cell culture

AT/RT cell lines CRL-3020, CRL-3036, and CRL-3038 and a medulloblastoma cell line CRL-3021 were cultured in DMEM/Hepes with 2% B27, 1% penicillin–streptomycin, 20 ng/ml FGF, and 20 ng/ml EGF. All the cells were cultured in an incubator at 37°C with 5% CO2.

RRBS

DNA was isolated from the frozen samples and cell lines using a QIAamp DNA Mini kit (QIAGEN) with RNAse treatment and from the FFPE sample with a turXTRAC FFPE DNA kit (Covaris). Library construction and sequencing were performed at the Beijing Genomics Institute (BGI), Hong Kong. Samples were sequenced using an Illumina HiSeq 2000 technology. Paired-end sequencing of 50 bp was used. The data goal was 5 Gb. Cell lines were sequenced at Novogene, with Illumina NovaSeq paired-end 150 bp. The data goal was 15 Gb.

RNA sequencing

RNA was isolated using a mirVana isolation kit (Invitrogen). Library construction and sequencing were performed at Novogene, Hong Kong. Samples were sequenced using an Illumina HiSeq 2000 technology. Paired-end sequencing of 150 bp was used. The data goal was 20 million raw reads per sample.

Western blotting

Samples from cell lines were scraped on ice with RIPA cell lysis buffer (Thermo Fisher Scientific) supplemented with protease inhibitors (Roche). Proteins were separated by a Ready Gel precast gel (Bio-Rad) and transferred onto a nitrocellulose membrane. Blocking was done with 5% milk, 1 h at RT. Primary antibody incubations were performed overnight at 4° for both polyclonal NEUROD1 (1:500, PA5-47381; Invitrogen) and monoclonal β-tubulin (1:20,000, SAP.4G5). Secondary antibodies (polyclonal rabbit anti-mouse immunoglobulins/HRP [1:1,000, P0260; Dako] and polyclonal rabbit anti-goat immunoglobulins/HRP [1:2,000, P0160; Dako]) were incubated for 1 h at RT. All antibody dilutions were made to 5% milk. Proteins were detected by enhanced chemiluminescence (Bio-Rad).

CUT&RUN sequencing

Each cell line was harvested by removing the media with cells from flasks and centrifuging the media at 175g for 7 min. The supernatant was removed, and cells were washed once with PBS followed by centrifugation at 175g for 7 min. Cells were counted, and 100,000 cells were collected for each reaction. CUT&RUN was performed using the Cell Signaling Technology CUT&RUN Assay kit (#86652). The experiment was done as described in the kit’s protocol. NEUROD1 antibody #PA5-47381 from Invitrogen was used. DNA concentrations were measured using Qubit 2.0. Library construction was done using DNA Library Prep Kit for Illumina Systems (ChIP-seq, CUT&RUN; #56795). Sequencing was performed at Novogene, UK. Samples were sequenced using an Illumina NovaSeq X technology. Paired-end sequencing of 150 bp was used. The data goal was 15 million raw reads per sample.

TrimGalore! version 0.6.7 was used to trim adapter reads. Reads with a length of less than 25 bp were filtered out. Reads were aligned to hg38 using bowtie2 version 2.4.5 and the –dovetail flag. Duplicated reads were removed using Picard MarkDuplicates version 2.27.1. Regions overlapping with the regions described by Nordin et al (2023) were removed using bedtools intersect version 2.29.1. Afterward, the reads with a length shorter or equal to 120 bp were extracted using samtools view (version 1.8). MACS3 version 3.0.0.bl was used to call the peaks for the fragments with the indicated length. Peaks located outside of the canonical chromosomes and on chrX or chrY were filtered out. Peaks that overlapped with peaks in the negative control samples (input DNA) were removed using bedtools intersect version 2.29.1. Replicates were unified by keeping peak regions called in both replicates. For all remaining peaks, the fasta sequence was extracted using bedtools getfasta. The NEUROD1 motif was downloaded from JASPAR 2020 (Fornes et al, 2020), and all peaks were scanned for the motif using FIMO (part of MEME suite v. 5.5.5) (Grant et al, 2011), keeping all peaks with a detected motif (P < 0.001). The peaks obtained with this method were combined with peaks overlapping with GTRD v.19.04 NEUROD1 binding sites and overlapped with obtained DMRs in R, using the GenomicRanges package (version 1.34.0).

Preprocessing and analysis of DNA methylation datasets

TF binding analyses

Differential expression analysis for expression datasets and integration with DNA methylation data

Linking DMRs to TADs for large-scale DNA methylation analysis

DMRs were overlapped with Cortex DLPFC and H1-NPC TADs using Bedtools GroupBy v2.25.0 and in-house scripts to calculate hyper- and hypomethylated records in each TAD (Quinlan & Hall, 2010; Wang et al, 2018). The input material was either a broadened (±500 bp) probe extracted from the i450k DMRs or the methylKit DMRs of the RRBS analysis. The resulting TADs for each comparison were filtered such that at least five similarly directed probes/tiles, each with a DNA methylation difference of at least 25% with both the i450k and RRBS approach, should hit the TAD. Simultaneously, less than 1/10 of the hits can undergo DNA methylation changes in the opposite direction. Moreover, the distance between the first and last probe or tile in the TAD had to be over 50 kb. The filtered TADs from different cell lines were pooled and reduced to produce single TAD coordinates for each comparison and then visualized with karyoploteR (Gel & Serra, 2017).

Comparing brain tumors with ESCs and primary fetal brain

Two PSC samples, nine primary FB samples (GSE116754) (Colunga et al, 2019), and 25 PSC samples (GSE31848) (Nazor et al, 2012) were downloaded from GEO and preprocessed using the Capper et al data (Capper et al, 2018). Beta-value density plots and tSNE clustering for the 10,000 most varying probes were used to check that samples were separated based on the tumor and normal sample type. The DMRs were called with DMRcate by comparing AT/RT, MB, and PLEX samples against PSC and FB cells. In addition, PSC samples were compared with FB.

From these seven comparisons, hypermethylated and hypomethylated DMRs were extracted (q < 0.05, abs. average beta fold change >= 0.25, and length between 100 and 5,000 bp). The region lists were overlapped with the same cancer-specific regions extracted for TFBS analysis (min 1 bp overlap). As a result, the experiment was summarized in a “developmental DMR matrix” of 2,301 non-overlapping regions showing which of the comparisons or cancers each region overlapped. The resulting region sets were categorized based on the differences detected in tumor–tumor, tumor–normal, and PSC-FB comparisons. TFBS enrichment analysis was performed for these categories by extracting and extending the probes hitting the regions and testing them against the GTRD using one-sided FET, as described.