scRNAseq analysis of medulloblastoma models reveals independent, metabolic subpopulations in Group 3 and SHH-specific extracellular matrix (ECM) subpopulations

In our study we used validated medulloblastoma cell lines representing the SHH and Group 3 medulloblastoma subgroups [8]. Three weeks after initiation, hydrogel models of ONS-76 (SHH) and HD-MB03 (Group 3) cells displayed heterogeneous growth and migration; to explore this behaviour at the gene expression level these models were then analysed using scRNAseq. We identified 9 different cell clusters in each model based on gene expression similarities and determined significantly up- and downregulated genes between the clusters (Additional file 1: Table S1, Fig. S1 a–b). Importantly, the SHH and Group 3 hydrogel models express different levels of subgroup-specific gene markers for the SHH pathway (SMO, PTCH1, SFRP2) and Group 3 (MYC, OTX2, GNB3) which is in line with previously performed bulk RNA sequencing data [8] (Additional file 2: Fig. S1 a–f; Additional file Methods). Based on gene expression each cell was assigned to a cell cycle phase and the cell cycle phase distribution was compared between clusters (Fig. 1 c–d, Additional file 2: Fig. S1g). Interestingly, both hydrogel models contained clusters with similar and very specific cell cycle distributions, such as clusters O7 and H3 that contained very few cells in G1. To identify general cluster similarities, all significantly up-and downregulated genes were compared between clusters, and clusters with shared up- or downregulated genes were connected in a network graph (Additional file 2: Fig. S2 a–f). In the SHH model all clusters are connected with at least one other cluster and clusters O5 and O7 display highest similarities (Additional file 2: Fig. S2c). In contrast, two Group 3 clusters, H4 and H7, are very distinct from other Group 3 clusters while cluster H9 shares similarities with 4 other clusters and could represent an intermediate state (Additional file 2: Fig. S2f). Interestingly, the comparison of cell clusters across models shows strong concordance between clusters O4, O2 and H4 as well as between O7 and H3 as was also predicted from the cell cycle analysis (Fig. 1c–e, Additional file 2: Fig. S2g–h). Notably, any similarities and differences observed between clusters were independent of their cell cycle status. KEGG pathway analysis of cluster-specific gene lists confirms some shared functional pathways between SHH and Group 3 MB models but also highlights subgroup-specific cellular functions (Fig. 1f, Additional file 2: Fig. S3). Some clusters share KEGG pathway similarities with neurological diseases such as Alzheimer’s disease (O4, O2, H4) or metabolism (H1, H3 and to a lesser degree O5, O7). Importantly, clusters H1 and H3 are the strongest metabolic clusters with the most gene expression changes related to metabolism and no comparable counterpart identified in the ONS-76 (SHH) model. Conversely, whilst enhanced metabolic clusters seem to be characteristic of the Group 3 model, clusters characterized by ECM components and adhesion pathways are exclusively found in the SHH model (O1, O3, O8; Fig. 1e–f).

Fig. 1

scRNA sequencing reveals tumour subpopulations characterized by specific gene expression patterns in SHH and Group 3 MB models. Analysis of scRNA sequencing data from ONS-76 (SHH) cells a and HD-MB03 (Group 3) cells b growing as 3D hydrogel models identified 9 different graph-based clusters for each cell line (t-SNE plot, see also Additional file 1: Table S1 for cluster characteristics). Gene expression of each cell was used to annotate a particular cell cycle phase. Relative proportion of cells in G1, S and G2/M phase for ONS-76 c and HD-MB03 d cells highlights heterogeneity between clusters (see Additional file 2: Fig. S1 for t-SNE plot of cell cycle annotation). e The network graph highlights similar clusters based on their gene expression (Intersection Tables with Jaccard Indices used to create the network graph are shown in Additional file 2: Fig. S2). Note that clusters O4 and H4 are most similar while cluster H7 is mostly unique in its gene expression compared to all other clusters. (The network graph shows cluster similarity based on the Jaccard indices of shared up-and downregulated genes with the number of connections indicating higher Jaccard values [1: 0.15 < Jaccard ≥ 0.1; 2: 0.2 < Jaccard ≥ 0.15; 3: 0.3 < Jaccard ≥ 0.2; 4: Jaccard ≥ 0.3]). f Extract of Top20 KEGG pathways of significantly up-(dark grey) and down-(light grey) regulated genes are listed for each cluster (full overview shown in Additional file 2: Fig. S3). Note the unique presence of ECM and adhesion subpopulations in ONS-76 clusters while metabolic clusters are dominant in HD-MB03

Taken together, the scRNAseq data show that although SHH and Group 3 models harbour some comparable subpopulations with shared functions, these exist alongside subgroup-specific clusters indicating heterogeneity in metabolism genes (Group 3 models) and tumour microenvironment interactions (SHH models).

Subgroup-specific features of the tumour composition are identified by metabolic imaging of 3D medulloblastoma hydrogel models

Although scRNAseq data delivers a great level of in-depth information on distinct cell clusters, their spatial location within the hydrogel is lost. In order to obtain metabolite distribution within each model the hydrogel samples were analysed by 3D OrbiSIMS to obtain mass spectrometry imaging (MSI) data. Three weeks after initiation, hydrogel models of ONS-76 (SHH) and HD-MB03 (Group 3) (Additional file 2: Fig. S4a–b) were snap frozen, sectioned, and analysed under cryogenic conditions [13]. The MSI data acquired enabled the detection of metabolite distribution both within and around the nodules and we were able to successfully resolve ions that exist at a similar m/z value (Additional file 2: Fig. S4c).

In order to identify further subgroup specific cluster differences, targeted analysis was undertaken of the acquired 3D OrbiSIMS data to specifically identify metabolites that have been previously shown to be differentially expressed in MB patients by MRI analysis [25]. Firstly, the expected high levels of choline, that typify high grade tumours including medulloblastoma [26], were observed within both SHH and Group 3 MB models (Additional file 2: Fig. S5). As observed in Group 3 patients, we also observed high levels of glucose, and significantly higher levels of lactate and taurine in the extracellular matrix surrounding Group 3 models relative to SHH[25] (Fig. 2a–d, Additional file 2: Fig. S5c). Interestingly, higher glutamate levels in SHH patients in comparison to Group 3 patients [25], were also replicated in the extracellular matrix of SHH 3D hydrogel models (Fig. 2c). Taken together, these findings indicate that the 3D hydrogel model, despite its simplicity, recapitulates some of the clinically-relevant metabolic differences between SHH and Group 3 tumours [25, 26] allowing their analysis in vitro. Inclusion of other tumour microenvironmental cell types will be required to more fully recapitulate metabolic programming in tumours.

Fig. 2

The 3D hydrogel models display subgroup-specific metabolic characteristics to those observed in MB patients. Several metabolites are known to be specifically high in MB patients [25]. 3D OrbiSIMS mass spectrometry imaging also confirms relatively high metabolite levels of glucose a. Lactate levels b are elevated in Group 3 nodules (D458, HD-MB03), while glutamate c is specifically high in SHH (DAOY, ONS-76). The Group 3-specific metabolite taurine d is also significantly higher in Group 3 nodules compared to SHH (a–d: mean ± SEM; n = 5; a, b: unpaired t-test, c: Mann Whitney test, d: unpaired t-test with Welch correction; *P < 0.05, **P < 0.01). Gel supernatants of 3 week-old SHH and Group 3 nodules were collected 6, 24, 48 and 72 h after medium change (containing 5000 µM glucose and 2000 µM glutamine) to measure consumption and secretion of glucose e, lactate f, glutamine g and glutamate h. Note that although both subgroups quickly consume glucose and glutamine, only lactate is secreted by both. Interestingly, glutamate secretion is characteristic for SHH nodules and below detection limit in Group 3 nodules (mean ± SEM; n = 3). i Gene expression for glycolysis (ALDOA, LDHB, GALM, AKR1A1, HKDC1, ENO2, ALDH2, ENO3, PGAM2, PFKM, PDHA1, PCK2, ALDH9A1), TCA cycle (SUCLG1, PDHA1, PCK2, MDH2, PC, MDH1, FH) and ECM genes (COL6A1, COL6A2, COL6A3, COL1A1, TNC) at the single cell level are displayed for the SHH (ONS 76) and Group 3 (HD MB03) model (purple: gene expressing cell; grey: non-expressing cell; percentage indicates relative proportion of total cells expressing the specified gene set). Cluster specific expression of these gene sets is shown in Additional file 2: Fig. S6

The next step, therefore, was to examine differences in energy metabolism by quantifying changes in glycolysis and glutaminolysis over time. To do this, supernatants of SHH and Group 3 models were collected over three days to measure glucose consumption and subsequent lactate secretion (Fig. 2e–f) as well as glutamine consumption and glutamate secretion (Fig. 2g–h). Starting glucose levels (5000 µM) reduced to half in all MB models after 6 h and glucose was almost completely consumed after one day (Fig. 2e). In parallel, lactate levels rose in all MB models, reaching a maximal concentration of approximately 5000 µM after 24 h and remaining stable until day three (Fig. 2f). In comparison glutamine (initial level 2000 µM) was consumed more slowly than glucose (halved after 24 h and completely consumed after 72 h in all apart from DAOY nodules, Fig. 2g). Importantly, secreted glutamate levels rose steadily over time only in the SHH models while remaining low or below detection limit in Group 3 models, thus further confirming subgroup differences in metabolic pathways between SHH and Group 3 MB (Fig. 2h).

The upregulation of metabolic pathways involved in glycolysis and the tricarboxylic acid (TCA) cycle was also strongly observed at the RNA level in Group 3 MB (Fig. 2i, Additional file 2: Fig. S6). In addition, many other pathways that can fuel the TCA cycle, such as fatty acid metabolism, are among the overall upregulated KEGG pathways in Group 3 models (Additional file 2: Fig. S6a). At the single cell level, 66% of HD-MB03 (Group 3) cells expressed genes involved in glycolysis in comparison to 28% of ONS-76 (SHH) cells (Fig. 2i, Additional file 2: Fig. S6b). For gene sets involved in the TCA cycle, metabolic differences became even more striking (52% of HD-MB03 cells expressed TCA genes relative to 18% of ONS-76 cells) (Fig. 2i, Additional file 2: Fig. S6c). In summary, there are clear metabolic differences found exclusively in Group 3 models that may represent therapeutic targets while on the other hand SHH tumours are characterized by a specific matrix components that are worthy of further analysis.

Small leucine rich proteoglycans and collagens are biomarkers of SHH MB

3D OrbiSIMS was also used to identify differences between SHH and Group 3 by unbiased mass analysis, revealing differences in the spatial distribution of several sulphur-containing compounds (Fig. 3a). Sulphur-containing species such as SO4− are specifically located at the outer shell of SHH nodules; conversely, these are more diffuse in Group 3 nodules, indicating a different functional role between MB subtypes (Fig. 3a). The location of these sulphur-containing species overlaps with collagens that are also specifically found at the outer shell of SHH nodules, validated by staining (Fig. 3b). Gene expression of these collagens, including type I- and type VI-collagens, is also significantly higher in SHH patients than Group 3 (Fig. 3c, Additional file 2: Fig. S7) and is associated with clusters O1 and O3, which have been identified as SHH-specific ECM receptor interaction and focal adhesion clusters using scRNAseq analysis (Fig. 1e and f). Interestingly, differential gene expression analysis of a genomic patient dataset revealed a strong correlation between the expression of both type I- and type VI-collagens (COL1A1 and COL6A1) and several small leucine rich proteoglycans (SLRPs), including lumican (LUM) and fibromodulin (FMOD) (Fig. 4a and b), as well as several of the other collagens that are more highly expressed in SHH patients than Group 3 (COL1A2, COL6A3, and COL3A1) (Fig. 3c, Fig. s8). Similarly, scRNAseq analysis revealed that within SHH nodules 70% of LUM positive cells co-expressed the genes that defined the ECM cluster (COL1A1, COL6A1 and TNC; Fig. 4c). Lumican and fibromodulin can both contain keratan sulfate (KS) chains and have overlapping but also unique roles in collagen fibrillogenesis, suggesting their KS may be the sulphated component of the shell-like structure that forms around SHH nodules as identified by 3D OrbiSIMS (Fig. 3a) [27,28,29]. In support of this, immunohistochemical staining of lumican in the ONS-76 and HD-MB03 nodules demonstrated that lumican does indeed form a shell-like structure around SHH nodules, similar to that of type I- and type VI- collagens. In contrast, lumican was absent from Group 3 nodules (Fig. 4d).

Fig. 3

SHH nodules are uniquely characterized by deposits of sulphur-containing species and a collagen-based outer shell. a 3D OrbiSIMS mass spectrometry imaging identified several sulphur-containing species within the outer shell of ONS 76 nodules that are non-specifically distributed in HD-MB03 nodules. b Similar to the localization of sulphur species, collagens (Col I and Col VI) are specifically located at the outer surface of SHH (DAOY, ONS-76) nodules compared to nonspecific or no expression in Group 3 nodules (D458, HD-MB03). c Expression of COL6A1, COL6A2, COL6A3, COL1A1, COL1A2 and COL3A1 is significantly higher in SHH patients compared to Group 3 patients in the Cavalli data set (Cavalli et al. [21], SHH: n = 223 and Group 3: n = 144; Kruskal–Wallis test with Dunn’s post hoc test; ***P < 0.001)

Fig. 4

Deposits of sulphur-containing species in SHH nodules are likely to be a combination of small leucine rich proteoglycans including lumican. a and b Small leucine rich proteoglycans (SLRPs) and other collagens are correlated with high COL1A1 and COL6A1 expression in SHH patients in the Cavalli data set (Cavalli et al. [21], SHH: n = 223, log-rank test used to establish ‘high’ and ‘low’ expression; one-way ANOVA, grey dots represent samples where p ≥ 0.05 and navy dots where p < 0.05) c scRNAseq analysis revealed that LUM (a SLRP) is expressed by 25% of cells within SHH nodules and that 70% of these cells also express ECM genes of interest (purple: gene expressing cell; grey: non-expressing cell; percentage indicates relative proportion of total cells expressing the specified gene and percentage ‘shared’ indicates the proportion of LUM expressing cells that also express the ECM genes COL6A1, COL6A2, COL6A3, COL1A1 and TNC. d Immunohistochemical staining of ONS76 and HD-MB03 hydrogel nodules for lumican reveals that lumican is specifically located on the outer edge of SHH nodules (similar to type-I and type-VI collagens) and is not present in Group 3 nodules

Hence, our scRNAseq, 3D OrbiSIMS and IHC data all point to the conclusion that SHH nodules are specifically characterized by a distinct ECM shell-like structure of which type I- and type VI- collagens and SLRPs (lumican) are components.

Fumarate accumulation during the TCA cycle as a novel target of Group 3 medulloblastoma

The upregulation of metabolic pathways involved in glycolysis and the TCA cycle, observed in Group 3 MB by scRNAseq analysis, then led us to interrogate all TCA cell cycle metabolites in the 3D OrbiSIMS mass spectrometry imaging data (Fig. 5a, b, Additional file 2: Fig. S9a) with the aim of identifying key players in Group 3 metabolism. While detected levels of pyruvate, isocitrate and glutamate were very low to low, intermediate levels of α-ketoglutarate, high levels of succinate and very high levels of fumarate were indicated (Fig. 5a, b). Importantly, malate appears to be absent in the model. This is believed to be an accurate representation of the chemistry of these systems (as opposed to a low ionizability of malate) since the 3D OrbiSIMS data acquired in a series of control experiments demonstrated that malate could be detected at concentrations as low as 10 nM in control gels (Additional file 2: Fig. S9b–e). Using an enzymatic detection assay approach, we then confirmed that increasing concentrations of fumarate are secreted by the MB hydrogel models with at least 5-times higher concentrations in the Group 3 model while malate was again completely absent within the assay’s detection limit (Fig. 5c). This indicates an accumulation of fumarate as an onco-metabolite. In other cancers fumarate accumulation has been shown to initiate protein succination in the cytosol [30]. Targets known to be upregulated as a result of enhanced succination include oxidative stress response and epithelial-to-mesenchymal transition markers [31,32,33,34]. Notably, several of these known targets (NRF2/NFE2L2, HIF1A, FOXM1, PDK1, ZEB1 and TWIST1) are significantly upregulated in MB patients compared to normal cerebellum (Additional file 2: Fig. S10 a–f).

Fig. 5

Fumarate accumulation is predominantly observed in Group 3 models. a 3D OrbiSIMS analysis of TCA metabolites after 16 h reveal low levels of malate (m/z 133.01), pyruvate (m/z 88.02), isocitrate (m/z 174.02) and glutamate (m/z 146.05), intermediate levels of α-ketoglutarate (m/z 145.01), high levels of succinate (m/z 117.02) and very high levels of fumarate (m/z 115.00). b Normalised intensity measurements confirm the accumulation of succinate and fumarate and the concurrent lack of malate after 16 h in ONS-76 and HD-MB03 cells. c Gel supernatants of 3 week-old SHH (DAOY, ONS 76) and Group 3 (D458, HD MB03) nodules were collected 6, 24, 48 and 72 h after medium change to quantify the secretion levels of fumarate and malate. Note the increasingly higher fumarate levels in the group 3 nodules compared to SHH nodules over time. Malate levels were below the assay detection limit at all time points (mean; n = 3)

These data suggest that the accumulation of fumarate and subsequent upregulation of target genes, including succination-inducible nuclear factor erythroid 2-related factor 2 (NRF2/NFE2L2), might be important cancer-driving mechanisms in Group 3 MB pathogenesis. This led us to hypothesise that pre-activation of oxidative stress response pathways in Group 3 MB could contribute to chemoresistance. To test this hypothesis, three-week hydrogel culture models were treated with 4 repeated doses of vincristine alone or in combination with an NRF2 inhibitor (ML385) over one week and cell viability was monitored over the following four weeks (Fig. 6a). We have previously shown that, in common with patients, Group 3 MB models are more resistant to chemotherapy than SHH models [8]. Here again, while vincristine alone was not sufficient to significantly decrease cell viability of the two Group 3 models, combination with the NRF2 inhibitor significantly enhanced the chemotherapy effect (Fig. 6b, c). In contrast, cell viability of the p53 wildtype ONS-76 SHH model was significantly decreased by vincristine. The NRF2 inhibitor also showed a modest decrease in growth (Fig. 6d). In contrast, the NRF2 inhibitor exhibited a comparable but not additional effect to vincristine on cell viability in the p53 mutant DAOY SHH model (Fig. 6e). NRF2 gene expression is upregulated in p53 mutated patients and cell lines independent of fumarate, thus explaining the effect of the NRF2 inhibitor on the DAOY cells (Additional file 2: Fig. S10g–h). In further support of this finding, whilst statistical significance is not quite reached, gene expression levels of NRF2 appear to better distinguish between patients with good and poor prognosis in Group 3 MB but not in SHH patients (p = 0.078 versus p = 0.924, Fig. 6f, g), although NRF2 gene expression is highest in the SHH alpha subtype which consists of mostly p53 mutated patients (Additional file 2: Fig. S10g–h; [21]).

Fig. 6

A combination of chemotherapy and NRF2 inhibition significantly improves long-term treatment of Group 3 models. a A scheme illustrates the course of the long-term drug treatment assay as previously established [8]. After 3 weeks of growth inside the HA hydrogels the SHH cell lines ONS-76 b, DAOY c and the Group 3 cell lines HD-MB-03 d, D458 e were treated four times with either 10 nM vincristine (green), 5 µM ML385 (NRF2 inhibitor; blue) or a combination (red) or vehicle (black) during one week and cell viability was monitored for the following four weeks in the absence of drug/vehicle present anymore. In the p53 wt SHH cell line ONS-76 only the chemotherapeutic reagent vincristine significantly decreases cell viability, while the NRF2 inhibitor (ML385) alone and in combination with vincristine are also effective in the p53mut SHH cell line DAOY. In contrast, in both Group 3 models the combination of vincristine and NRF2 inhibitor significantly reduced cell viability. (mean ± SEM, n = 3; Two-way ANOVA and Dunnett’s post hoc test, *P < 0.05, **P < 0.01 and ***P < 0.001 all relative to DMSO according to colour code). Analysis of the biggest publicly available MB data base [21] shows that NRF2 (gene name: NFE2L2) gene expression does not predict survival of SHH (f; logrank test, p = 0.924) patients, but is associated with worse survival in Group 3 patients (g; logrank test, p = 0.078). Note the exclusive effect of NRF2 expression in Group 3 patients

Fumarate and lumican are biomarkers that can predict outcome in aggressive Group 3 and SHH MB

Having shown that accumulation of the TCA cycle metabolite fumarate is associated with worse survival of Group 3 MB patients, and that blocking NRF2, a fumarate-mediated oxidative stress pathway member, can increase the efficacy of chemotherapeutic agents in our model, we then sought to corroborate these findings with patient data. HR-MAS NMR spectroscopy was used to detect fumarate levels in 29 SHH and Group 3 MB patient samples and the cohort was stratified into samples with no quantifiable fumarate and those with detectable fumarate levels. The presence of detectable fumarate was neither associated with metastatic disease (p = 0.298), large cell anaplastic histology (p = 1) nor MYC oncogene amplification (p = 1). When we analysed overall survival (OS) data, however, we found a strong trend between fumarate detection and poor outcome in Group 3 MB patients compared to no association in SHH patients (p = 0.1793 versus p = 0.661, Fig. 7a, b).

Fig. 7

HR-MAS NMR spectroscopy detected fumarate levels group 3 and high LUM gene expression levels in SHH MB patients predict better overall survival. HR-MAS NMR spectroscopy analysis of fumarate in SHH (a, n = 13) and Group 3 (b, n = 16) patients indicates no survival difference for SHH patients (no fumarate detected n = 8; fumarate detected n = 5) but a strong trend for worse overall survival in Group 3 patients (no fumarate detected n = 12; fumarate detected n = 4) with detectable fumarate levels. c Genomic analysis of LUM expression in SHH (c; logrank test, p < 0.001) and Group 3 (d; logrank test, p = 0.056) patients identified that high LUM expression correlates with better overall survival than low LUM expression, whereas the opposite trend is true in Group 3 patients. e A graphical abstract of this study illustrating that by using by using a combination of state-of-the-art techniques (scRNAseq and 3D OrbiSIMS), in a realistic 3D model, we have identified sub-group-specific tumour phenotypes in SHH and Group 3 medulloblastoma. In the SHH sub-group we have demonstrated the formation of an ECM shell-like structure composed of laminin, type I- and VI-collagens and lumican that can be used to identify low-risk SHH tumours. This could facilitate more accurate risk stratification of SHH medulloblastoma. In the Group 3 sub-group we have shown that fumarate accumulation identifies very-high risk patients that could benefit from a chemotherapy combination with NRF2 targeted therapy approaches

Based on our findings that SHH nodules form a sub-group specific ECM shell-like structure, we then sought to establish whether these newly identified ECM shell components could predict survival of SHH patients. Analysis of a large genomic patient dataset revealed that the SLRPs lumican, fibromodulin and decorin as well as COL1A1 and COL6A1 all correlate with good overall survival in SHH MB patients; suggesting that the presence of this ECM shell in SHH tumours may be a reliable biomarker of good outcome (Fig. 7c, Additional file 2: Fig. S11a and b) [21]. Notably, in Group 3 patients lumican expression is associated with a much poorer survival, highlighting that capturing lumican as part of a complex with other SLRPs and COL1A1 and COL6A1 in an ECM ‘shell’is key to this effect on outcome in the SHH sub-group (Fig. 7d).