Zongze He1†
Zheng Liu2†Qi Wang1†
Xingjian Sima3†Wei Zhao4†Chunmei He5Wenjie Yang1Han Chen1Bo Gong6,7*
Siyuan Song8*
Yi Wang4,9*- 1Department of Neurosurgery, Sichuan Provincial People’s Hospital, University of Electronic Science and Technology of China, Chengdu, China
- 2Department of Medicine, Baylor College of Medicine, Houston, TX, United States
- 3Medical School, Tongji Medical College of Huazhong University of Science and Technology, Wuhan, China
- 4Center of Critical Care Medicine, Sichuan Provincial People’s Hospital, School of Medicine, University of Electronic Science and Technology of China, Chengdu, China
- 5Department of Otolaryngology, Chongqing General Hospital of the Chinese People’s Armed Police Force, Chongqing, China
- 6Department of Health Management, Sichuan Provincial People’s Hospital, University of Electronic Science and Technology of China, Chengdu, Sichuan, China
- 7The Key Laboratory for Human Disease Gene Study of Sichuan Province and Institute of Laboratory Medicine, Sichuan Provincial People’s Hospital, University of Electronic Science and Technology of China, Chengdu, Sichuan, China
- 8Department of Neuroscience, Baylor College of Medicine, Houston, TX, United States
- 9Clinical Immunology Translational Medicine Key Laboratory of Sichuan Province, Sichuan Academy of Medical Science and Sichuan Provincial People’s Hospital, Chengdu, China
Background: Glioma is a highly heterogeneous malignancy of the central nervous system. This heterogeneity is driven by various molecular processes, including neoplastic transformation, cell cycle dysregulation, and angiogenesis. Among these biomolecular events, inflammation and stress pathways in the development and driving factors of glioma heterogeneity have been reported. However, the mechanisms of glioma heterogeneity under stress response remain unclear, especially from a spatial aspect.
Methods: This study employed single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (ST) to explore the impact of oxidative stress response genes in oligodendrocyte precursor cells (OPCs). Our analysis identified distinct pathways activated by oxidative stress in two different types of gliomas: high- and low- grade (HG and LG) gliomas.
Results: In HG gliomas, oxidative stress induced a metabolic shift from oxidative phosphorylation to glycolysis, promoting cell survival by preventing apoptosis. This metabolic reprogramming was accompanied by epithelial-to-mesenchymal transition (EMT) and an upregulation of stress response genes. Furthermore, SCENIC (Single-Cell rEgulatory Network Inference and Clustering) analysis revealed that oxidative stress activated the AP1 transcription factor in HG gliomas, thereby enhancing tumor cell survival and proliferation.
Conclusion: Our findings provide a novel perspective on the mechanisms of oxidative stress responses across various grades of gliomas. This insight enhances our comprehension of the evolutionary processes and heterogeneity within gliomas, potentially guiding future research and therapeutic strategies.
Introduction
Gliomas are the most frequent form of central nervous system (CNS) tumor (1). China is one of the top three countries in terms of CNS cancer cases (including glioma) and fatalities (2, 3). In the United States, the average incidence of malignant brain tumors was 7.08 per 100,000 people, with glioblastoma (GBM, WHO IV high-grade glioma) accounting for around 15% of all CNS tumors (4). Despite advances in surgery, chemotherapy, and radiation, these treatments have not significantly improved the overall survival rates for glioma patients, particularly those with high-grade gliomas such as GBM (5, 6). The blood-brain barrier poses a significant challenge by preventing most antitumor drugs from reaching the brain, thus limiting the effectiveness of systemic therapies (7). Therefore, investigating novel mechanisms and pathway alternations in gliomagenesis may improve diagnostic accuracy, prognosis and overall treatment outcomes for glioma patients.
Glioblastomas are the most malignant and aggressive type of gliomas, characterized with high heterogeneity and mortality. Researchers have showed that gliomas have various subtypes expressing different marker genes including typical NPC-like, OPC-like, AC-like and MES-like profiles using multiple datasets (8). Stress management is critical for cancer cell survival and evolution. Glioma cells generate more reactive oxygen species (ROS) than normal cells due to oncogenic activation, increased metabolic activity and mitochondrial dysfunction, driving heterogeneity. On the other hand, cancer cells may adapt to oxidative stress through mechanisms such as boosting antioxidant status to promotes ROS-driven proliferation while avoiding ROS levels that would cause senescence, apoptosis, or ferroptosis (9–13). Such adaptability leads to malignant transformation, metastasis, and resistance to anticancer treatments (14–17). It is crucial to examine how glioma cells adapt to various stress during the heterogeneity process.
Glioma development and heterogeneity is characterized by gene deregulation through genetic and epigenetic events, as well as cellular damage caused by ROS production (18). Glioma cells maintain metabolic balance by altering the expression of key genes via redox homeostasis pathways such as mitochondrial respiration, stabilizing pro-apoptotic proteins for survival, and preventing apoptosis when challenged. For example, in GBM stem cells (GSCs), activation of AMPK (AMP-activated protein kinase), a master regulator of metabolism, increases glioblastoma bioenergetics and tumor development (19). Another study showed that mitochondrial PKM2 (Pyruvate Kinase M2) regulates oxidative stress-induced apoptosis by stabilizing Bcl2 pro-apoptotic protein (20). Furthermore, the stemness of glioma stem cells is maintained through the interaction of TRAP1 and Sirtuin-3, which modulates mitochondrial respiration and oxidative stress (21). Studies have also found that PTPN2 induced by the inflammatory response and oxidative stress contributed to glioma progression (22). JNK activation increases intracellular ROS production, leading to oxidative stress-induced glioma cell parthanatos (23). Oncostatin M (OSMR) increases mitochondrial respiration by interacting with complex I’s NADH ubiquinone oxidoreductase 1/2 (NDUFS1/2). OSMR deletion reduces spare respiratory capacity, increases ROS, and makes glioma stem cells more susceptible to IR-induced cell death (24). 2-Hydroxyglutarate inhibits transaminase, impairing glutamate biosynthesis and redox homeostasis in glioma (25). Heat shock proteins also protect glioma cells from apoptosis through the endoplasmic reticulum (ER) stress and unfolded protein response pathways. For instance, the Hsp70 chaperone rescues C6 rat glioblastoma cells from oxidative stress by sequestering aggregating GAPDH (26). Recent single-cell RNA-sequencing (scRNA-seq) studies have explored the heterogeneity of gliomas. highlighting the importance of stress responses in heterogeneity and subsequent malignant transformation (27, 28). However, a comprehensive understanding of oxidative stress and its consequences across various gliomas subtypes, particularly at the whole-tissue level, is lacking.
Recent advances in scRNA-seq and spatial transcriptomics (ST) provided ways to visualize and study gene expression in tissue slices, enabling the exploration of transcriptional activity at the single-cell or spatial level (29, 30). scRNA-seq studies have revealed the cellular plasticity and environmental stress response in gliomas (31, 32). Additionally, scRNA-seq has been used to develop comprehensive prognostic signature for glioblastoma patients (33). ST has been applied to various of cancers, including prostate cancer (34), breast cancer (35), and pancreatic ductal adenocarcinomas (36). Our study combines scRNA-seq and ST data from four primary glioma tissue samples to evaluate oxidative stress and its effects. This study aims to elucidate regional and transcriptome-wide expression patterns, enabling that will enable us to deconvolve molecular processes in the oxidative stress response and pathway switching.
Materials and methods
Patient recruitment and sample collections
Patients were recruited at the Sichuan Provincial People’s Hospital Neurosurgery Clinic from July 2019 to September 2020. Detailed verbal and written information about the study was provided to all participants, who gave written informed consent before participating. from the study included two untreated patients: one female with WHO grade IV glioblastoma in the left temporal lobe (high-grade glioma, HG) and one male with WHO grade II oligodendrocyte astrocytoma in the right temporal lobe (low-grade glioma, LG). Detailed basic clinical and pathological information for the participants is provided in Supplementary Table 1.
Single-cell transcriptomic profiling of glioma samples
Four fresh samples were collected prior to surgery and processed into single-cell suspensions. These samples were analyzed using droplet-based single-cell RNA sequencing (scRNA-seq) technology of the 10 x Genomics Chromium system. The Cell Ranger software (version 3.1.0) was employed for initial data processing, including alignment and quality control. Post-quality control, the number of high-quality cells in each sample ranged from 2,419 to 10,733. Low-quality cells (doublets, multiplets, and apoptotic cells) were excluded, resulting in final cell counts ranging from 1,606 to 9,744. The average number of Unique molecular identifiers (UMIs) per cell ranged from 7,325 to 15,424, the average number of genes per cell is 2,424 to 3,858, and the mitochondrial genes ratio per cell from 5.69% to 11.90%. We used the R package Seurat (version 3.1.1) to process the filtered unique molecular identifier (UMI) count matrix (33), with variable gene identification using the algorithm provided by Macosko et al. (37) Principal component analysis (PCA) was used in Seurat to reduce dimensionality using the RunPCA function (PC num = 15) (33). Clusters were analyzed using FindClusters function and visualized via t-distributed stochastic neighbor embedding (t-SNE) (RunTSNE). Marker genes for each cluster were identified using the FindAllMarkers function in Seurat (33). Then, cell types were inferred using the R package SingleR, an automated annotation method for unbiased scRNA-seq cell type detection, with Human Primary Cell Atlas from Mabbott et al. as the reference transcriptome datasets, to infer the cell types (38, 39). Differentially expressed genes (DEGs) were identified with FindMarkers function in Seurat. The criterion for substantially different expression was established at P Value < 0.05 and |log2foldchange| > 0.58. Gene Ontology (GO) and KEGG Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on of DEGs using hypergeometric distribution.
10x Visium spatial RNA-seq data preprocessing
Spatial transcriptomic analysis was conducted using the 10x Genomics Visium platform. The Space Ranger software (version 1.0.0) was used for tissue identification, read alignment via the STAR aligner, feature-spot matrix construction, clustering, gene expression analysis, and spatial visualization slide images (40). This analysis contains the spatial transcriptome sequencing of 4 samples with 4 slices. The number of high-quality spots for Space Ranger quantitative quality control of each sample ranges from 2,791 to 3,821. The average UMI counts per spot ranges from 7,113 to 22,143, and the average number of genes per spot from 3,093 to 6,002, and mitochondrial genes ratios per spot ranges from 8.24% to 13.21%. To process the filtered unique molecular identifier (UMI) count matrix, the R package Seurat (version 3.1.1) was utilized (33). We initially standardized the data to account for differences in sequencing depth across data points before identifying high-variance features and storing the data in the SCT assay (41). FindClusters function was applied to analyze cell groups based on their gene expression profiles using graph-based clustering in Seurat. We used the RunTSNE function to display clusters using a 2-dimensional t-distributed stochastic neighbor embedding (t-SNE) method and the FindAllMarkers function to find marker genes in each cluster (33).
Copy number variation analysis
Genomic copy number and subcolonal structure of human tumors were predicted using the CopyKat and Infercnv packages, enabling high-throughput genomic analysis from scRNA-seq data (42–44).
Multimodal intersection analysis
Multimodal intersection analysis (MIA) was used to annotate spatial transcriptome spot clusters with single-cell transcriptome data, enhancing cell type resolution at the spot cluster level.
Histopathological examination
The collected glioma tissues were paraffin-embedded, sectioned, and stained with hematoxylin and eosin (H&E) for morphological analysis. For immunohistochemistry (IHC), the sections were routinely incubated with the c-Jun antibody (10586-1-AP from Proteintech, 1:200), and processed using a rabbit IgG-immunohistochemical SABC kit (Boster, Wuhan, China). After PBS washing, sections were counterstained with hematoxylin for nuclei staining and examined under an Olympus microscope (Japan).
Results
ScRNA-seg results reveal the links between heterogeneity and malignancy in gliomas
To explore the cell compositions in glioma tissues, scRNA-seq was conducted on all living cells isolated from four fresh primary glioma samples from two untreated patients with the 10x Genomics Chromium platform for parallel scRNA-seq and ST analysis (Figure 1A and method). The scRNA-seq data consisted of cells with approximately 7,325–15,424 unique molecular identifiers (UMIs) and approximately 2,424–3,858 uniquely expressed genes per cell (Supplementary Figure 1). The average percent of mitochondrial, ribosome, and hemoglobin subunit genes were also calculated and used for further data normalization and batch effects between two samples removed using harmony package (45) (Supplementary Figure 1). We found 16 cell clusters in this integrated scRNA-seq data (Figure 1B). We identified them as endothelial cells (ECs), natural killer T cells (NKTs), mural cells, microglia, oligodendrocytes, and malignant cells utilizing cell-type-specific marker genes (Figure 1C). The visualization plot was generated using the dimensionality reduction method UMAP (uniform manifold approximation and projection). The heatmap of the top 30 genes from each cluster is shown in Figure 1D. Moreover, the tSNE (t-distributed Stochastic Neighbor Embedding plot) for the expression of marker genes from each cell type is shown in Supplementary Figure 2A. Two immune-cell clusters (clusters 4, 11) were initially identified based on the expression level of PTPRC (CD45), a pan leukocyte marker (Supplementary Figures 2B; 3A). According to the presence of the typical marker CSF1R (Colony Stimulating Factor 1 Receptor, CD115), cluster 4 was mononuclear phagocytes (microglia) (Supplementary Figure 2D). Cluster 11 contains Natural Killer/T cells (NK/T), which express typical markers such as CD8A, CD8B, NKG7, CCL4, CCL5, but not CCR7, CD4, CTLA4, FOXP3, and IL2RA (Supplementary Figures 2C; 3A). In our glioma samples, this observation suggested the existence of CD8+ T cells rather than naive or CD4+ T cells. Cluster 9 has been identified as endothelial cells expressing specific markers such as FLT1, ICAM1, CDH5, CD34, and ENG (CD105) (Supplementary Figures 2E; 3B). Mural cells, as an integral part of the neurovascular unit, make up the majority of cluster 16, expressing specific markers such as ACTA2, PDGFRB, ZIC1, ABCC9, RGS5, CSPG4, and MCAM (CD146) (Supplementary Figures 2F; 3C). Clusters 7 and 8 are mature oligodendrocytes that produce myelinating oligodendrocyte markers, including SOX10, CNP, MAG, MOG, MBP, and others (Supplementary Figures 2J; 3D). However, they do not express proliferative factors such as EGFR, PDGFRA, or L1CAM. Cluster 1 and 12 expresses GFAP, an astrocyte marker, and are therefore classified as astrocytes. All other clusters express OLIG2, OLIG2, PDGFRA, and NG2 (CSPG4), indicating that they are OPC groups. In another cell cluster, we also discovered the high expression of well-known cancer stem cell markers SOX2 (Figure 1E) and EGFR (Epidermal Growth Factor Receptor), indicating abnormal proliferation and malignancy in these cell clusters (Figure 1F).
Figure 1. Isolation and cluster analysis of single-cell transcriptomes from four glioma samples derived from two patients. (A) Workflow of experimental strategy: isolation of human glioma tissues during clinical surgery, followed by parallel scRNA-seq and Spatial Transcriptomic (ST) analysis profiling using with the 10x Genomics Chromium platform and subsequent validation by H&E and IHC staining. (B) Uniform Manifold Approximation and Projection (UMAP) plot showing 16 major clusters. (C) Clusters are annotated for their cell types as predicted using canonical markers and signature-based annotation using Garnett. (D) Heatmap showing clustering with top 30 highly expressed genes. (E, F) Feature UMAP plots depicting cluster-specific expression of cell clusters markers including SOX2 (SRY-Box Transcription Factor 2) and EGFR (Epidermal growth factor receptor) to indicate the major malignant cell clusters.
Malignant cell heterogeneity and copy number variation analysis
Previous research suggested that four cellular states, NPC-like, OPC-like, AC-like, and MES-like, drive glioblastoma malignant cell heterogeneity. These malignant cells were identified OPC-like cell expression markers like OLIG1 OLIG1, OLIG2, PDGFRA, and NG2 (CSPG4), indicating that they are OPC groups. Infercnv package in R to further validate the malignant cells to calculate the copy number variations (CNVs) (42, 46) (Supplementary Figure 5). The malignant cell group showed a great heterogeneity. As we expected, the OPC-likes have the most significant alternation in CNV levels (Supplementary Figures 5B–D). According to the heatmap, we observed amplification of CNV of chromosome 7/10, a molecular marker of glioma, while a loss of chromosome 11 from the predicted results (Supplementary Figure 5A).
Spatial transcriptomics reveals distinct tumor microenvironments in high-grade and low-grade gliomas
For ST data, we performed the integrated clustering analysis of ST and revealed 9 different spatial clusters (Figure 2A). A significant variation in spatial spot type was observed between HG and LG to explain the heterogeneity. The spatial characteristics and cluster distribution of HG and LG are distinct (Figure 2B). Spatial cluster1 and cluster5 are primarily spread on the HG2 edge and the core portions of LG1 and LG2. While spatial cluster 1 was considered as mixed cell types, spatial cluster 5 was labeled as oligodendrocytes due to its MAG, MOG, GFAP, and MBP expressions. It was mostly found in LG2, although it was also discovered in the HG2 marginal region. Spatial cluster 8 was found in both HG1 and HG2, and it was apparently dispersed around the core tumor region (Figure 2B). The heatmap for the top 30 genes from each ST cluster was generated (Figure 2C). Interestingly, the extracellular matrix protein-encoding genes such as FN1, TIMP1, COL3A1, COL4A1, and COL2A1 were discovered in spatial cluster 8, regarded as tumor stroma tissues. Endothelial markers such as PECAM1, CD93, and ENG were also strongly expressed in spatial cluster 8. Cluster 4 was considered tumor stroma tissues mainly represented by fibroblasts. The correlations between ST clusters and scRNA-seq cell types were shown by integrative analysis of cell data and spatial transcriptome data using Multimodal intersection analysis (MIA). As shown in the heatmap of Figure 2D, clusters 2 and 4 are most likely to be malignant OPC-likes, while clusters 5 and 6 could be oligodendrocytes or astrocytes. Clusters 3 and 8 might be a mixture of immune cells, endothelial cells, and mural cells (Figure 2D).
Figure 2. Cluster analysis and annotations of spatial transcriptomes from four glioma samples derived from two patients. (A) t-distributed stochastic neighbor embedding (tSNE) projection of spots showing 9 major ST clusters. (B) The spatial relationship among 9 major ST clusters in four different samples. (C) Heatmap showing clustering with top 30 highly expressed genes from each ST clusters. (D) Heatmap showing the correlation matrix generated from 8 major ST Clusters and 7 major cell types from scRNA-seq calculated through MIA (Multimodal intersection analysis).
We next annotated the ST spot using SingleR and examined each sample individually to validate the results, considering tumor heterogeneity and individual differences among patients. The most frequent types of malignant cells in HG gliomas are MES-like (orange), OPC-like (yellow), and unidentified malignant type (grass green) (Supplementary Figures 4A, B). Oligodendrocytes (red) and AC-like malignant cells (dark blue) were found near the periphery of HG2, which can be differentiated from the primary malignant cells in HG1 (Supplementary Figure 4B). Normal immune cell infiltration also differs significantly between HG and LG gliomas. For example, monocyte infiltration (light purple) is more clearly observed in HG glioma tissue (HG1 and HG2), likely due to the looser structure of high-grade tumors. In contrast, LG gliomas (LG1 and LG2) have a tighter arrangement of tumor cells, and limiting immune cell infiltration. The great majority of malignant cells in LG glioma spatial spots are categorized as OPC-like (yellow), unidentified (grass green), or AC-like (dark blue), with NPC-like (purple) and MES-like (orange) malignant cells being uncommon (Supplementary Figure 4C). OPC-like-malignant is primarily found on the interior of the sample, whereas AC-like-malignant is more prevalent on the exterior. Simultaneously, LG2 contains more OPC-like-malignant components, while LG1 has more AC-like-malignant components (Supplementary Figure 4D).
OPC clusters revealed activation of stress response genes
We performed a PCA dimensionality reduction analysis on the OPC subpopulation single-cell data and used UMAP projection to obtain the expression matrix (Figure 3A). This analysis revealed 16 subclusters of OPC-like cells (Figure 3A). Using the CopyKat package, we assessed malignant and normal cells based on predicted aneuploid and diploid states to detect aberrant cells in single-cell clusters. The copy number variation (CNV) data show differences in OPC-like subclusters, with normal diploid cells primarily identified in subclusters 3, 4, 8, 10, 11, and aneuploid cells mainly discovered in other subclusters (Figures 3B, C). Both LG and HG samples shared Clusters 3 and 8 (Figure 3D). A heatmap displaying the top 30 genes from each OPC-like subcluster is shown in Figure 3E.
Figure 3. Single-cell RNA-seq of OPC cell populations. The OPC subclusters were extracted from the whole scRNA-seq data and processed for dimensionality reduction. (A) Uniform Manifold Approximation and Projection (UMAP) plot showing 16 major clusters of 7221 OPC single cells. (B-D) UMAP plot with CNV (copy number variation) annotations according to CopyKat package calculations. Aneuploid was depicted as red while diploid was shown as blue. Stacked bar charts showing the distribution of 16 subclusters of OPC with its sample origin and predicted CNV status. (E) Heatmap showing clustering with top 30 highly expressed genes from each OPC subclusters.
Using Monocle3, we set the start points of the pseudotime intervals to OPC-like subclusters 3 and 8 (Figure 4A). According to Monocle3 data, three primary lineages emerge from the starting point (Figures 4C, D). We created 28 distinct modules with varying gene expression patterns (Figure 4B) and found that several stress response genes exhibit distinct expression patterns. Oxidative stress-related genes with specific expression patterns include BTG1, EGR1, ETV1, FOS, HSPA1B, IER2, JUNB, JUN, DNAJB1, and HSPH1, all showing very similar trajectories (Figure 4E). Supplementary Figure 6 demonstrates the elevated levels of oxidative stress response genes and the elevation in pseudotime in an OPC-like subcluster UMAP plot, including stress response transcription regulators (BTG1, EGR1, FOS, IER2, JUNB, JUN) (Supplementary Figure 6A) and heat shock proteins (HSPH1, HSPA1B), as well as DNA damage response gene DNAJB1 (Supplementary Figure 6B). Additionally, tumor proliferation regulator and metabolic genes implicated in glycolysis and gluconeogenesis pathways, including VEGFA, CCN1, VIM, LDHA, LDHB, PKM, ENO1, and PGK1, were up-regulated in HG OPC-like subclusters (data not shown), indicating an increase in stress response during the gliomagenesis.
Figure 4. Pseudo-time analysis of OPC subclusters using Monocle3 (A) Monocle3 was used to track cells over pseudotime on the different OPC subclusters. (B) The heatmap of 39 modules calculating the gene expressions after prediction of trajectory. (C, D) Visualization of clusters (from Figure 1C) onto the pseudotime map. (E) Pseudo-time plots of gene expression for stress response genes (BTG1, EGR1, ETV1, FOS, HSPA1B, HSPH1, IER2, JUN and JUNB) with cells colored by OPC subcluster group indicating the similar patterns among those genes.
Specific expression patterns of stress response genes near the necrotic core
We validated the expression and cell-type distribution of these stress response transcription factors (TFs) using the Ivy Glioblastoma Atlas Project (Supplementary Figure 7). The peri-necrotic zone and pseudo-palisading cells around necrosis showed significantly higher mRNA levels of JUN, JUNB, FOS, and FOSL2 compared to other regions of the tumor, such as the leading-edge, infiltrating tumor, and cellular tumor (Supplementary Figures 7A-D). Pseudo-palisading cells around central degeneration have been recognized for nearly a century as a hallmark of GBM and indicate aggressive tumor behavior. The enrichment of FOSL2 protein in the peri-necrotic zone and pseudopalisading cells area was further confirmed by in situ hybridization (ISH) and H&E staining with tumor feature annotations (Supplementary Figures 7E, F).
Additionally, H&E staining of primary glioblastoma tissues from three patients (A020-21621, A021-12863, A021-15304) during surgery illustrated the representative necrotic areas and pseudo-palisading regions (Supplementary Figure 8A). The elevation of JUN protein expression around these areas was confirmed by immunohistochemical staining (IHC) using the specific antibody against JUN (Supplementary Figure 8B). Data from GEPIA2 also showed frequently increased FOS, JUNB, JUN, and ATF3 in glioma, especially GBM samples, compared to normal brain tissues (Supplementary Figure 9). These findings suggest that JUN family members are important regulators in glioma stress response, especially in GBM.
Oxidative stress response in high-grade glioma
For both HG1 and HG2, we used the ST transcriptome data and an enrichment study in R studio using the SPATA package (version 1.0.1). We discovered that the hypoxia response (Hall-Mark hypoxia pathway) is enriched in twisted bands or spots around the central areas (Figures 5A; 6A). Following that, we utilized the Monocle3 package (version 1.0.1) to determine the trajectory of gene alterations in the process of transitioning from a low hypoxic response to a high hypoxic response in both HG1 and HG2. Dimensionality reduction plots indicate cell population clustering in UMAP and t-SNE for HG1 and HG2 (Figures 5B, C; 6B, C).
Figure 5. The crucial pathway altered in the trajectory under an increasing hypoxia stress in HG1 (A) The trajectory in the HG1 with a increasing hypoxia score. (B, C) tSNE and UMAP plot showing distribution of different ST clusters. (D-H) The significantly altered pathways either descending or ascending in the trajectory. (I) The expressions of marker genes (DNAJB1, HSPA1A, HSPH1, JUNB and JUND) in the trajectory that response to hypoxia stress.
Figure 6. The crucial pathway altered in the trajectory under an increasing hypoxia stress in HG2 (A) The trajectory in the HG2 with an increasing hypoxia score. (B, C) tSNE and UMAP plot showing distribution of different ST clusters. (D-H) The significantly altered pathways either descending or ascending in the trajectory. (I) The expressions of marker genes (JUN, MT1X, PKM and STMN2) in the trajectory that response to hypoxia stress.
As predicted, the hypoxia response increases substantially throughout this phase (Figure 5F). In HG1, increasing routes include the initial immune response pathway (BP: GO) (Figure 5D) and the epithelial-mesenchymal transition (Hallmark) (Figure 5E), while decreasing pathways include glial cell differentiation (BP: GO) (Figure 5G) and nerve cell maturation (BP: GO) (Figure 5H). This suggests that the tumor cells advanced from a somewhat benign condition to a dedifferentiated or even invasive status in the central high-grade tumor sample. The rapid elevation of oxidative stress response gene expression, including JUNB, JUN, DNAJB1, HSPA1A, and HSPH1, was shown under oxidative stress (Figure 5I). On the other hand, hypoxia has a greater impact on cell cycle regulation (Figures 6D, F), DNA replication (Figures 6E, F), and cell proliferation in HG2. We also found that immune responses, such as the response to bacteria (Figures 6F, G) and interferon-gamma signaling (Figures 6F, H), are highly activated in the hypoxia regions in HG2, located in the marginal areas of high-grade glioma. This might be triggered by external stimuli in the brain.
Thus, it is clear that during gliomagenesis, cells first undergo uncontrolled proliferation and then respond to stress by abandoning the essential features of glial cells, leading to malignant transformation. Consistent with these findings, we discovered that tumor-proliferative factors such as VEGFA (Vascular Endothelial Growth Factor A), the representative downstream regulator of epithelial-mesenchymal transition, VIM (Vimentin), and the inflammatory cytokine SPP1, which regulates endoplasmic reticulum stress, are all significantly up-regulated in peripheral hypoxic regions in both HG1 and HG2 samples (Supplementary Figures 12, 14). The subcluster of OPC-like was projected using the ST data; cluster 12 was regarded as malignant, and cluster 3 was less malignant (Supplementary Figures 13, 15). Meanwhile, genes that have functions in the normal brain or neurons, such as CaMK2N1 (Calcium/Calmodulin Dependent Protein Kinase II Inhibitor 1), MEG3 (Maternally Expressed 3), SNCB (Synuclein Beta), SYT1 (Synaptotagmin 1), DKK3 (Dickkopf WNT Signaling Pathway Inhibitor 3), are almost not expressed in the tumor core. Additionally, stress-response transcription factors (JUNB, FOSL2, IER2, EGR1), DNA damage response proteins (DNAJB1), and heat shock proteins (HSPA1A) are dramatically up-regulated in the hypoxia regions (Figures 7, 8). They are concentrated in certain punctate areas, which overlap with hypoxia response areas in HG1 (Figure 7). However, on the other hand, those genes in HG2 are not as active as they are in HG1 regions. While the hypoxia score is growing in the trajectory, there is only a rapid increase in JUNB expression in HG2 (Figure 8).
Figure 7. Expression patterns of stress response genes in HG1 spatial transcriptomics data. (A-F) The expressions of stress response gene including JUN, DNAJB1, FOSL2, HSPA1A, IER2 and EGR1 was shown in the ST profile from HG1 tissues.
Figure 8. Expression patterns of stress response genes in HG2 spatial transcriptomics data. (A-F) The expressions of stress response genes including JUN, DNAJB1, FOSL2, HSPA1A, IER2 and EGR1 was shown in the ST profile from HG2 tissues.
Oxidative stress response in low-grade glioma
We use the Seurat tool to divide the spatial subgroups into 12 groups for both LG1 and LG2. After enrichment analysis by SPATA package in R studio (Supplementary Figures 10A; 11A). We used the Monocle3 package analysis of the SPATA toolkit to establish the trajectory of gene changes in the process from a hypoxic response to a high hypoxic response (from 0.00 to 1.00) (Supplementary Figures 10D; 11D). Unlike high-grade gliomas, the hypoxic response process increases slowly and steadily in low-grade gliomas. They are distributed on the edge of low-grade gliomas in a diffuse state. And the accompanying pathways are mainly normal nerve functions, such as synaptic endocytosis (Supplementary Figure 10E), neural projection (Supplementary Figure 11H), neurotransmitter secretion (Supplementary Figure 11E), and presynaptic activity area. The accompanying pathways are mainly normal nerve functions, such as synaptic endocytosis, neuroprotection, neurotransmitter secretion, and presynaptic activity area. The differentiation of OPC from core areas into the oligodendrocytes triggered the hypoxia areas in LG. Also, increased expression of stress response genes such as JUNB and FOSL2 is observed in the trajectory of LG1 (Supplementary Figure 11I) but not in LG2 (not shown). Neuron functional genes such as CaMK2N1, MEG3, SNCB, SYT1, DKK3 are almost expressed in the marginal area of the tumor (not shown). According to previous MIA analysis, it is believed that it is mainly oligodendrocytes. We believe that the increased oxidative stress is caused by the differentiation of OPC cells in glioma.
The metabolic shift in different types of glioma spatial subtypes
Using spatial transcriptome data, we uncovered an unexpected fact that mitochondrial gene expression differs in high grades of gliomas. In HG1, we found that all 13 human mitochondrial genes in HG1 were primarily expressed in normal oligodendrocytes or AC-like malignant cells rather than OPC-like or MES-like malignant cells (MT-ND1, MT-ND2, MT-ND3, MT-ND4L, MT-ND4, MT-ND5, MT-ND6, MT-CYB, MT-CO1, MT-CO3, MT-ATP6, MT-ATP8, MT-RNR2) except for MT-CO2 (Supplementary Figure 20). In the tumorous sections of HG2, the mitochondrial genes MT-ND1, MT-ND2, MT-ND4, MT-ATP6, MT-CYB, MT-CO3 were completely absent as compared to the normal parts, as shown in Supplementary Figure 21. A similar situation also occurred in HG1; the hypoxic area of HG1 nearly overlaps with the regions that have a poor expression of mitochondrial genes such as MT-CO3, MT-ATP6, MT-CO1, MT-ND3, MT-CYB, and MT-ND1. This finding suggests that mitochondrial dysfunction might be a hypoxia-induced cellular disorder. However, the results will need to be double-checked in further study. We saw a fast elevation in glycolysis-related genes such as PGK1 (Phosphoglycerate Kinase 1), LDHA (Lactate Dehydrogenase A), PKM2 (Pyruvate Kinase M2), and ENO1 (Enolase 1) in both HG1 and HG2 at the same time that mitochondrial genes were absent. When taken together, these findings imply that the tumor metabolic model might have changed under hypoxia stress, possibly shifting from mediated by mitochondria to glycolysis or gluconeogenesis. In The Cancer Genome Atlas (TCGA)-GBM/low grade glioma (LGG) (Supplementary Figure 16) and CGGA datasets (Supplementary Figure 18), we observed a decrease of MT genes in both datasets. Interestingly, mitochondrial gene expression is related to the overall survival of patients with glioma in the CGGA and TCGA-GBM/LGG datasets (Supplementary Figures 17, 19). The percentage of mitochondrial genes was decreased in HG1 and HG2 samples; however, the nUMI and nGene were increased (Supplementary Figure 22). These results indicated a stress response alternation in metabolism for HGs. The metabolic shift allows glioma cells to generate energy more rapidly and supports biosynthetic processes essential for rapid growth. By relying on glycolysis, glioma cells reduce oxidative stress and ROS production, which helps them evade apoptosis and resist cell death. Additionally, this metabolic reprogramming aids in immune evasion and correlates with increased tumor aggressiveness, contributing to a poorer prognosis for patients with high-grade gliomas (47, 48). In summary, this metabolic shift from oxidative phosphorylation to glycolysis is crucial for glioma progression as it enhances tumor proliferation, survival, and adaptation to hypoxic conditions (48, 49).
Discussions
Gliomas are heterogeneous tumors originating from glial cells and remain the deadliest form of brain cancer. Despite significant research, the mechanisms by which glioma stem cell contribute to the diverse cellular composition of gliomas remain inadequately define (50). Although previous studies have revealed the involvement of multiple pathway alternations and genes, such as those involved in the Wnt pathway, the precise mechanisms driving glioma pathogenesis are still not fully understood (51). Advances in single-cell sequencing technologies have provided new insights into tumorigenesis, including glioma, from different aspects. For example, Ochocka et al. identified sex-specific gene expression in glioma-activated microglia that may influence the incidence and outcomes of glioma patients (52). Müller et al. profiled the human gliomas using single-cell RNA-seq and revealed macrophage ontogeny as a factor contributing to regional differences in macrophage activation within the tumor microenvironment (53). Venteiche et al. proposed a unifying model for IDH-mutant gliomas and a comprehensive framework for analyzing the differences across human tumor subclasses (54). The interaction pattern between GSC and immune cells during carcinogenesis demonstrated by Zhai et al. provides a theoretical basis for GSC-targeted immunotherapy (55). However, the ST analysis of primary glioma remained mainly unexplored.
Various subtypes of gliomas have been identified, including OPC-like, AC-like, MES-like and NPC-like. By integrating single-cell and spatial transcriptomics. we observed the distribution of these subtypes in glioma samples. We discovered that oxidative stress response genes play a vital role in the generation of OPC-like subtype. We observed increased expressions of oxidative stress response genes, heat-shock proteins, and DNA damage response genes throughout the evolutionary trajectory of OPC subclusters. These findings suggest that stress responses are activated and potentially enhance the gliomagenesis. Notably, the situation was different between HG and LG gliomas. In HG gliomas, the hypoxic state is associated with uncontrolled proliferation and subsequent malignancy in the tumor core. Especially for HG1, the EMT was observed in the core region. For HG2, the pathways altered were more about the proliferation and cell cycle regulation.
We especially focused on JUN (AP1) and related TFs that regulate apoptosis and cell survival in various types of cancer. Our results revealed that the expression of AP1 in scRNA-seq data was associated with the EMT and expression of abnormal proliferation markers, such as VEGFA and VIM. In ST data, AP1 is enriched in more malignant OPC subtype MES-like and correlated with hypoxia stress. Histological analysis, including H&E and IHC results, confirmed that these genes were expressed in the necrotic areas and glioma pseudo-palisade regions typically surround necrotic foci (53). Pseudo-palisades in glioblastoma are characterized by hypoxia, extracellular matrix protease expression, and an actively migrating cell population (54). These observations were corroborated by IVY, an anatomic transcriptional atlas of human glioblastoma. We used the RNA-Seq as well as the image data from a cellular resolution of ISH tissue sections and adjacent H&E stained sections annotated for anatomic structures and found that JUN (AP1) was highly expressed in necrotic and pseudopalisade regions. These findings implied that oxidative stress-induced expression of JUN-AP1 might impact gliomagenesis and malignant transformation. AP1, a well-known stress-responsive TF, responds to various stimuli, including cytokines, growth factors, and stress signals. It plays a key role in various cancers, including glioma, by regulating its downstream targets. For instance, inhibition of c-Jun N-terminal kinase enhances temozolomide-induced cytotoxicity in human glioma cells (55). FRK controls the migration and invasion of human glioma cells by regulating JNK/c-Jun signaling (56). Similar situations were found for c-FOS, a major component of the AP1 complex. Silencing c-Fos sensitized glioma cells to radiation by increasing radiation-induced DNA double-strand breaks (DSBs), disturbing the DNA damage repair process, promoting G2/M cell cycle arrest, and enhancing apoptosis (57). Our results implied that these TFs might activate in response to oxidative stress sounded by the necrotic areas and trigger the protective role in cancer cells.
Additionally, we observed a loss of mitochondrial genes in HG glioma samples, particularly in the malignant MES-like subtype. Mitochondria are critical for apoptosis activation in mammalian cells; thus, mitochondrial dysfunction is often in gliomas as a mechanism to evade apoptosis (58). Previous studies showed a pathway-based classification of glioblastoma associated with a mitochondrial subtype with therapeutic vulnerabilities (59). These results suggested that mitochondrial subtype in OPC-like might associated with its therapeutic results. For instance, oncogene AIF-2 regulates proliferation, migration, and invasion of human glioma cells via mitochondrial dysfunction (60). HIF-1α, a key regulator of metabolism, affects Tregs migration and immunosuppression by switching between glycolytic and oxidative phosphorylation pathways (61). Moreover, failed apoptosis has a specific transcriptional signature regulated by JNK (an upstream kinase of AP1) and is enriched in metastatic melanoma (62). Our results also showed that mitochondrial gene lost is associated with high hypoxia score in HG gliomas, accompanied by increased expression of EGFR and other proliferative markers. This suggests that stress response TFs might regulate mitochondrial genes.
Furthermore, we observed distinct heterogeneities between LG1 vs LG2 and HG1 vs HG2. LG2 contains more OPC-like malignant components, whereas LG1 has more AC-like malignant components. This suggests that different subtypes of gliomas may have distinct cellular hierarchies and molecular characteristics, potentially influenced by intrinsic variations in glioma stem cells or microenvironmental factors. Additionally, genes in HG2 are less active compared to HG1 regions, likely due to varying microenvironmental conditions and stress responses. HG1 is characterized by a hypoxic state and EMT, driving uncontrolled proliferation, while HG2 focuses more on proliferation and cell cycle regulation pathways. The rapid increase in JUNB 1 expression observed in HG2 regions could be linked to differential stress responses and microenvironmental factors, indicating that glioma cells in these regions might experience unique stress conditions that promote specific gene regulatory mechanisms to enhance survival and proliferation.
Taken together, for the first time based on the scRNA-seq and ST data, this study provides novel insights into oxidative stress-induced distinct pathways in HG and LG, JUN (AP1) might be a critical regulator in the tumorigenesis processes of OPC-like glioma under stress conditions.
Our study, while comprehensive, has several limitations. Firstly, the combination of single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (ST) provided valuable insights, but the spatial resolution might not fully capture the complex cellular interactions and heterogeneity within gliomas. Additionally, our findings are based on a limited number of high-grade (HG) and low-grade (LG) glioma samples, which restricts the generalizability of the results across different glioma subtypes and patient demographics. The cross-sectional nature of our data further limits the understanding of temporal changes during gliomagenesis and the progression from low-grade to high-grade gliomas, necessitating longitudinal studies to capture these dynamic processes. Furthermore, our study primarily relies on transcriptomic data to infer gene function and pathway activation, underscoring the need for functional validation through in vitro and in vivo experiments to confirm the roles of identified genes and pathways in glioma progression. Lastly, while we identified potential therapeutic targets such as JUN (AP1), the translational aspect of these findings requires further investigation, including drug efficacy studies and clinical trials.
To build on our findings, several future directions can be pursued. Enhancing spatial resolution using advanced technologies can provide more detailed insights into the cellular architecture and interactions within the glioma microenvironment. Integrative multi-omics approaches that combine transcriptomics with proteomics, metabolomics, and epigenomics can offer a more holistic view of the molecular mechanisms driving gliomagenesis and progression. Conducting longitudinal studies to monitor the temporal evolution of gliomas and the impact of therapeutic interventions will be crucial for understanding disease dynamics and treatment responses. Functional studies, including gene knockdown or overexpression experiments and the use of animal models, are necessary to validate the roles of key genes and pathways identified in our study. Investigating the clinical relevance of identified biomarkers and therapeutic targets, including their predictive value for patient prognosis and response to treatment, will be a critical step towards translating these findings into clinical practice. Exploring the interactions between glioma cells and their microenvironment, including immune cells, stromal cells, and extracellular matrix components, will provide insights into the mechanisms of tumor progression and resistance to therapy. Lastly, developing targeted therapies that modulate stress response pathways and metabolic shifts identified in our study could offer new treatment strategies for glioma patients.
Data availability statement
The data, presented in the study are deposited in the GEO repository, accession number GSE270355.
Ethics statement
All procedures performed in studies involving human participants were following the ethical standards of Sichuan Provincial People’s Hospital and with the 1964 Helsinki Declaration. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.
Author contributions
ZH: Data curation, Writing – original draft. ZL: Formal analysis, Writing – original draft. QW: Methodology, Writing – original draft. XS: Methodology, Software, Writing – original draft. WZ: Writing – original draft, Methodology, Software. CH: Formal analysis, Writing – original draft. WY: Software, Writing – original draft. HC: Software, Writing – original draft. BG: Funding acquisition, Writing – review & editing. SS: Validation, Visualization, Writing – review & editing. YW: Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Validation, Visualization, Writing – review & editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This study is supported by the Natural Science Foundation of China (81802504), the grant from Health Commission Of Sichuan Provincial, China (No. 23LCYJ035), Sichuan Science and Technology Program, China (2023YFS0107, 2023YFS0106).
Acknowledgments
The authors would like to thank Tingjing Gong from Shanghai OEbiotech Co., Ltd (China) for performing the scRNA-seq analysis.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1452172/full#supplementary-material
Supplementary Table 1 | Baseline information of two patients in this study.
References
1. Wesseling P, Capper D. WHO 2016 classification of gliomas. Neuropathol Appl Neurobiol. (2018) 44:139–50. doi: 10.1111/nan.12432
2. Brain G.B.D. and Other C.N.S.C.C. Global, regional, and national burden of brain and other CNS cancer, 1990-2016: a systematic analysis for the Global Burden of Disease Study 2016. Lancet Neurol. (2019) 18:376–93. doi: 10.1016/S1474-4422(18)30468-X
3. Chen W, Zheng R, Baade PD, Zhang S, Zeng H, Bray FD, et al. Cancer statistics in China, 2015. CA Cancer J Clin. (2016) 66:115–32. doi: 10.3322/caac.21338
4. Ostrom QT, Cioffi G, Gittleman H, Patil N, Waite K, Kruchko C, et al. CBTRUS statistical report: primary brain and other central nervous system tumors diagnosed in the United States in 2012-2016. Neuro Oncol. (2019) 21:v1–v100. doi: 10.1093/neuonc/noz150
5. Nabors LB, Portnow J, Ammirati M, Baehring J, Brem H, Butowski N, et al. NCCN guidelines insights: central nervous system cancers, version 1.2017. J Natl Compr Canc Netw. (2017) 15:1331–45. doi: 10.6004/jnccn.2017.0166
6. Yalamarty SSK, Filipczak N, Li X, Subhan A, Parveen F, Ataide JA, et al. Mechanisms of resistance and current treatment options for glioblastoma multiforme (GBM). Cancers (Basel). (2023) 15(7):2116. doi: 10.3390/cancers15072116
7. Xu S, Tang L, Li X, Fan F, Liu Z. Immunotherapy for glioma: Current management and future application. Cancer Lett. (2020) 476:1–12. doi: 10.1016/j.canlet.2020.02.002
8. Neftel C, Laffy J, Filbin MG, Hara T, Shore ME, Rahme GJ, et al. An integrative model of cellular states, plasticity, and genetics for glioblastoma. Cell. (2019) 178:835–849.e21. doi: 10.1016/j.cell.2019.06.024
9. Reuter S, Gupta SC, Chaturvedi MM, Aggarwal BB. Oxidative stress, inflammation, and cancer: how are they linked? Free Radic Biol Med. (2010) 49:1603–16. doi: 10.1016/j.freeradbiomed.2010.09.006
10. Hayes JD, Dinkova-Kostova AT, Tew KD. Oxidative stress in cancer. Cancer Cell. (2020) 38:167–97. doi: 10.1016/j.ccell.2020.06.001
11. Reczek CR, Birsoy K, Kong H, Reyes IM, Wang Y, Gao P, et al. A CRISPR screen identifies a pathway required for paraquat-induced cell death. Nat Chem Biol. (2017) 13:1274–9. doi: 10.1038/nchembio.2499
12. Sies H, Jones DP. Reactive oxygen species (ROS) as pleiotropic physiological signalling agents. Nat Rev Mol Cell Biol. (2020) 21:363–83. doi: 10.1038/s41580-020-0230-3
13. Dodson M, Castro-Portuguez R, Zhang DD. NRF2 plays a critical role in mitigating lipid peroxidation and ferroptosis. Redox Biol. (2019) 23:101107. doi: 10.1016/j.redox.2019.101107
14. Rojo de la Vega M, Chapman E, Zhang DD. NRF2 and the hallmarks of cancer. Cancer Cell. (2018) 34:21–43. doi: 10.1016/j.ccell.2018.03.022
15. Shen J, Chen M, Lee D, Law CT, Wei L, Tsang FH, et al. Histone chaperone FACT complex mediates oxidative stress response to promote liver cancer progression. Gut. (2020) 69:329–42. doi: 10.1136/gutjnl-2019-318668
16. Kuehne A, Emmert H, Soehle J, Winnefeld M, Fischer F, Wenck H, et al. Acute activation of oxidative pentose phosphate pathway as first-line response to oxidative stress in human skin cells. Mol Cell. (2015) 59:359–71. doi: 10.1016/j.molcel.2015.06.017
17. Norouzi-Barough L, Sarookhani MR, Sharifi M, Moghbelinejad S, Jangjoo S, Salehi R, et al. Molecular mechanisms of drug resistance in ovarian cancer. J Cell Physiol. (2018) 233:4546–62. doi: 10.1002/jcp.26289
18. Sanchez-Perez Y, Soto-Reyes E, García-Cuéllar CM, Cacho-Díaz B, Santamarı́a A, Rangel-López E, et al. Role of epigenetics and oxidative stress in gliomagenesis. CNS Neurol Disord Drug Targets. (2017) 16:1090–8. doi: 10.2174/1871527317666180110124645
19. Chhipa RR, Fan Q, Anderson JL, Muraleedharan R, Huang Y, Ciraolo GM, et al. AMP kinase promotes glioblastoma bioenergetics and tumour growth. Nat Cell Biol. (2018) 20:823–35. doi: 10.1038/s41556-018-0126-z
20. Liang J, Cao R, Wang X, Zhang Y, Wang P, Gao H, et al. Mitochondrial PKM2 regulates oxidative stress-induced apoptosis by stabilizing Bcl2. Cell Res. (2017) 27:329–51. doi: 10.1038/cr.2016.159
21. Park HK, Hong JH, Oh YT, Kim SS, Yin J, Lee A-J, et al. Interplay between TRAP1 and sirtuin-3 modulates mitochondrial respiration and oxidative stress to maintain stemness of glioma stem cells. Cancer Res. (2019) 79:1369–82. doi: 10.1158/0008-5472.CAN-18-2558
22. Wu L, Wang F, Jiang X, Chen Z. PTPN2 induced by inflammatory response and oxidative stress contributed to glioma progression. J Cell Biochem. (2019) 120:19044–51. doi: 10.1002/jcb.29227
23. Zheng L, Wang C, Luo T, Lü B, Ma H, Zhou Z, et al. JNK activation contributes to oxidative stress-induced parthanatos in glioma cells via increase of intracellular ROS production. Mol Neurobiol. (2017) 54:3492–505. doi: 10.1007/s12035-016-9926-y
24. Sharanek A, Burban A, Laaper M, Heckel É, Joyal JS, Soleimani VD, et al. OSMR controls glioma stem cell respiration and confers resistance of glioblastoma to ionizing radiation. Nat Commun. (2020) 11:4116. doi: 10.1038/s41467-020-17885-z
25. McBrayer SK, Mayers JR, DiNatale GJ, Shi DD, Khanal J, Chakraborty AA, et al. Transaminase inhibition by 2-hydroxyglutarate impairs glutamate biosynthesis and redox homeostasis in glioma. Cell. (2018) 175:101–116.e25. doi: 10.1016/j.cell.2018.08.038
26. Lazarev VF, Nikotina AD, Mikhaylova ER, Nudler E, Полоник СГ, Guzhova IV, et al. Hsp70 chaperone rescues C6 rat glioblastoma cells from oxidative stress by sequestration of aggregating GAPDH. Biochem Biophys Res Commun. (2016) 470:766–71. doi: 10.1016/j.bbrc.2015.12.076
27. Combs SE, Schmid TE, Vaupel P, Multhoff G. Stress response leading to resistance in glioblastoma-the need for innovative radiotherapy (iRT) concepts. Cancers (Basel). (2016) 8(1):15. doi: 10.3390/cancers8010015
28. Zhang ZQ, Wang X, Xue B, Zhao Y, Xie F, Wanget SD. Chronic stress promotes glioma cell proliferation via the PI3K/Akt signaling pathway. Oncol Rep. (2021) 46(3):202. doi: 10.3892/or.2021.8153
29. Stahl PL, Salmén F, Vicković S, Lundmark A, Navarro JF, Magnusson JP, et al. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics. Science. (2016) 353:78–82. doi: 10.1126/science.aaf2403
30. Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol. (2018) 36:411–20. doi: 10.1038/nbt.4096
31. He C, Sheng L, Pan D, Jiang S, Ding L, Ma XJ, et al. Single-cell transcriptomic analysis revealed a critical role of SPP1/CD44-mediated crosstalk between macrophages and cancer cells in glioma. Front Cell Dev Biol. (2021) 9:779319. doi: 10.3389/fcell.2021.779319
32. Johnson KC, Anderson KJ, Courtois ET, Gujar AD, Barthel FP, Varn FS, et al. Single-cell multimodal glioma analyses identify epigenetic regulators of cellular plasticity and environmental stress response. Nat Genet. (2021) 53:1456–68. doi: 10.1038/s41588-021-00926-8
33. Fan F, Zhang H, Dai Z, Zhang Y, Xia Z, Cao H, et al. A comprehensive prognostic signature for glioblastoma patients based on transcriptomics and single cell sequencing. Cell Oncol (Dordr). (2021) 44:917–35. doi: 10.1007/s13402-021-00612-1
34. Berglund E, Maaskola J, Schultz N, Friedrich S, Marklund M, Bergenstråhle J, et al. Spatial maps of prostate cancer transcriptomes reveal an unexplored landscape of heterogeneity. Nat Commun. (2018) 9:2419. doi: 10.1038/s41467-018-04724-5
35. Bartoschek M, Oskolkov N, Bocci M, Lövrot J, Larsson C, Sommarin M, et al. Spatially and functionally distinct subclasses of breast cancer-associated fibroblasts revealed by single cell RNA sequencing. Nat Commun. (2018) 9:5150. doi: 10.1038/s41467-018-07582-3
36. Moncada R, Barkley D, Wagner F, Chiodin M, Devlin JC, Baron M, et al. Integrating microarray-based spatial transcriptomics and single-cell RNA-seq reveals tissue architecture in pancreatic ductal adenocarcinomas. Nat Biotechnol. (2020) 38:333–42. doi: 10.1038/s41587-019-0392-8
37. Macosko EZ, Basu A, Satija R, Nemesh J, Shekhar K, Goldman M, et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell. (2015) 161:1202–14. doi: 10.1016/j.cell.2015.05.002
38. Aran D, Looney A, Liu L, Xie T, Fong V, Hsu A, et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat Immunol. (2019) 20:163–72. doi: 10.1038/s41590-018-0276-y
39. Mabbott NA, Baillie JK, Brown H, Freeman TC, Hume DA. An expression atlas of human primary cells: inference of gene function from coexpression networks. BMC Genomics. (2013) 14:632. doi: 10.1186/1471-2164-14-632
40. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. (2013) 29:15–21. doi: 10.1093/bioinformatics/bts635
41. Hafemeister C, Satija R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. (2019) 20:296. doi: 10.1186/s13059-019-1874-1
42. Gao R, Bai S, Henderson YC, Lin Y, Schalck A, Yan Y, et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol. (2021) 39:599–608. doi: 10.1038/s41587-020-00795-2
43. De Falco A, Caruso FP, Su XD, Iavarone A, Ceccarelli M, et al. A variational algorithm to detect the clonal copy number substructure of tumors from scRNA-seq data. Nat Commun. (2023) 14:1074. doi: 10.1038/s41467-023-36790-9
44. Patel AP, Tirosh I, Trombetta JJ, Shalek AK, Gillespie S, Wakimoto H, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. (2014) 344:1396–401. doi: 10.1126/science.1254257
45. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. (2019) 16:1289–96. doi: 10.1038/s41592-019-0619-0
46. Slyper M, Porter C, Ashenberg O, Waldman J, Drokhlyansky E, Wakiro I, et al. A single-cell and single-nucleus RNA-Seq toolbox for fresh and frozen human tumors. Nat Med. (2020) 26:792–802. doi: 10.1038/s41591-020-0844-1
47. Jang M, Kim SS, Lee J. Cancer cell metabolism: implications for therapeutic targets. Exp Mol Med. (2013) 45:e45. doi: 10.1038/emm.2013.85
48. Burns JS, Manda G. Metabolic pathways of the warburg effect in health and disease: perspectives of choice, chain or chance. Int J Mol Sci. (2017) 18(12):2755. doi: 10.3390/ijms18122755
49. Agnihotri S, Zadeh G. Metabolic reprogramming in glioblastoma: the influence of cancer metabolism on epigenetics and unanswered questions. Neuro Oncol. (2016) 18:160–72. doi: 10.1093/neuonc/nov125
50. Laug D, Glasgow SM, Deneen B. A glial blueprint for gliomagenesis. Nat Rev Neurosci. (2018) 19:393–403. doi: 10.1038/s41583-018-0014-3
51. Takebe N, Miele L, Harris PJ, Jeong W, Bando H, Kahn M, et al. Targeting Notch, Hedgehog, and Wnt pathways in cancer stem cells: clinical update. Nat Rev Clin Oncol. (2015) 12:445–64. doi: 10.1038/nrclinonc.2015.61
52. Ochocka N, Segit P, Walentynowicz KA, Wojnicki K, Cyranowski S, Swatler J, et al. Single-cell RNA sequencing reveals functional heterogeneity of glioma-associated brain macrophages. Nat Commun. (2021) 12:1151. doi: 10.1038/s41467-021-21407-w
53. Muller S, Kohanbash G, Liu SJ, Alvarado B, Carrera D, Bhaduri A, et al. Single-cell profiling of human gliomas reveals macrophage ontogeny as a basis for regional differences in macrophage activation in the tumor microenvironment. Genome Biol. (2017) 18:234. doi: 10.1186/s13059-017-1362-4
54. Venteicher AS, Tirosh I, Hebert C, Yizhak K, Neftel C, Filbin MG, et al. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science. (2017) 355(6332):eaai8478. doi: 10.1126/science.aai8478
55. Zhai Y, Li G, Li R, Chang Y, Feng Y, Wang D, et al. Single-cell RNA-sequencing shift in the interaction pattern between glioma stem cells and immune cells during tumorigenesis. Front Immunol. (2020) 11:581209. doi: 10.3389/fimmu.2020.581209
56. Zhou X, Liu H, Zhang W, Zhu M, Shi Q, Li F, et al. FRK controls migration and invasion of human glioma cells by regulating JNK/c-Jun signaling. J Neurooncol. (2012) 110:9–19. doi: 10.1007/s11060-012-0933-1
57. Liu ZG, Jiang G, Tang J, Wang H, Feng GK, Chen FR, et al. c-Fos over-expression promotes radioresistance and predicts poor prognosis in Malignant glioma. Oncotarget. (2016) 7:65946–56. doi: 10.18632/oncotarget.v7i40
58. Guntuku L, Naidu VG, Yerra VG. Mitochondrial dysfunction in gliomas: pharmacotherapeutic potential of natural compounds. Curr Neuropharmacol. (2016) 14:567–83. doi: 10.2174/1570159X14666160121115641
59. Garofano L, Migliozzi S, Oh YT, Angelo FD’, Najac RD, Ko A, et al. Pathway-based classification of glioblastoma uncovers a mitochondrial subtype with therapeutic vulnerabilities. Nat Cancer. (2021) 2:141–56. doi: 10.1038/s43018-020-00159-4
60. Chen W, Líu H, Wang T, Bao G, Wang N, Li RC. Downregulation of AIF-2 inhibits proliferation, migration, and invasion of human glioma cells via mitochondrial dysfunction. J Mol Neurosci. (2019) 68:304–10. doi: 10.1007/s12031-019-01306-y
61. Miska J, Lee-Chang C, Rashidi A, Muroski ME, Chang AL, Lopez-Rosas A, et al. HIF-1alpha is a metabolic switch between glycolytic-driven migration and oxidative phosphorylation-driven immunosuppression of tregs in glioblastoma. Cell Rep. (2019) 27:226–237.e4. doi: 10.1016/j.celrep.2019.03.029
Keywords: glioma, oligodendrocyte precursor cells (OPCs), oxidative stress response, single-cell RNA sequencing, spatial transcriptomics, heterogeneity
Citation: He Z, Liu Z, Wang Q, Sima X, Zhao W, He C, Yang W, Chen H, Gong B, Song S and Wang Y (2024) Single-cell and spatial transcriptome assays reveal heterogeneity in gliomas through stress responses and pathway alterations. Front. Immunol. 15:1452172. doi: 10.3389/fimmu.2024.1452172
Received: 20 June 2024; Accepted: 31 July 2024;
Published: 27 August 2024.
Edited by:
Zuhui Pu, Shenzhen Second People’s Hospital, ChinaReviewed by:
Tuzhen Xu, MyFamily Hospice LLC, United StatesJing Ming Lim, University of Houston, United States
Wenyuan Li, Albert Einstein College of Medicine, United States
Qianwen Zhu, Johns Hopkins University, United States
Copyright © 2024 He, Liu, Wang, Sima, Zhao, He, Yang, Chen, Gong, Song and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yi Wang, w_yi2022@163.com; Siyuan Song, si-yuan.song@bcm.edu; Bo Gong, gongbo2007@hotmail.com
†These authors have contributed equally to this work