Increased Expression of Homeobox 5 Predicts Poor Prognosis: A Potential Prognostic Biomarker for Glioma

Introduction

Gliomas are the most widespread and malignant type of primary brain tumors in human adults and have a mortality rate of approximately 30%.^1,2 Although surgical resection followed by combined chemo-radiotherapy 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.^3,4 Although several studies have investigated the molecular mechanisms of glioma malignancy and aimed to identify therapeutic targets,² the factors involved in promoting malignant proliferation and metastasis have not been completely elucidated. A comprehensive exploration of the oncology-related molecular mechanism might help to identify novel glioma-predictive markers.

The HOX gene family is a family of homeobox genes that encode transcription factors with key roles in both tumor development and malignancy.^5,6 There are 39 HOX genes in mammals, which are divided into four clusters named 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 specific DNA-binding domain. Several researchers have demonstrated that some members of the HOXB cluster are dysregulated in glioma tissue, which contributes to oncogenesis. For example, HOXB2, HOXB3, and HOXB9 have been reported to be upregulated in glioma tissue and were shown to promote the proliferation and migration of glioma cells.^7–9 Conversely, HOXB1 expression in glioma tissue is significantly downregulated and might be strongly associated with the degree of malignancy.¹⁰ In our previous study, we constructed an endogenous RNA network that might have an effect on the survival rate of glioblastoma patients; HOXB5 is a member of this network and serves as the binding target of miR-7.¹¹ 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 its potential prognostic value 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. Moreover, we validated HOXB5 protein expression in glioma cells and clinical samples with Western blotting. Furthermore, we investigated the prognostic value of HOXB5 and built a nomogram for risk quantification. The functional mechanism and tumor microenvironment were explored subsequently.

Materials and Methods

RNA-Sequencing Data and Processing

This study included 1709 glioma samples and 225 non-tumor tissue samples. Samples were summarized from five datasets as follows: TCGA Lower Grade Glioma (LGG); TCGA Glioblastoma (GBM) cohort; The Genotype Tissue Expression (GTEx) datasets; Chinese Glioma Genome Atlas (CGGA) mRNAseq_325 dataset; CGGA mRNA-array_301 dataset and REMBRANDT microarray dataset. We merged five non-tumor samples from TCGA and 200 normal brain cortex samples from GTEx into the normal tissue group of discovery sets for comparison. Moreover, the other three datasets were defined as external validation sets. Furthermore, we used CGGA_mRNA693 dataset as additional validation and put the corresponding results in Figure S1.

The gene expression data (HTSeq-counts and TPM) of TCGA and GTEx datasets 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 RNA-seq data, the normalized count reads from the pre-processed data were log₂ transformed after adding a 0.001 pseudo count (to avoid an infinite value upon log transformation). All clinical information was obtained from Gliovis (http://gliovis.bioinfo.cnio.es/). Using the median HOXB5 expression profile, the cases were divided into high and low HOXB5 expression groups. The overall survival (OS) was estimated from data of diagnosis to death or final follow-up. Unavailable or unknown clinical features were regarded as missing values, and the data are summarized in Table S1.

Patient Samples

Patient samples for Western blot were collected at the First Hospital of China Medical University and included seven samples (nine glioma tissues: two cases each for grades II, III, and IV, and one glioma peri-tumor tissue as a control). This study complied with the Declaration of Helsinki and was approved by the Ethics Committee of the First Hospital of China Medical University. All experiments were performed in compliance with the relevant regulations, and all patients provided written informed consent. The patient information are presented in Table S2.

Cell Culture

The human glioma cell line U251 and LN229 cells were purchased from the Chinese Academy of Sciences cell bank (Shanghai, China) in September 2017. U87 was purchased from Gene-Chem (Shanghai, China) in March 2018. HS683 cells were obtained as a gift from Professor Shaowu Ou (Department of Neurosurgery, First Affiliated Hospital of China Medical University). The human neuroblastoma cell line Sh-SY5Y was purchased from the American Type Culture Collection (ATCC, Manassas, VA, USA) in December 2018. Glioma cells were maintained in Dulbecco’s modified Eagle’s medium (DMEM, Bioind, Israel) containing 10% fetal bovine serum (Bioind, Israel) and 1% penicillin/streptomycin (Gibco, USA) at 37°C with 5% CO₂. Sh-SY5Y cells were cultured in DMEM/F12 (Gibco, USA).

Protein Extraction and Western Blotting

Total protein was isolated from the samples using whole cell lysis buffer (Beyotime). Next, 20 μg of protein from each sample was loaded into the lanes and subjected to SDS-PAGE, followed by PVDF membrane (Millipore) transfer. Then, the membranes were blocked with 5% skimmed milk for 1 hour at room temperature, followed by primary antibody incubation overnight at 4°C (HOXB5, 1:1000, Bioss, bs-6668R; or GAPDH, 1:2000, Bioss, bs-2188R). Peroxidase-conjugated Affinipure coated anti-rabbit IgG (Bioss; 1:5000) was employed as the secondary antibody. Protein bands were visualized with chemiluminescence ECL reagents (Beyotime) and quantified using Image J software.

Survival Analyses and Model Construction

The Kaplan–Meier method followed by a Log rank test was used to assess survival distribution between the high HOXB5 expression group and the low HOXB5 expression group.¹² Receiver operating characteristic (ROC) analysis was performed to estimate the accuracy in predicting 1-, 3-, and 5-year survival. Multivariate cox analysis was applied to investigate independent risk factors. Schoenfeld’s Residuals Test analyses were used to test the assumptions of proportionality of the hazard ratios. Variance inflation factors were tested to ensure that all features had no multicollinearity. A nomogram was generated with the R package ‘rms’. Calibration plots and the concordance index (C-index) were used to evaluate the performance of the model.

Bioinformatics Analysis

The ClusterProfiler package was applied to perform GO and KEGG functional enrichment analysis on HOXB5-correlated genes. Gene set enrichment analysis (GSEA) was performed to identify the enrichment of specific gene sets from MSigDB (http://www.gsea-msigdb.org/gsea/msigdb/index.jsp).¹³ The criteria were |NES > 1.5|, p-value < 0.05, and FDR < 0.05. The stromal score and glioma purity were computed using the ESTIMATE R package.¹⁴ Endothelial cell infiltration scores were calculated with Microenvironment Cell Populations (MCP)-counter and estimate the proportion of immune and cancer (EPIC) package.^15,16

Statistical Analysis

All statistical analyses and plots were performed using the R v4.02 software. The Wilcoxon rank sum test was used to assess differences in gene expression, and Pearson correlation analysis was used to calculate correlations. Permutation tests were performed with corrections for Western blot results of clinical samples. Other statistical computations and figure generating were performed using R packages (ggplots2, corrplot, coin, and pheatmap). All reported P-values were two-sided and considered significant for values < 0.05.

Results

Expression of HOXB5 in Peri-Tumor and LGG and GBM Tumors

Based on RNA-sequencing (RNA-seq) data from TCGA and GTEx, we discovered that HOXB5 expression was elevated in gliomas compared to that in non-tumor brain tissues (Figure 1A; p < 0.001). Moreover, the expression level gradually increased with WHO grade, with Grade IV GBM having the highest level (Figure 1B and C; p < 0.001). We also confirmed this observation with CGGA_325 data (Figure 1D–F; p < 0.001).

Figure 1 Expression of HOXB5 in glioma from TCGA and CGGA datasets. (A) HOXB5 expression level between glioma and normal brain samples. (B) HOXB5 expression level between lower grade glioma (LGG) and glioblastoma (GBM). (C) HOXB5 expression level with different grades. (D–F) Validation in CGGA data (Wilcoxon, *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001).

HOXB5 Protein Expression Based on Western Blot Assays

Next, we proved the protein expression pattern of HOXB5 in glioma samples and several central nervous system cancer cell lines. As shown in (Figure 2A and B), six glioma samples showed the upregulation of HOXB5 expression compared with levels in one peri-tumor normal brain tissue (p = 0.02). In addition, we found that GBM U251 presented with the highest expression of HOXB5 protein, whereas HS683 had the least (p < 0.001, Figure 2C and D). These outcomes were consistent with the preceding results. Collectively, these data suggest that HOXB5 protein is overexpressed in aggressive glioma.

Figure 2 Western blot results of HOXB5 protein. (A) Western blot result of HOXB5 levels in central nervous system tumors and cells and quantification of expression (B). (C) Western blot analysis of HOXB5 in the indicated peri-tumor and glioma samples and quantification of expression (D). Peri-tumor, n = 1; grade II, n = 2; grade III, n = 2 grade IV, n = 2; Wilcoxon permutation test, *p < 0.05; **p < 0.01).

Association Between HOXB5 Protein Expression and Glioma Prognosis

As shown (Figure 3A–D), 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 based on four datasets (p < 0.0001). The time-dependent ROC (tROC) curves indicated that the HOXB5 expression profile had moderate accuracy in predicting 1-,3-, and 5- year survival (Figure 3E). Furthermore, multivariate Cox regression analysis indicated that HOXB5 protein expression level was an independent prognostic risk factor in glioma (HR = 1.12, p = 0.01; Table 1), which should arouse the attention of researchers.

Table 1 Multivariate Analysis of Clinical Prognostic Parameters in TCGA

Figure 3 Survival analysis and risk model with HOXB5. (A–D) Kaplan–Meier survival plots for HOXB5 in four cohorts. (E) Respective time-ROC (receiver operating characteristic) curves of each survival analysis for 1, 3, and 5 years. (F) Nomogram combined with HOXB5 and clinical characteristics, (G) calibration curves (1-, 3-, and 5-year survival), and C-index of nomogram.

Construction and Evaluation of a HOXB5-Associated Nomogram

To better quantify the risk assessment and survival probability for glioma patients, we extracted independent prognostic factors and built a nomogram (including age, WHO grade, and HOXB5 expression level (z-score normalized); Figure 3F). Then, we also verified the prediction accuracy of the line chart. The bias-corrected lines of 1-, 3-, and 5-year survival probability in the calibration plot were found to be quite close to the ideal curve (the 45-degree line). The C-index for the nomogram was 0.827 with 1000 bootstrap replicates (Figure 3G). All these results indicated that the nomogram is a promising model of survival prediction for glioma patients.

GO/KEGG Enrichment and GSEA Analysis of HOXB5 in Gliomas

To comprehensively explore the molecular mechanism underlying the effects of HOXB5 in glioma, we investigated 545 genes that were strongly correlated with this marker (Pearson |r| ≥ 0.4; Figure 4A). These genes were enriched in the regulation of cell growth, endothelial cell proliferation, and JAK−STAT signaling pathway, which are involved in tumor progression (Figure 4B–C). GSEA was used to verify the biological function of HOXB5, showing that the high expression group had an activated angiogenesis and IL6/JAK-STAT3 pathway phenotype (Figure 4D–F). These findings indicated that HOXB5 is potentially involved in several tumor-promoting processes.

Figure 4 Correlated genes and functional enrichments of HOXB5. (A) The global HOXB5-correlated genes were identified by a Pearson test. (B and C) GO and KEGG enrichment analysis. (D–F) GSEA analysis of highly expressed HOXB5 group to the low expression group.

Glioma Microenvironment Status with HOXB5

Next, we investigated association between HOXB5 and genes involved in angiogenesis among the four datasets. As shown in Figure 5A and Table S3, serval members, including VGEFA, TIPM1, and S100A4, exhibited positive correlations with HOXB5. Furthermore, we examined the microenvironment infiltration status. With ESTIMATE algorithms, we discovered that HOXB5 was positively correlated with the stromal score (r = 0.471) but was negatively correlated with glioma purity (r = −0.428; Figure 5B and C). Additionally, MCP-counter and EPIC algorithms both showed that high expression HOXB5 samples exhibited more endothelial cell infiltration (Figure 5D and E). Together with the preceding results, we indicated that HOXB5 might promote angiogenesis.

Figure 5 Microenvironmental enrichment and HOXB5 expression. (A) Genes involved in angiogenesis correlated with HOXB5 in four cohorts. (B) HOXB5 correlation with stromal score, (C) tumor purity, and (D–E) endothelial cell infiltration level.

Discussion

Glioma represents the most common primary malignancy occurring in the brain, and it has a high mortality rate. Several biomarkers have been reported to be associated with the malignant phenotype of glioma and other disease.^17,26,27 However, these are still not adequate. In our study, we comprehensively explored the expression and prognostic role, as well as the functional mechanism, of HOXB5, according to public datasets. Moreover, we validated the protein expression of HOXB5 in glioma cells and clinical tumor samples and respective peri-tumor samples with Western blot experiments. Our analyses proved that HOXB5 shows higher expression in glioma and could serve as a promising prognostic biomarker.

There has been little research on the expression pattern of HOX genes. Under normal conditions, HOX genes regulate the vertebrate central nervous system development at specific time points.¹⁸ This pattern is called “spatial and temporal specificity”. An aberrant expression pattern is generally found in poorly differentiated samples and is associated with oncogenic effects.¹⁹ However, the expression pattern of HOXB5 in glioma remained to be elucidated. In our study, RNA-seq data from TCGA and CGGA, with protein data from clinical samples, all demonstrated that the expression level of HOXB5 was higher in glioma than in normal brain tissue and that it was increased with increasing grades of glioma.

HOXB5 was reported to be overexpressed and promote poor prognosis in prostate cancer and acute myeloid leukemia.^20,21 Furthermore, we comprehensively invested the prognostic value of HOXB5 as a biomarker. Results of four independent cohorts showed that high HOXB5 expression was highly related to an unfavorable OS. Multivariate Cox analysis further proved that HOXB5 is an independent risk factor. Additionally, a nomogram combining HOXB5 and other meaningful clinical features was constructed. Our model showed accuracy and consistency with predicted values for the 1- and 3-year OS. Considering all these results together, we believe that HOXB5 might serve as a promising biomarker for glioma patients.

Angiogenesis is one of the characteristics of malignant glioma, contributing to tumor proliferation and an unfavorable prognosis.²² Activation of the inflammatory response and IL6/JAK-STAT3 signaling pathway has been shown to play important roles in promoting this process.^23,24 In the present study, multiple functional analysis of HOXB5-correlated genes revealed a significant association between HOXB5 and the regulation of angiogenesis, inflammatory response, and the IL6/JAK-STAT3 pathway. It is worth noting that Fessner et al reported that the overexpression of HOXB5 can enhance blood vessel remodeling, as well as leukocyte infiltration, in vivo by upregulating IL6.²⁵ In our research, HOXB5 was correlated with endothelial cell proliferation and angiogenesis. Consistently, the population of endothelial cells recruited to the glioma microenvironment was positively associated with HOXB5 expression. Hence, based on these findings, we could conceivably hypothesize that HOXB5 might participate in the regulation of neovascularization. These observations might explain why glioma patients overexpressing HOXB5 have aggressive phenotypes and devastating outcomes. However, further studies are required to better understand the mechanisms through which HOXB5 regulates the glioma microenvironment.

Although our research improved our understanding on the relationship between HOXB5 and glioma, there were a few limitations worth mentioning. First, owing to the limitations of storage conditions, although we validated the expression with cells and clinical samples, there were only seven samples included in the Western blot experiment. Thus, additional glioma samples are required to confirm our results in the future. Furthermore, the study of the direct mechanisms of the involvement of HOXB5 in the development of gliomas requires further functional experimental evaluation.

Conclusions

In summary, the expression level of HOXB5 was evaluated comparing glioma and normal brain samples, and it was found to be generally increased with WHO grade. Moreover, HOXB5 might serve as a promising prognostic biomarker for glioma patients. A nomogram combined with HOXB5 could be beneficial for clinicians in clinical practice. Further, HOXB5 might be associated with endothelial cell growth and could promote angiogenesis.