From the The Cancer Genome Atlas (TCGA) (https://xenabrowser.net/), Chinese Glioma Genome Atlas (CGGA) (http://www.cgga.org.cn/), and Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/) data sets, we collected HOXA5 data from LGG and GBM samples. RNA-seq data for specific tumor anatomic structure in GBM was from Ivy Glioblastoma Atlas Project (http://glioblastoma.alleninstitute.org/).

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.

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).

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.

The siRNA of HOXA5, siRNA-112 (sense 5′-3′: GGACUACCAGUUGCAUAAUTT; antisense 5′-3′: AUUAUGCAACUGGUAGUCCTT), siRNA-610 (sense 5′-3′: GCACAUAAGUCAUGACAACTT; antisense 5′-3′: GUUGUCAUGACUUAUGUGCTT) and siRNA-726 (sense 5′-3′: GCAGAAGGAGGAUUGAAAUTT; antisense 5′-3′: AUUUCAAUCCUCCUUCUGCTT) were purchased from HonorGene (Changsha, China). 95 µl serum-free DMEM medium, 5 µl HOXA5 siRNA, and 5 µl Lip2000 were added into the centrifuge tubes in turns. HOXA5 siRNA-112 was discarded because of its low efficiency according to the Western blot assay.

HOXA5 primary antibody (Abcam, Cat# ab140636, RRID: AB_2877721) was diluted to 1:1000. β-actin (Proteintech Cat# 66009-1-Ig, RRID: AB_2687938) was diluted to 1:5000 and utilized as the loading control and internal standard. Secondary antibody (Proteintech Cat# SA00001-2, RRID: AB_2722564; Proteintech Cat# SA00001-1, RRID: AB_2722565) was diluted to 1:5000. The total protein concentration was determined using the BCA Protein Assay Kit (Solarbio, China), according to the manufacturer’s instructions. The blots were subjected to three 5-min washes with TBST prior to 1-h incubation with the secondary antibody, repeated washing, and signal development with Western Lightning Plus-ECL.

The logarithmic growth phase transfected U251-MG and U87-MG GBM cells were obtained and digested for CCK8 assay. 1 × 10³ glioma cells and 100 μl 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₂.

U87-MG cells and U251 cells were digested and plated in six-well plates (300 cells per well) and cultured with 5% CO₂ 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.

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⁶ 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.

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.

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.

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 ).

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 pan-glioma 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 ).

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 ).

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⁻⁶).

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 ). In HOXA5^high samples, frequently amplified genomic regions included oncogenic driver genes, such as EGFR (7p11.2), IK3C2B (1q32.1), PDGFRA (4q12), and CDK4 (12q14.1), whereas deleted regions contained tumor suppressor genes, including CDKN2A/CDKN2B (9p21.3), PARK7 (1p36.23), and PTEN (10q23.3). In HOXA5^low samples, significant amplifications showed peaks in 8q24.21, 12p32.32, and 19p13.3, whereas the frequently deleted genomic regions were 2q37.3, 4q32.3, 9p21.3, and 19q13.43.

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 ).

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.

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 ).