
==== Front
Commun Biol
Commun Biol
Communications Biology
2399-3642
Nature Publishing Group UK London

39079981
6596
10.1038/s42003-024-06596-6
Article
A MYCN-driven de-differentiation profile identifies a subgroup of aggressive retinoblastoma
Ryl Tatsiana 1
Afanasyeva Elena 1
http://orcid.org/0000-0002-6993-347X
Hartmann Till 2
Schwermer Melanie 1
http://orcid.org/0000-0001-6634-4306
Schneider Markus 1
Schröder Christopher 2
Wagemanns Maren 1
Bister Arthur 1
Kanber Deniz 3
Steenpass Laura 4
Schramm Kathrin 5678
Jones Barbara 5678
http://orcid.org/0000-0002-2036-5141
Jones David T. W. 5678
Biewald Eva 9
http://orcid.org/0000-0002-7872-1997
Astrahantseff Kathy 10
Hanenberg Helmut 1
Rahmann Sven 11
Lohmann Dietmar R. 3
http://orcid.org/0000-0001-7670-7529
Schramm Alexander 12
http://orcid.org/0000-0002-8138-0441
Ketteler Petra petra.ketteler@uk-essen.de

13
1 grid.410718.b 0000 0001 0262 7331 Department of Pediatric Hematology and Oncology, University Hospital Essen, Essen, Germany
2 https://ror.org/04mz5ra38 grid.5718.b 0000 0001 2187 5445 Algorithms for Reproducible Bioinformatics, Genome Informatics, Institute of Human Genetics, University Hospital Essen, University of Duisburg-Essen, Essen, Germany
3 https://ror.org/04mz5ra38 grid.5718.b 0000 0001 2187 5445 Institute of Human Genetics, University Hospital Essen, University Duisburg Essen, Essen, Germany
4 https://ror.org/02tyer376 grid.420081.f 0000 0000 9247 8466 Human and Animal Cell Lines, Leibniz Institute DSMZ German Collection of Microorganisms and Cell Cultures, 38124 Braunschweig, Germany
5 grid.510964.f Division of Pediatric Glioma Research, Hopp Children’s Cancer Center (KiTZ), Heidelberg, Germany
6 grid.5253.1 0000 0001 0328 4908 National Center for Tumor Diseases (NCT), NCT Heidelberg, a partnership between DKFZ and Heidelberg University Hospital, Heidelberg, Germany
7 https://ror.org/04cdgtt98 grid.7497.d 0000 0004 0492 0584 German Cancer Research Center (DKFZ), Heidelberg, Germany
8 https://ror.org/038t36y30 grid.7700.0 0000 0001 2190 4373 Department of Pediatric Hematology and Oncology, Heidelberg University Hospital, University of Heidelberg, Heidelberg, Germany
9 https://ror.org/04mz5ra38 grid.5718.b 0000 0001 2187 5445 Department of Ophthalmology, Medical Faculty, University of Duisburg-Essen, 45147 Essen, Germany
10 https://ror.org/001w7jn25 grid.6363.0 0000 0001 2218 4662 Department of Pediatric Oncology and Hematology, Charité – University Medicine Berlin, Berlin, Germany
11 https://ror.org/01jdpyv68 grid.11749.3a 0000 0001 2167 7588 Algorithmic Bioinformatics, Center for Bioinformatics Saar and Saarland University, Saarland Informatics Campus, Saarbrücken, Germany
12 https://ror.org/04mz5ra38 grid.5718.b 0000 0001 2187 5445 Laboratory for Molecular Oncology, Department of Medical Oncology, West German Cancer Center, University Hospital Essen, University of Duisburg-Essen, Essen, Germany
30 7 2024
30 7 2024
2024
7 91913 12 2023
17 7 2024
© The Author(s) 2024, corrected publication 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Retinoblastoma are childhood eye tumors arising from retinal precursor cells. Two distinct retinoblastoma subtypes with different clinical behavior have been described based on gene expression and methylation profiling. Using consensus clustering of DNA methylation analysis from 61 retinoblastomas, we identify a MYCN-driven cluster of subtype 2 retinoblastomas characterized by DNA hypomethylation and high expression of genes involved in protein synthesis. Subtype 2 retinoblastomas outside the MYCN-driven cluster are characterized by high expression of genes from mesodermal development, including NKX2-5. Knockdown of MYCN expression in retinoblastoma cell models causes growth arrest and reactivates a subtype 1-specific photoreceptor signature. These molecular changes suggest that removing the driving force of MYCN oncogenic activity rescues molecular circuitry driving subtype 1 biology. The MYCN-RB gene signature generated from the cell models better identifies MYCN-driven retinoblastoma than MYCN amplification and can identify cases that may benefit from MYCN-targeted therapy. MYCN drives tumor progression in a molecularly defined retinoblastoma subgroup, and inhibiting MYCN activity could restore a more differentiated and less aggressive tumor biology.

DNA methylation analysis identifies a MYCN-driven subgroup of aggressive retinoblastoma that may benefit from MYCN-targeted therapy.

Subject terms

Paediatric cancer
Cancer epigenetics
Eye cancer
Deutsche Forschungsgemeinschaft, TE 1187/2-1 Kinderaugenkrebsstiftung, KAKS20220422f DKTK, IFORES 2011issue-copyright-statement© Springer Nature Limited 2024
==== Body
pmcIntroduction

Retinoblastomas are retinal tumors originating from postmitotic cone photoreceptor precursor cells1–3. Most retinoblastomas are characterized by a biallelic inactivation of RB1, the RB transcriptional corepressor 1, and ~45% of patients with retinoblastoma carry a constitutional, pathogenic RB1 variant4. Early diagnosis of intraocular retinoblastoma leads to high survival rates in high-income countries, while metastatic retinoblastoma is rare because the natural boundaries of the eye prevent metastatic spread in early disease stages5. In low and middle income country, metastatic spread of retinoblastoma is more common and leads to significantly lower survival rates5. Intraocular retinoblastoma can be treated with enucleation of the eye or, especially in patients with bilateral disease, with a range of therapies to preserve the eye globe and vision. Biopsies of intraocular retinoblastoma are obsolete, so that samples from early disease that received eye-preserving therapy are scarce. Advances in and increasing use of eye-preserving therapies complicate retinoblastoma molecular genetic characterization and development of innovative targeted therapies, although there is hope that this may be overcome by serial liquid biopsies in the future6.

Recurrent genetic alterations other than RB1 are rare in retinoblastoma7–10, while copy number aberrations of chromosomes 1, 6, and 16 are characteristic11. DNA methylation-based clustering combined with gene expression profiling have distinguished two retinoblastoma subtypes with differences in clinical characteristics, genomic alterations, and, potentially, prognoses12–14. Subtype 1 retinoblastomas (equivalent to cluster 2 in Kooi et al.13) are characterized by few genetic and chromosomal alterations, mostly limited to RB1 loss and 6p gain, and show gene expression resembling signatures of maturing photoreceptor cells that Kooi et al. termed a “photoreceptorness signature”13 and Liu et al. termed a “cone marker signature”12. Subtype 2 retinoblastomas (equivalent to cluster 1 described by Kooi13) are genetically more heterogeneous, carry more genetic and chromosomal alterations than subtype 1 retinoblastomas and express high levels of neuronal markers and genes for stemness with low levels of photoreceptor-related genes12. High TFF1 expression has been described as a surrogate marker for subtype 2 retinoblastoma12,15. Preliminary clinical data suggest that subtype 2 retinoblastomas are more aggressive and at higher risk for metastasis12,15 and, consequently, may require higher treatment intensity. To tailor treatment accordingly in the future, further characterization of the molecular and biological characteristics of each subtype is warranted.

MYCN oncogene amplification has been observed in a small subgroup of retinoblastomas grouped into subtype 2 retinoblastoma12,13,16. The MYC protein family are basic helix-loop-helix (bHLH) transcription factors regulating proliferation, differentiation, cell-cycle progression, protein synthesis, metabolism, and apoptosis17–21. Amplifications of MYCN or MYC occur, respectively, in the childhood cancers, neuroblastoma, and medulloblastoma, and the resulting protein overexpression is associated with aggressive growth and poor clinical outcome22–24. In retinoblastoma, MYCN amplification (MYCNA) occurs in RB1−/−MYCNA retinoblastomas and the rare, RB1-proficient MYCNA retinoblastomas. RB1-proficient MYCNA retinoblastomas are currently considered a separate disease entity with very aggressive behavior, in which MYCN is thought to exclusively drive tumorigenesis25,26. RB1-proficient MYCNA retinoblastoma could potentially arise from a different cell of origin or at a different maturation stage than RB1−/− retinoblastomas.

Here we aimed to define the impact of MYCN activity in the context of retinoblastoma subtypes 1 and 2 in a cohort of 61 retinoblastoma samples enriched by six RB1-proficient MYCNA and two extraocular metastatic relapse samples. Genome-wide DNA methylation profiling was used to molecularly classify the 61 primary retinoblastomas. We sought circuitry defining the molecular phenotypes through bioinformatics analyses on RNA sequencing datasets from 52 samples. Our aim was to precisely define the molecular patterns of retinoblastoma with oncogenic MYCN activity within subtype 2. We generated MYCN-knockdown retinoblastoma cell models to assess functional effects and define a MYCN-RB signature, as a potential application to identify patients with retinoblastomas driven by oncogenic MYCN activity who could benefit from MYCN-directed treatment.

Results

A distinct hypomethylated DNA pattern characterizes MYCN-driven retinoblastomas

We examined epigenetic and genomic levels for a detailed molecular view of subtypes 1 and 2 in a cohort of 61 primary retinoblastoma samples enriched by 6 RB1-proficient MYCNA retinoblastoma samples (4 previously described25) and 2 samples of relapsed, extraocular retinoblastomas (Supplementary Data 1, 2). RB1 status was derived from routine diagnostic sequencing and MYCN amplifications were derived from global DNA methylation profiles using whole-genome Illumina 450k or EPIC (850k) arrays (61 samples). We applied an unsupervised DNA methylation-based cluster assembly and increased the robustness of this clustering using a beta-value discretization approach27 in a consensus clustering of 8857 unsupervised clusterings (k-means and agglomerative or hierarchical), affinity propagation, mean-shift, DBSCAN, and spectral clusterings, each with various sets of hyperparameters (Fig. 1a). In contrast to the previously described 2 retinoblastoma subtypes, our consensus clustering distinguished 3 retinoblastoma groups: cluster A (n = 35, included 2 datasets from 1 sample), cluster B (n = 17) and cluster C (n = 10, Fig. 1b). The Liu et al. 9-CpG classifier12 designed for 450k methylation array data (and based on 8 probes also used in the EPIC array) showed a correspondence between cluster A and subtype 1 retinoblastoma, and between clusters B and C and subtype 2 retinoblastoma (Fig. 1c, Supplementary Fig. S1a). In line with this, the CpG subsets defined by Li et al.28 for each cluster corresponded with our clusters (Supplementary Fig. S1b).Fig. 1 MYCN-driven retinoblastomas have distinct molecular (epi)genetic and clinical features.

a Flowchart displaying method of DNA methylation data clustering. The flowchart summarizes the steps performed to build a threshold graph reflecting consensus clustering of the DNA methylation data from 61 retinoblastomas. b Separation of 3 retinoblastoma clusters by global DNA methylation-based consensus clustering (62 datasets of retinoblastomas), cluster A (n = 35, green), cluster B (n = 17, turquoise) and cluster C (n = 10, magenta). MYCN-driven cluster retinoblastoma included RB1−/−MYCNA (pink, n = 2), RB1-proficient MYCNA (purple, n = 6), and RB1−/−non-MYCNA (maroon, n = 2). c How molecular tumor and clinical characteristics in the 61 patients with retinoblastoma differed in between the 3 clusters and correlation of clustering with the 2 subtypes grouped by the 8-CpG classifier12 (Supplementary Data 1). d The numbers of differentially methylated CpG sites (Welsh t-test; abs.diff 0.2; Benjamini–Hochberg adjusted p < 0.001) in and outside of CpG islands (CpGi) are listed. e Differentially methylated CpGs in cluster C compared to cluster A|B and in cluster B compared to cluster A|C are depicted in the 2 heatmaps. f Density plots for DNA methylation levels (ß-values) of differentially methylated probes in cluster C vs. cluster A|B (left) and cluster B vs. cluster A|C (right). Inlet plots represent median density plots per cluster.

More CpGs were differentially methylated between clusters C and A (33462 probes) than between clusters B and A (2776 probes, Fig. 1d), emphasizing a distinct methylation pattern in cluster C retinoblastomas. Hypomethylation was characteristic for differentially methylated CpGs in cluster C (Fig. 1e), while most differentially methylated CpGs in cluster B retinoblastomas were hypermethylated (Fig. 1e, Supplementary Fig. S1c). Cluster C-specific hypomethylation and cluster B-specific hypermethylation occurred both inside and outside CpG islands (Fig. 1f). Only 8.4% of CpGs differentially methylated in cluster C were within islands, while 33.5% of differentially methylated CpGs in cluster B were in CpG islands. All MYCNA retinoblastomas (6 RB1-proficient MYCNA, 2 RB1−/−MYCNA) clustered together with 2 RB1−/−retinoblastomas lacking MYCN amplifications as cluster C (Fig. 1c). The 2 retinoblastomas lacking MYCN amplifications had aberrations that could trigger oncogenic MYCN activity through other routes. The focal 13q31.3 amplification in 1 retinoblastoma harbored MIR17HG, a microRNA known to activate MYC(N) signaling, and part of GPC5 (Supplementary Fig. S2a, c, Supplementary Data 3). The other retinoblastoma had multiple genomic abnormalities on chromosomes 6, 7, and 8, in line with chromothripsis, chromoanasynthesis or chromoplexy including an amplification containing the downstream MYC(N) signaling component, BRAF (Supplementary Fig. S2b, c, Supplementary Data 3).

Children with cluster A retinoblastomas were significantly younger at diagnosis (median: 0.87 years) than patients with cluster B retinoblastomas (median: 2.36 years, p-value (p)Wilcoxon rank test = 5.785e−08; Fig. 1c, Supplementary Fig. S2d, Supplementary Data 1). The age at which children were diagnosed with RB1−/− retinoblastomas in cluster C (median: 2.86 years) did not statistically differ from children with cluster B retinoblastomas (median: 2.36 years, pWilcoxon rank Test = 0.83; Supplementary Fig. S2d), while the 6 patients with RB1-proficient retinoblastomas in cluster C were significantly younger (median: 0.38 years, vs cluster B pWilcoxon rank test = 0.006961; Supplementary Fig. S2d). Classification separated the 2 extraocular relapsed retinoblastomas, grouping 1 in cluster A (subtype 1) and the other in cluster B (subtype 2). Gain of 6p was observed in all 3 clusters, while 1q gain and 16q loss were observed predominantly in cluster B, and 2p gain was exclusive to cluster C (Fig. 1c, Supplementary Fig. S2e). Genetic mutations (whole-exome sequencing) other than RB1 were rare, apart from somatic mutations in BCOR and the ARID family genes, ARID1A and ARID4A, as previously reported7,12 (Supplementary Data 4, 5). Consensus clustering subdivided subtype 2 retinoblastomas into 2 clusters, B and C, with distinct patterns in MYCN amplification, genomic rearrangements, and methylation profiles.

NKX2-5 and MYCN are the transcription factors best demarcating clusters B and C

We applied RNA sequencing to the 52 samples from the 61-retinoblastoma cohort (Supplementary Data 2). To characterize the activated signaling pathways in each cluster, we compared the RNA sequencing datasets from the 52 primary retinoblastomas cluster-wise (31 cluster A, 15 cluster B, 6 cluster C retinoblastomas) using sleuth29. Comparing gene expression from retinoblastomas in cluster A (subtype 1) with clusters B|C (subtype 2) retrieved 4887 differentially expressed transcripts corresponding to 3109 protein-coding genes (2403 upregulated genes in cluster A, 706 upregulated genes in clusters B|C, |b|> 0.3; Fig. 2a, Supplementary Data 6). Gene-set enrichment analysis (GSEA) showed that more strongly expressed genes in cluster A retinoblastomas were involved in inflammatory and interferon-gamma responses (Supplementary Fig. S3a, Supplementary Data 7). Clusters B and C shared high expression of genes from signatures defined for photoreceptor-poor, undifferentiated retinoblastoma (Liu score12: pAvsB = 3.1e−09, pAvsC = 2.2e−05, pCvsB = 0.78; Kooi Score13: pAvsB = 1.3e−08, pAvsC = 4.6e−05, pCvsB = 0.78; Wilcoxon rank-sum test; Fig. 2b). In line with this, pairwise comparisons of clusters B and C with cluster A (Supplementary Data 8, 9) showed that the majority of upregulated genes in cluster B (73%) and cluster C (82%) belonged to genes previously shown to be representative of photoreceptor-poor, undifferentiated retinoblastoma12,13 (Supplementary Fig. S3b). Clusters B and C shared 41 upregulated genes and 113 downregulated genes that included TFF1, EBF3, BCL11A, and SOX11 (Supplementary Data 10, Supplementary Fig. S3c). High-level MYCN expression characterized cluster C, while strong expression of the NKX2-5 and GATA4 transcription factors characterized cluster B (Fig. 2c, Supplementary Fig. S3d, Supplementary Data 8, 9). However, transcription factors co-expressed with MYCN or NKX2-5 (Supplementary Data 11) were not mutually exclusive. Only MYC expression showed a trend to anti-correlate with MYCN itself and factors co-expressed with MYCN (Fig. 2d).Fig. 2 Transcriptomic landscape in clusters B and C retinoblastomas (including RB1-proficient retinoblastoma).

a Volcano plot of differentially expressed transcripts between cluster A and cluster B|C. Up- and downregulated transcripts are plotted in violet and green, respectively. Key transcripts are plotted with the corresponding annotated gene symbols (Supplementary Data 6). b Box plots showing distribution of expression scores from genes representative of photoreceptor-poor, undifferentiated retinoblastoma (Liu score12: 3105 genes positively expressed; 3088 genes downregulated, adjusted p ≤ 0.05; or Kooi-score13: 3425 genes upregulated; 3477 genes downregulated, adjusted p ≤ 0.05) in 52 retinoblastomas. c Box plots showing aggregated log1p-transformed transcript tpm (transcripts per kilobase million) values and combined Storey-Tibshirani-adjusted p for MYCN and NKX2-5 expression (left) and associated gene sets (right) in 52 retinoblastomas. d Expression of transcription factors correlating with NKX2-5 and MYCN is illustrated in the circle plot. Color and size of the circles represent the correlation coefficients. e Schemes showing expression of genes defined for cluster B and C gene signatures. The top 15 upregulated and downregulated genes within these signatures are listed with the corresponding cytobands (complete outputs, Supplementary Data 15, 16). Genes encoding transcription factors are in bold. f Scatter plot demonstrating the differentially expression of genes from cluster B (524 transcripts, x-axis) and cluster C (117 transcripts, y-axis) signatures in 52 retinoblastomas. The average of log1p-transformed tpm values for a gene set was used as a signature score of each retinoblastoma sample. Boxes show upper and lower quantile and median. Whiskers extend from the hinge to ±1.5 times the interquartile range.

High MYC(N) target gene expression characterized cluster C, as well as heightened gene activity for ribosomal function and protein synthesis and low expression of genes involved in cilia function (GSEA, Supplementary Data 12). Cluster B was characterized by low immune-related gene activity and high target expression for a component of the RB-like, E2F and multi-vulval class B (DREAM) complex, hydroxysteroid dehydrogenase and polycomb repressive complex 2 (PRC2, Supplementary Data 13). Compared to cluster C, a subset of immune response-specific signatures was depleted and genes involved in transcriptional regulation were overrepresented in cluster B (Supplementary Fig. S3d, Supplementary Data 14). These data implicate distinct pathways in clusters C and B, with both routes leading to retinoblastoma de-differentiation.

We selected differentially expressed gene signatures defining either cluster C or B to better understand differences in cluster expression circuitry. A step-wise statistical filtering approach was applied to genes differentially expressed among clusters A through C to obtain signatures specific to cluster B or C. The signature defining cluster C contained fewer differentially expressed transcripts (cluster C: 372, cluster B: 1886, Supplementary Data 15, 16). Among the top 15 genes strongly expressed in only cluster C were several encoding bHLH transcriptions factors, including MYCN, NHLH1, and NEUROG1, while the developmentally regulated homeobox transcription factors, NKX2-5, ONECUT1, and GSC were among the genes defining cluster B (Fig. 2e). No other transcription factor was more abundant (transcript level) than MYCN in MYCNA retinoblastomas. The signature defining cluster C is expressed at a low level in cluster B (Fig. 2f, Supplementary Fig. S3e), and vice versa. GSEA recapitulated the expression signatures in pairwise comparisons of clusters B and C against A. Notably, retinoblastomas in both clusters B and C strongly expressed genes from a stemness signature and genes associated with neuroblasts and neuronal progenitors. Retinoblastomas in clusters B and C appear to share transcriptional characteristics, but be driven by distinct molecular circuitry indicated by their specific transcriptional hallmarks.

CpGs differentially methylated in clusters B and C are localized to genes defining these clusters

We detected hypermethylation in cluster B and hypomethylation in cluster C retinoblastomas (Fig. 1e, f). Differentially methylated genes in cluster B were 97% hypermethylated (compared to A│C) and in cluster C were 99% hypomethylated (compared to A│B). The presence of specific transcriptional hallmarks of clusters B and C suggested a correlative connection between our DNA methylation-based clustering and defining molecular circuitry in the retinoblastomas. To explore processes that DNA methylation could regulate, we analyzed enrichment of gene ontology (GO) terms corresponding to differentially methylated CpGs in each cluster. Hypermethylated genes in cluster B were involved in camera-eye development, regulating neuronal function and synapse organization (Fig. 3a). Cluster C was enriched with hypomethylated genes involved in smell perception (Fig. 3a, Supplementary Fig. S4a). Minimal regions around hypermethylated CpGs, distinguishing cluster B, were enriched for binding motifs specific for homeodomain and zinc-finger transcription factors (HOMER analysis, Fig. 3b, Supplementary Data 17). In contrast, the E-box (bHLH-binding motifs) variants, CACCTG and CATCTG, were enriched in minimal regions around hypomethylated CpGs in cluster C (HOMER analysis, Fig. 3b, Supplementary Data 18), indicating a potentially more open chromatin configuration for MYCN-driven gene regulation. Our results suggest that bHLH factors, such as MYCN, are highly expressed in cluster C retinoblastoma.Fig. 3 Differentially methylated CpGs in retinoblastoma clusters B and C.

a Biological processes retrieved from gene ontology enrichment analysis of the genes matching to differentially methylated CpGs between clusters B and A|C (left) and differentially methylated CpGs between clusters C and A|B (right). Dot size and color represent the number of genes and enrichment significance, respectively. x-axis indicates the gene enrichment ratio (GeneRatio) of a biological process GO term. b Selected motifs identified in 50-bp sequences surrounding CpGs that are hypermethylated in clusters B vs. A|C (4109 sequences, left) or hypomethylated in clusters C vs. A|B (25073 sequences, right) identified by HOMER known DNA-motif enrichment (complete outputs, Supplementary Data 17, 18).

We integrated our RNA sequencing and DNA methylation data to explore local regulation of gene expression via autosomal CpGs (in and outside islands) annotated to the protein-coding genes in each signature (Fig. 4a). The proposed surrogate marker for subtype 2 retinoblastoma, TFF1, was hypomethylated and strongly expressed in clusters B and C (Fig. 4a, b). The bHLH transcription factor, NHLH1, was hypomethylated and strongly expressed in cluster C (Fig. 4a, b). The majority of genes with differentially methylated CpGs (6320 CpGs of 2465 genes) negatively correlated with expression (R ≤ 0.4, Fig. 5a), while 2261 CpGs of 992 genes, including NKX2-5, positively correlated with gene expression (R ≥ 0.4; Fig. 4a, b, Supplementary Data 19–23). GO enrichment analysis revealed enrichment of neuronal development and function for negatively correlated hypomethylated genes in clusters B and C (Fig. 4c). Negatively correlated hypermethylated genes in clusters B and C were enriched for cell-to-cell communication and neuronal development terms (Fig. 4c). Positively correlated genes in cluster B and cluster C were enriched for the GO terms, muscle tissue development, mesenchyme development and mesenchymal cell differentiation (Fig. 4d). Our comparisons of differentially methylated CpGs between clusters (Fig. 1d) also identified 28% of CpGs analyzed in correlation analysis. Genes associated with enhancer elements previously annotated in developing retina or retinoblastoma30,31 were hypomethylated and highly expressed in our cohort (Supplementary Data 24, Supplementary Fig. S4b). Yet, enhancer elements, with the exception of hypomethylated in cluster C enhancers associated with MYCN and NHLH1, lacked cluster-specific differences in methylation (Supplementary Data 24). Our methylation and expression data confirm the previous association of TFF1 with subtype 2 retinoblastoma, provide further support for the defining role of MYCN activity in cluster C.Fig. 4 Integrative analysis of DNA methylation and gene expression in retinoblastoma clusters B and C.

a Volcano plot for Pearson correlations between mRNAs and methylation values of corresponding local CpGs in CpG island (left) and non-CpG island (right) contexts in 52 retinoblastomas (Supplementary Data 19, 20). Y-axes represent logarithm of p-value of each correlation coefficient. X-axes represent the Pearson correlation coefficient, r. Selected highly correlated CpG-mRNA pairs are labeled in blue (negative correlation) and red (positive correlation). b Relationship between gene expression in log1p-transformed tpm values and DNA methylation of CpGs (β-values) for selected genes in 52 primary retinoblastomas. Y-axes represent mRNA expression. X-axes represent DNA methylation. Each dot corresponds to a retinoblastoma sample and is color-coded by retinoblastoma cluster (green, cluster A; turquoise, cluster B; magenta, cluster C). Pearson (R) and Spearman (rho) coefficients with corresponding p-values, linear regression line and confidence interval are indicated for each comparison. Box plots showing average of log1p-transformed tpm values of negatively (c) and positively correlated (d) genes in retinoblastoma clusters (upper panels) and dot plots showing biological processes associated with these genes (lower panels). Boxes show upper and lower quantile and median. Whiskers extend from the hinge to ±1.5 times the interquartile range or the highest/lowest value.

Fig. 5 MYCN knockdown reduces clonogenicity and tumor formation in MYCN-knockdown retinoblastoma cell models.

a MYCN mRNA expression in retinoblastoma cell lines determined by RNAseq. b Flow cytometry analysis of MYCN protein expression in RB355, RB522, RB3823, and WERI-Rb1 cell lines transduced with doxycycline-inducible lentiviral shMYCN expression vector. MYCN protein was measured 48 h after adding doxycycline to the media. The results represent 3 experiments. c Soft agar colony formation analysis in retinoblastoma cell lines. Plots (upper panel) show ratios of colony numbers in shNTC-expressing (doxycycline-treated cells vs non-treated cells) and shMYCN-expressing (doxycycline-treated cells vs non-treated cells) retinoblastoma cell lines from three independent experiments. Data are from 3 independent experiments, with average values (blue bots) indicated and p-values calculated with Welch t-test reported. Density plots (lower panel) for colony size in shNTC-expressing and shMYCN-expressing retinoblastoma cell lines. Each line corresponds to one technical replicate. Representative images from shMYCN-expressing non-treated and doxycycline-treated cells are shown for each graph (shNTC-expressing non-treated and doxycycline-treated cells are in Supplementary Fig. S6d). Combined p-values calculated with one-sided Wilcoxon rank-sum test for the arbitrary difference in clone size (10% of maximum clone size) are shown. d Effects of MYCN knockdown in RB355 cells and the chick chorioallantoic model. Photos and plots, respectively, show the tumors in eggs at developmental stage E17 and the fraction of tumor cells and tumor area in engrafted eggs at stage E17. The results represent 3 replicates with 10 eggs for each condition.

MYCN knockdown rescues molecular circuitry driving subtype I in retinoblastoma cell models

We profiled 9 previously described retinoblastoma cell lines32 (Supplementary Fig. S5a, Supplementary Data 2) with our multi-omics approach. The RB1-proficient cell lines, RB522 and RB3823, and the RB1−/−MYCNA cell line, Y79, expressed the highest MYCN mRNA and MYCN protein levels and showed the highest proliferation rates (Fig. 5a, Supplementary Fig. S5b). Applying the Liu et al. 9-CpG classifier12 to these data grouped all 9 cell lines with subtype 2 retinoblastoma (Supplementary Fig. S5c). The 4 MYCNA cell lines comprised a separate branch closer to cluster C than B based on hierarchical clustering of mRNA expression of transcription factors expressed in these 2 retinoblastoma clusters (Supplementary Fig. S5d, e). The highest expression scores for cluster C signature genes were detected in the fast-growing, MYCNA RB3832, RB522 and Y79 cell lines (Supplementary Fig. S5f). Since no cell lines are derived from subtype 1 retinoblastoma, we explored the functional effects of MYCN in the subtype 2 background by generating MYCN knockdown models from these cell lines.

To define molecular and functional changes driven by MYCN in subtype 2 retinoblastomas, inducible knockdown models were generated from 4 retinoblastoma cells lines. We selected the 2 RB1 proficient MYCNA cell lines (shRNA-based MYCN knockdown in RB522, RB3832) because of high MYCN RNA and protein expression and two RB1−/− cell lines with and without MYCN amplification (shRNA-based MYCN knockdown in RB355 MYCNA cell lines and the MYCN-diploid WERI-Rb1 cell line). Knockdown reduced MYCN levels by 50–65% in cell models (Fig. 5b, Supplementary Fig. S6a). The proportion of cycling cells and the ability to form colonies was reduced in culture (Supplementary Fig. S6b, c, Fig. 5c, Supplementary Fig. S6d, e) and tumor growth was reduced in the in vivo chorioallantoic membrane model (Fig. 5d). Analysis of RNA sequencing data across all 4 cell models identified a common set of 477 genes that were differentially expressed (93 upregulated, 385 downregulated, Storey-Tibshirani p ≤ 0.1; Supplementary Data 25). More than 50% of commonly downregulated genes (n = 244) were known MYC(N) targets23,33–35 and define a MYCN-RB signature (Fig. 6a, Supplementary Data 26). MYCN knockdown downregulated expression of genes involved in protein synthesis and cell metabolism (Supplementary Data 27, second most common functional group).Fig. 6 Gene re-expression representative of photoreceptor-rich retinoblastoma in MYCN-knockdown cell models.

a Volcano plot shows the MYCN-RB signature of 244 genes derived from the overlap of previously identified MYCN targets and sleuth-modeled MYCN-knockdown expression pattern (Supplementary Data 25). Red dots correspond to significantly regulated MYCN targets (Storey-Tibshirani adjusted p ≤ 0.1). Orange dots correspond to other MYCN targets (Storey-Tibshirani adjusted p > 0.1). Black dots represent non-MYCN targets significantly downregulated in sleuth-modeled MYCN-knockdown expression pattern. b GSEA showing genes overexpressed in cluster A, genes representative of photoreceptor-rich retinoblastoma12,13 and the cluster B signature in the sleuth-modeled MYCN-knockdown expression pattern. Y-axis indicates enrichment score (ES). X-axis shows pathway genes. The dual-colored band represents the degree of correlation in the expression of these genes with MYCN-knockdown (red, MYCN-knockdown; blue, control). c Box plot (left) showing aggregated log1p-transformed transcript tpm values for OPTN expression and combined Storey-Tibshirani-adjusted p in 52 retinoblastomas. OPTN protein expression (right) in MYCN-knockdown retinoblastoma cell lines 48 h after doxycycline treatment. Dot plots show ratios of median fluorescence intensity (MFI) for OPTN in shNTC-expressing (doxycycline-treated cells vs non-treated cells) and shMYCN-expressing (doxycycline-treated cells vs non-treated cells) retinoblastoma from three independent experiments. Data are from 3 independent experiments, with average values (blue bots) indicated and p-values calculated with Welch t-test reported. Representative plots of flow cytometry for OPTN are shown in Supplementary Fig. S7c. Venn diagrams showing the pathways (FDR adjusted p ≤ 0.05) associated with downregulated genes (d) and upregulated genes (e) in MYCN-knockdown retinoblastoma cell lines. The 9 pathways depleted in all 4 cell lines are listed. f List of genes upregulated in both RB1−/− cell lines. Box plots showing the expression scores of the MYCN-RB signatures (g), and upregulated genes (h) after MYCN-knockdown in each of the 3 retinoblastoma clusters. Boxes show upper and lower quantile and median. Whiskers extend from the hinge to ±1.5 times the interquartile range.

Genes known to be markers for photoreceptor cells (NRL, OPTN, SLC17A7, PLEKHB1, GNGT2, ATP2B1, and BAZ2B) were among those upregulated in models after MYCN knockdown (Supplementary Data 25). MYCN knockdown upregulated genes descriptive of the photoreceptor-rich cluster A retinoblastoma (Fig. 6b) and downregulated genes representative of the photoreceptor-poor, undifferentiated and cluster C retinoblastoma (Supplementary Fig. S7a, b) in cell models. None of the genes defining cluster B were regulated by MYCN knockdown (sleuth-modeled MYCN-knockdown expression profile, Fig. 6b; pairwise non-induced/induced cell model comparisons, Supplementary Data 28), indicating MYCN plays no role in the molecular circuitry driving cluster B retinoblastoma. We confirmed upregulation of cluster A expressed OPTN on protein level using flow cytometry analysis (Fig. 6c, Supplementary Fig. S7c). Notably, pathways associated with genes downregulated by MYCN knockdown were consistent across all 4 cell models (Fig. 6d). The pathways associated with upregulated genes were similar within each of the 2 RB1−/− cell models (RB355 and WERI-Rb1) and RB1-proficient cell models (RB522 and RB3832), but differed between the 2 groups (Fig. 6e). In the absence of functional RB1, MYCN knockdown upregulated more genes encoding photoreceptor or ciliary proteins (compared with RB1-proficient cell lines; Fig. 6f, Supplementary Data 29), suggesting a stronger redifferentiation of the cells.

We applied the MYCN-RB signature to our RNA sequencing dataset from 52 primary retinoblastomas. MYCN-RB expression scores were lowest in cluster A and highest in cluster C (pAC = 0.00058, pAB = 0.00016, pBC = 0.02; Fig. 6g), peaking in MYCNA retinoblastomas (pMYCNnon-A-MYCNAWilcoxon rank sum test = 3.694e−06) and correlating with MYCN expression (R = 0.69, RMYCNnon-A = 0.60, p < 0.00001). We also assessed the opposite picture of expression (genes upregulated by MYCN knockdown) from our sleuth-modeled MYCN-knockdown expression profile in the 52 primary retinoblastomas. Retinoblastomas lacking MYCN amplifications strongly expressed this gene group (pMYCNnon-A-MYCNAWilcoxon rank sum test = 0.0002623), and expression was highest in cluster A, followed by cluster B (pAC = 0.00058, pBC = 0.00227, pAB = 0.37503; Wilcoxon rank-sum test; Fig. 6h). In our cohort of primary retinoblastoma, the MYCN-RB signature generated in the cell line model was able to identify MYCN-driven retinoblastomas with and without MYCNA, which cluster separately based on methylation as cluster C retinoblastomas. To validate this finding, we reanalyzed the DNA methylation array data of the Liu et al. cohort12 using our subsets of cluster B- and cluster C-specific CpGs, which revealed two cluster C retinoblastomas (Supplementary Fig. S8a, b). Based on Liu et al. RNA expression data, these two cluster C retinoblastoma and one retinoblastoma without methylation data showed the highest MYCN-RB signature scores (Supplementary Fig. S8c). These three tumors were identified as MYCNA, supporting the validity of MYCN-RB signature to identify MYCN-amplified retinoblastomas.

Likewise, the expression of genes descriptive of the photoreceptor-rich retinoblastoma in the cell line models after MYCN knockdown shows that MYCN is involved in maintaining undifferentiated state of retinoblastoma cells.

Discussion

Here, we first define a MYCN-driven group of subtype 2 retinoblastoma characterized by hypomethylation and high expression of MYCN target genes. This subtype 2-MYCN (cluster C) comprises not only the rare RB1-proficient MYCNA retinoblastoma but also RB1−/− retinoblastomas with overactive components of MYCN signaling. Knockdown of MYCN in cell models reactivated the differentiated subtype 1 photoreceptor signature, suggesting a transition between the differentiated subtype 1 retinoblastoma and the more aggressive subtype 2 retinoblastoma that can potentially be reversed via targeting MYCN. We also defined and validated a MYCN-RB gene signature to identify retinoblastomas that could benefit from MYCN-targeting treatment.

In contrast to the previously described hypermethylation pattern12,16 of subtype 2 retinoblastoma, this subtype 2-MYCN (cluster C) was characterized by hypomethylation and showed distinct patterns of genomic rearrangements, gene expression profiles and clinical characteristics to cluster B, underlining biological differences between these two clusters within subtype 2 (Fig. 7). Our data identify cluster C as a subgroup with different defining molecular circuitry and hallmarks in subtype 2 retinoblastomas12. Yet, the absence of anti-correlated gene expression between cluster B- and cluster C-specific transcription factors indicate that expression programs in cluster B and C retinoblastomas are probably not mutually exclusive. However, correlation of gene expression and methylation clearly confirmed MYCN as the transcription factor best demarcating cluster C from B.Fig. 7 Overview of characteristics of 3 retinoblastoma clusters (including RB1-proficient retinoblastoma).

The MYCN-driven cluster C included not only MYCN-amplified retinoblastomas, but two samples with high MYCN target expression that lacked MYCN amplifications, suggesting alternative ways of activating MYCN-driven tumor progression in retinoblastoma. Gene amplification is one of several mechanisms that lead to MYCN deregulation and activation of MYCN targets36. Several microRNAs, including MIR17HG, are involved in the complex regulation of MYC(N) signaling37, so that it is plausible that MIR17HG amplification drives MYCN signaling in the selected cluster C retinoblastoma. The high expression of MYC and other bHLH transcription factors (NHLH1 and NEUROG1) was similar among cluster C retinoblastomas, regardless of MYCN amplification status. The enrichment of CACCTG and CATCTG E-box variants, which are the preferred by TWIST, NEURO, and TCF/LEF family transcription factors and can bind MYCN under deregulated conditions38, emphasize the importance of bHLH transcription factors for cluster C retinoblastomas. Remarkably, hypomethylation of several olfactory receptor genes was characteristic for subtype 2-MYCN (cluster C). The olfactory receptor gene family is the largest in the human genome39. Ectopically expressed olfactory receptor genes have been linked with differentiation of cancer cells and prognosis in various cancers40. Activation of olfactory receptor genes may reflect de-differentiation of subtype 2-MYCN (cluster C) retinoblastoma. As MYCN amplification alone was insufficient to identify MYCN-driven retinoblastomas, we defined the MYCN-RB signature using data from MYCN-knockdown cell models. The MYCN-RB signature identified cluster C retinoblastomas and can be used to identify retinoblastoma cases that may benefit from MYCN-directed therapies.

MYCN is a known oncogene in various embryonal tumors and MYCN amplification is a well-known prognostic factor for patients with neuroblastoma41. A previous study demonstrated that MYCN oncogenic programs differed between retinoblastoma and neuroblastoma42. In contrast to neuroblastoma43, the MYCN-associated methylation pattern was characterized by hypomethylation in our retinoblastoma cohort, but results from our MYCN-knockdown cell models support a significant overlap between MYCN targets in retinoblastoma and known MYCN and MYC targets. Genes involved in protein synthesis appeared to be controlled by MYCN in retinoblastoma, consistent with previous findings in other cancers44,45. This supports the idea of MYCN cooperation with translation processes in retinoblastoma and opens up perspectives for new drug targets in MYCN-driven retinoblastomas. In retinoblastoma, MYCN amplification occurs in both the more common RB1−/− and rare RB1-proficient backgrounds25,46. RB1-proficient retinoblastoma has been described as a separate retinoblastoma entity with very aggressive behavior, in which MYCN is thought to be the exclusive driver of tumorigenesis25,46. Our consensus clustering based on genome-wide methylation did not distinguish between RB1-proficient and RB1−/− MYCN-amplified retinoblastomas, as other groups have reported16. While this is a relatively large group of the exceedingly rare RB1-proficient retinoblastomas, the comparison is limited by small cohort size, such that batch effects cannot completely be excluded. Apart from the previously reported very young age of diagnosis for RB1-proficient retinoblastoma, we did not observe differences in gene expression, chromosomal or genetic features between the RB1−/− and RB1-proficient MYCN-amplified retinoblastomas. Only in our cell models did differences become apparent after MYCN knockdown, when RB1−/− cell lines showed more pronounced expression of the photoreceptor signature than RB1-proficient retinoblastomas. DNA methylation and gene expression profiling in further RB1-proficient MYCN-amplified retinoblastomas will help to better characterize differences in molecular circuitry in RB1−/− and RB1-proficient MYCN-amplified retinoblastoma.

MYCN knockdown in our cell models caused a partial growth arrest and reactivated the photoreceptor gene signature of subtype 1 retinoblastoma in the subtype 2 cell lines. These observations can be explained by differentiation of retinoblastoma cells along a cone photoreceptor cell trajectory, and can be interpreted as a dynamic connection between subtype 2 and subtype 1 retinoblastoma. MYCN blocks differentiation pathways and maintains pluripotency during development and cancer pathogenesis47–49. Enforced reduction of MYCN expression in neuroblastoma cells is also associated with cell differentiation50. Kapatai et al.14 has suggested subtype 2 retinoblastoma originates from an earlier retinal lineage or an early uncommitted cell type, but subtype 2 retinoblastoma could alternatively originate from subtype 1 retinoblastoma via acquiring additional genetic and chromosomal alterations3,13. All available retinoblastoma cell lines grouped as subtype 2 retinoblastoma. This could either reflect that the more differentiated subtype 1 retinoblastoma is not amenable to culturing or that subtype 1 retinoblastoma acquire additional genetic and epigenetic alterations in culture and advance to subtype 2 retinoblastoma. The differentiation observed upon MYCN knockdown could hint that subtype 1 retinoblastoma is able to transition to subtype 2-MYCN (cluster C) retinoblastoma upon MYCN activation.

The lack of serial tissue sampling during retinoblastoma treatment leaves the question whether subtype 1 transitions into subtype 2 retinoblastoma by acquiring additional mutations3,13 or whether both subtypes arise independently from the same progenitor cell at different stages of maturation unresolved12,14. Subtype 2 occurs in older children, shows a more undifferentiated gene expression profile and predominates in metastatic cases12. One of the two extraocular retinoblastomas that we analyzed was unexpectedly a cluster A, subtype 1 retinoblastoma without elevated TFF1 expression. Remarkably, this patient is a long-term survivor despite localized central nervous system metastasis.

Few extraocular cases have been subtyped because molecular genetic analysis is conducted almost exclusively on cases treated in high-income countries, where nearly all retinoblastomas are diagnosed and treated during early intraocular stages. Accordingly, most primary retinoblastomas in our cohort were subtype 1 retinoblastomas (35/61, 57%). Countries with limited resources have high mortality rates from retinoblastoma5 because cases most often progress to extraocular metastasis. This implies that either geographic factors affect molecular retinoblastoma subtype distribution or that subtype 1 retinoblastoma can advance to subtype 2 retinoblastoma or take an aggressive clinical course given time. Collaborative international research projects are needed to answer this question and better define molecular retinoblastoma subtypes.

Given its correlation with rapid tumor progression and poor prognosis in several embryonal cancers MYCN is considered an ideal therapeutic target, but many direct or indirect MYCN modulators failed to result in an efficient MYCN specific therapy51,52. New approaches including histone deacetylases inhibitors, Aurora Kinase-A Inhibitors and Bromodomain and extra-terminal domain family Inhibitors show promising interaction with MYCN and its pathways51. As biopsies are obsolete in retinoblastoma, all tissue samples in our study derived from patients with advanced intraocular disease treated with enucleation. Therefore, molecular genetic characteristics of small intraocular retinoblastoma remain undetermined. This potential bias may be overcome in the future by the use of liquid biopsies. Results of molecular genetic characterization of cfDNA derived from aqueous humor at diagnosis, during treatmentis promising to provide information on smaller retinoblastoma, serial samples and correlation with response to therapy53. Analysis of cfRNA from aqueous humor could in the future also allow to identify MYCN-RB signature in retinoblastoma treated with eye-preserving therapy.

Our data suggest that de-differentiation from subtype 1 to subtype 2-MYCN (cluster C) retinoblastoma can be reversed by inhibiting MYCN. The newly defined MYCN-RB signature supports identification of patients who may benefit from MYCN-directed therapies, and presents a biomarker-guided option for patient selection in clinical trials incorporating agents targeting MYCN. Our data may also be applicable to MYCN signaling pathways in other cancers in which RB1 inactivation occurs, including lung, ovarian and breast cancers.

Methods

Patients and samples

We selected a cohort of 61 intraocular primary retinoblastoma samples from patients diagnosed and treated between 2007 and 2014 at the University Hospital Essen (written informed consent for use in research available). We enriched our cohort for the rare retinoblastoma types, RB1-proficient MYCNA (DNA and RNA from 4 retinoblastomas that were first described25 was reanalyzed) and extraocular metastatic relapses (patient and tumor characteristics in Supplementary Data 1). The remaining 2 RB1-proficient MYCNA samples lacked detected RB1 variants in routine diagnostics, raising the suspicion this case is RB1-proficient MYCNA. The data from the 2 patients with relapsed extraocular retinoblastoma evaluated in this study was produced and kindly provided by the INFORM program54–57. Selection criteria for all 59 primary retinoblastoma samples were: (1) consent for research projects available, (2) adequate DNA purity and quality, (3) no treatment prior to enucleation and (4) availability of sufficient tumor and blood DNA. Samples were obtained after primary enucleation or resection of the metastasis, directly frozen and stored at −80 °C until DNA and RNA preparation. RB1 genotyping was performed within genetic testing on matched blood and tumor samples at the request of individuals or their legal guardians with the aim to identify oncogenic alterations in RB1.

Cell lines

WERI-Rb1 (RRID:CVCL_1792) and Y79 (RRID:CVCL_1893) cell lines were purchased from the German Collection of Microorganisms and Cell Cultures (DSMZ, Braunschweig, Germany). RB1021 (RRID:CVCL_S624), RB247C (RRID:CVCL_2704), RB355 (RRID:CVCL_S611), RB383 (RRID:CVCL_S626), RB3823 (RRID:CVCL_ZF07), RB522 (RRID:CVCL_ZF04) were kindly provided by Brenda Gallie (Department of Ophthalmology and Vision Sciences, Hospital for Sick Children, Toronto). Rbl30 (RRID:CVCL_S621) was kindly provided by Ralf Küppers (University Hospital Essen). Cell line identity was validated in-house (Institute for Human Genetics) using RB1 mutational status and short tandem repeats. All retinoblastoma cell lines are well established, and were previously molecularly, cellularly and functionally characterized32. Cell lines were cultivated in Dulbecco’s modified Eagle medium (ThermoFisher Scientific) supplemented with 15% fetal calf serum (Thermo Fisher Scientific), 11 µM β-mercaptoethanol (MilliporeSigma), and 100 U/ml penicillin/streptomycin (Thermo Fisher Scientific) in a humidified atmosphere with 5% CO2 at 37 °C. All cell lines except RB355 (adherent monolayer) were cultivated in suspension. WERI-RB1, RB247C3, and RB1021 are RB1−/− cell lines lacking MYCN amplifications. Y79, Rb383 and RB355 are RB1−/−MYCNA cell lines. RB3823 and RB522 are RB1-proficient MYCNA cell lines with functional RB1.

Generating doxycycline-inducible MYCN-knockdown cell models

MYCN knockdown was achieved using the Dharmacon doxycycline-inducible lentiviral SMARTvector system with PGK promoter, TurboGFP and 3 human MYCN knockdown target sequences (CGAGCTGATCCTCAAACGA targeting 3′-UTR and open reading frame, ACGTCCGCTCAAGAGTGTC targeting 3′-UTR and open reading frame and CCACATAAGGGGTTTGCCA targeting only the 3′-UTR). The SMARTvector Inducible lentiviral NTC PGK-TurboGFP served as non-targeting control. The MYCN-knockdown vector targeting only the 3′-UTR performed best in all cell lines, and was used to generate all models. The polyethylenimine (PEI) method was used to create and transduce Vesicular Stomatitis Virus Glycoprotein-G (VSV-G)-pseudotyped replication-deficient lentiviral particles into HEK293T cells as described previously58. Briefly, 6 µg each of envelope plasmid (pczVSV-G), helper plasmid (pCD/NL-BH) and pSMART plasmids were mixed with 45 µg PEI and incubated with HEK293T cells. After 24 h, medium was changed to Iscove’s modefied Dulbecco’s medium (MilliporeSigma) supplemented with 10% FCS, 1% penicillin/streptomycin. Lentiviral supernatant was harvested at 24 h, filtered (0.45 µm filters) and added (1:2 dilution) to retinoblastoma cell cultures with protamine phosphate (final concentration 5 µg/mL medium, Millipore Sigma). Medium was aspirated off cells after 24 h, and cells washed once with phosphate-buffered saline (PBS, ThermoFisher Scientific) before culturing in complete medium for 7 days. Transduced cells were selected with 1 µg/mL puromycin (ThermoFisher Scientific). Cells were cultured in medium containing 0.5 µg/mL doxycycline (STEMCELL Technologies) for 24–72 h to establish a MYCN-knockdown population prior to further analysis.

Protein staining for flow cytometry

MYCN protein staining was performed using AlexaFluor 647-labeled antibodies against MYCN (sc-53993, Santa Cruz Biotechnology) and mouse IgG2a (sc-24637, Santa Cruz Biotechnology). OPTN protein staining was performed using non-conjugated antibodies (#70928, Cell Signaling) and secondary Alexa Fluor 647-labeled antibody (A-21244, Rabbit IgG (H + L) Cross-Adsorbed Secondary Antibody, Thermo Fisher Scientific, and Rabbit mAb IgG XP®Isotype Control (#2975 Cell Signaling). Cells (106) were prefixed with 4% paraformaldehyde (Morphisto) for 15 min at room temperature, resuspended in ice-cold methanol and fixed overnight at −20 °C. For immunostaining, cells were washed in PBS, incubated for 30 min in 1% bovine serum albumin (Roth) and 0.1% Triton X-100 in PBS (MilliporeSigma) then incubated 1 h at room temperature with the antibody (0.5 mg per sample). Then washed and if required incubated for 1 h with secondary antibody (dilution 1:400).

Assaying cell proliferation and viability

The proportion of cells in S phase was determined using EdU incorporation and cell cycle analysis. Cells were incubated with 10 µM EdU (Lumiprobe) for 1 h prior to collection, then 106 cells fixed in 4% paraformaldehyde (Morphisto) at room temperature for 15 min and in 90% ice-cold methanol before storing at −20 °C. Fixed cells were permeabilized for 30 min in 1 M Tris buffer (pH 7.4) containing 1% bovine serum albumin and 0.1% Triton X-100, then stained for 30 min with 3 µM Cy5 azide (Lumiprobe) in 1 M Tris (pH 7.4) with 2 mM CuSO4 and 20 mg/ml ascorbic acid. Cells were subsequently washed with PBS and counterstained with DAPI (Millipore, Sigma) for flow cytometric cell cycle analysis. Flow cytometry was conducted on a Beckman Coulter CytoFLEX instrument. Data analysis used FlowJo software V10 (Becton Dickinson). Viability was assessed in retinoblastoma cell lines and models using the 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) assay (Roth) according to manufacturer’s instructions. Cells (105 cells for suspension cell lines, 5 × 104 for adherent cells) were seeded onto 96-well plates in quadruplicate for MTT assay, and if MYCN knockdown was to be induced, incubated 24 h before 0.5 µg/ml doxycycline treatment. For soft agar assays, cells were trypsinised, and 0.5 × 104 cells were resuspended in a top agar solution containing 0.35% agar, then poured onto a base layer containing 0.7% agar in 12-well plates. Each well was covered with 0.7 ml of media. If MYCN knockdown was to be induced, cells were supplemented with 0.5 µg/ml doxycycline in top agar solution and upper media. Plates were incubated under standard culture conditions for 2 weeks. Colonies were fixed and stained with 10% methanol, 10% acetic acid and 0.1% of crystal violet. For clonogenicity assays, cells (0.5 × 104) were seeded onto 12-well plates pre-coated with poly-L-lysine (Thermo Fisher Scientific). MYCN knockdown was induced by adding doxycycline to a final concentration of 0.5 µg/ml, and plates were incubated under standard culture conditions for 2 weeks. Cells were fixed in 75% methanol, 25% acetic acid and stained with 10% methanol, 10% acetic acid, and 0.5% of crystal violet. The colony counting, calculation of colony areas was performed using Ilastik59-processed segmented images and ImageJ software (Rasband, WS, ImageJ, U.S. National Institutes of Health, Bethesda, MD, USA). Proliferation was assessed in the MYCN-knockdown RB355 and RB355 cell models using the chick chorioallantoic membrane (CAM) assay. Knockdown was induced in model cells (0.5 µg/ml doxycycline) 48 h before engrafting. Induced or uninduced model cells (106/50 µl PBS) were grafted onto the chorioallantoic membrane of E10 stage chick embryos after cutting a window into the eggshell. Proliferation was assessed 7 days after engraftment, at stage E17, by enlarging the windows to photograph the entire anterior chorioallantoic membrane region and carefully extract the tumors. Tumors were weighed to assess volumetric growth, and tumor areas were assessed using ImageJ software. The results represent 3 replicates with 10 eggs per test condition.

Genetic, epigenetic, and genomic analyses

RB1 variants were identified using DNA from fresh-frozen tumor samples or DNA from blood and one or more of the following in routine diagnostics (Dept. of Human Genetics, Essen) as previously described60–63 analysis of allele loss in tumors, cytogenetic analysis, denaturing high performance liquid chromatography, exon-by-exon sequencing, multiplex ligation-dependent probe amplification, methylation-sensitive PCR, quantitative fluorescent multiplexed PCR, quantitative real-time PCR, real-time PCR and single-strand conformation polymorphism analysis. RB1 mutational mosaicism was diagnosed when the signal ratio for the mutant:normal allele in blood DNA < 50%, but detecS§. If available, DNA from an additional tissue source (i.e., buccal mucosa) was analyzed for RB1 mosaic cases. Genome-wide methylation was analyzed using the Illumina Infinium Human Methylation 450k array (21 retinoblastoma samples) or EPIC 850k BeadChip array (39 retinoblastoma samples). Sample RB_33 was analyzed on both arrays to control for batch effects (Supplementary Data 1). Whole-exome sequencing (WES) was performed in 37 retinoblastomas (Supplementary Data 1).

RNA sequencing

RNAseq analysis was performed in 52 of the 61 intraocular retinoblastomas (Supplementary Data 2). Total RNA from primary tumor samples was isolated using the RNeasy Mini Kit (Qiagen) or Monarch Total RNA Miniprep kit (New England Biolabs GmbH). RNAseq was performed on retinoblastoma parental cell lines and models. Models were preincubated with/without 0.5 µg/ml doxycycline for 48 h before analysis. Total RNA from cells was prepared according to manufacturer’s instructions using RNeasy mini kit (Qiagen) with a separate DNAse I (Qiagen) digestion step. Whole-transcriptome sequencing on parental cell lines and tumor samples was performed by CeGat (Tübingen, Germany) using the SMART-Seq Stranded kit (Takara Bio Europe) for library preparation and NovaSeq 6000, 2x100bp (Illumina) for sequencing. For cell models, RNA amount was determined by Qubit fluorometer (Life Technologies), RNA quality was assessed by TapeStation 4200 with RNA ScreenTapes (Agilent) and the Lexogen QuantSeq 3′ mRNA-Seq Library Prep Kit FWD (Lexogen) was used to produce the library. Library quality were assessed on the TapeStation 4200 (High Sensitivity ScreenTapes, Agilent) and concentrations were measured using the Qubit system. RNA sequencing of MYCN-knockdown models were sequenced using either the Illumina NextSeq 2000 System (P2 reagents, 76 cycles, single-read, dual indices) or MiSeq using MiSeq Reagent Kit v3 (76 cycles, single-read, dual indices) using 3 biological replicates for each model and condition.

DNA methylation data processing and consensus clustering

Two retinoblastomas that were previously profiled (450k) within the INFORM study54–57 to identify molecular targets for therapeutic options were reused for our analyses (Supplementary Data 1). Datasets generated on 450k and EPIC arrays were reduced to the 450 K common probes and processed using minfi (version 1.36), then normalized using minfi’s funnorm function. Batch correction was conducted for each array type using comBat (from sva version 3.42). The resulting methylation levels were discretized into 3 classes: unmethylated (levels close to 0), methylated (levels close to 1), and unknown or semi-methylated (remainder). Exact cut-off thresholds were determined by fitting a 3-component mixture model of beta distributions to the histogram of methylation levels and identifying points of equal probability density, as described previously27. Probes were then filtered to remove probes without information content, retaining only probes that showed the “methylated” or “unmethylated” state in at least 10 different samples. Consensus clustering was performed on the reduced discretized sample-probe matrix by applying different clustering algorithms from the Python scikit-learn package (K-means, affinity propagation, agglomerative or hierarchical with single, average and complete linkage and cosine, Euclidean and Manhattan distance; spectral clustering; DBSCAN; MeanShift) in combination with different preprocessing methods (robust scaling, normalization and principal component analysis or UMAP with 1 to 4 dimensions) and different output cluster numbers (2 to 6). From the resulting 8857 successful individual clusterings, those with a silhouette score of zero or less were discarded. Remaining clusterings were used to define a similarity matrix between samples, describing in how many clusterings two samples were assigned to the same cluster. The matrix was transformed into a similarity graph, such that each sample is a node and weighted edges between nodes indicate relative similarity (between 0 and 1; proportion of clusterings for which the 2 samples are in the same cluster). As a last step, edges with a weight <0.5 were removed resulting in 3 separated components defining the final 3 clusters.

Differential methylation analysis

Retinoblastoma subtype-specific clustering was carried out using the published 8-probe classifier12, then confirmed using a modified 8-probe classifier supplemented with EPIC array probe cg04786667, which is closest to 450 array probe cg12750745 that is absent from the EPIC array. Hierarchical clustering (dplyr R package) was performed for the 8-probe classifiers with the parameters (distance: Manhattan, linkage: Ward.D2). Heatmaps were produced using the gplots R package. Differential methylation in retinoblastoma clusters was defined by differentially methylated probe sets for the batch-corrected (ComBat) whole-genome methylation data from tumors (Welch t-test with Benjamini–Hochberg false discovery rate adjustment; cut-off: β value mean difference of at least 0.2, adjusted p < 0.001). Cell lines were matched to cluster methylation profiles by extending the matrix for tumor methylation with β values from cell lines. Differential methylation in specified genomic regions (super-enhancers, CpG islands, etc) utilized the Bioconductor regioneR package on the usegalaxy.org public server64. CpG island coordinates were retrieved from the UCSC server (2020 update, http://hgdownload.cse.ucsc.edu/goldenpath/hg19/database/cpgIslandExt.txt.gz;). BED files corresponding to the enhancers identified in developing retinoblastoma and retina30 (GSE86981) were extracted from the SEDb2.0 database31. The Homer script, findMotifs.pl (http://homer.ucsd.edu/homer/, v4.10.3), was used to compare 50-bp sequences (containing known motifs) around CpGs (51 bp) that were differentially methylated with a strand-specific input.

RNA sequencing data analysis

Differential gene expression was analyzed among the 52 retinoblastomas profiled by RNA sequencing. RNA sequencing reads were preprocessed and aligned to the H. sapiens normal transcriptome construct (Ensembl v96 release) and relative transcript abundance was quantified for each retinoblastoma using kallisto (-k 31 parameter, -b 100 parameter)65. To identify differential gene expression, p-values were computed using sleuth (Wald test), then used to generate false discovery rates (q values, using the qvalue package) that were adjusted for multiple comparisons using the Storey-Tibshirani method (bootstrap method to calculate π0 values)66. p- or q-values were combined using Edgington’s67 method and metap package. Box plot visualization for intertumor comparisons of differential expression data from the matrix created by kallisto required batch correction, which was performed using removeBatchEffect from limma package (Bioconductor) using log1p-transformed tpm values. Differential expression in cell lines and models was analyzed using DESeq2-based Quantseq 2.3.6 FWD pipeline and Quantseq DE 1.4.0 pipeline on the BlueBee® Genomics Platform (Lexogen). Differentially expressed genes were ranked in cell models and tumors using π-values68. Gene sets from MSigDB69 were used for GSEA with parameter settings (statistic: weighted, number of permutations: 1000, minimal gene set size: 10, maximal gene set size: 500 or 5000 for customized gene sets). Spearman coefficients for correlation were visualized with corrplot. Hierarchical clustering of RNA expression data was performed using the dplyr R package with the parameter settings (distance: Euclidean, linkage: Ward.D2). Heatmaps were produced using gplots R package. Higher-level chromosome rearrangements (translocations, duplications) were analyzed using the Arriba pipeline for RB_52 and RB_8 (RNA-STAR settings in Supplementary Data 2).

Bioinformatic analyses

Chromosomal copy number variations (CNVs) were identified for the 61 samples by calling CNVs from the whole-genome methylation datasets for each sample using the Bioconductor conumee package70. Whole-exome sequencing data were analyzed using Varlociraptor v3.6.3 (https://github.com/snakemake-workflows/dna-seq-varlociraptor). Consensus clustering was performed using the R package ConsensusClusterPlus (settings: maxK = 4, reps = 1000, pItem = 0.8, pFeature = 1, clusterAlg = “hc”, distance = “Pearson”). Correlations between DNA methylation and gene expression were assessed using the Pearson correlation coefficient (cutoff: variance for β-values ≥ 1 10−4, Bonferroni-Hochberg-adjusted p < 0.05, R ≥ 0.4) and the batch-corrected matrix for normalized expression values of the corresponding genes. If multiple CpGs in a gene were differentially methylated, the Pearson correlation coefficient was calculated for each CpG-gene pair. Enrichment of expression or methylation gene signatures were analyzed using the clusterProfiler Bioconductor package71 with parameter settings (ontology: ALL, p-value cutoff: 0.05, method of p-value adjustment: Benjamini-Yekutieli, adjusted p-value cutoff: 0.25, minimal gene set size: 10, maximal gene set size: 500). Methylation probe mapping utilized RefSeq gene features. Venn diagrams for gene expression were generated using the VennDiagram and BioVenn R packages72. Box and volcano plots for gene expression and methylation β-value visualizations were generated using the ggplot2 and reshape2 R packages. Bioinformatic analyses were conducted using the R statistical package, v4.2.2 (2022) and Bioconductor libraries v3.16 (R Core Team, 2020). P values were adjusted using the Benjamini–Hochberg method73 for multiple comparisons, and unless otherwise stated, considered statistically significant if p < 0.05.

Statistics and reproducibility

Data were expressed as mean ± SD or mean ± SE. Unless otherwise stated, all statistical tests are two-sided. Unless otherwise stated, p < 0.05 was regarded as statistically significant. Every experiment was repeated three times independently.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Supplementary information

Supplementary Information

Description of Additional Supplementary Files

Supplementary Data 1-2

Supplementary Data 3-5

Supplementary Data 6-16

Supplementary Data 17-24

Supplementary Data 25-29

Supplementary Data 30

Reporting summary

Supplementary information

The online version contains supplementary material available at 10.1038/s42003-024-06596-6.

Acknowledgements

The authors thank all patients who participated in this study and all physicians for performing study examinations. We thank Brenda Gallie’s Lab for providing the cell lines representative for RB1-proficient retinoblastoma. This work was supported, in part, by funding from the Medical Research School Düsseldorf, DSO, Heinrich-Heine-Universität Düsseldorf, Düsseldorf. The INFORM program is financially supported by the German Cancer Research Center (DKFZ), several German health insurance companies, the German Cancer Consortium (DKTK), the German Federal Ministry of Education and Research (BMBF), the German Federal Ministry of Health (BMG), the Ministry of Science, Research and the Arts of the State of Baden-Württemberg (MWK BW); the German Cancer Aid (DKH), the German Childhood Cancer Foundation (DKS), RTL television, the aid organization BILD hilft e.V. (Ein Herz für Kinder) and the generous private donation of the Scheu family. This research was funded by the Deutsche Forschungsgemeinschaft (TE 1187/2-1 to P.K.), Kinderaugenkrebsstiftung (KAKS20220422f to P.K.), the German Cancer Consortium (DKTK) joint funding project “task force Paediatric Oncology”, IFORES program University Duisburg-Essen (IFORES 2011 to P.K.). Welcome Back grant by University Duisburg-Essen (2023 to T.R.).

Author contributions

Conceptualization: P.K., D.L., and A.S.; methodology: P.K., M.Schw., A.S., H.H., E.A., and T.R.; bioinformatics and statistics: S.R., E.A., T.H., and C.S.; cell model generation and functional experiments: T.R., M.Schw., M.W., D.K., M.Schn., A.B., and H.H.; manuscript preparation and revision: T.R., E.A., P.K., and K.A.; contributing discussions about data and findings: P.K., D.L., A.S., K.A., D.K, L.S., E.B., S.R., K.S., B.J., and D.T.W.J.; supervision: P.K., D.L., S.R., and A.S., project administration: P.K.; funding acquisition: P.K., D.L., and T.R. All authors have read and agreed to the submitted version of the manuscript.

Peer review

Peer review information

Communications Biology thanks Jesse Berry, Siyuan Zheng and Timothy Corson for their contribution to the peer review of this work. Primary Handling Editors: Derrick Ong and David Favero.

Funding

Open Access funding enabled and organized by Projekt DEAL.

Data availability

Raw DNA methylation data and RNA sequencing data have been deposited in GEO under study accession no. GSE267015 and GSE268136. GSE5998313, GSE5878512, and GSE8698130 datasets were used for the analyses. The source data for graphs and charts are available in Supplementary Data 30. All other data (raw WES data and the code for consensus clustering) will be provided upon reasonable request.

Competing interests

The authors declare no competing interests. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Ethical approval

The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the University Duisburg-Essen Ethics Committee (17-7384-BO from 9.2.2017 and 13-5345-BO from 29.1.2013). All ethical regulations relevant to human research participants were followed.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

These authors contributed equally: Tatsiana Ryl, Elena Afanasyeva.
==== Refs
References

1. Xu XL Retinoblastoma has properties of a cone precursor tumor and depends upon cone-specific MDM2 signaling Cell 2009 137 1018 1031 10.1016/j.cell.2009.03.051 19524506
Xu, X. L. et al. Retinoblastoma has properties of a cone precursor tumor and depends upon cone-specific MDM2 signaling. Cell 137, 1018–1031 (2009).19524506 10.1016/j.cell.2009.03.051
2. Bremner R Sage J Cancer: the origin of human retinoblastoma Nature 2014 514 312 313 10.1038/nature13748 25252972
Bremner, R. & Sage, J. Cancer: the origin of human retinoblastoma. Nature 514, 312–313 (2014).25252972 10.1038/nature13748
3. Xu XL Rb suppresses human cone-precursor-derived retinoblastoma tumours Nature 2014 514 385 388 10.1038/nature13813 25252974
Xu, X. L. et al. Rb suppresses human cone-precursor-derived retinoblastoma tumours. Nature 514, 385–388 (2014).25252974 10.1038/nature13813
4. Reschke M Eye tumors in childhood as first sign of tumor predisposition syndromes: insights from an observational study conducted in Germany and Austria Cancers 2021 13 1876 10.3390/cancers13081876 33919815
Reschke, M. et al. Eye tumors in childhood as first sign of tumor predisposition syndromes: insights from an observational study conducted in Germany and Austria. Cancers 13, 1876 (2021).33919815 10.3390/cancers13081876
5. Global Retinoblastoma Study G. The Global Retinoblastoma Outcome Study: a prospective, cluster-based analysis of 4064 patients from 149 countries Lancet Glob. Health 2022 10 e1128 e1140 10.1016/S2214-109X(22)00250-9 35839812
Global Retinoblastoma Study G. The Global Retinoblastoma Outcome Study: a prospective, cluster-based analysis of 4064 patients from 149 countries. Lancet Glob. Health 10, e1128–e1140 (2022).35839812 10.1016/S2214-109X(22)00250-9
6. Berry JL Genomic cfDNA analysis of aqueous humor in retinoblastoma predicts eye salvage: the surrogate tumor biopsy for retinoblastoma Mol. Cancer Res. 2018 16 1701 1712 10.1158/1541-7786.MCR-18-0369 30061186
Berry, J. L. et al. Genomic cfDNA analysis of aqueous humor in retinoblastoma predicts eye salvage: the surrogate tumor biopsy for retinoblastoma. Mol. Cancer Res. 16, 1701–1712 (2018).30061186 10.1158/1541-7786.MCR-18-0369
7. Kooi IE Somatic genomic alterations in retinoblastoma beyond RB1 are rare and limited to copy number changes Sci. Rep. 2016 6 25264 10.1038/srep25264 27126562
Kooi, I. E. et al. Somatic genomic alterations in retinoblastoma beyond RB1 are rare and limited to copy number changes. Sci. Rep. 6, 25264 (2016).27126562 10.1038/srep25264
8. Norrie JL Retinoblastoma from human stem cell-derived retinal organoids Nat. Commun. 2021 12 4535 10.1038/s41467-021-24781-7 34315877
Norrie, J. L. et al. Retinoblastoma from human stem cell-derived retinal organoids. Nat. Commun. 12, 4535 (2021).34315877 10.1038/s41467-021-24781-7
9. Zhang J A novel retinoblastoma therapy from genomic and epigenetic analyses Nature 2012 481 329 334 10.1038/nature10733 22237022
Zhang, J. et al. A novel retinoblastoma therapy from genomic and epigenetic analyses. Nature 481, 329–334 (2012).22237022 10.1038/nature10733
10. McEvoy J RB1 gene inactivation by chromothripsis in human retinoblastoma Oncotarget 2014 5 438 450 10.18632/oncotarget.1686 24509483
McEvoy, J. et al. RB1 gene inactivation by chromothripsis in human retinoblastoma. Oncotarget 5, 438–450 (2014).24509483 10.18632/oncotarget.1686
11. Grasemann C Gains and overexpression identify DEK and E2F3 as targets of chromosome 6p gains in retinoblastoma Oncogene 2005 24 6441 6449 10.1038/sj.onc.1208792 16007192
Grasemann, C. et al. Gains and overexpression identify DEK and E2F3 as targets of chromosome 6p gains in retinoblastoma. Oncogene 24, 6441–6449 (2005).16007192 10.1038/sj.onc.1208792
12. Liu J A high-risk retinoblastoma subtype with stemness features, dedifferentiated cone states and neuronal/ganglion cell gene expression Nat. Commun. 2021 12 5578 10.1038/s41467-021-25792-0 34552068
Liu, J. et al. A high-risk retinoblastoma subtype with stemness features, dedifferentiated cone states and neuronal/ganglion cell gene expression. Nat. Commun. 12, 5578 (2021).34552068 10.1038/s41467-021-25792-0
13. Kooi IE Loss of photoreceptorness and gain of genomic alterations in retinoblastoma reveal tumor progression EBioMedicine 2015 2 660 670 10.1016/j.ebiom.2015.06.022 26288838
Kooi, I. E. et al. Loss of photoreceptorness and gain of genomic alterations in retinoblastoma reveal tumor progression. EBioMedicine 2, 660–670 (2015).26288838 10.1016/j.ebiom.2015.06.022
14. Kapatai G Gene expression profiling identifies different sub-types of retinoblastoma Br. J. cancer 2013 109 512 525 10.1038/bjc.2013.283 23756868
Kapatai, G. et al. Gene expression profiling identifies different sub-types of retinoblastoma. Br. J. cancer 109, 512–525 (2013).23756868 10.1038/bjc.2013.283
15. Busch M Metz K Beier M Biewald E Dunker N Trefoil factor family 1 expression correlates with clinical outcome in patients with retinoblastoma Retina 2018 38 2422 2428 10.1097/IAE.0000000000001881 29016460
Busch, M., Metz, K., Beier, M., Biewald, E. & Dunker, N. Trefoil factor family 1 expression correlates with clinical outcome in patients with retinoblastoma. Retina 38, 2422–2428 (2018).29016460 10.1097/IAE.0000000000001881
16. Roohollahi K de Jong Y van Mil SE Fabius AWM Moll AC Dorsman JC High-level MYCN-amplified RB1-proficient retinoblastoma tumors retain distinct molecular signatures Ophthalmol. Sci. 2022 2 100188 10.1016/j.xops.2022.100188 36245757
Roohollahi, K., de Jong, Y., van Mil, S. E., Fabius, A. W. M., Moll, A. C. & Dorsman, J. C. High-level MYCN-amplified RB1-proficient retinoblastoma tumors retain distinct molecular signatures. Ophthalmol. Sci. 2, 100188 (2022).36245757 10.1016/j.xops.2022.100188
17. Alborzinia H MYCN mediates cysteine addiction and sensitizes neuroblastoma to ferroptosis Nat. Cancer 2022 3 471 485 10.1038/s43018-022-00355-4 35484422
Alborzinia, H. et al. MYCN mediates cysteine addiction and sensitizes neuroblastoma to ferroptosis. Nat. Cancer 3, 471–485 (2022).35484422 10.1038/s43018-022-00355-4
18. Trumpp A c-Myc regulates mammalian body size by controlling cell number but not cell size Nature 2001 414 768 773 10.1038/414768a 11742404
Trumpp, A. et al. c-Myc regulates mammalian body size by controlling cell number but not cell size. Nature 414, 768–773 (2001).11742404 10.1038/414768a
19. Evan GI Vousden KH Proliferation, cell cycle and apoptosis in cancer Nature 2001 411 342 348 10.1038/35077213 11357141
Evan, G. I. & Vousden, K. H. Proliferation, cell cycle and apoptosis in cancer. Nature 411, 342–348 (2001).11357141 10.1038/35077213
20. Ryl T Cell-cycle position of single MYC-driven cancer cells dictates their susceptibility to a chemotherapeutic drug Cell Syst. 2017 5 237 250.e238 10.1016/j.cels.2017.07.005 28843484
Ryl, T. et al. Cell-cycle position of single MYC-driven cancer cells dictates their susceptibility to a chemotherapeutic drug. Cell Syst. 5, 237–250.e238 (2017).28843484 10.1016/j.cels.2017.07.005
21. Meyer N Penn LZ Reflecting on 25 years with MYC Nat. Rev. Cancer 2008 8 976 990 10.1038/nrc2231 19029958
Meyer, N. & Penn, L. Z. Reflecting on 25 years with MYC. Nat. Rev. Cancer 8, 976–990 (2008).19029958 10.1038/nrc2231
22. Schwab M Westermann F Hero B Berthold F Neuroblastoma: biology and molecular and chromosomal pathology Lancet Oncol. 2003 4 472 480 10.1016/S1470-2045(03)01166-5 12901961
Schwab, M., Westermann, F., Hero, B. & Berthold, F. Neuroblastoma: biology and molecular and chromosomal pathology. Lancet Oncol. 4, 472–480 (2003).12901961 10.1016/S1470-2045(03)01166-5
23. Westermann F Distinct transcriptional MYCN/c-MYC activities are associated with spontaneous regression or malignant progression in neuroblastomas Genome Biol. 2008 9 R150 10.1186/gb-2008-9-10-r150 18851746
Westermann, F. et al. Distinct transcriptional MYCN/c-MYC activities are associated with spontaneous regression or malignant progression in neuroblastomas. Genome Biol. 9, R150 (2008).18851746 10.1186/gb-2008-9-10-r150
24. Northcott PA The whole-genome landscape of medulloblastoma subtypes Nature 2017 547 311 317 10.1038/nature22973 28726821
Northcott, P. A. et al. The whole-genome landscape of medulloblastoma subtypes. Nature 547, 311–317 (2017).28726821 10.1038/nature22973
25. Rushlow DE Characterisation of retinoblastomas without RB1 mutations: genomic, gene expression, and clinical studies lancet Oncol. 2013 14 327 334 10.1016/S1470-2045(13)70045-7 23498719
Rushlow, D. E. et al. Characterisation of retinoblastomas without RB1 mutations: genomic, gene expression, and clinical studies. lancet Oncol. 14, 327–334 (2013).23498719 10.1016/S1470-2045(13)70045-7
26. Blixt MKE MYCN induces cell-specific tumorigenic growth in RB1-proficient human retinal organoid and chicken retina models of retinoblastoma Oncogenesis 2022 11 34 10.1038/s41389-022-00409-3 35729105
Blixt, M. K. E. et al. MYCN induces cell-specific tumorigenic growth in RB1-proficient human retinal organoid and chicken retina models of retinoblastoma. Oncogenesis 11, 34 (2022).35729105 10.1038/s41389-022-00409-3
27. Schroder C Rahmann S A hybrid parameter estimation algorithm for beta mixtures and applications to methylation state classification Algorithms Mol. Biol. 2017 12 21 10.1186/s13015-017-0112-1 28828033
Schroder, C. & Rahmann, S. A hybrid parameter estimation algorithm for beta mixtures and applications to methylation state classification. Algorithms Mol. Biol. 12, 21 (2017).28828033 10.1186/s13015-017-0112-1
28. Li HT Characterizing DNA methylation signatures of retinoblastoma using aqueous humor liquid biopsy Nat. Commun. 2022 13 5523 10.1038/s41467-022-33248-2 36130950
Li, H. T. et al. Characterizing DNA methylation signatures of retinoblastoma using aqueous humor liquid biopsy. Nat. Commun. 13, 5523 (2022).36130950 10.1038/s41467-022-33248-2
29. Pimentel H Bray NL Puente S Melsted P Pachter L Differential analysis of RNA-seq incorporating quantification uncertainty Nat. Methods 2017 14 687 690 10.1038/nmeth.4324 28581496
Pimentel, H., Bray, N. L., Puente, S., Melsted, P. & Pachter, L. Differential analysis of RNA-seq incorporating quantification uncertainty. Nat. Methods 14, 687–690 (2017).28581496 10.1038/nmeth.4324
30. Aldiri I The dynamic epigenetic landscape of the retina during development, reprogramming, and tumorigenesis Neuron 2017 94 550 568.e510 10.1016/j.neuron.2017.04.022 28472656
Aldiri, I. et al. The dynamic epigenetic landscape of the retina during development, reprogramming, and tumorigenesis. Neuron 94, 550–568.e510 (2017).28472656 10.1016/j.neuron.2017.04.022
31. Jiang Y SEdb: a comprehensive human super-enhancer database Nucleic Acids Res. 2019 47 D235 D243 10.1093/nar/gky1025 30371817
Jiang, Y. et al. SEdb: a comprehensive human super-enhancer database. Nucleic Acids Res. 47, D235–D243 (2019).30371817 10.1093/nar/gky1025
32. Schwermer M Comprehensive characterization of RB1 mutant and MYCN amplified retinoblastoma cell lines Exp. Cell Res. 2019 375 92 99 10.1016/j.yexcr.2018.12.018 30584916
Schwermer, M. et al. Comprehensive characterization of RB1 mutant and MYCN amplified retinoblastoma cell lines. Exp. Cell Res. 375, 92–99 (2019).30584916 10.1016/j.yexcr.2018.12.018
33. Schuhmacher M The transcriptional program of a human B cell line in response to Myc Nucleic Acids Res. 2001 29 397 406 10.1093/nar/29.2.397 11139609
Schuhmacher, M. et al. The transcriptional program of a human B cell line in response to Myc. Nucleic Acids Res. 29, 397–406 (2001).11139609 10.1093/nar/29.2.397
34. Wei JS The MYCN oncogene is a direct target of miR-34a Oncogene 2008 27 5204 5213 10.1038/onc.2008.154 18504438
Wei, J. S. et al. The MYCN oncogene is a direct target of miR-34a. Oncogene 27, 5204–5213 (2008).18504438 10.1038/onc.2008.154
35. Zeller KI Jegga AG Aronow BJ O’Donnell KA Dang CV An integrated database of genes responsive to the Myc oncogenic transcription factor: identification of direct genomic targets Genome Biol. 2003 4 R69 10.1186/gb-2003-4-10-r69 14519204
Zeller, K. I., Jegga, A. G., Aronow, B. J., O’Donnell, K. A. & Dang, C. V. An integrated database of genes responsive to the Myc oncogenic transcription factor: identification of direct genomic targets. Genome Biol. 4, R69 (2003).14519204 10.1186/gb-2003-4-10-r69
36. Zimmerman MW MYC drives a subset of high-risk pediatric neuroblastomas and is activated through mechanisms including enhancer hijacking and focal enhancer amplification Cancer Discov. 2018 8 320 335 10.1158/2159-8290.CD-17-0993 29284669
Zimmerman, M. W. et al. MYC drives a subset of high-risk pediatric neuroblastomas and is activated through mechanisms including enhancer hijacking and focal enhancer amplification. Cancer Discov. 8, 320–335 (2018).29284669 10.1158/2159-8290.CD-17-0993
37. Misiak D The microRNA landscape of MYCN-amplified neuroblastoma Front Oncol. 2021 11 647737 10.3389/fonc.2021.647737 34026620
Misiak, D. et al. The microRNA landscape of MYCN-amplified neuroblastoma. Front Oncol. 11, 647737 (2021).34026620 10.3389/fonc.2021.647737
38. Murphy DM Global MYCN transcription factor binding analysis in neuroblastoma reveals association with distinct E-box motifs and regions of DNA hypermethylation PLoS ONE 2009 4 e8154 10.1371/journal.pone.0008154 19997598
Murphy, D. M. et al. Global MYCN transcription factor binding analysis in neuroblastoma reveals association with distinct E-box motifs and regions of DNA hypermethylation. PLoS ONE 4, e8154 (2009).19997598 10.1371/journal.pone.0008154
39. Olender T Lancet D Nebert DW Update on the olfactory receptor (OR) gene superfamily Hum. Genomics 2008 3 87 97 10.1186/1479-7364-3-1-87 19129093
Olender, T., Lancet, D. & Nebert, D. W. Update on the olfactory receptor (OR) gene superfamily. Hum. Genomics 3, 87–97 (2008).19129093 10.1186/1479-7364-3-1-87
40. Kalra S Analysis of single-cell transcriptomes links enrichment of olfactory receptors with cancer cell differentiation status and prognosis Commun. Biol. 2020 3 506 10.1038/s42003-020-01232-5 32917933
Kalra, S. et al. Analysis of single-cell transcriptomes links enrichment of olfactory receptors with cancer cell differentiation status and prognosis. Commun. Biol. 3, 506 (2020).32917933 10.1038/s42003-020-01232-5
41. Brodeur GM Seeger RC Schwab M Varmus HE Bishop JM Amplification of N-myc in untreated human neuroblastomas correlates with advanced disease stage Science 1984 224 1121 1124 10.1126/science.6719137 6719137
Brodeur, G. M., Seeger, R. C., Schwab, M., Varmus, H. E. & Bishop, J. M. Amplification of N-myc in untreated human neuroblastomas correlates with advanced disease stage. Science 224, 1121–1124 (1984).6719137 10.1126/science.6719137
42. Singh HP An immature, dedifferentiated, and lineage-deconstrained cone precursor origin of N-Myc-initiated retinoblastoma Proc. Natl Acad. Sci. USA 2022 119 e2200721119 10.1073/pnas.2200721119 35867756
Singh, H. P. et al. An immature, dedifferentiated, and lineage-deconstrained cone precursor origin of N-Myc-initiated retinoblastoma. Proc. Natl Acad. Sci. USA 119, e2200721119 (2022).35867756 10.1073/pnas.2200721119
43. Henrich KO Integrative genome-scale analysis identifies epigenetic mechanisms of transcriptional deregulation in unfavorable neuroblastomas Cancer Res. 2016 76 5523 5537 10.1158/0008-5472.CAN-15-2507 27635046
Henrich, K. O. et al. Integrative genome-scale analysis identifies epigenetic mechanisms of transcriptional deregulation in unfavorable neuroblastomas. Cancer Res. 76, 5523–5537 (2016).27635046 10.1158/0008-5472.CAN-15-2507
44. Boon K N-myc enhances the expression of a large set of genes functioning in ribosome biogenesis and protein synthesis EMBO J. 2001 20 1383 1393 10.1093/emboj/20.6.1383 11250904
Boon, K. et al. N-myc enhances the expression of a large set of genes functioning in ribosome biogenesis and protein synthesis. EMBO J. 20, 1383–1393 (2001).11250904 10.1093/emboj/20.6.1383
45. Hald OH Inhibitors of ribosome biogenesis repress the growth of MYCN-amplified neuroblastoma Oncogene 2019 38 2800 2813 10.1038/s41388-018-0611-7 30542116
Hald, O. H. et al. Inhibitors of ribosome biogenesis repress the growth of MYCN-amplified neuroblastoma. Oncogene 38, 2800–2813 (2019).30542116 10.1038/s41388-018-0611-7
46. Price EA MYCN amplification levels in primary retinoblastoma tumors analyzed by Multiple Ligation-dependent Probe Amplification Ophthalmic Genet. 2021 42 604 611 10.1080/13816810.2021.1923038 34003079
Price, E. A. et al. MYCN amplification levels in primary retinoblastoma tumors analyzed by Multiple Ligation-dependent Probe Amplification. Ophthalmic Genet. 42, 604–611 (2021).34003079 10.1080/13816810.2021.1923038
47. Kang JH Rychahou PG Ishola TA Qiao J Evers BM Chung DH MYCN silencing induces differentiation and apoptosis in human neuroblastoma cells Biochem. Biophys. Res. Commun. 2006 351 192 197 10.1016/j.bbrc.2006.10.020 17055458
Kang, J. H., Rychahou, P. G., Ishola, T. A., Qiao, J., Evers, B. M. & Chung, D. H. MYCN silencing induces differentiation and apoptosis in human neuroblastoma cells. Biochem. Biophys. Res. Commun. 351, 192–197 (2006).17055458 10.1016/j.bbrc.2006.10.020
48. Cotterman R Knoepfler PS N-Myc regulates expression of pluripotency genes in neuroblastoma including lif, klf2, klf4, and lin28b PloS one 2009 4 e5799 10.1371/journal.pone.0005799 19495417
Cotterman, R. & Knoepfler, P. S. N-Myc regulates expression of pluripotency genes in neuroblastoma including lif, klf2, klf4, and lin28b. PloS one 4, e5799 (2009).19495417 10.1371/journal.pone.0005799
49. Henriksen JR Conditional expression of retrovirally delivered anti-MYCN shRNA as an in vitro model system to study neuronal differentiation in MYCN-amplified neuroblastoma BMC Dev. Biol. 2011 11 1 10.1186/1471-213X-11-1 21194500
Henriksen, J. R. et al. Conditional expression of retrovirally delivered anti-MYCN shRNA as an in vitro model system to study neuronal differentiation in MYCN-amplified neuroblastoma. BMC Dev. Biol. 11, 1 (2011).21194500 10.1186/1471-213X-11-1
50. Reddy CD Anticancer effects of the novel 1alpha, 25-dihydroxyvitamin D3 hybrid analog QW1624F2-2 in human neuroblastoma J. Cell Biochem. 2006 97 198 206 10.1002/jcb.20629 16200638
Reddy, C. D. et al. Anticancer effects of the novel 1alpha, 25-dihydroxyvitamin D3 hybrid analog QW1624F2-2 in human neuroblastoma. J. Cell Biochem. 97, 198–206 (2006).16200638 10.1002/jcb.20629
51. Bartolucci D MYCN impact on high-risk neuroblastoma: from diagnosis and prognosis to targeted treatment Cancers 2022 14 4421 10.3390/cancers14184421 36139583
Bartolucci, D. et al. MYCN impact on high-risk neuroblastoma: from diagnosis and prognosis to targeted treatment. Cancers 14, 4421 (2022).36139583 10.3390/cancers14184421
52. Fletcher JI Ziegler DS Trahair TN Marshall GM Haber M Norris MD Too many targets, not enough patients: rethinking neuroblastoma clinical trials Nat. Rev. Cancer 2018 18 389 400 10.1038/s41568-018-0003-x 29632319
Fletcher, J. I., Ziegler, D. S., Trahair, T. N., Marshall, G. M., Haber, M. & Norris, M. D. Too many targets, not enough patients: rethinking neuroblastoma clinical trials. Nat. Rev. Cancer 18, 389–400 (2018).29632319 10.1038/s41568-018-0003-x
53. Xu L Establishing the clinical utility of ctDNA analysis for diagnosis, prognosis, and treatment monitoring of retinoblastoma: the aqueous humor liquid biopsy Cancers 2021 13 1282 10.3390/cancers13061282 33805776
Xu, L. et al. Establishing the clinical utility of ctDNA analysis for diagnosis, prognosis, and treatment monitoring of retinoblastoma: the aqueous humor liquid biopsy. Cancers 13, 1282 (2021).33805776 10.3390/cancers13061282
54. van Tilburg CM The pediatric precision oncology INFORM registry: clinical outcome and benefit for patients with very high-evidence targets Cancer Discov. 2021 11 2764 2779 10.1158/2159-8290.CD-21-0094 34373263
van Tilburg, C. M. et al. The pediatric precision oncology INFORM registry: clinical outcome and benefit for patients with very high-evidence targets. Cancer Discov. 11, 2764–2779 (2021).34373263 10.1158/2159-8290.CD-21-0094
55. Heipertz AE Outcome of children and adolescents with relapsed/refractory/progressive malignancies treated with molecularly informed targeted drugs in the pediatric precision oncology registry INFORM JCO Precis. Oncol. 2023 7 e2300015 10.1200/PO.23.00015 37364231
Heipertz, A. E. et al. Outcome of children and adolescents with relapsed/refractory/progressive malignancies treated with molecularly informed targeted drugs in the pediatric precision oncology registry INFORM. JCO Precis. Oncol. 7, e2300015 (2023).37364231 10.1200/PO.23.00015
56. Peterziel H Drug sensitivity profiling of 3D tumor tissue cultures in the pediatric precision oncology program INFORM NPJ Precis Oncol. 2022 6 94 10.1038/s41698-022-00335-y 36575299
Peterziel, H. et al. Drug sensitivity profiling of 3D tumor tissue cultures in the pediatric precision oncology program INFORM. NPJ Precis Oncol. 6, 94 (2022).36575299 10.1038/s41698-022-00335-y
57. Worst BC Next-generation personalised medicine for high-risk paediatric cancer patients—the INFORM pilot study Eur. J. Cancer 2016 65 91 101 10.1016/j.ejca.2016.06.009 27479119
Worst, B. C. et al. Next-generation personalised medicine for high-risk paediatric cancer patients—the INFORM pilot study. Eur. J. Cancer 65, 91–101 (2016).27479119 10.1016/j.ejca.2016.06.009
58. Wiek C Identification of amino acid determinants in CYP4B1 for optimal catalytic processing of 4-ipomeanol Biochem. J. 2015 465 103 114 10.1042/BJ20140813 25247810
Wiek, C. et al. Identification of amino acid determinants in CYP4B1 for optimal catalytic processing of 4-ipomeanol. Biochem. J. 465, 103–114 (2015).25247810 10.1042/BJ20140813
59. Berg S ilastik: interactive machine learning for (bio)image analysis Nat. Methods 2019 16 1226 1232 10.1038/s41592-019-0582-9 31570887
Berg, S. et al. ilastik: interactive machine learning for (bio)image analysis. Nat. Methods 16, 1226–1232 (2019).31570887 10.1038/s41592-019-0582-9
60. Albrecht P Ansperger-Rescher B Schuler A Zeschnigk M Gallie B Lohmann DR Spectrum of gross deletions and insertions in the RB1 gene in patients with retinoblastoma and association with phenotypic expression Hum. Mutat. 2005 26 437 445 10.1002/humu.20234 16127685
Albrecht, P., Ansperger-Rescher, B., Schuler, A., Zeschnigk, M., Gallie, B. & Lohmann, D. R. Spectrum of gross deletions and insertions in the RB1 gene in patients with retinoblastoma and association with phenotypic expression. Hum. Mutat. 26, 437–445 (2005).16127685 10.1002/humu.20234
61. Lohmann DR Constitutional RB1-gene mutations in patients with isolated unilateral retinoblastoma Am. J. Hum. Genet. 1997 61 282 294 10.1086/514845 9311732
Lohmann, D. R. et al. Constitutional RB1-gene mutations in patients with isolated unilateral retinoblastoma. Am. J. Hum. Genet. 61, 282–294 (1997).9311732 10.1086/514845
62. Zeschnigk M Bohringer S Price EA Onadim Z Masshofer L Lohmann DR A novel real-time PCR assay for quantitative analysis of methylated alleles (QAMA): analysis of the retinoblastoma locus Nucleic Acids Res. 2004 32 e125 10.1093/nar/gnh122 15353561
Zeschnigk, M., Bohringer, S., Price, E. A., Onadim, Z., Masshofer, L. & Lohmann, D. R. A novel real-time PCR assay for quantitative analysis of methylated alleles (QAMA): analysis of the retinoblastoma locus. Nucleic Acids Res. 32, e125 (2004).15353561 10.1093/nar/gnh122
63. Zhang K Nowak I Rushlow D Gallie BL Lohmann DR Patterns of missplicing caused by RB1 gene mutations in patients with retinoblastoma and association with phenotypic expression Hum. Mutat. 2008 29 475 484 10.1002/humu.20664 18181215
Zhang, K., Nowak, I., Rushlow, D., Gallie, B. L. & Lohmann, D. R. Patterns of missplicing caused by RB1 gene mutations in patients with retinoblastoma and association with phenotypic expression. Hum. Mutat. 29, 475–484 (2008).18181215 10.1002/humu.20664
64. Afgan E The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update Nucleic Acids Res. 2018 46 W537 W544 10.1093/nar/gky379 29790989
Afgan, E. et al. The Galaxy platform for accessible, reproducible and collaborative biomedical analyses: 2018 update. Nucleic Acids Res. 46, W537–W544 (2018).29790989 10.1093/nar/gky379
65. Bray NL Pimentel H Melsted P Pachter L Near-optimal probabilistic RNA-seq quantification Nat. Biotechnol. 2016 34 525 527 10.1038/nbt.3519 27043002
Bray, N. L., Pimentel, H., Melsted, P. & Pachter, L. Near-optimal probabilistic RNA-seq quantification. Nat. Biotechnol. 34, 525–527 (2016).27043002 10.1038/nbt.3519
66. Storey JD Tibshirani R Statistical significance for genomewide studies Proc. Natl Acad. Sci. USA 2003 100 9440 9445 10.1073/pnas.1530509100 12883005
Storey, J. D. & Tibshirani, R. Statistical significance for genomewide studies. Proc. Natl Acad. Sci. USA 100, 9440–9445 (2003).12883005 10.1073/pnas.1530509100
67. Edgington ES An additive method for combining probability values from independent experiments J. Psychol. 1972 80 351 363 10.1080/00223980.1972.9924813
Edgington, E. S. An additive method for combining probability values from independent experiments. J. Psychol. 80, 351–363 (1972).10.1080/00223980.1972.9924813
68. Xiao Y A novel significance score for gene selection and ranking Bioinformatics 2014 30 801 807 10.1093/bioinformatics/btr671 22321699
Xiao, Y. et al. A novel significance score for gene selection and ranking. Bioinformatics 30, 801–807 (2014).22321699 10.1093/bioinformatics/btr671
69. Liberzon A Subramanian A Pinchback R Thorvaldsdottir H Tamayo P Mesirov JP Molecular signatures database (MSigDB) 3.0 Bioinformatics 2011 27 1739 1740 10.1093/bioinformatics/btr260 21546393
Liberzon, A., Subramanian, A., Pinchback, R., Thorvaldsdottir, H., Tamayo, P. & Mesirov, J. P. Molecular signatures database (MSigDB) 3.0. Bioinformatics 27, 1739–1740 (2011).21546393 10.1093/bioinformatics/btr260
70. Hovestadt, V. & Zapat, M. conumee: Enhanced Copy-Number Variation Analysis using Illumina DNA Methylation Arrays. R Package Version 1.9.0 (2015).
71. Yu G Wang LG Han Y He QY clusterProfiler: an R package for comparing biological themes among gene clusters Omics 2012 16 284 287 10.1089/omi.2011.0118 22455463
Yu, G., Wang, L. G., Han, Y. & He, Q. Y. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287 (2012).22455463 10.1089/omi.2011.0118
72. Hulsen T de Vlieg J Alkema W BioVenn—a web application for the comparison and visualization of biological lists using area-proportional Venn diagrams BMC Genom. 2008 9 488 10.1186/1471-2164-9-488
Hulsen, T., de Vlieg, J. & Alkema, W. BioVenn—a web application for the comparison and visualization of biological lists using area-proportional Venn diagrams. BMC Genom. 9, 488 (2008).10.1186/1471-2164-9-488
73. Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B 57, 289–300 (1995).
