Skip to main content

Metabolic and transcriptomic profiles of glioblastoma invasion revealed by comparisons between patients and corresponding orthotopic xenografts in mice

Abstract

The invasive behavior of glioblastoma, the most aggressive primary brain tumor, is considered highly relevant for tumor recurrence. However, the invasion zone is difficult to visualize by Magnetic Resonance Imaging (MRI) and is protected by the blood brain barrier, posing a particular challenge for treatment. We report biological features of invasive growth accompanying tumor progression and invasion based on associated metabolic and transcriptomic changes observed in patient derived orthotopic xenografts (PDOX) in the mouse and the corresponding patients’ tumors. The evolution of metabolic changes, followed in vivo longitudinally by 1H Magnetic Resonance Spectroscopy (1H MRS) at ultra-high field, reflected growth and the invasive properties of the human glioblastoma transplanted into the brains of mice (PDOX). Comparison of MRS derived metabolite signatures, reflecting temporal changes of tumor development and invasion in PDOX, revealed high similarity to spatial metabolite signatures of combined multi-voxel MRS analyses sampled across different areas of the patients’ tumors. Pathway analyses of the transcriptome associated with the metabolite profiles of the PDOX, identified molecular signatures of invasion, comprising extracellular matrix degradation and reorganization, growth factor binding, and vascular remodeling. Specific analysis of expression signatures from the invaded mouse brain, revealed extent of invasion dependent induction of immune response, recapitulating respective signatures observed in glioblastoma. Integrating metabolic profiles and gene expression of highly invasive PDOX provided insights into progression and invasion associated mechanisms of extracellular matrix remodeling that is essential for cell–cell communication and regulation of cellular processes. Structural changes and biochemical properties of the extracellular matrix are of importance for the biological behavior of tumors and may be druggable. Ultra-high field MRS reveals to be suitable for in vivo monitoring of progression in the non-enhancing infiltration zone of glioblastoma.

Introduction

Management of patients suffering from glioblastoma (GBM, WHO grade IV), the most common and most malignant primary brain tumor in adults remains a challenge. Even with the latest treatment options, the median survival remains below two years [52]. This poor outcome has been attributed to the hallmarks of GBM that comprise a plethora of altered pathways, intra-tumoral (epi)genetic and metabolic heterogeneity, the interaction with the tumor microenvironment, and characteristic properties of tumor stem like cells [6, 13, 39]. These confer properties relevant for the invasive behavior and cell plasticity that render GBM particularly resistant to treatments [4, 33]. The diffuse infiltration of gliomas into the surrounding brain precludes complete resection and gives rise to tumor recurrence even in macroscopically fully resected tumors.

Importantly, the extent of infiltration is not visible using conventional T1 and T2-weighted Magnetic Resonance Imaging (MRI). Hence, it is difficult to target treatment to this “invisible” part that has a distinct microenvironment and is protected by an intact blood brain barrier (BBB).

Patient derived orthotopic xenograft (PDOX) models, implanting glioma stem-like cells or glioma derived spheres into the brain of immune-compromised mice, are frequently used to study tumor relevant molecular mechanisms and to evaluate novel therapies. The implanted cells give rise to tumors that recapitulate many relevant characteristics of the parental tumors, including diffuse infiltration [17, 24, 51, 59]. Glioma cells are metabolically reprogrammed in order to support proliferation, migration and invasion, and to survive in a hostile environment lacking oxygen and nutrients [5, 26, 35, 36]. In vivo detection of altered metabolic profiles in GBM, which can be achieved with ultra-high field 1H MR spectroscopy (1H MRS) yielding highly resolved spectra provides an opportunity to gain insights into the mechanisms of tumor progression [8, 22, 27]. This allows quantification of molecules involved in energy metabolism, myelination, neurotransmission, antioxidation, and osmoregulation, as determined in rodent xenograft models of glioma [23, 29, 37, 47]. The non-invasiveness of the technique allows longitudinal measurements of the metabolism during glioma development in the natural environment, while simultaneously probing the metabolism of the surrounding “non-tumoral” brain [29].

In this study we characterize molecular features associated with highly invasive growth of GBM that may serve as markers of early relapse or response to therapy, and may allow to shed light on underlying biological processes of tumor-host interaction. To this end, we investigated GBM in patients, paired with follow-up of their respective PDOX upon transplantation into the brains of immunocompromised mice. Analyses and comparisons included radiologic evaluation using ultra-high field MRI and 1H MRS of the patients (7 Tesla) and longitudinal follow-up of the respective developing PDOX (14.1 Tesla). Subsequently, we associated the metabolic profiles obtained by 1H MRS with the corresponding transcriptomes to gain insights into molecular mechanisms of tumor-host interaction of GBM invasion. Integrating metabolite profiles with the associated transcriptome allowed novel insights into the biological interactions, while distinguishing contributions originating from the human tumors and the invaded mouse brain, respectively.

Methods

Patient selection and 1H MRS of patients

Patients planned for surgery of a suspected GBM at the Lausanne University Hospital (CHUV) were enrolled (clinicaltrials.gov, NCT02904525) with written informed consent. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the local ethics committee CER-VD (F-25/99, 268/14).

Patients underwent 1H MRS/I in a 7 Tesla/68 cm MR-scanner (Siemens Medical Solutions, Erlangen, Germany). B0 field homogeneity was optimized using FASTMAP [18]. A 32-channel receive coil (NOVA Medical Inc., MA) with a single channel volume transmit coil or a 1H two loops surface coil were used depending on the location of the gliomas. 3D T1-weighted MR images acquired using MP2RAGE (TE/TR = 3.37/5000 ms, TI1/TI2 = 700/2200 ms, slice thickness = 1 mm, FOV = 176 × 256 mm2, matrix size = 176 × 256) [34] were used to position the volume of interest (VOI) for MRS measurements. Single voxel 1H MR spectra (SVS) were obtained using the semi-adiabatic SPECIAL localization sequence [37, 64] with TE = 16 ms, TR = 5.5–8.5 s (depends on the SAR restriction), and number of averages = 48–96. 2D multi-voxel 1H magnetic resonance spectroscopic imaging (MRSI) was measured by a sEmi-Adiabatic Spin-Echo MRSI sequence (EASE) [63] for P1, P12 and P14, using the following parameters: TE = 16 ms, TR = 4.3 s, NA = 1, FOV = 200 × 200mm2, VOI = 60 × 60mm2, slice thickness = 15 mm, matrix = 16 × 16, elliptical k-space sampling. Water and lipid suppression techniques were applied prior to the localization using VAPOR with OVS according to consensus recommendations [56].

1H MRS spectra were analyzed by LCModel [40] using a basis set with simulated metabolite spectra and an experimentally measured macromolecule baseline [48]. LCModel simulated macromolecule and lipid components were used during the analysis to allow fitting for potential lipid and MM resonances that arose from tumors. Due to time restriction in patient scans, additional water acquisition for normalization purposes could not always be performed. Metabolite ratios to total creatine were calculated for both SVS and MRSI measurement. Metabolites that were quantified with Cramer-Rao lower bound (CRLB) > 50% were excluded in the analysis.

Orthotopic mouse glioma model

Tumor tissue obtained at surgery was split in two, one part was frozen for subsequent analyses and the second part was dissociated into single cells and re-suspended in stem cell media (DMEM/F12 supplemented with B27 and growth factors) as described [49]. The next day 105 cells were transplanted stereotactically in a volume of 5 µl (Hanks’ balanced salt solution, HBBS, with phenol red, no calcium, no magnesium; Thermo Fisher Scientific) into the striatum (coordinates: bregma 0.5 mm anterior, 2 mm lateral and 3 mm ventral)[58] of male immunocompromised mice (n = 3–6/patient; age, 6–8 weeks; NOD-SCID gamma knock-out mice, NOD.Cg-Prkdcscid II2rgtm1Wjl/SzJ, bred in-house) using a micro pump (injection rate 5 µl/1 min, Stoelting). For the procedure, anesthetized mice were placed into a stereotactic frame, and fixed with a mouth piece (Stoelting) as previously described [19]. Mice were checked daily and sacrificed at first signs of neurologic symptoms (lethargy, ataxia and seizures) or body weight loss (> 10%). All animal procedures were performed under anesthesia/ analgesia, and protocols were approved by the concerned Swiss authorities (VD-1181_6; VD-2777).

In vivo 1H MRS of orthotopic xenografts in the mouse

1H MRS experiments were carried out in a 14.1 Tesla animal scanner with a 28-cm horizontal bore (Agilent Technologies, Palo Alto, CA, USA) using the SPECIAL sequence (TR = 4 s, TE = 2.8 ms, 160 or 240 scans) [37] as previously described [29]. An axial T2-weighted image (fast spin-echo sequence, TR/TE = 5000/13 ms, FOV = 18 × 18 mm, slice thickness 0.6 mm, 6 averages) was acquired before the 1H-MRS for voxel positioning (VOI, 2 × 2x2 mm3), centered in the striatum of the injected, and symmetrically, in the contralateral hemisphere. SVS was preferred over MVS due to the location of the tumors and the shorter acquisition time, thus allowing repeated measurement in the same animal. Modifications in the BBB were assessed in selected cases with T1-weighted coronal fast spin-echo images with Gadolinium contrast (Gadovist 1.0, Bayer, Leverkusen, Germany).

The first scan was performed 6 weeks after injection or at onset of symptoms (earliest week 4), and was repeated every 1–2 weeks. Spectra were quantified using the LCModel [29, 40] using a simulated basis set of brain metabolites combined with an experimental spectrum of macromolecules (Mac) acquired in a healthy subject and simulated lipids (at 0.9, 1.3 and 2.0 ppm) when present according expert recommendations [12]. Data shown are selected for accurate quantification, following the criterion CRLB < 40%. Metabolites were normalized to tCr unless stated otherwise. Metabolite concentrations using water as internal reference were also computed.

Tissue processing and immunohistochemistry

At final sacrifice, the fresh mouse brain was cut using a brain mold. The central coronal brain slice (~ 5 mm) encompassing the injection-site was frozen in OCT (Tissue-Tek), and the rest was formalin fixed and paraffin embedded for histology. Serial frozen sections were cut on a cryostat from the central slice of the mouse brains. The first and last sections were used to determine tumor location and estimate tumor cell content, based on hematoxylin and eosin (H&E) staining and nuclear immune-positivity for human-specific NCL (ab13541, Abcam, Cambridge, UK) as previously described [29]. Xenografts were macro-dissected, guided by H&E and NCL staining, and collected from the injected and contralateral side separately. In absence of apparent spread to the contralateral side, the tissue of the symmetric area on the contralateral side was collected. Immunohistochemistry was performed for EGFR (E30), Ki67 (MIB1), GFAP (G-A-5), and TP53 (DO-7) (platform Ventana, Roche).

DNA/RNA isolation and PCR

RNA/DNA was isolated from the macro-dissected xenografts (injected/contralateral side, separately) and human GBM (AllPrep DNA/RNA Mini Kit, Qiagen). The ratio of human/mouse cells in the xenografts was estimated by species specific PCR (DNA) [2]. Library preparation and RNA-sequencing was performed at the Lausanne Genomic Technologies Facility (LGTF, University of Lausanne; TruSeq Stranded Total RNA Library Kit; Illumina HiSeq 2500). The samples were barcoded and randomized between lanes for sequencing. For samples with a tumor/mouse cell ratio of ~ 50%, we aimed at 60 × 106 reads, and for all other samples, including the human GBM 30 × 106 reads.

Data preparation and analysis of RNA sequencing

Preprocessing of RNA sequencing (RNAseq) data was performed following the standard pipeline and recommendations from bcbio-nextgen (version 1.0.4, http://bcbio-nextgen.readthedocs.org/en/latest/). Reads were aligned to the human (GRCh37) and mouse reference genome assembly (mm10) by hisat2 aligner (version 2.1.0), and classified into three classes (mouse, human, ambiguous) by the Disambiguate algorithm [1]. Transcripts with low read counts, or classified as ambiguous were removed. The gene expression data were summarized by trimmed means of M-values (TMM) of normalized counts (R package edgeR) [44, 45], including log-transformation and using read counts and full library size. The Variant calling analysis was performed with the software VarDict [30] and the genomic variant annotations were obtained by SNPeff [9]. Single nucleotide variants (SNVs) with low quality estimation and silent and synonymous mutations were excluded. The genes listed in the COSMIC database as mutated in glioma and glioblastoma [54] were used for our analysis. In addition, the SNVs identified were compared with the mutation database from The Cancer Genome Atlas (TCGA) for Glioma and Glioblastoma (R package RTCGA.mutations). For personal privacy reasons the RNA-sequencing raw data will be made available upon request.

Data analyses and variable selection

The metabolite profiles and expression data were examined by principal component analysis (PCA) and heatmap representation based on Euclidean distance and Ward’s algorithm for clustering. Missing values were imputed by regularized iterative PCA algorithm [25]. The global differences between groups were tested by a Monte-Carlo test (permutation) on the between-groups inertia percentage [46].

The MVS in patients and longitudinal metabolic profiles in mice were analyzed by the STATIS prcedure [31, 57] that allows simultaneous analysis of different data arrays, matched by common columns (same variables) based on principal component analysis (PCA). Briefly, in using inertia operators and RV-coefficients [52], STATIS compares the “images” (interstructure) of each dataset, to find a consensus (compromise) and simultaneous representations of each dataset on the same factorial plane (intrastructure).

A flowchart details the strategies for gene selection and procedures used for data integration with metabolites and gene ontology (Additional file 1: Fig. S1). Gene expression was analyzed by sparse PCA (SPCA) [50] for gene selection, using singular value decomposition and lasso regularization (Additional file 1: Fig. S1A-B). The gene signature was consolidated by bootstrap (50 repetitions). At each iteration, two components were retained to describe data organization, and 50 variables (genes) were kept in each loading vector [50]. Thus, the most frequently selected genes were retained (cut-off ≥ 0.1).

The correlation structure between metabolite profiles and gene expression, was investigated by sparse Partial Least Squares analysis (SPLS) [3, 32] after removing unwanted effects of tumor origin by within-group PCA (WCA) [43] (Additional file 1: Fig. S1C). The SPLS approach combines both integration and additional variable selection (lasso regularization) simultaneously on two data sets in a one-step strategy. The selected gene-set was consolidated by bootstrap (50 repetitions). At each iteration, the association between metabolite profiles and gene expression was summarized by two components, for which all metabolites were retained and 50 variables (genes) were kept in each loading vector. Finally, the most frequently selected genes were retained (cut-off ≥ 0.1). The Coinertia analysis [14], a multivariate method for coupling two tables, summarizes the correlation structure between metabolite profiles (1H MRS) and expression of the selected genes (R packages ade4 and mixOmics). The common correlation structures between gene expression and metabolite profiles were investigated by permutation test and reported as RV-coefficient (vectorial correlation coefficient) [21].

Gene set enrichment analysis (GSEA) was performed with the molecular signature database (MSigDB v7.0, updated August 2019, all 8 collections) [53] using hypergeometric tests (R packages msigdbr and ClusterProfiler). Gene-sets with Bonferroni adjusted P-values ≤ 0.1 were considered significant. The conversion of mouse genes into the corresponding human homologs was performed with R package biomaRt [15]. All analyses and graphical representation were performed with R version 3.6.1 (URL http://www.R-project.org) [42] and the R packages survival [55], missMDA and ade4 [7].

Results

To determine metabolite profiles across distinct parts of unresected GBM, patients with suspected GBM were enrolled into the study and underwent scans at 7 Tesla prior to planned surgery. Highly resolved spectra were obtained using single-voxel and multi-voxel MRS (MVS) (Table 1).

Table 1 Baseline information of patients and PDOX

Depending on the patient’s clinical condition and tumor location, not all analyses were possible. For nine patients enrolled, fresh tumor tissue became available at surgery for orthotopic transplantation into mice. Six of the patients’ tumors were subsequently diagnosed as GBM, and three as astrocytoma grade III, whereof one was diagnosed as IDHmt, 1p/19q non-codeleted. The patient baseline description is summarized in Table 1.

Metabolic profiles of tumor development and invasion in PDOX

After stereotactic transplantation of the patient derived GBM cells into the striatum of 4 to 6 immunocompromised mice, the tumor development was followed longitudinally by MRI and 1H-MRS at ultra-high magnetic field (14.1 T). First scans were performed 4 to 6 weeks after transplantation, and measurements were repeated thereafter at one to two week intervals, placing the voxel in the presumed injection site and symmetrically in the contralateral side of the brain. The acquired highly resolved spectra, allowed measurements of 21 metabolites (Fig. 1a; Additional File 2). Temporal changes of the metabolite profiles monitored on the injected and contralateral sides, were in general more sensitive for early detection of tumor growth than standard clinical MRI acquisitions performed in parallel. Eventually tumor development was revealed for all patients’ samples, except for the IDH-mutant astrocytoma grade III (patient 10, P10). The growth kinetics of the tumors and the spread to the contralateral side was reflected in the evolution of the metabolite profiles over time (Fig. 1c–d). In contrast, the metabolite profiles of the mice transplanted with cells from P10 remained at base line, over the observation time of up to 250 days, in line with the lack of tumor development (Fig. 1d). The metabolite profiles from the injected and contralateral side are visualized for each mouse and measurement in function of time from injection, stratified by patient (Fig. 1d). These metabolite profiles are represented by their coordinates on the first axis of the STATIS compromise analysis (similar to PCA) including all measurements. The corresponding representation of the metabolites on the first axis is depicted in Fig. 1e. Growth and invasion were associated with decreasing neuronal metabolism (N-acetyl aspartate, NAA; glutamate, Glu; and gamma aminobutyric acid, GABA) and increase in metabolite markers specific for high cellular turnover such as choline compounds (total choline, tCho; glycerophosphocholine, GPC) and myo-inositol (Ins) as visualized by their coordinates on the first axis of STATIS (Fig. 1e). High values on this axis correspond to high “tumoral properties”, while low values are associated with more normal features, resembling the profiles of normal brain. Accordingly, small or no changes in the metabolite slopes registered in the contralateral side of xenografts reflected absence or lower (e.g. P12, P14) invasive capacity as confirmed by histology (Fig. 1c–d). In line, no changes from baseline (“normal”) were observed on either side of the brain in mice orthotopically injected with cells from P10 that did not form any histologically detectable tumors over the observation period.

Fig. 1
figure1

Longitudinal metabolite changes indicate tumor development and invasion. a The spectra of the last scans of the injected (inj) and the non-injected contralateral side (contra) are displayed, labeled for main metabolites and their changes (indicated by arrows). The corresponding MRIs b are annotated for patient origin (P) of the xenografts (X), and mouse number (m). The histology c shows the invasive growth patterns of the corresponding representative xenografts, with close-up for characteristic areas (squares). Of note, XP14 displays migration of GBM cells across the corpus callosum and scattered infiltrated cells on the contralateral side. The human GBM cells are visualized by immunostaining against human specific nucleolin (hNCL), P53, MIB-1, or EGFR. d The longitudinal measurements of the metabolite profiles of all mice were analyzed together using STATIS. The metabolite profiles are displayed for all mice by patient against time from injection (days). Each point is a measurement of an individual mouse, identified by a specific symbol as indicated, and corresponds to the respective coordinate of the metabolite profile on the first axis of the STATIS compromise. The 21 metabolites of the profiles are represented on the first axis of STATIS as indicated in panel (E). Metabolite profiles are indicated in blue for the injected, and in pink for the contralateral side for each PDOX-series. The temporal trends are visualized by loess regression. The colors of the metabolites e correspond to their major function as indicated (energy metabolism, red; myelination, yellow; macromolecules, black; neurotransmission, green; anti-oxidation & osmoregulation, blue). Abbreviations: N-acetyl aspartate (NAA), N-acetylaspartylglutamate (NAAG), glutamate (Glu), glutamine (Gln), glycerophosphocholine (GPC), phosphocholine (PCho), glucose (Glc), glycine (Gly), myo-inositol (Ins), creatine (Cr), phosphocreatine (PCr), lactate (Lac), glutathione (GSH), taurine (Tau), alanine (Ala), aspartate (Asp), ascorbate (Asc), phosphorylethanolamine (PE), acetate (Ace), gamma aminobutyric acid (GABA), macromolecules (Mac)

The measurements preformed at ultra-short echo-time combined with the increased spectral resolution at 14.1 T allowed reporting additional changes of metabolites with low concentrations and overlapping signals that include metabolites like GABA, glutamine (Gln), GPC (dominating the increase seen in tCho), glycine (Gly) in a single measurement without the need of specific editing acquisition sequences (Fig. 1a, e, Additional file 1: Fig. S2). Of note, these metabolites have been rarely reported in previous studies mainly due to technical limitations.

The profiles of the individual metabolite concentrations, averaged over all PDOX by patient, are displayed in Additional file 1: Fig. S2. This illustrates the differences in individual metabolite changes, induced by tumor growth in the injected side of the brain and invasion of the contralateral side, between patients, and the similarity among mice injected with tumor cells of the same patient. Absence or minimal increase of lactate (Lac) was observed for xenografts derived of most patients’ tumors, in line with the invasive nature and absence of a necrotic tumor core in the mouse xenografts [29]. The growth pattern of the tumors was patient-dependent, and ranged from diffuse growth with migration across the corpus callosum and infiltration of the non-injected side to well delineated tumors with only local invasion (Fig. 1c). The PDOX from P12 yielded the most compact tumors, associated with the highest increase of Lac, tCho, Gly and decrease in Gln, Cr, PCr, Ins among all PDOX (Additional file 1: Fig. S2), and displayed some Gadolinium uptake in T1w-imaging (not shown).

Spatial metabolite patterns of human GBM resemble temporal evolution of metabolite profiles of PDOX development

To evaluate the metabolite patterns across different areas of the patients’ GBM, MVS was performed at high magnetic fields (7 T). For three patients, MVS was acquired sampling an array of 16 to 30 voxels covering distinct regions of the GBM, encompassing parts of the tumor core, the presumed infiltration zone, and adjacent brain (Fig. 2a). This allowed the measurement of 16 metabolites with their corresponding spatial information (Additional File 2), as visualized for nine individual metabolites in distinct heatmaps projected onto the respective MRI of P1 (MRSI, Additional file 1: Fig. S3). Benefiting from the increases in sensitivity and spectral resolution at 7 T together with the short-TE used, we could characterize metabolites beyond those commonly reported at 3 T (i.e. NAA, tCho, tCr, Glu + Gln and Ins), including a separate measure of Glu and Gln, and low concentration metabolites such as glycine (an important glioma marker), NAAG, PE and GSH. To determine the spatial pattern of the metabolite profiles, the MVS data of the three patients’ GBM was analyzed simultaneously by STATIS, for which the corresponding representation of the metabolites on the first axis is depicted in Fig. 2b. The metabolite profile attributed to each individual voxel, defined by their coordinates on the first axis of STATIS, was then projected as a heatmap onto the MRIs of the three patients to visualize the spatial organization of the metabolite profiles in the tumor (Fig. 2a). Globally, the spatial pattern of the metabolite profiles was characterized by low levels of NAA, and neurotransmitters, and high concentrations of GPC.PCho, Gln and lipids for voxels located in the tumor core, corresponding to high values on this first axis of STATIS (Fig. 2b) (red in the heatmap, Fig. 2a). In contrast, the probable infiltration zones and “normal” brain were reflected in low values on this first axis (blue in in the heatmap, Fig. 2a) with NAA as a prominent marker (Fig. 2b).

Fig. 2
figure2

Spatial metabolite profiles determined by multi-voxel analyses. a Multi-voxel spectroscopy (MVS) data was acquired for patients P1, P12 and P14. The metabolite spectra acquired were simultaneously analyzed by STATIS, and the respective coordinates from the first axis were then projected as a heatmap onto the corresponding voxel on the respective MRI. The color gradient corresponds to the coordinates projected on the first axis of STATIS, as indicated. The most malignant parts of the tumors are located in the “orange-red” areas. b The first axis of STATIS shows the organization of the compromise of the 16 metabolites from MVS analyses of the three patients. The color code of the metabolites corresponds to their major function as indicated. c Comparison of the first axis of the 13 common metabolites between human MVS (spatial organization) and mouse metabolite data (temporal organization) is shown in a scatter plot, and displays a remarkable similarity (Spearman correlation = -0.68, p < 0.013). Abbreviations as in Fig. 1

To elucidate the resemblance of the spectral patterns of the patients’ GBM, composed of necrotic and infiltrative tumor regions, with the spectra from the PDOXs we analyzed the 13 common metabolites (less than 50% missing values). To this end, the correlation of the coordinates of the common metabolites of the respective STATIS analyses was determined as depicted in a scatter plot (Fig. 2c). This revealed a remarkable similarity (Spearman correlation = 0.68, p < 0.013) between the spatial metabolite profiles derived from MVS across the different areas of the human glioblastoma and the temporal metabolite profiles of the PDOX, reflecting early and late stages of tumor development and invasive growth. Hence, the observed alterations of the metabolite profiles across the sampled (MVS) areas of the patients that encompassed ”normal” brain, invasion zone, and the tumor core, resembled the temporal changes of the metabolite profiles in the mouse brains during tumor development, from normal brain, invasive growth, to the full blown orthotopic xenografts at end-stage. This reflected the progression of the metabolite profiles from high NAA/low GPC.PCho (more “normal”, low values on the first axes of human MV and PDOX), to low NAA, GABA, Glu and high GPC.PCho, and Gln (more “tumoral”, high values on the first axes of human MV and PDOX), and remarkably respecting the gradient of the metabolites.

Gene expression profiles integrating human and mouse reads

To investigate the molecular underpinnings of tumor invasion the human glioblastoma (n = 8, excluding P10) and the macro-dissected PDOX of the injected and corresponding region of the contralateral side, were subjected to RNA-sequencing (Table 1, Additional File 3). The samples selected for this analysis included those from P1, 12, 14, and 7, for which we were able to obtain MVS of the patients, or for which a large number of samples were available (Table 1), The human and mouse reads, classified using the Disambiguate algorithm [1] were proportional to the presence of human cells in the mouse brain, as estimated on the DNA level with species-specific PCR and the ratio of human tumor cells/mouse cells (Spearman correlation = 0.867, p < 0.001) as determined by immunohistochemistry with a human-specific nuclear marker (NCL). Analyzing the human reads only, the genomic variant analysis (SNVs) established the filiation of the original patient tumors and the corresponding PDOX [30]. Characteristic molecular features of the parental tumors, such as mutations or previously described expression signatures [38] (e.g. associated with EGFR overexpression), were mostly retained in the PDOX (Additional file 1: Fig. S4).

Integration of mouse and human reads, expectedly revealed that the first axis of the PCA of all samples, based on a gene set selected by sparse PCA (SPCA, 247 genes, comprising 134 human and 113 mouse genes, selected with lasso regression, consolidated bootstrap; see flowchart of analysis in Additional file 1: Fig. 1a) was dominated by the species origin, which at the same time reflected the tissue type (tumor vs mouse brain) (Additional file 1: Fig. 5, Additional File 4). In line, pathway analysis, using gene set enrichment analysis (GSEA) and the molecular signature database (MSigDB), for which the mouse genes were converted into the human homologs, revealed that the top 25 pathways were dominated by cell cycle, proliferation, EGFR signaling, and tumor progression associated gene signatures, and were expectedly mainly described by the human/tumor derived genes (Additional File 5). Interestingly, a tumor invasion related signature was among these top 25 pathways, with similar contributions from both mouse/host-derived and human/tumor-derived genes.

Molecular pathways associated with changes in metabolite profiles (1H MRS) of invasive tumors.

The main interest in this study was to investigate molecular changes underlying tumor development and invasion by integrating gene expression and metabolite information. The latter is amenable to non-invasive analyses allowing longitudinal monitoring of tumor progression and invasion with potential for clinical use.

We determined gene expression profiles associated with the metabolite profiles in the PDOXs acquired from the injected and contralateral side at the last MRI/S scan (end-stage) (see flowchart of analysis in Additional file 1: Fig. S1C). The patient specific effects were removed in the expression data using within-group PCA (WCA) [43] in order to focus on features of tumor host interaction. A set of 185 genes associated with the metabolite profiles was selected using sparse Partial Least Squares analysis (SPLS) [3, 32], consolidated by bootstrap) that allows combination of different data types linked to the same samples (common column of both tables). The gene set comprised 60 mouse and 125 human genes (Additional File 6). A heatmap illustrates the strong association between the 13 metabolites obtained by MRS and the 185 selected genes in the cross platform comparison (cross table coinertia analysis; coefficient of vectorial correlation, RV 0.73, p = 0.01, Fig. 3a). Representing the PDOX samples in function of their metabolite and gene expression profiles (defined by their coordinates on the first axes of the coinertia analysis) revealed a clear gradient following tumor progression, from normal brain (high in NAA) to tumor (low in NAA, and high in choline compounds and enhanced Lac), and change of gene expression (Fig. 3b, c). The extreme features were most prominent for XP12 that yielded the most compact xenografts with some contrast enhancement, while in the contralateral side tumor spread was neither detectable by MRS nor subsequent histology. Of note, the detection of metabolites by MRS is agnostic to the species origin. The RNA coverage revealed a drift from mouse to human reads during tumor development, as expected.

Fig. 3
figure3

Correlation structure between metabolite profiles (1H-MRS) and associated transcriptome (human and mouse). The heatmap a illustrates the correlations between 13 metabolites from the last scans of the xenograft bearing mice (injected and contralateral side) and the 185 metabolite associated genes selected by SPLS and retained by the bootstrap procedure (≥ 0.1). The correlation matrix (coinertia) between genes (rows) and metabolites (columns) is ordered by the 1st axes obtained from coinertia. For each gene the species origin (human or mouse) and the frequency of selection by bootstrap are annotated on the right. b The relation between metabolite profiles and gene expression is visualized for all PDOX samples per patient on the vectorial plane defined by the coordinates on the first axes of gene expression and the metabolite profiles from the coinertia analysis, respectively. The samples from the injected side are represented by squares, and by circles for the contralateral side, the color gradient of the symbols indicates the percentage of human reads. c The metabolites projected on the first axis of the coinertia analysis. d Similarly, the correlation structure between metabolite profiles (1H-MRS) and averaged expression of the significant pathways emerging from GSEA (p ≤ 0.1) using the MSigDB is illustrated in a heatmap (d). The pathways are annotated with the adjusted p-value (p ≤ 0.1), and the proportion of human genes contributing to the pathway

Pathway analysis of the 185 selected genes using the MSigDB database, identified a set of 24 significant pathways (p ≤ 0.1, Bonferroni adjusted). The majority of the pathways were associated with extracellular matrix, extracellular structure, tissue remodeling, adhesion, morphogenesis and remodeling of blood vessels, multi-cancer signature(s) of invasion, metastasis, epithelial mesenchymal transition, and including a mesenchymal glioblastoma signature (Additional File 7). Interestingly, expression of both, mouse and human genes contributed to the selection of all these pathways, ranging from 20 to 80% human genes, depending on the pathway. A heatmap visualizes the correlation of the pathways with the metabolite profiles (Fig. 3d). Interestingly, the gene set linked to epithelial to mesenchymal transition was associated with early stages, while gene sets such as those for matrix remodeling, and mesenchymal glioblastoma were associated with metabolite profiles of more advanced stages of tumor growth.

Pathway analysis of the invaded brain

In order to determine the molecular effects of glioblastoma invasion on the mouse brain, we specifically analyzed the mouse reads that originate from mouse derived cells enclosed in the macro-dissected xenografts, and the mouse cells of the respective mirrored region from the contralateral side, invaded or not by GBM cells. Two samples with < 106 mouse reads were excluded. Both samples were from the injected side of the PDOX, derived from P12 that are highly compact and therefore comprise only few mouse cells (Fig. 1). The heatmap of genes selected by SPCA (n = 208, consolidated by bootstrap; Additional File 4; see flowchart of analysis, Additional file 1: Fig. S1B) yielded 3 major gene clusters (Fig. 4a). The organization of the samples by gene expression seemed to be driven by the extent and pattern of the tumor invasion the mouse brains were exposed to. This intriguing observation seems impacted by the patient origin of the injected tumor cells, despite the fact that GBM-derived human reads were excluded from the analysis. A PCA illustrates the samples projected onto the first 2 axes of gene expression (Fig. 4b–d). The first axis explains the difference between the injected and the contralateral side dominated by cluster 3 genes. The gradient of expression of these genes are explained by the extent of invasion, as measured by combination of the percentage of human reads in the samples, and the PDOX type (injected or contralateral side) and explains more than 50% of the variance (combined variation fraction > 0.5; Fig. 4e). The pathways associated with gene expression of cluster 3 are related to immune response, regulation of cytokine production, including interferon-α interleukin 6, and cytokine signaling (e.g. via C-X-C Motif Chemokine Receptor 1, CXCR1), and inflammatory response such as response to interferon-γ (Fig. 4f, Additional File 5). Interestingly, it includes the mesenchymal GBM signature as significant pathway that we also identified as significant, when selecting metabolism associated gene expression signatures, as described above (Fig. 3). Of note, this signature, included in the MSiGDB, corresponds to the signature associated with the expression-based classification of the mesenchymal subtype of GBM [60]. The second axis of the PCA seems to be driven by the impact of PDOX of patient P7 (plus 1 PDOX of P14), with a particular expression pattern captured in cluster 2 that is independent of the extent of tumor invasion. A characteristic feature of these PDOX was rapid growth and/or massive invasion of both hemispheres (Fig. 1). The pathways selected with cluster 2 genes were dominated by ribosomal genes that were highly redundant among the numerous selected pathways. They represented generic protein and RNA related processes, such as regulation of translation, cellular localization of proteins, catabolic processes, and pathways linked to infectious disease.

Fig. 4
figure4

Effect of tumor invasion on the mouse brain transcriptome. Mouse gene expression profiles included PODX samples and one mock injected mouse brain (mB). The heatmap (a) illustrates the normalized expression of the 208 mouse genes selected by SPCA and consolidated by bootstrap (≥ 0.1). The genes were classified in 3 clusters (consensus k-means clustering for 100 repetitions). (B-D) The samples are projected onto the first two axes of the PCA for the selected genes, where the human read proportion b, the sample type c and tumor origin d were added as supplementary variables, annotated with the color code. The contribution of each gene cluster to explain the variance is evaluated by variation partitioning e represented by a Venn diagram of variation fractions (percentage) for the three supplementary variables (frequency of human reads, FreqHum; sample origin, Origin; type if tissue, Type). (F) The averaged gene expression of the significant pathways emerging from GSEA (p ≤ 0.1) of cluster 3 is illustrated in a heatmap. The pathways are annotated with the adjusted p-value and the gene ratio (number of selected genes in pathway/number of selected genes in the analysis)

Discussion

The invasive capacity of GBM plays a key role in the aggressiveness of the tumor, its resistance to treatment, recurrence and poor prognosis. The invasive component is typically shielded by the BBB, and it remains a challenge to monitor tumor infiltration of the brain parenchyma using standard MRI techniques [16, 62]. It is therefore of crucial importance to develop adequate tumor models featuring this infiltrative part that is considered highly relevant for tumor recurrence.

We evaluated molecular features of infiltrative PDOX that may serve as proxy for studying the infiltrative front of GBM. To this aim, we present the first in vivo comparison of human tumors and respective PDOX in the mouse brain on the levels of radiological behavior and metabolism, as determined by ultra-high field 1H-MRS of the patients (7 T) and the respective PDOX (14.1 T). We found a good concordance between the temporal changes of the metabolite profiles observed during invasive growth of the PDOX, and the metabolite profiles reflecting the spatial organization of MVS originating from different parts of the human tumors, comprising normal brain, infiltration zone, and the tumor core (Fig. 2c). Hence, metabolic changes may inform on alterations in the infiltration zone that is not visible on routine MRI evaluation. Importantly, highly resolved 1H-MRS was more sensitive to detect early development of the highly infiltrative PDOX than conventional structural MRI acquisitions. This may open the possibility to test and non-invasively monitor early treatment effects in the infiltrative zone shielded by an intact BBB.

Interrogating the molecular mechanisms related to tumor invasion and infiltrative growth by means of multi-dimensional analysis allowed simultaneous exploration of metabolic and transcriptomic changes of the samples and respective associated pathways. This combined analysis provided insights into the biological mechanisms that are underlying the metabolic changes that can be followed non-invasively. Most importantly, both tumor and host contributed to gene expression profiles associated with the biological pathways uncovered, supportive of their biological relevance. These signatures indicated active processes associated with changes of the extracellular matrix, tissue and blood vessel remodeling, along with signatures attributed to tumor invasion, metastasis, and mesenchymal GBM. Lately, the tumor matrix, and the associated cell–matrix interaction have received more attention in solid extra cerebral tumors [11], while little remains known in brain tumors [41]. Interestingly, gene expression annotated for hallmarks of epithelial/mesenchymal transition was associated with more normal brain-like metabolic features, while the gene expression related to matrisome, metastasis, and mesenchymal GBM were more strongly correlated with metabolite profiles of more advanced tumors (Fig. 3d).

Taking a different view, and evaluating the transcriptome originating from the invaded brain only (mouse reads, Fig. 4), revealed inflammatory signatures and profiles of cytokine mediated regulation of immune response. The mesenchymal GBM signature that was again associated with a more aggressive extent of tumor invasion/tumor burden as estimated by the proportion of human reads, despite the fact that they were excluded from the analysis. The mesenchymal GBM subtype has been associated with more prominent recruitment of macrophages [61] for which the functional interaction with the GBM cells has been described recently, evoking that this interaction drives the mesenchymal-like state of GBM [20]. Along the same lines, testing our previously reported GBM gene expression signatures [38] with GSEA, we found two matching clusters that were significantly associated with the expression profiles of the invaded brain: an interferon-induced gene signature (G12; adjusted p-value, 0.0002), and an innate immune response signature (G24; adjusted p-value, 0.0001). This suggests some concordance of the brain reacting to tumor invasion in the present PDOX model and human GBM.

These molecular insights suggest that the tumor/host interaction of tumor invasion can be modeled by monitoring metabolic changes in PDOX and that these changes present a good match with their human counterparts as described in this study. The insights derived from temporal metabolic changes associated with diffuse tumor invasion and progression may be applicable to help monitor the invasion zone in patients to identify early responses to treatment or on the contrary, early tumor progression. The obvious limitations are the lack of an intact immune system in this mouse model that is obviously expected to play an important role in tumor progression and treatment [28].

The PDOX may represent an attractive perspective to develop and evaluate drugs aimed at treating tumor cells in the invasion zone. Other invasive orthotopic brain tumor models may also be suitable for metabolic monitoring. For instance, we have reported similar longitudinal metabolic changes using established GBM derived sphere (GS) lines yielding highly infiltrative orthotopic xenografts [29]. The changes in the metabolic profiles were sensitive enough to follow the distinct evolution when transducing the GS cells with a tumor suppressor gene (Wnt inhibitory factor 1) [29]. This is important, as the median latency of freshly resected patient derived tumor cells to fully develop into PDOX is generally too long (5–35 weeks in this study) to serve as faithful avatars for patient specific testing of targeted drugs to guide treatment choices.

Taken together, invasive PDOX models exert spectroscopic and transcriptomic features of brain infiltration that show similarities with the presumed infiltration zone of GBM. This supports their suitability as relevant models for studying the non-enhancing part of GBM. The possibility of non-invasive in vivo monitoring of invasive growth in this difficult to evaluate compartment, may allow early detection of relapse and monitoring of treatment effects of novel drugs that eventually may be translated into the human setting.

Availability of data and materials

The datasets produced, used, and analyzed during the current study available from the corresponding author on reasonable request.

Abbreviations

MRI:

Magnetic resonance imaging

PDOX:

Patient derived orthotopic xenografts

1H MRS:

1H magnetic resonance spectroscopy

GBM:

Glioblastoma

BBB:

Blood brain barrier

VOI:

Volume of interest

SVS:

Single voxel spectroscopy

MVS:

Multi-voxel spectroscopy

MRSI:

Magnetic resonance spectroscopic imaging

MVS:

Multi voxel speczroscopy

CRLB:

Cramer-Rao lower bound

NAA:

N-acetyl aspartate

NAAG:

N-acetylaspartylglutamate

Glu:

Glutamate

Gln:

Glutamine

GPC:

Glycerophosphocholine

PCho:

Phosphocholine

Glc:

Glucose

Gly:

Glycine

Ins:

Myo-inositol

Cr:

Creatine

PCr:

Phosphocreatine

Lac:

Lactate

GSH:

Glutathione

Tau:

Taurine

Ala:

Alanine

Asp:

Aspartate

Asc:

Ascorbate

PE:

Phosphorylethanolamine

Ace:

Acetate

GABA:

Gamma aminobutyric acid

Mac:

Macromolecules

H&E:

Hematoxylin and eosin

NCL:

Nucleolin

EGFR:

Epidermal growth factor receptor

RNAseq:

RNA sequencing

SNVs:

Single nucleotide variants

COSMIC:

Catalogue of Somatic Mutations in Cancer

TCGA:

The Cancer Genome Atlas

PCA:

Principal component analysis

SPCA:

Sparse PCA

SPLS:

Sparse partial least squares

WCA:

Within-group PCA

GSEA:

Gene set enrichment analysis

MSigDB:

Molecular signature database

IDHmt:

Isocitrate dehydrogenase mutant

1p/19q non-codeleted:

Chromosome arms 1p and 19q non-codeleted

References

  1. 1.

    Ahdesmaki MJ, Gray SR, Johnson JH, Lai Z (2016) Disambiguate: An open-source application for disambiguating two species in next generation sequencing data from grafted samples. F1000Res 5:2741. https://doi.org/10.12688/f1000research.10082.2

    CAS  Article  PubMed  Google Scholar 

  2. 2.

    Alcoser SY, Kimmel DJ, Borgel SD, Carter JP, Dougherty KM, Hollingshead MG (2011) Real-time PCR-based assay to quantify the relative amount of human and mouse tissue present in tumor xenografts. BMC Biotechnol 11:124–124. https://doi.org/10.1186/1472-6750-11-124

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  3. 3.

    Allen GI, Peterson C, Vannucci M, Maletic-Savatic M (2013) Regularized partial least squares with an application to NMR spectroscopy. Stat Anal Data Min 6:302–314. https://doi.org/10.1002/sam.11169

    Article  PubMed  Google Scholar 

  4. 4.

    Bao S, Wu Q, McLendon RE, Hao Y, Shi Q, Hjelmeland AB, Dewhirst MW, Bigner DD, Rich JN (2006) Glioma stem cells promote radioresistance by preferential activation of the DNA damage response. Nature 444:756–760. https://doi.org/10.1038/nature05236

    CAS  Article  PubMed  Google Scholar 

  5. 5.

    Bougnaud S, Golebiewska A, Oudin A, Keunen O, Harter PN, Mader L, Azuaje F, Fritah S, Stieber D, Kaoma T et al (2016) Molecular crosstalk between tumour and brain parenchyma instructs histopathological features in glioblastoma. Oncotarget 7:31955–31971. https://doi.org/10.18632/oncotarget.7454

    Article  PubMed  PubMed Central  Google Scholar 

  6. 6.

    Ceccarelli M, Barthel FP, Malta TM, Sabedot TS, Salama SR, Murray BA, Morozova O, Newton Y, Radenbaugh A, Pagnotta SM et al (2016) Molecular profiling reveals biologically discrete subsets and pathways of progression in diffuse glioma. Cell 164:550–563. https://doi.org/10.1016/j.cell.2015.12.028

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  7. 7.

    Chessel D, Dufour AB, Thioulouse J (2004) The ade4 package-I: One-table methods. R News 4:5–10

    Google Scholar 

  8. 8.

    Choi C, Ganji SK, DeBerardinis RJ, Hatanpaa KJ, Rakheja D, Kovacs Z, Yang XL, Mashimo T, Raisanen JM, Marin-Valencia I et al (2012) 2-hydroxyglutarate detection by magnetic resonance spectroscopy in IDH-mutated patients with gliomas. Nat Med 18:624–629. https://doi.org/10.1038/nm.2682

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  9. 9.

    Cingolani P, Platts A, le Wang L, Coon M, Nguyen T, Wang L, Land SJ, Lu X, Ruden DM (2012) A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin) 6:80–92. https://doi.org/10.4161/fly.19695

    CAS  Article  Google Scholar 

  10. 10.

    Coulibaly B, Nanni I, Quilichini B, Gaudart J, Metellus P, Fina F, Boucard C, Chinot O, Ouafik L, Figarella-Branger D (2010) Epidermal growth factor receptor in glioblastomas: correlation between gene copy number and protein expression. Hum Pathol 41:815–823. https://doi.org/10.1016/j.humpath.2009.09.020

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Cox TR (2021) The matrix in cancer. Nat Rev Cancer 21:217–238. https://doi.org/10.1038/s41568-020-00329-7

    CAS  Article  PubMed  Google Scholar 

  12. 12.

    Cudalbu C, Behar KL, Bhattacharyya PK, Bogner W, Borbath T, de Graaf RA, Gruetter R, Henning A, Juchem C, Kreis R et al (2021) Contribution of macromolecules to brain (1) H MR spectra: Experts’ consensus recommendations. NMR Biomed 34:e4393. https://doi.org/10.1002/nbm.4393

    Article  PubMed  Google Scholar 

  13. 13.

    Dirkse A, Golebiewska A, Buder T, Nazarov PV, Muller A, Poovathingal S, Brons NHC, Leite S, Sauvageot N, Sarkisjan D et al (2019) Stem cell-associated heterogeneity in glioblastoma results from intrinsic tumor plasticity shaped by the microenvironment. Nat Commun 10:1787. https://doi.org/10.1038/s41467-019-09853-z

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Doledec S, Chessel D (1994) Co-Inertia analysis - an alternative method for studying species environment relationships. Freshw Biol 31:277–294. https://doi.org/10.1111/j.1365-2427.1994.tb01741.x

    Article  Google Scholar 

  15. 15.

    Durinck S, Moreau Y, Kasprzyk A, Davis S, De Moor B, Brazma A, Huber W (2005) BioMart and Bioconductor: a powerful link between biological databases and microarray data analysis. Bioinformatics 21:3439–3440. https://doi.org/10.1093/bioinformatics/bti525

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Eidel O, Burth S, Neumann JO, Kieslich PJ, Sahm F, Jungk C, Kickingereder P, Bickelhaupt S, Mundiyanapurath S, Baumer P et al (2017) Tumor infiltration in enhancing and non-enhancing parts of glioblastoma: A correlation with histopathology. PLoS ONE 12:e0169292. https://doi.org/10.1371/journal.pone.0169292

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  17. 17.

    Golebiewska A, Hau AC, Oudin A, Stieber D, Yabo YA, Baus V, Barthelemy V, Klein E, Bougnaud S, Keunen O et al (2020) Patient-derived organoids and orthotopic xenografts of primary and recurrent gliomas represent relevant patient avatars for precision oncology. Acta Neuropathol 140:919–949. https://doi.org/10.1007/s00401-020-02226-7

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  18. 18.

    Gruetter R (1993) Automatic, localized in vivo adjustment of all first- and second-order shim coils. Magn Reson Med 29:804–811. https://doi.org/10.1002/mrm.1910290613

    CAS  Article  PubMed  Google Scholar 

  19. 19.

    Gusyatiner O, Bady P, Pham MDT, Lei Y, Park J, Daniel RT, Delorenzi M, Hegi ME (2021) BET inhibitors repress expression of Interferon-stimulated genes and synergize with HDAC inhibitors in glioblastoma. Neuro Oncol. https://doi.org/10.1093/neuonc/noab115

    Article  PubMed  Google Scholar 

  20. 20.

    Hara T, Chanoch-Myers R, Mathewson ND, Myskiw C, Atta L, Bussema L, Eichhorn SW, Greenwald AC, Kinker GS, Rodman C et al (2021) Interactions between cancer cells and immune cells drive transitions to mesenchymal-like states in glioblastoma. Cancer Cell 39(779–792):e711. https://doi.org/10.1016/j.ccell.2021.05.002

    CAS  Article  Google Scholar 

  21. 21.

    Heo M, Gabriel KR (2010) A permutation test of association between configurations by means of the rv coefficient. Commun Statist Simul Comput 27:843–856. https://doi.org/10.1080/03610919808813512

    Article  Google Scholar 

  22. 22.

    Horska A, Barker PB (2010) Imaging of brain tumors: MR spectroscopy and metabolic imaging. Neuroimaging Clin N Am 20:293–310. https://doi.org/10.1016/j.nic.2010.04.003

    Article  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Hulsey KM, Mashimo T, Banerjee A, Soesbe TC, Spence JS, Vemireddy V, Maher EA, Bachoo RM, Choi C (2015) (1)H MRS characterization of neurochemical profiles in orthotopic mouse models of human brain tumors. NMR Biomed 28:108–115. https://doi.org/10.1002/nbm.3231

    CAS  Article  PubMed  Google Scholar 

  24. 24.

    Jacob F, Salinas RD, Zhang DY, Nguyen PTT, Schnoll JG, Wong SZH, Thokala R, Sheikh S, Saxena D, Prokop S et al (2020) A patient-derived glioblastoma organoid model and biobank recapitulates inter- and intra-tumoral heterogeneity. Cell 180(188–204):e122. https://doi.org/10.1016/j.cell.2019.11.036

    CAS  Article  Google Scholar 

  25. 25.

    Josse J, Husson F (2012) Handling missing values in exploratory multivariate data analysis methods. Journal de la Société Francaise de Statistique 153:1–21

    Google Scholar 

  26. 26.

    Kathagen-Buhmann A, Schulte A, Weller J, Holz M, Herold-Mende C, Glass R, Lamszus K (2016) Glycolysis and the pentose phosphate pathway are differentially associated with the dichotomous regulation of glioblastoma cell migration versus proliferation. Neuro Oncol 18:1219–1229. https://doi.org/10.1093/neuonc/now024

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Keunen O, Taxt T, Gruner R, Lund-Johansen M, Tonn JC, Pavlin T, Bjerkvig R, Niclou SP, Thorsen F (2014) Multimodal imaging of gliomas in the context of evolving cellular and molecular therapies. Adv Drug Deliv Rev 76:98–115. https://doi.org/10.1016/j.addr.2014.07.010

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Klemm F, Maas RR, Bowman RL, Kornete M, Soukup K, Nassiri S, Brouland JP, Iacobuzio-Donahue CA, Brennan C, Tabar V et al (2020) Interrogation of the microenvironmental landscape in brain tumors reveals disease-specific alterations of immune cells. Cell 181(1643–1660):e1617. https://doi.org/10.1016/j.cell.2020.05.007

    CAS  Article  Google Scholar 

  29. 29.

    Lai M, Vassallo I, Lanz B, Poitry-Yamate C, Hamou MF, Cudalbu C, Gruetter R, Hegi ME (2018) In vivo characterization of brain metabolism by (1) H MRS, (13) C MRS and (18) FDG PET reveals significant glucose oxidation of invasively growing glioma cells. Int J Cancer 143:127–138. https://doi.org/10.1002/ijc.31299

    CAS  Article  PubMed  Google Scholar 

  30. 30.

    Lai Z, Markovets A, Ahdesmaki M, Chapman B, Hofmann O, McEwen R, Johnson J, Dougherty B, Barrett JC, Dry JR (2016) VarDict: a novel and versatile variant caller for next-generation sequencing in cancer research. Nucleic Acids Res 44:e108. https://doi.org/10.1093/nar/gkw227

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Lavit C, Escoufier Y, Sabatier R, Traissac P (1994) The Act (Statis Method). Comput Stat Data Anal 18:97–119. https://doi.org/10.1016/0167-9473(94)90134-1

    Article  Google Scholar 

  32. 32.

    Le Cao KA, Rossouw D, Robert-Granie C, Besse P (2008) A sparse PLS for variable selection when integrating omics data. Stat Appl Genet Mol Biol 7, Article 35. https://doi.org/10.2202/1544-6115.1390

  33. 33.

    Liau BB, Sievers C, Donohue LK, Gillespie SM, Flavahan WA, Miller TE, Venteicher AS, Hebert CH, Carey CD, Rodig SJ et al (2017) Adaptive chromatin remodeling drives glioblastoma stem cell plasticity and drug tolerance. Cell Stem Cell 20(233–246):e237. https://doi.org/10.1016/j.stem.2016.11.003

    CAS  Article  Google Scholar 

  34. 34.

    Marques JP, Kober T, Krueger G, van der Zwaag W, Van de Moortele PF, Gruetter R (2010) MP2RAGE, a self bias-field corrected sequence for improved segmentation and T1-mapping at high field. Neuroimage 49:1271–1281. https://doi.org/10.1016/j.neuroimage.2009.10.002

    Article  PubMed  Google Scholar 

  35. 35.

    Mashimo T, Pichumani K, Vemireddy V, Hatanpaa KJ, Singh DK, Sirasanagandla S, Nannepaga S, Piccirillo SG, Kovacs Z, Foong C et al (2014) Acetate is a bioenergetic substrate for human glioblastoma and brain metastases. Cell 159:1603–1614. https://doi.org/10.1016/j.cell.2014.11.025

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Mishkovsky M, Gusyatiner O, Lanz B, Cudalbu C, Vassallo I, Hamou MF, Bloch J, Comment A, Gruetter R, Hegi ME (2021) Hyperpolarized (13)C-glucose magnetic resonance highlights reduced aerobic glycolysis in vivo in infiltrative glioblastoma. Sci Rep 11:5771. https://doi.org/10.1038/s41598-021-85339-7

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  37. 37.

    Mlynarik V, Cudalbu C, Xin L, Gruetter R (2008) 1H NMR spectroscopy of rat brain in vivo at 14.1Tesla: improvements in quantification of the neurochemical profile. J Magn Reson 194:163–168. https://doi.org/10.1016/j.jmr.2008.06.019

    CAS  Article  PubMed  Google Scholar 

  38. 38.

    Murat A, Migliavacca E, Gorlia T, Lambiv WL, Shay T, Hamou MF, de Tribolet N, Regli L, Wick W, Kouwenhoven MC et al (2008) Stem cell-related “self-renewal” signature and high epidermal growth factor receptor expression associated with resistance to concomitant chemoradiotherapy in glioblastoma. J Clin Oncol 26:3015–3024. https://doi.org/10.1200/JCO.2007.15.7164

    CAS  Article  PubMed  Google Scholar 

  39. 39.

    Neftel C, Laffy J, Filbin MG, Hara T, Shore ME, Rahme GJ, Richman AR, Silverbush D, Shaw ML, Hebert CM et al (2019) An integrative model of cellular states, plasticity, and genetics for glioblastoma. Cell 178(835–849):e821. https://doi.org/10.1016/j.cell.2019.06.024

    CAS  Article  Google Scholar 

  40. 40.

    Provencher SW (2001) Automatic quantitation of localized in vivo 1H spectra with LCModel. NMR Biomed 14:260–264. https://doi.org/10.1002/nbm.698

    CAS  Article  PubMed  Google Scholar 

  41. 41.

    Quail DF, Joyce JA (2017) The microenvironmental landscape of brain tumors. Cancer Cell 31:326–341. https://doi.org/10.1016/j.ccell.2017.02.009

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  42. 42.

    R Core Team (2020) R: A language and environment for statistical computing. https://www.R-project.org2020

  43. 43.

    Rao CR (1964) The use and interpretation of principal component analysis in applied research. Sankhya Ser A 26:329–358

    Google Scholar 

  44. 44.

    Robinson MD, McCarthy DJ, Smyth GK (2010) edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26:139–140. https://doi.org/10.1093/bioinformatics/btp616

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  45. 45.

    Robinson MD, Oshlack A (2010) A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol 11:R25. https://doi.org/10.1186/gb-2010-11-3-r25

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  46. 46.

    Romesburg HC (1985) Exploring, confirming, and randomization tests. Comput Geosci 11:19–37. https://doi.org/10.1016/0098-3004(85)90036-6

    Article  Google Scholar 

  47. 47.

    Rosol M, Harutyunyan I, Xu J, Melendez E, Smbatyan G, Finlay JL, Krieger MD, Gonzalez-Gomez I, Reynolds CP, Nelson MD et al (2009) Metabolism of orthotopic mouse brain tumor models. Mol Imaging 8:199–208. https://doi.org/10.2310/7290.2009.00019

    Article  PubMed  Google Scholar 

  48. 48.

    Schaller B, Xin L, Gruetter R (2014) Is the macromolecule signal tissue-specific in healthy human brain? A (1)H MRS study at 7 Tesla in the occipital lobe. Magn Reson Med 72:934–940. https://doi.org/10.1002/mrm.24995

    CAS  Article  PubMed  Google Scholar 

  49. 49.

    Sciuscio D, Diserens AC, van Dommelen K, Martinet D, Jones G, Janzer RC, Pollo C, Hamou MF, Kaina B, Stupp R et al (2011) Extent and patterns of MGMT promoter methylation in glioblastoma- and respective glioblastoma-derived spheres. Clin Cancer Res 17:255–266. https://doi.org/10.1158/1078-0432.CCR-10-1931

    CAS  Article  PubMed  Google Scholar 

  50. 50.

    Shen HP, Huang JHZ (2008) Sparse principal component analysis via regularized low rank matrix approximation. J Multivar Anal 99:1015–1034. https://doi.org/10.1016/j.jmva.2007.06.007

    Article  Google Scholar 

  51. 51.

    Stieber D, Golebiewska A, Evers L, Lenkiewicz E, Brons NH, Nicot N, Oudin A, Bougnaud S, Hertel F, Bjerkvig R et al (2014) Glioblastomas are composed of genetically divergent clones with distinct tumourigenic potential and variable stem cell-associated phenotypes. Acta Neuropathol 127:203–219. https://doi.org/10.1007/s00401-013-1196-4

    CAS  Article  PubMed  Google Scholar 

  52. 52.

    Stupp R, Taillibert S, Kanner A, Read W, Steinberg D, Lhermitte B, Toms S, Idbaih A, Ahluwalia MS, Fink K et al (2017) Effect of tumor-treating fields plus maintenance temozolomide vs maintenance temozolomide alone on survival in patients with glioblastoma: a randomized clinical trial. JAMA 318:2306–2316. https://doi.org/10.1001/jama.2017.18718

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  53. 53.

    Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES et al (2005) Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA 102:15545–15550. https://doi.org/10.1073/pnas.0506580102

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Tate JG, Bamford S, Jubb HC, Sondka Z, Beare DM, Bindal N, Boutselakis H, Cole CG, Creatore C, Dawson E et al (2019) COSMIC: the catalogue of somatic mutations in cancer. Nucleic Acids Res 47:D941–D947. https://doi.org/10.1093/nar/gky1015

    CAS  Article  PubMed  Google Scholar 

  55. 55.

    Therneau TM, Grambsch PM (2000) Modeling survival data: extending the cox model. Springer, Berlin

    Book  Google Scholar 

  56. 56.

    Tkac I, Deelchand D, Dreher W, Hetherington H, Kreis R, Kumaragamage C, Povazan M, Spielman DM, Strasser B, de Graaf RA (2021) Water and lipid suppression techniques for advanced (1) H MRS and MRSI of the human brain: Experts’ consensus recommendations. NMR Biomed 34:e4459. https://doi.org/10.1002/nbm.4459

    CAS  Article  PubMed  Google Scholar 

  57. 57.

    Van Deun K, Smilde AK, van der Werf MJ, Kiers HA, Van Mechelen I (2009) A structured overview of simultaneous component based data integration. BMC Bioinf 10:246. https://doi.org/10.1186/1471-2105-10-246

    CAS  Article  Google Scholar 

  58. 58.

    Vassallo I, Zinn P, Lai M, Rajakannu P, Hamou MF, Hegi ME (2016) WIF1 re-expression in glioblastoma inhibits migration through attenuation of non-canonical WNT signaling by downregulating the lncRNA MALAT1. Oncogene 35:12–21. https://doi.org/10.1038/onc.2015.61

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Vaubel RA, Tian S, Remonde D, Schroeder MA, Mladek AC, Kitange GJ, Caron A, Kollmeyer TM, Grove R, Peng S et al (2020) Genomic and phenotypic characterization of a broad panel of patient-derived xenografts reflects the diversity of glioblastoma. Clin Cancer Res 26:1094–1104. https://doi.org/10.1158/1078-0432.CCR-19-0909

    CAS  Article  PubMed  Google Scholar 

  60. 60.

    Verhaak RG, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, Miller CR, Ding L, Golub T, Mesirov JP et al (2010) Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell 17:98–110. https://doi.org/10.1016/j.ccr.2009.12.020

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  61. 61.

    Wang Q, Hu B, Hu X, Kim H, Squatrito M, Scarpace L, deCarvalho AC, Lyu S, Li P, Li Y et al (2017) Tumor evolution of glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. Cancer Cell 32(42–56):e46. https://doi.org/10.1016/j.ccell.2017.06.003

    CAS  Article  Google Scholar 

  62. 62.

    Warntjes M, Blystad I, Tisell A, Larsson EM (2018) Synthesizing a contrast-enhancement map in patients with high-grade gliomas based on a postcontrast MR Imaging quantification only. AJNR Am J Neuroradiol 39:2194–2199. https://doi.org/10.3174/ajnr.A5870

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  63. 63.

    Xin L, Mlynarik V, Gruetter R (2012) A new approach to short-TE full-sensitivity MRSI of human brain at 7T. Proc Intl Soc Mag Reson Med 20th Annual Meeting & Exhibition, Melbourne, Australia

  64. 64.

    Xin L, Schaller B, Mlynarik V, Lu H, Gruetter R (2013) Proton T1 relaxation times of metabolites in human occipital white and gray matter at 7 T. Magn Reson Med 69:931–936. https://doi.org/10.1002/mrm.24352

    CAS  Article  PubMed  Google Scholar 

Download references

Acknowledgements

We thank the patients and their families for their support and participation. Information about TCGA can be found at “http://cancergenome.nih.gov”. We acknowledge access to the facilities, technical platforms and expertise of the CHUV, UNIL, SIB, and the CIBM.

Funding

The study was funded by the Swiss Bridge Award, the Swiss Cancer Ligue (KFS-3998–08-2016, KFS-4461–02-2018), the Swiss National Science Foundation (SNF-3103–163297, SNF-3103–182821), Harry’s ride and Cadot foundations.

Author information

Affiliations

Authors

Contributions

C.C., A.F.H, and M.E.H conceived and designed the study; P.B. designed, performed the biostatistical analyses, and interpretation; L.X. performed MRI/MRS (7 T) of patients; C.C. and M.La. performed MRI/MRS of mice; O.G., and M.-F.H. performed transplantations, and sample analyses; M.Le. assisted animal work; R.T.D. provided surgical specimen and J.-P.B. neuropathology expertise; A.F.H. lead the patient study; M.E.H. coordinated data integration, wrote the manuscript with P.B.. All authors contributed to data interpretation and manuscript writing. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Monika E. Hegi.

Ethics declarations

Ethics approval and consent to participate

Patients planned for surgery of a suspected GBM at the Lausanne University Hospital (CHUV) were enrolled (clinicaltrials.gov, NCT02904525) with written informed consent. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the local ethics committee CER-VD (F-25/99, 268/14). Animal protocols were designed under 3R considerations and all procedures were performed under anesthesia/ analgesia, and protocols were approved by the concerned Swiss authorities (VD-1181_6; VD-2777).

Consent for publication

Not applicable.

Competing interests

The authors declare to have no competing interests.

Additional information

Publisher's Note

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

Andreas F. Hottinger and Monika E. Hegi share senior authorship

Supplementary Information

Additional file 1.

 Fig. S1 Flow charts for gene selection procedure. Fig. S2. Longitudinal variations of metabolite concentrations of PDOX. Fig. S3. Multivoxel 1H MRS at 7 Tesla of patient 1. Fig. S4. Retention of characteristic molecular features in PDOX. Fig. S5. Expression analyses including all PDOX and original human GBMs.

Additional file 2.

Metabolite datasets.

Additional file 3.

Description samples (RNAseq, NB reads, categories, etc.).

Additional file 4. 

List of candidate genes (n=247) for aggregated human and mouse RNAseq data; Candidate genes (n=208) mouse only.

Additional file 5. 

Pathways associated with gene expression using all samples. Pathways associated with gene expression restricted to mouse sequences.

Additional file 6. 

Candidate genes (n=185) associated with metabolites.

Additional file 7. 

Pathways linked with metabolite associated gene expression.

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

Verify currency and authenticity via CrossMark

Cite this article

Cudalbu, C., Bady, P., Lai, M. et al. Metabolic and transcriptomic profiles of glioblastoma invasion revealed by comparisons between patients and corresponding orthotopic xenografts in mice. acta neuropathol commun 9, 133 (2021). https://doi.org/10.1186/s40478-021-01232-4

Download citation

Keywords

  • Glioblastoma
  • Invasion
  • 1H MRS at ultra-high fields (UHF)
  • Transcriptome
  • Patient-derived orthotopic xenografts (PDOX)
  • Tumor host interaction