A high-risk retinoblastoma subtype with stemness features, dedifferentiated cone states and neuronal/ganglion cell gene expression

Abstract

Retinoblastoma is the most frequent intraocular malignancy in children, originating from a maturing cone precursor in the developing retina. Little is known on the molecular basis underlying the biological and clinical behavior of this cancer. Here, using multi-omics data, we demonstrate the existence of two retinoblastoma subtypes. Subtype 1, of earlier onset, includes most of the heritable forms. It harbors few genetic alterations other than the initiating RB1 inactivation and corresponds to differentiated tumors expressing mature cone markers. By contrast, subtype 2 tumors harbor frequent recurrent genetic alterations including MYCN-amplification. They express markers of less differentiated cone together with neuronal/ganglion cell markers with marked inter- and intra-tumor heterogeneity. The cone dedifferentiation in subtype 2 is associated with stemness features including low immune and interferon response, E2F and MYC/MYCN activation and a higher propensity for metastasis. The recognition of these two subtypes, one maintaining a cone-differentiated state, and the other, more aggressive, associated with cone dedifferentiation and expression of neuronal markers, opens up important biological and clinical perspectives for retinoblastomas.

Fig. 1. Multi-omics-based molecular subtypes of retinoblastoma and clinical characteristics.

a Consensus clustering of retinoblastomas based on transcriptomic, DNA methylation, and copy-number alteration data (top panel). Unsupervised cluster-of-clusters analysis (middle panel). Supervised centroid-based classification (bottom panel). Final omics subtype: subtype 1, n = 31 (gold); subtype 2, n = 38 (blue); unclassified, n = 3 (gray). b Heatmap showing methylation values (methylome arrays) for the nine-CpG-based classifier (left panel). Correlation between methylation values assessed by pyrosequencing and by methylome array, for 17 tumors (middle panel). A two-sided Pearson’s correlation test was used. The nine-CpG-based classifier applied to a subset of 17 tumors of the initial series, led to the same classification as obtained by the -omics approach in 16 cases (one case being not classified by the nine-CpG-based classifier). Subtype assignment of 30 additional tumors based on the nine-CpG-based classifier (right panel). c Final molecular classification of 96 retinoblastomas and their key clinical and pathological characteristics. p ≥ 0.05 (ns), p < 0.05 (*), p < 0.01 (**), p < 0.001 (***), p < 0.0001 (****). For comparisons of RB1 germline mutation, laterity, growth pattern, tumor diameter, and necrosis between two subtypes, Chi² tests were used. For comparisons of age at diagnosis and tumor diameter between two subtypes, two-sided Kruskal–Wallis rank tests were used. For comparisons of optic nerve invasion and choroid and sclera invasion between two subtypes, two-sided Fisher’s exact tests were used. Exact p-values are provided in Table 1.

Fig. 2. Genomic characterization, somatic mutational landscape, and DNA methylation profiles of the two retinoblastoma subtypes.

a Pattern of somatic copy-number alterations in subtype 1 (top, n = 38) and subtype 2 (bottom, n = 58) retinoblastomas. b Boxplots comparing genomic instability between subtype 1 tumors (n = 38) and subtype 2 tumors (n = 58). Among the subtype 2 tumors, non-MYCN-amplified (n = 48) and MYCN-amplified (n = 10) tumors are also shown. Significant differences were tested by two-sided Wilcoxon tests for Subtype 1 vs Subtype 2: p = 3.3 × 10⁻⁷; Subtype 1 vs Subtype 2 non-MYCN: p = 1.2 × 10⁻⁷; Subtype 1 vs Subtype 2 MYCN-amplified: p = 0.147; and Subtype 2 non-MYCN-amplified vs Subtype 2 MYCN-amplified: p = 0.014. c Boxplots comparing the number of somatic mutations between subtype 1 tumors (n = 25) and subtype 2 tumors (n = 41). Among the subtype 2 tumors, non-MYCN-amplified (n = 33) and MYCN-amplified (n = 8) tumors are also shown. Significance differences were tested by two-sided Wilcoxon tests for Subtype 1 vs Subtype 2: p = 8.1 × 10⁻⁷; Subtype 1 vs Subtype 2 non-MYCN-amplified: p = 3.5 × 10⁻⁶; Subtype 1 vs Subtype 2 MYCN-amplified: p = 0.001; and Subtype 2 non-MYCN-amplified vs Subtype 2 MYCN-amplified: p = 0.775. b, c In the boxplots, the central mark indicates the median and the bottom and top edges of the box the 25th and 75th percentiles. The whiskers are the smaller of 1.5 times the interquartile range or the length of the 25th percentiles to the smallest data point or the 75th percentiles to the largest data point. Data points outside the whiskers are outliers. Note: p ≥ 0.05 (ns), p < 0.05 (*), p < 0.01 (**), p < 0.001 (***), p < 0.0001 (****). d Somatic mutations of the three genes recurrently altered by tumor subtype. For RB1 are indicated the germline mutations. MYCN amplifications, 1q gains, and 16q losses are also shown. e Heatmap of the 6607 differentially methylated CpGs (difference of methylation level >0.2, adjusted p < 0.05, two-sided Wilcoxon test and BH correction) between subtype 1 and subtype 2. f Distribution, in subtype 2 as compared to subtype 1, of hypomethylated CpGs (upper panel) and hypermethylated CpGs (lower panel), by CpG content and neighborhood context. g Density plots showing the distribution of methylation levels of the differentially methylated CpGs located in CpG islands (upper panel) and outside CpG islands (lower panel).

Fig. 3. Transcriptomic differences between the two retinoblastoma subtypes.

a Volcano plot with genes significantly upregulated in subtype 1 (n = 26) (gold) and subtype 2 (n = 31) (blue). The genes related to cone-cell and neuronal/ganglion-cell differentiation are indicated (in gold and blue, respectively), together with the most highly differentially expressed genes in each subtype. b Hierarchical clustering of the significantly differentially expressed genes identified three main gene clusters. c Upper panels: Gene sets from the GOBP collection enriched in clusters 1.1, 1.2, 2 in hypergeometric tests. Results are presented as networks of enriched gene sets (nodes) connected based on their overlapping genes (edges). Node size is proportional to the total number of genes in the gene set concerned. The names of the various GOBP terms are given in Supplementary Data 3. Bottom panels: Top 5 Gene sets from the HALLMARK collection enriched in clusters 1.1, 1.2, 2. d Upper panel: Boxplots of stemness indices, determined as in Malta et al., in the two subtypes of retinoblastoma (subtype 1 tumors: n = 26, subtype 2 tumors: n = 31). In the boxplots, the central mark indicates the median and the bottom and top edges of the box the 25th and 75th percentiles. Whiskers are the smaller of 1.5 times the interquartile range or the length of the 25th percentiles to the smallest data point or the 75th percentiles to the largest data point. Data points outside the whiskers are outliers. Significance was tested by a two-sided Wilcoxon test, p = 1.9 × 10⁻⁷. Bottom panel: Heatmap of stemness indices and meta-score of the most correlated and anti-correlated HALLMARK (HM) pathways and MCP-score of the most anti-correlated immune cells. Spearman’s rho and p-value are shown in the figure. p < 0.0001 (****). e Heatmap representing expression pattern of cone- and ganglion-associated genes in the two subtypes of retinoblastoma. Statistical significance and log2 fold-change in expression between subtype 2 and subtype 1 are also shown. Adjusted.p ≥ 0.05 (ns), adjusted.p < 0.05 (*), adjusted.p < 0.01 (**), adjusted.p < 0.001 (***), adjusted.p < 0.0001 (****). Limma moderated two-sided t-tests and BH correction were used. Exact p-values are provided in Supplementary Data 3.

Fig. 4. Expression of cone and neuronal/ganglion cell markers in retinoblastoma and retinal organoids.

a Heatmap showing the expression of cone and ganglion markers in retinal organoids at different differentiation time points, and in subtype 1 and subtype 2 tumors assessed by NanoString technology. Differences in gene expression between the two subtypes were assessed by two-sided t-tests with BH correction. Exact p-values are provided in Supplementary Data 4. b Pearson’s correlation of the expression of 8 cone markers, between the centroids of the 2 retinoblastoma subtypes and retinal organoids at different time points in differentiation. C1: centroid of subtype 1; C2: centroid of subtype 2. c Phylogenetic tree based on cone marker expression, for retinal organoids at different differentiation time points and for retinoblastoma samples. d Immunohistochemical staining of CRX, ARR3, EBF3, and Ki-67 in normal retina and retinoblastoma. For RB617, the black arrows indicate the mutually exclusive patterns for ARR3 and EBF3. Immunohistochemistry experiments were performed on 34 samples (subtype 1, n = 9; subtype 2, n = 25). Two representative images are shown for each subtype. e Boxplots showing the quick score (QS) for the differentiation markers used in the immunohistochemical analysis: CRX, ARR3, and EBF3. In the boxplots, the central mark indicates the median and the bottom and top edges of the box the 25th and 75th percentiles. The whiskers are the smaller of 1.5 times the interquartile range or the length of the 25th percentiles to the smallest data point or the 75th percentiles to the largest data point. Data points outside the whiskers are outliers. Two-sided Wilcoxon tests were used.

Fig. 5. Intratumor heterogeneity at the single-cell level of a subtype 2 retinoblastoma (RBSC11).

a 2D t-SNE plot of 1198 single retinoblastoma cells from one patient. Each dot represents one cell. b Heatmap of top cluster markers (top 20 most upregulated genes per cluster according to fold-change). Representative cluster markers and enriched gene sets are shown. Cluster marker p-values were calculated by hypergeometric tests with BH correction. c Expression of selected genes shown in 2D t-SNE plot (early photoreceptor markers: CRX, OTX2; late cone markers: ARR3, GUCA1C; neuronal/ganglion markers: EBF3, GAP43, DCX; proliferation marker: MKI67; pro-apoptotic marker: BNIP3; macrophage marker: CD14; T-cell marker: CD3D). d CNV profiles inferred from single-cell gene expression. Each row represents the profile of one individual cell. The genes on chromosome 6p overexpressed in the non-malignant cells monocyte/microglia correspond to HLA complex genes and should not be interpreted as CNV in cluster 5. e Upper panel: Diagram summarizing the interpretation of the different clusters of the 2D t-SNE plot. Lower panel: A progression model for this retinoblastoma case based on genomic alterations.

Fig. 6. Subtype 2 tumors are associated with a higher risk of metastasis.

a Immunostaining of CRX, ARR3, and TFF1 in normal retina and retinoblastoma. Immunohistochemistry experiments were performed on 55 samples (subtype 1, n = 18; subtype 2, n = 37) from the initial series of 102 retinoblastomas. Representative images are shown: one subtype 1 tumor (RB1) and two subtype 2 tumors (RB635, RBsjd8). The subtype 2 tumors presented either a co-staining (RB635) or a mirror pattern (RBsjd8) for ARR3 and TFF1. b Boxplots showing the quick score (QS) for TFF1 in 55 tumors of the initial series (subtype 1, n = 18; subtype 2, n = 37), and in 112 tumors of the HRPF series. In the boxplots, the central mark indicates the median and the bottom and top edges of the box the 25th and 75th percentiles. The whiskers are the smaller of 1.5 times the interquartile range or the length of the 25th percentiles to the smallest data point or the 75th percentiles to the largest data point. Data points outside the whiskers are outliers. Two-sided Wilcoxon tests were used to assess the difference of the QS for Subtype 1 vs Subtype 2, p = 1.1 × 10⁻⁷, and metastatic vs non-metastatic, p = 0.007. c Immunostaining of TFF1 for primary tumors of metastatic retinoblastoma (left) and their metastatic sites (right), at low and high magnification. TFF1 expression could be assessed by immunohistochemistry for 6 of 7 available primary-metastasis tumor pairs. Representative images of four are shown.