{'fraction_non_alphanumeric': 0.03238546603475513, 'fraction_numerical': 0.027606635071090048, 'mean_word_length': 5.097038285576692, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Genome-wide methylation profiling identifies novel methylated genes in neuroblastoma tumors

ABSTRACT Neuroblastoma is a very heterogeneous tumor of childhood. The clinical spectra range from very aggressive metastatic disease to spontaneous regression, even without therapy. Aberrant DNA methylation pattern is a common feature of most cancers. For neuroblastoma, it has been demonstrated both for single genes as well as genome-wide, where a so-called methylator phenotype has been described. Here, we present a study using Illumina 450K methylation arrays on 60 neuroblastoma tumors. We show that aggressive tumors, characterized by International Neuroblastoma Risk Group (INRG) as stage M, are hypermethylated compared to low-grade tumors. On the contrary, INRG stage L tumors display more non-CpG methylation. The genes with the highest number of hypermethylated CpG sites in INRG M tumors are TERT, PCDHGA4, DLX5, and DLX6-AS1. Gene ontology analysis showed a representation of neuronal tumor relevant gene functions among the differentially methylated genes. For validation, we used a set of independent tumors previously analyzed with the Illumina 27K methylation arrays, which confirmed the differentially methylated sites. Top candidate genes with aberrant methylation were analyzed for altered gene expression through the R2 platform (http://r2.amc.nl), and for correlations between methylation and gene expression in a public dataset. Altered expression in nonsurvivors was found for the genes B3GALT4 and KIAA1949, CLIC5, DLX6-AS, TERT, and PIRT, and strongest correlations were found for TRIM36, KIAA0513, and PIRT. Our data indicate that methylation profiling can be used for patient stratification and informs on epigenetically deregulated genes with the potential of increasing our knowledge about the underlying mechanisms of tumor development.


Introduction
Neuroblastoma (NB) is a pediatric extracranial solid tumor of the peripheral sympathetic nervous system showing vast clinical heterogeneity. The tumors are characterized by patterns of genetic aberrations with unbalanced translocations and focal changes in ploidy. The proto-oncogene MYCN is amplified in 20% of the primary tumors, 1 and is associated with a more aggressive disease. 2 Rare but recurrent somatic mutations have been identified in the anaplastic lymphoma receptor tyrosine kinase (ALK) gene in both sporadic and familial cases. [3][4][5][6][7] The International Neuroblastoma Risk Group (INRG) classification system classifies tumors into localized tumors (stage L1 and L2), metastatic tumors (M), and metastatic tumors with particular features in children younger than 18 months (MS). 8,9 Aberrant DNA methylation has emerged as an important feature of both development and progression in many cancers. Although the full scope of regulatory mechanisms exerted from DNA methylations remain to be elucidated, certain modes of the regulatory effect on expression patterns have been discovered. For instance the hypermethylation of CpG islands in promoter regions are generally associated with reduced transcriptional activity. Thus, this adds to yet another mechanism that can silence a tumor suppressor gene, apart from deletional events and mutational gene silencing, and thus contribute to induction of malignancy. However, not only gene promoter regions harbor deregulated methylation sites with transcriptional regulatory impact. In other malignancies, the combined analysis of global DNA methylation patterns and global transcription patterns have revealed that hypomethylation of intragenic methylation sites are associated with increased gene expression, but for some genes also the opposite is true. In neuroblastoma the transcriptional regulation of DNA methylation has also been shown to be valid for a smaller subset of genes. 10,11 Furthermore, genomic regions such as enhancers and silencers are known to have cis-regulatory importance but can be positioned distant from its target gene.
Recently, species of non-coding RNAs have emerged as important regulators of cancer, and the dysregulation of these RNAs have been reported to be associated with regional methylation alterations in neuroblastoma. 12 Most analyses of DNA methylation in NB tumors have been performed on single genes or with limited sample sizes. Several genes in various cellular pathways (apoptosis, cell cycle, differentiation, invasion, and metastasis) have been identified as aberrantly methylated. The study by Alaminos et al. was one of the first to demonstrate that the clustering of NB tumors based on the methylation profile of 10 genes could divide NB into clinical risk groups. 13 Since then, genome-wide analyses of DNA methylation have revealed a DNA methylator phenotype in NB with poor prognosis, characterized by the methylation of a set of multiple CpG islands. 14 A study of 7 genes in 2 NB cell lines (tumorigenic LA1-55n and non-tumorigenic LA1-5s) by Yang et al. reports that tumorigenic properties can be inhibited by reversing epigenetic changes with a demethylating agent, suggesting that epigenetic aberrations contribute to the NB phenotype. 15 Recent advances in high-throughput technology in the field of epigenetics now enable large-scale studies using bisulfite modified DNA hybridized to DNA methylation arrays. We have previously used the Illumina HumanMe-thylation27 DNA analysis BeadChip to determine methylation profiles of NB primary tumors and identified aberrantly methylated genes that can be used as biomarkers. 16 Extending on that study we now used the Illumina 450K arrays, profiling more than 450,000 CpG sites, for a large set of NB tumors.

Results
Clustering the samples based on methylation data identifies a NB methylator phenotype Genome-wide methylation analysis using the Infinium HumanMethylation450 BeadChips was performed on a set of 60 primary NB tumors. The sample set contained 17 tumors with MYCN amplification, 9 samples with 11q-deletion and 12 with ALK mutations, and is further described in Supplementary Table 1 and in Fig. 1. After quality control analysis and filtering of methylation data, a set of 364,774 probes remained for further analyses. To get an overview of the methylation patterns in the data, we applied unsupervised clustering of the 1% most variable probes (nD3,648). The analysis separated samples into 2 clusters (Fig. 1A) that discriminated also by INRG classification (PD0.00017), deletion of chromosome 1p (PD0.015) and on MYCN-amplification (MNA) (PD0.034) (Fig. 1A). The two INRG MS cases in the cohort clustered together (Fig 1A). Furthermore, the probes that were used for clustering differed in mean methylation levels (Fig. 1B, PD2.2£10 ¡6 , t-test), and in 5-year overall survival (OS) (PD0.014; Fig 1A and 1C). The CpG site with largest difference in methylation b value between the 2 clusters was cg17878351, located in the gene CLIC5. As the samples separated most strongly based on INRG stage, we focused our further analyses on this grouping.

Sites with variable methylation are hypermethylated in INRG stage M tumors
Using the ChAMP package 17 we compared the groups of INRG metastatic (M, 26 samples) and INRG localized (L, 32 samples). In total, 37,891 probes gave an adjusted P-value below 0.05, and 4,557 of these probes showed a difference in b-value methylation between the INRG M and L groups that was larger than 0.2 (20%). We defined this set of probes as the Methylation Variable Positions (MVPs) ( Fig. 2A and Supplementary Table 2). More than 80% of these MVPs were hypermethylated in INRG M ( Fig. 2A; 3,738 hyper-and 819 hypo-methylated). The genomic distribution of the MVPs in relation to CpG islands (CGI) was compared to the reference distribution of probes in relation to CGIs for the data set of 364,774 probes for each locational feature analyzed. MVPs were overrepresented in regions of low CpGs ("ocean") and underrepresented in CGIs (P<0.001 for both, Pearson x 2 test). Also CGI shelves and shores showed a non-random distribution between the identified MVPs and the genomic CGI representation on the analyzed platform ( Fig. 2B; P<0.01 for shelves and shores, respectively). When the MVPs were divided into hyper-or hypo-MVPs, both the direction of distributional change and significance remained for probes located in oceans and CGI ( Fig. 2B upper panel). In hyper-MVPs, shores were overrepresented compared to the reference distribution (P<0.0001, Fig. 2B), whereas localization of MVPs in shelves were underrepresented compared to reference (P<0.001, Fig. 2B). On the contrast, hypo-MVPs were underrepresented at shores compared to the reference (P<0.001; Fig. 2B).

INRG differentially methylated sites are enriched in enhancer regions
When MVPs were mapped in relation to gene positions we found a significant overrepresentation of MVPs in intergenic regions (IGR; P<0.0001) and this significance remained also when subdividing into hyper-and hypo-methylated MVPs ( Fig. 2B lower panel). An increased representation of MVPs in gene body (P<0.01) for MVPs was not seen when subdividing into hyper-and hypo-methylated MVPs. An overall underrepresentation in MVPs located in the first exon (P<0.0001, Fig. 2B) was present in both hyper-and hypo-MVPs (P<0.0001 and P<0.0012, respectively; Fig. 2B). Furthermore, there was an underrepresentation of MVPs in the upstream region of transcription start site (TSS) that remained when analyzed separately for in hyper-and hypo-MVPs, both in the region 200 bases upstream of TSS and in the region 1,500 bases upstream of TSS (P<0.0001 for all, Fig. 2B lower panel). In line with this, there was an underrepresentation of MVPs in 5 0 untranslated region (UTR). However, the representation of MVPs in 3 0 UTR was similar as in the reference data, both when hyper-and hypo-MVPs were analyzed together and separately. Both hypermethylated and hypomethylated MVPs were enriched in enhancer sites (P<0.0001), suggesting a regulatory function of the differential methylation.

Location of MVPs
The distributions of the 3 most discriminatory hyper-MVPs and 3 hypo-MVPs for INRG are shown in Fig. 2C. The cg13441112 is located on chromosome region 1q32 within the cluster of complement factor H related genes. Its genomic feature is classified as intergenic, although its position can be intronic by alternative splicing of the nearby complement factor H-related 4 (CFHR4) gene. The INRG M hypermethylated cg07469926 and cg23195028 are located in gene bodies of the genes tripartite motif containing 36 (TRIM36; chr5q22.3) and protein phosphatase 2, regulatory subunit B, gamma (PPP2R2C; chr4p16.1), respectively. The three most selective hypomethylated probes were located in intergenic regions on chromosome region 1p36 (cg02971581 and cg09510180), and on chromosome 13 (cg25233139). Kaplan-Meier survival curves based on the methylation of these sites are shown in We next explored the predictive capacity of using tumor methylation signatures for improved patient stratification for the 23 patients who died from the disease and the 18 patients with >60 months of no evidence of disease. We performed multivariate analysis (logit regression) using INRG stage and MNA status as predictors of outcome, which gave a predictive accuracy of AUC D 0.86. When using the methylation cluster and the 6 top discriminating CpG sites mentioned above, the predictive accuracy increased to 0.95 (P-value for difference between the models P<0.05). Combining the methylation factors with INRG stage and MNA status did not further increase the accuracy of the model (AUC D 0.96).
The 4,557 identified MVPs mapped to 1,516 genes; out of these, 1,265 genes were hypermethylated and 290 genes hypomethylated (Supplementary Table 2) in the INRG M group. The genes with the highest number of MVPs (Table 1) were the hypermethylated genes protocadherin gamma subfamily A, 4 (PCDHGA4; chr5q31) and telomerase reverse transcriptase (TERT; chr5p15). PCDHGA4 is located in a cluster of multiple overlapping protocadherin genes and has previously been recognized as hypermethylated in NB. 14 The TERT gene had 24 MVPs, the majority of which were located in the gene body (Fig. 3A). The gene list also contains the 2 adjacent genes DLX5  (Table 1). Genes with both hyper-and hypo-MVPs were also identified (nD39; detailed in Supplementary Table 3). The gene PR domain containing 16 (PRDM16), located at chromosome 1p36, a region frequently deleted in neuroblastoma, was the only gene showing at least 2 of each hyper-and hypo-MVPs.

Gene ontology analysis identifies NB relevant terms
To evaluate the functions of the genes with associated MVPs we used gene ontology classification through the DAVID bioinformatics resource. 18 Genes containing either hyper-or hypomethylated MVPs were submitted to default settings in DAVID. This analysis showed a representation of neuronal tumor relevant gene functions, such as cell adhesion, neuronal development, neuronal differentiation, and neuron projection development (Fig. 3B).

Differentially methylated regions
To identify coherent chromosomal areas of aberrant methylation the dataset was analyzed for differentially methylated regions (DMRs). Using the default P-value of 0.01 in the ChAMP package we identified 459 DMRs (spanning 21-7,185 base pairs) distributed over 2,425 probes. One fourth (26%) of these probes was classified as MVPs in the above analysis. The DMRs represented 283 annotated genes out of which 204 genes were included among the identified MVP genes above. The DMRs with the lowest P-values are listed in Table 2 (B3GALT4,   KIAA1949, DLX6AS, HLA-DOA, intergenic region at 12q24.33, GABBR1). As we noted an uneven distribution among the chromosomes for the DMRs, we visualized the genomic locations of DMRs located to genes (Fig. 3C). Interestingly, multiple DMRs were located at chromosomal regions 1p36 and 17q, 2 regions that frequently show copy number alterations in neuroblastoma, and on chromosome 6 several genes with DMRs were located in the MHC region.

Analysis of other prognostic factors in NB and overall survival
As a next step, we looked for methylation discriminating between samples with or without MYCN-amplification, 5-year OS, and 11q deletion. We applied the same strategy for MVP calling (volcano plots are shown in Supplementary Fig. 2 Supplementary  Table 7. Only 8 probes were called MVPs in all 4 analyses and these were located in the genes HORMAD2 (2 probes), SUSD5, PPAPDC1A, NKPD1 and 2 intergenic probes located on chromosome 2 and 9 ( Supplementary Fig. 4, Supplementary Table 8).

Non-CpG methylation separates NB patients
The Illumina arrays also contain 3,091 non-CpG sites, out of which 1,678 remained after filtering (see Materials & Methods). Cluster analysis of the 10% most variable sites (nD168) separated the patients into separate clusters, of which also here INRG was the strongest clinical variable (Fig. 4A). For non-CpG sites, however, the highest methylation frequency was found in the INRG stage L (PD0.001; compared to INRG M, Fig. 4B).
Validation of the obtained methylation profiles using published data sets  Supplementary Fig 5A). The distribution of the 3 most top ranked hyper-and hypo-MVPs in this data set showed a similar trend as in our dataset (Supplementary Fig 5B).

Expression of genes with differential methylation in prognostic groups
To further look into the biological relevance of the genes with identified MVPs, we used the R2: Genomics Analysis and Visualization Platform (http://r2.amc.nl). We used the Kocak set of 649 NB tumors (GEO ID: GSE45547). 20 Many of the genes identified in the current study also showed significant differences in gene expression when it comes to overall survival of the patients (Supplementary Fig 6). For example, the genes B3GALT4 and KIAA1949, identified as hypermethylated in INRG stage M showed lower gene expression in patients who died from their disease (Bonferroni corrected P-values 5.4 Ã 10 ¡14 and 1.6 Ã 10 ¡17 , respectively). The gene DLX6AS1, which showed hypermethylation in the gene body (indicative of stable gene expression), showed high expression in patients who died from their disease (Bonferroni corrected P-values, 5.7 Ã 10 ¡13 ). The CLIC5 gene, with the top site that discriminated the 2 methylation clusters, showed lower expression in non-survivors (Bonferroni corrected P-values, 2.1 Ã 10 ¡5 ). We also verified that the same was true using other gene expression data sets (using the SEQC_NB data set composed of 498 samples, GEO ID: GSE62564). The PIRT gene, with the CpG sites with the

Correlation of DNA methylation and gene expression in neuroblastoma tumors
To further explore if the identified genes with differential methylation and differential expression were also correlated to each other in individual tumors, we used a recently published smaller set of NB tumors for which paired methylation and gene expression data was available for 19 samples (GEO accession numbers GSE54719 and GSE54720). 19 Correlation between DNA methylation and gene expression with nominal significance P<0.05 were found for 16 of 39 genes tested (Supplementary Table 9). The three strongest correlations were found for TRIM36 (rho=-0.83), KIAA0513 (rhoD0.80), and PIRT (rhoD-0.79); for TRIM36 and KIAA0513, the tested methylation site was located in the corresponding gene body; for PIRT, the methylation site was located within 200 bases upstream of the transcription start site (Supplementary Table 9 and Supplementary Fig. 5C).

Discussion
Current risk stratification for NB tumors is based on MYCNamplification status and presence of chromosomal segmental breaks. Attempts have been made to improve classification based on methylation of single-genes, or genome-wide profiles. 16,[21][22][23][24][25][26] The studies are often limited to small cohorts without validation sets. Here we present a larger study of 60 patients using Illumina methylation arrays, which profile 485,000 CpG sites; as validation, we use our previously published study using the 27,000 CpG site platform. 11 Cluster analysis identified 2 clusters with different methylation levels: the hypermethylated cluster had a mean methylation of 0.64 and the hypomethylated cluster had a mean methylation of 0.43. The patients in the hypermethylated cluster had a poorer survival than the other group of patients. The methylation profiles separated the NB tumors into clinically relevant groups, of which INRG stage was the strongest clinical factor. Gene ontology analysis of the identified genes identified relevant terms as cell adhesion, neuronal development, neuronal differentiation, and neuron projection development. One of the top differentially methylated genes in INRG stage M was the telomerase reverse transcriptase gene TERT. Telomerase expression plays a role in cellular senescence, as it is normally repressed in postnatal somatic cells resulting in progressive shortening of telomeres. Deregulation of telomerase expression in somatic cells may be involved in oncogenesis. Gene expression of TERT is significantly higher in NB patients who die from their disease (R2: Genomics Analysis and Visualization Platform). Methylation of TERT has recently been proposed as a biomarker for risk stratification of childhood brain tumors, as hypermethylation of a region upstream of the transcription start site of TERT was found to be associated with tumor progression and poor prognosis. 27 A high-risk neuroblastoma subgroup defined by genomic rearrangements causing increased expression of TERT was recently presented. 28,29 In our study, TERT had the highest number of hypermethylated sites in INRG M (Table 1), suggesting altered epigenetic regulation, which is consistent with a high-risk TERT subgroup, although the genomic rearrangements are not known for our samples. TRIM36 was hypermethylated in INRG stage M and has lower expression in stage M tumors. Overexpression of TRIM36 has been described to decelerate the cell cycle and limit cell growth 28 and TRIM-NHL proteins are important regulators of developmental transitions, for example, by promoting differentiation and repressing cell growth and proliferation in stem and progenitor cells. 29 Also, the PIRT gene shows the same pattern in NB. PIRT functions as a regulatory subunit of the TRPV1 protein with roles in cell migration, cytoskeleton reorganization, and in neuronal guidance. TRPV1 could be important during neuronal differentiation as its expression is increased during all-trans-retinoic acid differentiation in NB cell lines. 30 As validation, we used a set of NB tumors previously analyzed on the Illumina 27K platform 11 and also here separate clusters were formed where the hypermethylated group of patients was associated with INRG stage M and a poor survival.
Our results fit well into and extend on previous studies. A hypermethylator phenotype linked to poor outcome has been suggested previously, which our study supports. 14 For single genes, for example the KRT19 gene has previously been suggested as a marker of clinical risk factors in NB and was also in the current study identified as differentially methylated in INRG stage M compared to stage L. 11,21 A recent study of a smaller number of NB tumors also reported on differences in non-CpG methylation where the presence of non-CpG methylation was mostly associated with tumors characterized by favorable clinicobiological features, 19 which the current study confirms.
Our data show that methylation profiles can separate NB tumors into relevant clinical groups. It also highlights genes that, in addition to the differential DNA methylation, also show differential gene expression of the respective gene, suggesting that these genes might have a role in tumor development or progression, which needs further functional studies to address.

Patients and samples
A panel of 60 neuroblastoma tumors was used in this study (Supplementary Table 1). Ethical permission was granted by the local ethics committee.

DNA extraction and bisulfite modification
DNA extraction was carried out with the DNeasy Blood & Tissue Kit (Qiagen, Hilden, Germany), including a RNA removal step, according to the protocol provided by the supplier. The DNA was quantified with the Nanodrop and 1 ug was used for bisulfite modification using the EZ DNA methylation kit (D5001, Zymo Research, Orange, CA) according to the protocol provided with the modification step according to the recommendations for array processing of the samples. Control PCR reactions were carried out before array analysis to confirm successful modification of the DNA.

DNA methylation arrays
The bisulfite-modified DNA (500 ng) was applied to the Infinium HumanMethylation450 BeadChips (Illumina), which determine the methylation levels of 485,000 CpG sites. 31 After bisulfite treatment of the DNA samples, the cytosines in the CpG sites were genotyped as C/T polymorphisms according to the manufacturer's protocol. The fluorescence signals were measured from the BeadArrays using an Illumina BeadStation GX scanner. The fluorescence data were then analyzed using the GenomeStudio software (Illumina). The software assigns a score called a "b value" to each CpG site, which corresponds to the ratio between the fluorescence signal of the methylated allele (C) and the sum of the fluorescent signals of the methylated (C) and unmethylated (T) alleles. The data generated by the BeadStudio software was exported and further analyses were performed in the R programming environment.
The R-package ChAMP 17 was used for data preprocessing, normalization and comparison between groups. Singular value decomposition analysis was performed to identify confounding factors.
Positions on the X and Y chromosome were filtered away using ChAMP and we applied a filter for SNPs from the IMA package (https://www.rforge.net/IMA/); hereafter, a set of 364,774 probes remained for analyses. Probes with highest variation were identified using median absolute deviation. ChAMP was used to calculate probes that were differentially methylated between groups. We defined methylated variable positions (MVPs). We set the criteria for MVP as calling at significance of a Benjamini-Hochberg adjusted P<0.05 and a difference in b-value between groups larger than 0.2. Hierarchical clustering of samples was done with the heatmap function in the R NMF package using Euclidean distances. Samples classified as noevidence-of-disease (NED) with less than 60 months follow-up time (nD19) were omitted from the survival analysis.
Analysis of differentially methylated regions (DMRs) was done using ChAMP. The web-based program Phenogram (http://visualization.ritchielab.psu.edu/) was used to map and visualize chromosomal positions of DMRs located to genes. When a DMR included more than one gene it was mapped to the gene with most probes or if equal to the gene with the strongest signal.
Differences between groups were tested with t-test, Pearson x 2 square test or by Fisher exact 2-tailed test as appropriate. Spearman rank correlation was used to test for correlation between expression and methylation. Log-rank test was used to compare survival distributions.

Disclosure of potential conflicts of interest
No potential conflicts of interest were disclosed.