Increased expression of homeobox 5 predicts poor prognosis: A potential therapeutic target for glioma


 Background: The homeobox gene 5 (HOXB5) encodes a transcription factor that regulates the central nervous system embryonic development. Of note, its expression pattern and prognostic role in glioma remain unelucidated. This study aimed to identify the relationship between HOXB5 and glioma by investigating the HOXB5 expression data from the The Cancer Genome Atlas (TCGA) and The Genotype Tissue Expression (GTEx) databases and validating the obtained data using the Chinese Glioma Genome Atlas (CGGA) database. Kaplan-Meier and univariate cox regression analyses were performed to assess the prognostic value of HOXB5. The key functions and signaling pathways of HOXB5 were analyzed using GSEA and GSVA. Immune infiltration was calculated using Microenvironment Cell Populations-counter (MCP-counter), single-sample Gene Set Enrichment Analysis (ssGSEA), and ESTIMATE algorithms.Result: HOXB5 expression was elevated in glioma tissues. The increased levels of HOXB5 were significantly correlated with a higher WHO grade and aggressive cancer phenotypes. HOXB5 overexpression represented a risk factor that was associated with shorter overall survival (OS) while exhibiting a moderate forecast efficiency in most clinical subgroups. These results were validated using the CGGA and Rembrandt datasets. Furthermore, the functional analysis showed enrichment of angiogenesis, the IL6/JAK-STAT3 pathway, and inflammatory response in the tissues that showed high expression of HOXB5. Lastly, the high expression of HOXB5 was associated with enrichment of Tregs and MDSCs, and HOXB5 expression was shown to play a role in several immune checkpoint genes.Conclusions: HOXB5 may serve as a predictive factor of glioma malignancy and prognostic status and represents potential as a molecular treatment candidate.


Abstract
Background: The homeobox gene 5 (HOXB5) encodes a transcription factor that regulates the central nervous system embryonic development. Of note, its expression pattern and prognostic role in glioma remain unelucidated. This study aimed to identify the relationship between HOXB5 and glioma by investigating the HOXB5 expression data from the The Cancer Genome Atlas (TCGA) and The Genotype Tissue Expression (GTEx) databases and validating the obtained data using the Chinese Glioma Genome Atlas (CGGA) database. Kaplan-Meier and univariate cox regression analyses were performed to assess the prognostic value of HOXB5. The key functions and signaling pathways of HOXB5 were analyzed using GSEA and GSVA. Immune in ltration was calculated using Microenvironment Cell Populationscounter (MCP-counter), single-sample Gene Set Enrichment Analysis (ssGSEA), and ESTIMATE algorithms.
Result: HOXB5 expression was elevated in glioma tissues. The increased levels of HOXB5 were signi cantly correlated with a higher WHO grade and aggressive cancer phenotypes. HOXB5 overexpression represented a risk factor that was associated with shorter overall survival (OS) while exhibiting a moderate forecast e ciency in most clinical subgroups. These results were validated using the CGGA and Rembrandt datasets. Furthermore, the functional analysis showed enrichment of angiogenesis, the IL6/JAK-STAT3 pathway, and in ammatory response in the tissues that showed high expression of HOXB5. Lastly, the high expression of HOXB5 was associated with enrichment of Tregs and MDSCs, and HOXB5 expression was shown to play a role in several immune checkpoint genes.
Conclusions: HOXB5 may serve as a predictive factor of glioma malignancy and prognostic status and represents potential as a molecular treatment candidate.

Background
Gliomas are the most widespread and malignant type of primary brain tumors in human adults and have a mortality rate of approximately 30% [1]. While surgical resection followed by combined chemoradiotherapy represents the standard treatment, the prognosis is not encouraging. Following diagnosis, patients with glioblastoma (GBM, the most aggressive type of glioma) have a median survival rate of 14 months [2,3]. Although several studies have investigated the molecular mechanisms of glioma malignancy and aimed to identify therapeutic targets [4], the factors involved in promoting malignant proliferation and metastasis have not been completely elucidated. A comprehensive exploration of the oncology-related molecular mechanism of glioma may help us identify novel glioma-predictive markers.
The HOX gene family is one of the families of homeobox genes that encode transcription factors that play key roles in both tumor development and malignancy [5,6]. There are 39 HOX genes in mammals, which are divided into four clusters named as HOXA, HOXB, HOXC, and HOXD . The HOXB cluster consists   of 11 genes (including, HOXB1, HOXB2, HOXB3, HOXB4, HOXB5, HOXB6, HOXB7, HOXB8, HOXB9, HOXB10, and HOXB13), which encode nuclear proteins containing a speci c DNA-binding domain.
Several researchers have demonstrated that some members of the HOXB cluster are dysregulated in glioma tissue, which contributes to oncogenesis. For instance, HOXB2, HOXB3, HOXB7, and HOXB9 have been reported to be upregulated in glioma tissue and shown to promote proliferation and migration of glioma cells [7][8][9][10]. Conversely, HOXB1 expression in glioma tissue is signi cantly downregulated and it might be strongly associated with the degree of malignancy [11]. In our previous study, we constructed an endogenous RNA network that might have an impact on the survival rate of glioblastoma patients, and HOXB5 is a member of the network and serves as the binding target of miR-7 [12]. However, the clinical and prognostic value of HOXB5 in glioma remains unknown, and its functional role in glioma has never been investigated.
This study aimed to identify the relationship between HOXB5 and glioma, as well as the potential prognostic value of HOXB5 in glioma. For this, we analyzed the expression level of HOXB5 from glioma and normal brain tissue data from The Cancer Genome Atlas (TCGA) and The Genotype Tissue Expression (GTEx) databases. Furthermore, we performed a survival analysis to identify the prognostic value of HOXB5 and validated its effectiveness in three datasets. Lastly, the biological functions of HOXB5 involved in glioma were investigated using multiple functional enrichment analyses.

Association between HOXB5 expression and clinicopathologic variables
Based on RNA sequencing (RNA-seq) data from TCGA, GTEx, and Chinese Glioma Genome Atlas (CGGA) databases, we discovered that HOXB5 expression was elevated in gliomas compared to that in non-tumor brain tissues ( Fig. 1a-b). Moreover, a high level of HOXB5 expression was signi cantly correlated with a high WHO grade, and according to histopathologic classi cations, HOXB5 was the most enriched in GBM ( Fig. 1f-g). We further explored the expression pattern of HOXB5 according to ve molecular features that previous reported [13][14][15]. As shown in Fig. 1c-e and g-h, increased HOXB5 expression was signi cantly correlated with IDH-1 wildtype (WT), unmethylated MGMT promoter, and non-1p19q codeletion. Moreover, mesenchymal gliomas, which are the most malignant glioma subtype, exhibited the most HOXB5 enrichment among gliomas. We obtained similar results when using the validation set (Fig s1) for analysis validation. These results suggested that HOXB5 might be speci cally enriched in gliomas with aggressive phenotypes.
As shown in Fig. 2b-c, the Kaplan-Meier survival analysis suggested that patients expressing a high level of HOXB5 had a worse prognosis than those expressing a low level of HOXB5. The receiver operating characteristic (ROC) curves indicated that the HOXB5 expression pro le had a moderate accuracy in predicting survival (Fig. 2d-e). To better understand the mechanism underlying the action of HOXB5 in glioma, we investigated the relationship between HOXB5 expression and clinicopathologic variables using univariate Cox analysis (Fig. 2a). While IDH-1 mutant, 1p19q codeletion, and young patients (≤ 42 years of age) did not show an association with high HOXB5 expression, the overexpression of HOXB5 represented a risk factor (HR > 1 and P-value < 0.05) in most subgroups. Consistent with this, we discovered analog data in the validation set (Table s4). These results indicated that HOXB5 expression levels can differently impact the prognosis of patients with different types of gliomas.
Identi cation of genes co-expressed with HOXB5 and their functional annotation To explore differential biological features between the groups with low and high HOXB5 expression, we performed a differential expression analysis using DESeq2 R package, followed by a Pearson correlation analysis. Based on the cutoff criteria (|log2 Foldchange| > 1.5, |r > 0.4|, adjusted P-value < 0.05), a total of 545 corelated differential expressed genes (co-DEGs) were detected (Table s3). Next, gene ontology (GO) enrichment analysis was performed using Metascape. While the co-DEGs engaged in several terms, we found that blood vessel development, wounding response, leukocyte migration, and growth factor binding were signi cantly in uenced by these genes (Fig. 3a). To capture the association between the enriched terms, we selected a subset of terms and built a network plot. The edges linked the terms with a similarity > 0.3. Each node stood for an enriched term colored by the cluster-ID (Fig. 3b).

Biological Pathways
To comprehensively explore the molecular mechanism of action of HOXB5 in glioma, gene set enrichment analysis (GSEA) of the differences between the groups with low and high HOXB5 expression was performed to identify the key biological processes and signaling pathways in which HOXB5 is involved. The thresholds used were |NES| > 1.5, P-value < 0.05, and False Discovery Rate (FDR) < 0.01. We also used the gene set variation (GSVA)[16] method to screen the variation in biological processes in the glioma patient data and identi ed their associations with HOXB5 expression using Pearson correlation and differential analyses (|r| >0.4, |log2Foldchange > 0.4| and P-value < 0.05). HOXB5 was positively correlated with the in ammatory response, angiogenesis, and IL6/JAK-STAT3 pathway (Fig. 3c-f). Combined with the previous function annotation results, our ndings indicate that HOXB5 is potentially involved in several tumor-promoting processes and immune response processes.

Microenvironment Immune Function
The tumor microenvironment of glioma contains noncancerous cell types including stromal and immune cells, which contribute actively to the regulation of tumor progression by cross talking with cancer cells in the microenvironment [17,18]. By using ESTIMATE, Microenvironment Cell Populations-counter (MCP-counter) [19,20], and ssGSEA algorithms [21], we discovered that HOXB5 was positively corelated with the stromal score (r = 0.471) and immune score (r = 0.386), while being negatively correlated with glioma purity (r = − 0.428) (Fig. 4a-e). Additionally, immune cell in ltration calculation showed that HOXB5 was positively associated with enrichments of Tregs, MDSCs ( Fig. 4f-g), and endothelial cells (Table s5). Consistently, GSVA scores for angiogenesis, in ammatory response, and the IL6/JAK-STAT3 pathway shared the same correlation tendency with HOXB5. These results indicate that gliomas expressing high levels of HOXB5 recruit more immune cells and promote angiogenesis.
Immune Checkpoint Blockade Therapy Immune checkpoint blockade therapy provides therapeutic targets for immune cells and has shown potential bene ts for cancer treatment [22]. We examined the association of HOXB5 and several classical immune checkpoint genes using TCGA and CGGA datasets. As illustrated in Fig. 4h

Discussion
Previous studies have shown that HOXB5 is involved in embryo development, immune cell differentiation, and vascular remodeling [23,24]. Aberrant HOXB5 expression has been demonstrated to have carcinogenic effects and to promote the development of breast cancer and pancreatic cancer [25,26].
However, the expression pattern of HOXB5 and its potential prognostic impact on glioma remains to be elucidated. In our study, high throughput RNA-seq data from TCGA and CGGA databases demonstrated that the HOXB5 expression was upregulated in glioma tissues and was associated with a series of clinicopathologic characteristics including WHO grade, histological type, and molecular type. Using strati ed survival analysis, we established that a high HOXB5 expression represents an important unfavorable factor for the prognosis of glioma patients. Furthermore, multiple functional analyses indicated that angiogenesis, in ammatory response, and the IL6/JAK-STAT3 signaling pathway were signi cantly enriched in samples expressing high levels of HOXB5. Finally, our study showed that several microenvironment-in ltrated immune as well as stromal cells and immune checkpoint markers were signi cantly correlated with HOXB5 expression in glioma. These ndings suggest that HOXB5 may serve as a potential indicator of malignancy and as a prognostic marker for glioma patients.
Under normal conditions, HOX genes regulate the vertebrate central nervous system development at speci c time points [27]. This pattern is called 'spatial and temporal speci city'. An aberrant expression pattern is generally found in poorly differentiated samples and is associated with oncogenic effects [28]. Our study showed that compared to that in normal brain tissues, the expression of HOXB5 was signi cantly higher in glioma tissues. Furthermore, overexpression of HOXB5 was signi cantly enriched in patients with aggressive features including WHO Grade IV, IDH-wildtype, unmethylated MGMT promoter, and GBM. Moreover, mesenchymal GBM, the most aggressive and least differentiated subtype of glioma, exhibited the highest HOXB5 levels amongst the subtypes. These ndings demonstrate that glioma tissues with elevated expression of HOXB5 have malignant pathological phenotypes. In addition, the Kaplan-Meier survival analysis showed that high HOXB5 expression was highly relevant with unfavorable overall survival (OS) among glioma patients. The strati ed univariate Cox regression analysis also indicated that HOXB5 expression remains a powerful forecaster of most of the glioma subgroup prognoses. We obtained similar results following the analysis of the validation datasets. Collectively, our analyses indicate that HOXB5 is an oncogenic gene and may serve as a promising biomarker for predicting the malignancy and prognosis.
Several studies have shown that both glioma purity and the in ltrated immune cell components have clinical and molecular implications [29,30]. While lower purity represents aggressive progression, worse prognosis, and overloaded immune activity, the in ltration of immune cells is associated with the biological behavior and survival of glioma patients. We found that HOXB5 was signi cantly corelated with the immune score, stromal score, and glioma purity. This indicated that the tissue microenvironment with high HOXB5 expression has more variations and complexities. Consistently, positive correlations were discovered between HOXB5 and MDSCs as well as regulatory T cells. MDSCs can mediate tumorinduced immune tolerance and secrete IL-6 [31]. Simultaneously, Treg cells can suppress the effector T cell responses through a variety of mechanisms [32]. High enrichment of MDSCs and Treg cells in glioma often indicates an immune suppressive microenvironment that promotes tumor progression and facilitates tumor immunity escape. These observations may explain why glioma patients overexpressing HOXB5 have aggressive phenotypes and devastating outcomes.
Immune checkpoint blockade is emerging as a promising strategy for neoplasm treatment. PD-1 as well as CTLA4, LAG3, and TIM3 have been investigated in a clinical setting as potential targets for cancer treatment [33]. We identi ed signi cant correlations between HOXB5 and most of these targets. Therefore, it can be hypothesized that glioma patients with an increased expression of HOXB5 might bene t more from immune checkpoint blockade treatments.
Angiogenesis is one of the characteristics of malignant glioma, contributing to tumor proliferation and unfavorable prognosis [34]. The activation of the in ammatory response and IL6/JAK-STAT3 signaling pathway has been shown to play important roles in promoting this process [35,36]. In the present study, the multiple functional analysis of HOXB5-corelated genes revealed a signi cant association between HOXB5 and the regulation of angiogenesis, in ammatory response, and the IL6/JAK-STAT3 pathway. It is worth noting that Fessner et al. have reported that the overexpression of HOXB5 can enhance blood vessel remodeling as well as leucocyte in ltration in vivo by upregulating IL6 [23]. Consistently, the population of endothelial cells recruited by the glioma microenvironment was positively associated with HOXB5 expression (r = 0.24). Hence, based on these ndings, we could conceivably hypothesize that HOXB5 might participate in the regulation of the in ammatory response and neovascularization.
Considering the role of HOXB5 in the prognosis of glioma, it may serve as a promising potential target for developing new therapeutic strategies targeting the angiogenic process in glioma. However, further studies are required to better understand the underlying mechanisms of HOXB5 through which it regulates the glioma microenvironment.
Although our research improved our understanding on the relationship between HOXB5 and glioma, there were a few limitations worth mentioning. The current study was performed by using RNA-seq data from TCGA database; therefore, the HOXB5 mRNA expression levels should be veri ed using cell lines and clinical samples. Secondly, due to the limitations of public databases, several clinical factors, such as details regarding the treatments of the patients, remained uncovered. Furthermore, the study of the direct mechanisms of HOXB5 involved in the development of gliomas requires further functional experimental evaluation.

Conclusions
In summary, by exploring glioma patient data from multiple databases, our study revealed the expression pattern, prognostic value, as well as functional mechanisms of HOXB5 in glioma for the rst time. The ndings clearly suggest that HOXB5 is a pro-tumorigenic factor in glioma and has the potential of predicting malignancy, OS time, and the bene t of immune checkpoint treatment. Second, tissues expressing high levels of HOXB5 exhibited signi cantly enriched immune-suppressed microenvironments, which might promote tumor progression and an ineffective immune response evasion. Finally, HOXB5 might be involved in regulating the IL6/JAK-STAT3 pathway and local angiogenic processes. Taken together, our results revealed HOXB5 as a potentially reliable forecaster of aggressive phenotype and overall survival in glioma along with its potential as a target gene candidate for novel therapeutic treatments.

Patient samples
The gene expression data (HTSeq-counts and TPM) of TCGA datasets (Tumor: 697, Normal: 5) and GTEx datasets (Normal: 200) were downloaded from the UCSC Xena database ( http://xena.ucsc.edu/ ). The normalized expression matrices of the CGGA microarray cohort, CGGA RNA-seq cohort, and REMBRANDT microarray cohort were obtained from the CGGA database ( http://www.cgga.org.cn/ ). For CGGA RNAseq data, the normalized count reads from the pre-processed data were log2 transformed after adding a 0.5 pseudo count (to avoid in nite value upon log transformation). All clinical information was obtained from Gliovis ( http://gliovis.bioinfo.cnio.es/ ). Using the median HOXB5 expression pro le, the cases were divided into high and low HOXB5 expression groups. The OS was estimated from data of diagnosis to death or nal follow-up. Unavailable or unknown clinical features were regarded as missing values, and the data are summarized in Table S1. Bioinformatics Analysis Expression pro les (HTSeq-counts) were compared between the high and low HOXB5 expression groups to identify DEGs using DESeq2 R package, and the thresholds were |log2-foldchange (FC)| > 1.5 and adjusted P-value < 0.05. Metascape ( https://metascape.org/gp/index.html ) was used to perform GO (Gene Ontology) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses on the DEGs. GSEA and GSVA were performed to identify the enrichment of speci c gene sets from MSigDB ( http://www.gsea-msigdb.org/gsea/msigdb/index.jsp ). Stromal and immune score as well as glioma purity were computed using the ESTIMATE R package, as previously described. The immune in ltration analysis was done using the ssGSEA method and the GSVA package. Based on the signature genes of 28 types of immunocytes, the MCP-counter was used to estimate eight classical immune cells and two types of stromal cells.

Statistical analysis
All statistical analysis and plots were conducted using the R v4.02 software. Wilcoxon rank sum test was used to assess differences in gene expression, and Pearson correlation analysis was used to calculate correlations. The Kruskal-Wallis and Wilcoxon signed-rank tests were used to evaluate the relationships between clinicopathologic features and HOXB5 expression. The Kaplan-Meier method was used to evaluate the survival distribution between the high HOXB5 expression group and the low HOXB5 expression group. A log-rank test was used to estimate the difference between the high and low HOXB5 expression groups. The Cox regression analysis was used to estimate the prognostic value of HOXB5 in the different clinical subgroups. ROC analysis was performed using the timeROC package. Other statistical computations and gures were performed using the R packages (ggplots2, corrplot, and pheatmap). All reported P-values were two-sided and considered signi cant for values < 0.05. Declarations Ethics approval and consent to participate All data used in this study were obtained from TCGA and CGGA, and hence ethics approval and informed consent were not required.

Consent for publication
Not applicable.

Competing interests
The authors have declared that no competing interest exists.

Funding
The project is supported by the Science and Technology Project of Shenyang (18-014-4-03), the Science and Technology Project of the Education Department of Liaoning Province (LFWK201705) and Liaoning BaiQianWan Talents Program.
Authors' contributions GL selected the research topic and conducted the guidance of the process of the topic. CX wrote the article, performed data processing, and produced the gure. JH produced parts of tables. LL and YY contributed technology support of the R software and contributed parts of the manuscript. All authors read and approved the nal manuscript.