HOXA9 is dramatically upregulated in hepatocellular carcinoma

To explore the expression of HoxA9 in different organ and tissues, we analyzed mouse single-cell transcriptome dataset (GSE109774) derived from Tabula Muris (Supplementary Figure 1A). The mRNA level of HoxA9 gene was highly expressed in bladder, kidney, muscle and marrow, and rarely expressed in heart, aorta, liver, thymus, tongue, trachea, spleen, mammary gland, and lung (Supplementary Figure 1B, C).

HoxA9 gene has been found to be aberrant expressed in several solid tumors, including colon cancer, nasopharyngeal carcinoma and breast cancer (Bhatlekar et al. 2014). The analysis of expression profiles in 33 tumors of TCGA (The Cancer Genome Atlas) database showed that HoxA9 mRNA was significantly upregulated in CHOL (cholangiocarcinoma), COAD (colon cancer), ESCA (Esophageal carcinoma),GBM (Glioblastoma multiforme), HNSC (head and neck squamous cell carcinoma), KIRP (Kidney renal papillary cell carcinoma), LIHC (hepatocellular carcinoma), LUSC (lung squamous cell carcinoma), PRAD (Prostate adenocarcinoma), READ (rectal cancer) and STAD (Stomach adenocarcinoma) and was significantly downregulated in BRCA (Breast invasive carcinoma), KICH (Kidney Chromophobe), KIRC (renal clear cell carcinoma), THCA (Thyroid carcinoma) (Figure 1A). These results are consistent with the previous studies that HoxA9 gene plays opposite roles in the development of distinct cancers (Bhatlekar et al. 2014). Surprisingly, the mRNA level of HoxA9 gene was upregulated in HCC tumors compared to adjacent non-tumor tissues, which was never reported before. Similar results were obtained by analyzing different HCC cohorts derived from TCGA+GTEx and Gene Expression Omnibus (GEO) datasets (Figure 1B), implying a potential and important role of HoxA9 in HCC development.

Fig. 1

HoxA9 RNA expression was significantly upregulated in HCC and upregulation of HOXA9 in HCC predicts poor prognosis. A TCGA datasets were used to analyze the RNA expression of HoxA9 in tumors from TIMER2.0. The RNA expression of HoxA9 was low in normal liver tissues and the expression of HoxA9 was significantly upregulated in hepatocellular carcinoma (red virtual box); B The RNA expression of HoxA9 was significantly upregulated in liver cancer tissues according to the TCGA + GTEx and multi-GEO datasets. The GEO liver cancer datasets include: GSE54236, GSE36376, GSE25097. C Kaplan–Meier curves of overall survival (OS) in HOXA9high and HOXA9low patients with LIHC using TCGA data from the HPA database; D the expression difference of HOXA9 at the protein level in liver tissues showed by immunohistochemical results from HPA public database (antibody number HPA-061982); Scale bars, 200 μm; E–K Box plots evaluating HOXA9 expression among different groups of patients based on clinical parameters in LIHC using the UALCAN database. Analysis is shown for sample types (E), tumor grades (F), cancer stage (G), sex (H), metastasis (I), TP53 mutation status (J) and histological subtypes (K). N0: no regional lymph node metastasis; N1: metastases in 1 to 3 axillary lymph nodes. All the datasets were collected from the Human Protein Atlas website (https://www.proteinatlas.org/) and the UALCAN database (http://ualcan.path.uab.edu/analysis-prot.html). ***P < 0.001, **P < 0.01, *P < 0.05

Upregulation of HOXA9 in HCC predicts poor prognosis

To investigate whether HOXA9 expression correlates with patients’ survival, we analyzed it in HOXA9high and HOXA9low patients with liver cancer from the Human Protein Atlas (HPA) database. As expected, patients with HOXA9high liver tumors had a worse median survival than HOXA9low tumors (Figure 1C). In terms of protein expression levels in LIHC and normal tissues, immunohistochemistry results in HPA database showed that HoxA9 was highly expressed in liver cancer (Figure 1D).

HOXA9 expression among groups of patients according to different clinical parameters in TCGA dataset demonstrated that a significant increase in HOXA9 expression was observed in sample types (Figure 1E), tumor grades (Figure 1F), tumor stages (Figure 1G), both sexes (Figure 1H), cancer stages (Figure 1I), TP53-mutant (Figure 1J), histological subtypes (Figure 1K) of LIHC patients compared to corresponding normal controls. These results suggest that HOXA9 is a prognostic indicator for the patients with liver cancer, which are highly aggressive and lethal.

HOXA9 expression in clinical HCC samples

To further confirm the upregulation of HoxA9 gene in HCC tumors, both mRNA and protein levels of HoxA9 were analyzed in 24 paired specimens of tumor and adjacent non-tumor tissues. To verify the expression of HOXA9 in HCC tissues, we investigated the mRNA level of HoxA9 in 12 pairs of tumor and adjacent non-tumor tissues by real-time PCR. Consistently, HOXA9 expression was significantly increased in HCC samples (Figure 2A, B). HOXA9 expression was also significantly increased in human HCC cell lines (HepG2, Hep3B, HCC-LM3, Huh7 and MHCC-97H) campared with normal human hepatocyte cell line (THLE-2) (Figure 2C, D). Then the expression of HOXA9 in tumor tissues and parallel non-tumor tissues resected from patients with HCC were assessed by IHC, which confirmed the protein level of HOXA9 in HCC was higher than that in adjacent tissues (Figure 2E). These results suggest that the expression of HOXA9 is upregulated in HCC.

Fig. 2

HOXA9 promotes tumorigenesis in HCC samples and cell lines. A The expression levels of HOXA9 in tumor tissues (T) and paracancer tissues (P) of 24 HCC patients were compared by Western blot. ImageJ was used to calculate the gray value and compare the HOXA9 Fold for statistical analysis (HOXA9 Fold = HOXA9/β-Actin; P < 0.05); The original blots are presented in Supplementary Fig. 4A; B The RNA expression of HoxA9 in HCC tumor and adjacent non-tumor tissues was detected by qPCR; C Western blotting confirmed that HOXA9 was highly expressed in HCC cell lines compared with Human normal liver cell line THLE-2; The original blots are presented in Supplementary Fig. 4B; D The RNA expression of HoxA9 in HCC cell lines compared with Human normal liver cell line THLE-2 was detected by qPCR. E Immunohistochemical staining and statistical analysis of HOXA9 in paired cancer tissues (T) and paracancer tissues (P) from (number) HCC patients. Scale bars, 200 μm; Results are means ± SD (n = 5 per group); ***P < 0.001, **P < 0.01, *P < 0.05. Statistical analysis was determined by Student’s t-test

HOXA9 promotes cell proliferation of HCC in vitro and in vivo

To explore the functions of HOXA9 in HCC tumorigenesis, we conducted overexpression and knockdown of HOXA9 in HCC cells. HCC cell lines were performed HOXA9 knockdown with lentivirus-based shRNA. On the other hand, we also transfected HOXA9 overexpressing plasmid (pLVX-Flag-HOXA9) into HCC cell lines (Figure 3A). Afterwards, CCK-8 and EdU assays demonstrated that knockdown of HOXA9 significantly inhibited cell proliferation in Hep3B and MHCC-97H cell lines (Figure 3B, C), which was consistent with the results of colony formation assays (Figure 3D). On the contrary, HOXA9 overexpressing in Hep3B and MHCC-97H cells significantly promoted cell proliferation (Figure 3B–D).

Fig. 3

HOXA9 promotes proliferation in HCC cell lines. A The efficiency of transfection of HOXA9 knockdown and overexpressing were detected by western blot. The original blots are presented in Supplementary Fig. 4C. B–C CCK8 assays (B) and EdU asseys (C) were performed to measure cell proliferation ability of HOXA9 overexpression and shRNA knockdown groups in MHCC-97H and Hep3B cells compared with the control groups. Data are shown as means ± SD (n = 3). D Colony formation assays were performed to measure cell proliferation ability of HOXA9 overexpression and shRNA knockdown groups in MHCC-97H and Hep3B cells compared with the control groups. Data are shown as means ± SD (n = 3). E Images of nude mice with tumors formed from MHCC-97H cells with shRNA knockdown and Hep3B cells with HOXA9 overexpression groups compared with the control groups. F The weight of tumor samples in each group formed from MHCC-97H and Hep3B cells mentioned above by the time of harvest. G The tumor growth curves were measured after the injection of MHCC-97H and Hep3B cells in each group once the tumor had formed, and the volume was calculated every 7 days. Data are shown as means ± SD (n = 5 per group); ***P < 0.001, **P < 0.01, *P < 0.05. Statistical analysis was determined by Student’s t-test

To investigate the effect of HOXA9 on HCC tumorigenesis in vivo, we delivered a xenograft tumor model in nude mice (Figure 3E). As observed, knockdown of HOXA9 significantly suppressed tumor growth and shHOXA9 groups gained lighter weights (Figure 3F) and smaller volumes (Figure 3G) in tumors formed from MHCC-97H cells. Besides, HOXA9 overexpressing groups showed opposite results in tumors formed from Hep3B cells. Overall, HOXA9 plays a promoting role in HCC proliferation both in vitro and in vivo.

HOXA9 propels HCC metastasis and prevents apoptosis in vitro

The wound healing and transwell assays of HCC cells confirmed that silencing of HOXA9 significantly weakened the cell migration ability in wound healing assay (Figure 4A) and migration and invasion cells in HOXA9 knockdown groups were fewer than that in the control group (Figure 4B). By contrast, overexpressing of HOXA9 showed an opposite trend. Overall, these results suggest that HOXA9 may play a promoting role in HCC metastasis in vitro.

Fig. 4

HOXA9 promotes migration, invasion and prevents apoptosis in HCC cell lines. A–B Representative images and quantification of the effect of HOXA9 overexpression and shRNA knockdown groups on cell migration or invasion through wound healing assay (A) and transwell assay (B) in MHCC-97H and Hep3B cells compared with the control group. Scale bar, 200 μm. Results are means ± SD (n = 3 per group); C Flow cytometry analysis was performed to measure cell apoptosis rates of HOXA9 overexpression and shRNA knockdown groups in MHCC-97H and Hep3B cells compared with the control group. Data are shown as means ± SD (n = 3); D The BAK, BCL-2, Caspase3, Cleaved-Caspase3 and β-Actin were detected by Western blot analysis in HOXA9 knockdown and overexpressing HCC cells (Hep3B and MHCC-97H). The original blots are presented in Supplementary Fig. 4D; ***P < 0.001, **P < 0.01, *P < 0.05. Statistical analysis was determined by Student’s t-test

Meanwhile, flow cytometric analysis in HCC cells confirmed that the apoptosis rates of shHOXA9 groups in HCC cells were significantly increased (Figure 4C). Besides, intended proteins were detected by Western blot analysis. Pro-apoptosis protein BAK was upregulated and anti-apoptosis protein BCL-2 was downregulated in shHOXA9 groups of HCC cells (Figure 4D). Meanwhile, a rise of protein levels of Caspase3 and Cleaved-Caspase3 was also observed in shHOXA9 groups (Figure 4D). By contrast, HOXA9 overexpressing groups showed opposite results. Thus, the data above suggested that HOXA9 prevents apoptosis in HCC cells.

HOXA9 expression is regulated by RPL38 in HCC

Ribosomal Protein L38 (RPL38) has been known to regulate Homeobox (Hox) genes’ expression via IRES dependent translation (Kondrashov et al. 2011; Xue et al. 2014). To investigate whether RPL38 regulates HOXA9 expression in HCC, RPL38 was ablated in two human HCC cell lines (Hep3B and HCC-LM3) or overexpressed in human HCC cell line MHCC-97H (Figure 5A, B). Whereas RPL38 ablation decreases the protein level of HOXA9, overexpression of RPL38 increases it (Figure 5A, B). Curiously, the mRNA level of HOXA9 was also regulated by RPL38 whose expression positively correlated with HOXA9 mRNA level (Figure 5C). Taken together, these results suggest that both mRNA and protein levels of HOXA9 were controlled by RPL38.

Fig. 5

RPL38 regulates the expression of several Homeobox genes, including HoxA9. A HCC cell lines (Hep3B and HCC-LM3) were performed RPL38 knockout with CRISPR/Cas9 plasmids compared with the control groups. The HOXA9, RPL38 and α-Tublin were detected by Western blot analysis in RPL38-silenced HCC cells (Hep3B and HCC-LM3). The original blots are presented in Supplementary Fig. 4E; B MHCC-97H cell line was transfected with the RPL38 overexpressing plasmid (pLVX-Flag-RPL38) compared with the control group. The HOXA9, RPL38 and α-Tublin were detected by Western blot analysis in MHCC-97H RPL38-overexpressing cell line. The original blots are presented in Supplementary Fig. 4E; C The RNA expression of HoxA9 in RPL38-silenced Hep3B and HCC-LM3 cells and MHCC-97H RPL38-overexpressing cell line were detected by qPCR; Data are shown as means ± SD (n = 3); D Hierarchical clustering analysis of differentially expressed genes between RPL38-silenced (Cas9-RPL38 knockout) and control of HCC-LM3 cells. A heat map of selected Homeobox genes with significent diffferences was made; E The RNA expression of HoxA7, HoxA10, HoxA13 in RPL38-silenced Hep3B and HCC-LM3 cells were detected by qPCR; Data are shown as means ± SD (n = 4); F The RNA expression of HoxA7, HoxA10, HoxA13 in MHCC-97H RPL38-overexpressing cell line were detected by qPCR; Data are shown as means ± SD (n = 4); G The RNA expression of HoxA1, HoxA2, HoxA5 in RPL38-silenced Hep3B and HCC-LM3 cells were detected by qPCR; Data are shown as means ± SD (n = 4); H The RNA expression of HoxA1, HoxA2, HoxA5 in MHCC-97H RPL38-overexpressing cell line were detected by qPCR; Data are shown as means ± SD (n = 4); ***P < 0.001, **P < 0.01, *P < 0.05. Statistical analysis was determined by Student’s t-test

Meanwhile, RNA-sequencing was conducted in parental and RPL38-ablated HCC-LM3 cells. Intriguingly, the mRNA levels of multiple Hox genes were dramatically upregulated or downregulated by RPL38 ablation (Figure 5D). qRT-PCR analysis of several Hox genes in human cancer cell lines or normal human hepatocytes showed that while RPL38 expression positively correlated with the mRNA levels of HOXA7, HOXA10, and HOXA13 genes (Figure 5E, F), it has a reversed effect on the expression of HOXA1, HOXA2 and HOXA5 genes (Figure 5G, H). These results suggest that in addition to upregulation of HOXA9 expression, RPL38 displays opposite effects on distinct Hox genes’ expression.

Identification of HOXA9-interacting genes and proteins analysis

To better explore the underlying pathogenic mechanism of HOXA9 in LIHC, the coexpressed genes of HOXA9 were screened out in TCGA-LIHC patients using the LinkedOmics database (Figure 6A). The top 50 genes that were positively and negatively correlated with HOXA9 in LIHC are shown (Figure 6B, C).

Fig. 6

HOXA9-related genes and enrichment analysis. A Coexpressed genes associated with HOXA9 in HCC, analyzed by using LinkedOmics and visualized by Volcano plot; B Heat maps showing the top 50 genes positively correlated with HOXA9 in HCC; C Heat maps showing the top 50 genes negatively correlated with HOXA9 in HCC; D GO enrichment analysis of genes correlated with HOXA9, including biological process, cellular component and molecular functions; E KEGG enrichment analysis of top genes correlated with HOXA9

Notably, in terms of cellular composition (CC), and molecular function (MF), HOXA9 was enriched in DNA methylation processes, such as methyltransferase complex and methyl-CpG binding in LIHC (Figure 6D). In addition, the top 10 KEGG pathways for HOXA9 and its correlated genes are shown in Figure 6E. Among these pathways, some immune-related pathways were highly associated with HOXA9, including T cell receptor signaling pathway, Th1 and Th2 cell differentiation and Natural killer cell mediated cytotoxicity in LIHC (Figure 6E). These results strongly imply that HOXA9 is involved in the regulation of DNA methylation and immune infiltration in liver cancer.

Correlations between HOXA9 and immune cell infiltration in the tumor microenvironment (TME)

The tumor microenvironment (TME) plays a critical role in the development and progression of HCC. Studies indicate that HCC employs multiple immunomodulatory mechanisms to disrupt the immune system (Flecken et al. 2012). We further explored the correlations between HOXA9 expression and immune cells in the TME that HOXA9 expression presented strong connections to CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells (DCs) in LIHC (Supplementary Figure 2A). Then, a heatmap showed correlation between HOXA9 expression and the above 6 immune cells in TME of TCGA-LIHC dataset (Supplementary Figure 2B). Besides, there was a statistically significant positive correlation between HOXA9 expression and the immune infiltration levels of three types of immune cells: CD4+ T cells (R = 0.21, P = 8.2e-5), neutrophils (R = 0.18, P = 4.0e−4) and macrophages (R = 0.55, P = 0.03) for LIHC cases (Supplementary Figure 2C). These results showed close correlations between HOXA9 and immune cell infiltration in the tumor microenvironment (TME).

To deepen our understanding of HOXA9 with the immune response, we validated the correlations between HOXA9 expression and diverse immune signatures in LIHC. After adjusting for tumor purity, HOXA9 expression was significantly associated with most immune markers of immune cells in LIHC (Supplementary Table 3). We also examined the correlation between HOXA9 expression and various functional T cells and HOXA9 expression was significantly correlated with 32 of 38 T-cell markers in LIHC after adjusting for tumor purity (Supplementary Table 4). These findings further support that HOXA9 expression is significantly related to immune infiltration and HOXA9 might be vital in regulating immune cell infiltration.

DNA methylation analysis of HoxA9 gene

The validation of the DNA methylation status from TCGA data showed that the promoter methylation level of HOXA9 was significantly elevated in LIHC tissues of patients in terms of sample types, tumor grades, tumor stages, both sexes, cancer stages, TP53-mutant and TP53 wild-type LIHC patients compared to normal controls (Supplementary Figure 3A–F). To summary, these results suggest that DNA methylation of HOXA9 may be a potienal prognostic indicator for LIHC patients.