Hypoxia Contributes to Poor Prognosis in Primary IDH-wt GBM by Inducing Tumor Cells MES-Like Transformation Trend and Inhibiting Immune Cells Activity

Aims: To reveal the influence of hypoxia on tumor cells and immune cells in primary IDH-wt glioblastoma patients.

Methods: Single-cell RNA-seq data and bulk RNA-seq data were acquired from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases, respectively. Hypoxia status and subtypes of tumor cells were identified based on single-sample Gene Set Enrichment Analysis (ssGSEA). Regulon network analysis of different subtypes under different conditions was conducted by SCENIC. Within tumor microenvironment, biological process activity analysis and cell–cell communication network were conducted to uncover the inner links between each cell subtype under different hypoxia status.

Results: Different types of tumor cell in GBM possessed different hypoxia status, and MES-like subtype was under a more severe hypoxia condition than other subtypes. Hypoxia also induced MES-like signature gene expression within each tumor cell, which could stimulate tumor cell proliferation and invasion by regulating cell–cell communication. Additionally, hypoxia inhibited immune cell activity in the tumor microenvironment by inducing macrophage phenotype polarization and upregulating immune-inhibited cell–cell interaction within immune cells. Interactions between tumor cells and immune cells under hypoxia status also promoted tumor progression.

Conclusions: Hypoxia was a poor prognostic marker for primary IDH-wt GBM patients. Meanwhile, it could induce tumor cells’ MES-like transformation trend and inhibit antitumor function of immune cells.

Introduction

Glioblastoma (GBM) is the most fatal tumor in the central nervous system (CNS) with the median overall survival of 12–15 months. A previous study has divided GBM into four subtypes, namely, mesenchymal (ME), classical (CL), proneural (PN), and neural (NE), based on bulk RNA-seq data (1) and later retain three subtypes including ME, CL, and PN (2). Among these subtypes, ME and CL are associated with poor prognosis compared with others (3). To evaluate the intratumoral heterogeneity, a recent study identified four cellular subtypes in GBM based on single-cell data, which are a bit different from the former Verhaak’s classification (4). It distinguishes GBM cells into various states named neural-progenitor-like (NPC-like), oligodendrocyte-progenitor-like (OPC-like), astrocyte-like (AC-like), and mesenchymal-like (MES-like) based on cell meta modules, which can coexist within the same tumor sample because of intratumoral heterogeneity. Among all tumor cell subtypes, MES-like subtype is linked to hypoxia, indicating the inner association between MES-like subtype and hypoxia.

Hypoxia, a microenvironment feature in solid tumor, is associated with many cancer “hallmarks,” involving impaired immune response, metabolic reprogramming, stimulation of tumor angiogenesis, and promotion of tumor invasion (5). It has been shown to be a factor of chemotherapy and radiation resistance in GBM (6), which is also linked to a proneural–mesenchymal transformation (PMT) in the transcriptome (7). Moreover, hypoxia is a strong inducer of PMT in GBM by upregulating specific transcription factors (TFs) like hypoxia-inducible factor 1-alpha (HIF1α) (8). In the process of PMT, MES-related gene signatures are associated with therapeutic resistance, while PN-related signatures are associated with therapeutic sensitivity (7). Several lines of evidence suggest that PMT occurs in glioma stem cells under hypoxia status (9), then promotes tumor invasion, angiogenesis, and therapeutic resistance.

Collectively, the vast majority of studies focusing on hypoxia effect are based on bulk RNA-seq datasets. Due to intratumoral heterogeneity, tumor purity, and single-cell sequencing development, hypoxia evaluation through bulk data is unable to reveal specific hypoxia status of each cell type, which, however, can be complemented by single-cell analysis. To study the influence of hypoxia on single cells, including tumor and non-tumor cells, especially immune cells in the tumor microenvironment, we performed the study on single-cell RNA-seq datasets derived from primary IDH-wild type (IDH-wt) GBM patients. This work provided an opportunity to investigate the influence of hypoxia on single cells and uncover new possible mechanisms for hypoxia as a poor prognostic marker.

Materials and Methods

Datasets and Data Processing

A total of 7,930 cells from primary IDH-wt GBM patients’ single-cell RNA-seq data according to Smart-seq2 protocol (GSE131928) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo) (4). We additionally downloaded 4,825 cells from oligodendroglioma patients’ and 6,341 cells from IDH-mutant astrocytoma patients’ single-cell RNA-seq data according to Smart-seq2 protocol (GSE70630, GSE89567) to test the efficiency of the biological process activity (BPA) transformation analysis (10–12). Six GBM primary cell lines RNA-seq data (three samples in normoxia and three samples in hypoxia) were downloaded from GSE108013 (13). One hundred forty-five samples of Illumina HiSeq RNA-seq data of primary IDH-wild-type GBM and corresponding phenotype data were downloaded from the TCGA database (https://tcga-data.nci.nih.gov/) via Xena Browser developed by University of California Santa Cruz (UCSC). The QC process, normalization, and cell-type identification of single-cell RNA-seq data were performed as described in published studies (4, 10, 11), which used the same criteria. RNA-Seq data of TCGA and GSE108013 were normalized by transcripts per kilobase million (TPM) method for further analysis. A computational pipeline for analyzing effect of hypoxia in this study is exhibited in Figure 1.

FIGURE 1

Figure 1 Schematic representation of the analysis pipeline in this study.

Subtypes and Glioma Stem-Like Cells (GSCs) Identification

Subtypes of GBM cells and glioma stem-like cells (GSCs) were identified by ssgsea.GBM.classification R package (2). First, we generated numerous virtual samples by randomly selecting expression values of the gene as virtual samples’ corresponding gene expression from datasets. Then, the ssGSEA scores of each category were calculated, respectively. We set 5,000 virtual random samples and correlated these samples with the real sample. At the same time, we counted the number of matches with random samples under each subtype. We defined the subtype as the one that had the fewest matches to the random sample. If there were more than one subtypes sharing the min matches in one sample, we defined the sample as MIX. Besides, we defined the cells with no random sample match as mesenchymal GSCs (mGSCs) or proneural GSCs (pGSCs). The four GBM subtypes (MES-like, AC-like, OPC-like, and APC-like) signatures were acquired from previous research (4). Additionally, mGSCs and pGSCs signatures were obtained from Wang’s study (14). Meanwhile, macrophage-subtype marker genes were obtained from the CellMarker database (15).

Hypoxia Score Calculation and Hypoxia Status Identification

According to the hypoxia 15-gene signature evaluated by published papers (16–20) in pan-cancer analysis, we used GSVA algorithm (21) to calculate the hypoxia score in TCGA bulk samples and single-cell samples, respectively. Unsupervised hierarchical clustering method was employed to classify hypoxia status of bulk samples based on 15 hypoxic gene expression profiles, same as the process in previous research (16). The top three clusters were designed as hypoxia, intermediate, and normoxia groups in TCGA GBM cohorts. Due to the dropout rate and sparse matrix of the single-cell data, we used quarter quantile to identify the hypoxia status of cells. The cells possessing the upper quarter hypoxia score were labeled as hypoxia status, and the cells with the lower quarter hypoxia score were identified as normoxia status; the other cells were defined as intermediate status.

MES-Like Score Calculation and Deconvolution Analysis

We used the same method as calculating hypoxia score to acquire the MES-like score of each GBM tumor cell based on MES-like signatures. The deconvolution of TCGA GBM samples were conducted by CIBERSORTx, in which we got the target cell proportion within the bulk samples (22).

Tumor Purity Calculation in Bulk Samples and Purity–Hypoxia Correlation

We adopted five tumor purity scores in the analysis, namely, the consensus purity estimation (CPE) score, ESTIMATE score, ABSOLUTE score, LUMP score, and IHC score. ESTIMATE was based on specific gene expression profiles consisting of 141 immune genes and 141 stromal genes (23); ABSOLUTE used somatic copy-number data (24); LUMP (leukocytes unmethylation for purity) averaged 44 non-methylated immune-specific CpG sites; and IHC was estimated by image analysis of hematoxylin and eosin stain slides produced by the Nationwide Children’s Hospital Biospecimen Core Resource (25). For CPE, it was the median purity level derived from the four methods mentioned above after normalization as previously described (25). We used Pearson correlation to calculate the coefficient of purity score and hypoxia score among different calculation methods.

Pseudo-Time Lineage Trajectory

Lineage trajectory was constructed by monocle R package (26). Briefly, we selected variant features identified by Seurat v3 as an input gene set. Then, the monocle used a machine learning method to order single cell into a trajectory in a two-dimension space.

Biological Process Activity Analysis and GSEA Enrichment

Biological process activity (BPA) analysis was based on gene list in msigdb database to evaluate each pathway’s activity in single cells (12). We used KEGG database, GO database, and msigdb immune-related C7 database to calculate pathways activity in cells, in which the parameters were set as advised (12, 27). To evaluate the efficiency of the method, we adopted another two single-cell datasets of glioma as formerly mentioned to test it and visualized BPA result by TSNE. Gene set enrichment analysis (GSEA) was performed by clusterProfiler R package based on gene list from msigdb (28).

Identification of Differentially Expressed Genes (DEGs) and Pathways

Differentially expressed genes (DEGs) and differentially expressed pathways between the hypoxia and normoxia groups were identified by limma R package. All parameters were set as default.

Transcription Factor Regulation and Cell Communication Network Identification

We used epigenetic landscape in silico deletion analysis (LISA), based on DEGs, to identify main transcription factors in hypoxia group (29). To analyze transcription factor regulons further, we adopted SCENIC R package (30) with default parameters. For visualization, we mapped the regulon activity (AUC) scores in heatmap plot. Additionally, intratumoral cell–cell communication network based on potential receptor–ligand interaction was inferred by CellPhoneDB from single-cell transcriptomic data (31).

Survival and Statistics Analysis

R packages survival and survminer were used for overall survival analysis. Univariate cox and multivariate cox regression analysis methods were conducted to identify independent prognostic features. In addition, single- and two-factor Kaplan–Meier curves were shown to compare the prognostic results between different groups. All statistical analyses were performed using R software, version 3.6.2 (The R Foundation for Statistical Computing, http://www.rproject.org/). Continuous variables that conform to normal distribution between groups were compared by the Student’s t-test with post-hoc pairwise Bonferroni tests. Variables that do not conform to the normal distribution were compared by the Wilcoxon rank-sum test with post-hoc pairwise Bonferroni tests.

Results

Hypoxia Contributes to Poor Prognosis in Primary IDH-wt GBM Bulk Samples

We divided 145 primary IDH-wt GBM patients from TCGA database into three hypoxia status: hypoxia (n = 78), intermediate state (n = 17), and normoxia (n = 50) based on 15-gene hypoxic signature, which performed the best in evaluating tumor hypoxic condition reported by former study (Supplementary Figure S1A) (16). To investigate the effect of hypoxia, we took hypoxia and normoxia groups into further analysis, and we observed that hypoxia was associated with worse prognosis (Figure 2A). Additionally, we conducted univariate Cox regression on hypoxia status and found that it was a risk factor [hazard ratio (HR) = 1.63; 95%CI = 1.09–2.42, p = 0.016].

FIGURE 2

Figure 2 The effect of hypoxia on tumor cells and prognosis. (A) Kaplan–Meier curve of different hypoxia status of primary IDH wild-type GBM in TCGA database. (B) The expression of 15 hypoxia-related genes in GBM tumor cells. (C) The distribution of hypoxia score within each tumor cell type. (D) Intercellular communications within each tumor subtype. (E) Distribution of MES-like cell proportion within TCGA primary IDH-wt GBM samples between hypoxia and normoxia group. (F) Two-factor Kaplan–Meier curve of TCGA primary IDH-wt GBM. The two factors were MES-like cell proportion and hypoxia status of each sample. The character ***, **** means p < 0.001 and p < 0.0001, respectively.

Identify Hypoxia Status in Primary IDH-wt GBM Single Cells

To thoroughly evaluate hypoxia effects on primary IDH-wt GBM, we calculated hypoxia score on six subtypes of tumor cells in GBM (MES-like, AC-like, OPC-like, NPC-like, mGSC, and pGSC). We first identified four subtypes of GBM cells (MES-like, AC-like, OPC-like, and NPC-like) and then classified all tumor cells into another two GSC subtypes. Interestingly, mGSCs were derived from MES- and AC-like subtypes, while pGSCs were derived mainly from OPC- and NPC-like cells (Supplementary Figure S1B). Afterwards, we divided each of them into three hypoxia status as previous depicted (Figure 2B). These 15 hypoxia marker genes were upregulated in hypoxia groups comparing to intermediate and normoxia groups; meanwhile, we performed LISA analysis based on differentially upregulated genes in the hypoxia group [logFC > 1 and false discovery rate (FDR) < 0.05], discovering that both transcription factors HIF1A and HIF2A (EPAS1) were upregulated in the hypoxia group (Supplementary Table S1). We then compared hypoxia score among cell subtypes to distinguish overall hypoxic condition of each subtype (Figure 2C). Apart from GSCs, MES-like subtype was the most hypoxic one, followed by AC-, OPC-, then NPC-like subtype. To uncover the relationship between hypoxia status and these four tumor cell subtypes, we analyzed the autocrine signaling under different hypoxia status and found that MES-like subtype with hypoxia sent out more epidermal growth factor receptor (EGFR), vascular endothelial growth factor receptor (VEGFR), and fibroblast growth factor receptor (FGFR) signal to itself, indicating that it was more proliferative and invasive than others (Figure 2D). Moreover, we conducted deconvolution analysis on TCGA samples to calculate each subtype proportion by CIBERSORTx and found that hypoxia groups in bulk dataset had more MES-like cells (Figure 2E). Consistent with hypoxia, MES-like cell proportion was identified as a risk factor by univariate Cox regression (HR = 1.48, 95%CI = 1.01–2.18, p = 0.046) and possessed synergistic effect with hypoxic condition on poor prognosis (Figure 2F). In order to explore potential role of the hypoxia and MES-like cell proportion for patients’ prognosis in detail, we included gender, age, O⁶-methylguanine DNA methyltransferase (MGMT) promoter methylation status, and purity to perform multivariate Cox analysis, which suggested that hypoxia was an independent risk factor (Table 1). However, we did not find obvious associations between tumor purity and hypoxia, which indicated that hypoxia was attributed to multiple cell types in the tumor microenvironment (Supplementary Figure S1C).

TABLE 1

Table 1 Statistic of prognostic related factors in TCGA cohort.

Hypoxia Changes Immune-Related Metabolism and Induced MES-Like Signature Expression in GBM Tumor Cells

We performed BPA analysis on primary IDH-wt GBM tumor cells based on KEGG, Msigdb C7, and GO databases. BPA based on single-cell RNA sequencing (scRNA-seq) scRNA-seq data, which could decrease the impact of dropout events and batch effects to provide a robust description on cellular states, had already been tested well in the additional two single-cell datasets (Supplementary Figure S1D) (12). We observed that hypoxia- and HIF1A-related pathways in KEGG database were upregulated in the hypoxia group (Figure 3A and Supplementary Table S2). In addition, there were more immune-related gene sets obviously upregulated in the hypoxia group than in the normoxia group such as interferon (IFN)-, interleukin-15 (IL15)-, and signal transducer and activator of transcription (STAT)-associated signatures, which showed a higher immune activity, whether pro- or anti-tumor, within hypoxia GBM microenvironment (Figure 3B and Supplementary Table S4). However, few GO pathways were obviously changed between hypoxia and normoxia groups (Supplementary Table S2). We then used SCENIC to distinguish regulons between hypoxia and normoxia groups. As a result, ATF3, JUNB, and MYC regulons were highly expressed in the hypoxia group (Supplementary Figure S2A). For more details about hypoxia and cell subtypes, we conducted SCENIC again on both cell subtypes and hypoxia status and found that MES-like subtype shared more common regulons with AC-like subtype revealing that these two subtypes were from the same progenitor that differentiated to different directions due to microenvironment (Supplementary Figure S2B) (32), while NPC-like subtype shared more common regulons with OPC-like subtype (Figure 3C). Within the MES- or AC-like cells, ATF3, a key inhibitory of transcriptional regulator in the inflammatory response, was higher in the hypoxia group. Since autophagy, which could be induced by hypoxia, was a significant catabolic mechanism for GBM tumor cells to survive and resist antitumoral therapy (33), we tested autophagy pathways’ activity between hypoxia and normoxia groups by GSEA. The results showed that these pathways were enriched in the hypoxia group, especially at MES-like subtype, with statistical significance comparing to other subtypes (Figures 3D, E).

FIGURE 3

Figure 3 The different pathway activities between hypoxia and normoxia groups. (A) Heatmap of KEGG metabolism pathway activity between hypoxia and normoxia groups within tumor cells. (B) Heatmap of msigdb C7 immune signatures between hypoxia and normoxia groups within tumor cells. (C) Heatmap of binarized regulon network activity of different tumor cell subtype under hypoxia or normoxia. (D) GSEA analysis plot of autophagy-related pathways in tumor cells under hypoxia condition with p < 0.05. (E) GSEA analysis plot of autophagy-related pathways in MES-like tumor cells under hypoxia condition with p < 0.05.

Meanwhile, we also discovered that CD44, an MES marker gene, was higher in the hypoxia group (Supplementary Figure S2C). This was consistent with the former analysis that MES-like cell proportion was higher in the hypoxia group in bulk tissue (Figure 2E). We then calculated MES-like score among subtypes of tumor cells based on two identified MES-like signatures (4) and found that with the exacerbation of hypoxia, the MES1- and MES2-like score all increased within each subtype (Figure 4A). To verify that hypoxia could induce MES-like signatures expressed in GBM tumor cells, we analyzed another GBM RNA-seq dataset from primary cultured cell lines of primary IDH-wt GBM patients (13). These cells were cultured under hypoxia condition (n = 3, 3% O₂) or normoxia condition (n = 3, 21% O₂) for 48 h, respectively. Consistently, two MES-like score of cell lines under hypoxia were higher than normoxia (Figure 4B). Because of the limited sample size, we did not observe p-values with statistical significance. However, the possibility of type I error between these two groups was 10%. Considering the limitation of sample size, the error rates were considered acceptable.

FIGURE 4

Figure 4 Cell–cell communication network in the tumor microenvironment. (A) The distribution of MES-like score in each tumor cell subtype under different hypoxia status. (B) The distribution of MES-like score in primary cell lines of primary IDH-wt GBM cultured under different hypoxia conditions. Hypoxia group was cultured under 3% O₂ for 48 h, and normoxia group was cultured under 21% O₂ for 48 h. (C) The expression of macrophage marker gene in hypoxia or normoxia group. M1 marker, IRF5 and M2 markers, and CD163 and CD206 were selected. (D) Intercellular interaction within immune cells under different hypoxia conditions. (E) Cell–cell communications between macrophages and tumor cells under different hypoxia conditions. (F) Cell–cell communications between T cells and tumor cells under different hypoxia conditions. The character **, ***, **** means p < 0.01, p < 0.001 and p < 0.0001, respectively. ns, no significance.

Hypoxia Weaken Antitumor Activity of Immune Cells in GBM Microenvironment

The immune microenvironment of primary IDH-wt GBM mainly consisted of tumor-associated macrophages (TAMs) (34). We selected TAMs identified as hypoxia or normoxia for further analysis (Supplementary Figure S2D). The TAMs under hypoxia status highly expressed M2 markers, CD163 and CD206, while TAMs in normoxia upregulated M1 marker, IRF5 (35, 36) (Figure 4C). Meanwhile, TAMs were more active in normoxia compare to those in hypoxia (Supplementary Figures S2E, F), and CCL-related cell–cell interactions were higher in normoxia status between TAMs and tumor cells (Figure 4E). We then tested autocrine signaling within TAMs under different hypoxia status and observed that CD44-related communication, which could induce macrophage M1 to M2 polarization (37), was more frequent and stronger in hypoxia (Figure 4D). Furthermore, we investigated T cells with hypoxia and normoxia status as well but found no obvious change in expression of inhibitory immune checkpoint genes (Supplementary Figure S3A). However, the TIGIT-related communication was a little stronger within hypoxia T cells (Figure 4D), while tumor necrosis factor (TNF)- and FASLG-related communications were stronger in normoxia between T cells and tumor cells (Figure 4F and Supplementary Figure S3B). Otherwise, tumor cells in hypoxia could inhibit TAMs activation via HLA-G/LILRB interactions (Supplementary Figure S3C), assist TAMs polarization via CD44-related communication (Figure 4E and Supplementary Figure S3C), and inhibit T-cell activation via PDCD1 and TIGIT (Figure 4F and Supplementary Figure S3B). Conversely, TAMs could promote tumor cells proliferation and invasion via EGFR-, FGFR2-, and CD44 (38)-related interaction in hypoxia (Figure 4E and Supplementary Figure S3C).

Discussion

The classification of GBM subtypes according to bulk RNA-seq and scRNA-seq possesses inner link as previously reported that TCGA-CL and TCGA-ME subtypes in bulk correspond to tumors enriched with AC- and MES-like cells, respectively, while TCGA-PN subtype corresponds to the combination of OPC- and NPC-like cells (4). TCGA GBM subtypes can also be distinguished by the cell composition that TCGA-CL contain high level of astrocytes, TCGA-ME are infiltrated highly by immune cells, and TCGA-PN are characterized by the highest level of oligodendrocytes and neurons (14). Nowadays, studies on glioma stem cells identify them into two main types, namely, mGSC and pGSC (14, 39). Similarly, the TCGA-ME and TCGA-CL samples are enriched by mGSCs, and TCGA-PN contains high level of pGSCs (14). In our study, we demonstrated that AC-like cells shared more common features with MES-like cells, while OPC-like cells were more similar to NPC-like cells in terms of regulon network and cell–cell communication. Considering all the evidence, we hypothesize that these four GBM cell subtypes come from two lineages derived from mGSC or pGSC, respectively. One is the mesenchymal lineage including mGSC, MES-like subtype, and AC-like subtype, and the other is proneural lineage containing pGSC, OPC-like subtype, and NPC-like subtype. Compared to proneural lineage, the hypoxia degree of mesenchymal lineage was more severe. Such phenomenon is partially associated with metabolism distinction that mGSC preferentially utilizes glycolysis, while pGSC employs oxidative phosphorylation (39). However, the metabolism of each cell subtype of GBM and the lineage hypothesis need further study.

One mechanism of hypoxia inducing GSCs PMT is altering DNA methylation pattern at promoter of mesenchymal-related genes, regulators, or enhancers (7). Nevertheless, we did not find any obvious methylation pattern of mesenchymal markers in methylation data from TCGA (Supplementary Table S3). This result can be influenced by tumor purity in that the methylation data contain not only the tumor cells but also non-tumor cells like immune cells and stromal cells. Unfortunately, single-cell methylation analysis on GBM, which would be a meaningful direction for studying hypoxia effect on single cells, does not exist. Preliminary works about PMT of GBM are focused on GSCs, but in our comprehensive analysis, non-GSCs cells of GBM tumor could also upregulate MES-like signature genes under hypoxia. Although hypoxic MES-like cells were more malignant in inner tumor cell interaction such as invasiveness mediated by CD44-related interaction (38), proliferation induced by EGFR- and MIF-related interactions (40), and angiogenesis activated by VEGF-related interaction, other tumor cell subtypes under hypoxia performed more or less similar manners (Figure 2D). This MES-like transformation trend under hypoxia environment could partially explain the role of hypoxia as the poor prognosis marker of GBM.

Hypoxia is a microenvironment feature that could influence tumor cells and other non-tumor cells simultaneously. In our results, hypoxia was identified as an independent risk factor for prognosis of GBM patients. Although the result of multivariate Cox regression analysis denied MES-like ratio as a risk factor, it could be explained that hypoxia had inner contact with MES-like ratio, thereby influencing MES-like ratio’s impact on prognosis. In addition, we found that hypoxia could inhibit immune activity. Since the immune microenvironment of primary GBM mainly consisted of TAMs (34), we analyzed TAMs in the datasets and found that TAMs exhibited M2 phenotype in hypoxia status, while they exhibited upregulated M1 markers when in normoxia status. The antitumor ability of macrophages was weakened by hypoxia so that it could benefit tumor growth. Not only cells’ autocrine signaling, like gene expression regulation, hypoxia also influences intercellular interaction between immune cells and tumor cells. In hypoxia status, tumor cells inhibited immune cells’ antitumor ability by interacting with inhibitory receptors in immune cells like PDCD1, TIGIT, and LILRB or weakening the immune-stimulation ligand–receptor interactions such as CCL, TNF, and FASLG. Furthermore, TAMs also promoted tumor progression under hypoxia by interacting with corresponding receptors, such as EGFR- and FGFR2-promoting tumor cell proliferation, CD44-promoting tumor invasiveness, and VEGF-inducing angiogenesis to increase invasion ability of tumor cells.

In summary, we explored in this study the effects of hypoxia on single cells and on both tumor cells and immune cells from primary IDH-wt GBM. Hypoxia was identified as a poor prognostic marker for primary IDH-wt GBM patients by analyzing bulk RNA-seq data. We also explored the potential mechanism of hypoxia contributing to poor prognosis on scRNA-seq data and found that hypoxia could induce tumor cells’ MES-like transformation trend and inhibit immune antitumor function such as inducing macrophage M1 to M2 polarization, thereby promoting tumor progression.

Data Availability Statement

All data was downloaded from the public database. The single-cell data was got from the GEO database and the bulk data was from the TCGA database.

Author Contributions

ZX designed the workflow, performed the data analysis, and wrote the manuscript. HL and CH assisted with the data analysis, designed the graphic image, and revised the manuscript. XL assisted with the study design and contributed to the manuscript. All authors contributed to the article and approved the submitted version.

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.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Nos. 81472594 and 81770781). Science Foundation of Hunan Province, China (2019JJ50978). Special funds for the construction of innovative provinces in Hunan Province (2020SK2062).

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021.782043/full#supplementary-material

Supplementary Figure 1 | (A) The expression of 15 hypoxia-related genes in TCGA GBM samples. (B) The composition ratio of 4 GBM tumor cell subtypes in glioma stem like cells. (C) The distribution of tumor purity score under different hypoxia conditions and correlation with hypoxia score. (D) TSNE plot based on BPA transformation values.

Supplementary Figure 2 | (A) Heatmap of binarized regulon network activity in tumor cells with hypoxia or normoxia. (B) Lineage trajectory among mGSCs, MES-like cells and AC-like cells. (C) Distribution of CD44, MES marker, between hypoxia and normoxia groups. (D) The expression of 15 hypoxia-related genes in GBM non-tumor cells. (E) GSEA analysis plot of activated pathways in macrophage under normoxia condition. (F) GSEA analysis plot of activated pathways in macrophage under hypoxia condition.

Supplementary Figure 3 | (A) The relative expression level about inhibitory immune checkpoints of T cells under different hypoxia conditions. (B) Cell-cell communications between tumor cells and T cells under different hypoxia conditions. (C) Cell-cell communications between tumor cells and macrophages under different hypoxia conditions.

Supplementary Table 1 | The expression rank of transcription factors in tumor cells under hypoxia condition identified by LISA.

Supplementary Table 2 | Differentially expressed BPA pathways between different hypoxia groups, hypoxia v.s. normoxia.

Supplementary Table 3 | Differentially methylation condition nearby MES-like signature genes between different hypoxia groups, hypoxia v.s. normoxia.

Supplementary Table 4 | Immune-related signatures unregulated in hypoxia group within tumor cells showed in Figure 3B.

References

1. Verhaak RG, Hoadley KA, Purdom E, Wang V, Qi Y, Wilkerson MD, et al. Integrated Genomic Analysis Identifies Clinically Relevant Subtypes of Glioblastoma Characterized by Abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell (2010) 17(1):98–110. doi: 10.1016/j.ccr.2009.12.020

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Wang Q, Hu B, Hu X, Kim H, Squatrito M, Scarpace L, et al. Tumor Evolution of Glioma-Intrinsic Gene Expression Subtypes Associates With Immunological Changes in the Microenvironment. Cancer Cell (2017) 32(1):42–56.e6. doi: 10.1016/j.ccell.2017.06.003

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Lin N, Yan W, Gao K, Wang Y, Zhang J, You Y. Prevalence and Clinicopathologic Characteristics of the Molecular Subtypes in Malignant Glioma: A Multi-Institutional Analysis of 941 Cases. PloS One (2014) 9(4):e94871. doi: 10.1371/journal.pone.0094871

PubMed Abstract | CrossRef Full Text | Google Scholar

4. 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(4):835–49.e21. doi: 10.1016/j.cell.2019.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Wigerup C, Påhlman S, Bexell D. Therapeutic Targeting of Hypoxia and Hypoxia-Inducible Factors in Cancer. Pharmacol Ther (2016) 164:152–69. doi: 10.1016/j.pharmthera.2016.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Bao S, Wu Q, McLendon RE, Hao Y, Shi Q, Hjelmeland AB, et al. Glioma Stem Cells Promote Radioresistance by Preferential Activation of the DNA Damage Response. Nature (2006) 444(7120):756–60. doi: 10.1038/nature05236

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Segerman A, Niklasson M, Haglund C, Bergström T, Jarvius M, Xie Y, et al. Clonal Variation in Drug and Radiation Response Among Glioma-Initiating Cells Is Linked to Proneural-Mesenchymal Transition. Cell Rep (2016) 17(11):2994–3009. doi: 10.1016/j.celrep.2016.11.056

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Joseph JV, Conroy S, Pavlov K, Sontakke P, Tomar T, Eggens-Meijer E, et al. Hypoxia Enhances Migration and Invasion in Glioblastoma by Promoting a Mesenchymal Shift Mediated by the HIF1α-ZEB1 Axis. Cancer Lett (2015) 359(1):107–16. doi: 10.1016/j.canlet.2015.01.010

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Bhat KPL, Balasubramaniyan V, Vaillant B, Ezhilarasan R, Hummelink K, Hollingsworth F, et al. Mesenchymal Differentiation Mediated by NF-κb Promotes Radiation Resistance in Glioblastoma. Cancer Cell (2013) 24(3):331–46. doi: 10.1016/j.ccr.2013.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Tirosh I, Venteicher AS, Hebert C, Escalante LE, Patel AP, Yizhak K, et al. Single-Cell RNA-Seq Supports a Developmental Hierarchy in Human Oligodendroglioma. Nature (2016) 539(7628):309–13. doi: 10.1038/nature20123

PubMed Abstract | CrossRef Full Text | Google Scholar

11. 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. Sci (New York NY) (2017) 355(6332). doi: 10.1126/science.aai8478

CrossRef Full Text | Google Scholar

12. Ding H, Blair A, Yang Y, Stuart JM. Biological Process Activity Transformation of Single Cell Gene Expression for Cross-Species Alignment. Nat Commun (2019) 10(1):4899. doi: 10.1038/s41467-019-12924-w

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Heiland DH, Gaebelein A, Börries M, Wörner J, Pompe N, Franco P, et al. Microenvironment-Derived Regulation of HIF Signaling Drives Transcriptional Heterogeneity in Glioblastoma Multiforme. Mol Cancer Res MCR (2018) 16(4):655–68. doi: 10.1158/1541-7786.mcr-17-0680

CrossRef Full Text | Google Scholar

14. Wang L, Babikir H, Müller S, Yagnik G, Shamardani K, Catalan F, et al. The Phenotypes of Proliferating Glioblastoma Cells Reside on a Single Axis of Variation. Cancer Discovery (2019) 9(12):1708–19. doi: 10.1158/2159-8290.cd-19-0329

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Zhang X, Lan Y, Xu J, Quan F, Zhao E, Deng C, et al. CellMarker: A Manually Curated Resource of Cell Markers in Human and Mouse. Nucleic Acids Res (2019) 47(D1):D721–d8. doi: 10.1093/nar/gky900

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ye Y, Hu Q, Chen H, Liang K, Yuan Y, Xiang Y, et al. Characterization of Hypoxia-Associated Molecular Features to Aid Hypoxia-Targeted Therapy. Nat Metab (2019) 1(4):431–44. doi: 10.1038/s42255-019-0045-8

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Fox NS, Starmans MH, Haider S, Lambin P, Boutros PC. Ensemble Analyses Improve Signatures of Tumour Hypoxia and Reveal Inter-Platform Differences. BMC Bioinf (2014) 15:170. doi: 10.1186/1471-2105-15-170

CrossRef Full Text | Google Scholar

18. Buffa FM, Harris AL, West CM, Miller CJ. Large Meta-Analysis of Multiple Cancers Reveals a Common, Compact and Highly Prognostic Hypoxia Metagene. Br J Cancer (2010) 102(2):428–35. doi: 10.1038/sj.bjc.6605450

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Thienpont B, Steinbacher J, Zhao H, D'Anna F, Kuchnio A, Ploumakis A, et al. Tumour Hypoxia Causes DNA Hypermethylation by Reducing TET Activity. Nature (2016) 537(7618):63–8. doi: 10.1038/nature19081

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Haider S, McIntyre A, van Stiphout RG, Winchester LM, Wigfield S, Harris AL, et al. Genomic Alterations Underlie a Pan-Cancer Metabolic Shift Associated With Tumour Hypoxia. Genome Biol (2016) 17(1):140. doi: 10.1186/s13059-016-0999-8

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Hänzelmann S, Castelo R, Guinney J. GSVA: Gene Set Variation Analysis for Microarray and RNA-Seq Data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

22. Newman AM, Steen CB, Liu CL, Gentles AJ, Chaudhuri AA, Scherer F, et al. Determining Cell Type Abundance and Expression From Bulk Tissues With Digital Cytometry. Nat Biotechnol (2019) 37(7):773–82. doi: 10.1038/s41587-019-0114-2

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring Tumour Purity and Stromal and Immune Cell Admixture From Expression Data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Carter SL, Cibulskis K, Helman E, McKenna A, Shen H, Zack T, et al. Absolute Quantification of Somatic DNA Alterations in Human Cancer. Nat Biotechnol (2012) 30(5):413–21. doi: 10.1038/nbt.2203

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Aran D, Sirota M, Butte AJ. Systematic Pan-Cancer Analysis of Tumour Purity. Nat Commun (2015) 6:8971. doi: 10.1038/ncomms9971

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Qiu X, Hill A, Packer J, Lin D, Ma YA, Trapnell C. Single-Cell mRNA Quantification and Differential Analysis With Census. Nat Methods (2017) 14(3):309–15. doi: 10.1038/nmeth.4150

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) Hallmark Gene Set Collection. Cell Syst (2015) 1(6):417–25. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. Omics J Integr Biol (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

29. Qin Q, Fan J, Zheng R, Wan C, Mei S, Wu Q, et al. Lisa: Inferring Transcriptional Regulators Through Integrative Modeling of Public Chromatin Accessibility and ChIP-Seq Data. Genome Biol (2020) 21(1):32. doi: 10.1186/s13059-020-1934-6

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Aibar S, Gonzalez-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: Single-Cell Regulatory Network Inference and Clustering. Nat Methods (2017) 14(11):1083–6. doi: 10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Vento-Tormo R, Efremova M, Botting RA, Turco MY, Vento-Tormo M, Meyer KB, et al. Single-Cell Reconstruction of the Early Maternal-Fetal Interface in Humans. Nature (2018) 563(7731):347–53. doi: 10.1038/s41586-018-0698-6

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Xiong Z, Yang Q, Li X. Effect of Intra- and Inter-Tumoral Heterogeneity on Molecular Characteristics of Primary IDH-Wild Type Glioblastoma Revealed by Single-Cell Analysis. CNS Neurosci Ther (2020) 26(9):981–9. doi: 10.1111/cns.13396

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Jawhari S, Ratinaud MH, Verdier M. Glioblastoma, Hypoxia and Autophagy: A Survival-Prone 'Ménage-À-Trois'. Cell Death Dis (2016) 7(10):e2434. doi: 10.1038/cddis.2016.318

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Klemm F, Maas RR, Bowman RL, Kornete M, Soukup K, Nassiri S, et al. Interrogation of the Microenvironmental Landscape in Brain Tumors Reveals Disease-Specific Alterations of Immune Cells. Cell (2020) 181(7):1643–60.e17. doi: 10.1016/j.cell.2020.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Krausgruber T, Blazek K, Smallie T, Alzabin S, Lockstone H, Sahgal N, et al. IRF5 Promotes Inflammatory Macrophage Polarization and TH1-TH17 Responses. Nat Immunol (2011) 12(3):231–8. doi: 10.1038/ni.1990

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Zhang F, Parayath NN, Ene CI, Stephan SB, Koehne AL, Coon ME, et al. Genetic Programming of Macrophages to Perform Anti-Tumor Functions Using Targeted mRNA Nanocarriers. Nat Commun (2019) 10(1):3974. doi: 10.1038/s41467-019-11911-5

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Kwon JH, Kim M, Bae YK, Kim GH, Choi SJ, Oh W, et al. Decorin Secreted by Human Umbilical Cord Blood-Derived Mesenchymal Stem Cells Induces Macrophage Polarization via CD44 to Repair Hyperoxic Lung Injury. Int J Mol Sci (2019) 20(19). doi: 10.3390/ijms20194815

CrossRef Full Text | Google Scholar

38. Mooney KL, Choy W, Sidhu S, Pelargos P, Bui TT, Voth B, et al. The Role of CD44 in Glioblastoma Multiforme. J Clin Neurosci Off J Neurosurgical Soc Australasia (2016) 34:1–5. doi: 10.1016/j.jocn.2016.05.012

CrossRef Full Text | Google Scholar

39. Yan K, Yang K, Rich JN. The Evolving Landscape of Glioblastoma Stem Cells. Curr Opin Neurol (2013) 26(6):701–7. doi: 10.1097/wco.0000000000000032

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Baron N, Deuster O, Noelker C, Stüer C, Strik H, Schaller C, et al. Role of Macrophage Migration Inhibitory Factor in Primary Glioblastoma Multiforme Cells. J Neurosci Res (2011) 89(5):711–7. doi: 10.1002/jnr.22595

PubMed Abstract | CrossRef Full Text | Google Scholar