Diagnostic and prognostic value of FOXD1 expression in head and neck squamous cell carcinoma

FOXD1 has been reported to function as an oncogene in several types of cancer. This study evaluated the expression of FOXD1 and its role in head and neck squamous cell carcinoma (HNSCC). We mined the Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases for expression profiles, clinical significance, and potential mechanisms of FOXD1in HNSCC. Our validation cohort consisted of FOXD1 mRNA expression in 162 paired HNSCC and adjacent normal tissues, as determined using quantitative real-time polymerase chain reaction. FOXD1 expression was upregulated in HNSCC in the public databases and in the validation cohort. The expression level of FOXD1 was associated with DNA amplification and methylation level. The areas under the curves (AUC) of TCGA cohort and the validation cohort were 0.855 and 0.843, respectively. Furthermore, higher FOXD1 expression was significantly associated with worse overall survival (hazard ratio [HR]: 1.849, 95% confidence interval [CI]: 1.280-2.670, P = 0.001) and a lower rate of recurrence-free survival (HR: 1.650, 95% CI: 1.058-2.575, P = 0.027) in patients with HNSCC. Moreover, gene set enrichment analysis showed that cases of HNSCC with FOXD1 overexpression were enriched in bladder cancer, cell cycle, DNA replication, glycosaminoglycan biosynthesis chondroitin sulfate, homologous recombination, glycan biosynthesis, nucleotide excision repair, p53 signaling pathway, pyrimidine metabolism, and spliceosome pathways. In summary, FOXD1 was significantly upregulated in HNSCC and was a good diagnostic biomarker and an independent predictor of poor survival and low rate of recurrence-free survival in patients with HNSCC.


Introduction
Head and neck cancer is the sixth most common cancer globally, and is comprised of malignant tumors in the oral cavity, oropharynx, hypopharynx, and larynx. Head and neck squamous cell carcinoma (HNSCC) is the predominant histological type (>90%) of head and neck cancer [1,2]. The International Agency for Research on Cancer reported that the global incidence of HNSCC is more than one million new cases annually, and HNSCC results in 543,000 deaths per year [3]. Many studies have shown that tobacco exposure and alcohol consumption are important risk factors for the development of HNSCC [4,5]. Recently, human papillomavirus (HPV) has also been shown to be a strong and independent risk factor for the development of HNSCC [6]. Management of HNSCC requires a multi-faceted approach that includes surgery, radiation, and/or chemotherapy. Although substantial progress has been made in treatment of HNSCC, outcomes remain poor, especially in patients with advanced disease. The 5-year survival rate is approximately 50%, which represents only a slight improvement over the last two decades [7]. This poor survival rate may be due to late diagnosis, low therapeutic response rates, and high rates of recurrence and metastasis [8]. Furthermore, reliable and specific biomarkers for Ivyspring International Publisher diagnosis and prognosis of HNSCC are lacking. Therefore, it is critical to determine the molecular correlates of HNSCC and to identify reliable biomarkers for diagnosis, prognosis, and monitoring of recurrence, which would allow for improved personalized treatment strategies for patients with HNSCC.
Forkhead box D1 (FOXD1), located on chromosome 5q12, belongs to the forkhead box transcription factor family, and is involved in numerous physiological processes and biological functions, such as embryonic development and organogenesis, cell cycle regulation, control of metabolism, stem cell niche maintenance, and signal transduction [9,10]. FOXD1 was identified and described for the first time in the forebrain neuroepithelium and is considered an important factor during retinal development [11]. Furthermore, FOXD1 was found to be a mediator of successful progression of cell reprogramming through selfrenewal and differentiation [12,13]. Recent studies showed that FOXD1 was associated with carcinogenesis, tumor progression, and metastasis in numerous cancers [14][15][16]. FOXD1 was reported to be overexpressed in colorectal cancer tissues, and expression levels correlated with tumor size, differentiation, tumor node and metastasis (TNM) stage, lymph node metastasis, and poor prognosis [17]. FOXD1 has also been shown to be highly expressed in non-small cell lung cancer, and promoted cell proliferation, migration, invasion, and metastasis through activation of Vimentin [18]. Zhao et al. found that FOXD1 was up-regulated in breast cancer tissues and promoted cell proliferation and chemotherapeutic drug resistance by targeting p27 expression [19]. These findings suggested that FOXD1 may function as an oncogene in several cancers. However, the role of FOXD1 in HNSCC has not been characterized.
In the present study, we evaluated FOXD1 mRNA expression in HNSCC using data from The Cancer Genome Atlas (TCGA) database and the Gene Expression Omnibus (GEO) database. We also investigated the association of FOXD1 expression with clinicopathological parameters and evaluated its potential as a diagnostic biomarker of HNSCC. Furthermore, we assessed the prognostic value of FOXD1 for overall survival (OS) and recurrence-free survival (RFS). We performed gene set enrichment analysis (GSEA) to identify FOXD1 related signaling pathways involved in tumorigenesis and progression of HNSCC. In addition, tissue samples from patients for whom clinicopathological and survival data were available were analyzed to validate the results of the bioinformatic analysis.

Mining analysis using the Cancer Genome Atlas (TCGA) dataset
Level-3 FOXD1 RNA-seq data consisting of 520 HNSCC tissues and 44 normal controls were downloaded from the University of Santa Cruz Xena browser (https://xenabrowser.net, up to June 11, 2019). Related clinicopathological data including sample type, age at initial pathologic diagnosis, gender, histological types, smoking history, alcohol history, anatomic neoplasm subdivision, HPV status by p16 testing, perineural invasion present, histologic grade pathologic T, pathologic N, pathologic stage, OS status, OS time, RFS status, RFS time, DNA methylation, and gistic2 threshold-processed copynumber alteration (-1: copy deletion; 0: no change; +1: amplification) were obtained for secondary analysis.
The β values of cg23454038 probes mapping 200 bp downstream of the transcription start sites of FOXD1 were defined as the FOXD1 promoter methylation level.

Mining analysis using the Gene Expression Omnibus (GEO) dataset
The keywords "head and neck squamous cell carcinoma" were used to search the GEO database (https://www.ncbi.nlm.nih.gov/gds/), and FOXD1 expression profiles from GSE6631 [19] were downloaded. The platform for GSE6631 was GPL8300, [HG_U95Av2] Affymetrix Human Genome U95 Version 2 Array, which contains twenty-two paired head and neck squamous cell carcinoma samples and adjacent normal tissue.

Clinical sample collection
One hundred sixty-two surgical specimens from patients with HNSCC, including tumor and adjacent normal tissues, were collected (125 male and 37 female, age 38-79 years; average age 61.7-years of age). The patients were diagnosed based on clinical features and histopathological examination at the Ningbo Medical Centre Lihuili Hospital (Ningbo, China) and the Affiliated Tumor Hospital of Xiangya Medical School (Changsha, China) from February 2014 to November 2018. Upon removal, specimens were immediately placed in RNA-fixer Reagent (Bioteke, Beijing, China) and stored at -80 °C until be used in experiments. None of the patients underwent radiation or chemotherapy prior to surgery. Histological type was identified independently by two experienced pathologists. Clinicopathological features were collected from medical records, and tumor stages were classified according to the 8th Edition HNSCC TNM staging system of the American Joint Committee on Cancer (AJCC) [20]. During the follow-up period, 11 patients were censored, and 72 patients died. The median patient follow-up time was 27.2 months. Overall survival time was defined as the period between pathological diagnosis and death. This study was approved by the Ethics Committee of the Ningbo Medical Centre Lihuili Hospital and the Affiliated Tumor Hospital of Xiangya Medical School. All patients provided written informed consent.

Survival analysis
Patients with HNSCC for whom FOXD1 expression and survival data were collected were divided into 2 groups (low and high FOXD1 expression) based on the maximum Youden index of the receiver operating characteristic (ROC) curves for death and recurrence. Overall survival and RFS were compared between the high and low FOXD1 expression groups using Kaplan-Meier analysis with the log-rank test. Univariate and multivariate Cox proportional hazards models were performed to evaluate the relative risk factors associated with OS or RFS, and hazard ratios (HR) with 95% confidence intervals (CI) were obtained for each variable. Only significant factors in the univariate analysis were included in the multivariate analysis.

Gene set enrichment analysis (GSEA) resulted in identification of FOXD1-related signaling pathways in HNSCC
Gene set enrichment analysis was performed to identify potential mechanisms of FOXD1 in development of HNSCC. Samples from TCGA were divided into high and low FOXD1 expression groups based on the median FOXD1 expression. MsigDB Collection (h.all.v6.2.symbols.gmt) was used as a reference for GSEA. Significantly enriched pathways related to tumor biological process were selected according to normalized enrichment score (NES) with 1000 permutations. A pathway was regarded as significantly enriched when the P-value was less than 0.05.

Statistical analysis
Statistical Program for Social Sciences (SPSS) 20.0 software (SPSS Inc., Chicago, IL, USA) and R 3.1.2 software (https://www.r-project.org/) were used to perform statistical analysis and to generate figures. Between-group comparisons of FOXD1 expression and the correlation between FOXD1 expression and clinicopathological features were performed using independent or paired Student's t-tests, as appropriate. Receiver operating characteristic curves were generated, and areas under the curves (AUC) were calculated to determine the diagnostic power of FOXD1 in HNSCC. The cut-off point was defined as the maximum Youden index. Pearson correlation coefficients were generated to evaluate the association between FOXD1 mRNA expression and FOXD1 DNA methylation levels. P ≤ 0.05 was considered statistically significant.

FOXD1 was upregulated in patients with HNSCC
Using Xena browser, we reviewed FOXD1 mRNA expression in 520 HNSCC and 44 normal tissues in TCGA database. The results indicated that FOXD1 mRNA expression was significantly higher in HNSCC tissues than that in normal tissues (P = 9.64E-22, Figure 1A and B). Similarly, FOXD1 was significantly upregulated in one dataset (containing 22 pairs of HNSCC and adjacent normal tissues) obtained from the GEO database (P = 8.04E-4, Figure  1C). To further investigate the expression of FOXD1 in HNSCC, we collected HNSCC tissue and adjacent normal tissue from 162 patients with HNSCC for use as a validation cohort. Quantitative RT-PCR analysis showed that FOXD1 expression was significantly higher in HNSCC tissues than that in normal tissues (P = 2.26E-27, Figure 1D).

Relationship between FOXD1 expression levels and clinicopathological factors
We examined the correlation between FOXD1 expression and clinicopathological features of patients with HNSCC. In the TCGA cohort, we found that patients with HNSCC who did not smoke, and had tumors located in the larynx, had significantly elevated FOXD1 expression. However, FOXD1 expression was not associated with age, gender, alcohol history, histologic grade, HPV status, perineural invasion, tumor category, nodal category, or pathologic stage (Table 1). Furthermore, FOXD1 expression was not significantly associated with any clinicopathological factors in the validation cohort ( Table 2).

The diagnostic value of FOXD1 for HNSCC
A receiver operating characteristic (ROC) curve was generated and the area under the ROC curve (AUC) was calculated to evaluate the diagnostic value of FOXD1. The AUC of TCGA cohort was 0.855, with sensitivity and specificity values of 0.746 and 0.864, respectively. In the validation cohort, the AUC was 0.843, with sensitivity and specificity values of 0.833 and 0.722, respectively (Figure 2).

High FOXD1 expression was an independent risk factor for poor OS in patients with HNSCC
Five hundred seventeen patients with complete FOXD1 expression and overall survival data in TCGA cohort were divided into high expression (289 patients) and low expression groups (228 patients) according to the maximum Youden index of the ROC curve for death. The Kaplan-Meier curve showed that high mRNA expression of FOXD1 was associated with poor overall survival rate ( Figure 3A, log-rank P = 2.51E-4). To further assess the prognostic value of FOXD1, we collected follow-up information from 162 patients with HNSCC. These results showed that patients with high FOXD1 expression had worse overall survival than patients with lower FOXD1 expression ( Figure 3B, log-rank P = 0.0098). Subsequently, a Cox proportional hazard regression model was used to screen prognostic factors for patients with HNSCC in TCGA database. As shown in the

High FOXD1 expression was an independent factor for poor RFS in patients with HNSCC patients
Using TCGA dataset, which contained follow-up data for HNSCC recurrence in 438 patients, we evaluated the association of FOXD1 expression with RFS in patients with HNSCC. The Kaplan-Meier curve for RFS showed that high FOXD1 expression was associated with poor RFS (Figure 4, log-rank P = 0.007). This finding was further confirmed by univariate Cox proportional hazard analysis (Table 4), in which patients with high FOXD1 expression in tumors were at significantly increased risk of recurrence compared to patients with low FOXD1 expression in tumors (HR = 1.725, 95% CI: 1.152-2.583, P = 0.008). Subsequently, multivariate analysis showed that high FOXD1 expression was an independent risk factor for RFS in patients with HNSCC (HR = 1.650, 95% CI: 1.058-2.575, P = 0.008) after adjustment for significant prognostic clinicpathological parameters (smoking history and pathologic stage).

FOXD1 expression was related to DNA copy number alteration and promoter methylation in HNSCC
To further investigate the mechanism of FOXD1 overexpression in HNSCC, we examined the association of FOXD1 expression with copy number alterations and promoter methylation using TCGA database. Among 514 patients for whom FOXD1 DNA copy number was determined, 33 cases (14.5%) showed DNA amplification (+1) and 216 cases (31.8%) showed copy deletion (-1). We showed that DNA amplification was associated with increased FOXD1 expression (P = 0.007, compared to the copy deletion group, Figure 5A). Subsequently, we observed a negative correlation (Person r = -0.295, P = 7.27E-12) between FOXD1 expression and promoter methylation ( Figure 5B).

GSEA identified FOXD1-related signaling pathways in HNSCC
We compared the low and high FOXD1 groups using GSEA analysis to identify FOXD1-related signaling pathways activated in HNSCC. The results showed significant differences in gene sets (false discovery rate [FDR] P-value < 0.05) in the enrichment of MSigDB Collection (h.all.v6.2.symbols.gmt). Detailed results are provided in Table 5. The results showed that bladder cancer, cell cycle, DNA replication, glycosaminoglycan biosynthesis chondroitin sulfate, homologous recombination, glycan biosynthesis, nucleotide excision repair, p53 signaling pathway, pyrimidine metabolism, and spliceosome were enriched in the high FOXD1 expression group (Figure 6).

Discussion
Recently, FOXD1 has been reported to be highly expressed in several cancers, such as colorectal cancer [17], lung cancer [18], breast cancer [19], and Hodgkin's lymphoma [22]. Increased FOXD1 expression has been shown to promote cell proliferation, migration, invasion, and tumorigenesis.
However, the role of FOXD1 in HNSCC has not been well characterized. The Cancer Genome Atlas (TCGA) was a large-scale effort to comprehensively characterize 33 major cancer types, including HNSCC, for which 528 cases were included [23]. Our study using TCGA database showed that FOXD1 mRNA expression was higher in HNSCC tissue than in normal tissue. The HNSCC mRNA microarray dataset GSE6631 was downloaded from the GEO database. The GEO database is a public functional genomics data repository that includes array-and sequencebased gene profiles and next-generation sequencing [24]. Analysis of the GEO dataset showed that FOXD1 was significantly upregulated in HNSCC tissues compared with corresponding adjacent normal tissues. Furthermore, we collected 162 paired HNSCC tissues and adjacent normal tissues to validate our bioinformatics analysis. Quantitative RT-PCR analysis confirmed that FOXD1 mRNA was significantly overexpressed in HNSCC tissues. These findings indicated that FOXD1 may function as an oncogene in HNSCC.
HNSCC is a heterogeneous group of tumors located in the upper aerodigestive tract with multifactorial etiologies. Alcohol consumption, betel nut chewing, and human papillomavirus infection are significant risk factors for development of HNSCC in the upper gastrointestinal tract [25][26][27][28]. In contrast, tobacco smoking is a significant risk factor for development of HNSCC in the upper respiratory tract [29,30]. In the present study, we found that tumors located in the upper respiratory tract (larynx) showed significantly higher FOXD1 expression levels than those in the upper gastrointestinal tract (oral cavity, oropharynx and hypopharynx), which suggested that laryngeal squamous cell carcinoma may be more susceptible to FOXD1 overexpression.    Early diagnosis is critical to successful management of cancer [31,32]. More than half of HNSCC patients are at advanced stages at the time of diagnosis because of concealment of the anatomical site and lack of specific and reliable indicators. Furthermore, the age of diagnosis is slowly decreasing [33]. Therefore, identification of biomarkers for early diagnosis of HNSCC is of great importance. In the current study, ROC curves were generated, and AUC values were calculated to evaluate FOXD1 expression as a potential diagnostic marker of HNSCC. The AUC values of TCGA cohort and the validation cohort were 0.855 and 0.843, respectively. Compared to traditional tumor markers, such as CEA (Carcino Embryonic Antigen), SCC (Squamous Cell Carcinoma Antigen), TPS (Tissue Polypeptide Specific Antigen), and CYFRA 21-1 [34,35], FOXD1 had a greater ability to discriminate patients with HNSCC from healthy individuals, which indicated that FOXD1 could serve as a potential early diagnostic biomarker for HNSCC, especially when combined with other efficient markers.
Molecular characteristics, pathogenesis, and prognosis of head and neck cancers are heterogeneous. Several studies have shown that HPVpositive patients with HNSCC showed significantly improved overall and disease-free survival compared with HPV-negative patients [6,36]. Development of new technologies, such as microarray technology and next-generation sequencing, have allowed for collection of large amounts of data for molecular subclassification of cancer, which has allowed for development of individualized treatment programs [37,38]. Recent studies showed that overexpression of FOXD1 were associated with poor OS in colorectal cancer [17], non-small cell lung cancer [18], and breast cancer [19]. Consistent with these studies, the log-rank test performed in our study showed that high FOXD1 expression was associated with significantly worse OS in both TCGA cohort and the validation cohort. In addition, multivariate proportional hazard Cox regression analysis showed that high FOXD1 expression could serve as an independent indicator of poor prognosis in patients with HNSCC, which suggested that FOXD1 may be a promising prognostic biomarker and therapeutic target for HNSCC. Recurrence rate is believed to be an important contributor to the poor prognosis associated with HNSCC [39]. Approximately 30-40% of patients with HNSCC suffer from recurrence or metastasis following treatment [40]. In our study, multivariate Cox regression analysis showed that FOXD1 was an independent predictor of recurrence. This result suggested that monitoring FOXD1 expression may improve outcomes in patients with HNSCC.
Cancer results from accumulation of genetic and epigenetic modifications of oncogenes and tumor-suppressor genes, resulting in metabolic dysfunction and uncontrolled proliferation [41,42]. Amplification of DNA is the major genetic change that results in cancer-specific expression of critical genes [43,44]. In addition, DNA methylation is an important epigenetic modification involved in the inactivation of numerous tumor suppressor genes [45,46]. Therefore, we explored the association of FOXD1 expression with copy number alterations and promoter methylation levels to determine the mechanism of FOXD1 overexpression in HNSCC. The results showed that DNA amplification was associated with elevated FOXD1 RNA expression. In addition, we also found that FOXD1 RNA expression was negatively correlated with promoter methylation level. These findings indicated that both genetic and epigenetic alterations contributed to dysregulation of FOXD1 in HNSCC. Gene set enrichment analysis showed that bladder cancer, cell cycle, DNA replication, glycosaminoglycan biosynthesis chondroitin sulfate, homologous recombination, glycan biosynthesis, nucleotide excision repair, p53 signaling pathway, pyrimidine metabolism, and spliceosome may be key pathways regulated by FOXD1 in HNSCC. These findings should be further validated by rigorous in vitro and in vivo experiments.

Conclusion
In summary, our analysis showed that FOXD1 expression was significantly elevated in HNSCC tissues relative to normal tissues in TCGA database, the GEO database, and the validation cohort. Genetic and epigenetic alterations contributed to upregulation of FOXD1 in HNSCC. In addition, elevated FOXD1 expression was a good diagnostic biomarker and independent predictor of poor OS and RFS in patients with HNSCC. Furthermore, FOXD1 overexpression was significantly associated with bladder cancer, cell cycle, DNA replication, nucleotide excision repair, and p53 signaling pathways. Future studies are needed to characterize the specific role of FOXD1 in HNSCC.