HOX Gene Family Denes Tumor Microenvironment and Drug Sensitive in Hepatocellular Carcinoma

Hepatocellular carcinoma (HCC) pathophysiology is prevalently related with HOX genes. However, the study on associations of extensive HOX genes with tumor microenvironment and drug sensitivity of HCC remains scarce. The data sets of HCC were downloaded from TCGA, ICGC and GEO by bioinformatics method and analyzed. RT-qPCR and IHC were used to verify the expression levels of 10 HOX genes in normal and HCC. Based on unsupervised clustering, two distinct HOX genes patterns has been identied, and a signicant shorter OS in clusterB from HCC patients was observed than those in clusterA. Furthermore, based on a computational frame, we divided HCC samples into high and a low HOXscore groups, and signicantly shorter survival time in the high HOXscore was observed relative to low HOXscore group using survival analysis. GSEA revealed that the high HOXscore group was more likely to be enriched in cancer-specic pathways. Furthermore, the high HOXscore group was associated with the inltration of inhibitory immune cells. Besides, the HOXscore has greater prognostic signicance than TMB alone and is capable of identifying intermediate subtype group existing with partial TMB functionality in low TMB patients. In response to anti-cancer drugs, the high HOXscore group was more sensitive to mitomycin and cisplatin. Importantly, the HOXscore was associated with the therapeutic ecacy of PD-L1 blockade, suggesting that the development of potential drugs targeting these HOX genes to aid the clinical benets of immunotherapy. In addition, RT-qPCR and immunohistochemistry showed 10 HOX genes mRNA expression was higher in HCC compared to the normal tissues. Finally, high HOXscore was closely related with poorer overall survival of HCC in an independent validation cohort. Our study provide a comprehensive analysis of HOX genes family in HCC. We revealed the potential function of these HOX genes family in TME and identied their therapeutic liability in targeted therapy and immunotherapy. This work highlights the cross-talk and potential clinical utility of HOX genes family in HCC therapy.


Introduction
Hepatocellular carcinoma (HCC) is one of the most common malignant tumors, and has aroused more and more attention with its high rate of recurrence and distant metastasis [1,2]. Its recurrence rate is up to 70%, and the 5-year survival rate is less than 60% with successful surgical resection or liver transplantation [3]. Due to the lack of speci c symptoms of early HCC, it is di cult to be detected, and 90% of patients are already in the advanced stage at the time of diagnosis [4]. Inspite of achieving great progress in prevention, diagnosis, and treatment for this disease recent years, the clinical outcome is still unsatisfactory [5]. There are limited therapeutic methods for HCC. Surgery and percutaneous ablation can achieve long-term control in some early HCC patients. However, since most of the early HCC patients have no obvious symptoms, most of the patients are in the late stage and cannot be resected by radical surgery [6]. Including arterial chemoembolization, radiotherapy, immunotherapy or hormone therapy, the treatment results are poor, and the ve-year survival rate of patients is extremely low [7,8]. Therefore, the search for HCC prognostic markers and drug sensitivity may provide relevant guidance for the survival and treatment of HCC patients.
The pathogenesis of HCC is heterogeneous and complex, which may be a process involving multiple genes and multiple signaling molecules, in which gene mutation, activation of oncogenes and inactivation of tumor suppressor genes play an important role [2]. HOX gene is a family of genes encoding transcription factor (TF), which plays an important role in embryonic development and cell differentiation.
The HOX gene family contains a common 120-base-pair DNA sequence that encodes a 61-amino acid peptide called the homologous domain. In the human genome, the HOX gene family contains 39 members, located on four different chromosomes [9]. In recent years, studies have found that the HOX gene family is involved in the development of a variety of tumors [10]. HOX gene can not only be upregulated as a transcriptional activator in cancer, but also play a role as a transcriptional suppressor. The up-regulation or down-regulation of TF encoded by HOX gene family in cancer is critical to the development and progression of cancer. This abnormal regulation of the HOX gene in cancer suggests that HOX plays a tissue-speci c role in tumors [11]. Based on the crucial role in modulating the initiation and progression of cancer, numerous investigators utilized these HOX genes to establish prognostic models for tumor diagnosis or prognosis,, even the guidance of therapeutic strategies. Similarly, HOX genes were recently used to develop an optimal prognostic gene signature for predicting the prognosis of cancer [12,13].
Thus far, the mechanisms underpinning the role of HOX genes in HCC remain mostly unknown. In this study, based on unsupervised clustering, two distinct HOX genes patterns has been identi ed, characterized under these two patterns were highly consistent with the three immune phenotypes of tumors. Next, based on a computational frame, we divided HCC samples into high and a low HOXscore groups, and signi cantly shorter survival time in the high HOXscore was observed relative to low HOXscore group using survival analysis. Besides, the HOXscore has greater prognostic signi cance than TMB alone and is capable of identifying intermediate subtype group existing with partial TMB functionality in low TMB patients. In response to anti-cancer drugs, the HOXscore group was more sensitive to mitomycin and cisplatin. Importantly, the HOXscore was associated with the therapeutic e cacy of PD-L1 blockade.

Data source and preprocessing
The work ow of our study was shown in Figure 1A. HCC RNA-seq data (FPKM value) were obtained from TCGA, ICGC, and GSE76427 database [14], and their clinical information also was downloaded simultaneously. We obtained 725 HCC tissue samples. Considering all data in this study was obtained from PDB Public DataBase, the ethics approval and informed consent were not provided. The detailed clinical characteristics of these patients are summarized in Table 1.

Selection of HOX genes in HCC
Limma packages in R language was used screening the differentially expressed HOX genes in HCC. The default Benjamini-Hochberg false discovery rate (FDR) was employed to adjust the p-value to remove false positive results. P<0.05 and and absolute log2FC>1 was considered a criterion for the selection of differentially expressed genes. Differently expressed HOX genes were processed by hierarchical clustering analysis with R package ''pheatmap'' [15].

Consensus clustering analysis
The cohort was clustered into subgroups according to the consensus expression of HOX genes using "ConsensusClusterPlus" in R. Subsample 80% of samples at each iteration and partition each subsample into up to k (max K = 9) groups by k-means algorithm upon Euclidean distance. This process was repeated for 1,000 repetitions [16]. Sequentially, we selected a perfect clustering result by considering the relative change in area under the cumulative distribution function curve, which was assessed by the proportion of ambiguous clustering (PAC) score, a low value of PAC implies a at middle segment, allowing conjecture of the optimal K by the lowest PAC. A HOXscore for each patient was calculated according to the following formula: HOXscore = ∑(PC1 i + PC2 i ), where i is the TPM value of each screened gene. The overall survival (OS) discrepancy among subclusters was assessed though the Kaplan-Meier method and log-rank test [17]. Principal component analysis (PCA) was used to exhibit the gene expression patterns in the three cohort subclusters using the R package [18]. The differences of clinic pathological features, including age, gender, fustat and stage, were calculated by Mann-Whitney-Wilcoxon test between the subgroups.
Gene set enrichment analysis (GSEA) In order to further investigate the potential role of HOXscore in HCC, we performed gene enrichment analysis. The HOXscore was divided into high HOXscore group and low HOXscore group based on the median value in samples. The go and kegg standardized pathway gene set obtained from the GSEA website which was used for enrichment analysis. We sequenced 1000 gene sequences and sequenced the enrichment pathways of each phenotype using nominal p value and normalized In reviewenrichment score (NES). The discovery rate |NES|>1, NOM p-val < 0.05, FDR q-val< 0.25 were considered to be signi cantly enriched [19].

Evaluation of immunity in strati ed three cohorts
We performed the ssGSEA method to calculate the enrichment scores on the basis of metagenes [20]. The reason why we considered the metagene a robust approach was the two main characteristics: (1) the use of a set of genes instead of single genes that represent one immune subpopulation, because the use of single genes as markers for immune subpopulations can be misleading as many genes are expressed in different cell types, and (2) the assessment of relative expression changes of a set of genes in relation to the expression of all other genes in a sample. Referring [21]. Immunohistochemistry 150 HCC patients tissues and corresponding adjacent tissues were collected to explore the expression of 10 HOX genes in the tissue samples by using immunohistochemical staining (IHC). IHC staining was performed according to the manufacturer's instructions.

Statistical analysis
HOX gene expression and clinical data obtained from TCGA, ICGC, and GSE76427 were combined through perl5.30.0.1. R4.0.3 was used to screen differential HOX genes of HCC, analyze independent prognosis, and search for clinical correlation. In addition, genetic and clinical related data were used to perform univariate Cox regression analysis. The data difference between two groups was calculated by Mann-Whitney-Wilcoxon test. P 0.05 was considered as statistical signi cance.

Result
Identi cation of HOX genes in TCGA-LIHC patients In order to visualize the expression pattern of HOX genes between the HCC and normal control, heatmap of 39 HOX genes were obtained in Figure 1A In order to predict the clinical outcomes of HCC with HOX genes, univariate Cox regression analysis was performed to identify the HOX genes that associated with OS in TCGA-LIHC cohort. 11 HOX genes were signi cantly correlated with OS ( Figure 1B). Unreasonably, HOXB1 has no signi cant between HCC samples and normal control samples, but its has signi cantly correlated with OS in the univariate Cox analysis, so it was excluded from further study ( Figure 1C). A total of 10 prognostic HOX-related DEGs were preserved, and Figure 1D-M show the expression distribution of 10 HOX genes between paired normal and HCC tissues ( Table 2). Unsupervised clustering of 10 HOX genes in the HCC meta-cohort Next, Three datasets (TCGA-LIHC, LIRI-JP, and GSE76427) with available OS data and clinical information were enrolled into one meta-cohort. Based on the expression similarity of the 10 HOX genes (HOXD13, HOXB1, HOXD11, HOXD10, HOXA6, HOXA7, HOXC8, HOXA11, HOXA9, HOXC6, and HOXD9), the 725 HCC patient cohort was divided into two clusters (cluster A and cluster B) when k = 2 was used for consensus clustering analysis ( Figure 3A). A signi cant shorter OS in clusterB from HCC patients was observed than those in clusterA ( Figure 3B). To explore functional role of the 2 distinct HOX genes patterns in HCC progression, we performed GSVA analysis, and the results indicated that clusterB was more likely to be enriched in cancer-speci c pathways, such as cell cycle, RNA polymerase, RNA degradation, spliceosome, and nucleotide excision repair, and clusterA was more likely to be enriched in cell adhesion molecules, ECM receptpr interction, dilated cardiomyopathy, and hypertrophic cardiomyopathy ( Figure 3C). Moreover, the principal component analysis (PCA) analysis of the transcriptional pro le between cluster A and cluster B denoted a markedly discrimination each other ( Figure 3F). Then, we assessed the associations between the clustering and clinicpathological features. There were signi cant differences between the cluster A and cluster B for the N, stage, gender and age (P < 0.05), and fustat (P < 0.01), and the signi cant overexpressions of 10 HOX genes in cluster B samples were observed relative to cluster A samples ( Figure 3D). It has been con rmed that the tumor microenvironment (TME) is correlated with the sensitivity to immunotherapy and HCC prognosis [22], and ssGSEA method was used to calculate the enrichment scores on the basis of metagenes, and was employed in the gene expression data from three datasets to calculate immune score. As shown in gures 3E, the proportion of immune cell in the clusterB group was signi cantly lower than that in the clusterA group (P <0.05), suggesting that the immune ability of clusterB is suppressed.

Generation of 10 HOX gene signatures and functional annotation
To further investigate the potential biological behavior of each HOX gene pattern, we determined 1417 HOX phenotype-related DEGs using limma package. Next, the GSEA was used to reveal the malignant hallmarks of HCC, and GSEA GO results indicated that immune regulation was signi cantly associated with HCC ( Figure 3G), GESA KEGG results indicated that metabolism pathway was signi cantly associated with HCC ( Figure 3H). To further validate this regulation mechanism, we consensus clustering of 1417 HOX phenotype-related DEGs, and the 725 HCC patient cohort was divided into three clusters (cluster A, cluster B, and cluster B) when k = 3 was used for consensus clustering analysis ( Figure 4A). A signi cant shorter OS in clusterC from HCC patients was observed than those in clusterA and clusterB ( Figure 4B). Then, we also observed that there were signi cant differences between the cluster A, cluster B, and cluster C for the N, stage, gender and age, and fustat (P < 0.05) ( Figure 4C), and the signi cant overexpressions of 10 HOX genes in cluster C samples were observed relative to cluster A and B samples ( Figure 4D). The above results showed suggested HOXgene played a role in TME, however, these analyses were only based on the patient population and could not accurately predict the pattern of HOX gene in individual patients.
Based on these DEGs, we constructed a set of scoring system to quantify the HOX gene pattern of individual patients with HCC, we termed as HOXscore, and a signi cant shorter OS in high HOXscore from HCC patients was observed than those in low HOXscore ( Figure 4E), and the positive correlation between HOXscore and immune cell was also observed ( Figure 4F). We also used alluvial diagram to showed the relationship between gene cluster, HOX cluster, HOXscore, and fustat, and con rmed this relationship ( Figure 4G). We further analyzed the HOXscore value between the HOX clusterA and B group, and signi cantly higher HOXscore value in the clusterB group was observed relative to clusterA group ( Figure  4H). In addition, we analyzed HOXscore value between the gene clusterA, B, and C group, and signi cantly higher HOXscore value in the clusterC group was observed relative to clusterA and B group ( Figure 4I).

Independent prognostic analysis of HOXscore
Next, multivariate Cox regression analysis was performed on clinical factors and HOXscore to determine the independent prognostic value of HOX ( Table 3). As shown in gures 5a and 5b, in both the young and old groups, signi cantly shorter survival time in the high HOXscore group was observed relative to low HOXscore group using survival analysis (P <0.05) ( Figure 5A and B). Subsequently, signi cantly shorter survival time in the high HOXscore group was observed relative to low HOXscore group using survival analysis in the male, female, stage1-2, and stage 3-4 (P <0.05) ( Figure 5C-G). Mitomycin and cisplatin are commonly used chemotherapy drugs for HCC, and we found that the IC50 to mitomycin and cisplatin in the high HOXscore group was lower than that in the low HOXscore group, suggesting that the high HOXscore group was more sensitive to mitomycin and cisplatin. Abovementioned evidence demonstrated that the prognostic signature based on HOXscore is reliable in HCC prognosis. To explore functional role of this prognostic signature in HCC progression, we performed GSEA analysis, and the results indicated that the high HOXscore group was more likely to be enriched in cancer-speci c pathways ( Figure 6A and B).
The HOXscore is better than tumor somatic mutation (TMB) can predict HCC prognosis It has been con rmed that the abnormality of TMB is correlated with the sensitivity to immunotherapy and HCC prognosis [23]. We further analyzed the HOXscore value between the gene clusterA, B, and C group, and signi cantly higher HOXscore value in the clusterC group was observed relative to clusterAand B group ( Figure 7A). In addition, we analyzed the HOXscore value between the HOX clusterA and B group, and signi cantly higher HOXscore value in the clusterB group was observed relative to clusterA group ( Figure 7B), and a signi cant shorter OS in high TMB from TCGA-LIHC patients was observed than those in low TMB (( Figure 7C). Next, to evaluate whether HOXscore is better than TMB status can predict HCC prognosis, we divide all patients into H-TMB/H-HOXscore, H-TMB/L-HOXscore, L-TMB/L-HOXscore, and L-TMB/L-HOXscore. The results showed that patients in L-TMB/H-HOXscore had a poor survival as compared with L-TMB/L-HOXscore in TCGA set ( Figure 7D), suggest the HOXscore has greater prognostic signi cance than TMB alone and is capable of identifying intermediate subtype group existing with partial TMB functionality in low TMB patients. Furthermore, we analyzed the somatic mutation between the high and low HOXscore group, and no signi cant was observed between the two group ( Figure 7E).
HOX patterns in the role of PD-1/L1 immunotherapy HCC is immunogenic and immunosuppressed, thus, the application of immunotherapy is of great signi cance for the development of THE treatment of HCC. We rst analyzed the expression of PD-L1 between high and low HOXscore group, and signi cantly higher PD-L1expressionin the high HOXscore group was observed relative to low HOXscore group ( Figure 8A). Next, we download the immunotherapy data of TCGA-LIHC patients, and the signi cant therapeutic advantages and clinical response to PD-1/L1 immunotherapy in patients with high HOXscore compared to those with low HOXscore were con rmed ( Figure 8B-E). It has been con rmed that the abnormality of MSI is correlated with the sensitivity to immunotherapy. We further analyzed the MSI between high and low HOXscore group, and signi cantly higher MSI in the high and low HOXscore group was observed relative to low HOXscore group ( Figure 8F). In addition, we analyzed the HOXscore value between the MSS, MSI-L, and MSI-H group, and signi cantly higher HOXscore value in the MSI-H group was observed relative to MSS and MSI-L group ( Figure 8G). Abovementioned evidence demonstrated that HOXscore is reliable in HCC PD-1/L1 immunotherapy.

Validation the expression and prognosis of HOXscore in GSE76427
To further validate HOX genes expression in HCC, GSE76427 was used to measure the expression level of HOX genes, and the result showed that compared with normal group, the HOXA6, HOXA7, HOXA9, HOXA11, HOXC6, HOXC8, HOXD10, HOXD11, HOXD13, and HOXD9 level were signi cantly higher in HCC group. (Figure 9A-J). In addition, signi cantly shorter survival time in the high HOXscore group was observed relative to low HOXscore group using survival analysis ( Figure 9K).
Next, to further validate HOXscore in HCC, we measured the 10 HOX gnes protein level by immunohistochemistry, and the result showed that compared with the normal group, the 10 HOX gnes protein level was signi cantly higher in HCC group ( Figure 10A). In addition, RT-qPCR was used to detect the 10 HOX gnes mRNA expression in HCC. Compared with the normal group, the 10 HOX gnes mRNA level was signi cantly higher in the HCC group ( Figure 10B). Additionally, Kaplan-Meier and Cox's proportional hazards regression model survival analysis revealed that patients with high HOXscore group had shorter overall survival in HCC ( Figure 10C).

Discussion
Increasing evidence demonstrated that HOX genes took on an indispensable role in embryonic development and cancer. While most studies have focused on a single type of HOX gene, the mutual relationships and functions of multiple types of HOX genes in cancer are not fully understood. Here, we selected 39 HOX genes in this study, and 33 aberrant expression of HOX genes in HCC were observed. 11 HOX genes were signi cantly associated with OS of TCGA-LIHC cohort. Of these, 10 HOX genes were used to divided HCC patient into two clusters, and a signi cant OS between the two clusters. The abundance of immune cells in the TME was signi cantly different between the two HCC subtypes, speci cally, the immune ability of clusterB is suppressed. Base on the DEGs between two HCC subtypes, the HCC patient divided into three clusters, and a signi cant OS between the three clusters.
Simultaneously, based on these DEGs, we constructed HOXscore, and HOXscore is an independent prognostic in HCC. Moreover, HOXscore has greater prognostic signi cance than TMB alone and is capable of identifying intermediate subtype group existing with partial TMB functionality in low TMB patients. Finally, HOXscore is reliable in HCC PD-1/L1 immunotherapy, and sensitive to mitomycin and cisplatin.
HCC is one of the most popular common forms of liver cancer with the four highest cancer-associated mortality worldwide. Its initiation and progress is a complicated pathomechanism, which involves the epigenetics alterations at the posttranscriptional level [2]. HOX is abnormally expressed in a variety of tumors and is closely related to tumor stage and prognosis [24]. The regulation of HOX gene expression is a major factor affecting the e cacy of HOX in tumors, involving chromosome epigenetic modi cation and signaling molecules, such as FGF, Wnt, RA and a variety of steroid small molecules [25][26][27]. Meanwhile, the transcription activity of HOX protein is regulated by cofactors and other transcription factors [28]. HOX regulates the proliferation, apoptosis, migration and invasion of tumor cells by promoting the transcription of target genes, and participates in the occurrence and development of tumors and drug tolerance [29,30]. The role of HOX in various HCC has been reported. HOXA7 upregulates the expression of cyclin E1 and CDK2 proteins, thereby promoting the proliferation of HCC cells [31]. HOTAIR gene is located in the HOXC gene cluster and inhibits the expression of HOXD cluster genes such as HOXD10 by recruiting methyltransferase [32]. The expression of HOXD10 was downregulated in liver cancer tissues and HepG2 cells and promoted the invasion and metastasis of cells. In HCC, the expression of HOXD10 is down-regulated due to promoter methylation, which promotes the survival, migration and invasion of HCC cells. The HOXD10 downstream gene IGFBP-3 induces apoptosis of HCC cells by activating caspase-3 and Caspase-8 [33].
Recently, immunotherapies, especially immune checkpoint inhibitor therapy through targeting the PDL1/PD-1 axis, have been focused in the eld of HCC therapy [34]. Bewilderingly, less than 20% of cancer patients obtain the bene cial effects following treatment with immune checkpoint inhibitor, whereas the remainders are frustrated. It may be explained that different subtypes in homogeneous cancer exhibit disparate sensitivity to immunotherapy [35]. In HCC, previous study has constructed a prognostic model on the basis of immune signatures, which suggested that the high immunity subtype shows longer OS and more sensitivity to immunotherapy [36]. In our study, the similar tendency was observed that the high HOXscore group with lower immune and stromal scores exhibited shorter OS. Meanwhile, the increased TMB and tumor purity were found in high HOXscore group, which indicated a greater likelihood of malignant behavior. These evidences support that HOXscore can be helpful for guiding the selection of immunotherapeutic strategies in HCC.
There have been many studies on the composition and prognosis of immune cells in the tumor microenvironment. M2 macrophages and proin ammatory M1 macrophages are associated with the growth and spread of tumors and the formation of immunosuppressive microenvironment [37,38]. High concentrations of CD3 + T cells, CD8 + cytotoxic T cells, and CD45RO + memory T cells were signi cantly associated with longer disease-free survival and overall survival [39,40]. The effect of CD4 + regulatory T cells on tumor has been debated, and studies have shown that there was no effect on survival while other studies showed a negative association with prognosis in head and neck cancer, colorectal cancer, bladder cancer, and ovarian cancer [41][42][43]. Studies have used CD57 as a phenotypic marker to detect NK cells in liver cancer tissues, showing that THE in ltration of CD57 + cells indicates a good prognosis [44]. The role of B cells in tumors is unclear. The spontaneous cancer model in mice suggests that B cells have a carcinogenic effect, which may activate M2 macrophages and promote the early stage of carcinogenesis, and may also promote metastasis by promoting the transformation of resting T cells into regulatory T cells [45,46]. However, other studies have shown that B cells can act as antigen-presenting cells, inducing CD4 + T cell-dependent CD8 + memory T cells, thereby contributing to the control of tumor invasion and metastasis [47]. In our study, we found most of the immune cell activity is suppressed in HCC with low survival rate, and the positive correlation between HOXscore and immune cell was also observed.
There are several limitations in this study. First, there is a lack of biological veri cation. Future molecular studies are needed to identify the interactions between 10 HOX genes from our prognostic model. Moreover, all data in this study is obtained from public databases, so future studies with prospective validation are still warranted.
In conclusion, this work is the rst to provide a comprehensive analysis of HOX genes family in HCC. We revealed the potential function of these HOX genes family in TME and identi ed their therapeutic liability in targeted therapy and immunotherapy. This work highlights the cross-talk and potential clinical utility of HOX genes family in HCC therapy, which also provide a wider insight for understanding the HCC.

Declarations
Authors' contributions Dan Li and Tao Yu designed the study, performed statistical analysis, and drafted the manuscript. Zhi Zeng, Jingjing Han, and Jie Wu helped to draft the manuscript. Benhong Zhou conceived the study, participated in its design and coordination, and helped to draft the manuscript. All authors read and approved the nal manuscript.