Transcriptomic, copy number variations (CNV), and clinical data were extracted from The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) databases, and 952 samples with lower grade glioma were finally included. 508 samples extracted from the TCGA database were defined as the training set, while 444 from the CGGA database were set as the validation set. Normal or glioblastoma (GBM, WHO grade IV glioma) samples were excluded.

We evaluated the expression distribution of HMOX1 between tumor and normal tissues in TCGA pan-cancers. HMOX1 expression between different clinic-related subgroups were explored in both cohorts. The associations between HMOX1 expression and patient outcomes, including overall survival (OS), disease specific survival (DSS), and progression free interval (PFI), in TCGA pan-cancer sets were evaluated using univariate Cox analysis and displayed using R package “forestplot.” And the survival curves were correspondingly established by Kaplan-Meier analysis to evaluate the relationship of LGG patients’ prognosis and HMOX1 expression level as well as mutation status.

The Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE; R package “estimate”) analysis was employed to measure the tumor purity. Meanwhile, we used the CIBERSORT and the Tumor Immune Estimation Resource (TIMER) algorithm to assess the fractions of human immune infiltrating cell types (Li et al., 2017; Chen et al., 2018). The abundance of these cells was compared between groups with either high or low HMOX1 expression levels (according to the median value) in WHO grade II/III gliomas. Besides, the correlations between HMOX1 and tumor stemness, tumor mutational burden, as well as microsatellite instability (MSI) were assessed based on Spearman’s correlation analysis, which were displayed as radar charts using R package “radar.”

We included 29 glioma samples from the Department of Neurosurgery, Xiangya Hospital from February 2019 to February 2021. Patients with recurrent gliomas, glioblastomas, or other cancers, or serious underlying diseases were excluded. Five fresh glioma samples were collected and then immediately stored in 4% paraformaldehyde in room temperature. Slides were sequentially incubated in graded ethanol after deparaffinization for 3 h at 60 degrees. Antigen was exposed using citrate buffer (pH = 6.0). After blocking, slides were treated with anti-HMOX1 (rabbit, AiFang biological AF300167, 1:100) in antibody diluent (abcam ab64211), and subsequent HRP Goat anti-rabbit IgG (H+L) secondary antibody (ABclonal AS014, 1:1,000).

We acquired scanned images on immunohistochemical sections using a digital scanner, and the area of tissue measurement was automatically read using an image analysis system (Servicebio). The positive grade was divided as follows: 0 for negative staining (without staining); 1 for weakly positive (light yellow staining); 2 for moderately positive (brownish yellow staining); and 3 for strongly positive (tan staining). The histochemistry score (H-score) was calculated to reflect the degree of positivity using the following formula:

In the formula, “p_i” indicates percentage of negative/weak/moderate/strong intensity area. “i” represents the positive grade. H-score ranges from 0 to 300, where a higher H-score means a stronger positivity.

Twenty-four samples were collected and then stored in liquid nitrogen for further sequencing on a BGISEQ-500 platform (BGI-Shenzhen, China). The gene expression levels were calculated using RSEM (v1.2.12). This study was approved by the Ethics Committee of Xiangya Hospital (No. 2017121019). The written informed consents were obtained in advance from all participants or their family representatives.

We calculated the correlations between HMOX1 and other genes to identify HRGs (| correlation coefficient| > 0.6; false discovery rate (FDR) < 0.001), which were recorded for subsequent functional enrichment analysis. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms were used to explore the biological functions of the gene set. Next, we developed a HRGs based prognostic signature for the LGG patients by performing the least absolute shrinkage and selection operator (LASSO) Cox regression analysis based on the R package “glmnet” (Friedman et al., 2020), The risk score calculating formula is:

Where “n” means the number of genes, “β_i” is the coefficient for each gene, “x_i’’ means the expression value (log transformed) of each gene. The predicted protein-protein interactions (PPI) among these model genes were achieved by GeneMANIA.¹ And their correlations were also calculated.

We then calculated the risk score for the patients, and divided them into low- and high-risk groups according to the median risk score. The relationships between risk signature and survival as well as other clinicopathological characteristics were also assessed.

The efficacies of four therapies (radiotherapy, chemotherapy, targeted therapy, and immunotherapy) in high-risk and low-risk groups were evaluated.

The Tumor immune dysfunction and exclusion (TIDE) was used to assess the response of LGG patients to the immunotherapies (Jiang et al., 2018). Three resources for therapeutic biomarker discovery in cancer cells, including Genomics of Drug Sensitivity in Cancer (GDSC), Cancer Therapeutics Response Portal (CTRP), and Cancer Cell Line Encyclopedia (CCLE), were employed to evaluate the relationship between drug sensitivity [half maximal inhibitory concentration (IC50)] and mRNA expression of genes.

Single-cell RNA sequencing data derived from the single-cell dataset GSE84465 which contains four cell groups and 3,590 cells (2,342 tumor cells and 1,248 periphery cells) (Darmanis et al., 2017). The correlations between gene expression and functional states including angiogenesis, apoptosis, cell cycle, differentiation, DNA damage, DNA repair, epithelial-mesenchymal transition (EMT), hypoxia, inflammation, invasion, metastasis, proliferation, quiescence, and stemness, were calculated based on the CancerSEA database² (Yuan et al., 2019).

We constructed a nomogram integrating clinical information and prognostic signature by R package “rms” using the variables screened out by Cox regression analysis, which could predict the 1-, 3-, and 5-year overall survival of LGG patients in an intuitionistic and easy-to-utilize manner. Calibration plots and C-indexes were obtained to access the predictive accuracy of the model.

We used Kaplan-Meier curve combined log-rank test to evaluate the patient survival. Subgroups were stratified on the basis of clinical, pathological, and molecular features including age (≤ 40 or > 40-year-old), gender (male or female), grade (II or III), histological type (astrocytoma, oligoastrocytoma, or oligodendroglioma), and IDH1 status (mutant or wildtype). When comparing variables between groups, we used Wilcoxon test or Kruskal-Wallis test. Most statistical analyses were achieved via R language (version 4.0.3). And P < 0.05 was considered statistically significant.

The flow chart of our study was displayed as Supplementary Figure 1. And the characteristics of included patients were summarized in Table 1. HMOX1 expression levels were significantly elevated in tumor samples compared to normal controls (Figures 1A,B). Importantly, this trend remained consistent at single-cell level (Figure 1C). Similarly, glioma tissue had elevated protein level of HMOX1 than normal cerebral cortex tissue (Figures 1D,E). In the TCGA cohort, youngers (≤40 years old) and males had relatively higher HMOX1 expression levels (Figure 1F), while no significant difference was observed in CGGA cohort (Figure 1G), which could be partially attributed to the heterogeneity between cohorts.

Furthermore, we comprehensively investigated the expression distribution of HMOX1 at tissue and cell levels in HPA database (Supplementary Figure 2). In normal tissues, HMOX1 was enriched in lymphoid tissues, especially spleen (Supplementary Figure 2A). Regarding cell types, HMOX1 demonstrated a high expression level in blood and immune cells including monocytes, macrophages, Kupffer cells, and Hofbauer cells (Supplementary Figure 2B). And in both tumor tissues and cell lines, HMOX1 was found to be elevated in glioma tissues and cells (U-138 MG and U-87 MG; Supplementary Figures 2C,D).

Moreover, the mutational profiles of HMOX1 have been illustrated in Supplementary Figure 3. Amplification stood as the most common mutation type in LGG cohort (Supplementary Figure 3A). R123H/C site alteration was frequently observed among HMOX1 mutations (Supplementary Figure 3B). Compared to the HMOX1-wildtpe group, HMOX1-mutant group shared significantly higher mutation frequencies of TP53, TTN, MCM5 in pan-cancers (Supplementary Figure 3C), and higher mutation frequencies of MCM5, APOL6, MB in lower grade gliomas (Supplementary Figure 3D).

Elevated HMOX1 expression was significantly correlated with an unfavorable OS, DSS, and PFI in both uveal melanoma (UVM) and LGG (Figures 2A–C). Subsequent Kaplan-Meier survival analysis showed that HMOX1^high patients had significantly poorer OS (p = 0.0020), poorer DSS (p = 0.0064), and poorer PFI (p = 0.0240) in TCGA cohort (Figures 2D–F). Furthermore, LGG patients with HMOX1 wildtype had better survival than those with HMOX1 mutation (Figures 2G–I). Specifically, patients with HMOX1 CNV deletion demonstrated the worst prognosis which can be explained by the fact that lower CNV (deletion when < 0) represented higher HMOX1 gene expression, as CNV was negatively correlated with expression for HMOX1 (Supplementary Figures 4A,B). Notably, when exploring the impact of the CNV status on the prognosis of glioma patients, the sample sizes varied considerably. Briefly, there were 416 individuals in the HMOX1 “wild type” group vs. only 26 in the “amplification” group. This between-group variation is highly likely to contribute to a misinterpretation of statistical results.

Positively correlations were observed between HMOX1 expression and regulatory T cells (Tregs), activated CD4+ memory T cells, as well as M2 macrophages (Figure 3A) in most cancer types. Regarding LGG, significant correlations were observed between higher HMOX1 expression with increased abundance of naïve B cells, resting CD4+ memory T cells, Tregs, and M2 macrophages, but with decreased abundance of memory B cells and naive CD4+ T cells (Figure 3B). Furthermore, a negative relationship existed between HMOX1 expression and tumor purity in lower grade gliomas, as it positively correlated with stromal, immune, and ESTIMATE scores (Figures 3C–E). Moreover, we found that the abundances of CD4+ T cells, neutrophils, macrophages, and dendritic cells were positively associated with HMOX1 expression based on TIMER algorithm (Figure 3F).

To further explore the specific correlations between HMOX1 and subpopulations of M2 macrophages. We summarized the functions, inducers, and surface markers of M2a∼d macrophages (Colin et al., 2014; Martinez and Gordon, 2014; Wang et al., 2019). Then we calculated the correlation between HMOX1 and each marker to reflect the correlations between HMOX1 mRNA expression and subclassified M2 macrophages (Table 2). The results showed that HMOX1 expression was strongly correlated with M2b macrophages, as it was significantly correlated with M2b markers including CD86 (R = 0.752, p < 0.0001; Supplementary Figure 4C), HLA-DR (R = 0.705, p < 0.0001; Supplementary Figure 4D), HLA-DP (R = 0.672, p < 0.0001; Supplementary Figure 4E), and HLA-DQ (R = 0.621, p < 0.0001; Supplementary Figure 4F).

Using our own samples, we validated that grade III gliomas had higher HMOX1 protein expression compared with grade II gliomas (Figure 4A). The detailed immunohistochemical staining results were summarized in Supplementary Table 2. Furthermore, we confirmed the positive relationships between HMOX1 expression levels and immune infiltrates including B cell and resting memory CD4+ T cell (Figures 4B,C), which was in consistent with the results obtained from the TCGA cohort. However, the poor correlations between HMOX1 expression and cancer stemness, TMB, and MSI (Figures 3G–I) indicated that HMOX1 was unlikely to influence oncogenic processes by engaging in genetic alterations or epigenetic modifications. Furthermore, HMOX1 was significantly correlated with several recognized immune checkpoints including leukocyte associated immunoglobulin like receptor 1 (LAIR1), CD80, programmed cell death 1 ligand 2 (PDCD1LG2), CD40, and CD86 (Supplementary Figure 4G).

We identified 505 HRGs (| Spearman R| > 0.6, FDR < 0.001). Mainly located on cytoplasmic vesicle membrane, HRGs participated molecular functions such as immune receptor activity and interactions. In terms of biological processes, the gene set was mainly involved in leukocyte activity involved in immune response, leukocyte activation, positive regulation of immune response, etc. (Figure 5A). Regarding KEGG pathways, HRGs were intimately involved in tumor immunity such as NOD-like receptor signaling pathway, cytokine-cytokine receptor interaction, and NK cell mediated cytotoxicity (Figure 5A).

A total of 377 HRGs were significantly prognostic in the training set following a univariate Cox analysis (p < 0.05). A LASSO Cox analysis was performed on training dataset and a risk signature was then constructed containing 27 genes (Figure 5B). The result of PCA analysis indicated that the signature could discriminate lower grade gliomas from GBM and normal brain tissues (Figure 5C). Moreover, the coefficients and univariate Cox analysis results were displayed in Figures 5D,E and Supplementary Table 3. Protein-protein interactions and transcriptomic correlations were also illustrated in Figures 5F,G.

Risk plots and survival diagrams were displayed as Figure 6A. Kaplan-Meier curve indicated that high-risk LGG patients had significantly poorer prognosis (p < 0.0001, Figure 6B). The time-dependent ROC curves exhibited a promising ability of the model to predict OS in the training cohort (1-year AUC = 0.83, 3-year AUC = 0.88, 5-year AUC = 0.91; Figure 6C). The results were similar in the external CGGA cohort (Figure 6D). Higher risk scores also indicated poorer OS (p < 0.0001, Figure 6E). The risk model still exhibited stable and high predication ability (1-year AUC = 0.75, 3-year AUC = 0.77, 5-year AUC = 0.78; Figure 6F). Furthermore, we performed a stratification analysis and found that the risk model retained the ability to predict OS in various subgroups in both cohorts (Figures 6G,H).

Sankey diagrams were displayed showing the distribution of risk score and clinicopathologic characteristics among LGG patients (Figures 7A,B). LGG patients with increasing age and higher WHO grade possessed higher risk scores. Besides, an individual patient would have a higher risk score if his pathologic type of glioma was an astrocytoma or if he carried a wildtype IDH1 (Figures 7C,D). As for the immune microenvironment of glioma, samples with higher risk demonstrated higher percentages of resting CD4+ memory T cells, M1, M2 macrophages, and activated dendritic cells in both cohorts (Figures 7E,F and Supplementary Table 4).

As for the association between risk signature and therapeutic benefit, most therapies showed poor efficacy in both low-risk (Figures 8A–C) and high-risk groups (Figures 8E–G) in the TCGA cohort. However, patients who had lower risk scores exhibited uniquely better prognosis if they had undergone immunotherapies (p = 0.0078; Figure 8D). But the condition was not observed in the high-risk group (p = 0.0620; Figure 8H).

We then decoded the different functional states of malignant cells and target molecules at single-cell resolution and discovered three functional states which were significantly related to HMOX1, including angiogenesis (correlation = 0.44, p < 0.001; Figure 8I), hypoxia (correlation = 0.40, p < 0.001; Figure 8J), and metastasis (correlation = 0.31, p < 0.01; Figure 8K). Furthermore, patients with lower HMOX1 expression levels had high TIDE scores and exhibited better responses to the immunotherapies (Figure 8L). This is in consistent with the results of prediction of therapy efficacy in different risk subgroups. In addition, three drug databases were accessed to find relationship between drug sensitivity and HMOX1 expression. TKI258 (dovitinib) showed significant sensitivity in the CCLE database (Supplementary Figure 5A). WZ3105 and Repligen 136 were selected from the CTRP (Supplementary Figure 5B) and GDSC (Supplementary Figure 5C) drug banks, respectively.

Multivariate Cox analyses showed that age at first pathological diagnosis, WHO grade, histological type, and risk score were independent predictors of OS in both cohorts (Figure 9A). We then constructed a nomogram model incorporating these prognostic variables as a clinically applicable quantitative tool to predict survival in patients with primary lower grade gliomas (Figure 9B). The calibration plots demonstrated satisfactory predictive performance of the nomogram in both training (Figures 9C–E) and validation (Figures 9F–H) cohorts.