Identification of HOXD10 as a Marker of Poor Prognosis in Glioblastoma Multiforme

Introduction

Glioblastoma multiforme (GBM) is the most common and aggressive primary intracranial tumor in adults, with a high fatality rate. The median survival time of GBM is only 14–15 months, with a 10% probability of 5-year survival.¹ Improving the survival time of GBM patients remains a difficult problem.

HOX genes encode proteins that function as critical master regulatory transcription factors during embryogenesis and also play key roles in the development of cancers.² In humans, 39 HOX genes are organized into four clusters (HOX-A, HOX-B, HOX-C, and HOX-D) on different chromosomes. At present, the expression of HOX family genes in a variety of tumors has been reported, such as lung cancer,^3–5 breast cancer,^6–8 ovarian cancer,^9–12 prostate cancer,^13–16 colorectal carcinoma^17–20 and glioma.^21–23

Although deregulated HOX gene expression was originally associated with oncogenic activities, it is now apparent that HOX genes might be lost as well as gained in cancer, and they may function as tumor suppressors as well as tumor promoters.²⁴ The aberrant expression of specific HOX genes tends to differ based on the tissue type and tumor site.²⁵ HOXD10, a tumor suppressor gene, is expressed at a low level in most tumors compared with normal tissues, such as in cholangiocellular carcinoma,²⁶ gastric carcinogenesis,²⁷ prostate cancer²⁸ and hepatocellular carcinoma.²⁹ However, the expression of HOXD10 in glioma is limited and controversial, especially in GBM. Fattah et al analyzed the expression profile of HOX family genes in neoplastic and non-neoplastic astrocytes using both cell lines and human brain tissues. That study showed that HOXD10 was silent or expressed at very low levels in nonneoplastic astrocytes, while they were over-expressed in malignant astrocytic cell lines and GBM specimens.³⁰ Additionally, Buccoliero et al used RT-PCR to examine the relative expression of HOX-D genes in a group of 14 pediatric low-grade gliomas, and the results reveled that HOXD10 expression was found in four (29%) out of 14 astrocytoma samples, whereas it was not found in any non-neoplastic brain tissues.³¹ In Chinese studies, HOXD10 is highly expressed in normal brain tissues and expressed at a low level in human glioma tissues, and it is negatively correlated with the malignant degree and prognosis of glioma.^32,33

Considering the inconsistency of the abovementioned research results and the small sample sizes included in these studies, with the aim to further clarify the significance of HOXD10 in GBM, we analyzed the mRNA expression of HOXD10 and its effect on prognosis using data derived from online databases such as GEPIA and The Cancer Genome Atlas (TCGA), the Chinese Glioma Genome Atlas (CGGA), and the Gene Expression Omnibus (GEO) (GSE74187 and GSE4412). Importantly, we examined the protein expression of HOXD10 in 71 GBM patients collected from our hospital, and their overall survival time (OS) and progression-free survival time (PFS) were analyzed. Additionally, based on the mRNA expression of HOXD10, GBM patients from TCGA were divided into a low-HOXD10-expression group and a high-HOXD10-expression group to analyze gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and construct a long non-coding RNA (lncRNA)-microRNA(miRNA)-mRNA network as well as a protein-protein interaction (PPI) network.

Materials and Methods

Databases

GEPIA is a newly developed interactive web server for analyzing the RNA sequencing expression data of tumors and normal samples from TCGA and the Genotype Tissue Expression (GTEx) projects using a standard processing pipeline, which was used to compare the mRNA expression of HOXD10 in GBM and normal tissues.

The mRNA-seq data in count format of the TCGA-GBM cohort was obtained from the TCGA dataset (https://portal.gdc.cancer.gov/). mRNA microarray data from 301 glioma patients and mRNA sequencing data from 693 glioma patients were obtained from the CGGA data portal (http://www.cgga.org.cn/). Patients histologically classified with primary GBM were screened for subsequent analyses. mRNA microarray data from GSE74187 and GSE4412 were obtained from GEO (https://www.ncbi.nlm.nih.gov/geo/).

At the same time, survival information matching the above datasets was obtained, and patients with a survival time greater than 30 days were used for survival analysis.

Clinicopathological Data Collection and Analysis

Clinical data such as age, gender and pathological results were retrospectively collected from 71 patients who underwent surgical resection and were pathologically confirmed with GBM from January 2016 to December 2018 at the Henan Provincial People’s Hospital. The patients included 43 males and 28 females aged between 10 and 81 years. The tumors of the patients were staged according to the 2016 WHO classification system. Survival time was calculated from the date of diagnosis to the end of follow-up. The patients were followed up until death or until the end date of the study (31 January 2021). The median follow-up was 11.4 months. The study was approved by the Ethics Committees of the Henan Provincial People’s Hospital and the ethical approval numbers was 2020-107. All the patients selected for our study were fully informed about our experiment protocols and signed an informed consent to participate in this study. The study was conducted in accordance with the Declaration of Helsinki.

Immunohistochemistry

All specimen tissues were formalin fixed, paraffin embedded and sectioned into 4-µm serial sections. Additionally, 15 normal brain tissue samples were assessed as the control group. The immunohistochemistry was performed according to our previously published article.³⁴ Briefly, after antigen retrieval (using a microwave at high power for 8 min, followed by middle-low power for 7 min in EDTA solution at pH 9.0), the sections were treated to quench endogenous peroxidase activity. Rabbit serum was used for blocking. The slides were then stained with HOXD10 antibody at 1:200 (ab85698, Abcam, Cambridge, UK) overnight at 4°C and secondary antibodies were incubated for 1 h at room temperature. After adding substrate and hematoxylin staining, the slides were covered and observed by microscope. The proportion of stained cells in each field was assessed as follows: 0% = 0; 1–24% = 1; 25–49% = 2; 50–74% = 3; and 75–100% = 4. The intensity of overall staining was graded as follows: no staining = 0; weak staining = 1; moderate staining = 2; and strong staining = 3. The overall immunostaining score for each specimen was calculated by multiplying the percentage score by the intensity score.³⁵ A final score of 0–3 indicated negative expression and 4–12 indicated positive expression.

Differentially Expressed Gene Analysis and Functional Enrichment Analysis

With the aim to understand the role of HOXD10 in GBM, we performed the functional enrichment analysis. Differentially expressed genes (DEGs) were identified with the “DEseq2 package” in R software. Fold change > 1.0 and adjusted P < 0.05 were set as the cutoffs. A volcano plot of the DEGs was drawn using the “ggplot2” package. Then, we transformed the gene symbols into Entrez IDs via the “org.Hs.eg.db” package and performed pathway enrichment analysis for the DEGs based on the GO database and KEGG by using the “cluster Profiler”, “ggplot2” and “enrich plot” packages. The P-values were adjusted by the False Discovery Rate (FDR) method and an FDR < 0.25 was considered significantly enriched pathways.

Construction of a lncRNA-miRNA-mRNA (ceRNA) Network

The “survival” package was used to identify six lncRNAs that are related to the prognosis of GBM. Relevant lncRNAs-target data (miRNAs) were obtained from miRCode. Relevant miRNA-target data (mRNAs) were obtained from miRTarBase. Then, a lncRNA-miRNA-mRNA network was constructed according to the hypothesis of ceRNA.³⁶ In total, we obtained 89 nodes (45 mRNAs, 39 miRNAs and 5 lncRNAs) and 227 edges to construct the lncRNA-miRNA-mRNA network. Nodes and edge files were exported and visualization was performed with Cytoscape 3.6.1.

PPI Network Analysis

The interaction of 138 HOXD10 correlated differentially expressed mRNAs (DEmRNAs) were analyzed using the online tool STRING (https://www.string-db.org/). A PPI network was constructed by setting medium confidence at 0.400 and excluding irrelevant genes. The network was exported and visualized using Cytoscape 3.6.1.

Statistical Analysis

Differences between groups were analyzed using the χ² test or Student’s t-test depending on the data type. Kaplan-Meier plots and Log rank tests were used for survival analysis. Multivariate analyses were based on the Cox proportional hazards regression model. All statistical analyses were performed in RStudio (version 3.6.8), and we considered a P-value < 0.05 to be statistically significant.

Results

High mRNA Expression of HOXD10 in GBM Predicted Poor Prognosis

The mRNA expression of HOXD10 in GBM and normal tissues was analyzed using GEPIA. This analysis revealed that the HOXD10 expression was much higher in GBM compared to normal tissues (Figure 1A). Based on the database analysis of TCGA, CGGA, GSE74187 and GSE4412, high HOXD10 expression was associated with shorter survival time than those with low HOXD10 expression in GBM (Figure 1B–F).

Figure 1 (A) The mRNA expression of HOXD10 is significantly upregulated in GBM compared with normal brain tissues. *P < 0.05. (B–F) High mRNA expression of HOXD10 in GBM predicted poor prognosis based on database analyses. (B), TCGA database; (C and D), CGGA database; (E and F), GEO database. Abbreviations: GBM, glioblastoma multiforme; TCGA, The Cancer Genome Atlas; CGGA, Chinese Glioma Genome Atlas; GEO, Gene Expression Omnibus.

Immunohistochemical analysis of HOXD10 protein expression and the association with clinicopathological characteristics.

Then, we collected and followed up 71 GBM patients and 15 normal brain tissues from the Henan Provincial People’s Hospital. IHC was used to examine the protein expression of HOXD10 in GBM. Immunostaining of HOXD10 was observed in the nucleus of cancer cells (Figure 2A–D), and the positive staining rate in GBM was 26.80% (19/71), while the positive rate in normal brain tissues was 73.3% (11/15). The expression of HOXD10 was much higher in normal brain tissue than that in GBM samples (χ² = 11.825, P = 0.001). The association between HOXD10 expression and the clinicopathological characteristics of GBM patients are listed in Table 1. The results revealed that the expression of HOXD10 was significantly associated with the status of IDH1.

Table 1 Associations Between HOXD10 Protein Expression and GBM Patient Clinicopathological Characteristics

Figure 2 Representative immunohistochemical staining of HOXD10 in GBM and normal brain tissues. (A), HOXD10 was negative in normal brain tissue (X400). (B), HOXD10 was positive in normal brain tissue (X400). (C), HOXD10 was negative in GBM (X400). (D), HOXD10 was positive in GBM (X400). Abbreviation: GBM, glioblastoma multiforme.

Univariate and Multivariate Survival Analyses

Survival analysis showed that complete surgical resection, postoperative radiotherapy, postoperative temozolomide (TMZ) chemotherapy, IDH1 mutation, and low HOXD10 expression were beneficial prognostic factors (Figure 3A–E). Further multivariate Cox regression analysis revealed that complete surgical resection and postoperative radiotherapy were independent prognostic factors (Table 2). Furthermore, GBM patients who received complete surgical resection, postoperative radiation therapy and had only one lobe involved at the initial treatment had a longer median PFS (mPFS) (Figure 3F–H), and these three factors were independent predictive factors (Table 3).

Table 2 Multivariate Analysis of Prognosis for GBM (N=71)

Table 3 Multivariate Analysis of PFS for GBM (N=71)

Figure 3 Univariate and multivariate survival analyses. (A), Complete surgical resection can improve the survival time of GBM patients. (B), Postoperative TMZ chemotherapy can improve the survival time of GBM patients. (C), Postoperative radiation therapy can improve the survival time of GBM patients. (D), GBM patients with IDH1 mutation had a longer survival time. (E), GBM patients with high HOXD10 protein expression had a shorter survival time. (F), Complete surgical resection can improve the PFS of GBM patients. (G), Tumors that invade multiple lobes were more likely to recur after surgery. (H), Postoperative radiation therapy can prolong the PFS of GBM patients. Abbreviations: GBM, glioblastoma multiforme; TMZ, Temozolomide; mOS, median overall survival time; mPFS, median progression-free survival time.

DEG Profiles and Functional Enrichment Analysis

To understand the biological functions related to the DEGs, we divided GBM patients from TCGA database into low and high expression groups based on the median value of HOXD10 expression. A total of 138 DEGs were identified (Table S1; Figure 4A). GO contains biological process (BP), cellular components (CC) and molecular function (MF). There was no item enriched in CC process and the enriched items in BP and MF were shown in Figure 4B. KEGG enrichment analyses showed that the DEGs based on HOXD10 expression were mainly involved in cytokine-cytokine receptor interactions, viral protein interactions with cytokines and cytokine receptors and neuroactive ligand-receptor interaction signaling pathways (Figure 4C).

Figure 4 Differential analysis of gene expression profiles in high- and low-HOXD10 expression groups and pathway enrichment analysis. (A), Volcano plot showing the differentially expressed genes. (B), Bubble chart of the results of GO enrichment analysis. (C), Bar plot of the results of KEGG enrichment analysis. Abbreviations: FDR, false discovery rate; BP, biological process; MF, molecular function.

ceRNA Network

The “survival” package was used to identify six survival-related lncRNAs (Figure 5A–F). LINC01539 was removed because there was no suitable miRNA match. Based on the expression level of each lncRNA, the ceRNA network was divided into a lncRNA high expression group and a lncRNA low expression group. In the up-regulated ceRNA network, there were four lncRNAs, 37 miRNAs and 42 mRNAs (Figure 6A). In the down-regulated ceRNA network, there was one lncRNA, three miRNAs and three mRNAs (Figure 6B).

Figure 5 Survival analysis of differentially expressed lncRNAs. (A–D and F), High expression of HOTAIR, HOXC-AS3, HOXC12-AS, LINC01305 and LINC01539 were poor prognostic factors. (E), High expression of MIR663AHG was good prognostic factor.

Figure 6 Construction of a ceRNA network. (A), The ceRNA network was constructed based on the up-regulation of lncRNAs screened by HOXD10 expression level and related to the survival of GBM patients. (B), The ceRNA network was constructed based on down-regulated lncRNAs screened by HOXD10 expression level and related to the survival of GBM patients. Note: Yellow, mRNA; Green, miRNA; Red, LncRNA.

PPI Network

To construct the PPI network, 138 HOXD10-correlated DEmRNAs were analyzed, and 71 genes were included in the network after the elimination of disconnected nodes. The network is shown in Figure 7. HOXD10 was closely related to other HOX family members, including HOXD8, HOXD9, HOXD11, HOXD13 and HOXB3.

Figure 7 Protein–protein interaction (PPI) network analysis. After eliminating disconnected nodes, 71 genes were involved in the PPI. The network was visualized with Cytoscape 3.6.1.

Discussion

Previous studies on HOXD10 in GBM are limited, and most have focused on gliomas of all grades. In the present study, to our knowledge, this is the first examination of HOXD10 expression specifically in GBM. We found that the protein expression of HOXD10 in GBM was much lower than that in normal brain tissues, contrary to the results of the online tool GEPIA. Sample size may go some way to explaining the difference. But, more importantly, we found it interesting that in the small number of studies by Chinese authors, the expression of HOXD10 in normal brain tissues was much higher than that in human glioma tissues.^32,33 While studies from the non-Chinese scholars, HOXD10 expression in glioma was slightly higher than that in non-neoplastic brain.^30,31 Whether racial difference is one reason for the inconsistency of our data and the online tool GEPIA remains to be further verified. Lastly, whether HOXD10 differs in the expression levels of RNA and protein remains to be further studied.

Both our own data and database-derived data suggest that GBM patients with high HOXD10 expression have significantly worse OS. However, a previous study of 40 glioma samples showed that high expression of HOXD10 was a factor that predicted better prognosis.³³ In addition to small sample sizes, the previous studies included not only GBM, but also other grades of glioma.

In fact, HOXD10 is recognized as a tumor suppressor and is known to be down-regulated in a variety of cancer types,^26–29,32 such as colon carcinomas,³⁷ early-stage head and neck squamous cell carcinoma (HNSCC)³⁸ and oral squamous cell carcinoma,³⁹ in which HOXD10 is up-regulated. In colon carcinomas, up-regulation of HOXD10 contributes to the crypt colonic stem cell overpopulation that drives colorectal cancer (CRC) development³⁷ and in HNSCC and cervical cancer, HOXD10 over-expression resulted in increased proliferation ability.^38,40 Based on our data analysis, although HOXD10 expression was low in GBM compared with normal brain tissue, patients with high HOXD10 expression have a worse prognosis. The underlying mechanism could be like that of TGF-β, a well-known bifunctional regulator that can either inhibit or stimulate cell proliferation.⁴¹ Typically, this cellular effect accompanies an induction of mesodermal cell fate and inhibition of neural cell differentiation. However, during pathological conditions, these defined effects of TGF-β can be reversed; for example, the growth inhibitory effect is replaced with its tumor promoting ability.⁴² Therefore, based on unique features with respect to its expression and function, HOXD10 is best described as a “tumor modulator” rather than an oncogene or tumor suppressor gene. HOXD10 may play different or even opposite roles in different stages of GBM occurrence and development. Further larger samples and deeper exploration need to be performed.

IDH1 mutation as a marker of prognosis has been widely recognized and it is regarded as a therapeutic target,⁴³ which was further confirmed in this study. Correlation analysis showed that HOXD10 expression was associated the status of IDH1, although the specific regulatory mechanism remains unclear.

GO enrichment analyses revealed the DEGs based on the expression of HOXD10 from TCGA database mainly acted as the DNA-binding transcription activators and involved in the particular tissue or organ development. In accordance with the previous studies that HOX proteins can act as transcription factors to enhance or repress gene expression during normal embryonic development.² HOX proteins are integral to normal temporospatial limb and organ development along the anterior-posterior axis.⁴⁴ KEGG analysis showed that the signaling pathways mainly involved in the cytokine-cytokine receptor interaction, which indicated that HOXD10 may regulate the inflammatory response and immune function in GBM.

In the ceRNA network, we found that HOTAIR and LINC01305 were enriched in a number of nodes. HOTAIR is one of the most extensively studied lncRNAs found to be dysregulated in human cancer.⁴⁵ Tan et al used qRT-PCR to measure the HOTAIR expression in serum from 43 GBM and 40 control samples, and the results revealed that HOTAIR expression was significantly correlated with high grade brain tumors and that serum HOTAIR can be used as a novel prognostic and diagnostic biomarker for GBM.⁴⁶ Previous research also indicated that HOTAIR promoted malignant progression and poor prognosis in glioma patients and exhibited pro-oncogenic activity.⁴⁷ LINC01305 has been widely explored in cervical cancer and is involved in the epithelial-mesenchymal transition by regulating the PI3K/AKT signaling pathway.^48,49 In the present study, the database analysis showed that LINC01305 was over-expressed and predicted poor prognosis in GBM. However, its clinical significance and prognostic value needed to be further verified in our tissue samples. Through gene enrichment and PPI network analysis of HOXD10-related genes, we found that HOXD10 is mainly involved in the signal transduction of cytokine interactions, and there was mutual regulation between HOXD10 and multiple cytokines, such as CCL20 and CXCL6. HOXD10 has been found to be involved in tumor proliferation, invasion or angiogenesis through multiple signaling pathways, such as the RHOC/AKT/MAPK and ERK signaling pathways.^{26,29,38,40,50,51} Whether HOXD10 interacts with cytokines to regulate the GBM tumor microenvironment is a future research direction.

Conclusion

Based on our experimental data, although HOXD10 expression is low in GBM compared with normal brain tissues, GBM patients with high HOXD10 expression have a worse prognosis. HOXD10 may play different or even opposite roles in different stages of GBM occurrence and development. For patients with GBM, HOXD10 may be a valid predictor of prognosis.