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