Transcriptome analysis reveals molecularly distinct subtypes in retinoblastoma

Identification of molecular subtypes in retinoblastoma

To investigate the different retinoblastoma molecular subtypes, we analysed the transcriptomes of 114 retinoblastoma. Consensus clustering of all samples were used for the subtyping of retinoblastoma. Consensus average linkage hierarchical clustering of 114 samples identified two clearly robust subtypes with clustering stability increasing. We also noticed there was a heterogeneity within the two subtypes and consensus matrix showed a clear separation between subgroups and there was a minimal overlap across subgroups when k = 6 (Supplementary Fig. 1). Notably, subtype 1 can be categorized into four subgroups (subtype 1a, subtype 1b, subtype 1c and subtype 1d) while subtype 2 was divided into two subgroups (subtype 2a and subtype 2b) (Supplementary Table 1). We then queried the differentially expressed genes (DEGs) between subtypes to uncover the biological differences. A total of 56 genes were identified as DEGs between subtype 1 and subtype 2 (Supplementary Table 2). Of them, 39 genes were up-regulated in the subtype 1. Of these 39 subtype 1 genes, we can still find the differences between subgroups. Photoreceptor genes (ARR3, GUCA1B) were relatively enriched in subtype 1a and subtype 1d, while WNT pathway related gene (WIF1) were relatively high in subtype 1b (Fig. 1). Subtype 1c had the moderate expression of all subtype 1 genes. We also noticed a high level of immune response genes (CD14, CD163, HLA-DMA, HLA-DRA and HLA-DPA1) were presented in the subtype 1b, suggesting this subgroup could have a higher immune response and immune cell infiltration. Within subtype 2, subtype 2a and 2b had a higher expression of ganglion genes (EBF3 and GAP43) but subtype 2b has a relative low expression of neurodevelopment genes (DCX, ROBO1, ST6GALNAC5, TFF1). Our subtyping was also consistent with the previous study, where one subtype 1 and 2 had a higher expression of cone and ganglion markers, respectively (Fig. 1).

Consensus clustering showing two major subtypes and six subgroups in retinoblastoma. (A) Heatmap showing the different expressed genes between subtypes. (B) Volcano plot showing the different expressed genes between two subtypes. Red dots, different expressed genes; Purple dots, ganglion related genes; Green dots, cone related genes. (C) Heatmap showing the cone and ganglion genes between subtypes. (D) Barplot showing the percentage of RB1 germline mutation and Growth type between two major subtypes.

We also noticed that higher percentage of RB1 germline mutation was observed in the Subtype1, especially in Subtype 1b (Fig. 1A,D), as compared with Subtype2. Exophytic growth was predominately found in Subtype1 while endophytic growth was the major growth type in subtype2 (Fig. 1D). We also noticed that young patients (< 24 months) were predominately in Subtype1, especially in Subtype1a and 1b. More female patients were classified into subtype 2a (Supplementary Table 3).

To further explore the difference of biological function between subtypes, we then investigated the biological pathways by DEGs. Dotplot for DEGs shows immune response, epithelial cell proliferation and visual perception were the most significantly different between subtype 1 and 2 (Fig. 2). This result was confirmed by the emap plot, where the significant pathways were categorized into three main clusters. These results suggested that the biological difference between subtypes of retinoblastoma were mainly in the immune response, visual perception and epithelial cell proliferation.

Gene Ontology analysis of the different expressed genes between subtypes. Dotplot (A) and Emap plot (B) showing the enriched pathways of the different expressed genes.

Subgroups of retinoblastoma shows distinct profiles regarding retinal markers and stemness

Since the different subtypes of retinoblastoma had distinct activation of biological pathways, we then examined whether these retinal markers maintained the same levels across subgroups. We investigated that the meta score for cone markers and ganglion markers in subtype 1 and 2, respectively. As expected, subtype 1 had a higher meta score of cone markers (Wilconxon test, p = 1.81 × 10⁻⁵) while subtype 2 had a higher level of ganglion meta score (Wilconcon test, p = 7.432 × 10⁻⁸) (Fig. 3). Interestingly, these cone and ganglion markers did not remain the same levels within the same subtypes. Subtype 1ahad the highest meta score for cone markers and subtype 1b showed the lowest meta score in all four subtype 1 subclasses. Subtype 2ahad the lowest cone meta score but subtype 2b had the similar score to subtype 1. For ganglion score, subtype 1d had the highest of all four subtype 1 subclasses and subtype 2ahad a significant higher meta score across all subgroups. Similarly, subtype 2b had the lower meta score of ganglion markers as compared with subtype 2abut the comparable meta score as compared with the subgroups in subtype 1, suggesting subtype 2b might originated from cone rather than ganglion. These results highlighted the heterogeneity of subtypes in retinoblastoma in terms of tissue of origin and classifying the subtypes of retinoblastoma by the tissue of origin need to be justified.

The cone and ganglion expression between subtypes. (A) Boxplot showing the expression of cone related genes in subtypes. (B) Box plot showing the expression of ganglion related genes in subtypes.

We then tested whether the stemness was different in subtypes of retinoblastoma. Subtype 2 had a higher stemness index than subtype 1, suggesting subtype 2 likely had more undifferentiated cells. In subtype 1, samples with low stemness index were enriched in the subtype 1b, suggesting subtype 1b likely contains more differentiated cells (Fig. 4). Then, we queried the RB and TP53 related genes expression to test the correlation of them with stemness feature across subtypes because a decrease of RBL2 and TP53 was observed in neuronal/ganglion subtype of retinoblastoma. Interestingly, we found RBL2 and TP53 relevant genes (MDM4, MDM2 and TP53) were slightly higher in subtype 2 retinoblastoma than these in subtype 1. In subtype 1, stemness index was negatively correlated with RBL2 and TP53 but showed a strong correlation of MDM4. In subtype 2, there was no strong correlation of these genes with stemness index (Fig. 4).

Stemness index between subtypes. (A) Boxplots showing the stemness index between subtypes. (B) Dot plots showing the correlation of stemness index with RBL2 and TP53 related genes between subtypes.

A large difference of immune infiltration patterns between subtypes of retinoblastomas

Since the pathway analysis of retinoblastoma subtypes indicate the immune response was the one of the key features, we next investigated the immune cell infiltration in each subtype of retinoblastoma. Through MCPcount algorithms, the percentages of different immune cell lineages across subtypes of retinoblastoma were predicted. Notably, a higher B cell lineage was only observed in subtype 1a and high proportion of neutrophils, cytotoxic lymphocytes, monocytic lineage and myeloid dendritic cells were presented in subtype 1b (Fig. 5). In subtype 2, the proportion of all immune cells lineages were relatively low as compared with subtype 1, suggesting subtype 2 retinoblastoma was likely immune “cold” tumor (Fig. 5). We also noticed that the percentage of B cell lineage are invertedly correlated with the percentage of cytotoxic lymphocytes, monocytic lineage and myeloid dendritic cells in subtype 1a and 1b, suggesting different immunity (innate immunity in subtype 1a versus adaptive immunity in subtype 1b) approaches should be taken into consideration when we design immunotherapy for retinoblastoma.

Immune cell infiltration between subtypes in retinoblastoma. (A) Heatmap showing the immune cell infiltration between subtypes. (B) Dot plots showing the correlation of B lineage proportion with other immune cells across subtypes.

For further investigating the inflammation occurred in the retinoblastoma, we then queried the cytokine profile across all subtypes of retinoblastoma. Interestingly, all subtypes of retinoblastoma shared the similar cytokine profiles (Supplementary Fig. 2).

Subtype 1b retinoblastoma had a high risk of tumor invasion and prognosis

Lastly, we examined whether the biological subtyping is indicated of the patients’ outcomes. Then, we compared the expression of retinoblastoma invasive biomarkers between subtypes. Based on our previous study, we have successfully identified a few retinoblastomas invasive signatures, providing the risk of tumor invasion. Interestingly, the retinoblastoma invasion markers (CLUL1, CNGB1, ROM1 and RDH12) were predominantly higher in subtype 1b, suggesting subtype 1b retinoblastoma prone to be more invasive as compared with other subtypes (Fig. 6). Kaplan–Meier Curve of patients’ overall survival also demonstrated that subtype 1b patients had worse outcomes while subtype 2a and 2b patients had relatively better outcomes (Log-rank test, p < 0.0001) (Supplementary Fig. 2).

Invasive feature between subtypes. Boxplots showing the expression level of the invasive gene panel of CLUL1 (A), CNGB1 (B), ROM1 (C) and RDH12 (D) across subtypes.