HOXA5 Is Recognized as a Prognostic-Related Biomarker and Promotes Glioma Progression Through Affecting Cell Cycle

Glioma is malignant tumor derives from glial cells in the central nervous system. High-grade glioma shows aggressive growth pattern, and conventional treatments, such as surgical removal and chemo-radiotherapy, archive limitation in the interference of this process. In this work, HOXA5, from the HOX family, was identified as a glioma cell proliferation-associated factor by investigating its feature in the TCGA and CGGA data set. High HOXA5 expression samples contain unfavorable clinical features of glioma, including IDH wild type, un-methylated MGMT status, non-codeletion 1p19q status, malignant molecular subtype. Survival analysis indicates that high HOXA5 expression samples are associated with worse clinical outcome. The CNVs and SNPs profile difference further confirmed the enrichment of glioma aggressive related biomarkers. In the meantime, the activation of DNA damage repair-related pathways and TP53-related pathways is also related to HOXA5 expression. In cell lines, U87MG and U251, by interfering HOXA5 expression significantly inhibit glioma progression and apoptosis, and cell cycle is arrested at the G2/M phase. Collectively, increased HOXA5 expression can promote glioma progression via affecting glioma cell proliferation.


INTRODUCTION
Glioma is a malignant tumor that derived from the central nervous system. Four grades were proposed for evaluating its malignancy according to its pathological features, including grades I, II, II, and IV. Grade IV glioma, also called as glioblastoma or glioblastoma multiforme (GBM), is the worst subtype of glioma with median survival time less than 14.6 months (1). High proliferative ability, disordered tumor angiogenesis, and massive tumor necrotic area were recognized as its feature in GBM (2). Recent study proposed different classification of GBM based on its transcriptomic signatures (3,4). Treatment like the STUPP protocol is able to prolong patient's survival time but cannot inhibit glioma progression. Therefore, identification of potential target to glioma may assist in improving its prognosis.
HOXA5 belongs to the gene family encoding transcription factors, which contained homeobox. This gene family expressed subsequently and affected tissue developmental and organogenesis (5). HOXA5 is located on chromosome 7 in human and encodes ANTP-class homeodomain protein consisting of 270 amino acids (5). HOXA5 widely participated in the development of organs such as the respiratory system, digestive system, lipid tissue, mammary gland, and so on (5,6). Therefore, HOXA5 acts a critical role in human development.
In this work, the expression profile of the HOXA5 in multiple tumor types was mapped. Its association with glioma progression and potential mechanisms were analyzed by using the GO/ KEGG enrichment analysis. We identified that HOXA5 as an unfavorable factor for glioma and its expression are highly connected with the activation of p53-related pathways. We also confirmed that HOXA5 affected tumor cells progression through regulating tumor cell proliferation.

Biological Function and Gene Set Enrichment Analyses
Correlation analysis of HOXA5 regarding different pathological features was performed in the TCGA and CGGA data sets with R language (https://www.r-project.org/). Differentially expressed genes between HOXA5 high and HOXA5 low groups with the adjusted p-value <0.05. and the absolute FC larger than 2.0 were considered to be statistically significant. Association between HOXA5 expression and gene sets from the Molecular Signatures Database (MSigDB) were analyzed using gene set enrichment analysis (GSEA). Gene ontology (GO), KEGG (Kyoto Encyclopedia of Genes and Genomes), and HALLMARK analysis was performed using the R package GSVA. Somatic mutations and somatic copy number alternations (CNAs) were downloaded from the TCGA database. Copy number alternations associated with HOXA5 expression were analyzed using GISTIC 2.0.

Survival Analysis
Patients were subdivided into high and low groups according to HOXA5 expression. The overall survival (OS), progression-free interval (PFI), and disease-specific survival (DSS) rates of patients in low and high groups were compared by the Kaplan-Meier method with log-rank test. ROC was performed to evaluate the prediction performance of HOXA5 expression in various aspects, including 3-year, 5-year OS, and subtype of GBM (classical, mesenchymal, neural, proneural).

Cell Culture
U87-MG cells and U251 cells, purchased from the Chinese Academy of Sciences, were cultured in DMEM medium with 10% Gibico FBS+1% penicillin-streptomycin, and the medium was changed every 2 to 3 days.

CCK8 Assay
The logarithmic growth phase transfected U251-MG and U87-MG GBM cells were obtained and digested for CCK8 assay. 1 × 10 3 glioma cells and 100 ml of medium were placed into 96-well plates. The absorbance at 450 nm was measured after hatched for 1 h under the condition of 37°C and 5 % CO 2 .

Colony-Forming Assay
U87-MG cells and U251 cells were digested and plated in six-well plates (300 cells per well) and cultured with 5% CO 2 at 37°C for 2 weeks. The colonies were then fixed with 4% methanol (1 ml per well) for 15 min and stained with crystal violet for 30 min at room temperature. After photograph, discoloration was performed with 10% acetic acid, and cells were measured absorbance at 550 nm.

Cell Cycle Assay
Transfected cells were digested, and cell suspension was obtained after centrifuging. Then, cells were washed two to three times with PBS, and we adjusted the number of cells to 1 × 10 6 cells/ml, 400 µl PBS was added, the cells were gently resuspended to separate as individual cells, 1.2 ml of pre-cooled 100% ethanol was added, and the cells were placed overnight at 4°C for fixation. PI was excited by a 488-nm argon ion laser and was received by a 630-nm pass filter. The percentage of each cell cycle was analyzed using the PI fluorescence histogram.

Cell Apoptosis Assay
Transfected cells were digested, and cell suspension was obtained after centrifuging. Then, cells were washed two to three times with PBS, 500-µl binding buffer was added, then the cells were gently resuspended to separate into individual cells, followed by staining with 5-µl Annexin V-APC and 5-µl PI solution for 10 min at room temperature and in dark place. The apoptotic cells were measured by the flow cytometer.

Statistical Analysis
Kaplan-Meier survival curves were generated and compared using the log-rank test. Wilcoxon rank test (nonnormally distributed variables), t test (normally distributed variables), and one-way analysis of variance were used to analyze the expression difference of HOXA5 in different clinical factors, including WHO grades, GBM subtypes, and treatment outcome. The Pearson correlation was applied to evaluate the linear relationship between gene expression levels. Statistical analyses of the colony-forming assay and the CCK8 assay were carried out by GraphPad Prism (version 8.0). The Kolmogorov-Smirnov test was used to assess the normal distribution of data. All tests were two-sided, and P values <0.05 were considered to be statistically significant.

HOXA5 Expression Is Elevated in Aggressive Gliomas
The mRNA expression levels of HOXA5 in different WHO grade gliomas were evaluated using expression data from publicly available databases: TCGA, n = 672; CGGA, n = 1013. HOXA5 was observed to be significantly up-regulated in GBM (WHO grade IV) compared with low-grade glioma (LGG) samples (WHO grade III and WHO grade II) in the TCGA and CGGA cohorts (P <0.001, respectively; Figure 1A). The expression of HOXA5 was also higher in WHO grade III than WHO grade II cases in the TCGA and CGGA cohorts (P <0.001, respectively; Figure 1A). The HOXA5 levels in common cancer types other than gliomas were further evaluated (P <0.001, respectively; Figure 1B). We also evaluated the expression levels of HOXA5 in different age groups in pan-glioma analysis in TCGA and CGGA data sets, and in GBM patients in TCGA microarray data set, where patients older than 45 years have higher expression of HOXA5 (P <0.001, respectively; Figure 1C). HOXA5 was upregulated in the IDH mutant gliomas in TCGA data set, in the IDH mutant GBM in TCGA microarray data set, and in both the IDH mutant gliomas and GBM alone in CGGA data set (P <0.001, respectively; Figure 1D). Furthermore, HOXA5 was upregulated in the 1p19q codeletion in pan-glioma analysis in both TCGA and CGGA data sets (P <0.001, respectively; Figure 1E). Additionally, HOXA5 was upregulated in the unmethylated gliomas in TCGA data set (P <0.001; Figure 1F).

Inter-Tumor and Intra-Tumor Heterogeneous Characteristics of HOXA5 in Gliomas
Human gliomas have been molecularly categorized into distinct sub-classes: classical (CL), mesenchymal (MES), proneural (PN), and neural (NE). CL and MES subtypes showed more aggressive growth pattern than PN or NE subtypes (21). We subsequently investigated the inter-tumor heterogeneity of HOXA5 among different molecular subtypes based on the VERHAAK 2010 classification scheme (22). HOXA5 was upregulated in CL and ME subtypes by comparing with NE and PN subtypes in panglioma analysis in TCGA data set, where CL subtype had the highest expression of HOXA5 (P <0.001; Figure 2A). Receiver operating characteristic (ROC) curve further indicated that HOXA5 might serve as a predictor for CL and MES subtypes in pan-gliomas analysis (AUC value = 0.898, P <0.001; Figure 2B). Based on the Ivy Glioblastoma Atlas Project data, HOXA5 was found to be abundant in peri-necrotic zones, pseudopalisading cells around necrosis and cellular tumor compared with other pathological areas (P <0.001; Figure 2C). Furthermore, the different expression level of HOXA5 in glioma in regard to histology was shown in Figure 2D. HOXA5 was also highly expressed in primary glioma (P <0.001, respectively; Figure 2E), whereas the analysis for first-course treatment outcome showed that progressive disease was correlated with higher expression of HOXA5 (P <0.001, respectively; Figure 2F).

HOXA5 Expression Is Associated With Poor Survival in Glioma Patients
We next assessed the prognostic value of HOXA5 expression in human gliomas using ROC curve analysis and Kaplan-Meier analysis. The ROC curve analysis showed the prognostic value of HOXA5 in survival in TCGA and CGGA cohorts (3-year AUC value = 0.848, 5-year AUC value = 0.813; 3-year AUC value = 0.774, 5-year AUC value = 0.778, respectively; Figure 2G). Kaplan-Meier survival curves were generated based on median values of HOXA5 expression in gliomas. In both TCGA and CGGA data sets, HOXA5 high patients exhibited significantly shorter overall survival (OS), disease-specific survival (DSS), progression free interval (PFI) than HOXA5 low patients in pan-glioma analysis, (P <0.001, respectively; Figure 3A, Figure  S1A). In addition, HOXA5 high patients exhibited significantly shorter OS, DSS, PFI compared with HOXA5 low patients in LGG alone in TCGA data set (P <0.001, respectively; Figure 3B). Poor OS, DSS, PFI was also associated with high expression of HOXA5 in GBM alone in TCGA data set (P <0.001, respectively; Figure 3C). The DSS, OS, and PFI survival probabilities of GBM patients were further testified in the TCGA microarray data set (P <0.001, respectively; Figure 4A), in which HOXA5 was correlated with poor survival. Further, HOXA5 predicted worse OS in GSE108474 data set ( Figure 4B). In TCGA microarray data set, HOXA5 had the highest expression in GBM patients without IDH mutation (P <0.001, respectively; Figure 4C), in GBM patients without radiotherapy (P <0.001, respectively; Figure 4D), and in GBM patients without chemotherapy (P <0.001, respectively; Figure 4E). Similar results were obtained in CGGA data set (P <0.001, respectively; Figures S1F, S1G, S1H). In addition, in pan-glioma analysis in TCGA microarray data set, HOXA5 also had the highest expression in patients without IDH mutation (P <0.001, respectively; Figure 4F), in patients without radiotherapy (P <0.001, respectively; Figure 4G). The results were further confirmed in CGGA data set (P <0.001, respectively; Figures  S1B-D). Patients without 1p19q codeletion also had the highest expression of HOXA5 in pan-glioma analysis and LGG alone in TCGA microarray data set (P <0.001, respectively; Figures 4H, I).
Similar results were obtained in pan-glioma analysis and LGG alone in CGGA data set (P <0.001, respectively; Figures S1E, I).

HOXA5 Expression Levels Are Associated With Distinct Genomic Alterations
As for the 15 methylation probes designed for HOXA5 from Infinium Human Methylation450 BeadChip, the mean value of methylation probes exhibited negative association with expression of HOXA5 (Pearson test, R = −0.61, P < 2.2 × 10-16). As the IDH mutation exerted great influence on the methylation of the whole genome, we separately analyzed the relationship between HOXA5 and methylation status for IDH-wild and IDH-mutant gliomas. As shown in Figure 5A, in IDH-wild gliomas, methylation was much lower than that of IDH-mutant tumors, as expected (Student's t test, P = 1.8 × 10 −6 ).
To determine whether HOXA5 expression levels were associated with specific genomic characteristics in gliomas, we performed copy number variation (CNV) and somatic mutation  analysis using the TCGA data set. A distinct overall CNV profile emerged from the comparison of the HOXA5 low (n = 158) versus the HOXA5 high (n = 158) cluster ( Figures 5B, C). Amplification of chr7 and deletion of chr10, which are both common genomic events in GBM, frequently occurred in the HOXA5 high cluster ( Figure 5B). Non-deletion of 1p and 19q, a genomic hallmark of oligodendroglioma, also more frequently appeared to be associated with the HOXA5 high cluster ( Figure 5C) Analysis of somatic mutation profiles based on HOXA5 expression levels revealed a high frequency of mutations in EGFR (30%), TP53 (28%), TTN (25%), and PTEN (23%) in the HOXA5 high group (n = 158), whereas IDH1 (89%), TP53 (35%), CIC (32%), and ATRX (24%) were more frequently mutated in the HOXA5 low group (n = 158; Figure 5D).

Potential Mechanism of HOXA5 in Regulating the Progression of Gliomas
To elucidate whether HOXA5 could play a role in promoting gliomas occurrence, GSVA analysis was performed. HOXA5 was found to be associated with regulation of DNA damage response signal transduction by p53 class mediator, signal transduction by p53 class mediator, nucleotide excision repair DNA gap filling, DNA synthesis involved in DNA repair, DNA damage response detection of DNA damage, signal transduction in response to DNA damage, mismatch repair, and base excision repair in TCGA and CGGA data sets ( Figure 6A, Figure S2A). The correlation analysis between HOXA5 and these GO pathways is shown in Figure 6B. As the threshold was set as logFC > 2 and adjust P <0.01, a total number of 3,446 differentially expressed genes (DEGs) were detected between high expression of HOXA5 sample and low expression of HOXA5 sample ( Figure 6C). We paid special attention to two pathways, DNA damage response signal transduction by p53 class mediator and signal transduction by p53 class mediator, in which the GO analysis revealed changes in gene sets related to these two pathways in patients with a higher expression of HOXA5 (Figures 6D, E). We further investigated whether HOXA5 might have a role in DNA damage response and p53 signal transduction in gliomas using GSEA analysis in TCGA and CGGA data sets ( Figure 6F, Figure  S2B). In KEGG pathway analysis, HOXA5 was found to be associated with mismatch repair, DNA replication, p53 signaling pathway, and base excision repair in TCGA and CGGA data sets (Figures 7A, B). The correlation analysis between HOXA5 and KEGG pathways in TCGA and CGGA data sets is shown in Figure S2C and Figure S2D, respectively. We also investigated the role that HOXA5 might play in mismatch repair and p53 signal transduction in gliomas using GSEA analysis in TCGA and CGGA data sets (Figures 7C, D). In HALLMARK pathway analysis, HOXA5 was found to be associated with g2m checkpoint and p53 pathway in TCGA and CGGA data sets ( Figures 7E, F), in which the correlation analysis is shown in Figure S2E and Figure S2F, respectively. These results all suggested that HOXA5 was associated with oncogenic processes.

HOXA5 Affect Glioblastoma Cell Proliferation, Viability, and Apoptosis
To further investigate the role that HOXA5 plays in the proliferation of GBM, colony-forming assay, CCK8 assay, and cell cycle analyses were performed. Western blot results verified the silence of HOXA5 by siRNA (Figures 8A, B). The cell colony forming assay (Figures 8C, D) revealed the remarkable suppression of cell clonality of GBM after silencing HOXA5 in U87 cell line and in U251 cell line ( Figure 8E). The CCK8 assay revealed that the cells proliferation ability is inhibited by silencing HOXA5 ( Figure 8F). Cell cycle analysis suggested that cells were blocked at the G2/M phase after silencing HOXA5 ( Figures 8G, H).
The potential role of HOXA5 in the apoptosis of GBM was also explored. The analysis of flow cytometry described that down-regulation of HOXA5 remarkably enhanced apoptosis in U87 cell line and U251 cell line ( Figures 9A, B).

DISCUSSION
Glioma is a central nervous system-derived tumor with highly aggressive growth pattern, treatment resistance, and poor prognosis. Nowadays, strategies show their limitations on improving patient's survival outcome, including surgery, chemotherapy, and radiotherapy (23). To reveal the mechanism concealed behind its aggressive growth pattern, we identify HOXA5 as prognostic-related genes through regulating cell proliferation. HOXA5, locate on chromosome 7, is involved in tissue development and organogenesis. Multiple studies reported that high HOXA5 expression was a promotor in tumor progression, including esophageal cancer (24), breast cancer (12), gastric cancer (25), and renal cancer (26). In this study, high HOXA5 expression is associated with malignant clinical features like high-grade glioma or IDH wildtype glioma, worse clinical treatment outcome, and strong cells proliferative ability. We also noticed genes like EGFR, PDGFRA were amplified in HOXA5 high samples; and genes like CDKN2A/CDKN2B, PARK7 were down-regulated. Those genes were also identified as glioma promotor-related factors (27)(28)(29)(30)(31)(32)(33). Besides, tumor cells cycle was arrested at G2/M phase when inhibited HOXA5 expression implying its role in mediating cell cycle. Along with results from other studies that HOXA5 modulated cell cycle in Jurkat cells and hematopoietic stem cells (34), cervical cancer, mesenchymal stem cells (35), and leukemia cells (36). Further, HOXA5 was found to decrease the apoptosis rates of GBM cells. Together, those evidences supported that HOXA5 promoted glioma progression through affecting cell cycle, cell proliferation, and cell apoptosis.
The product of HOXA5 is a DNA-binding transcription factor and able to regulate multiple gene expression, including TP53. Previous study reported that HOXA5 can bind to the promoter of TP53 to activate its transcription and affect corresponded pathways (16,37). HOXA5 affected tumor progression through influencing TP53 homeostasis in breast cancer (10,38), TP53-dependent apoptosis in liposarcomas (39), and TP53-mediated cell proliferation in cervical cancer   (16). Therefore, HOXA5 can affect TP53-mediated pathway by regulating TP53 expression. Next, we investigated biofunction difference between HOXA5 high and HOXA5 low samples. Results suggested that DNA repair-related pathways were activated in HOXA5 high samples. Because DNA damage repair is also associated tumor resistance to treatment, such as chemo-or radio-therapy, HOXA5 might also contribute to tumor resistance to treatments (40,41). In this work, survival analysis based on radio-and chemo-therapy suggested that high HOXA5 samples showed worse clinical outcome than low HOXA5 samples in this study. Besides, previous study confirmed that HOXA5 can affect glioblastoma cells sensitivity to radiation therapy (9). Collectively, HOXA5 may target TP53-mediated DNA damage repair to modulate tumor resistance to treatments (42)(43)(44). Taken together, HOXA5 may affect glioma response to chemoor radio-therapy by targeting TP53.
In conclusion, this work proved that high HOXA5 expression positively correlated with glioma malignancy's clinical features and was associated with an unfavorable survival outcome. Moreover, HOXA5 may affect glioma progression and apoptosis by modulating TP53 expression and corresponding pathways. Critically, HOXA5 can arrest cell cycle at G2/M phase to inhibit tumor cell proliferation.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
FD and PC analyzed the data and performed in vitro environments. PB assisted in manuscript writing and revising. WP and QC. designed and supervised the study, and offered funding. All authors contributed to the article and approved the submitted version.