Altered lipid metabolism marks glioblastoma stem and non-stem cells in separate tumor niches

Abstract

Glioblastoma (GBM) displays marked cellular and metabolic heterogeneity that varies among cellular microenvironments within a tumor. Metabolic targeting has long been advocated as a therapy against many tumors including GBM, but how lipid metabolism is altered to suit different microenvironmental conditions and whether cancer stem cells (CSCs) have altered lipid metabolism are outstanding questions in the field. We interrogated gene expression in separate microenvironments of GBM organoid models that mimic the transition between nutrient-rich and nutrient-poor pseudopalisading/perinecrotic tumor zones using spatial-capture RNA-sequencing. We revealed a striking difference in lipid processing gene expression and total lipid content between diverse cell populations from the same patient, with lipid enrichment in hypoxic organoid cores and also in perinecrotic and pseudopalisading regions of primary patient tumors. This was accompanied by regionally restricted upregulation of hypoxia-inducible lipid droplet-associated (HILPDA) gene expression in organoid cores and pseudopalisading regions of clinical GBM specimens, but not lower-grade brain tumors. CSCs have low lipid droplet accumulation compared to non-CSCs in organoid models and xenograft tumors, and prospectively sorted lipid-low GBM cells are functionally enriched for stem cell activity. Targeted lipidomic analysis of multiple patient-derived models revealed a significant shift in lipid metabolism between GBM CSCs and non-CSCs, suggesting that lipid levels may not be simply a product of the microenvironment but also may be a reflection of cellular state. CSCs had decreased levels of major classes of neutral lipids compared to non-CSCs, but had significantly increased polyunsaturated fatty acid production due to high fatty acid desaturase (FADS1/2) expression which was essential to maintain CSC viability and self-renewal. Our data demonstrate spatially and hierarchically distinct lipid metabolism phenotypes occur clinically in the majority of patients, can be recapitulated in laboratory models, and may represent therapeutic targets for GBM.

Introduction

Glioblastoma (GBM) is the most common primary malignant brain tumor in adults. Despite aggressive standard treatment strategies including surgical resection followed by radiation and chemotherapy, the median survival for patients with GBM is approximately 15 months from the time of diagnosis [1]. A key challenge to GBM treatment is the intratumoral heterogeneity at both the cellular and microenvironmental levels [2, 3]. Maintenance of heterogeneity may be driven by a population of cells within the tumor termed cancer stem cells (CSCs), which are highly plastic and responsive to their environment and hold self-renewal and tumor initiation capacity [4, 5]. At the single cell level, GBM is highly heterogeneous with a spectrum of stem cell and metabolic phenotypes [6, 7], and contains both fast-cycling and slow-cycling cells that have distinct metabolisms and cancerous phenotypes [8]. Critically for treatment, this diversity within the cell population means that while cells from one microenvironment or cellular state may respond to a therapy, others may not, resulting in therapeutic resistance of the overall tumor.

Lipid metabolism is abnormally regulated in gliomas compared to normal cells, with changes in the expression of lipid-related genes such as SREBP1 and FAS, which results in altered lipid composition and lipogenesis to keep up with energy demands [9,10,11]. GBM tumors also accumulate more fatty acids than surrounding normal brain tissue [9]. These lipid stores can be used as an energy reservoir [12], can fuel GBM cell proliferation [13], and must be maintained to avoid oxidative damage and lipotoxicity [14]. Recently, lipid metabolism has emerged as a potential therapeutic target to treat gliomas, including GBM [9, 14, 15], and brain metastases [16].

Lipid droplets are cytosolic organelles that, among other functions, serve as a storage medium for the fatty acids, protecting the cell from oxidative damage and providing an energy source to maintain proliferation in unfavorable microenvironments. Lipid droplet formation particularly occurs under stressful conditions such as hypoxia and nutrient deprivation [17]. Accumulation of lipid droplets has been observed in a variety of cancers, including hepatic cancer, lung cancer, breast cancer, and gliomas, and is an important regulator of critical facets of cancer including angiogenesis, inflammatory responses, apoptosis and cell death, and hypoxia-mediated alterations of lipid metabolism [18]. Although lipid droplet accumulation has been observed in hypoxic cells within several tumor types [19, 20], differential lipid metabolism between heterogeneous GBM cell types and microenvironments has remained unexplored. In this study, we analyzed spatial alterations in gene expression and lipid content to determine the metabolic alterations present in GBM CSCs.

Methods

Human cell and organoid culture

Glioblastoma samples were obtained either directly from patients undergoing resection following written informed consent in accordance with protocol #2559 approved by the Cleveland Clinic Institutional Review Board or from collaborators as previously established patient-derived tumorsphere cultures. Patient tissue samples were either finely minced prior to organoid formation or were dissociated into single-cell suspensions, red blood cells were then removed by brief hypotonic lysis, and cells were counted for number and viability using trypan blue. Cells were cultured as tumorspheres in Neurobasal medium supplemented with 10 ng/mL EGF (R&D Systems, Minneapolis MN), 10 ng/mL bFGF (R&D Systems), B27 (Invitrogen, Carlsbad CA), glutamine (CCF media core), sodium pyruvate (Invitrogen), and antibiotics (Antibiotic–Antimycotic, Invitrogen) (“NBM complete”). All cells used in this work were patient-derived primary cultures, and all specimens were verified by comparison of short tandem repeat (STR) analysis performed periodically during the course of experimentation. Tumorspheres were used to form xenografts and harvested for analysis as previously described [21]. All animal experiments were approved by the Cleveland Clinic Institutional Animal Care and Use Committee. Organoids were formed as previously described [22] by suspending tumor cells in 80% Matrigel (BD Biosciences, San Jose, CA) and forming 20 µL pearls on parafilm molds prior to culture. Organoids were cultured in 6-well or 10-cm plates with shaking in NBM complete media. For prospective stem cell sorting, subcutaneous xenografts were minced and digested with papain (Worthington) as previously described [23], and dissociated cells were allowed to recover overnight prior to use. Following overnight recovery from papain digestion, dissociated xenograft GBM cells were magnetically sorted based on CD133 expression using magnetic beads (CD133/2 beads, Miltenyi). This approach has previously been validated to show differences in tumorigenic potential between CD133-positive and CD133-negative fractions [21, 23,24,25,26].

Regional isolation of GBM organoids

Regional labeling and subsequent cell isolation from organoid layers was achieved using a 20 µM final concentration of CellTracker Blue CMAC Dye (Invitrogen, #C2110) in media. Mature organoids were incubated with dye for 2 h at 37 °C with shaking to allow outer layer labeling. Desired labeling depth was verified using a confocal microscope and a compatible imaging dish (MatTek #P35GC-1.0-14-C). After labeling, organoids were finely chopped and dissociated using Accutase (FisherSci, #ICN1000449) at 4 °C for 15 min and then heated to 37 °C for another 10 min. Cells were then single-cell filtered, and live cells were isolated by fluorescence-activated cell sorting (FACS) using 1 µM Calcein-AM (Invitrogen, #C3099MP) and 1:2000 TO-PRO3 (Invitrogen, #T3605) according to the manufacturer’s protocols. Cell sorting and analysis were performed using a BD FACS ARIA II flow cytometer.

RNA-sequencing analysis

Total RNA was extracted using TRIzol reagent (Life Technologies) and purified with phenol–chloroform extraction including the use of PhaseLock tubes (5PRIME). Samples were prepared for RNA-seq according to the manufacturer’s instructions (Illumina) and sequenced using 101-bp paired-end chemistry on a HiSeq-2000 machine in the UTHSCSA Genomic Facility. FASTQ files were trimmed with TrimGalore, which implements Cutadapt [27] and FastQC [28], to remove low quality reads and trim adapters. Reads were aligned to Gencode v29 using Salmon [29] with correction for GC bias, positional bias and sequence-specific bias. The R/Bioconductor package tximport [30] was used to generate TPM values. Low-expressed genes without at least 1 count in 4 samples were excluded from analysis. Comparisons were performed using DESeq2 [31] on raw counts. PCA plots were generated using the ‘prcomp’ function from the R/Bioconductor package ‘stats’ with default arguments. GSEA analysis was performed using mSigDB (Broad Institute).

Gene expression measurement by RT-qPCR

Sorted cells (above) were expanded in permissive media (NBM complete for CSCs or DMEM with 10% FBS (CCF media core) and antibiotics (Antibiotic–Antimycotic, Invitrogen) for non-CSCs). Cell pellets were collected by centrifugation, flash frozen and stored at − 80 °C for later RNA isolation. Total RNA was extracted using TRIzol reagent (Life Technologies) and Phase Lock Gel Heavy tubes (5PRIME). RNA was quantified using Nanodrop 2000 spectrophotometer (Thermo Scientific), and single strand complementary DNA (cDNA) was prepared using SuperScript III Reverse Transcriptase (Invitrogen #18080093) using 1 µg of RNA each. RT-qPCR was performed using TaqMan gene expression assays (Invitrogen #4331182) for FADS1 (Assay ID #Hs00203685_m1) and FADS2 (Assay ID #Hs00927433_m1) in an Applied Biosystems 7500 real-time PCR system, with 18S as the housekeeping gene. Analysis was performed using the ΔΔCt method.

FACS sorting and analysis

Cells were acutely sorted from subcutaneous xenografts as described above and incubated with BODIPY 505/515 (Invitrogen, D3921) or Nile Red (Invitrogen, N1142) for 15 min in the dark. Cells were washed multiple times with sterile PBS and resuspended in Neurobasal medium with 0.5% BSA prior to FACS sorting and analysis using a BD FACSAria II cytometer in the Cleveland Clinic Flow Cytometry core.

Knockdown of FADS1 and FADS2 gene expression

To knockdown gene expression, independent human FADS1 and FADS2 specific shRNA clone sequences were purchased, along with non-targeting control shRNA (Sigma). Recombinant lentiviruses were produced by transfection of 293 T cells plated in 10 cm dishes (4.5 × 10⁶ cells/dish) cultured in DMEM (with 10% FBS and Antibiotic–Antimycotic) medium. The FADS1/2 lentiviral vector (20 µg), psPAX2 (10 µg), and pMD2.G (5 µg) were co-transfected into 293 T cells using calcium phosphate mediated transient transfection. The next morning, cells were fed with fresh NBM complete medium and virus containing media (VCM) were collected twice over 48 h, filtered using 0.22 micron filter and freeze stored at – 80 °C for later use. For infection of GBM CSCs, cells were plated in Geltrex matrix coated 10 cm dishes (1 × 10⁶ cells/dish) overnight in NBM complete medium which was replaced by VCM the next day. After a day of infection, cells were allowed to recover before undergoing antibiotic selection using 2 µg/mL of puromycin in NBM medium. Cells surviving antibiotic selection were then dissociated and collected. Knockdown efficiency of all FADS1 and FADS2 shRNAs were determined by performing RT-qPCR as described above, and 2 shRNAs for each gene with the highest knockdown efficiency were chosen for downstream experiments.

CellTiter-Glo luminescent cell viability assay

Cells collected after FADS1/2 knockdown were resuspended in NBM complete medium, counted and plated in Geltrex matrix coated 96-well plates (2000 cells/well). Cells were then grown in a 37 °C tissue culture incubator in culture media for 5 days. Cell viability was determined using CellTiter-Glo luminescent cell viability assay (Promega, G7570) as per manufacturer’s instructions.

Limiting dilution assay

To determine tumorsphere-propagating potential, live cells were FACS sorted into wells of 96-well plates at concentrations ranging from 1 to 32 cells per well. For cells from FADS1/2 knockdown, cells were counted and plated into 96-well plates at concentrations ranging from 1 to 128 cells per well. Cells were then grown in a 37 °C tissue culture incubator in culture media for 14 days. The presence or absence of spheres in each well was then assessed and analyzed using the ELDA analysis tool (http://bioinf.wehi.edu.au/software/elda/) to calculate stem cell frequencies [32].

Oil Red O histochemistry

For organoid Oil Red O staining, whole organoids were fixed in 4% paraformaldehyde, cryoprotected in 30% sucrose, and snap frozen in OCT using an isopentane bath chilled with dry ice. Tissue sections with a thickness of 10 µm were cut on a cryomicrotome and mounted on glass slides. Oil Red O staining was performed by the Cleveland Clinic Lerner Research Institute imaging core following standard core protocols using commercial control slides. Slides were digitized with a Leica Aperio digital slide scanner (Leica). For primary patient GBM samples, freshly resected GBM tissue specimens from the Department of Neurosurgery, Odense University Hospital, Odense, Denmark, were frozen at – 40 °C using the MCC cryoembedding compound and PrestoCHILL device (Milestone). Tissue sections with a thickness of 8 µm were cut on a cryomicrotome and mounted on glass slides. The slides were left to dry for 15 min at room temperature and fixed for one hour. Following fixation, slides were washed three times with deionized water and incubated with Oil Red O staining solution (Fluka, CI26125, dissolved in 60% triethyl phosphate) for 30 min. After staining, slides were washed with deionized water three times and counterstained with Mayer’s hematoxylin (Merck). Slides were then rinsed with deionized water for 5 min and mounted with a coverslip using Aquatex mounting medium. Finally, slides were digitalized with the NanoZoomer 2.0HT digital image scanner (Hamamatsu, Japan). The use of tissue specimens was approved by the Danish Data Inspection Authority (approval number 16/11065) and the Regional Scientific Ethical Committee of the Region of Southern Denmark (approval number S-20150148).

Targeted lipidomic profiling

Quantification of neutral lipids and glycerophospholipids was conducted as previously described [33, 34]. Briefly, approximately 10 mg of frozen mouse liver was homogenized in 800 mL ice-cold 0.1 N HCl:CH₃OH (1:1) using a tight-fit glass homogenizer (Kimble/Kontes Glass, Vineland, NJ) for ~ 1 min on ice. The suspension was then transferred to cold 1.5 mL microfuge tubes (Laboratory Product Sales, Rochester, NY) and vortexed with 400 mL cold CHCl₃ for 1 min. The extraction proceeded with centrifugation (5 min, 4 °C, 18,000 g) to separate the two phases. The lower organic layer was collected, and the solvent was evaporated. The resulting lipid film was dissolved in 100 mL isopropanol:hexane:100 mmol/L NH₄CO₂H(aq) (58:40:2) (mobile phase A). Quantification of glycerophospholipids was achieved by the use of a liquid chromatography–mass spectrometry technique using synthetic (non-naturally occurring) diacyl and lysophospholipid standards. Typically, 200 ng of each odd-carbon standard was added per 10–20 mg tissue. Glycerophospholipids were analyzed on an Applied Biosystems/MDS SCIEX 4000 Q TRAP hybrid triple quadrupole/linear ion trap mass spectrometer (Applied Biosystems, Foster City, CA) and a Shimadzu high-pressure liquid chromatography system with a Phenomenex Luna Silica column (2, 3, 250 mm, 5 mm particle size) using a gradient elution. The identification of the individual species, achieved by liquid chromatography-tandem mass spectrometry, was based on their chromatographic and mass spectral characteristics. This analysis allows identification of the two fatty acid moieties but does not determine their position on the glycerol backbone (sn-1 vs. sn-2). Triacylglycerol (TAG), diacylglycerol (DAG), and monoacylglycerol (MAG) from frozen mouse liver tissue (10–15 mg) were extracted by homogenizing tissue in the presence of internal standards (500 ng each of 14:0 MAG, 24:0 DAG, and 42:0 TAG) in 2 mL PBS and extracting with 2 mL ethyl acetate: trimethylpentane (25:75). After drying the extracts, the lipid film was dissolved in 1 mL hexane:isopropanol (4:1) and passed through a bed of 60 Å Silica gel to remove the remaining polar phospholipids. Solvent from the collected fractions was evaporated, and lipid film was redissolved in 100 mL CH₃OH:CHCl₃ (9:1) containing 10 mL of 100 mmol/L CH₃COONa for mass spectrometry analysis as described previously [33, 34].

Measurement of de novo lipogenesis flux

Bulk GBM tumors were sorted to enrich populations in cancer stem cells (CD133+) or non-stem (CD133-) populations which were seeded in 35-mm plates. Measurement of de novo lipogenesis rates was accomplished by tracing [14C]-acetate or [3H]-oleate into triacylglycerol and total phospholipids as described previously [33, 35]. Sorted cell populations were simultaneously incubated with 0.5 µCi [14C]-acetate (substrate for de novo fatty acid synthesis) or 1 µCi [3H]-oleate (substrate for direct esterification into complex lipids), and then cells were harvested at various time points (30 min, 60 min, 120 min and 240 min) post substrate addition. From each time point cells were rinsed with PBS (twice), lipids were extracted using a Folch extraction [36], and separated by thin layer chromatography (TLC) using hexane:diethyl ether:acetic acid (70:30:1) as a solvent system. Total phospholipids and Triacylglycerol spots were scraped off of the plate, and the incorporation of [14C]-acetate and [3H]-oleate into each lipid class was determined by liquid scintillation counting. Radiation count was normalized to amount of protein, as quantified by BCA assay (Pierce).

Results

GBM organoids mimic the pathologic transition zones and molecular heterogeneity of GBM patient tumors

GBM tumors have a complex microenvironment defined by histologic hallmarks of angiogenesis and pseudopalisading necrosis, and these two anatomic features harbor distinct GBM cell populations. CSCs are particularly enriched in the perivascular niche of glioblastoma tumors [37]. Current culture models fail to replicate the complex microenvironments of CSCs, limiting our ability to study and therapeutically target GBM. We previously developed 3D patient-derived GBM organoids that show a heterogeneous mix of cells forming zones comparative to the pathologic transition zones in patient tumors [22]. GBM organoids can be conceptually divided into two zones—an outer cell-dense rim consists of dividing cells in a high-oxygen and high-nutrient environment provided by the nearby media, mimicking the conditions of the tumor perivascular niche, and a hypoxic core with necrotic cells relatively deprived of the nutrient media components, which phenotypically mimics the hypoxic and perinecrotic regions of GBM tumors (Fig. 1A).

Fig. 1

GBM organoids mimic the pathologic transition zones and molecular heterogeneity of GBM patient tumors. A H&E staining of GBM 3D organoids (right panel) reveals histological zones comparable to GBM primary patient tumors (left panel). The perivascular region and hypoxic core in primary patient tumors (left panel) are mimicked by the organoid proliferative rim and hypoxic core regions (right panel), respectively. B To compare the molecular signature of these histological regions, the organoids were stained whole to label the entire outer rim region, and single cells were isolated. C Limiting-dilution assays showed that the organoid proliferative rim is functionally enriched for stem cells compared to the hypoxic core. Calculated stem cell frequencies and 95% confidence intervals are shown. D Upon RNA-seq analysis of single cells, distinct cell-type signatures were found to be enriched within spatially separate niches in the organoids. E Mapping expression in organoids to the regional Ivy GAP database showed region-specific enrichment. Scale bar = 100 µm

To determine whether these ex vivo culture regions molecularly mimic the corresponding patient GBM tumor regions, we developed a method to three-dimensionally label and sort live GBM cells from our organoid cultures (Fig. 1B). Whole organoids were incubated with a blue lipophilic dye for an optimized period that we empirically determined to be sufficient to label only the outer proliferative niche as viewed by live confocal z-section imaging. Organoids were then dissociated and sorted via FACS to separate these spatially distinct cell populations. Cells in the outer rim display enhanced self-renewal, a CSC hallmark, compared to cells within the core, as determined by limiting-dilution assay (Fig. 1C).

We performed spatially defined RNA-sequencing of the different regions of GBM organoids to investigate the gene expression in each niche population. Recently, Neftel et al. [7] defined single-cell heterogeneity of patient GBM cell populations as being dominated by clusters of genes that change in relationship to each other, called meta-modules. We compared these gene expression signatures to these meta-module gene signatures from single cell RNA sequencing (scRNA-seq) of clinical GBM. We found distinct cell-type signatures enriched within spatially separate niches in our organoids (Fig. 1D). As anticipated, the cells from the organoid core were enriched for hypoxia hallmark genes as determined by gene set enrichment analysis (GSEA; FDR q value = 3.22 × 10⁻¹⁶) and displayed a signature of the hypoxia-dependent mesenchymal-like 2 (MES2) meta-module defined by scRNA-seq of patient tumor cells. In contrast, cells from the organoid rim were enriched for genes found in the astrocyte-like (AC) and oligodendrocyte precursor cell (OPC)-like meta-modules (Fig. 1D). Additionally, the cells from the organoid rim highly expressed genes related to the G1/S meta-module, indicating high proliferation in this niche. Overall, the expression profiles within our spatially segregated organoid microenvironments represent 3 out of 4 categories of discrete GBM cell types identified in clinical tumors by Neftel et al. [7].

While these data provide resolution to connect the diverse cell types present in GBM organoids with individual cell signatures from clinical tumors, the single-cell sequencing data cannot be directly traced to spatial tumor subregions. To correlate the niche-specific gene expression in organoids and the gene expression of different regions of primary GBM tumors, we reflected our sequencing data upon subregion sequencing data from 41 patient samples in the Ivy Glioblastoma Atlas Project (Ivy GAP) database [38] (Fig. 1E). We found that genes significantly enriched in the organoid core were specifically highly expressed in the pseudopalisading and perinecrotic tumor regions of primary GBM. Conversely, the organoid rim was found to be enriched for genes expressed in the cellular tumor regions and depleted of those expressed in regions of hypoxia. Taken together, these data show that the 3D environments within GBM organoids recapitulate at least part of the cellular and microenvironmental diversity within primary GBM tumors at both a histologic and molecular level.

Genes that are differentially overexpressed in tumor compared to normal brain tissue may suggest a possible therapeutic window for a candidate therapeutic target in brain tumors. To determine genes from our analysis that may be upregulated in tumor compared to normal brain, and upregulated in increasing grades of tumor, we compared the TCGA, Gravendeel, Bao and Ivy GAP datasets using the GlioVis portal [39]. First, to validate our findings, we chose to look at differential expression of specific genes that represent the different regions, particularly highlighting representative genes from the three different meta-modules, including a well-established hypoxia marker, Carbonic Anhydrase 9 (CAIX, Additional file 1: Fig. S1). We found that CAIX and vimentin (VIM) (Mes2 meta-module) were highly expressed in GBM brain tumor compared to non-tumor tissue, due in part to their increase in pseudopalisading regions. This was reflected in GBM organoids where expression of the CAIX and VIM genes was likewise increased in the core region compared to the rim. Although the expression of the OLIG1 (OPC meta-module) and AQP4 (AC meta-module) genes was not significantly increased in GBM compared to non-tumor tissue, these genes are differentially expressed within GBM patient tumor regions and these differences are again recapitulated in the corresponding organoid regions. Taken together, the above findings demonstrate that we can recapitulate tumor-relevant cellular heterogeneity and maintain microenvironmentally regulated GBM cell behavior ex vivo.

Lipid droplets accumulate in the core of GBM organoids and the corresponding perinecrotic and pseudopalisading zones of GBM patient tumors

In addition to the above genes, we identified the gene hypoxia-inducible lipid droplet-associated (HILPDA), which encodes a protein necessary for lipid trafficking in cytosolic lipid droplets, to be consistently differentially expressed in the different regions of 3D GBM organoids and GBM clinical datasets (Fig. 2A). HILPDA expression was significantly higher in brain tumor tissue compared to normal brain and higher in GBM compared to lower-grade brain tumors. Moreover, in primary GBM, HILPDA expression was significantly higher in the pseudopalisading region of the tumor. The increase in HILPDA expression in only GBM and not lower-grade brain tumors, combined with the specific increase in pseudopalisading regions, is consistent with brain tumor pathology as pseudopalisading necrosis is a defining diagnostic feature of GBM. We observed similar differences in the corresponding regions of GBM organoids: HILPDA expression was significantly higher in the core region of organoids compared to the rim (Fig. 2A).

Fig. 2

Lipid droplet accumulation in perinecrotic and pseudopalisading tumor regions and corresponding GBM organoid cores. A Publicly available databases show that HILPDA is consistently increased in GBM tumors and specifically enriched in hypoxic pseudopalisading cells, which is recapitulated by the organoid core. * p < 0.01; ** p < 0.001; ns, p > 0.05. Lipid droplet staining with Oil Red O shows higher staining in the B organoid core and C pseudopalisading and perinecrotic regions of primary tumors. Scale bar for wide field images = 100 µm and 50 µm for other images. D Pie charts representing the findings in 11 patient tumors and 9 organoids

Although we had not initially sought to investigate lipid metabolism, the increase in HILPDA expression, combined with an increase in hallmark adipogenesis genes (p = 2.45 × 10⁻⁵, FDR q-value = 1.36 × 10⁻⁴) and cholesterol homeostasis genes (p = 2.16 × 10⁻⁵, FDR q-value = 1.36 × 10⁻⁴) in the organoid core compared to the rim, prompted us to investigate whether there was any lipid accumulation phenotype that would indicate an overall alteration in lipid metabolism and storage. We therefore stained sections from primary GBM samples and lab-grown GBM organoids with a lysochrome diazo dye, Oil Red O, for histological visualization of lipid droplets. We observed notable differences in lipid droplet staining within the two regions of organoids. Oil Red O staining was concentrated in the cells of the core region of GBM organoids, indicating accumulation of lipid droplets in these cells, whereas the cells in the rim region were devoid of the stain (Fig. 2B, Additional file 1: Fig. S2). To determine whether this ex vivo phenotype is representative of human tumors, we further analyzed multiple primary patient GBM sections. Similar to the GBM organoids, cells in the corresponding pseudopalisading and perinecrotic regions of primary tumors specifically stained for the Oil Red O dye, while the cells in the cellular tumor region lacked the stain (Fig. 2C, Additional file 1: Fig. S3). This staining pattern held true for the vast majority (73%) of samples (Fig. 2D). These results demonstrate that our initial in vitro findings represent a clinically relevant phenomenon of altered lipid metabolism and storage between different regions of GBM tumors.

Differential lipid accumulation marks GBM CSC and non-CSC populations

As the rim region of GBM organoids is enriched for CSCs and the hypoxic core has limited CSCs (Fig. 1C), we investigated whether lipid droplet accumulation is associated with stem cell phenotype. Since GBM organoids are a relatively new technology and can be strongly influenced by the choice of media and culture conditions, we chose to utilize in vivo patient-derived xenograft (PDX) models for this purpose. In vivo tumors are a gold standard for tumor biology and represent perhaps the most realistic and cellularly diverse recapitulation of human tumor microenvironments and stem cell populations. Although stem-like cells within a tumor are controversial and likely exist along a spectrum as opposed to a dichotomous population, prospective sorting into CSC-enriched and -depleted populations requires selection of a surface marker. We utilized the CD133 epitope which is the most longstanding such marker in GBM and has previously been shown effective in enriching CSCs in the cell models used in this study [21, 23,24,25,26]. We collected and dissociated GBM tumor xenografts to obtain a heterogeneous mix of single tumor cells, then magnetically sorted the cells into populations of CD133-positive cancer stem cells (CSCs) and CD133-negative non-stem cancer cells (non-CSCs) (Fig. 3A). Upon flow cytometric analysis, we observed increased fluorescence in the non-CSC populations of all three PDX specimens for both lipid dyes tested (Fig. 3C). Oil Red O dye staining of fixed sorted cells also confirmed this result (Fig. 3B).

Fig. 3

GBM CSCs and non-CSCs have differential levels of lipid accumulation. A Dissociated cells from GBM PDX models were stained for the stem cell marker CD133 and magnetically sorted. B Oil Red O staining was higher in the cultured non-CSC population than in the CSC population. C The CD133-negative non-CSC population showed increased fluorescence compared to CD133-positive CSCs for both Nile red and BODIPY lipid-specific fluorescent dyes. D Magnetically sorted CSC and non-CSC populations were subsequently FACS sorted into lipid-high and lipid-low populations based on BODIPY fluorescence (top and bottom 20%). Sphere-forming capability was determined by limiting-dilution assay. Non-CSCs showed uniformly low sphere-forming capability. BODIPY-low CSCs had significantly higher self-renewal capacity than their BODIPY-high counterparts (p = 0.002). Calculated stem cell frequencies and 95% confidence intervals are shown

To further validate this observation, we asked whether lipid content can enrich for CSCs with increased sphere-forming capability. We sorted dissociated GBM cells for CD133 as above, stained the CSCs and non-CSCs with BODIPY lipid dye, FACS sorted BODIPY-high and BODIPY-low populations, and tested their sphere-forming capability by limiting-dilution assay. Non-CSCs had uniformly low sphere-forming capability, as expected. However, BODIPY-low CSCs were functionally enriched for sphere-forming behavior compared to CSCs with high lipid content (Fig. 3D). We therefore conclude that CSCs are enriched among cells with reduced lipid droplet accumulation and that lipid accumulation may be an indicator of the CSC/non-CSC cell state.

To ask whether the enhanced lipid droplet accumulation seen in non-CSCs was due to increased de novo lipid synthesis by these cells, and to shed light on the separate lipid biosynthetic pathways that may be utilized, we used radiocarbon labeled acetate or tritiated oleate and measured their incorporation into complex lipids over time in either GSCs or non-GSCs. The acetate incorporation is indicative of the de novo lipogenic pathway, while incorporation of tritiated oleate is indicative of direct esterification of the existing labeled fatty acid into a complex lipid. Lipids were extracted and separated into total phospholipids (PL) or triacylglycerols (TG, stored in cytosolic lipid droplets) then analyzed by liquid scintillation counting (Additional file 1: Fig. S4). We found a dramatic increase in radiolabeled phospholipids in GSCs compared to non-GSCs, and this was due to both de novo synthesis and esterification pathways. We did not observe a difference in the contribution of de novo synthesis to TG in either cell population, suggesting that increased de novo fatty acid synthesis in non-CSCs is not driving this phenotype. Taken together, these data show that GSCs prefer to shuttle de novo synthesized (from acetate) or exogenous (oleate) fatty acids into bulk phospholipids and away from triacylglycerols. We did find a notable increase in esterified oleate in non-GSCs at the later timepoint (Additional file 1: Fig. S4D), suggesting that esterification of existing fatty acids play a contributing role to the higher lipid droplet content maintained in these cells.

Lipidomic profiling of CSCs and non-CSCs from patient-derived samples reveals increased neutral lipid species in non-CSC populations

Non-CSCs have higher lipid accumulation in comparison to CSCs, but we cannot resolve individual lipid species by dye. We therefore investigated the differences in lipid metabolism between CSCs and non-CSCs using targeted lipidomic approaches. We isolated cells from 5 patient-derived PDX models, sorted for CD133 ± cell populations, and analyzed these as pools of CD133 + CSCs and CD133- non-CSCs (Fig. 4A) via targeted lipidomic profiling as previously described [33, 34]. We found that the high lipid content in the non-CSC population is due to a broad increase in lipids known to be preferentially stored in cytosolic lipid droplets. Neutral lipid species including diacylglycerol (DAG) and triacylglycerol (TAG) were significantly enriched in the non-CSCs (Fig. 4B, C). These data are consistent with the marked accumulation of cytosolic lipid droplets in non-CSCs.

Fig. 4

Lipid profiling of GBM stem and non-stem cells show differences in global lipid profiling. A Matched pairs (n = 5) of magnetically sorted CD133-positive CSCs and CD133-negative non-CSCs were pooled and analyzed together for different classes of lipids using mass spectroscopy. GBM CSCs have decreased levels of the neutral lipid species B DAG and C TAG. Specific classes of phospholipids are altered in CSCs vs non-CSCs: D Overall phospholipid distribution and E Distribution of minor phospholipid classes. *p < 0.05

GBM CSCs from patient-derived samples exhibit species-specific alterations in glycerophospholipids

Most, but not all (Additional file 1: Fig. S5), molecular species of neutral lipids (DAGs and TAGs) were enriched in non-CSC populations. However, we observed numerous species-specific alterations in glycerophospholipid levels (Figs. 4, 5 and 6). Analysis of total levels of glycerophospholipids revealed that CSCs exhibit modest decreases in minor phospholipid classes including phosphatidic acid (PA), phosphatidylethanolamine (PE), phosphatidylglycerol (PG), phosphatidylinositol (PI) and phosphatidyl (Fig. 4D, E). However, the most abundant class of glycerophospholipids, phosphatidylcholines (PC), was unaltered (Fig. 4D). When we examined molecular species within each class, CSC populations exhibited a general decrease in the levels of longer-chain polyunsaturated fatty acid (PUFA) species of PS (Fig. 5A), PG (Fig. 5B), PtdOH (Fig. 5C), PE (Fig. 5D), PC (Fig. 5E), and LPC (Fig. 5F) lipid classes, indicating that the esterification of PUFAs into these complex lipids may be selectively impaired (summarized in Fig. 5G).

Fig. 5

Lipidomic profiling of GBM stem and non-stem cells shows differences in specific lipid species. GBM CSCs and non-CSCs have differentially expressed phospholipid species, with most of the increase in lipid species in non-CSCs. A Phosphatidylserine, B Phosphatidylglycerol, C Phosphatidic acid, D Phosphatidylethanolamine, E Phosphatidylcholine, and F Lysophosphatidylethanolamine species. G Table highlighting statistically significant differences in lipid groups. * p < 0.05; ns, p > 0.05

Fig. 6

Fatty acid desaturases are upregulated in CSCs and required for CSC survival. A GBM CSCs and non-CSCs have notable differences in phosphatidylinositol species. In particular, CSCs have reduced arachidonic acid (AA) levels but increased levels of its precursor DGLA. B and C qRT-PCR shows that FADS1 and FADS2 levels are higher in CSCs compared to non-CSCs from PDX models. D and E RNA-seq shows higher expression of FADS1 and FADS2 in cells of the GBM organoid proliferative rim and patient cellular tumor regions. CellTiter-Glo cell viability assay shows decreased viability of GBM cells upon FADS1 (F) or FADS2 (G) knockdown. H Proposed mechanism of high accumulation vs high flux of lipids in GBM non-CSCs and CSCs, respectively. * p < 0.01; ** p < 0.001; *** p < 0.0001

Fatty acid desaturases are upregulated in CSCs

Within the PI class, we observed a marked decrease in 38:4 PI and a reciprocal increase in 38:3 PI (Fig. 6A). It is generally accepted that 38:4 PI contains arachidonic acid (AA; 20:4; n-6) in the sn-2 position, while 38:3 PI is expected to harbor di-homo-γ-linolenic acid (DGLA; 20:3; n-6) in the sn-2 position. This reciprocal alteration prompted us to examine the expression levels of the delta-5 desaturase enzyme fatty acid desaturase 1 (FADS1), which converts DGLA to AA. Interestingly, the expression of FADS1 and that of the delta-6 desaturase FADS2 was elevated in CSC populations (Fig. 6B, C). We also found that both FADS1 and FADS2 were significantly upregulated in the organoid rim region and corresponding cellular tumor regions in patient tumors (Fig. 6D, E). A recent report showed that another PUFA-synthesizing enzyme, fatty acid elongase 2 (ELOVL2), exhibits a similar distribution [40]. Increased ELOVL2 and FADS2 further process AA into long-chain PUFAs, thus preventing high AA accumulation in CSCs despite high PUFA pathway activity. Given these findings, there is growing evidence for a role for PUFA synthesis in supporting tumorigenesis within the GBM microenvironment.

FADS1 and FADS2 are essential for GBM CSC survival and maintenance

To functionally test the requirement for FADS1 and FADS2 activity in CSCs, we identified 2 independent shRNA sequences each for FADS1 and FADS2 with validated effective RNA knockdown, achieving approximately 87–97% gene expression knockdown by qPCR (Additional file 1: Fig. S6). We found that upon FADS1 or FADS2 knockdown, GBM CSCs were unable to normally proliferate and survive. Proliferation assays indicated significantly fewer cells in wells containing FADS1 or FADS2 shRNA knockdown cells compared to control shRNA (Fig. 6F, G). Additionally, limiting dilution assays performed to determine the tumorsphere forming capability of GBM CSCs showed that functionally stem-like cell behavior was almost non-existent after FADS1 or FADS2 knockdown (Additional file 1: Fig. S7).

The current evidence suggests that the generation of PUFA-enriched glycerophospholipids appears to be favored in the nutrient-rich CSC microenvironment and supported by high FADS1 and FADS2 expression, whereas lipid accumulation and storage are favored in the nutrient-low non-CSC microenvironment (Fig. 6H). Taken together, our results demonstrate a striking degree of metabolic diversity in GBM depending on each cell’s microenvironment and CSC state, and this heterogeneity must be taken into account in both basic research and in therapeutic targeting of GBM.

Discussion

To meet energy demands in resource-sparse tumor microenvironments, tumor cells undergo metabolic reprogramming, which is known as a hallmark of GBM in addition to many other cancers [41, 42]. Metabolic targeting has been proposed as a therapy for many tumor types, and inhibition of DGAT1 has recently been proposed to alter fat metabolism and increase oxidative stress in GBM [14]. However, cellular heterogeneity and plasticity are features of GBM and drive therapeutic resistance [43]. Recently, intratumoral heterogeneity has become a highly researched focus of both pre-clinical and clinical GBM studies aiming to develop targeted treatment methodologies [44], and it is critically important to mimic this feature through in vitro cultures for more relevant study outcomes. Our overall findings show that lipid droplets accumulate in the hypoxic core of GBM organoids and also in perinecrotic and pseudopalisading regions of GBM patient tumors. This was accompanied by overall increased accumulation of fatty acid species in the CD133-negative non-CSC population versus matched CD133-positive CSCs obtained from patient-derived xenografts. In short, we show that intratumoral lipid metabolism heterogeneity exists and must be considered at the pathologic, cellular and molecular levels.

Pioneering studies by Patel et al. [3], and more recently by Neftel et al. [7], used single-cell RNA sequencing (scRNA-seq) to show that GBM cells vary in their expression of different transcriptional programs, including those associated with proliferation and hypoxia. Consistent with these studies, we found that distinct meta-module signatures from patient tumor cells are enriched within our distinct GBM 3D organoid spatial regions. While cells from the organoid core highly express genes belonging to the mesenchymal-like (MES) meta-module, cells from the rim have higher expression of the astrocyte-like (AC) and oligodendrocyte-progenitor (OPC)-like meta-modules. Although spatial information is lost when tumors are dissociated for scRNA-seq, it is assumed that these different populations from patient tumors derive from different microenvironments within the tumor. Here, we further showed that the gene expression in organoids corresponds to the gene expression in the regional Ivy GAP database [38]. The genes enriched in the organoid core reflect the patient pseudopalisading and perinecrotic tumor regions in the Ivy GAP data, whereas the genes enriched in the organoid rim were associated with the gene enriched in the cellular tumor region. Thus, we successfully ascertained that molecular heterogeneity in the 3D organoid culture corresponds to primary GBM tumors at both the single-cell and spatial microenvironmental levels. This spatial heterogeneity is especially important to appreciate when looking at datasets derived from small patient samples as the sample site selection may have a profound influence upon the obtained data (for instance, whether a pseudopallisade was present within the sample site when investigating HILPDA expression).

We show that HILPDA expression, which is canonically driven by HIF-1α, is consistently increased in primary GBM tumors and specifically in the pseudopalisading/perinecrotic region of GBM (Fig. 2A). This supports a finding by Mao et al. showing that HILPDA was upregulated in GBM compared to normal brain tissue or lower-grade gliomas and in GBM cells cultured in hypoxic conditions [45]. Additionally, studies show that HILPDA is involved in triglyceride fatty acid secretion [46] and regulates lipid metabolism and hypoxia-induced lipid droplet biogenesis [47]. When we investigated lipid droplet accumulation in our 3D GBM organoid model, we discovered that lipid droplets were exclusively enriched in the hypoxic core region of the 3D organoids. We further traced this finding back to patients, showing that all patients in our panel with clear pseudopalisading necrosis by pathology, and most patients overall, accumulate lipid droplets in the hypoxic regions of their tumors. To our knowledge, this is the first concrete demonstration of lipid droplet accumulation as a marker of cells surrounding pseudopalisading necrosis.

It is an outstanding open question in the field [42] whether there is a link between self-renewal and lipid metabolism in GBM. Here, we addressed this question using GBM patient tissue, fresh patient-derived 3D ex vivo cultures, and patient-derived xenografts. We obtained consistent findings at the single molecular, transcriptional, cellular, and tissue scales, all of which show increased lipid content in nutrient-poor and non-stem GBM cells. Our results support findings showing that slower-cycling GBM cells (as found in GBM organoid cores (Fig. 1D and previous work [22])) have increased lipid droplet content [8] and that targeting lipid homeostasis in GBM has antiproliferative effects [14]. However this also contrasts with findings that slow-cycling tumor cells may be CSCs [8] and that cultured colorectal cancer CSCs have increased lipid content compared to non-CSC populations [48]. Hypoxia is indeed known to promote the GBM stem cell phenotype [49,50,51] and we observe rare functionally verified CSCs within organoid cores (Fig. 1C). These hypoxic CSCs may have a different lipid profile, or at least a different place within the tumor cell spectrum, than hypoxic non-CSCs or non-hypoxic CSCs and warrant individualized studies.

Along with lipogenesis, lipolysis (through oxidation of fatty acids) is critical for the renewal of stem cells and maintenance of stemness [52], which could explain the decreased lipid accumulation and overall decrease of lipids in CSCs despite high lipid synthesis rates (Fig. 3D). One initial paradox in our data was the increased 20:4 AA species in non-CSCs despite lower DGLA (20:3) species compared to CSCs. Upon examining the expression levels of the fatty acid desaturase genes involved in AA/DGLA catabolism, FADS1 and FADS2, we found that both are enriched in the CSC-rich organoid rim. Additionally, ELOVL2, a critical enzyme downstream of AA, has been shown to be enriched in SOX2- and OLIG2-positive GBM cells [40]. We propose that ELOVL2 and FADS2 in CSCs facilitate turnover/utilization of lipid species in high nutrient conditions, unlike in non-CSCs and GBM cells in nutrient-poor environments, where the fatty acids accumulate, forming lipid droplets. We believe that this potential mechanism makes intuitive sense from a standpoint of survival, where energy storage is favored in environments where nutrient resources are scarce and uncertain. Also we have demonstrated that CSCs have a high dependence on FADS1 and FADS2 function to proliferate and maintain self-renewal. There are multiple commercial or published FADS1/2 or PUFA inhibitors including: CP 24,879 [53, 54]; 8,11-Eicosadiynoic Acid [55]; and compound-326 [56, 57]. Our work suggests that targeting PUFA synthesis in GBM with these or other improved FADS inhibitors may have therapeutic benefit in GBM.