Stage 4S neuroblastoma tumors show a characteristic DNA methylation portrait

ABSTRACT

Stage 4S neuroblastoma (NB) is a special type of NB found in infants with metastases at diagnosis and is associated with an excellent outcome due to its remarkable capacity to undergo spontaneous regression. As genomics have not been able to explain this intriguing clinical presentation, we here aimed at profiling the DNA methylome of stage 4S NB to better understand this phenomenon. To this purpose, differential methylation analyses between International Neuroblastoma Staging System (INSS) stage 4S, stage 4 and stage 1/2 were performed, using methyl-CpG-binding domain (MBD) sequencing data of 14 stage 4S, 14 stage 4, and 13 stage 1/2 primary NB tumors (all MYCN non-amplified in order not to confound results). Stage 4S-specific hyper- and hypomethylated promoters were determined and further characterized for genomic localization and function by cytogenetic band enrichment, gene set enrichment, transcription factor target enrichment and differential RNA expression analyses. We show that specific chromosomal locations are enriched for stage 4S differentially methylated promoters and that stage 4S tumors show characteristic hypermethylation of specific subtelomeric promoters. Furthermore, genes involved in important oncogenic pathways, in neural crest development and differentiation, and in epigenetic processes are differentially methylated and expressed in stage 4S tumors. Based on these findings, we describe new biological mechanisms possibly contributing to the stage 4S-specific tumor biology and spontaneous regression. In conclusion, this study is the first to describe the highly characteristic stage 4S DNA methylome. These findings will open new avenues to further unravel the NB pathology in general and stage 4S disease specifically.

KEYWORDS:

• DNA methylation
• methyl-CpG-binding domain (MBD) sequencing
• neuroblastoma
• stage 4S (MS)
• spontaneous regression

Abbreviations

DMR	=	differentially methylated region
GSEA	=	gene set enrichment analysis
INRGSS	=	International Neuroblastoma Risk Group Staging System
INSS	=	International Neuroblastoma Staging System
MBD	=	methyl-CpG-binding domain
NB	=	neuroblastoma
NGF	=	nerve growth factor

Introduction

Neuroblastoma (NB) is a childhood tumor that originates from precursor cells of the sympathetic nervous system and is a heterogeneous disease with prognosis ranging from long-term survival to fatal outcome. NB tumor staging systems have been established, based on clinical and biological parameters (age of the patient, tumor stage, MYCN amplification, DNA index, and histopathology), and are being used for therapeutic stratification. The special stage IV (IV-S) was first described in the staging system of Evans et al., introduced in 1971, and defined patients who would otherwise have been stage I or II (with localized primary tumors), but who had disseminated disease confined to liver, skin or bone marrow, or any combination of these, without involvement of the skeleton. Remarkably, these patients did not have the dismal prognosis of other metastatic NB, as most of these tumors are characterized by spontaneous regression without any therapeutic intervention.^1,2 Approximately 20 y later, the International Neuroblastoma Staging System (INSS) refined the criteria of this tumor stage, renamed as stage 4S, and limited it to children younger than 12 months of age at diagnosis.³ In 2005, this age limit was extended to 18 months in the International Neuroblastoma Risk Group Staging System (INRGSS), as it was demonstrated that an age cut-off of 18 months resulted in a stratification that optimized the prognostic contribution of the age at diagnosis.⁴ In the INRGSS, this special tumor stage is indicated as MS.⁵

Several genomic alterations have been identified in NB, but so far no specific changes, apart from hyperdiploidy and near triploidy with numerical chromosome alterations common to localized tumors without MYCN amplification, have been associated with stage 4S NB. This indicates that other molecular mechanisms, such as epigenetic alterations, may play a role in stage 4S NB.⁶ In this study, we therefore specifically focused on identifying patterns of DNA methylation that are characteristic of this particular NB tumor stage.

Results

Stage 4S tumors demonstrate a unique promoter methylation portrait

In order to profile the stage 4S NB methylome, we re-analyzed methyl-CpG-binding domain (MBD) sequencing data of 41 MYCN non-amplified primary NBs,^7,8 including INSS stage 1 (n = 7), stage 2 (n = 6), stage 4 (n = 14), and stage 4S (n = 14) tumors. As patients with stage 4S NB, in sharp contrast to stage 4 patients, show excellent outcome despite widespread metastases at diagnosis, we first performed a differential methylation analysis between these tumor stages. In total, 5,914 and 10,974 promoters were identified as being hyper- or hypomethylated in stage 4S tumors, respectively (Supplemental Table 1). Next, stage 4S tumors were also compared to stage 1/2 tumors, which are also associated with a favorable patient outcome, but patients with stage 1/2 tumors do not present with metastatic disease. Here, 1,288 and 3,854 promoters were hyper- and hypomethylated in stage 4S tumors, respectively (Supplemental Table 1). Subsequently, we determined the number of stage 4S-specific hyper- and hypomethylated promoters by determining the overlap between these 2 differential methylation analyses and by excluding promoters differentially methylated in the stage 1/2 vs. stage 4 comparison. As such, we identified 393 hypermethylated and 1,150 hypomethylated promoters that typify the stage 4S methylome (Fig. 1A; Supplemental Table 1). Similarly, also for stage 4 and stage 1/2 NB, specific promoter methylation patterns were determined (Fig. 1A; Supplemental Table 1). Principal component analysis using the normalized read counts of the promoters included in Fig. 1A visualizes the NB stage-specific promoter methylation patterns (Fig. 1B). A dense clustering is observed for stage 4S tumors and this tumor stage seems most similar to stage 1/2.

Figure 1. MBD sequencing analysis of MYCN non-amplified stage 4S, stage 4, and stage 1/2 tumors identifies stage-specific promoter methylation portraits. A. Diagram of the promoters identified as hypermethylated in the respective sample groups. For example, 393 and 1,150 promoters are specifically hyper- and hypomethylated in stage 4S, respectively. B. Principal component (PC) analysis using the normalized read counts of the promoters included in the diagram. Shown are the data for the first 2 PCs. PC1 explains 29.9% of the variability of the data and PC2 10.1%.

Specific chromosomal locations are enriched in stage 4S differentially methylated promoters

Next, we determined the chromosomal distribution of the promoters contributing to the specific stage 4S methylation portrait. The percentages of hyper- and hypomethylated promoters located on the different chromosomes is indicated in Fig. 2A and Fig. 2B, respectively. Importantly, as chromosomes differ in length and gene density, these chromosome percentages were calculated based on the number of identified promoters corrected for the total number of promoters located on the corresponding chromosome. Stage 4S hypermethylated promoters are overrepresented (P < 0.05) on chr17 and chr19, and underrepresented (P < 0.05) on chr1, chr3, chr4, chr5, chr7, chr10, chr11, chr14, chr15, chr18, chr20, and chr21 (Fig. 2A). Stage 4S hypomethylated promoters are overrepresented (P < 0.05) on chr3, chr7, chr8, chr13, and chr18, and underrepresented (P < 0.05) on chr2, chr6, chr12, chr16, chr17, chr19, and chr22 (Fig. 2B). To further pinpoint specific chromosomal regions enriched in stage 4S differentially methylated promoters, cytogenetic band enrichment analyses in WebGestalt^9,10 were performed (Fig. 2C and Supplemental Table 2). These analyses indicate that stage 4S hyper- and hypomethylated promoters are not randomly distributed across the genome, but that specific chromosomal cytogenetic bands are significantly enriched in stage 4S differentially methylated promoters (Fig. 2).

Figure 2. Stage 4S differentially methylated promoters are not randomly distributed across the genome. A. Percentage of stage 4S hypermethylated promoters per chromosome, based on promoter density corrected numbers. Significantly under- and overrepresented chromosomes are indicated with a dark dot. The vertical line represents the percentage if the stage 4S hypermethylated promoters would be randomly distributed across the genome (null hypothesis). B. Percentage of stage 4S hypomethylated promoters per chromosome, based on promoter density corrected numbers. Significantly under- and overrepresented chromosomes are indicated with a dark dot. The vertical line represents the percentage if the stage 4S hypomethylated promoters would be randomly distributed across the genome (null hypothesis). C. Significantly enriched cytogenetic locations of stage 4S hyper- and hypomethylated promoters. A list of genes located in the enriched regions can be found in Supplemental Table 2.

Specific subtelomeres are enriched in stage 4S hypermethylated promoters

Telomerase activity upregulation and alternative lengthening of telomeres are key mechanisms in counteracting telomere shortening in NB and contribute to tumor cell immortalization.^11-14 As DNA methylation of subtelomeres represents an additional way to suppress telomere elongation^15,16 and as tumor regression (and thus loss of immortality) is a hallmark of stage 4S NB, we specifically analyzed the subtelomeric methylation patterns of the different NB tumor stages. To this purpose, we first determined the proportions of promoters of the different stage-specific methylation portraits that are located in subtelomeres (1 Mb or 2 Mb regions proximal to the telomeres). Fig. 3A and B illustrate that stage 4S hypermethylated and stage 4 hypomethylated promoters have a significantly larger proportion of 1 Mb and 2 Mb subtelomeric promoters (P < 0.001), compared to stage 1/2 and stage 4 hypermethylated and stage 4S hypomethylated promoters. Stage 1/2 hypomethylated promoters also comprise a larger proportion of 2 Mb subtelomeric promoters (P < 0.001; Fig. 3B), compared to stage 1/2 and stage 4 hypermethylated and stage 4S hypomethylated promoters. Next, using the gene lists of the differential methylation analyses, gene set enrichment analyses (GSEA) were performed to evaluate whether specific subtelomeres are enriched in differentially methylated promoters (Table 1). These analyses firmly demonstrate that multiple subtelomeres are enriched in stage 4S hypermethylated promoters compared to stage 4 and stage 1/2 tumors. Remarkably, not a single subtelomere is enriched in stage 4 hypermethylated promoters. Table 1 further shows that the chr5p 1 Mb and 2 Mb subtelomeres are enriched in promoters hypermethylated in stage 4S compared to both stage 1/2 and stage 4, and in stage 1/2 compared to stage 4. TERT, the gene encoding the catalytic subunit of telomerase, is located on the chr5p 2 Mb subtelomere, and its promoter is hypermethylated in stage 4S compared to stage 4 (Supplemental Table 1). Using publically available mRNA expression data,¹⁷ we could further demonstrate that TERT is expressed at lower levels in stage 4S compared to stage 4 (P < 0.001; log₂ fold change = 0.82).

Figure 3. Stage 4S hypermethylated, and stage 4 and 1/2 hypomethylated promoters are frequently located on subtelomeres. Depicted are the percentages of the stage-specific differentially methylated promoters located in the 1 Mb (A) and 2 Mb (B) subtelomeres. Stage 4S hypermethylated and stage 4 hypomethylated promoters have a significantly larger proportion of 1 Mb and 2 Mb subtelomeric promoters, compared to stage 1/2 and stage 4 hypermethylated and stage 4S hypomethylated promoters. Stage 1/2 hypomethylated promoters also comprise a larger proportion of 2 Mb subtelomeric promoters, compared to stage 1/2 and stage 4 hypermethylated and stage 4S hypomethylated promoters.

Table 1. One third of subtelomeres are enriched in stage 4S hypermethylated promoters.

		Enriched 1Mb subtelomeres				Enriched 2 Mb subtelomeres
Comparison		Subtelomere	NES	Nom P-value	FDR q-value	Subtelomere	NES	Nom P-value	FDR q-value
Stage 4S vs. stage 4	Hypermethylated stage 4S	chr4p	1.58	0.037	0.040	chr1p	2.21	<0.001	<0.001
		chr5p	2.29	<0.001	<0.001	chr4p	1.92	<0.001	0.002
		chr10q	1.87	<0.001	0.002	chr4q	1.47	0.041	0.045
		chr13q	2.12	<0.001	<0.001	chr5p	2.60	<0.001	<0.001
		chr16p	2.32	<0.001	<0.001	chr5q	1.49	0.011	0.042
		chr17q	1.68	0.006	0.022	chr8q	1.35	0.022	0.083
		chr19p	2.28	<0.001	<0.001	chr9q	1.58	<0.001	0.029
						chr10q	1.78	<0.001	0.003
						chr13q	1.97	<0.001	0.001
						chr16p	2.50	<0.001	<0.001
						chr16q	2.17	<0.001	<0.001
						chr17q	1.65	<0.001	0.020
						chr19p	2.80	<0.001	<0.001
						chr20q	1.83	<0.001	0.004
	Hypermethylated stage 4			None				None
Stage 4S vs. stage 1/2	Hypermethylated stage 4S	chr4q	2.04	<0.001	<0.001	chr4q	2.12	<0.001	<0.001
		chr5p	1.59	0.022	0.030	chr5p	1.81	<0.001	0.006
		chr9q	2.30	<0.001	<0.001	chr5q	1.96	<0.001	<0.001
		chr12q	1.60	0.027	0.031	chr8q	1.79	<0.001	0.006
		chr13q	1.72	<0.001	0.013	chr10q	1.73	<0.001	0.008
		chr17p	1.81	0.008	0.008	chr13q	1.75	<0.001	0.009
		chr17q	2.22	<0.001	<0.001	chr20q	3.56	<0.001	<0.001
		chr20q	2.19	<0.001	<0.001	chr21q	1.64	<0.001	0.018
		chr22q	1.53	<0.001	0.044	chr22q	1.28	0.047	0.109
	Hypermethylated stage 1/2	chr1p	1.38	0.035	0.141			None
Stage 1/2 vs. stage 4	Hypermethylated stage 1/2	chr5p	2.29	<0.001	<0.001	chr4p	1.49	0.023	0.054
						chr5p	1.92	<0.001	<0.001
						chr10p	1.49	0.035	0.059
	Hypermethylated stage 4			None				None

For each of the 3 sample group comparisons, the subtelomeres significantly enriched for promoters hypermethylated in each of the sample groups are indicated, as well as their corresponding normalized enrichment score (NES), nominal P-value (nom P-value < 0.05), and false discovery rate q-value (FDR q-value < 0.25).

Genes involved in important oncogenic pathways, neural crest development and differentiation, and epigenetic processes are differentially methylated and expressed in stage 4S tumors

To functionally characterize the stage 4S differentially methylated promoters, we first examined whether these promoters share common DNA motifs, by analyzing transcription factor target enrichment in WebGestalt.^9,10 In total, 24 and 320 motifs were significantly enriched for stage 4S hyper- and hypomethylated promoters, respectively (Supplemental Table 3), of which the top 10 (ranked based on the ratio of enrichment) are shown in Table 2. Motifs enriched in stage 4S differentially methylated promoters represent binding sites for transcription factors implicated in controlling cell cycle, apoptosis, neural crest development, and neural differentiation (Table 2). Further, we performed GSEA on the gene lists of the differential methylation analyses, using the oncogenic signatures database containing gene sets that represent signatures of cellular pathways dysregulated in cancer.¹⁸ Enriched gene sets were only identified for promoters hypomethylated in stage 4S compared to stage 4, in stage 4S compared to stage 1/2, and in stage 1/2 compared to stage 4 (Supplemental Table 4). Of note, 19 gene sets (Table 3) were specifically enriched in promoters hypomethylated in stage 4S compared to both stage 1/2 and stage 4, including gene sets involved in downstream signaling of important oncogenes and epigenetic actors. Transcription factor target gene sets (Table 2) and the leading edge genes of the oncogenic signatures gene sets (Table 3) significantly enriched in stage 4S differentially methylated promoters were further tested for differential expression, using the global test¹⁹ and publically available mRNA expression data.¹⁷ In total, 25 of the 39 (64.10%) gene sets were identified as differentially expressed in stage 4S compared to both stage 4 and stage 1/2 (Table 2 and Table 3). Finally, using covariates plots¹⁹ and differential gene expression analyses on the individual gene level, we evaluated which genes contributed the most to the significant global test results (Supplemental Table 5 and Table 4).

Table 2. Specific transcription factor targets are enriched in stage 4S differentially methylated promoters.

Hypermethylated in stage 4S
		Enrichment analysis						Differential expression analysis
Transcription factor target	Transcription factor	C	E	O	Genes	R	q	Stage 4S vs. stage 4	Stage 4S vs. stage 1/2
								q	q
CCGNMNNTNACG		75	0.28	3	HMG20B, HSPD1, FOSB	10.58	0.035	0.198	0.777
NTNGCGTGNNN	AHR	209	0.79	6	TSSK6, SLC9A5, FHOD1, FOSB, CIRBP, LPHN1	7.60	0.012	0.051	0.292
TTCNRGNNNNTTC		150	0.57	4	MKNK2, HSPD1, FOSB, CCT4	7.06	0.034	0.225	0.944
TAAWWATAG	MEF2A	163	0.62	4	EIF5A, SLC9A5, ALDOA, JSRP1	6.49	0.036	< 0.001	< 0.001
SGCGSSAAA	E2F1 TFDP1 RB1	167	0.63	4	MCM2, MCM7, SLC9A5, FHOD1	6.34	0.038	< 0.001	0.001
TGAMCTTTGMMCYT	HNF4A	260	0.98	6	PRKRA, NDUFA3, PRKCSH, RTN2, C4A, ALDOA	6.11	0.017	0.033	0.147
TTTSGCGS	E2F1	228	0.86	5	GPN3, MCM2, MCM7, SLC9A5, FHOD1	5.80	0.034	< 0.001	0.001
NCSCGCSAAAN	E2F1 TFDP1	234	0.88	5	EIF5A, TNPO2, GPN3, MCM2, MCM7	5.65	0.034	0.053	0.044
GCCATNTTN	YY1	238	0.90	5	NCOR1, PIGL, PRKCSH, PTBP1, LPHN1	5.56	0.034	0.121	0.003
NKCGCGCSAAAN	E2F1 TFDP1	239	0.90	5	EIF5A, TNPO2, MCM2, MCM7, PRKCSH	5.54	0.034	0.053	0.044
Hypomethylated in stage 4S
		Enrichment analysis						Differential expression analysis
Transcription factor target	Transcription factor	C	E	O	Genes	R	q	Stage 4S vs. stage 4	Stage 4S vs. stage 1/2
								q	q
NNTGTTACTAAAAATAGAAMNN		25	0.32	3	SLCO5A1, PPARGC1A, MGST3	9.27	0.011	0.045	0.213
NTGCGTGGGCGK	EGR3	86	1.11	7	CELF4, GPR176, ATP6V1C1, EPHB1, KCNQ5, UBE2W, CLTC	6.29	< 0.001	< 0.001	< 0.001
CARAACTAGGNCAAAGGTCA	PPARA	37	0.48	3	SS18, MARCH10, PPARGC1A	6.27	0.026	0.046	0.379
KMCATNNWGGA		87	1.13	7	CELF4, PDE4D, ADAMTSL1, RUNX1T1, OSBPL6, ZFAT, ARNTL	6.22	< 0.001	0.530	< 0.001
TGATTTRY	GFI1	291	3.77	23	RUNX1T1, FOXN3, WDPCP, AGR2, ALPK2, MYOT, SND1, CREB5, CDK14, IGF2BP3, FOXP2, NRXN3, MIR137HG, PTEN, C12orf42, ALDH1A2, JAZF1, ZNF521, MYBPC1, ZBTB20, PPARGC1A, MAPK10, RIMS2	6.11	< 0.001	< 0.001	< 0.001
NNGAATATKCANNNN	POU2F1	211	2.73	16	ADAMTSL1, NPAS3, RUNX1T1, PHLPP1, NRXN3, FOXP2, MARCH10, PPFIA2, ZNF521, ZBTB20, LPL, PDE4D, CACNA1D, PPARGC1A, RGS8, ARPP21	5.86	< 0.001	< 0.001	< 0.001
CTBATTTCARAAW	CEBPG	252	3.26	19	BRAF, RUNX1T1, NR6A1, FOXN3, PTK2, CREB5, WDR64, KCNQ5, CXCL9, FOXP2, SPAG9, PPFIA2, ALDH1A2, ST18, ZBTB20, C8orf46, SLIRP, PPARGC1A, PABPC1	5.83	< 0.001	< 0.001	0.007
NNNNNNKCTAWAAATAGMNNNN		226	2.92	17	LARS, SLC8A3, CELF4, ALPK2, CKMT2, KCNQ5, CACNA2D3, ARHGEF38, ARNTL, FBXO40, TSC22D1, TNNT2, PPARGC1A, C8orf46, ADRA1A, MLLT3, ARPP21	5.81	< 0.001	0.040	0.067
CNGTAWNTG	MSX1	173	2.24	13	FOXP2, BCAR3, KCNAB1, ADAMTSL1, ELMO1, C12orf42, RUNX1T1, NR6A1, ARNTL, CBLN2, SND1, INHBA, IGF2BP3	5.81	< 0.001	0.340	0.010
AGATAAGATAA	EVI1	55	0.71	4	CREB5, CACNA1D, PPARGC1A, RUNX1T1	5.62	0.014	< 0.001	0.324

Shown are the top 10 significantly (q < 0.05) enriched transcription factor targets (motifs) in stage 4S hyper- and hypomethylated promoters. For each enriched motif, matching transcription factor annotation is shown (if available), as well as the number of genes in the category (C), the expected number in the category (E), the number of genes hyper- or hypomethylated in stage 4S in the category (O), the annotation for the genes in O, the ratio of enrichment (R), and the q-value. Global test q-values of the differential expression analyses are also indicated.

Table 3. Oncogenic signatures are enriched in stage 4S hypomethylated promoters.

	Enrichment analysis						Differential expression analysis
	Stage 4S vs. stage 4			Stage 4S vs. stage 1/2			Stage 4S vs. stage 4	Stage 4S vs. stage 1/2
Gene set	NES	Nom P-value	FDR q-value	NES	Nom P-value	FDR q-value	q	q	Description of gene set
BCAT_BILD_ET_AL_DN	1.38	0.033	0.053	1.47	0.002	0.026	0.002	0.205	Genes downregulated in primary epithelial breast cancer cell culture overexpressing activated CTNNB1 (PMID: 16273092)
BMI1_DN.V1_DN	1.30	0.025	0.088	1.30	0.002	0.102	0.002	< 0.001	Genes downregulated in DAOY medulloblastoma cells upon knockdown of BMI1 by RNAi (PMID: 17452456)
CAHOY_NEURONAL	1.49	0.001	0.042	1.27	0.021	0.131	< 0.001	< 0.001	Genes upregulated in neurons (PMID: 18171944)
CRX_NRL_DN.V1_DN	1.38	0.006	0.056	1.23	0.030	0.182	< 0.001	< 0.001	Genes downregulated in retina cells from CRX and NRL double knockout mice (PMID: 17653270)
EGFR_UP.V1_DN	1.40	0.001	0.056	1.26	0.007	0.131	< 0.001	< 0.001	Genes downregulated in MCF7 breast cancer cells positive for ESR1 and engineered to express ligand-activated EGFR (PMID: 11888208)
EIF4E_DN	1.31	0.035	0.076	1.42	< 0.001	0.036	< 0.001	< 0.001	Genes downregulated in HMEC primary mammary epithelium cells upon overexpression of EIF4E (PMID: 17638893)
ERB2_UP.V1_DN	1.25	0.041	0.145	1.41	< 0.001	0.038	< 0.001	0.005	Genes downregulated in MCF7 breast cancer cells positive for ESR1 and engineered to express ligand-activated ERBB2 (PMID: 11888208)
KRAS.AMP.LUNG_UP.V1_UP	1.37	0.006	0.063	1.30	0.006	0.113	0.003	< 0.001	Genes upregulated in epithelial lung cancer cell lines overexpressing KRAS (PMID: 19847166)
KRAS.KIDNEY_UP.V1_UP	1.66	< 0.001	0.004	1.50	< 0.001	0.021	< 0.001	< 0.001	Genes upregulated in epithelial kidney cancer cell lines overexpressing an oncogenic form of KRAS
MEK_UP.V1_DN	1.47	< 0.001	0.044	1.37	< 0.001	0.054	< 0.001	< 0.001	Genes downregulated in MCF7 breast cancer cells positive for ESR1; MCF7 cells stably overexpressing constitutively active MAP2K1 (PMID: 16585219)
MEL18_DN.V1_DN	1.36	0.004	0.062	1.20	0.040	0.244	< 0.001	0.006	Genes downregulated in DAOY medulloblastoma cells upon knockdown of PCGF2 gene by RNAi (PMID: 17452456)
PIGF_UP.V1_UP	1.31	0.009	0.076	1.38	< 0.001	0.053	< 0.001	0.005	Genes upregulated in HUVEC endothelium cells by treatment with PIGF (PMID: 15516835)
PRC2_EDD_UP.V1_DN	1.30	0.014	0.088	1.32	< 0.001	0.106	< 0.001	< 0.001	Genes downregulated in TIG3 fibroblast cells upon knockdown of EED (PMID: 16618801)
PRC2_EZH2_UP.V1_DN	1.33	0.003	0.071	1.21	0.034	0.209	0.021	< 0.001	Genes downregulated in TIG3 fibroblast cells upon knockdown of EZH2 (PMID: 16618801)
RAF_UP.V1_DN	1.39	< 0.001	0.056	1.30	< 0.001	0.105	0.002	< 0.001	Genes downregulated in MCF7 breast cancer cells positive for ESR1; MCF7 cells stably overexpressing constitutively active RAF1 (PMID: 16585219)
RAF_UP.V1_UP	1.33	0.007	0.068	1.25	0.009	0.139	0.002	< 0.001	Genes upregulated in MCF7 breast cancer cells positive for ESR1; MCF7 cells stably overexpressing constitutively active RAF1 (PMID: 16585219)
STK33_NOMO_UP	1.21	0.040	0.204	1.27	< 0.001	0.127	0.007	< 0.001	Genes upregulated in NOMO1 acute myeloid leukemia cells after knockdown of STK33 by RNAi (PMID: 19490892)
STK33_SKM_UP	1.27	0.007	0.115	1.28	< 0.001	0.120	0.002	< 0.001	Genes upregulated in SKM1 acute myeloid leukemia cells after knockdown of STK33 by RNAi (PMID: 19490892)
STK33_UP	1.33	< 0.001	0.070	1.31	< 0.001	0.109	0.004	0.004	Genes upregulated in NOMO1 and SKM1 acute myeloid leukemia cells after knockdown of STK33 by RNAi (PMID: 19490892)

For each gene set significantly enriched in promoters hypomethylated in stage 4S compared to stage 4 and stage 1/2, the corresponding normalized enrichment score (NES), nominal P-value (nom P-value < 0.05) and false discovery rate q-value (FDR q-value < 0.25) are shown. Global test q-values of the differential expression analyses are also indicated, as well as a description of the gene set (with PubMed ID of the source publication, as described in the Molecular Signatures Database).¹⁸

Table 4. Multiple stage 4S differentially methylated genes are also differentially expressed.

Hypermethylated in stage 4S
		Stage 4S vs. stage 4			Stage 4S vs. stage 1/2
Enriched gene set	Gene	q	logFC	Overexpression in	q	logFC	Overexpression in
E2F1, E2F1_TFDP1_RB1 and MEF2A	SLC9A5	<0.001	0.833	stage 4	0.002	0.500	stage 1/2
Hypomethylated in stage 4S
		Stage 4S vs. stage 4			Stage 4S vs. stage 1/2
Enriched gene set	Gene	q	logFC	Overexpression in	q	logFC	Overexpression in
CAHOY_NEURONAL	DYNC1I1	0.010	1.035	stage 4	< 0.001	1.517	stage 1/2
	HS3ST2	0.021	0.680	stage 4	0.031	0.589	stage 1/2
	NSF	0.006	0.470	stage 4	< 0.001	0.538	stage 1/2
	GPR22	< 0.001	1.272	stage 4S	0.016	0.847	stage 4S
	CALB1	< 0.001	1.718	stage 4S	0.048	0.850	stage 4S
CRX_NRL_DN.V1_DN	GPR124	0.011	0.471	stage 4S	0.050	0.344	stage 4S
	PVRL3	< 0.001	1.153	stage 4S	0.009	0.719	stage 4S
EGFR_UP.V1_DN	C9orf3	< 0.001	0.833	stage 4S	0.011	0.574	stage 4S
	DERA	< 0.001	0.897	stage 4S	0.023	0.451	stage 4S
EGR3	ATP6V1C1	0.002	0.380	stage 4	< 0.001	0.444	stage 1/2
	CLTC	0.008	0.391	stage 4	0.002	0.419	stage 1/2
	GPR176	0.031	0.184	stage 4	0.018	0.185	stage 1/2
	UBE2W	0.048	0.177	stage 4	0.031	0.177	stage 1/2
EIF4E_DN	CLTC	0.008	0.391	stage 4	0.002	0.419	stage 1/2
	NPEPPS	0.002	0.348	stage 4	0.004	0.305	stage 1/2
	TP53BP1	0.027	0.287	stage 4	0.010	0.304	stage 1/2
	LMO7	0.001	0.490	stage 4S	0.028	0.323	stage 4S
	SH3YL1	0.004	0.428	stage 4S	0.049	0.282	stage 4S
ERB2_UP.V1_DN	NPEPPS	0.002	0.348	stage 4	0.004	0.305	stage 1/2
	C9orf3	< 0.001	0.833	stage 4S	0.011	0.574	stage 4S
	LAMB1	< 0.001	0.690	stage 4S	0.022	0.354	stage 4S
KRAS.AMP.LUNG_UP.V1_UP	DYNC1I1	0.010	1.035	stage 4	< 0.001	1.517	stage 1/2
	C8B	0.038	0.340	stage 4S	0.006	0.408	stage 4S
	GCKR	0.027	0.240	stage 4S	0.032	0.215	stage 4S
KRAS.KIDNEY_UP.V1_UP	DYNC1I1	0.010	1.035	stage 4	< 0.001	1.517	stage 1/2
	ESRRG	0.021	0.475	stage 4	0.018	0.449	stage 1/2
	CALB1	< 0.001	1.718	stage 4S	0.048	0.850	stage 4S
MEK_UP.V1_DN	NPEPPS	0.002	0.348	stage 4	0.004	0.305	stage 1/2
	C9orf3	< 0.001	0.833	stage 4S	0.011	0.574	stage 4S
	TCF12	0.001	0.370	stage 4S	0.015	0.261	stage 4S
MEL18_DN.V1_DN	DIO2	0.048	0.460	stage 4	0.019	0.497	stage 1/2
PIGF_UP.V1_UP	LRRC17	< 0.001	0.568	stage 4S	0.031	0.283	stage 4S
PRC2_EDD_UP.V1_DN	CBX8	< 0.001	0.815	stage 4	0.002	0.588	stage 1/2
RAF_UP.V1_DN	C9orf3	< 0.001	0.833	stage 4S	0.011	0.574	stage 4S
	PALLD	0.042	0.330	stage 4S	0.019	0.347	stage 4S
	TCF12	0.001	0.370	stage 4S	0.015	0.261	stage 4S
RAF_UP.V1_UP	CABYR	0.003	0.462	stage 4	< 0.001	0.589	stage 1/2
	GATA2	0.032	0.711	stage 4S	0.012	0.760	stage 4S
STK33_NOMO_UP	KIAA1467	< 0.001	0.560	stage 4	0.002	0.447	stage 1/2
	NPEPPS	0.002	0.348	stage 4	0.004	0.305	stage 1/2
	UBE2W	0.048	0.177	stage 4	0.031	0.177	stage 1/2
STK33_SKM_UP	UBE2W	0.048	0.177	stage 4	0.031	0.177	stage 1/2
	PALLD	0.042	0.330	stage 4S	0.019	0.347	stage 4S
STK33_UP	UBE2W	0.048	0.177	stage 4	0.031	0.177	stage 1/2

For each gene set differentially expressed in stage 4S compared to stage 4 and stage 1/2 (Table 2 and Table 3), differentially expressed individual genes are shown. Q-values are indicated, as well as the log₂ fold change (logFC) and the group in which the gene is overexpressed.

Discussion

Genome-wide characterization studies in NB often combine patients with stage 4S or other low-risk tumors into one group for comparison to high-risk patients, as low-risk patients generally have an excellent outcome, which is in sharp contrast to high-risk patients who often show fatal progression. Patients with stage 4S tumors demonstrate a typical pattern of metastatic spread at diagnosis and have the capacity to undergo spontaneous regression, and these characteristics are not exhibited by stage 1/2 tumors. As genomics alone have not been able to explain this intriguing clinical presentation,⁶ we aimed to characterize the DNA methylation pattern of stage 4S NB, in order to reveal new biological insights into this special type of NB.

By performing differential methylation analyses between stage 4S and stage 1/2 and stage 4, we could clearly define a unique stage 4S promoter methylation portrait. We further showed that stage 4S-specific hyper- and hypomethylated promoters are located on specific chromosomal locations. Non-random chromosomal distribution of regions differentially methylated between tumor and normal samples,^20,21 and between different tumor stages,^21,22 has been previously reported for other cancer types. In this study, we show that specific chromosomal cytogenetic bands are enriched in stage 4S differentially methylated promoters and, most importantly, that stage 4S tumors show characteristic hypermethylation of specific subtelomeric promoters. Telomere length and its regulatory mechanisms are well-studied in NB.^11-14,23-25 Most stage 4S tumors have low telomerase activity or short telomeres,^11,12,23 and it has been suggested that the association between stage 4S and spontaneous regression might be related to differences in telomere length regulation.⁶ While subtelomeric DNA methylation has been shown to suppress telomere elongation in other cancer types,^15,16 such telomere length regulatory mechanism is unexplored in NB. In this study, we have shown that specific subtelomeric promoters are hypermethylated in stage 4S tumors compared to stage 1/2 and stage 4 tumors, which suggests that subtelomeric DNA methylation may represent an additional mechanism by which telomere length and spontaneous regression are controlled in NB. Therefore, the role of DNA methylation in telomere length regulation in stage 4S NB should definitely be further explored.

Further, we have shown that multiple transcription factor target gene sets are enriched in stage 4S differentially methylated promoters. Methylation of cytosines within transcription factor recognition motifs has been shown to cause structural changes influencing DNA-protein interactions (either in a positive or negative manner) and transcription factor activity.²⁶ As such, DNA methylation serves as a regulatory mechanism for transcription factor activity. Therefore, the transcription factor target gene sets enriched in stage 4S differentially methylated promoters may indicate differential activity of the respective transcription factors in stage 4S tumors compared to stage 1/2 and stage 4 tumors. Of note, this assumption is reinforced by the fact that several of the identified transcription factors of which multiple target genes are differentially methylated in stage 4S show proven dependency on the DNA methylation status of the binding site for transcription factor-DNA interaction. For example, DNA methylation inhibits the binding of E2F1 and CEBPG.^26-28 Additionally, the differentially methylated target gene sets are differentially expressed in stage 4S NB, further strengthening the possibility that DNA methylation affects the activity of these transcription factors. Of note, E2F1 also binds to the TERT promoter,²⁹ which is hypermethylated in stage 4S compared to stage 4, and TERT is expressed at lower levels in stage 4S compared to stage 4. This indicates that TERT DNA methylation might also regulate telomerase activity in stage 4S NB.

Furthermore, our transcription factor target enrichment analyses particularly demonstrate that the DNA methylation portrait of stage 4S NB is dominated by differential methylation of target genes of transcription factors involved in neural crest development and (sympathetic) neural differentiation (MSX1,^30,31 EVI1,^32,33 EGR3,³⁴ AHR,^35,36 MEF2A,^37,38 YY1,³⁹ PPARA,⁴⁰ POU2F1,⁴¹ and GFI1^42,43). These findings suggest that differences in the pathogenic mechanisms between stage 4S and stage 1/2 and stage 4 tumors are at least in part related to differences in transcriptional activity during neural crest development and differentiation, which might indicate that stage 4S, stage 1/2 and stage 4 tumors arise during different phases of neural crest cell development. This is in line with the previous postulation that stage 4S NB originate from defective premigratory neural crest cells, which explains its clinical symptomology, as stage 4S tumor sites mirror developmental neural crest cell migration.⁴⁴

One of the genes contributing the most to the significant differential expression of the target genes of E2F1 and MEF2A is SLC9A5. The SLC9A5 promoter is hypermethylated in stage 4S and these tumors show low SLC9A5 expression levels compared to stage 4 and stage 1/2 tumors. SLC9A5 is a member of the Na(+)/H(+) exchanger family and is abundantly expressed in neurons. In rat pheochromocytoma, a neuroendocrine tumor of the adrenal medulla, Diering et al. have shown that SLC9A5 modulates trafficking of TrkA between recycling endosomes and the plasma membrane, neurotrophin signaling, and neural differentiation.⁴⁵ The authors demonstrate that upon knockdown of SLC9A5 in these tumor cells, there is no difference in the total abundance of TrkA compared to control cells, suggesting that reduced expression of SLC9A5 has little effect on TrkA degradation or synthesis. However, although the overall TrkA protein levels did not change, the cell surface abundance of TrkA was reduced upon knockdown of SLC9A5.⁴⁵ Stage 4S NB tumors generally show high levels of TrkA and, depending on the presence or absence of nerve growth factor (NGF) in the tumor microenvironment, the tumors undergo neural differentiation or spontaneous regression, respectively. Therefore, neurotrophin deprivation (TrkA without NGF) and activation of apoptosis have been proposed to explain the phenomenon of spontaneous regression in NB.⁶ Given the crucial role of SLC9A5 in regulating the cellular distribution of TrkA and NGF signaling,⁴⁵ our findings suggest that, over time, stage 4S tumor cell surfaces might get depleted of TrkA (with an accompanying intracellular accumulation of TrkA) due to the low expression levels of SLC9A5, ultimately leading to apoptosis and regression.

Our gene set enrichment analyses also highlight differential methylation of genes involved in important oncogenic pathways, especially of genes implicated in EGFR and RAS signaling. Although the exact mechanisms by which DNA methylation contributes to differential regulation of these pathways in stage 4S specifically cannot be determined in our study, our findings indicate that it is crucial to separate stage 4S tumors from other low-stage (stage 1/2) tumors in studies focusing on NB biology, as this might reveal new insights into the stage 4S NB pathogenesis. The importance hereof is highlighted by the fact that all but one differentially methylated oncogenic signatures are also differentially expressed in stage 4S NB, confirming differential regulation of these oncogenic pathways in stage 4S on the one hand and stage 4 and stage 1/2 on the other hand. Finally, the enrichment analyses also illustrate that stage 4S tumors present differential activity of important epigenetic players, such as BMI1, PCGF2, EED, and EZH2. Potentially, differential regulation of epigenetic mechanisms might also contribute to the regression phenotype of stage 4S NB, as epigenetic alterations can more easily be reversed. Clearly, these epigenetic mechanisms need to be further studied in NB.

In conclusion, this study is the first to describe the stage 4S DNA methylome and demonstrates that stage 4S NB tumors show characteristic DNA methylation patterns compared to stage 1/2 and stage 4 tumors. Genomic localization and functional characterization analyses of this stage 4S DNA methylation portrait suggested new mechanisms that may contribute to the stage 4S-specific tumor biology and spontaneous regression. Future studies must confirm these observations and hypotheses.

Materials and methods

Neuroblastoma primary tumors

Sample collection for MBD sequencing (n = 41) and expression profiling (n = 180) is described in Decock et al.⁷ and Oberthuer et al, ¹⁷ respectively. All tumor samples were assigned to one of the 3 following subgroups, based on the INSS tumor stage³ and MYCN amplification status: (1) MYCN non-amplified stage 4S (MBD sequencing: n = 14; expression profiling: n = 28), (2) MYCN non-amplified stage 4 (MBD sequencing: n = 14; expression profiling: n = 46), or (3) MYCN non-amplified stage 1 or 2 (MBD sequencing: n = 13; expression profiling: n = 106). Detailed criteria for patient inclusion in the study are described in Decock et al.^8,46 MYCN amplified tumor samples were excluded. Detailed clinicobiological characteristics of the different sample cohorts are shown in Supplemental Table 6. The study was approved by the ethical committee of the Ghent University Hospital (approval number: B67020109912).

DNA methylome and transcriptome profiling

A detailed description of the generation and raw data analysis of the MBD sequencing data is provided in Decock et al.⁷ Data can be accessed through the Gene Expression Omnibus numbers GSE69243 and GSE69268. For the mRNA expression analyses, microarray data of 180 NB tumors were used.¹⁷

Bioinformatics and statistical analyses

Differential methylation analyses on the MBD sequencing data were performed as described in Decock et al.,⁸ using DESeq⁴⁷ and promoter region count matrices with data on all Ensembl transcripts (release 68). Promoters located on X and Y chromosomes were excluded to avoid gender-related confounding. To estimate the background signal and to be able to exclude possible falsely identified differentially methylated regions (DMRs) in the DESeq analyses, a filtering procedure based on 2 correction steps was applied on the resulting DMR lists. A first correction step consisted of excluding DMRs, which were also significantly identified in the same comparison on the input (non-MBD-enriched) sample data (input correction).⁸ Secondly, computational simulation analyses were performed, in which 100 random patient annotation swaps were executed to create 2 artificial sample groups on which the same statistics were applied as on the biologically relevant groups. Only biologically relevant DMRs with a π-value (= −ln pval * log₂ fold change)⁴⁸ that equals or outperforms the π-value of the corresponding simulation DMR were considered for further evaluation (swap correction). DESeq normalized read counts of the stage-specific DNA methylation portraits were used for principal component analysis in R (package factoextra).

Cytogenetic band enrichment analysis and transcription factor target enrichment analysis were performed using WebGestalt^9,10 with gene symbols as input list. The entire human genome was used as reference list. Benjamini-Hochberg multiple testing correction was applied⁴⁹ and only statistically enriched terms (q < 0.05) with at least 2 genes were considered.

Non-corrected transcript lists, including all transcripts included in the differential methylation analyses (preranked on the π-value) were used for GSEA.^18,50 To this purpose, Ensembl transcript IDs were converted to gene symbols. For genes with multiple transcripts, the data of the Ensembl transcript with the highest absolute π-value was retained. Tested gene set databases are the oncogenic signatures (version 5.1) of the Molecular Signatures Database¹⁸ and 2 in-house created subtelomere gene set lists. These two subtelomere gene set lists contain all genes in the 1 Mb or 2 Mb subtelomeres, which are defined as the 1 Mb or 2 Mb genomic regions²⁴ proximal to the telomeres (10 kb chromosomal ends; hg19). GSEA was also performed on the input data to exclude falsely identified enriched gene sets.

The sets of differentially methylated transcription factor targets and the leading edge genes of the enriched oncogenic signature gene sets were further tested for differential expression using the globaltest R package¹⁹ and publically available mRNA expression data.¹⁷ For significantly differentially expressed gene sets, covariates plots¹⁹ were generated to inspect the contributions of each gene to the test result. Differential expression analyses on the individual gene level were performed using the R package limma.⁵¹ Benjamini-Hochberg multiple testing correction was applied.⁴⁹

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Funding

This study is supported by the Fournier-Majoie Foundation (FFM), the Belgian National Lottery, the Flemish League Against Cancer (VLK), UGent GOA (grant number 01G01910) and the ISCIII Spain & FEDER (grant number PI14/01008). AD was a recipient of a grant of the Research Foundation Flanders (FWO) and an Emmanuel Van der Schueren research grant (VLK).