Next-generation sequencing identifies HOXA6 as a novel oncogenic gene in low grade glioma

Analysis the expression, clinical information and prognosis of HOXAs

We utilized TCGA database (https://www.cancer.gov/) and CGGA (http://www.cgga.org.cn/) [8–10] to examine the expression, clinical information and prognosis of HOXAs in LGG. HPA (https://www.proteinatlas.org/) database utilized to examine the protein of HOXAs in LGG [11].

Function analysis for HOXAs in LGG

In the present research, we utilized the LinkedOmics database (http://www.linkedomics.org/login.php) obtained the co-expression genes of HOXAs in LGG. We were using clusterProfiler package to examine the function of HOXAS in LGG [12, 13]. GeneMANIA database (http://genemania.org/) used to constructed gene-gene interaction network of HOXAs [14].

Cox regression analysis and Kaplan-Meier survival analysis

We utilized cox regression analysis to examine the correlation between HOXAs expression and overall survival and disease-specific survival of patients using the TCGA databases. The Kaplan-Meier method was used to assess the difference between high and low risk groups based on the best separation of HOXAs expression, employing R packages of survminer and survival.

Immune infiltration analysis

TIMER (https://cistrome.shinyapps.io/timer/) [15], an interactive web portal, could perform comprehensive analysis on the infiltration levels of different immune cells. In this study, we were using the TIMER database to explore the correlation between HOXAs and diverse immune cell infiltration in LGG. We also were using GSVA R package to quantify the LGG immune infiltration of 24 tumor-infiltrating immune cells in tumor samples via ssGSEA [16]. The TISIDB (http://cis.hku.hk/TISIDB/) database utilized analysis the expression of HOXAs in different immune subtype of LGG [17].

Analysis the correlation between the HOXAs expression and drug sensitivity

The Genomics of Drug Sensitivity in Cancer (GDSC) database (http://www.cancerRxgene.org) is the largest public resource for information on drug sensitivity in cancer cells and molecular markers of drug response. Data are freely available without restriction. GDSC currently contains drug sensitivity data for almost 75 000 experiments, describing response to 138 anticancer drugs across almost 700 cancer cell lines. In this study, we utilized the Genomics of Drug Sensitivity in Cancer (GDSC) (http://www.cancerRxgene.org) and the Cancer Therapeutics Response Portal (http://www.broadinstitute.org/ctrp) databases to analysis the correlation between HOXAs expression and drug sensitivity [18, 19]. Spearman’s rank correlation was employed to measure the correlation between HOXAs expression and diverse drug sensitivity.

Gene mutation analysis

The cBioPortal database (http://www.cbioportal.org/) and GSCA (http://bioinfo.life.hust.edu.cn/web/GSCALite/) used to analysis the gene mutation, DNA methylation of HOXAs in LGG [20, 21].

LncRNA/miRNA/mRNA network construction

To explore the potential upstream regulatory of HOXAs in LGG, we used the LnCeVar database (http://www.bio-bigdata.net/LnCeVar/index.jsp)constructed a competing endogenous RNA (ceRNA) network [22].

Cell culture conditions

GBM cells lines (including A172, SF295, T98G and U251 cells) were purchased from cell bank of Kunming Institute of Zoology, and cultured in DMEM medium (Corning) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin at 37° C in atmosphere containing 95% air and 5% CO2.

Quantitative real-time PCR

The qRT-PCR assay was performed as documented [23]. The primer sequences are list follows: HOXA6-F: TCCCGGACAAGACGTACAC, HOXA6-R: CGCCACTGAGGTCCTTATCA. B-actin-F: CTTCGCGGGCGACGAT, β-actin-R: CCATAGGAATCCTTCTGACC. The expression quantification was obtained with the 2^−ΔΔCt method.

Cell proliferation assay

Cell proliferation was performed as previously documented [24]. Briefly, for cell proliferation assay, indicated cells were plated into 12-well plates at a density of 1.5×10⁴, the cell numbers were subsequently counted each day using an automatic cell analyzer Countstar (Shanghai Ruiyu Biotech Co., China, IC1000). (ns, p ≥ 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001).

Cell migration assay

Cell migration assays was performed as previously documented [24]. Briefly, to produce a wound, the monolayer cells in 6-well plate were scraped in a straight line with pipette tips. Plate was then washed with PBS to remove detached cells. Photographs of the scratch were taken at indicated time points using Nikon inverted microscope (Ti-S).

Statistical analysis

Correlation analysis was performed using Pearson correlation test. Kaplan-Meier survival curves were plotted to exhibit the overall survival for LGG patients. Univariate and multivariate Cox regression analyses were used to examine the independent prognostic significance of each variable enrolled in this finding. The significance of the data between two experimental groups was determined by Student’s t-test and multiple group comparisons were analyzed by one-way ANOVA. P < 0.05 (*), P < 0.01 (**) and P < 0.001 (***), were considered significant.

HOXAs were highly expressed in LGG

We used TCGA and GTEx databases explored the expression of 11 HOXAs in GBM and LGG, the results showed that 11 HOXAs were significantly higher in lower-grade glioma (LGG) and glioblastoma multiforme (GBM) than those in normal samples (Figure 1A, 1B). Furthermore, we uncover that the protein of HOXA4, HOXA6, HOXA9 and HOXA10 were highly expressed in LGG tissue compared with control group based on the Human Protein Atlas datasets (Figure 1C). Collectively, these results suggested that HOXAs were overexpressed in LGG tissues.

Correlation between HOXAs expression level and clinicopathological characteristics of patients with LGG

We further explored the correlation between HOXAs and clinical characteristics of LGG. The results showed that up-regulated 11 HOXAs were significantly associated with a higher tumor stage, IDH mutation status, 1p/19q chromosome co-deletion, histological type, primary therapy outcome, and age based on TCGA dataset (Figures 2, 3 and Supplementary Figure 1, Supplementary Table 1). This result was verified by CGGA database (Supplementary Figures 2, 3). The Logistic regression analysis of results also confirmed that HOXAs expression was remarkably associated with poor clinicopathological characteristics (Figure 4). Given that HOXAs were up-regulated in LGG and its higher expression associated with poor clinical characteristics. Therefore, we further explored the correlation among distinct HOXAs and the survival rate of LGG patients. The results indicated that higher expression of 11 HOXA was correlated with shorter overall survival (OS), disease-specific survival (DSS) and progression-free survival (PFS) in LGG patients (Figures 5–7). These results were verified by CGGA datasets (Supplementary Figure 4). Given that HOXAs differentially expression and associated with poor clinic-pathologic features and prognosis in LGG. The univariate cox proportional hazards regression analyses results confirmed that HOXA members (HOXA1, HOXA2, HOXA3, HOXA4, HOXA5, HOXA6, HOXA7, HOXA9, HOXA10 and HOXA11) expression and six clinical features (WHO grade, 1p/19q co-deletion, primary therapy outcome, IDH status, histological type and age) were correlated with adverse clinical outcomes of LGG patients. Additionally, the multivariate analyses revealed that HOXA1, HOXA6 expression and tumor grade, age, primary therapy outcome and age were independent factors affecting the prognosis of LGG patients (Table 1).

Diagnostic values of HOXAs in LGG

Given that HOXAs were significantly up-regulated in LGG and its higher expression was associated with adverse clinical features and adverse clinical outcomes. Therefore, we decided explored the diagnosis values of HOXAs in LGG. ROC curve analysis of HOXAs showed that HOXAs had a high accuracy (AUC > 0.80) in predicting LGG (Figure 8). These results confirmed that HOXAs has the potential to act as a diagnosis biomarker for the diagnosis of LGG with high sensitivity and specificity.

Construction and validation of HOXAs based nomogram

In order to provide clinicians with predicted prognosis of LGG patients quantitatively, we constructed the nomogram combining HOXA1, HOXA6 expression and independent clinical risk factors (tumor grade, age, primary therapy outcome and age). Higher total points on the nomogram for overall survival (OS), disease-specific survival (DSS) and progression-free survival (PFS), respectively, indicated a worse prognosis (Figure 9). In summary, these results indicate that the nomogram can well predict short- or long-term survival of LGG patients.

Gene mutation analysis

We utilized cBioPortal analysis the genetic alteration of differentially expressed HOXAs family. Results confirmed that HOXA1, HOXA2, HOXA3, HOXA4, HOXA5, HOXA6, HOXA7, HOXA9, HOXA10, HOXA11 and HOXA13 were all altered, with 1.2, 0.8, 1, 0.6, 0.8, 0.8, 1, 1, 0.6, 1 and 0.6alterations in the LGG samples, respectively (Figure 10A). We also summarize the incidence of CNV and somatic mutations of 11 HOXAs in LGG. The results confirmed that missense mutation was the most common variant classification, SNPs were the most common variant type and C > T ranked as the top SNV class. The results also showed HOXA1 as the gene with the highest mutation frequency, followed by HOXA11 and HOXA6, among the 11 HOXAs (Figure 10B, 10C). CNV analysis results confirmed that copy number variations of HOXAs was significantly positive correlated with its expression in LGG (Figure 10D). Furthermore, HOXAs CNV affect the overall survival (OS), disease-specific survival (DSS) and progression-free survival (PFS) in LGG patients (Figure 10E). Finally, we uncover that the DNA methylation level of HOXA13 negatively associated with its expression (Figure 10F). More importantly, lower DNA methylation level of HOXA1, HOXA3, HOXA9 and HOXA10 were correlated with poor DSS, OS and PFS (Figure 10G).

Functional analysis of HOXAs in LGG

To explore the potential functions of HOXAs in LGG progression, we were using the GeneMANIA database constructed the gene-gene interaction network for HOXAs and the altered neighboring genes (Figure 11A). Furthermore, we perform the GO and KEGG enrichment to explore the functions of HOXAs and their neighboring genes in LGG. The enrichment results confirmed that the biological processes for these genes were mostly neutrophil activation, neutrophil mediated immunity, neutrophil activation involved in immune response, T cell activation, positive regulation of cytokine production, regulation of immune effector process, regulation of lymphocyte activation, positive regulation of cell activation, lymphocyte differentiation, mononuclear cell proliferation, lymphocyte proliferation, leukocyte cell−cell adhesion, regulation of T cell activation, regulation of leukocyte cell−cell adhesion, regulation of mononuclear cell proliferation, regulation of lymphocyte proliferation, response to interferon−gamma, T cell differentiation and T cell proliferation (Figure 11B). The molecular functions for these genes were mostly secretory granule membrane, external side of plasma membrane, lysosomal membrane, membrane microdomain, phagocytic vesicle, specific granule membrane, MHC protein complex and MHC class II protein complex (Figure 11C). KEGG enrichment results confirmed that these genes mainly involved in the Cytokine−cytokine receptor interaction, Osteoclast differentiation, Cell adhesion molecules, Kaposi sarcoma−associated herpes, virus infection, Hematopoietic cell lineage, Chemokine signaling pathway, Th17 cell differentiation, NOD−like receptor signaling pathway, JAK−STAT signaling pathway, Natural killer cell mediated cytotoxicity, Th1 and Th2 cell differentiation, Inflammatory bowel disease, Leukocyte transendothelial migration, B cell receptor signaling pathway, Toll−like receptor signaling pathway, NF−kappa B signaling pathway, Viral protein interaction with cytokine and cytokine receptor and TNF signaling pathway (Figure 11D). Collectively, these data suggest that HOXAs plays an essential regulatory role in the inflammation response and immune regulation.

Correlation between HOXAs and immune subtypes of LGG

For determine the expression of HOXAs in diverse immune subtypes of LGG, we utilized the TISIDB database analysis and found that HOXAs mainly highly expressed in C3 and C4 subtype (Supplementary Figure 5).

Immune infiltration analysis of the HOXAs in LGG

Given that HOXAs was significantly correlated with immune regulation related signaling pathway. Therefore, we decide to examine the correlation between HOXAs expression and immune cell infiltration. In this finding, we used TIMER database analysis and found that somatic copy number alterations of HOXAs were significantly correlated with diverse immune cell infiltration levels in LGG (Supplementary Figures 6, 7). Furthermore, we also utilized the TIMER database analysis the correlation between HOXAs and diverse immune cell. The results confirm that HOXAs were positively correlated with the immune infiltration of B cells, CD4+ T cells, CD8+ T cells, Neutrophils, Macrophages and Dendritic cells in LGG (Figure 12 and Supplementary Figure 8).

Finally, ssGSEA with Spearman’s rank correlation was employed to measure the correlation between HOXAs expression and infiltration levels of 24 immune cell types in glioma. The results showed that HOXA1 expression was correlated with the immune infiltration of macrophages and eosinophils, negatively correlated with the immune infiltration of NK CD56bright cells and pDC. HOXA2 positively correlated with the immune infiltration of Macrophages, Eosinophils and T cells, negatively correlated with the immune infiltration of Tcm and pDC. HOXA3 positively correlated with the immune infiltration of NK cells, Th2 cells and T cells, negatively correlated with the immune infiltration of NK CD56bright cells and pDC. HOXA4 positively correlated with the immune infiltration of macrophages, eosinophils and neutrophils, negatively correlated with the immune infiltration of pDC and NK CD56bright cells. HOXA5 positively correlated with the immune infiltration of Macrophages, Eosinophils, NK cells and Neutrophils, negatively correlated with the immune infiltration of NK CD56bright cells and pDC. HOXA6 positively correlated with the immune infiltration of Eosinophils, Th2 cells and Macrophages, negatively correlated with the immune infiltration of CD8 T cells and pDC. HOXA7 positively correlated with the immune infiltration of Th2 cells, neutrophils and NK cells, negatively correlated with the immune infiltration of NK CD56bright cells and pDC. HOXA9 positively correlated with the immune infiltration of Neutrophils and Eosinophils, negatively correlated with the immune infiltration of NK CD56bright cells and TReg. Furthermore, HOXA10 expression associated with the immune infiltration of several subsets of myeloid cells, including macrophages, Eosinophils and Neutrophils, HOXA10 expression negatively correlated with the immune infiltration of TReg and NK CD56bright cells. There are positive correlations between HOXA11 expression and the immune infiltration of Macrophages, T helper cells and Neutrophils, negatively correlated with the immune infiltration of TRegNK CD56bright cells. There is also a strong positive correlation between HOXA13 expression and immune infiltration of Tgd, T helper cells and Macrophages, negatively correlated with the immune infiltration of TReg and NK CD56bright cells (Figure 13 and Supplementary Figure 9). Taken together, these results confirmed that HOXAs plays a pivotal role in the immune regulation.

Correlation between HOXAs and immune modulator in LGG

Considering immune modulator plays a crucial role in the immune response and development and progression of LGG. We explored the correlation between HOXAs and the immune checkpoint-related genes by using Pearson correlation analysis in LGG, including CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG12, TIGIT, SIGLEC15. The result showed that HOXAs expression was significantly positively correlated with CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG12 and SIGLEC15 (Figure 14). These results demonstrate that HOXAs expression was significantly correlated with the expression of immune checkpoint-related genes in LGG.

Correlation between HOXAs expression and diverse drug sensitivity

We downloaded the IC50 values of anti-cancer drugs and gene expression profiles in the relative cell lines from the GDSC and CTRP database. To explore the influence of HOXAs expression on the sensitivity of anti-cancer drugs, we divided the tumor cells into high- and low- HOXAs groups and compared their IC50 values. For example, we found that the IC50 values of several anti-cancer drugs decreased in the high- HOXA3 group, including abiraterone, Trametinib, VAF-347, BRD-K09344309, BRD-K34099515 and erlotinib (Supplementary Figure 10), indicating that patients exhibiting high HOXA3 expression levels are relatively sensitive to these anti-cancer drugs. We also uncover that HOXAs were positively correlated with the diverse drug sensitivity. For instance, HOXA3 was positively correlated with the diverse drug sensitivity of tacedinaline, BRD-A94377914, LRRK2-IN-1, PX-12, Apicidin, PRIMA-1, Necrosulfonamide, tipifarnib-P2, vorinostat, avrainvillamide, serdemetan, cerulenin, CCT036477 and PL-DI (Supplementary Figure 10). Taken together, these results confirmed that HOXAs was positively or negatively correlated with the diverse drug sensitivity in the cancer therapeutic response portal database.

Construction of a network of HOXAs related ceRNA

Emerging reports have revealed that lncRNA /miRNA regulatory axis play a vital role in the gene expression [25]. To uncover that the upstream regulatory mechanism of HOXAs in LGG, we utilized the LnCeVar database identification of HOXAs related ceRNA events in LGG of patients from TCGA. The results confirmed that ZHF337-AS1 or SNHG12/miR-210-3p / regulatory axis may regulate the expression of HOXA1. LINC00174/miR-26a-5p regulatory axis may regulate the expression of HOXA5. HCG18 or TRG- AS1/miR-196a-5p regulatory axis may regulate the expression of HOXA7. LNC00665/miR-320a regulatory axis may regulate he expression of HOXA10. KDM4A-AS1/ miR-27a-5p regulatory axis may regulate he expression of HOXA13 (Supplementary Figure 11). These results confirmed that HOXAs related ceRNA network regulatory axis may play a vital role in the progression of LGG.

Knock down of HOXA6 inhibits the cell proliferation and migration of LGG cells

Given that HOXA1 expression was independent factors affecting the prognosis of LGG patients and previous no study reported the function of HOXA6 in LGG. Therefore, we decide explored function of HOXA6 in LGG. We first examined the expression of HOXA6 in LGG cells lines by qRT-PCR assay. The results confirmed that HOXA6 was up-regulated in LGG cells lines, especially in SF295 and A172 cells (Figure 15A). Next, we utilized siRNA knock down of HOXA6 in SF295 and A172 cells, the knock down efficiency was verify by qRT-PCR assay (Figure 15B, 15C). The growth curve assay showed that depletion of HOXA6 significantly inhibits LGG cells proliferation ability (Figure 15D, 15E). Furthermore, the transwell assay confirmed that HOXA6 knock down significantly inhibits LGG cells migration ability (Figure 15F). Collectively, these results demonstrate that HOXA6 was highly expressed in LGG cells lines and significantly affected their proliferation and migration abilities.