FOXD1 expression in head and neck squamous carcinoma: a study based on TCGA, GEO and meta-analysis

Head and neck cancer (HNC) is an aggressive malignant tumor arising in oral cavity, oropharynx, hypopharynx and larynx, representing the sixth most common cancer and fifth leading cause of cancer-related deaths worldwide [1]. It is estimated that more than 51000 new cases of HNC were diagnosed in 2018, resulting in an annual incidence of 10000 deaths in United States [2]. The vast majority of HNC are squamous carcinomas (head and neck squamous cancer, HNSC), accounting for more than 90% of HNC. Over the past 30 years, even if the combination of surgery, chemotherapy, radiotherapy and novel immune checkpoint inhibitors have been applied, the 5-year survival rate for HNSC, especially for advanced HNSC patients, has not yet been remarkably improved. Growing evidence has demonstrated that HNSC is a highly complex and heterogeneous disease, involving distinct histological types, different anatomical sites and varied genetic and molecular alterations, thus, hampering the clinicians to accurately diagnose, guide individualized treatment, and improve the prognosis of HNSC patients. Therefore, there is still an urgent need to explore the effective biomarkers and underlie the molecular mechanism in the development and progression of HNSC.

Forkhead box D1 (FOXD1), also known as FKHL8, FREAC-4 and FREAC4, is a new member of FOX transcription factor family and located on 5q13.2. FOXD1 has demonstrated multilevel roles during normal development, adult physiology, and several diseases’ pathogenesis. Previous studies reported that FOXD1 was a mediator and indicator in the process of cell programming, especially in kidney development [3–5]. Herrera et al. reported that foxd1 plays a dual role in retinal ganglion cell axon migration through the optic chiasm [6]. Moreover, accumulating evidence has demonstrated that FOXD1 is implicated in the carcinogenesis, including lung cancer [7], colorectal cancer [8], ovarian cancer [9], breast cancer [10], nasopharyngeal carcinoma [11], gastric cancer [12], and prostate cancer [13]. Gao et al. also found that FOXD1 was up-regulated and directly correlated with the glioma grade, and regulated glioblastoma cell behaviors [14]. However, little is known about the role of FOXD1 in the progression of HNSC.

In the present study, we analyzed FOXD1 expression pattern using The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets. First, we investigated the FOXD1 expression between HNSC tissues and normal tissues, analyzed the correlation between FOXD1 expression and clinical characteristics, and evaluated the prognostic value of FOXD1. Moreover, we assessed the relationship between FOXD1 expression and tumor microenvironment (TME) and immune cell infiltration using Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) and Cell-type Identification By Estimating Relative Subsets Of known RNA Transcripts (CIBERSORT) algorithms. Accordingly, the relationship between FOXD1 expression and underlying biological functions and signal pathways in HNSC was analyzed using Gene Set Enrichment Analysis (GSEA) tool. Together, our results could provide a new insight into the potential mechanism contributing to HNSC development and progression, and highlight new targets for clinical treatment.

FOXD1 expressions were compared between cancer tissues and normal tissues across all TCGA tumors using TIMER online tool (http://timer.cistrome.org/). In addition, we analyzed the significant FOXD1-correlated genes using LinkedOmics online database (http://linkedomics.org). The LinkedOmics database contains multiomics data for 32 cancer types and a total of 11158 patients from TCGA project. Then, we evaluated the biological functions and pathways for FOXD1 and FOXD1-correlated genes through STRING database (https://string-db.org/).

The RNA-sequencing data and related clinical information were downloaded from TCGA database (https://www.cancer.gov/tcga). The sample inclusion criteria: completed FOXD1 sequencing data and detailed clinical information, including age, gender, TNM stage, follow-up information.

The microarray HNSC datasets were searched from GEO database (http://www.ncbi.nlm.nih.gov/geo/) up to December 2020. The search terms included the following keywords: (head OR neck OR oral) AND (tumor OR cancer OR carcinoma OR neoplasm). Each GEO dataset should meet the following criteria: (1) completed FOXD1 expression data; (2) study including HNSC and control group; (3) the number of each group was greater than three. The flow diagram of the GEO dataset selection is shown in Figure 1. All characteristics of GEO datasets are listed in Table 1.

FOXD1 expressions were compared between HNSC tissues and normal tissues in TCGA and GEO datasets. We also compared FOXD1 expression between paired normal and tumor tissues. Moreover, the associations between FOXD1 expression and clinical parameters were evaluated using TCGA dataset, including age, gender, distant metastasis, clinical stage, T stage, N stage and HPV status. A P<0.05 was considered significant.

Kaplan–Meier curves were used to evaluate the correlation between FOXD1 expression and overall survival (OS), tested by Log-Rank test. Univariate and multivariate Cox regression analyses were performed to evaluate the independent prognostic value of FOXD1, as well as age, gender, distant metastasis, clinical stage, N stage, and T stage in HNSC. All independent prognostic parameters identified by multivariate Cox regression analyses were used to construct a nomogram to investigate the probability of 3- and 5-year OS of HNSC. The concordance index (C-index) was calculated to quantify the discrimination performance of this nomogram. The calibration curve was plotted to assess whether the predict probability was in agreement with actual rate in the nomogram.

TME provides a more complex matrix around tumor cells, which contributes to tumor development, progression, metastasis, and drug resistance [15]. We adopted ESTIMATE method to calculate the stromal score and immune score, inferring the tumor purity for TCGA dataset. The stromal score and immune score were compared between high- and low- FOXD1 expression group. Moreover, the CIBERSORT algorithm was performed to analyze the transcript data from TCGA dataset and calculate the proportions of immune infiltrating cell subtypes. Moreover, the correlation between FOXD1 expression and the proportion of immune cell subtype was evaluated by Pearson correlation coefficient.

GSEA was employed to identify the underlying biological functions and pathways of high- and low-FOXD1 expression groups. The reference gene sets ‘c2.cp.kegg.v7.2.symbols.gmt gene sets’ and ‘C5.go.bp.v7.2.symbols.gmt gene set’, were obtained from the Molecular Signatures Database (MSigDB). The nominal P-value <0.05 and normalized enrichment score (NES) < 25% were considered significant.

The human normal oral epithelial cells HOEC were obtained from Keygen Company Co., Ltd. HNSC cell lines CAL-27, SCC-9, and TCA-8113 were purchased from ATCC. HOEC and HNSC cell lines were cultured in 1640 and DMEMs containing 10% fetal bovine serum (Thermo Fisher Scientific, Waltham, MA) at 37°C and 5% CO₂ in an incubator, separately.

Total RNA was extracted by TRIzol reagent (Invitrogen) and cDNA was synthesized using PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa, Japan). The real-time polymerase chain reaction (RT-qPCR) was performed to determine RNA level using SYBR Green PCR Kit (TaKaRa, Japan). The primer sequences for the detection of FOXD1 were 5′-TGAGCACTGAGATGTCCGATG-3′ (forward primer) and 5′-CACCACGTCGATGTCTGTTTC-3′ (reverse primer). Glyceraldehyde‐3‐phosphate dehydrogenase (GAPDH) was amplified to normalize FOXD1 levels. The expression level was calculated by 2^−ΔΔC_t method. Each sample was measured in triplicates.

Cells were lysed with RIPA Buffer (Beyotime Biotechnology Institute). Denatured proteins were separated on 10% SDS/PAGE and transferred on to PVDF membranes. The blocked membranes were incubated with primary FOXD1 antibody (Absin Bioscience, Inc) at 4°C overnight, and then incubated with horseradish peroxidase-conjugated goat anti-rabbit secondary antibody (Cell Signaling Technology, Inc) at room temperature for 2 h. Protein band was visualized by enhanced chemiluminescence reagent (Thermo Fisher Scientific, Inc.).

Primary HNSC tissues were collected from HNSC patients, who underwent surgical resection in the School and Hospital of Stomatology, China Medical University. All patients received no preoperative chemotherapy, immune therapy, and radiation therapy. The study was approved by the Ethics Committee of China Medical University. All study involving human participants were in accordance with the Declaration of Helsinki. All the enrolled participants provided written informed consent.

Paraffin-embedded tissue blocks were sectioned into 4-µm-thick sections for immunohistochemistry. All samples were deparaffinized and rehydrated. The citrate buffer (pH 6.0) and 3% hydrogen peroxide were used for antigen retrieval and endogenous peroxidase activity blocking, respectively. Then, all slides were incubated with goat anti-FOXD1 polyclonal antibody (1:100, Abcam Company, #ab129324) overnight at 4°C, followed by incubating with rabbit anti-goat horseradish peroxidase (HRP)-conjugated secondary antibody at room temperature for 30 min. Finally, all slides were visualized using DAB Horseradish Peroxidase Color Development Kit (Maixin Co., Fuzhou, China).

All data were analyzed using R software (version 3.6.3) and SPSS 22.0 (SPSS Inc., Chicago, IL, U.S.A.). All data were expressed as mean ± standard deviation (SD). Statistical differences were evaluated using Student’s t tests or one-way ANOVA. Survival analyses were evaluated using Kaplan–Meier curve and Log-rank test. The univariate and multivariate Cox regression analyses were performed to identify the potential prognostic values of FOXD1 and clinical parameters. Pearson correlation analysis was used to assess the relationship between FOXD1 expression and the proportion of immune cell subtype. P-value <0.05 was considered significant.

We analyzed tumor samples from TCGA database to identify FOXD1 expression characteristics (Figure 2A). FOXD1 expression was up-regulated in breast invasive carcinoma (BRCA), cholangiocarcinoma (CHOL), colon adenocarcinoma (COAD), esophageal carcinoma (ESCA), head and neck squamous cell carcinoma (HNSC), liver hepatocellular carcinoma (LIHC), lung squamous cell carcinoma (LUSC), prostate adenocarcinoma (PRAD), and stomach adenocarcinoma (STAD). FOXD1 expression was down-regulated in kidney renal clear cell carcinoma (KIRC), kidney renal papillary cell carcinoma (KIRP), thyroid carcinoma (THCA), and uterine corpus endometrial carcinoma (UCEC). Then, we used LinkedOmics tool to analyze the FOXD1-correlated genes (Supplementary Table S1) based on TCGA dataset and analyzed the GO functions and KEGG pathways with STRING database. The main GO functions and gene distributions, including biological process (BP), cell component (CC), and molecular function (MF), were shown in Figure 2B and Supplementary Table S2. The enrichment signal pathways were mainly distributed on proteoglycans in cancer, Rap1 signaling pathway, Ras signaling pathway, and T-cell receptor signaling pathway (Figure 2C).

According to inclusion criteria, a total of 449 HNSC samples and 44 normal samples from TCGA database were finally included in the study. The clinical characteristics were list in Supplementary Table S3. Compared with normal tissues, FOXD1 expression was significantly up-regulated in HNSC tissues in TCGA dataset (P<0.001, Figure 3A). To validate the result reliability, we also compared FOXD1 expression between paired HNSC tissues and adjacent normal tissues (n=43 pairs). FOXD1 was also up-regulated in HNSC tissues than normal tissues (P<0.001, Supplementary Figure S1). Moreover, in GSE30784 dataset, FOXD1 expression in HNSC tissues was significantly higher than benign dysplasia and normal tissues (P<0.001, Figure 3B).

Then, GEO datasets were as external validation. A total of 13 GEO datasets were included in the meta-analysis based on the inclusion criteria. These datasets were submitted from 2007 to 2020. FOXD1 expression in tumor tissues was significantly higher than normal tissues in GSE6631, GSE19089, GSE23558, GSE25099, GSE30784, GSE31056, GSE74530, GSE78060, and GSE85195 datasets (all P<0.05, Figure 3C). Meta-analyses indicated that the pooled standard mean difference (SMD) of FOXD1 was 0.91 (P<0.001, 95% CI: 0.47–1.35) using random-effects model, indicating that FOXD1 expression in HNSC was significantly higher than normal tissues. The funnel plot of SMD for the included studies appeared to be symmetric, and displayed no publication bias (Figure 3D).

In TCGA dataset, FOXD1 expression was significantly associated with tumor site (P=0.010) and HPV infection (P=0.010). But, FOXD1 expression was not associated with age (P=0.354), gender (P=0.483), distant metastasis (P=0.205), clinical stage (P=0.189), N stage (P=0.496), T stage (P=0.081), and grade (P=0.572) (Figure 4).

In TCGA dataset, Kaplan–Meier curves showed that high FOXD1 expression was significantly correlated with poor OS (OS: P<0.001, Figure 5A). GSE41613 and GSE65858 datasets were used as external validation. Although the difference was not significant in GSE41613 and GSE65858 datasets, high FOXD1 expression had a tendency of low OS (Figure 5B,C). Univariate and multivariate Cox regression analyses showed that FOXD1 expression was an independent prognostic factor for OS (P=0.001, Table 2). According to the result of multivariate Cox regression, the nomograms were generated to predict the 3- and 5-year OS probabilities, respectively. We introduced three independent factors into the nomogram for OS, including distant metastasis, T stage and FOXD1 expression, and each factor was assigned a score in proportion to its contribution to the risk of survival, thus obtaining the corresponding predicted survival rate (Figure 5D). The C-index of OS prediction was 0.62. Calibration of the nomograms displayed good agreement between the predicted 3‐ and 5‐year OS rates and actual observations (Figure 5E,F).

To understand the immune characteristics of FOXD1, ESTIMATE was applied to calculate the stromal score and immune score in 449 HNSC patients. We found that stromal scores and immune scores displayed no significant difference between high- and low-FOXD1 expression groups (P>0.05, Figure 6A). Then, the gene expression matrix of the HNSC dataset was analyzed using CIBERSORT algorithm to estimate the proportions of 22 immune infiltrating cells (Figure 6B). The proportions of naïve B cells (P=0.023), plasma cells (P=0.019), and resting dendritic cells (DCs) (P=0.018) in high FOXD1 expression were significantly lower than low FOXD1 expression group, while the proportion of activated mast cells in high FOXD1 expression group was significantly higher than low FOXD1 expression group (P=0.014). Moreover, we found that the proportions of naïve B cells (r = −0.141, P=0.006, Figure 6C), plasma cells (r = −0.170, P=0.001, Figure 6D), and resting DCs (r = −0.118, P=0.020, Figure 6E) were negatively correlated with FOXD1 expression, otherwise, the proportion of activated mast cells was positively correlated with FOXD1 expression (r = 0.189, P<0.001, Figure 6F).

FOXD1‐related signaling pathways were analyzed between high- and low-FOXD1 expression phenotypes through GSEAs. We list the top ten up-regulated signal pathways in the high-expression group of FOXD1, including metabolism of xenobiotics by cytochrome P450, arachidonic acid metabolism, steroid hormone biosynthesis, PPAR signaling pathway, cell adhesion molecule (Figure 7A). Moreover, the biological functions were enriched in keratinocyte differentiation, epidermal cell differentiation, complement activation, B-cell receptor signaling pathway, humoral immune response mediated by circulating immunoglobulin in high FOXD1 phenotype (Figure 7B).

To further validate the mRNA expression of FOXD1, we assessed the expression levels of FOXD1 in HNSC cell lines by qRT-PCR. The mRNA levels of FOXD1 were significantly higher in SCC-9 and CAL-27 cell lines compared with those of HOEC cell line (Figure 8A). FOXD1 protein levels were also significantly up-regulated in SCC-9 cell lines (Figure 8B,C). Moreover, FOXD1 expression in HNSC occurred mainly in the nucleus, and FOXD1 expression was significantly increased in cancer tissues compared with corresponding adjacent normal tissues (Figure 8D).

FOX family is an important and complex gene family, comprising diverse cell- and tissue-specific ‘winged-helix’ transcriptional factors [16]. Previous studies demonstrate that FOX family members participate in DNA damage repair, cell proliferation, differentiation, metabolism, apoptosis, tissue homeostasis, and development of immune system [17]. As a consequence, the alterations of FOX expression can influence the tumorigenesis as well as cancer progression. To date, several key FOX gene subfamilies, such as FOXA, FOXC, FOXM, FOXO and FOXP, have been found to be strongly linked to cancer initiation, invasion, metastasis, and drug resistance [18–21]. In the present study, we found that FOXD1 was up-regulated in HNSC tissues at both mRNA and protein levels, and significantly associated with primary tumor site and HPV infection. Moreover, FOXD1 expression was significantly correlated with OS in HNSC patients. In order to explore the function and molecular mechanism of FOXD1, we analyzed the correlation of FOXD1 expression and TME, and found that naïve B cells, plasma cells, resting DCs, and activated mast cells were significantly different between low- and high-FOXD1 expression groups. Moreover, GSEA revealed that FOXD1 was mainly involved in cancer-related signaling pathway and metabolism-related pathways.

The biological function of FOXD1 in human cancers has not been completely explored. FOXD1, as a newly discovered FOX family transcriptional factor, played a controversial role in varied cancer types. Some studies suggested that FOXD1 was involved in carcinogenesis and functioned as a tumor promoter in many cancer types [8]. Overexpression of FOXD1 could enhance the ability of cell proliferation and chemoresistance in MCF-7 breast cells, whereas silencing of FOXD1 expression attenuate cell proliferation and chemoresistance in MDA-231 cells [10]. Zhang et al. reported that FOXD1, as a novel oncogene, promoted the proliferation, migration, invasion, and radio-resistance in nasopharyngeal carcinoma cells [11]. In colorectal cancer patients, FOXD1 expression was up-regulalted and correlated with more invasive phenotype, such as lymphatic metastasis and TNM stage [8]. Tao et al. demonstrated that inhibitior of FOXD1 suppressed cell proliferation, migration and invasion in osteosarcoma cells, while overexpression of FOXD1 promoted osteosarcoma cell proliferation and migration [22]. In prostate cancer, down-regulation of FOXD1 affected the expression of cell cycle control genes and suppressed the androgen-independent growth of 22RV1 cells [13]. However, another study in ovarian cancer indicated that up-regulated FOXD1 could inhibit cell proliferation by inducing cell cycle arrest at G₁ phase in ovarian cancer cells, and high FOXD1 expression was significantly correlated with favorable prognosis [9]. The results obtained from breast cancer, lung cancer, nasopharyngeal carcinoma, colorectal cancer and prostate cancer are evident, but, the function of FOXD1 in ovarian cancer remain controversial. We performed a pan-cancer anlysis for FOXD1 at mRNA levels, and found that FOXD1 was up-regulated in BRCA, CHOL, COAD, ESCA, HNSC, LIHC, LUSC, PRAD and STAD, and down-regulated in KIRC, KIRP, THCA and UCEC. Due to non-available data of normal tissues, there was no comparison in ovarian carcinoma (OV). In view of the up-regulation of FOXD1 at mRNA levels in TCGA HNSC samples, we validated the FOXD1 expression in TCGA and GEO datasets. Thus, we speculated that FOXD1 was a tumor promoter in HNSC. Simialr with our result, a recent study reported that FOXD1 up-regulation was related to metastasis and poor clinical outcomes in oral squamous carcinoma [23]. Lin et al. also found that FOXD1 knockdown dramatically suppressed the colony-forming ability and confered radioresistance by down-regulating the JAK-STAT pathway in oral cancer cells [24]. In the present study, we found that mRNA and protein levels of FOXD1 was up-regulated in HNSC cell lines, and FOXD1 expression was increased in cancer tissues compared with corresponding adjacent normal tissues.

In the last decade, it is clear that FOX family members have potential roles in various aspects of immune regulation, immune homeostasis, and development and differentiation of immune cells [25,26]. FOXO subfamily members participated in the DC activity, CD8⁺ T-cell response, and macrophage activation [27]. FOXJ1 restrains B-cell activation and the maturation of humoral responses to regulate B lymphocyte homeostasis [28]. FOXP1 repressed immune signaling in the central nervous system and contributed to the immune dysfunction in central nervous system development [29]. FOXP3, a molecular marker of Tregs, is highly expressed in lymphoid tissue, and plays a key role in the development and function of Treg cells, representing an important therapeutic target for cancer [30]. FOXN1 is indispensable for thymic epithelial cell differentiation, growth, function, and thymic epithelium cell homeostasis [31,32]. However, the immune function of FOXD1 is little known and its role in TME has not been fully elucidated yet. Lin et al. reported that Foxd1 coordinated the regulation of the activity of NF-AT and NF-κB, and FOXD1 deficiency could result in multiorgan and systemic inflammation, and exaggerated Th cell-derived cytokine production and T-cell proliferation [33]. HNSC is an immunosuppressive disease, which is characteristic of decreased lymphocyte counts [34], impaired natural killer (NK) cell activity [35], down-regulated antigen-presenting function [36], and abnormality of tumor-infiltrating T lymphocytes [37]. Our results suggested that immune score decreased in HNSC tissues, although the difference was not significant (P=0.067). Moreover, we adopt CIBERSORT algorithms to analyze the proportions of immune cell subtypes and found that the naïve B cells, plasma cells, and resting DCs were significantly decreased with FOXD1 expression increase. Wouters et al. searched Pubmed databse, covering 69 studies and 19 cancer types, and provided the evidence to support a positive role for B cells and plasma cells in antitumor immunity [38]. Moreover, the preoperative DC counts and DC surface molecular expression were impaired in HNSC patients in comparison with healthy controls [39]. Furthermore, we investigated the correlation of FOXD1 expression and changed immune cells, and we found that the proportions of naïve B cells, plasma cells, and resting DCs were negatively correlated with FOXD1 expression; otherwise, the proportion of activated mast cells was positively correlated with FOXD1 expression.

However, our study has some limitation. First, the study was mainly based on TCGA and GEO datasets to investigate the prognostic value of FOXD1, while histological validation was conducted only in a few patients without completed follow-up information. Second, in vivo function experiments were lacking. Third, we analyzed the FOXD1-related pathways through GSEA, and there was a lack of relevant mechainsm research.

Taken together, we identified that up-regulation of FOXD1 was a potential prognostic marker in HNSC, and was associated with HPV status and HNSC sites. Moreover, FOXD1 expression played an important role in TME and immune cell infiltration. Hence, further studies are required to validate the role of FOXD1 in HNSC and to improve the understanding of the underlying mechanisms.

The datasets generated and analyzed during the current study are available from the TCGA and GEO database, and the data presented in this manuscript are available from the corresponding author on reasonable request.

The authors declare that there are no competing interests associated with the manuscript.

This work was supported by the National Natural Science Foundation of China [grant number 81301835].

T.j.W. and B.L. conceived and designed the experiments. B.L. and J.j.H. analyzed the data. T.j.W. and B.L. wrote the drafts of he paper, and T.j.W. approved the final draft.

The present study was approved by the Ethics Committee of China Medical University. All study involving human participants were in accordance with the Declaration of Helsinki.

We would like to thank TCGA and GEO database for sharing the data with the public.

BRCA

breast invasive carcinoma

CHOL

cholangiocarcinoma

CIBERSORT

Cell-type Identification By Estimating Relative Subsets Of known RNA Transcripts

COAD

colon adenocarcinoma

C-index

concordance index

DC

dendritic cell

ESCA

esophageal carcinoma

ESTIMATE

Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data

FOXD1

Forkhead box D1

GEO

Gene Expression Omnibus

GO

gene ontology

GSEA

Gene Set Enrichment Analysis

HNC

head and neck cancer

HNSC

head and neck squamous cancer

HPV

human papillomavirus

KEGG

Kyoto Encyclopedia of Genes and Genomes

KIRC

kidney renal clear cell carcinoma

KIRP

kidney renal papillary cell carcinoma

LIHC

liver hepatocellular carcinoma

LUSC

lung squamous cell carcinoma

NF-AT

nuclear factor of activated T cells

OS

overall survival

PPAR

peroxisome proliferator-activated receptor

STAD

stomach adenocarcinoma

TCGA

The Cancer Genome Atlas

THCA

thyroid carcinoma

TME

tumor microenvironment

TNM

tumor, nodes, metastasis

Treg

regulatory T cells

UCEC

uterine corpus endometrial carcinoma