Computational and Mathematical Methods in Medicine
Volume 2022 (2022), Article ID 3268386, 16 pages
https://doi.org/10.1155/2022/3268386
Research Article

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

Buze Chen,1,2 Cui Gao,3 Haihong Wang,1 Jieyun Sun,1 and Zhengxiang Han4

1Department of Gynecology, The Affiliated Hospital of Xuzhou Medical University, Xuzhou, 221000 Jiangsu, China
2Xuzhou Medical University, Xuzhou, 221000 Jiangsu, China
3Department of Obstetrics, Jinhu County People’s Hospital, Huai’an, 223000 Jiangsu, China
4Department of Oncology, The Affiliated Hospital of Xuzhou Medical University, Xuzhou, 221000 Jiangsu, China

Correspondence should be addressed to Buze Chen and Zhengxiang Han

Received 15 May 2022; Revised 29 August 2022; Accepted 2 September 2022; Published 29 September 2022

Academic Editor: Hongxun Tao

Copyright © 2022 Buze Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

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; ) and progression-free survival (PFS) (HR: 0.69; 95% CI: 0.54–0.87; ). 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.

1. 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 summary, 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.

In this study, we used a comprehensive bioinformatics analysis to analyze the potential of HOXDs in OC as a predictor of OC prognosis, possible regulatory mechanisms, and relationship with immune infiltration. We hope that our study will be useful for the prognosis of biomarker and treatment of OC.

2. Materials and Methods

2.1. cBioPortal Analysis

The cBio Cancer Genomics Portal (cBioPortal) (http://cbioportal.org) was applied to study mutations in HOXD genes in OC [16]. Queries for visualization and analysis were performed by entering (1) cancer type: ovarian cancer; (2) 3 selected studies: ovarian serous cystadenocarcinoma (TCGA, Nature 2011), ovarian serous cystadenocarcinoma (TCGA, PanCancer Atlas), and ovarian serous cystadenocarcinoma (TCGA, Firehose Legacy); (3) molecular profile: mutations, structural variants, and copy number alterations; (4) selection of patients/case sets: all samples (1365); and (5) input genes: HOXD1 (ENSG00000128645), HOXD3 (ENSG00000128652), HOXD4 (ENSG00000170166), HOXD8 (ENSG00000175879), HOXD9 (ENSG00000128709), HOXD10 (ENSG00000128710), HOXD11 (ENSG00000128713), HOXD12 (ENSG00000170178), and HOXD13 (ENSG00000128714). After submission of queries, accessions were made including origin studies, mutation profiles, mutation number, overall survival (OS) status, OS (months), disease-free status, and disease-free period (months) tracks.

2.2. Differential Expression of HOXDs

ONCOMINE (https://www.oncomine.org/resource/login.html) was used to analyze the levels of HOXD mRNAs in OC tissues and normal tissues [17]. Screening criteria are as follows: , , and top 10% of gene rank [18].

The analysis was carried out according to the reference [19, 20]. Software: R (version 3.6.3). R package: mainly ggplot2 (for visualization). Molecules: HOXD1/3/4/8/9/10/11/12/13. Data: UCSC XENA (https://xenabrowser.net/datapages/) RNAseq data in TPM (transcripts per million reads) format for TCGA and GTEx processed uniformly by the Toil process [21]. Extracted TCGA (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) OC and corresponding normal tissue data in GTEx. Data filtering: none. Data transformation: RNAseq data in TPM format and log2 transformed for sample-to-sample expression comparisons. Significance markers: ns, ; , ; and , .

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

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

2.5. 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 molecular expression. Additional data: prognostic data from the reference [24].

The survival curves of HOXD12 were plotted using the online Kaplan-Meier plotter database [25].

2.6. Univariate and Multivariate Cox Regression Analysis

Software: R (version 3.6.3). R package: survivor package (version 3.2-10). Statistical methods: Cox regression module. Prognosis type: OS and PFS. Included variables: HOXD1/3/4/8/9/10/11/12/13. Data: RNAseq data in level 3 HTSeq-FPKM format from TCGA OC project. Data conversion: 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.

2.7. ROC Curve Analysis

The analysis was carried out according to the reference [18]. Software: R (version 3.6.3). R packages: mainly the pROC package (for analysis) || ggplot2 package. Molecules: HOXD1/3/4/8/9/10/11/12/13. Clinical variables: tumor vs. normal. Data: UCSC XENA RNAseq data in TPM format for TCGA and GTEx processed uniformly by the Toil process [21]. Extracted OC for TCGA and corresponding normal tissue data in GTEx. Data filtering: none. Data transformation: RNAseq data in TPM format and log2 transformed for analysis.

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

2.9. GO and KEGG Analyses

DAVID database was used to do Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses for the coexpression genes of HOXDs, including BP (biological process), MF (molecular function), CC (cellular component), and pathway analysis [23, 27].

2.10. Correlation between the Expression of HOXD Genes in OC and Immune Cells

The analysis was carried out according to the reference [20]. Software: R (version 3.6.3). R package: GSVA package (version 1.34.0) [28]. Immunoinfiltration algorithm: ssGSEA (built-in algorithm of the GSVA package). Molecules: HOXD1/3/4/8/9/10/11/12/13. Immune cells: aDC (activated DC); B cells; CD8 T cells; cytotoxic cells; DC; eosinophils; iDC (immature DC); macrophages; mast cells; neutrophils; NK CD56bright cells; NK CD56dim cells; NK cells; pDC (plasmacytoid DC); T cells; T helper cells; Tcm (T central memory); Tem (T effector memory); Tfh (T follicular helper); Tgd (T gamma delta); Th1 cells; Th17 cells; Th2 cells; and Treg. Data: RNAseq data in level 3 HTSeq-FPKM format from TCGA OC project. Data conversion: RNAseq data in FPKM format were converted to TPM format and log2 transformed. Data filtering: control/normal removed (not all items have control/normal). Other data: markers for 24 immune cells were obtained from the reference [29].

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. were considered statistically significant.

3. Results

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

Figure 1: The genetic alteration of HOXD genes in OC by cBioPortal.
Figure 2: Percentage of HOXD genes in OC cases calculated using the cancer type summary in cBioPortal.

As shown in Figure 3 and Table 1, the expression of HOXD4/10/11 mRNA in OC tissues was significantly lower than that in normal ovarian tissues (). Among them, HOXD10 had the highest expression change (, ), and 2 data sets confirmed this [30, 31]. As shown in Figure 4, the HOXD1 expression in OC tissues was significantly higher than that in normal ovarian tissues ( vs. , ), the HOXD4 expression in OC tissues was significantly lower than that in normal ovarian tissues ( vs. , ), the HOXD8 expression in OC tissues was significantly lower than that in normal ovarian tissues ( vs. , ), the HOXD9 expression in OC tissues was significantly lower than that in normal ovarian tissues ( vs. , ), the HOXD10 expression in OC tissues was significantly lower than that in normal ovarian tissues ( vs. , ), the HOXD11 expression in OC tissues was significantly lower than that in normal ovarian tissues ( vs. , ), the HOXD12 expression in OC tissues was significantly higher than that in normal ovarian tissues ( vs. , ), the HOXD13 expression in OC tissues was significantly higher than that in normal ovarian tissues ( vs. , ), and there was no significant difference in HOXD3 between OC tissues and normal ovarian tissues ( vs. , ). 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 (, ), HOXD4 was significantly lower expressed in OC tumor tissues (, ), HOXD8 was significantly lower expressed in OC tumor tissues (, ), HOXD9 was significantly lower expressed in OC tumor tissues (, ), and HOXD12 was significantly higher expressed in OC tumor tissues (, ). 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.

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.
Table 1: Differential expression of HOXD mRNA in OC and ovarian tissues (ONCOMINE database).
Figure 4: The expression of HOXDs in normal ovarian tissue was compared with that of the OC tissues. Significance markers: ns, ; , ; and , .
Table 2: Differential expression of HOXD mRNA in OC and ovarian tissues (GEO29450).

We examined the correlation between HOXD genes using the Pearson correlation analysis. As shown in Figure 5, there was no significant correlation between HOXD12 and HOXD1/3/4. Other HOXD genes were significantly positively correlated with each other.

Figure 5: Correlation between every two genes of HOXD genes.
3.2. 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. The clinical characteristics included FIGO stage, primary therapy outcome, race, age, histologic grade, anatomic neoplasm subdivision, venous invasion, lymphatic invasion, tumor residual, tumor status, and age. As shown in Table S2, high expression of HOXD1 was associated with FIGO stage (), low expression of HOXD3 was associated with FIGO stage () and histological grade (), low expression of HOXD4 was associated with FIGO stage (), low expression of HOXD8 was associated with FIGO stage (), low expression of HOXD9 was associated with FIGO stage (), primary therapy outcome (), and tumor status (), low expression of HOXD10 was associated with anatomic neoplasm subdivision (), low expression of HOXD11 was associated with age (), and high expression of HOXD13 was associated with FIGO stage ().

Figure 6: Expression of HOXDs in OC patients at different tumor stages (GEPIA). (a) HOXD1, (b) HOXD3, (c) HOXD4, (d) HOXD8, (e) HOXD9, (f) HOXD10, (g) HOXD11, (h) HOXD12, and (i) HOXD13 were analyzed in this study.

As shown in Figure 7 and Figure S1A, low HOXD9 expression was associated with shorter overall survival (OS) (HR: 0.75; 95% CI: 0.58–0.98; ) and progression-free survival (PFS) (HR: 0.69; 95% CI: 0.54–0.87; ). Low HOXD12 expression was associated with shorter OS (HR: 0.86; 95% CI: 0.75–1; ) and PFS (HR: 0.78; 95% CI: 0.68–0.89; ). The OS of the HOXD9/12 mRNA lower expression group was lower than that of the HOXD9/12 low expression group at all time points. It suggested that HOXD9 were risk factors of OC. As shown in Figure 8 and Figure S1B, HOXD9 and HOXD12 low expressions were associated with PFS shortening. The mRNA of HOXD9 can be used as indicators for predicting OC/PFS progression. As shown in Table 3, HOXD9 (HR: 0.754; 95% CI: 0.581-0.978; ) was independently correlated with OS in multivariate analysis. As shown in Table 4, HOXD9 (HR: 0.690; 95% CI: 0.544-0.875; ) was independently correlated with OS in multivariate analysis.

Figure 7: The expression of HOXDs is associated with poor OS in patients with OC. (a) HOXD1, (b) HOXD3, (c) HOXD4, (d) HOXD8, (e) HOXD9, (f) HOXD10, (g) HOXD11, and (h) HOXD13 were analyzed in this study.
Figure 8: The expression of HOXDs is associated with poor PFS in patients with OC. (a) HOXD1, (b) HOXD3, (c) HOXD4, (d) HOXD8, (e) HOXD9, (f) HOXD10, (g) HOXD11, and (h) HOXD13 were analyzed in this study.
Table 3: Univariate and multivariate Cox regression analyses of HOXDs associated with OS.
Table 4: Univariate and multivariate Cox regression analyses of HOXDs associated with PFS.
3.3. 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 () in OC patients and healthy individuals, indicating that HOXD1/8/9 can be used as biomarkers for OC.

Figure 9: ROC curve showed the efficiency of HOXD expression level in distinguishing OC tissue from nontumor tissues.
3.4. 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.

Figure 10: GO analysis of HOXD coexpression genes predicted by DAVID. BP: biological process; MF: molecular function; CC: cellular component.
Figure 11: KEGG analysis of HOXD coexpression genes predicted by DAVID.
3.5. Correlation of HOXD Gene Expression and Immune Cells in OC

As shown in Figure 12, there was a correlation between HOXD gene expression and immune cells in OC. HOXD1 gene expression was positively correlated with some TIICs, including aDC, cytotoxic cells, DC, iDC, macrophages, neutrophils, NK CD56bright cells, NK CD56dim cells, T cells, T helper cells, Tcm, Tem, TFH, Th1 cells, Th17 cells, and TReg, and negatively correlated with NK cells. HOXD3 gene expression was positively correlated with some TIICs, including aDC, DC, NK CD56bright cells, NK CD56dim cells, T helper cells, Tcm, Tem, and Th1 cells. HOXD4 gene expression was positively correlated with some TIICs, including aDC, cytotoxic cells, DC, NK CD56bright cells, NK CD56dim cells, pDC, T helper cells, Tcm, Tem, Th1 cells, and TReg. HOXD8 gene expression was positively correlated with some TIICs, including aDC, DC, NK CD56dim cells, and Tem. HOXD9 gene expression was positively correlated with some TIICs, including aDC, DC, NK CD56dim cells, Tem, TFH, Th2 cells, and TReg. HOXD10 gene expression was positively correlated with some TIICs, including DC, iDC, macrophages, neutrophils, T helper cells, Tem, and TFH. HOXD11 gene expression was positively correlated with some TIICs, including B cells, CD8 T cells, DC, eosinophils, iDC, macrophages, mast cells, neutrophils, NK CD56dim cells, T cells, T helper cells, Tem, TFH, Tgd, Th1 cells, Th2 cells, and TReg. HOXD12 gene expression was positively correlated with some TIICs, including macrophages, T helper cells, Th2 cells, and TReg. HOXD13 gene expression was positively correlated with some TIICs, including iDC, macrophages, neutrophils, NK CD56dim cells, NK cells, T helper cells, Tcm, Tem, TFH, Tgd, Th1 cells, Th2 cells, and TReg, and negatively correlated with Th17 cells.

Figure 12: Correlation between the expression of each HOXD gene and the 24 TIICs of OC (lollipop plot). (a) HOXD1, (b) HOXD3, (c) HOXD4, (d) HOXD8, (e) HOXD9, (f) HOXD10, (g) HOXD11, (h) HOXD12, and (i) HOXD13. In the color bar, the darker the color, the smaller the value, indicating a higher statistical significance. The bubble size represents the correlation value; the larger the bubble, the larger the correlation value.

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

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

Authors’ Contributions

BC and ZH designed the research. BC, CG, and HW collected the data. BC, CG, HW, JS, and ZH did the data analysis. BC, CG, and HW wrote the manuscript. BC and ZH reviewed the manuscript. All authors read and approved the final manuscript. Buze Chen, Cui Gao, and Haihong Wang contributed equally to this work.

Acknowledgments

The authors thank TCGA and GEO databases for providing the data. This work was supported by the Six Talent Peaks Project in Jiangsu Province (LGY202006) and Xuzhou Key R&D Programme (ZYSB20210489).

Supplementary Materials

Supplementary 1. Figure S1: the prognostic value of mRNA level of HOXD12 in OC patients (Kaplan-Meier plotter). (A) Overall survival of HOXD12 and (B) progression free interval of HOXD12. The HR and log-rank values are indicated at the right corner of the plot. indicates that the difference is statistically significant.

Supplementary 2. Table S1: patient characteristics in TCGA for OC.

Supplementary 3. Table S2: correlation of HOXD genes with clinical characteristics of OC patients.

Supplementary 4. Table S3: the list of HOXD coexpression genes.

Supplementary 5. Table S4: GO and KEGG analyses of HOXD coexpression genes.