A MYCN-driven de-differentiation profile identifies a subgroup of aggressive retinoblastoma

Abstract

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.

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^−/−MYCN^A (pink, n = 2), RB1-proficient MYCN^A (purple, n = 6), and RB1^−/−non-MYCN^A (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 classifier (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.

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 score: 3105 genes positively expressed; 3088 genes downregulated, adjusted p ≤ 0.05; or Kooi-score: 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.

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).

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.

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 retinoblastoma^, 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.

Fig. 7

Overview of characteristics of 3 retinoblastoma clusters (including RB1-proficient retinoblastoma).