Heme Oxygenase-1 Predicts Risk Stratification and Immunotherapy Efficacy in Lower Grade Gliomas

Background: Gliomas are the most common tumors in human brains with unpleasing outcomes. Heme oxygenase-1 (HMOX1, HO-1) was a potential target for human cancers. However, their relationship remains incompletely discussed. Methods: We employed a total of 952 lower grade glioma (LGG) patients from TCGA and CGGA databases, and 29 samples in our hospital for subsequent analyses. Expression, mutational, survival, and immune profiles of HMOX1 were comprehensively evaluated. We constructed a risk signature using the LASSO Cox regression model, and further generated a nomogram model to predict survival of LGG patients. Single-cell transcriptomic sequencing data were also employed to investigated the role of HMOX1 in cancer cells. Results: We found that HMOX1 was overexpressed and was related to poorer survival in gliomas. HMOX1-related genes (HRGs) were involved in immune-related pathways. Patients in the high-risk group exhibited significantly poorer overall survival. The risk score was positively correlated with the abundance of resting memory CD4+ T cells, M1, M2 macrophages, and activated dendritic cells. Additionally, immunotherapy showed potent efficacy in low-risk group. And patients with lower HMOX1 expression were predicted to have better response to immunotherapies, suggesting that immunotherapies combined with HMOX1 inhibition may execute good responses. Moreover, significant correlations were found between HMOX1 expression and single-cell functional states including angiogenesis, hypoxia, and metastasis. Finally, we constructed a nomogram which could predict 1-, 3-, and 5-year survival in LGG patients. Conclusion: HMOX1 is involved in immune infiltration and predicts poor survival in patients with lower grade glioma. Importantly, HMOX1 were related to oncological functional states including angiogenesis, hypoxia, and metastasis. A nomogram integrated with the risk signature was obtained to robustly predict glioma patient outcomes, with the potential to guide clinical decision-making.


INTRODUCTION
Gliomas, originating from intrinsic constituent cells of the brain, are the most common primary tumors there (Sanai et al., 2005). Before genome-wide molecular profiling researches revealed the comprehensive genomic landscape of human gliomas (Brat et al., 2015;Suzuki et al., 2015), glioma classification was largely based on their microscopic and immunohistochemical features. Development at the molecular biology level identifies novel biomarkers for optimized classification strategy and promising treatment targets. Predictive biomarkers recognized and used in clinics mainly included isocitrate dehydrogenase (IDH) mutation, the discovery of which constituted a key breakthrough in the understanding of WHO grade II/III gliomas (Yan et al., 2009). Besides, the presence of O6-methylguanine-DNA methyltransferase (MGMT) promoter methylation predicts benefit from temozolomide-based chemotherapy in patients with IDH-wildtype glioma (Wick et al., 2012). Furthermore, 1p/19q codeletion is predictive for benefit from combined radiotherapy and chemotherapy (with procarbazine, lomustine, and vincristine) in two phase III trials (Cairncross et al., 2013;van den Bent et al., 2013). Novel pathogenesis-based treatments targeting oncogenic signaling pathways such as BRAF mutation (Robinson et al., 2014), epidermal growth factor receptor (EGFR) amplification (Phillips et al., 2016), and fibroblast growth factor receptor (FGFR)-TACC fusion (Stefano et al., 2015) harbored the potential for LGG treatment.
Although these molecular subclassifications deepened our understanding of tumorigenesis and personalized therapeutics, a certain LGG population still acquired resistance to these targeted therapies. Furthermore, gliomas are not considered highly immunogenic due to the low mutational loads, besides they are featured by severe immunosuppression mediated by immune-inhibitory factors, such as programmed celldeath 1 ligand 1 (PDCD1LG1) and secreted transforming growth factor β (Nduom et al., 2015). Therefore, we hope that newly identified biomarkers can overcome immunosuppression, exploit antitumor immune responses, and guide individualized treatments.
Heme oxygenase is an essential enzyme in heme catabolism as it cleaves cellular heme to form biliverdin. HMOX1 overexpression is observed in various solid malignancies, including bladder (Miyata et al., 2014), breast (Noh et al., 2013), colon (Yin et al., 2014), glioma (Gandini et al., 2014), lung (Degese et al., 2012), prostate (Maines and Abrahamsson, 1996), and gastric (Yin et al., 2012), cancers. Although HMOX1 prevents DNA damage under normal conditions, HMOX1 overexpression paradoxically promotes cancer cell proliferation and invasiveness at late phase of tumorigenesis Ph et al., 2013). Targeting HMOX1 was effective for hormone-refractory prostate cancer (Alaoui-Jamali et al., 2009) and it has been shown to reverse imatinib resistance in myeloid leukemia (Mayerhofer et al., 2008). Several imidazole-based non-porphyrin HMOX1 inhibitors were recently developed, which exhibited selectivity toward HMOX1 . Meanwhile, they showed potent anti-tumor activities both in vitro and in vivo, with the potential for clinical applications . In a nutshell, these findings shed a light on the future for targeting HMOX1 to promote cancer immunotherapy.
Currently, several clinical trials are attempting to explore the clinical benefits of targeting HMOX1 (or related molecules) in the treatment of lower grade gliomas and other solid cancers (Supplementary Table 1). We comprehensively analyzed two independent glioma cohorts, as well as samples from our institution, to explore the HMOX1 profiling in the context of gliomas. In addition, an HRG-based risk signature was established to predict the outcome of patients diagnosed with primary lower grade gliomas. The multifaceted performance of the IIRS was also examined to reveal its superior predictive ability of response to immunotherapy.

Data Extraction
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.

Expression, Mutational, and Survival Analysis
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.

Immune Infiltration Analysis in Lower Grade Gliomas
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."

Immunohistochemical Staining and RNA Sequencing of Glioma Samples
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.

Construction and Validation of Risk Signature Based on HMOX1-Related Genes
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. 1 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.

Prediction of Therapy Efficacy, Correlation With Pathways in Single Cell Landscape, and Drug Response
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.

Establishment of Nomogram
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.

Statistical Analyses
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.

Heme Oxygenase-1 Is Overexpressed in Lower Grade Glioma
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

Heme Oxygenase-1 Predicts Poorer Survival Outcome in Lower Grade Glioma
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.

Heme Oxygenase-1 Participates in Immune Infiltration and Tumor Microenvironment
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; 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).

Risk Signature Construction in the Training Set
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)
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).

Risk Signature Correlates With Clinicopathological Characteristics, Immune Microenvironment, and Therapeutic Efficacy
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; The abundance of immune infiltrates between HMOX1 high and HMOX1 low groups in 24 glioma samples. (C) The correlations between HMOX1 expression and immune infiltrates in 24 glioma samples. *p < 0.05, **p < 0.01, ***p < 0.001, ns, not significant. 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   The risk score in different subgroups stratified by age, gender, WHO grade, histological type, and IDH1 status in the TCGA and CGGA cohorts. (E,F) Correlations between risk score and immune infiltration in the TCGA and CGGA cohorts. ***p < 0.001, ****p < 0.0001, ns, not significant. A, astrocytoma; O, oligodendroglioma; OA, oligoastrocytoma; IDH1, isocitrate dehydrogenase 1.

Construction and Validation of the Nomogram
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.

DISCUSSION
Molecular profiling of lower grade gliomas, common tumors in the human brain, has enhanced insights of molecular oncology and identified few prognostic and therapeutic targets. Over the past decades, surgical section following radiotherapy and chemotherapy still serve as the mainstream means for glioma treatment. Due to the unique microenvironment and low mutational burden, gliomas acquire immunosuppressive phenotypes and poorly response to established immunotherapies. Clinicians thus need renewed molecular targets for better response rate as well as improved prognosis. Comprehensive bioinformatic analyses were used at multiple levels of expression, survival, and biological function to demonstrate the essential role of HMOX1. HMOX1 and the risk signature are reliable predictors of LGG prognosis, and targeting HMOX1 will help clinicians optimize the management of lower grade glioma patients in future preclinical and clinical practice.
The mechanism can be explained by the antioxidant role of HMOX1 in malignancies. Heme oxygenase-1 is a phase II enzyme that responds to adverse conditions, such as oxidative stress, cellular injury, and diseases. Its diverse roles in tumorigenesis have been well reviewed (Nitti et al., 2017;Chiang et al., 2018). Furthermore, HMOX1 has been widely recognized to play a cytoprotective role in tumor cells to overcome the assault of enhanced oxidative stress in the tumor microenvironment, thereby preventing the cancer cells from apoptosis and autophagy (Chiang et al., 2018). On the other hand, our study found that HMOX1 was enriched in blood and immune cells, especially monocytes and macrophages, which were thought to be risk factors for the overall survival of glioma patients (Zhang et al., 2021a,b). Therefore, a higher HMOX1 expression level represented a stronger resistance to oxidative stress and a higher abundance of risk immune cell types, resulting an unfavorable outcome for LGG patients.
We observed that HMOX1 expression was increased in lower grade gliomas at both tissue and single-cell levels, and such increase was associated with worse OS, DSS, and PFI in patients with LGG. Thereafter we suspect that HMOX1 might influence the malignant properties and patient outcomes by intervening immune infiltration in the TME, as it was significantly associated with tumor purity, immune checkpoints' levels, immune infiltration abundance, as well as immuneoncological pathways. These findings were partially validated by our own samples.
In this study, 27of 377 prognostically valuable HRGs were screened out to establish a risk signature. Among them, ARHGDIB (Su et al., 2019), LRRC25 , PLAUR (Tan et al., 2020), and TREM1 (Kong et al., 2020) were selected as prognostic hub genes in lower grade glioma in previous bioinformatic analyses. Besides, GMFG (Lan et al., 2021), GNG5 , S100A10 (Zhang et al., 2021c), TNFRSF12A (Tran et al., 2006), and TYROBP (Lu et al., 2021) were found to be upregulated and correlated with worse prognosis in LGG. Importantly, CAPG was identified as a prognosis factor correlated with macrophages in gliomas (Wei et al., 2020;Prescher et al., 2021). Moreover, a recent study reported that targeting CLIC1 reduced gliomagenesis in tumoral stem or progenitor cells, indicating CLIC1 as a potential target and prognostic biomarker there Setti et al., 2013). Briefly, previous literature indicated that the model genes and encoded proteins were related to the progression of gliomas, which were perfectly consistent to our risk signature. The risk signature we constructed may reflect the distribution of immune infiltrates involved in the TME, as well as reliably and independently predicts prognosis of patients with primary gliomas. In addition, combining the constructed signature with clinicopathological features, we obtained a nomogram with strong predictive power which is of potential in clinical utilization.
Intriguingly, both the risk signature and HMOX1 expression are correlated with fractions of infiltrating cells, especially resting memory CD4+ T cells and M2 macrophages. Resting memory CD4+ T cells were indicated to play an vital role in latent HIV-1 infection (Siliciano and Siliciano, 2015), its functions is yet warranted to be discovered in glioma. Additionally, abundance of M2 macrophages was associated with immune suppressive phenotype and short-term relapse after radiation therapy . Moreover, M2b macrophages are the most closely related to HMOX1 in M2 macrophage subsets. Considering its involvement in type 2 helper T (Th2) activation, immunomodulation, and tumor progression, combined with the fact that Th2/Th1 ratio is a poor prognostic marker in gliomas (Kumar et al., 2006), we have good reasons to believe that HMOX1 is involved in immunoregulation in glioma microenvironment, driving tumors toward more malignant properties as well as a worse prognosis. Further work is certainly imperative for the validation of HMOX1 functions in lower grade glioma. Besides, selectively targeting HMOX1 was effective for therapy resistance in various cancers (Mayerhofer et al., 2008;Alaoui-Jamali et al., 2009). Several recently developed imidazole based nonporphyrin HMOX1 inhibitors are highly selective for HMOX1 without affecting other heme-containing proteins , and show potent antitumor activities in preclinical models , exhibiting high-quality prospects for clinical use.
A key finding in our study is that immunotherapy demonstrated unique superiority in the low-risk group as patients with vs. without immunotherapy exhibited significantly better survival. Considering that patients with lower HMOX1 expression had lower TIDE score and showed higher response rate to immunotherapy, we conclude that HMOX1, as well as HMOX1-based risk signature, is capable for predicting the immunotherapy efficacy in lower grade gliomas. Regarding the mechanisms, we found that both HMOX1 and risk signature were correlated with M2 macrophages. An important and promising observation for cancer immunotherapy is that macrophages can direct T or B cell responses in the presence or absence of specific antigens (Mills et al., 2016). In particular, M2 macrophages induce T cells into Treg and other T cell type responses without anticancer activity via innate signals such as transforming growth factor-β (TGFβ) and interleukin 10 (IL-10) (Mills et al., 2000;Ruffell and Coussens, 2015). Therefore, we believe that the key to the sensitivity of immunotherapies to the risk score is that the HMOX1-based risk score well reflects the distribution of M2 macrophages in gliomas.
There were several limitations in our study. First, our sample size is too small to be completely convinced. Second, only one cohort with single-cell sequencing data was employed to evaluated the functions of HMOX1 in glioma. Third, we included 27 genes in our risk signature, proposing a great challenge for experimental validation. Finally, The TCGA dataset does not provide details of adjuvant therapy and, in particular, immunotherapy, which reduces the generalizability of our conclusions. Unfortunately, we do not have such a LGG dataset to study the immunotherapy response. Although we used patients from the CGGA database (validation set) to validate the results obtained from the TCGA database (training set) and showed good concordance, we believe that a larger multi-center cohort of glioma patients who undergo immunotherapy will provide insight into this issue.
In our study, we employed data from independent cohorts to explore the expression profile of heme oxygenase-1 in glioma. HMOX1 influences immune infiltration as well as survival prognosis of LGG patients. Importantly, HMOX1 were related to oncological functional states including angiogenesis, hypoxia, and metastasis. A risk signature and nomogram based on HRGs have robust predict power as well as the potential for clinical applications.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of Xiangya Hospital. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
WY performed the data analysis, interpreted the data, and prepared the draft. FL collected the samples in our cohort and they were responsible for the subsequent RNA sequencing of them. CL performed the visualization. ZL and FL revised the manuscript. CL and FL designed the research and supervised all the work. All authors have read and approved the final manuscript, and agree to be accountable for the content of the work.