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 .

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

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.

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

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.

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

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.

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

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.

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].

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

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

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.

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