The Expression and Survival Significance of FOXD1 in Lung Squamous Cell Carcinoma: A Meta-Analysis, Immunohistochemistry Validation, and Bioinformatics Analysis

Accumulating evidence demonstrated that FOXD1 dysregulation was correlated with a broad spectrum of malignancies. However, litter is known about the role of FOXD1 in the progression of lung squamous cell carcinoma (LUSC). We conducted the comprehensive bioinformatics analysis to investigate FOXD1 expression in LUSC from TCGA and GEO datasets, and validated the FOXD1 expression pattern in clinical samples using immunohistochemistry method. ESTIMATE and CIBERSORT algorithms were performed to assess the relationship of FOXD1 and tumor microenvironment and immune cell infiltration. Our study showed that FOXD1 expression was significantly upregulated in LUSC tissues in TCGA dataset, validated by GEO datasets and clinical samples. In TCGA dataset, Kaplan-Meier curves showed that high FOXD1 expression was significantly correlated with favorable prognosis in LUSC patients. Moreover, FOXD1 expression has an impact on immune score and the proportions of immune cell infiltration subgroups. Finally, we predicted FOXD1 may be involved in many immune-related biological functions and cancer-related signaling pathways. Taken together, FOXD1 was upregulated in LUSC tissues, and FOXD1 expression could be a potential prognostic marker. FOXD1 might be associated with tumor microenvironment and perhaps a potential target in the tumor immunotherapy.


Introduction
According to the last GLOBOCAN estimates, the incidence and mortality of lung cancer ranked first worldwide in 2018, accounting for 2,093,876 new cases (11.6% of all new cases of malignant tumors) and 1,761,007 deaths (18.4% of all tumor deaths) [1]. Lung cancers are divided into two main subtypes, non-small cell lung cancer (NSCLC) and small cell lung cancer (SCLC). Lung squamous cell carcinoma (LUSC), a subtype of NSCLC, represents of 30% of NSCLC [2]. Previous studies demonstrated that the molecular mechanism of LUSC is associated with gene mutation, multiple altered expression of genes and pathways, and chromosomal instabilities in the progression of LUSC [3]. Moreover, smoking is one of the major risk factors in the development and progression of LUSC, and the rate of smoking exposure in LUSC patients exceeds 90% [4]. Although diagnosis and treatment have made remarkable progress in recent years, no specific biomarkers or relatively optimal targeted therapies have been identified for LUSC patients [5]. In LUSC patients, the 5-year overall survival rate for stage I and stage II patients is approximately 40%, while the overall survival (OS) rate is less than 5% for stage III and stage IV patients [6]. Therefore, it is of importance to unveil the molecular mechanism of LUSC and explore the target therapy strategies for LUSC patients.
The forkhead box (Fox) transcription factors (TFs) play important roles in biological processes, including cell growth, differentiation, proliferation, apoptosis, and longevity [7,8]. Accumulating evidence has demonstrated that their mutation and dysregulation have been correlated to a broad spectrum of malignancies [9]. FOXD1, as a new member of the FOX family, is mainly located on 5q13.2 and encodes a DNA-binding protein containing 465 amino acids. Recently, many studies indicated that FOXD1 was involved in the development and progression of different types of human cancers, and its dysregulation was mainly associated with cell proliferation, migration, invasion, radioresistance, and epithelial-to-mesenchymal transition (EMT) [10][11][12]. Sohei Nakayama, et al. demonstrated that knockdown of FOXD1 could suppress cell growth in lung cancer cell lines, and high FOXD1 mRNA level was correlated with poor prognosis [13]. Li D, et al. indicated that FOXD1 has oncogenic characteristics by activating vimentin in NSCLC cell lines [14]. In LUSC, miR-30a-5p may inhibit the proliferation and migration of LUSC cells by inhibiting FOXD1 expression [15]. Although FOXD1 is involved in the pathological development of lung cancer, especially in LUSC, the related functions of FOXD1 contributing to LUSC progress have remained largely unknown.
In the present study, we analyzed FOXD1 expression in LUSC using The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) datasets, and validated the result by immunohistochemistry method. Then, we evaluated the association of FOXD1 expression and clinical parameters. Subsequently, we analyzed the relation between FOXD1 expression and tumor microenvironment (TME), as well as the immune cell infiltration.

Data Collection.
The transcriptome data and clinical information of LUSC were downloaded from TCGA database (https://www.cancer.gov/tcga). A total of 495 LUSC samples with detailed clinical information were enrolled in the study, as shown in Table 1. We searched GEO database (http://www.ncbi.nlm.nih.gov/geo/) to collect LUSC datasets up to December, 2020.The search terms included the following keywords: (lung) AND (tumor OR cancer OR carcinoma OR neoplasm). Each GEO dataset should meet the following criteria: (1) The sample includes cancer tissue and normal tissue; (2) The number of each group was greater than three. FOXD1 expression data was extracted from TCGA and GEO datasets. Then, the flow diagram of the study process was shown in Figure 1.

FOXD1
Expression in LUSC Samples from TCGA and GEO Datasets. FOXD1 expression was compared between LUSC tissues and adjacent noncancerous tissues in TCGA dataset. Subsequently, a total of 23 GEO datasets were included in the meta-analysis. The features of GEO datasets were listed in Table 2. We analyzed FOXD1 expression in all GEO datasets using Review Manager software (RevMan). The standard mean difference (SMD) and 95% confidence interval were used to estimate the expression pattern of FOXD1. Evidence of bias was assessed by visual funnel plots.
2.3. Immunohistochemistry Validation. Primary LUSC tissues were collected from LUSC patients, who underwent surgical resection in Shengjing Hospital, China Medical University. The study was approved by the Ethics Committee of Shengjing Hospital, China Medical University. All studies involving human participants were in accordance with the Declaration of Helsinki. All enrolled participants provided written informed consent.
First, we sectioned paraffin-embedded tissue blocks into 4-μm thick sections for immunohistochemistry (IHC). All samples were deparaffinized and rehydrated in a series of xylene and graded ethanol solutions. 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 FOXD1 antibody (1 : 100, Abcam Company, #ab129324) overnight at 4°C. Tissue sections were then incubated with rabbit anti-goat horseradish peroxidase (HRP)-conjugated secondary antibody at room temperature for 1 h. Finally,      BioMed Research International all slides were visualized using DAB Horseradish Peroxidase Color Development Kit (Maixin Co., Fuzhou, China).

Association of FOXD1 Expression and Clinical
Parameters. The association of FOXD1 expression and clinical parameters were evaluated using TCGA dataset. The clinical parameters included age, gender, distant metastasis, clinical stage, tumor size, node metastasis, and smoking. P <0.05 was considered as statistically significant.
2.5. Prognosis of FOXD1 Expression in LUSC Patients. The prognosis analysis was conducted in TCGA (n =495 cases), GSE3141 (n =53 cases), GSE13213 (n =117 cases), GSE14814 (n =52 cases), GSE37745 (n =66 cases), and GSE50081 (n =43 cases) datasets. We divided LUSC patients into two groups according to mean of FOXD1 expression, high FOXD1 expression group and low FOXD1 expression group. Kaplan-meier curve was plotted to evaluate the prognostic value of FOXD1, examined by Log-rank test. In TCGA dataset, the univariate and multivariate Cox proportional hazards models were performed orderly to screen the prognostic factors among the clinical parameters, as well as FOXD1 expression. P <0.05 was considered as statistically significant. Then, a nomogram was constructed on the basis of FOXD1 expression and clinical parameters, which were the independent prognostic factors in multivariate Cox regression analysis. Moreover, we used C-index and calibration plots to evaluate the performance of the prediction model. The model was visualized using R programming (R 3.6.2; http://www.r-Project.org).

Correlation of FOXD1 Expression and Immune Cell
Infiltration. The stromal cells and immune cells are two major components in TME. The levels of stromal cells and immune cells in LUSC tissues were computed by Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) tool (https:// bioinformatics.mdanderson.org/estimate/) based on gene expression data. The stromal and immune scores were compared between high FOXD1 expression group and low FOXD1 expression group.
To further explore the influence of FOXD1 in the immune microenvironment, we evaluated the proportion of immune cell infiltration for each LUSC sample using the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm. CIBERSORT is a deconvolution algorithm, which could infer the abundance of 22 immune cell types based on gene expression data [16]. Then, we compared the difference of immune cell infiltration between high FOXD1 expression group and low FOXD1 low expression group. Furthermore, we analyzed the correlation of FOXD1 expression and immune cell subtypes. A P <0.05 was considered as statistically significant.

Function Enrichment Analysis for FOXD1-Related
Genes. FOXD1-related genes were screened using LinkedOmics online tool (http://www.linkedomics.org/login.php). The screening threshold was set as: |r| >0.4 and P <0.05. Gene ontology (GO), including biological process (BP), molecular function (MF), and cell component (CC) were used to depict the function of FOXD1-related genes using STRING online tool (http://www.string-db.org). Then, we    BioMed Research International analyzed the potential miRNAs and transcription factors, which could regulate FOXD1 expression.
2.8. GSEA Analysis. We performed gene set enrichment analysis (GSEA) to explore the potential altered signaling pathways between high FOXD1 expression group and low FOXD1 expression group. The significantly enriched gene sets were identified with a nominal P value <0.05 and the |enrichment score| >0.4.
2.9. Statistical Analysis. All statistical analyses were conducted using SPSS 22.0 (SPSS, IL, United States) and R. 3.6.2 (https://www.r-project.org/). The correlation between FOXD1expression and various clinical features was analyzed by t-test and analysis of variance (ANOVA). Survival analysis was estimated by Kaplan-Meier curve, and tested by Logrank test. Cox proportional risk regression models was used to assess independent risk factors for overall survival (OS) of LUSC patients. A P <0.05 was considered statistically significant.

FOXD1 Expression Pattern in LUSC Samples.
We first analyzed FOXD1 expression between LUSC tissues and adjacent noncancerous tissues in TCGA dataset, and found that FOXD1 was upregulated in LUSC tissues (P <0.05, Figure 2(a)). In order to verify the difference in FOXD1 expression in TCGA dataset, a comprehensive metaanalysis was conducted based on GEO datasets. A total of 23 GEO datasets were included in the study. The result showed that FOXD1 was upregulated in LUSC tissues in GEO datasets (SMD =1.05, 95% CI: 0.54-1.56, P <0.001, Figure 2(b)). The funnel plot of SMD for the included studies appeared to be symmetric, and displayed no publication bias (Figure 2(c)). We also performed IHC to validate FOXD1 protein expression in clinical samples. IHC results indicated that FOXD1 in LUSC tissues was significantly upregulated than normal lung tissues (Figure 2(d)).

Association of FOXD1 Expression and Clinical
Parameters. Moreover, we analyzed the correlations between FOXD1 expression and clinical parameters, including age, gender, TNM stage, clinical stage, and smoking history. FOXD1 expression displayed no significant differences in all clinical characteristics (all P >0.05, Figures 3(a)-3(g)). In stage IV and N3 stages, FOXD1 showed a tendency of upregulation, but the difference was not significant (P >0.05).

FOXD1 Expression Was Positively Correlated with
Overall Survival in LUSC Patients. In TCGA dataset, we plotted Kaplan-meier curves, examined by Log-rank test, and found that high FOXD1 expression was positively correlated with OS in LUSC patients (P =0.036, Figure 4(a)). Moreover, we found that LUSC patients with FOXD1 expression had a longer OS in comparison with those with low FOXD1 expression in GSE3141, GSE14814, GSE37745, GSE13213, and GSE50081 datasets (Figures 4(b)-4(f)). Due to the small sample size, the differences were not significant.
As shown in Table 3

Function Enrichment Analysis of FOXD1-Related Genes.
A total of 360 FOXD1-related genes were identified through transcriptome co-expression analysis, including 308 positively correlated genes and 52 negatively correlated genes. We list the top 50 positive and negative FOXD1-related genes using heatmaps (Figure 6(a) and 6(b)). GO analyses showed that BP were mainly involved in immune response, positive regulation of T cell proliferation, cell adhesion, inflammatory response, and leukocyte migration; and MF were mainly involved in receptor activity, cytokine receptor activity, MHC class II receptor activity, SH3/SH2 adaptor activity, and chemokine receptor activity (Figure 6(c)). Moreover, we predicted the potential miRNAs (miR-328, miR-485-5P, and miR-346) and transcription factors (ZR5, VDR, NF-κB, PEA3, ELF1, et al.), which could regulate FOXD1 expression (Figures 6(d) and 6(e)).

GSEA Analysis.
The differentially expressed genes between high FOXD1 expression and low FOXD1 expression groups were conducted to perform KEGG pathway analyses. The results showed that enrichment pathways included glycosaminoglycan biosynthesis keratin sulfate, hedgehog signaling pathway, Wnt signaling pathway, and ERBB signaling pathway in high FOXD1 expression group (Figures 7(a)-7(d)). Otherwise, a series of immune-related pathways, such as antigen processing and presentation, cell adhesion molecules cams, cytokine cytokine receptor interaction, complement and coagulation cascades, natural killer cell mediated cytotoxicity, JAK-STAT signaling pathway, primary immunodeficiency, and chemokine signaling pathway, were mainly enriched in low FOXD1 expression group (Figures 7(e)-7(l)).

Discussion
FOX proteins are a superfamily of evolutionarily conserved transcription factors, which share 'winged helix' DNAbinding domain (FOX domain) [9]. Recent studies identified that FOX domains have been found in over 100 proteins, ranging from FoxA to FoxQ subclasses [7,17]. FOX family genes are involved in carcinogenesis as oncogenes and/or cancer suppressor genes. FOXD subfamily consists of 9 members, including FOXD1, FOXD2, FOXD3, FOXD4, FOXD5, FOXD6, FOXDL4, FOXDL5, and FOXDL6. Initially, the function of FOXD1 was identified in renal morphogenesis and the development of optic chiasm [18,19]. Recently, the roles of FOXD1 have begun to be elucidated in the development and progression of tumors. In breast cancer, FOXD1 could induce G1 to S phase transition, thus promoting cell proliferation and chemoresistance [20]. In lung cancer, FOXD1 coupled with Gal-3 increased tumor growth and motility, whereas depletion of Gal-3 attenuated FOXD1-mediated tumorigenesis [21]. Moreover, FOXD1 promoted cell proliferation, migration and invasion in colorectal cancer cells by regulating the phosphorylation of ERK 1/2 pathway [22]. In melanoma, overexpression of FOXD1    [24]. Therefore, FOXD1 plays important roles in many cancers, such as proliferation, metastasis, and drug resistance, but some cancers are on the contrary. In the present study, we found that the levels of FOXD1 mRNA were upregulated in LUSC tis-sues in TCGA and GEO datasets. Moreover, we collected clinical LUSC samples and validated FOXD1 protein upregulation using IHC experiments.

22
BioMed Research International expression had longer survival time in our study. As every coin has two sides, one possible reason may be that FOXD1 was an inducer in the progression of LUSC, activating and/or inhibiting some potential signal pathways.
Another explanation could be the existence of multiple influential factors that may lead to increased expression of FOXD1 under different circumstances. Other explanation was that the prognosis of LUSC was not determined   Regulation of the immune system is a critical part of anticancer therapies including immunotherapy, chemotherapy, and radiotherapy [25]. Our results indicated that high FOXD1 expression group showed a significantly lower immune score in comparison with low FOXD1 expression. Growing evidence has demonstrated that the immune response has antitumor effects, which immunologically mediated elimination of transformed cells has been widely accepted in the context of cancer for many decades [26]. FOXD1 low expression was companied by the low immune score, implying that FOXD1 was involved in the modulation of immune response in LUSC. Although still in debate, amounting evidence suggests that CD4 + T cells play an important role as a mediator in the maintenance and control of protective immune responses [27,28]. Monocytes are one of the most abundant cells in the solid tumor bulk. Initially, monocytes contribute antitumor functions, and finally become tumor-supportive and immunosuppressive undergoing a phenotypic switch [29]. Macrophages are versatile immune cells, which are polarized to two opposite types, classically activated M1 macrophages and alternatively activated M2 macrophages. M1 macrophages, as an antitumor phenotype, exert an immune protective role by producing chemokines and cytokines to destroy tumor cells, whereas M2 macrophages protect cancer cells from antitumor immune responses and contribute to tumor progression [30,31]. Our results showed that resting memory CD4 + cells, monocytes, and M1 macrophages in high FOXD1 expression group were significantly lower than low FOXD1 expression group, indicating that the immune protective function was dampened with FOXD1 downregulation.
The present study had some limitations. First, most samples were downloaded from TCGA and GEO datasets, and clinical samples in IHC experiments were absent from detailed clinical information, which may affect the results to some extent. Second, some important treatment information, such as chemotherapy and immunotherapy, were incomplete, so that we did not fully analyze the association of FOXD1 expression and immune therapy. In addition, functional experiments should be performed to further elucidate the molecular mechanisms of FOXD1.

Conclusion
In summary, our study illustrated that FOXD1 was upregulated in LUSC samples and could predict the prognostic outcome in LUSC patients. Moreover, FOXD1 expression was correlated with immune infiltration. Therefore, FOXD1 could be a new target gene, which provides a new therapeutic target in LUSC. Further studies are required to investigate its molecular function.

Data Availability
The raw data can be obtained from public TCGA and GEO database. All data are available from the authors on reasonable request.

Conflicts of Interest
All authors declare there are no competing interests.