Skip to main content

Identifying new biomarkers of aggressive GroupĀ 3 and SHH medulloblastoma using 3D hydrogel models, single cell RNA sequencing and 3D OrbiSIMS imaging

Abstract

The most common malignant brain tumour in children, medulloblastoma (MB), is subdivided into four clinically relevant molecular subgroups, although targeted therapy options informed by understanding of different cellular features are lacking. Here, by comparing the most aggressive subgroup (Group 3) with the intermediate (SHH) subgroup, we identify crucial differences in tumour heterogeneity, including unique metabolism-driven subpopulations in Group 3 and matrix-producing subpopulations in SHH. To analyse tumour heterogeneity, we profiled individual tumour nodules at the cellular level in 3D MB hydrogel models, which recapitulate subgroup specific phenotypes, by single cell RNA sequencing (scRNAseq) and 3D OrbiTrap Secondary Ion Mass Spectrometry (3D OrbiSIMS) imaging. In addition to identifying known metabolites characteristic of MB, we observed intra- and internodular heterogeneity and identified subgroup-specific tumour subpopulations. We showed that extracellular matrix factors and adhesion pathways defined unique SHH subpopulations, and made up a distinct shell-like structure of sulphur-containing species, comprising a combination of small leucine-rich proteoglycans (SLRPs) including the collagen organiser lumican. In contrast, the Group 3 tumour model was characterized by multiple subpopulations with greatly enhanced oxidative phosphorylation and tricarboxylic acid (TCA) cycle activity. Extensive TCA cycle metabolite measurements revealed very high levels of succinate and fumarate with malate levels almost undetectable particularly in Group 3 tumour models. In patients, high fumarate levels (NMR spectroscopy) alongside activated stress response pathways and high Nuclear Factor Erythroid 2-Related Factor 2 (NRF2; gene expression analyses) were associated with poorer survival. Based on these findings we predicted and confirmed that NRF2 inhibition increased sensitivity to vincristine in a long-term 3D drug treatment assay of Group 3Ā MB. Thus, by combining scRNAseq and 3D OrbiSIMS in a relevant model system we were able to define MB subgroup heterogeneity at the single cell level and elucidate new druggable biomarkers for aggressive Group 3 and low-risk SHH MB.

Introduction

The most common malignant brain tumour in children, medulloblastoma (MB), is subdivided into 4 subgroups which differ in metastasis patterns and related prognoses [1,2,3]. Despite these significant subgroup differences, therapy options are still mainly non-targeted cytotoxins and novel approaches are urgently required [4]. Group 3Ā MB is frequently associated with metastases and has the worst outcome of all MB subgroups with metastasis being an independent predictor of survival in this group of patients [5,6,7]. Recently, we have shown that MB cell lines grown in 3D hydrogel models mimicking the human brain extracellular matrix can be used to recapitulate the laminar metastatic phenotype of GroupĀ 3Ā MB patients as well as the exclusively nodular phenotype of sonic hedgehog (SHH) MB patients observed in the clinic [8, 9]. Whilst Group 3 tumours are universally high risk, SHH tumours are highly heterogeneous in their outcomes; hence, an ability to delineate risk using biomarkers would be beneficial, particularly if these biomarkers are druggable. Importantly, the 3D hydrogel model recapitulated subgroup-specific adhesion, invasion, chemoresistance and extracellular matrix (ECM) remodelling phenotypes, which make it a tool well-suited to exploring subgroup-specific therapy targets.

Tumour heterogeneity is a major obstacle in cancer treatment and is an area of active research in brain tumours [10, 11]. In our 3D hydrogel models, we also observed that laminar invasion and resistance to chemotherapy were features of distinct subpopulations of GroupĀ 3 cells rather than a feature of the global cell population [8]. Consequently, we hypothesized that contrasting therapy responses of the different MB subgroups could be due to both subgroup-specific molecular signalling dependencies and heterogeneity in the underlying cell populations. Novel technologies, such as single cell RNA sequencing (scRNAseq) and label free metabolite imaging (3D OrbiSIMS) allow the analysis of cell phenotypes and functions on the single cell level and therefore visualisation of the full spectrum of brain cancer cell states [12]. This in-depth understanding is important since recurrence and metastasis can only be prevented if therapy approaches target all cancer subpopulations.

Here, we use a combination of state-of-the-art scRNAseq and mass spectrometry imaging (3D OrbiSIMS) to compare 3D MB hydrogels of the highly aggressive GroupĀ 3 to their less aggressive SHH counterpart. The 3D OrbiSIMS combines the label free imaging attributes of secondary ion mass spectrometry (SIMS) with the high mass-resolving power of an Orbitrap. Uniquely (amongst mass spectrometry imaging technologies) 3D OrbiSIMS offers the opportunity to probe organic molecules, metabolites in this instance, at high spatial resolution (2Ā Ī¼m) and with high mass resolving power and accuracy (within 5Ā ppm) [13,14,15]. The combination of the spatial resolution imaging data obtained from the 3D OrbiSIMS with the single cell gene expression subpopulation data, allowed unique characteristics of both Group 3 and SHH medulloblastoma to be identified. In the SHH sub-group we demonstrated the presence of sub-group specific ECM subpopulations and structures that correlate with good overall survival in patient data, representing a potential biomarker for risk stratification in SHH tumours. In Group 3, we identified metabolism-driven subpopulations by scRNAseq and 3D OrbiSIMS analyses where we also observed fumarate accumulation in GroupĀ 3 nodules. NMR spectroscopy confirmed an association between fumarate and poorer overall survival in patients, a feature that could be readily exploited by targeted therapy combinations.

Materials and methods

Medulloblastoma cell lines and reagents

DAOY cell line was obtained from ATCC (Manassas, USA), ONS-76 from Dr. Annette KĆ¼nkele (CharitĆ© UniversitƤtsmedizin Berlin, Germany), HD-MB03 from Dr Till Milde (DKFZ Heidelberg, Germany) and D458 from Dr. John R. Silber (University of Washington, Seattle, USA). DAOY and D458 were grown in low glucose DMEM (Thermo Fisher) with 10% FBS (HyClone, Thermo Fisher), ONS-76 and HD-MB03 were grown in low glucose RPMI 1640 (Sigma) with 10% FBS (Thermo Fisher). All cells were maintained at 37Ā Ā°C in a humidified atmosphere containing 5% CO2. During the course of this study, all cell lines were confirmed as negative for mycoplasma contamination (MycoAlert, Lonza).

HA hydrogel preparation and long-term cell culture

Hydrogels for MB 3D culture were prepared as previously described [8]. Briefly, hyaluronan cross-linked hydrogels were prepared according to the manufacturerā€™s recommendations (Hystem; BioTime Inc) with 1% Hyaluronan and 2% Extralink (PEGDA) in order to achieve a matrix stiffness of around 1.5Ā kPa and free diffusion of particles less than 75Ā kDa [16, 17]. For long-term embedding experiments, a pure gel layer (50Ā ĀµL) was pipetted and equally distributed in a 24-well plate insert (Greiner). A central cell-containing gel layer was then added (DAOY, ONS-76, D458, HD-MB03: 10,000 cells per 150Ā ĀµL gel). The ā€œgel sandwichā€ was finished by adding another pure gel layer on top (50Ā ĀµL). Each gel layer was allowed to set properly before the next layer was added. Finally, each ā€œgel sandwichā€ was covered with 100Ā Āµl of the appropriate cell medium and the lower well was filled with 300Ā Āµl of the appropriate cell medium. The cells were allowed to grow over three weeks with cell medium being replaced every 2ā€“3Ā days. All gels were monitored daily by visual inspection and photomicroscopy.

Immunohistochemistry of HA hydrogels

HA hydrogels were removed from the well plate and embedded in Histogel (Thermo Fisher) as described by Pinto et al. [18]. The Histogel blocks were incubated on ice for 10Ā min, fixed in 4% paraformaldehyde and embedded in paraffin. Sections were further processed for immunohistochemistry. Primary antibodies against Ki67 (1:400; CS9449, Cell Signalling), collagen I (1:100, ab7463), collagen VI (1:100, ab45139) were used in combination with a goat Anti-Rabbit IgG H&L secondary antibody (ab214880; all Abcam except for Ki67). The primary antibody against lumican (1:10, AF2846, R and D systems) was used in combination with a donkey Anti-Goat IgG H&L secondary antibody (1:1000, ab6885).

Measurement of glucose, lactate, glutamate, glutamine, malate and fumarate in hydrogel supernatants

Secreted metabolites were measured using the Glucose-Glo Assay, the Lactate-Glo Assay, the Glutamine/Glutamate-Glo Assay (all Promega) as well as the Fumarate Detection Assay and the Malate Assay Kit (both Abcam) following the manufacturerā€™s instructions. To avoid additional background signals, the hydrogels were cultured in phenol-red free medium supplemented with dialyzed FBS (A3382001, Thermo Fisher) for 19Ā days. Hydrogels were washed three times with PBS and fresh medium was added. Supernatant of hydrogels was collected and snap-frozen 6Ā h, 24Ā h, 48Ā h and 72Ā h after the medium change. For glucose, lactate and glutamine/glutamate detection samples were diluted at least 1:50 for measurement. Samples or sample dilutions, media controls and standard dilutions were prepared and mixed with assay buffer as described by the manufacturers. The resulting luminescence (Glucose-Glo Assay, Lactate-Glo Assay, Glutamine/Glutamate-Glo Assay) or absorbance (Fumarate Detection Assay, Malate Assay Kit) were measured on a FLUOstar Omega plate reader (BMG Labtech Ltd). For quantification, standard titration curves were prepared for glucose, lactate, glutamine, glutamate, fumarate and malate and included on each plate. Medium samples were measured and subtracted from all samples as background.

scRNA sequencing and data analysis

Four 3-week hydrogel cultures of ONS-76 and HD-MB03 cells (pools of three gels per cell lines) were transferred into gentleMACS C tubes and processed using the Tumour Dissociation kit and the GentleMACS Dissociation programme (all Miltenyi Biotec) according to the manufacturerā€™s instructions. Dead cells and cell debris were removed from the resulting suspension of single cells using MS columns and the Dead Cell Removal Kit (both Miltenyi Biotec) according to the manufacturerā€™s instructions.

Library preparation and sequencing were performed by the Deep Sequencing facility of the University of Nottingham. Briefly, single cell 3ā€™ whole transcriptome sequencing libraries were prepared from dissociated cell suspensions using the Chromium Single Cell 3ā€™ Library and Gel Bead Kit v2, the Chromium Single Cell A Chip Kit and the Chromium Multiplex Kit, 96 rxns (10X Genomics, Pleasanton, California, USA; PN-120267, PN-1000009 and PN-120262). Cell counts and viability estimates were obtained using the LUNA-II Automated Cell Counter (Logos Biosystems), Trypan Blue Stain, 0.4% and Luna Cell Counting Slides (Logos Biosystems; T13001 and L12001). Live cell counts were used to calculate cell input, rather than total cell count, as visual inspection of cell field on the LUNA II and the gating histogram, showed thatā€‰>ā€‰95% of cells looked viable and that the cell counter appeared to be counting some extracellular debris as non-viable cells. Five thousand input cells were targeted per sample, with the aim of generating sequencing libraries fromā€‰~ā€‰2,500 single cells. All steps, including GEM Generation and Barcoding, Post GEM-RT Cleanup and cDNA Amplification and Library Construction were performed according to the Chromium Single Cell 3ā€™ Reagents Kits v2 User Guide, Rev D. Variable steps of this protocol included using 12 cycles of cDNA amplification and 10 cycles of library amplification. Amplified cDNA was quantified using the Qubit Fluorometer and the Qubit dsDNA HS Assay Kit (ThermoFisher Scientific; Q32854) and fragment length profiles were assessed using the Agilent 4200 TapeStation and Agilent High Sensitivity D5000 ScreenTape Assay (Agilent; 5067ā€“5592 and 5067ā€“5593). Completed sequencing libraries were quantified using the Qubit Fluorometer and the Qubit dsDNA HS Assay Kit and fragment length distributions assessed using the Agilent 4200 TapeStation and the High Sensitivity D1000 ScreenTape Assay (Agilent; 5067ā€“5584, 5067ā€“5585).

Libraries were pooled in equimolar amounts and the final library pool was quantified using the KAPA Library Quantification Kit for Illumina Platforms (Roche; KK4824). Libraries were sequenced on the Illumina NextSeq 500 using a NextSeq 500 High Output v2.5 150 cycle kit (Illumina; 20,024,907) to generateā€‰>ā€‰50,000 raw reads per cell for each sample, using custom sequencing run parameters. Cellranger V3.1.0 (10X Genomics; Genome Reference: hg38; Annotation Source: Ensembl build 94) was used to create fastq files of all cells, to count UMIs for each gene in each cell and create the loupe files. Then Loupe Browser V4.2.0 (10X Genomics) was used to generate t-SNE plots. The Seurat package in R [19] was used to classify cells based on their cell cycle stages. This was performed for both samples following the process outlined by the Satija Lab at (https://satijalab.org/seurat/v3.2/cell_cycle_vignette.html) and as performed in [20].

To determine similarities and differences between cell subpopulations, lists of differentially expressed genes were compared using Multiple List Comparator (MolBioTools) and pairwise intersections were visualized according to the Jaccard index for the nine clusters of ONS-76 (O1-O9) and HD-MB-03 (H1-H9). Since clusters O6 and H5 did not contain significantly differently expressed genes compared to other clusters they were excluded from the further analysis. For each remaining cell cluster pair, the mean Jaccard index was calculated (mean value of Jaccard index for upregulated genes and Jaccard index for downregulated genes). Based on the mean Jaccard value, which indicates similarities in up-and down-regulated genes between two clusters, the number of network connection between these clusters was determined (mean Jaccardā€‰>ā€‰0.3: 4 connections; mean Jaccardā€‰>ā€‰0.2: 3 connections; mean Jaccardā€‰>ā€‰0.15: 2 connections; mean Jaccardā€‰>ā€‰0.1: 1 connection). The resulting network matrix was visualized in R.

Cell viability assay and Vincristine and ML385 drug treatment

Three-week old gels were treated with either 10Ā nM vincristine (VCR) (S1241; Selleckchem), 5Ā ĀµM NRF2 inhibitor ML385 (6243, Bio-Techne) or both for one week. DMSO was used as a vehicle control. Drugs were renewed 24, 72 and 144Ā h after the first dose. In order to assess the cellā€™s recovery potential all gels were washed and covered with fresh, drug-free medium one week after the final VCR treatment and monitored for a further four weeks (total experiment time: 8Ā weeks). For cell viability assays, Prestoblue reagent (Thermo Fisher) was used according to the manufacturerā€™s instructions and fluorescence intensity was measured after 40Ā min incubation time using a microplate reader (FLOUstar Omega; BMG Labtech, Ortenberg, Germany). All gels were washed four times with HBSS buffer (Thermo Fisher) before fresh medium was added to each well. Cell viability within each gel was measured weekly after drug or vehicle treatment.

Analysis of publicly available data bases

Analysis of available microarray mRNA expression and survival data was conducted using the ā€œR2: Genomics Analysis and Visualization Platformā€ (http://r2.amc.nl). The Cavalli Medulloblastoma data set was analysed [21]. In this study 763 fresh frozen primary MB samples were hybridised with the Affymetrix Human Gene 1.1 ST Array. For the comparison between MB samples and normal cerebellum the R2: MegaSampler tool (Pfister data set [22] compared to the 9 cerebellum samples from the Roth353 dataset [23]) was used.

3D OrbiSIMS imaging

Three-week hydrogel cultures of ONS-76 and HD-MB03 cells were frozen in liquid nitrogen, sectioned (30Ā Āµm thick atā€‰āˆ’ā€‰20Ā Ā°C) and arranged on indium tin oxide coated glass slides. Samples were freeze dried overnight before introduction to the mass spectrometer. 3D OrbiSIMS analysis was performed on a Hybrid SIMS ToF and OrbiTrap instrument (IONTOF GmbH). Calibration of the Orbitrap analyser was performed on a silver sample plate, using silver clusters following the method described by Passarelli et al. [13] The Bi3+ liquid metal ion gun in spectrometry mode and the ThermoFisher Tune software were employed for calibration. For the acquisition of 3D OrbiSIMS images, a 20Ā keV Ar3000+ analysis beam was used as primary ion source. Samples were analysed at room temperature across a 200ā€‰Ć—ā€‰200Ā Āµm area in positive and negative polarity. Ar3000+ primary ions were used with a target current of 0.2 nA with charge compensation performed using a low energy (20Ā eV) electron flood gun. Duty cycle was set to 38.5% and cycle time to 200Ā Āµs. Argon gas flooding was in operation to aid with charge compensation, which led to a pressure of 9.0ā€‰Ć—ā€‰10āˆ’7Ā bar in the main chamber. Optimal target potential varied for different samples. Images were acquired from selected samples at a beam size of either 2Ā Āµm or 20Ā Āµm diameter, pixel size of 2Ā Āµm and each measurement lasted one scan. Mass spectra were recorded at a resolution of 240 000 at m/z 200 in the mass range of 50ā€“750Ā m/z. The AGC target was off with the maximum injection time set at 500Ā ms. Assignments were determined by accurate mass within 3Ā ppm error of the calculated mass. Both data acquisition and the subsequent data processing were performed using SurfaceLab 7 software (IONTOF GmbH). Orbitrap data was acquired using a Thermo Fisher Orbitrap HF mass spectrometer. Visualisation of the nodules was assisted by microscopy (Olympus, CKX41 with attached Canon camera, DS126431).

HR-MAS NMR spectroscopy of fumarate levels in MB patients

Frozen diagnostic tumour tissue for 60 medulloblastoma cases were acquired from the CCLG tissue bank (CCLG 2015 BS 05). Tissue was snap frozen in liquid nitrogen shortly after surgical resection and stored at āˆ’Ā 80Ā Ā°C. Tissue had been removed prior to the patient receiving any chemo- or radiotherapy.

HR-MAS was performed at the Henry Wellcome Building Biomolecular NMR Facility at the University of Birmingham. Tissue was cut with a scalpel over dry ice to fit into a 12Ā Ī¼l or 50Ā Ī¼l zirconium rotor before being weighed. The internal standard 3-(Trimethylsilyl)propionic-2,2,3,3-d4 acid sodium salt (TMSP) (Cambridge Biosciences, Cambridge, UK) was added to the sample in a rotor dependent manner. 3Ā Ī¼l of standard was added to 12Ā Ī¼l rotors whilst 5Ā Ī¼l of standard was added to 50Ā Ī¼l rotors. D2O (Sigma Aldrich, Dorset, UK) was added to completely fill the rotor before it was fully assembled. The sample and rotor were kept cold over dry ice during preparation to prevent tissue degradation.

Spectra were acquired using a Bruker Avance spectrometer (Bruker, Coventry, UK) with a magnetic field strength of 500Ā MHz fitted with a 4Ā mm three channel HCD HRMAS z-PFG band probe. The rotor was spun at a temperature of 4Ā Ā°C to prevent metabolic activity and a frequency of 4800Ā Hz to remove spinning sidebands from the spectra. A NOESY pulse sequence was used with 2Ā s pre-saturation to suppress the water signal and a repetition time of 4Ā s. 256 or 512 averages were acquired for 50Ā Ī¼l and 12Ā Ī¼l rotors respectively. Using this protocol, an experiment with 256 averages was completed in 17Ā min.

Free Induction Decays were Fourier transformed in Topsin 2.0 (Bruker, Coventry, UK) and the resulting spectra were imported into MestReNova 9.0.1 software suite (Mestrelab Research, Spain). To ensure quality control, the spectra were visually examined for a high signal to noise ratio, a well-defined TMSP peak and clear discrimination of the choline, phosphocholine and glycerophosphocholine peaks at 3.2ā€“3.3Ā ppm. The spectra were phased, baseline corrected and the chemical shift referenced with respect to TMSP at 0Ā ppm. Fumarate was fit as a singlet peak at 6.5Ā ppm, and quantified by comparing the area of the fumarate peak to the area of the internal standard peak, taking into account the number of protons contributing to both fumarate and the standard signal.

Fumarate concentration was normalised by dividing by sample mass. The normalised fumarate concentration was compared between medulloblastoma subgroups using a non-parametric Kruskalā€“Wallis test. The cohort was stratified into samples with no quantifiable fumarate and those with detectable fumarate and subjected to Kaplanā€“Meier tests for the whole cohort and each subgroup. Hazard Ratios for the two groups were calculated using Cox Proportional Hazards models. Fishers Exact Tests were performed to test for associations between the presence of detectable fumarate and metastatic disease, large cell anaplastic histology and MYC oncogene amplification.

Generation of diagrams

The diagram in Fig.Ā 4d was generated using the online software BioRenderĀ® (https://biorender.com/).

Statistics

Data were expressed as meansā€‰Ā±ā€‰SEM. The biological number of samples and corresponding statistical test and significance level is indicated in each figure legend. All statistical analyses and plots were carried out with GraphPad Prism 9.00 (GraphPad Software Inc., La Jolla, CA, USA) unless stated otherwise. Kaplan Meier survival data analysis was performed using the R2: Genomics Analysis and Visualization Platform (http://r2.amc.nl) and statistical significance was tested using the logrank test as described in Bewick et al. [24].

Results

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

Discussion

In this study we have demonstrated, using single cell RNA Seq and high mass resolved metabolic imaging, that biological differences between medulloblastoma molecular subgroups also exist at the level of ECM structure and cellular metabolites. Importantly, these differences can be used to identify clinically-relevant targetable differences between subgroups.

Our study identified relevant extracellular matrix factors and adhesion pathways in SHH subpopulations, while extensive TCA cycle activation and in particular fumarate accumulation is a feature of Group 3 medulloblastoma. Fumarate accumulation is well described as a driver of tumour aggressiveness in renal cancers that are fumarate hydratase-deficient, due to loss of heterozygosity of the enzyme or upregulation of the mammalian target of rapamycin (mTOR) [35,36,37]. As a consequence of high intracellular fumarate levels, fumarate can react with the sulphydryl moiety of cysteine residues to produce 2-succinyl-cysteine, known as protein succination [38]. Protein succination can modify the activity of several cysteine-containing proteins, including the NRF2 inhibitor Kelch-like ECH-associated protein 1 (KEAP1) thus activating pathways involved in epithelial-to-mesenchymal-transition as well as oxidative stress responses [31,32,33,34]. Since oxidative stress response pathways are activated despite normal oxygen conditions this phenomenon has been described as a pseudo-hypoxic state [39]. Several chemotherapeutic reagents, including vincristine, are known to induce oxidative stress [40,41,42]. Therefore, the permanent upregulation of electrophilic and oxidative stress response pathways, for example mediated by the master regulator NRF2, can cause chemoresistance in several cancer types [43,44,45,46]. NRF2 inhibition or knockdown restored chemosensitivity in breast, colorectal, pancreatic and gallbladder cancer and is being increasingly discussed as a novel chemotherapy combination strategy [47,48,49]. In glioma, the most aggressive type of brain tumour in adults, overall survival as well as disease-free survival were significantly shorter in patients with high NRF2 expression and NRF2 expression also correlated with tumour grade [50].

Our study highlights that a similar NRF2-mediated pre-activation of oxidative stress response pathways could cause chemoresistance and impair treatment of the most aggressive MB subgroup. Interruption of NRF2 signalling could therefore become an interesting new target for MB therapy, especially in combination with existing chemotherapy regimens. In addition, NRF2 inhibitors could be beneficial in p53 mutated SHH cases. Although we did not detect increased levels of fumarate in the p53 mutated SHH model, NRF2 gene expression is elevated as result of the p53 mutation as described in other cancer types [51, 52]. Therefore, future studies, including in vivo studies, will be needed to investigate whether GroupĀ 3 and p53 mutated SHH MB tumour treatment can be improved by NRF2 inhibitors.

Apart from effects on NRF2, it is important to mention that protein succination causes multiple, complex intracellular changes. Therefore, broader therapy approaches that prevent fumarate accumulation inside cancer cells might be more efficient in GroupĀ 3Ā MB. In addition to new therapy options, fumarate levels also need to be tested as potential markers for patient prognosis and indicators of therapy failure. In our study, we detected fumarate levels in tumour samples of Group 3 patients at first resection before chemotherapy and showed a strong trend for poorer overall survival in these patients. Future larger studies should include the measurement of fumarate levels before and after treatment to test its potential as a prognostic marker for aggressive GroupĀ 3Ā MB.

In addition, the combination of single cell sequencing and mass spectrometry imaging techniques identified a specific ECM phenotype in SHH tumour models. Recently, we showed that laminin is a core component of the SHH shell-like structure and that laminin expression predicts good overall survival for SHH patients [8]. Here we show that collagens and sulphur-containing species are also part of this shell -like structure. SHH tumours are highly heterogeneous ranging from very high risk p53 mutated tumours to low risk desmoplastic tumours found in infants [53], in addition to standard risk tumours with a classical phenotype. The shell-like structure that we observe may be a more subtle form of nodularity, not detectable in standard histological analyses, but able to predict good outcome in tumours that would otherwise be classified as standard risk classical tumours.

Through the analysis of genomic patient data and immunohistochemical staining of SHH nodules grown in hydrogels, we have demonstrated that these sulphur-containing species are most probably a combination of SLRPs. Our data supports the proposal that lumican is one of these SLRPs, since it is located in the shell structure. Lumican is a proteoglycan that can contain multiple KS chains and plays a crucial role in collagen fibrillogenesis. It therefore may facilitate the formation and organisation of the collagen component of the SHH shell-like structure [27,28,29]. In addition, survival analysis of the same patient dataset revealed that lumican (plus fibromodulin and decorin) as well as COL1A1 and COL6A1 predict good overall survival for SHH patients in the same manner as laminin, further suggesting that they are all part of the same ECM shell-like structure as depicted in Fig.Ā 7 (Fig.Ā 7c and e, Additional file 2: Fig. S11aā€“d). In fact, during the course of this study, lumican, COL1A1 and COL6A1 have all been validated independently by Trombetta-Lima et al. as components of the medulloblastoma secreted extracellular matrix, by proteomic and immunohistochemical analysis of patient tumours [54]. Furthermore, the expression of lumican has previously been shown to have a restrictive role in the progression of several other cancer types as well as being correlated with better survival, illustrating the potential anti-tumorigenic effect of this ECM shell-like structure on SHH MB [55,56,57,58,59,60,61].

Conclusions

This study highlights that a combination of state-of-the-art single cell sequencing and mass spectrometry imaging techniques in a realistic 3D model allows the identification of sub-group-specific tumour cell populations whose unique functions can be exploited for future prognostic and therapeutic strategies. In this first example of such an approach, we show that fumarate accumulation in aggressive Group 3Ā MB tumours identifies very-high risk patients, potentially identifiable using HR-MAS NMR spectroscopy, who could benefit from a chemotherapy combination with NRF2 targeted therapy approaches. Conversely, we also demonstrate that low-risk SHH MB tumours can be identified by the formation of an ECM shell-like structure composed of laminin, type I- and VI-collagens and lumican; these patients may benefit from less harsh treatment regimens in order to ensure the best possible quality of life post-treatment. Patients without these ECM proteins could be classed as high-risk and may require more intensive treatment.

In paediatric cancer a crucial aim is reliable risk stratification including individual therapy approaches to maximise therapy effects while minimising quality of life reducing long-term side effects. Therefore, we recommend including fumarate and a panel of the SHH ECM shell proteins in future trials to test if SHH and Group 3 patients can be stratified based on these biomarkers for personalized treatments.

Availability of data and materials

All data generated or analysed during this study are included within the article (and its Additional Information files). The scRNA and RNA sequencing data have been deposited in the ArrayExpress database at EMBL-EBI (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-12591 and E-MTAB-9823 respectively. Genes differentially expressed between clusters can be found in Additional File 1.

Abbreviations

MB:

Medulloblastoma

SLRPs:

Small leucine-rich proteoglycans

TCA:

Tricarboxylic acid

scRNAseq:

Single cell RNA sequencing

NRF2:

Nuclear Factor Erythroid 2-Related Factor 2

SHH:

Sonic hedgehog

ECM:

Extracellular matrix

3D OrbiSIMS:

A Time of Flight Secondary Ion Mass Spectrometer (ToFā€“SIMS) with hybrid OrbiTrapTM functionality

KS:

Keratan sulfate

References

  1. Taylor MD, Northcott PA, Korshunov A et al (2012) Molecular subgroups of medulloblastoma: the current consensus. Acta Neuropathol 123:465ā€“472

    ArticleĀ  CASĀ  Google ScholarĀ 

  2. Northcott PA, Korshunov A, Witt H et al (2011) Medulloblastoma comprises four distinct molecular variants. J Clin Oncol 29:1408ā€“1414

    ArticleĀ  Google ScholarĀ 

  3. Gibson P, Tong Y, Robinson G et al (2010) Subtypes of medulloblastoma have distinct developmental origins. Nature 468:1095ā€“1099

    ArticleĀ  CASĀ  Google ScholarĀ 

  4. Ramaswamy V, Taylor MD (2017) Medulloblastoma: from myth to molecular. J Clin Oncol 35:2355ā€“2363

    ArticleĀ  CASĀ  Google ScholarĀ 

  5. Kool M, Korshunov A, Remke M et al (2012) Molecular subgroups of medulloblastoma: an international meta-analysis of transcriptome, genetic aberrations, and clinical data of WNT, SHH, Group 3, and Group 4 medulloblastomas. Acta Neuropathol 123:473ā€“484

    ArticleĀ  CASĀ  Google ScholarĀ 

  6. Gajjar A, Pfister SM, Taylor MD et al (2014) Molecular insights into pediatric brain tumors have the potential to transform therapy. Clin Cancer Res 20:5630ā€“5640

    ArticleĀ  CASĀ  Google ScholarĀ 

  7. Shih DJH, Northcott PA, Remke M et al (2014) Cytogenetic prognostication within medulloblastoma subgroups. J Clin Oncol 32:886ā€“896

    ArticleĀ  Google ScholarĀ 

  8. Linke F, Aldighieri M, Lourdusamy A et al (2021) 3D hydrogels reveal medulloblastoma subgroup differences and identify extracellular matrix subtypes that predict patient outcome. J Pathol 253:326ā€“338

    ArticleĀ  CASĀ  Google ScholarĀ 

  9. Zapotocky M, Mata-Mbemba D, Sumerauer D et al (2017) Differential patterns of metastatic dissemination across medulloblastoma subgroups. J Neurosurg Pediatr 21:1ā€“8

    Google ScholarĀ 

  10. Friedmann-Morvinski D (2014) Glioblastoma heterogeneity and cancer cell plasticity. Crit Rev Oncog 19:327ā€“336

    ArticleĀ  Google ScholarĀ 

  11. Hovestadt V, Smith KS, Bihannic L et al (2019) Resolving medulloblastoma cellular architecture by single-cell genomics. Nature 572:74ā€“79

    ArticleĀ  CASĀ  Google ScholarĀ 

  12. Qazi MA, Bakhshinyan D, Singh SK (2019) Deciphering brain tumor heterogeneity, one cell at a time. Nat Med 25:1474ā€“1476

    ArticleĀ  CASĀ  Google ScholarĀ 

  13. Passarelli MK, Pirkl A, Moellers R et al (2017) The 3D OrbiSIMS - Label-free metabolic imaging with subcellular lateral resolution and high mass-resolving power. Nat Methods 14:1175ā€“1183

    ArticleĀ  CASĀ  Google ScholarĀ 

  14. Meurs J, Scurr DJ, Lourdusamy A et al. (2021) Sequential 3D OrbiSIMS and LESA-MS/MS-based metabolomics for prediction of brain tumor relapse from sample-limited primary tissue archives. bioRxiv. 2020.07.15.182071

  15. Kotowska AM, Trindade GF, Mendes PM et al (2020) Protein identification by 3D OrbiSIMS to facilitate in situ imaging and depth profiling. Nat Commun 11:5832

    ArticleĀ  CASĀ  Google ScholarĀ 

  16. Hynes WF, Doty NJ, Zarembinski TI et al (2014) Micropatterning of 3D microenvironments for living biosensor applications. Biosensors 4:28ā€“44

    ArticleĀ  Google ScholarĀ 

  17. Vanderhooft JL, Alcoutlabi M, Magda JJ et al (2009) Rheological properties of cross-linked hyaluronan-gelatin hydrogels for tissue engineering. Macromol Biosci 9:20ā€“28

    ArticleĀ  CASĀ  Google ScholarĀ 

  18. Pinto MP, Jacobsen BM, Horwitz KB (2011) An immunohistochemical method to study breast cancer cell subpopulations and their growth regulation by hormones in three-dimensional cultures. Front Endocrinol 2:1ā€“6

    ArticleĀ  Google ScholarĀ 

  19. Stuart T, Butler A, Hoffman P et al (2019) Comprehensive integration of single-cell data. Cell 177:1888-1902.e21

    ArticleĀ  CASĀ  Google ScholarĀ 

  20. Nestorowa S, Hamey FK, Pijuan Sala B et al (2016) A single-cell resolution map of mouse hematopoietic stem and progenitor cell differentiation. Blood 128:e20ā€“e31

    ArticleĀ  CASĀ  Google ScholarĀ 

  21. Cavalli FMG, Remke M, Rampasek L et al (2017) Intertumoral heterogeneity within medulloblastoma subgroups. Cancer Cell 31:737-754.e6

    ArticleĀ  CASĀ  Google ScholarĀ 

  22. Northcott PA, Buchhalter I, Morrissy AS et al (2017) The whole-genome landscape of medulloblastoma subtypes. Nature 547:311ā€“317

    ArticleĀ  CASĀ  Google ScholarĀ 

  23. Roth RB, Hevezi P, Lee J et al (2006) Gene expression analyses reveal molecular relationships among 20 regions of the human CNS. Neurogenetics 7:67ā€“80

    ArticleĀ  CASĀ  Google ScholarĀ 

  24. Bewick V, Cheek L, Ball J (2004) Statistics review 12: survival analysis. Crit Care 8:389ā€“394

    ArticleĀ  Google ScholarĀ 

  25. BlĆ¼ml S, Margol AS, Sposto R et al (2016) Molecular subgroups of medulloblastoma identification using noninvasive magnetic resonance spectroscopy. Neuro Oncol 18:126ā€“131

    ArticleĀ  Google ScholarĀ 

  26. Panigrahy A, Krieger MD, Gonzalez-Gomez I et al (2006) Quantitative short echo time 1H-MR spectroscopy of untreated pediatric brain tumors: preoperative diagnosis and characterization. Am J Neuroradiol 27:560ā€“572

    CASĀ  Google ScholarĀ 

  27. Lopez SG, Bonassar LJ (2022) The role of SLRPs and large aggregating proteoglycans in collagen fibrillogenesis, extracellular matrix assembly, and mechanical function of fibrocartilage. Connect Tissue Res 63:269ā€“286

    ArticleĀ  CASĀ  Google ScholarĀ 

  28. Kalamajski S, Oldberg ƅ (2010) The role of small leucine-rich proteoglycans in collagen fibrillogenesis. Matrix Biol 29:248ā€“253

    ArticleĀ  CASĀ  Google ScholarĀ 

  29. Schaefer L, Iozzo RV (2008) Biological functions of the small leucine-rich proteoglycans: from genetics to signal transduction. J Biol Chem 283:21305ā€“21309

    ArticleĀ  CASĀ  Google ScholarĀ 

  30. Fuhler GM, Eppinga H, Peppelenbosch MP (2017) Fumarates and cancer. Trends Mol Med 23:3ā€“5

    ArticleĀ  CASĀ  Google ScholarĀ 

  31. Merkley ED, Metz TO, Smith RD et al (2014) The succinated proteome. Mass Spectrom Rev 33:98ā€“109

    ArticleĀ  CASĀ  Google ScholarĀ 

  32. Zhang DD, Hannink M (2003) Distinct cysteine residues in Keap1 are required for Keap1-dependent ubiquitination of Nrf2 and for stabilization of Nrf2 by chemopreventive agents and oxidative stress. Mol Cell Biol 23:8137ā€“8151

    ArticleĀ  CASĀ  Google ScholarĀ 

  33. Kerins MJ, Vashisht AA, Liang BX-T et al. (2017) Fumarate mediates a chronic proliferative signal in fumarate hydratase-inactivated cancer cells by increasing transcription and translation of ferritin genes. Mol Cell Biol 37

  34. Sciacovelli M, GonƧalves E, Johnson TI et al (2016) Fumarate is an epigenetic modifier that elicits epithelial-to-mesenchymal transition. Nature 537:544ā€“547

    ArticleĀ  CASĀ  Google ScholarĀ 

  35. Alam NA, Rowan AJ, Wortham NC et al (2003) Genetic and functional analyses of FH mutations in multiple cutaneous and uterine leiomyomatosis, hereditary leiomyomatosis and renal cancer, and fumarate hydratase deficiency. Hum Mol Genet 12:1241ā€“1252

    ArticleĀ  CASĀ  Google ScholarĀ 

  36. Drusian L, Nigro EA, Mannella V et al (2018) mTORC1 upregulation leads to accumulation of the oncometabolite fumarate in a mouse model of renal cell carcinoma. Cell Rep 24:1093-1104.e6

    ArticleĀ  CASĀ  Google ScholarĀ 

  37. Pollard PJ, Spencer-Dene B, Shukla D et al (2007) Targeted inactivation of Fh1 causes proliferative renal cyst development and activation of the hypoxia pathway. Cancer Cell 11:311ā€“319

    ArticleĀ  CASĀ  Google ScholarĀ 

  38. Alderson NL, Wang Y, Blatnik M et al (2006) S-(2-Succinyl)cysteine: a novel chemical modification of tissue proteins by a Krebs cycle intermediate. Arch Biochem Biophys 450:1ā€“8

    ArticleĀ  CASĀ  Google ScholarĀ 

  39. Bratslavsky G, Sudarshan S, Neckers L et al (2007) Pseudohypoxic pathways in renal cell carcinoma. Clin Cancer Res 13:4667ā€“4671

    ArticleĀ  CASĀ  Google ScholarĀ 

  40. Kaur G, Jaggi AS, Singh N (2010) Exploring the potential effect of Ocimum sanctum in vincristine-induced neuropathic pain in rats. J Brachial Plex Peripher Nerve Inj 5:3

    Google ScholarĀ 

  41. Carozzi VA, Canta A, Chiorazzi A (2015) Chemotherapy-induced peripheral neuropathy: what do we know about mechanisms? Neurosci Lett 596:90ā€“107

    ArticleĀ  CASĀ  Google ScholarĀ 

  42. Parums DV (2022) Editorial: the world health organization (WHO) fungal priority pathogens list in response to emerging fungal pathogens during the COVID-19 pandemic. Med Sci Monit 28:BR22ā€“BR29

    ArticleĀ  Google ScholarĀ 

  43. Jeddi F, Soozangar N, Sadeghi MR et al (2017) Contradictory roles of Nrf2/Keap1 signaling pathway in cancer prevention/promotion and chemoresistance. DNA Repair 54:13ā€“21

    ArticleĀ  CASĀ  Google ScholarĀ 

  44. Rushworth SA, Zaitseva L, Murray MY et al (2012) The high Nrf2 expression in human acute myeloid leukemia is driven by NF-ĪŗB and underlies its chemo-resistance. Blood 120:5188ā€“5198

    ArticleĀ  CASĀ  Google ScholarĀ 

  45. Jiang T, Chen N, Zhao F et al (2010) High levels of Nrf2 determine chemoresistance in type II endometrial cancer. Cancer Res 70:5486ā€“5496

    ArticleĀ  CASĀ  Google ScholarĀ 

  46. Vomund S, SchƤfer A, Parnham M et al (2017) Nrf2, the master regulator of anti-oxidative responses. Int J Mol Sci 18:2772

    ArticleĀ  Google ScholarĀ 

  47. Kim SK, Yang JW, Kim MR et al (2008) Increased expression of Nrf2/ARE-dependent anti-oxidant proteins in tamoxifen-resistant breast cancer cells. Free Radic Biol Med 45:537ā€“546

    ArticleĀ  CASĀ  Google ScholarĀ 

  48. Wang XJ, Li Y, Luo L et al (2014) Oxaliplatin activates the Keap1/Nrf2 antioxidant system conferring protection against the cytotoxicity of anticancer drugs. Free Radic Biol Med 70:68ā€“77

    ArticleĀ  CASĀ  Google ScholarĀ 

  49. Shibata T, Kokubu A, Gotoh M et al (2008) Genetic alteration of Keap1 confers constitutive Nrf2 activation and resistance to chemotherapy in gallbladder cancer. Gastroenterology 135:1358-1368.e4

    ArticleĀ  CASĀ  Google ScholarĀ 

  50. Zhao M, Xu H, Zhang B et al (2015) Impact of nuclear factor erythroid-derived 2ā€“like 2 and p62/sequestosome expression on prognosis of patients with gliomas. Hum Pathol 46:843ā€“849

    ArticleĀ  CASĀ  Google ScholarĀ 

  51. Tung M-C, Lin P-L, Wang Y-C et al (2015) Mutant p53 confers chemoresistance in non-small cell lung cancer by upregulating Nrf2. Oncotarget 6:41692ā€“41705

    ArticleĀ  Google ScholarĀ 

  52. Lisek K, Campaner E, Ciani Y et al (2018) Mutant p53 tunes the NRF2-dependent antioxidant response to support survival of cancer cells. Oncotarget 9:20508ā€“20523

    ArticleĀ  Google ScholarĀ 

  53. Orr BA (2020) Pathology, diagnostics, and classification of medulloblastoma. Brain Pathol 30:664ā€“678

    ArticleĀ  Google ScholarĀ 

  54. Trombetta-Lima M, Rosa-Fernandes L, Angeli CB et al (2021) Extracellular matrix proteome remodeling in human glioblastoma and medulloblastoma. J Proteome Res 20:4693ā€“4707

    ArticleĀ  CASĀ  Google ScholarĀ 

  55. Coulson-Thomas VJ, Coulson-Thomas YM, Gesteira TF et al (2013) Lumican expression, localization and antitumor activity in prostate cancer. Exp Cell Res 319:967ā€“981

    ArticleĀ  CASĀ  Google ScholarĀ 

  56. De Wit M, Belt EJT, Delis-Van Diemen PM et al (2013) Lumican and versican are associated with good outcome in stage ii and iii colon cancer. Ann Surg Oncol 20:348ā€“359

    ArticleĀ  Google ScholarĀ 

  57. Li X, Truty MA, Kang Y et al (2014) Extracellular lumican inhibits pancreatic cancer cell growth and is associated with prolonged survival after surgery. Clin Cancer Res 20:6529ā€“6540

    ArticleĀ  CASĀ  Google ScholarĀ 

  58. Troup S, Njue C, Kliewer EV et al (2003) Reduced expression of the small leucine-rich proteoglycans, lumican, and decorin is associated with poor outcome in node-negative invasive breast cancer. Clin Cancer Res 9:207ā€“214

    CASĀ  Google ScholarĀ 

  59. Li X, Kang Y, Roife D et al (2017) Prolonged exposure to extracellular lumican restrains pancreatic adenocarcinoma growth. Oncogene 36(38):5432ā€“5438

    ArticleĀ  CASĀ  Google ScholarĀ 

  60. Karamanou K, Franchi M, Vynios D et al (2020) Epithelial-to-mesenchymal transition and invadopodia markers in breast cancer: lumican a key regulator. Semin Cancer Biol 62:125ā€“133

    ArticleĀ  CASĀ  Google ScholarĀ 

  61. Li X, Lee Y, Kang Y et al (2018) Hypoxia-induced autophagy of stellate cells inhibits expression and secretion of lumican into microenvironment of pancreatic ductal adenocarcinoma. Cell Death Differ 26(2):382ā€“393

    ArticleĀ  Google ScholarĀ 

Download references

Acknowledgements

The authors would like to thank Dr Victoria Wright, Dr Nadine Holmes and Dr Fei Sang from the DeepSeq facility of the University of Nottingham for their support with the scRNA sequencing.

Funding

The research presented in this paper was supported by The Stoneygate Trust (ā€œCareer catalyst awardā€ awarded to Dr Beth Coyle in 2016), Biotechnology and Biological Sciences Research Council (BB/M008770/1) funding awarded to James Johnson, Engineering and Physical Sciences Research Council (EPSRC) EP/R025282/1 funding via the NanoPrime scheme (awarded to Dr Franziska Linke, Dr Rian L Griffiths and Dr David J Scurr), Wellcome Trust Joint Senior Investigator Award 103882 awarded to Prof Morgan Alexander, an independent Anne McLaren Fellowship (University of Nottingham) awarded to Dr Rian Griffiths and Andrea Burkillā€™s charity ā€œClaudiaā€™s causeā€ (awarded to Dr Beth Coyle). The 3D OrbiSIMS was funded by an EPSRC EP/P029868/1 strategic equipment award to Prof Morgan Alexander and Dr David J. Scurr.

Author information

Authors and Affiliations

Authors

Contributions

FL performed the in vitro hydrogel experiments with JJ contributing to some experiments. DJS and RLG helped with the design, then DJS, and RLG performed and optimised the analyses of the 3D OrbiSIMS mass spectrometry imaging data. MRA provided the 3D OrbiSIMS facility and contributed to editing of the manuscript. SK assisted with these analyses. DL and FL analysed scRNAseq data with AL assisting with bulk RNAseq data analysis. CDB and ACP performed and analysed patientā€™s patient NMR spectroscopy data. SC provided subgroup data for patients. All authors were involved in research design, manuscript writing and the final approval. FL, JJ and BC were responsible for the overall design of the research, analysis and interpretation of data. FL, JJ and BC wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Beth Coyle.

Ethics declarations

Ethics approval and consent to participate

This study has Research Ethical Committee approval (NRES East Midlands-Derby, 04/MRE04/41) and CCLG Biological Studies Committee approval.

Consent for publication

Not applicable.

Competing Interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

scRNAseq data ONS76 HD-MB 03 - ALL SIGNIFICANT GENES PER CLUSTER. This is and Excel spreadsheet that contains data relating to all the significantly altered genes with their associated Log2 fold change and p-values.

Additional file 2.

Supplementary Figures and Methods. This is a word file that contains additional figures and their legends, alongside additional methods and references relating to these.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Linke, F., Johnson, J.E.C., Kern, S. et al. Identifying new biomarkers of aggressive GroupĀ 3 and SHH medulloblastoma using 3D hydrogel models, single cell RNA sequencing and 3D OrbiSIMS imaging. acta neuropathol commun 11, 6 (2023). https://doi.org/10.1186/s40478-022-01496-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40478-022-01496-4

Keywords