Molecular Analysis of Prognosis and Immune Infiltration of Ovarian Cancer Based on Homeobox D Genes

Background Homeobox D (HOXD) genes were associated with cancer pathogenesis. However, the role of HOXD genes in ovarian cancer (OC) and the possible mechanisms involved are unclear. In this study, we analyzed the function and regulatory mechanisms and functions of HOXD genes in OC based on comprehensive bioinformatics analysis. Methods Expression of HOXD1/3/4/8/9/10/11/12/13 mRNA was analyzed between OC tissue and normal tissue using ONCOMINE, GEO, and TCGA databases. The relationship between HOXD expression and clinical stage was studied by GEPIA. The Kaplan-Meier plotter was used to analyze prognosis. cBioPortal was used to analyze the mutation and coexpression of HOXDs. GO and KEGG analyses were performed by the DAVID software to predict the function of HOXD coexpression genes. Immune infiltration analysis was used to evaluate the relationship between the expression of HOXD genes and 24 immune infiltrating cells. Results The expression of HOXD3/4/8/9/10/11 was significantly lower in OC tissues than in normal ovarian tissues, while the expression of HOXD1/12/13 was significantly higher in OC tissues. The expression of HOXD genes was associated with FIGO stage, primary therapy outcome, tumor status, anatomic neoplasm subdivision, and age. The expression levels of HOXD1/3/4/8/9/10 correlated with tumor stage. HOXD1/8/9 could be served as ideal biomarkers to distinguish OC from normal tissue. Low HOXD9 expression was associated with shorter overall survival (OS) (HR: 0.75; 95% CI: 0.58–0.98; P = 0.034) and progression-free survival (PFS) (HR: 0.69; 95% CI: 0.54–0.87; P = 0.002). The HOXD coexpression genes were associated with pathways including cell cycle, TGF-beta signaling pathway, cellular senescence, and Hippo signaling pathway. HOXD genes were significantly associated with immune infiltration. Conclusion The expression of HOXD genes is associated with clinical characteristics. HOXD9 is a new biomarker of prognosis in OC, and HOXD1/4/8/9/10 may be potential therapeutic targets. The members of the HOXD genes may be the response to immunotherapy for OC.


Introduction
Ovarian cancer (OC) is one of the most common gynaecological tumors, ranking fourth in incidence and third in mortality worldwide [1,2]. In China, OC has the second highest mortality rate among gynaecological tumors and is on the rise, while the incidence is declining [3]. It is difficult to detect at an early stage, and most patients are diagnosed at an advanced stage [4]. Despite advances in the treatment of OC with chemotherapy, radiotherapy, surgery, and targeted therapies, the 5-year OS rate for patients with advanced OC is around 30% [5,6]. Therefore, there is a need to explore the genetic signature of prognostic prediction associated with the underlying mechanisms of OC progression.
Homeobox genes are regulatory genes that share a common 180-183 bp sequence and encode a 61-amino acid structural domain known as the homeodomain. This homeodomain is a DNA binding domain that functions as a transcription factor [7]. In humans, the HOX genes are divided into four clusters (HOXA, HOXB, HOXC, and HOXD) on different chromosomes [8]. HOXDs contain 9 members, including HOXD1, HOXD3, HOXD4, HOXD8, HOXD9, HOXD10, HOXD11, HOXD12, and HOXD13. HOXD1 inhibited cell proliferation, cell cycle, and TGF-β signaling in kidney renal clear cell carcinoma (KIRC) [9]. By identifying the YY1-HOXD3-ITGA2 regulatory axis as a potential therapeutic target for hepatocellular carcinoma (HCC) treatment, a new and complete pathway for HCC treatment is offered [10]. Increased expression of HOXD3 was an independent and important predictor of poor prognosis in breast cancer (BRCA) patients [11]. Overexpression of HOXD4 is significantly associated with poorer prognosis in patients with gastric cancer (GC), suggesting the potential of HOXD4 as a novel clinical predictive biomarker and drug target [12]. HOXD8 may be associated with cisplatin resistance and metastasis in advanced OC [13]. Downregulation of miR-142-5p induced resistance to gefitinib in lung cancer PC9 cells through upregulation of HOXD8 [14]. In sum-mary, some members of HOXDs are closely associated with clinical features and drug resistance of tumors, and their expression levels can be used as predictors of tumor prognosis, metastasis, and response to chemotherapy and targeted therapy. HOXD genes play a role in the pathogenesis of pediatric low-grade gliomas [15]. However, the role of HOXDs in OC is unclear. Studying the prognostic value of HOXDs for patients with OC may help to improve the prediction of clinical prognosis in OC and inform personalized treatment.

Correlation Heat Map.
Correlation between every two genes of HOXDs was assessed using a Pearson's correlation coefficient [16]. Software: R (version 3.6.3). R package: mainly ggplot2 (version 3.3.3). Data: RNAseq data in level 3 HTSeq-FPKM format from the TCGA OC project. Data conversion: RNAseq data in FPKM (fragments per kilobase per million) format were converted to TPM format and log2 transformed. Data filtering: remove control/normal (not all items have control/normal).

The Relationship between HOXDs and Clinical
Characteristics of OC. Software: R (version 3.6.3). R package: basic R package. Molecules: HOXD1/3/4/8/9/10/11/12/13. The grouping condition is the median. Data were obtained from the TCGA OC project for RNAseq data in level 3 HTSeq-FPKM format. RNAseq data in FPKM format were converted to TPM format and then log2 transformed.
Expression and correlation analyses of HOXDs were carried out on the GEPIA website (http://gepia.cancer-pku.cn/) [22]. The expression of HOXDs at different clinical stages was generated online.

The Relationship between HOXDs and Prognosis of OC.
Using the Kaplan-Meier method, the analysis was carried out according to the reference [18,23]. Software: R (version 3.6.3). R package: survminer package (for visualization) and survival package (for statistical analysis of survival data). Molecules: HOXD1/3/4/8/9/10/11/12/13. Subgroups: 0-50 vs. 50-100. Prognosis type: OS and progression-free survival (PFS). OS is defined as the time from the beginning to death from any cause. PFS is defined as the time from initiation to the onset of arbitrary tumor progression or the onset of death. Data: RNAseq data and clinical data in level 3 HTSeq-FPKM format from the TCGA OC project. Data filtering: retain data with clinical information. Data conversion: RNAseq data in FPKM format were converted to TPM format and analyzed by grouping them according to     Figure 3: Changes in HOXD mRNA expression between different types of cancer and normal tissues using the ONCOMINE database. Cell color is determined by the best gene rank percentile for the analyses within the cell. Red indicates an increase in expression, blue indicates a decrease in expression, and white indicates that the copy number is neutral. The data in the middle of the square represents the number of data sets. RNAseq data in FPKM format were converted to TPM format and log2 transformed. Supplementary data: prognostic data from the reference [24]. Data filtering: remove control/normal (not all items have control/normal) + keep clinical information.

Correlation
Analysis for Genes Associated with HOXDs in OC. cBioPortal was also used to analyze the relationship between the mutation of HOXDs and survival in OC. Coexpression levels were calculated according to the online instructions of "Similar Genes" part of GEPIA2 (http:// gepia2.cancer-pku.cn/index.html#example#e3). The first 100 coexpressed genes were kept separately for each gene, and finally, all coexpressed genes were summarized. To further validate the accuracy of the ONCOMINE and TCGA databases, OC samples from the GEO database were downloaded for analysis [26]. 10 ovarian cancer tumor tissues and 10 normal ovarian tissues contained in GSE29450 were used for differential gene expression analysis.   Figure 4: The expression of HOXDs in normal ovarian tissue was compared with that of the OC tissues. Significance markers: ns, P ≥ 0:05; * , P < 0:05; and * * * , P < 0:001.
2.11. Statistical Analysis. The methodology of our analysis follows the previous literature [18]. The expression of HOXDs between OC tissue and normal ovarian tissue was analyzed using the Wilcoxon rank sum test. P < 0:05 were considered statistically significant.

HOXD Gene Alterations and mRNA Expression in OC.
The cBioPortal online tool was used to analyze the gene expression of HOXD genes in OC patients. Alterations in the HOXD genes in OC ranged from 4% to 5% (Figure 1). The structural variation data, mutation data, and CNA (copy number alteration) data from the 3 studies are depicted in Figure 2.
The mRNA expression levels of HOXD1/12/13 in OC tissues were significantly higher than that in normal ovarian tissues, and the mRNA expression levels of HOXD4/8/9/10/11 in OC tissues were significantly lower than those in normal ovarian tissues. There was no significant difference in HOXD3. As shown in Table 2, compared with normal ovarian tissues, HOXD3 was significantly lower expressed in OC tumor tissues (fold change = 0:331, P = 0:026), HOXD4 was significantly lower expressed in OC tumor tissues (fold change = 0:155, P < 0:001), HOXD8 was significantly lower expressed in OC tumor tissues (fold change = 0:280, P < 0:001), HOXD9 was significantly lower expressed in OC tumor tissues (fold change = 0:373, P = 0:009), and HOXD12 was significantly higher expressed in OC tumor tissues (fold change = 6:720, P < 0:001). There was no significant difference in HOXD1/10/11/13. The above results from different databases showed that the expression of HOXD4/8/9/10/11 was significantly lower in OC tissues than in normal ovarian tissues, while the expression of HOXD1/12/13 was significantly higher in OC tissues than in normal ovarian tissues. We examined the correlation between HOXD genes using the Pearson correlation analysis. As shown in Figure 5, there        Computational and Mathematical Methods in Medicine was no significant correlation between HOXD12 and HOXD1/3/4. Other HOXD genes were significantly positively correlated with each other.

Relationship between HOXD mRNA Expression and the
Clinical Stage and Prognosis of OC. As shown in Figure 6, HOXD1/3/4/8/9/10 were negatively correlated with the clinical stage of OC. HOXD1/3/4/8/9/10 may be closely related to the development of OC. As shown in Table S1, in the TCGA database, the clinical information of 379 OC patients was used for prognostic analysis of HOXD genes.

Diagnostic Value of HOXD Gene Expression in OC.
As shown in Figure 9, the area under curve (AUC) of HOXD1 was 0.890, the AUC of HOXD3 was 0.538, the AUC of HOXD4 was 0.615, the AUC of HOXD8 was 0.700, the AUC of HOXD9 was 0.748, the AUC of HOXD10 was 0.666, the AUC of HOXD11 was 0.575, the AUC of HOXD12 was 0.537, and AUC of HOXD3 was 0.666. The above results suggest that the expression of HOXD1/8/9 showed good classification efficiency (AUC > 0:7) in OC patients and healthy individuals, indicating that HOXD1/8/ 9 can be used as biomarkers for OC.

The Function of Genes Associated with HOXD Genes.
Some proteins were closely related to the HOXDs (Table S3). These results suggested that changes in the expression profile of HOXDs contributed to the development of OC. The results contained 139 biological processes, mainly including positive regulation of stem cell differentiation, apoptotic process involved in development, kidney mesenchyme development, epithelial tube morphogenesis, and urogenital system development ( Figure 10 and Table S4). The 3 enriched molecular functions included DNA-binding transcription repressor activity, RNA polymerase II-specific, DNA-binding transcription activator activity, RNA polymerase II-specific, and heparin binding ( Figure 10 and Table S4). The results contained 2 cell components, which were mainly related to perinuclear endoplasmic reticulum and transcription factor complex ( Figure 10 and Table S4). The analysis of these functions provides further insight into the cellular localization, geometric distribution, and functional classes of the HOXDs. KEGG analysis showed that 9 pathways, including cell cycle, TGF-beta signaling pathway, gastric cancer, chronic myeloid leukemia, bladder cancer, cellular senescence, Hippo signaling pathway, hepatitis C, and hepatocellular carcinoma, in OC were associated with HOXDs ( Figure 11 and Table S4). These results contributed to the study of the mechanism of action of HOXDs in the development of OC and the possibilities for clinically targeted therapy.

Discussion
HOXDs play an important role in the development, metastasis, and prognosis of various tumors, but the mechanisms are complex. This study used bioinformatics tools to investigate the relationship between HOXDs and the development and prognosis of OC. The results suggested that members of HOXDs could be used as new therapeutic targets and predictive markers for OC. HOXD dysregulation has been reported in many cancers.
HOXD4 protein expression may be associated with poorer prognosis in ovarian serous carcinoma [32]. miR-5692a has oncogenic effects in HCC by targeting HOXD8, which may shed new light on new therapeutic targets and biomarkers for HCC [33]. HOXD1 was lowly expressed in KIRC and correlates with patient OS, DFS, and advanced tumor stage [9]. HOXD9 is upregulated in cervical cancer species, is strongly associated with metastasis rate and poor prognosis in cervical cancer patients, and stimulates the migration and invasive ability of cervical cancer cells by positively regulating HMCN1 levels [34]. Reduced HOXD10 expression promotes a proliferative and aggressive phenotype of prostate cancer [35]. The miR-224/HOXD10 axis may be useful as a promising biomarker and therapeutic approach for the control of NSCLC cell metastasis [36]. HOXD11 may be used as a candidate biomarker for the clinical application of targeted drugs and prognostic assessment therapy for glioma [37]. Progesterone receptor positive cancer tissues have higher levels of HOXD12 and D13 than negative cancer tissues in BRCA [38]. Downregulation of HOXD13 may be a potentially useful prognostic marker for BCRA patients [39]. In this study, the mRNA expression    13 Computational and Mathematical Methods in Medicine levels of HOXD1/12/13 in OC tissues were significantly higher than that in normal ovarian tissues, and the mRNA expression levels of HOXD3/4/8/9/10/11 in OC tissues were significantly lower than that in normal ovarian tissues. The expression of HOXD genes was associated with FIGO stage, primary therapy outcome, tumor status, anatomic neoplasm subdivision, and age. HOXD1/3/4/8/9/10 was negatively correlated with the clinical stage of OC. ROC analysis results suggested that HOXD1/8/9 could be served as ideal biomarkers to distinguish OC from normal tissue. The HOXD9 low expression was associated with OS/PFS shortening.
The lncRNA insulin-like growth factor 2 antisense RNA (IGF2-AS) is predicted to exert a tumor suppressive effect by HOXD1 [40]. HOXD3 plays a key role in BRCA stemness and drug resistance through integrin β3-mediated Wnt/βcatenin signaling [41]. HOXD3 promotes the growth of colorectal cancer (CRC) cells and plays a key role in the development and survival of malignant human CRC cells [42]. miRNA-10a inhibits the expression of HOXD4 in human BRCA cells [43]. HOXD8 upregulates caspases 6 and 7 and cleaves PARP, thereby inducing apoptotic events in CRC cells [44]. HOXD9-RUFY3 axis was associated with the development and progression of GC [45]. HOXD9 promotes a malignant biological process in GC, which could be a potential therapeutic target for GC [46]. HOXD10 was inhibited in colon adenocarcinoma cells, thereby downregulating the RHOC/AKT/MAPK pathway to enhance apoptosis and restrain proliferation, migration, and invasion of colon cancer cells [47]. Downregulation of HOXD10 expression by miR-10b overexpression may induce an increase in prometastatic gene products, such as MMP14 and RHOC, and contribute to the acquisition of a metastatic phenotype by epithelial ovarian cancer cells [48]. POU2F1 activity regulates HOXD10 and HOXD11 to promote proliferative and invasive phenotypes in head and neck cancer [49]. GALNT10 can regulate the proliferation and migratory capacity of GC cells by enhancing the expression of HOXD13 and decreasing the sensitivity to 5-Fu [50]. miR-7156-3p regulates stemness, invasion, and growth of glioma cells by mediating HOXD13 [51]. In this study, KEGG analysis showed that HOXDs were related to pathways including cell cycle, TGF-beta signaling pathway, chronic myeloid leukemia, Hippo signaling pathway, HTLV-I infection, microRNAs in cancer, and signaling pathways regulating pluripotency of stem cells in OC.
Immune infiltration and antitumor immune evasion are key mechanisms of tumor progression [19]. HOXD13 was negatively associated with Th17 cells. HOXD1/3/4/8/9/10/ 11/12/13 were positively associated with other T cells. The emergence of adaptive T cell-based oncology therapies, such as chimeric antigen receptor T cell therapy, may be a promising paradigm for OC, and a better understanding of HOXDs could improve treatment strategies.
This study integrates information on expression levels, mutations, and immune responses to identify potential biomarkers and alterations of HOXD genes in OC. The results promote the understanding of the complex impact of HOXD genes on OC and may help improve clinical decision making. The present study has some limitations in that no in vitro or in vivo experiments were performed to validate the identified role of HOXD genes in OC, which should be attempted in future studies.

Conclusion
The expression of HOXD genes is associated with clinical characteristics. Downregulation of HOXD9 is an independent factor in the poorer prognosis of OC. HOXDs were key players in mediating OC development and progression through multiple pathways, including regulating immune cells and cell cycle, TGF-beta signaling pathway, cellular senescence, and Hippo signaling pathway. The findings suggested that HOXD9 was a new marker of OC prognosis, while HOXD1/4/8/9/10 may be potential targets for the treatment of OC. The members of the HOXD genes may be the response to immunotherapy for OC.

Data Availability
The data analyzed during the current study are available in the TCGA database with the accession number TCGA-OC (Ovarian Serous Cystadenocarcinoma). The data analyzed during the current study are available in the GEO database with the accession number GSE29450. The data used to support the findings of this study are included within the article.

Ethical Approval
Since the resources of TCGA, ONCOMINE, GEO, and other databases are freely available, the patients involved in the databases have received ethical approval and users can download the relevant data for their studies and publish relevant articles free of charge. Our study is based on open data, so there are no ethical issues or other conflicts of interest.

Consent
Consent is not applicable for this study.

Conflicts of Interest
The authors declare that they have no competing interests.