Dysregulation of FOXD2-AS1 promotes cell proliferation and migration and predicts poor prognosis in oral squamous cell carcinoma: a study based on TCGA data

High FOXD2-AS1 expression in OSCC

To verify the FOXD2-AS1 expression in cancer, we first explored its expression using pan-cancer analysis, which included tumor tissues and paracancerous tissues. FOXD2-AS1 was upregulated in most cancer types (Figure 1A, 1B), including head and neck squamous cell carcinoma (HNSCC). The Wilcoxon rank-sum test was used to compare the expression of FOXD2-AS1 in normal, genotype-tissue expressed (GTEx) tumor samples with HNSCC samples from TCGA (Figure 1C). We excluded non-oral samples (Hypopharynx, Larynx, Oropharynx, Tonsil) in this analysis. A total of 331 OSCC patients with FOXD2-AS1 expression were obtained from TCGA. FOXD2-AS1 was significantly upregulated in tumor tissues relative to nontumor tissues (Figure 1D). We also compared the expression of FOXD2-AS1 in tumor tissues with their corresponding paracancerous samples and found that FOXD2-AS1 was highly expressed in the OSCC samples (Figure 1E). The receiver operating characteristic analysis revealed that FOXD2-AS1 can be used as a marker to distinguish tumor from non-tumor cells (Figure 1F).

Figure 1. Expression level of FOXD2-AS1 was upregulated in most cancer types including OSCC. (A) FOXD2-AS1 expression in pan-cancer from normal TCGA samples of GTEx combined with samples of TCGA head and neck squamous cell carcinoma. (B) FOXD2-AS1 expression in pan-cancer of tumor tissues and paracancer tissues from TCGA. (C) FOXD2-AS1 expression of normal samples and tumor samples in TCGA head and neck squamous cell carcinoma. (D) FOXD2-AS1 expression of normal samples and tumor samples in OSCC patients. (E) FOXD2-AS1 expression of tumor tissues and paracancer tissues in OSCC patients. (F) ROC analysis of FOXD2-AS1 expression showing promising discrimination power between tumor and normal tissues.

Association of FOXD2-AS1 expression with clinicopathological parameters

To further elucidate the FOXD2-AS1 involvement in the development of OSCC, we analyzed the correlation of FOXD2-AS1 expression with clinicopathological parameters. In this study, the clinical information of 331 OSCC patients from TCGA were analyzed; among which, 228 male and 103 female patients with a median age of 62 years old (range: 52 to 71) were enrolled. Other clinicopathological features are shown in Table 1. The expression of FOXD2-AS1 in OSCC tissues was labeled as “low” or “high expression”, based on median values. High FOXD2-AS1 expression was strongly associated with radiation therapy, anatomic localization, high histologic grade, and lymphovascular invasion (P values were 0.016, 0.033, <0.001, and 0.011 respectively) (Figure 2A–2D). Other characteristics had no correlation with FOXD2-AS1 expression; detailed information is shown in Table 1. Logistic regression analysis of differential FOXD2-AS1 expressions is shown in Table 2.

Figure 2. The association between FOXD2-AS1 expression and clinicopathological parameters. (A) The history of radiation therapy. (B) The anatomic localization of tumor. (C) The histologic grade of tumor. (D) The status of lymphovascular invasion.

Prognostic value of FOXD2-AS1 in OSCC

Kaplan-Meier curves of OS and disease-specific survival (DSS) were plotted to analyze the prognosis of OSCC patients using different expressions of FOXD2-AS1. Log-rank test of OS and DSS revealed that a high expression of FOXD2-AS1 was significantly associated to shorter survival (Figure 3A, 3B). We then performed a subgroup analysis to further explore the correlation between FOXD2-AS1 expression and patient prognosis with different characteristics. (Figure 3C–3H). The results showed that in non-radiated, histologic grades I and II, and non-lymphovascular invasive groups, a high FOXD2-AS1 expression predicted poor survival. Univariate Cox regression analysis revealed that perineural invasion, lymphovascular invasion, radiation therapy, and FOXD2-AS1 expression were prognostic factors for overall survival (OS) and DSS, while multivariate Cox regression showed that FOXD2-AS1 expression was an independent risk factor in OSCC patients (Tables 3, 4).

Figure 3. Prognostic value of FOXD2-AS1 in OSCC. (A, B) Kaplan-Meier curves of OS and DSS with different expression level of FOXD2-AS1. (C, D) Kaplan-Meier curves of OS and DSS with different expression level of FOXD2-AS1 in non-radiation therapy group. (E, F) Kaplan-Meier curves of OS and DSS with different expression level of FOXD2-AS1 in histologic grade I and grade II group. (G, H) Kaplan-Meier curves of OS and DSS with different expression level of FOXD2-AS1 in non-lymphovascular invasion group.

Establishment of survival prognostic models for OSCC

Since the results mentioned above suggested that FOXD2-AS1 was an independent prognostic factor in OSCC, we tried to establish a prediction model of OS and DSS by fitting the expression of FOXD2-AS1 with the clinicopathological parameters. We constructed a nomogram to integrate FOXD2-AS1 and other prognostic factors, including the tumor (T) stage, perineural invasion, radiation therapy, and lymphovascular invasion (Figure 4A, 4B). A worse prognosis was represented by a higher point on the nomogram. The nomogram performance with FOXD2-AS was evaluated using a calibration curve, and the C-index of OS and DSS were 0.715 and 0.716, respectively (Figure 4C, 4D), suggesting that this nomogram might be a better model for predicting survival in OSCC patients than individual prognostic factors.

Figure 4. Establishment and evaluation of survival prognostic models for OSCC. (A, B) Nomogram for predicting the probability of 1-, 3- and 5-year OS and DSS for OSCC patients. (C, D) Calibration curve for evaluation the accuracy of nomogram of OS and DSS for OSCC patients.

The identification of differentially expressed genes (DEGs) in High FOXD2-AS1 expression and Low FOXD2-AS1 expression samples via functional cluster analysis

To explore the potential mechanism of FOXD2-AS1 in causing tumor progression, we analyzed DEGs in high and low FOXD2-AS1 expression samples. A total of 478 DEGs were identified, with 348 upregulated genes and 130 downregulated genes. The expression of DEGs was demonstrated using a Heatmap and a Volcano Plot (Figure 5A, 5B). The functions of co-expression in OSCC patients were then predicted by GO and KEGG enrichment analyses. The top GO terms in the biological process (BP), molecular function (MF), and cellular component (CC) groups were keratinization, structural constituent of muscle, and cornified envelope, respectively. KEGG analysis revealed that the retinol metabolism was enriched (Figure 5C–5F). We also performed a GSEA to identify the key pathway related to FOXD2-AS1. The most significantly enriched pathways based on the NES were the E2F targets and G2/M checkpoint (Figure 5G, 5H).

Figure 5. Differentially expressed genes (DEGs) in High FOXD2-AS1 expression and Low FOXD2-AS1 expression samples and functional cluster analysis. (A) Heat map of top 5 DEGs in high and low FOXD2-AS1 expression groups. (B) Volcano map of DEGs with |log2FoldChange >2| and adjusted P value<0.05. (C–F) Significantly enriched GO and KEGG annotations of FOXD2-AS1 related genes. (G, H) Enrichment of genes in E2F targets and G2/M checkpoint by GSEA.

The correlation between FOXD2-AS1 expression and immune infiltration

We analyzed the correlation between FOXD2-AS1 expression and immune cell infiltration level (generated by ssGSEA) using Spearman correlation, as shown in Figure 6A. The results indicated that FOXD2-AS1 expression was strongly negatively correlated with mast cells. The Spearman R value was -0.395 with P<0.001. The Wilcoxon rank-sum test also indicated that the enrichment score of mast cells was higher in low FOXD2-AS1 expression group than in the high expression group (Figure 6B, 6C). Other immune cells with statistically significant enrichments were neutrophils, DCs, iDCs, and B cells (Figure 6D–6K).

Figure 6. The expression level of FOXD2-AS1 was associated with the immune infiltration in tumor environment. (A) The plots showed the correlation between FOXD2-AS1 expression and immune cells subsets. (B) Spearman correlation between expression of FOXD2-AS1 and mast cells. (C)The plots of mast cells expression in low and high FOXD2-AS1 samples. (D–K) Spearman correlation and expression distribution of Neutrophils, DC cells, iDC cells and B cells in low and high FOXD2-AS1 samples.

FOXD2-AS1 is upregulated in OSCC tissues and cell lines and promotes cell proliferation and migration in vitro and in vivo

Twenty-five paired OSCC tissue samples and 3 OSCC cell lines were used to investigate the expression of FOXD2-AS1 via qPCR. FOXD2-AS1 expression was elevated in both tumor tissues and cell lines (Figure 7A, B). The subcellular localization was performed by qPCR and the results indicated that FOXD2-AS1 was localized predominantly in the nucleus (Figure 7C). Since the Cal27 cell lines had the highest expression of FOXD2-AS1, we performed a loss of function test in Cal27 using an siRNA and the Smart Silence Kit (Figure 7D). First, we validated the expression of the top 3 upregulated and downregulated genes filtered by TCGA. The results showed that the expression of PNCK, HOXB8, and HOXB9 was also decreased after the knockdown of FOXD2-AS1, while KRT6C and DSG1 were upregulated; this result was consistent with a previous TCGA analysis (Figure 7E). Next, we investigated the function of FOXD2-AS1 in the Cal27 cell line. Knockdown of FOXD2-AS1 remarkably impaired cell migration ability, as demonstrated in the wound healing and Transwell migration assays (Figure 7F, 7G). The mRNA and protein levels of N-cadherin and Snail1 were decreased, while the expression of E-cadherin was upregulated, as measured using qPCR and western blotting (Figure 7H, 7I). The silencing of FOXD2-AS1 also inhibited the proliferation of cancer cells according to the results of the CCK-8 and colony formation assays (Figure 7J, 7K). We also performed in vivo knockout experiments by injecting FOXD2-AS1 antisense oligonucleotides (ASOs) into nude mice. A xenograft tumor model was established via subcutaneous injection of Cal27 cells into nude mice. The injection of ASOs inhibited tumor growth (Figure 7L) and reduced tumor volume (P value < 0.001) (Figure 7M) and tumor weight (P value < 0.05) (Figure 7N). These results showed that FOXD2-AS1 was highly expressed in OSCC and can promote the proliferation and migration of cancer cells.

Figure 7. Knockdown of FOXD2-AS1 decreased OSCC cell proliferation and migration in vitro and in vivo. (A) The relative expression of FOXD2-AS1 in 25 paired tumor and adjacent normal tissues. The expression of FOXD2-AS1 was significantly up-regulated in Tumor tissues. (B) The relative expression of FOXD2-AS1 in 3 OSCC cell lines and Normal oral primary keratinocytes (NOK) cells. (C) Cell nucleus/cytoplasmic fractionation and qRT-PCR showed the subcellular localization of FOXD2-AS1 in Cal27 cells (MALAT1 and β actin were used as separation quality standards and endogenous controls). (D) The relative expression of FOXD2-AS1 after transfected with SiRNA and Smart Silence in Cal27 cell line. (E) The relative expression of PNCK, HOXB8, HOXB9, KRT6C, EPG and DSG1 after FOXD2-AS1 knockdown by SiRNA. (F) FOXD2-AS1 knockdown inhibited the migration capacity of Cal27 cells as detected by wound-healing assay after 48hs. (G) The cell migration abilities of Cal27 cells transfected with SiRNA and Smart silence were determined by transwell assays. (H) The relative mRNA expression of E-cadherin, N-cadherin and Snail in Cal27 cells transfected with SiRNA and Smart silence. (I) The protein expression of E-cadherin, N-cadherin and Snail in Cal27 cells transfected with SiRNA and Smart silence. FOXD2-AS1 knockdown inhibited the proliferation capacity of Cal27 cells as detected by CCK-8 assay (J) and colony formation assay (K). (L) Tumors were collected from different groups after 8 times ASO-FOXD2-AS1/ASO-NC injection. (M) Volume curve of tumors from different groups at the time of injection. (N) Weight of collected tumors from different groups. *p < 0.05, **p < 0.01, ***p < 0.001.

RNA-sequencing data from TCGA

Gene expression data and corresponding clinical information of OSCC patients (331 cases, Workflow type: HTseq-FPKM) were downloaded from the HNSCC projects of TCGA (https://genome-cancer.ucsc.edu/). Patients diagnosed with OSCC with complete follow-up information were included. Next, the level 3 HTseq-FPKM data were transformed into transcripts per million reads (TPM) for further analyses. Unavailable and unknown clinical features were regarded as missing values. This study met the publication guidelines stated by TCGA. All data used in this study were obtained from TCGA, thus, ethics approval and informed consent were not required.

DEG analysis

The samples were divided into high- and low-expression groups according to the level of FOXD2-AS1 expression (median as cut-off value). The expression data (TPM) of the high and low FOXD2-AS1 expression groups were then compared to identify DEGs using the limma package in R [41]. log₂(Fold Change) >2 and adjusted P value <0.05 were set as threshold values for DEGs.

Metascape analysis

Metascape, a free, well-maintained, and user-friendly gene list analysis tool for gene annotation and analysis [42], was used to conduct pathway and process enrichment analyses of DEGs. The conditions for differences to be deemed statistically significant included a P value <0.01, minimum count of 3, and an enrichment factor >1.5.

Gene set enrichment analysis

GSEA was used to elucidate the significant function and pathway difference between high- and low- FOXD2-AS1 groups. Gene set permutations were performed 1000 times for each analysis. The expression level of FOXD2-AS1 was used as a phenotype label. Adj. P value< 0.05 and FDR< 0.25 were identified as significantly related genes. Statistical analysis and graphical plotting were conducted using the R package clusterProfiler [43].

Analysis of immune infiltration and its correlation with FOXD2-AS1 expression

The immune infiltration analysis of OSCC was conducted via the single sample GSEA (ssGSEA) method using the GSVA package (http://www.bioconductor.org/packages/release/bioc/html/GSVA.html) in R. Based on the signature genes of 24 types of immunocytes in literature [44], the relative enrichment score of every immunocyte was quantified from gene expression profiles of each tumor sample. Wilcoxon rank-sum test and Spearman correlation were used to evaluate the association between FOXD2-AS1 expression and the 24 types of immune cells.

Tissue samples and cell culture

Paired, fresh OSCC tissues and their corresponding adjacent normal tissues were collected from the Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine. Written informed consent was obtained from each patient. The patients were well-informed and have consented to this study, and the process was approved by Ethics Committee of the Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine.

Cal27, HN6, and SCC9 cells and normal oral primary keratinocytes (NOK) were used in this study. Cal27, SCC9 and HN6 were purchased from the Type Culture Collection of Chinese Academy of Sciences (Shanghai, China). NOK were cultured from gingival tissues after tooth extraction from healthy patients. These cell lines were maintained in DMEM (Basalmedia, Shanghai) supplemented with 10% fetal bovine serum, 1% glutamine, and 1% penicillin/streptomycin. Cells were cultured in a standard humidified atmosphere with 5% CO₂ at 37° C.

RNA extraction and qRT-PCR

Total RNA was extracted using TRIzol (Invitrogen, CA, United States). An equal amount of RNA was reverse-transcribed using the HiScript II Q RT Supermix and was quantified by qPCR using SYBR Green (Bimake). The primer sequences are shown in Supplementary Table 1.

Western blotting

Cells transfected with SiRNA or the Smart Silence Kit for 48 h in 6-well plates were rinsed with PBS. Then, 150 μL of 1× SDS Protein Loading Buffer containing protease and phosphatase inhibitors were added to each well. Cell lysates were denatured at 98° C for 10 min and were subjected to 10% SDS-PAGE. The protein bands were then transferred to PVDF membranes, which were then blocked with 5% skimmed milk for 1 h and incubated with antibodies at 4° C overnight. The antibodies used in western blotting included those targeting β-actin (ABclonal AC026), E-cadherin (ABclonal A11492), N-cadherin (ABclonal A3045), and Snail (ABclonal A11794). HRP-conjugated secondary antibodies were used for detection.

Isolation of nuclear and cytoplasmic RNA

Nuclear and cytoplasmic RNA was isolated using the PARIS™ kit (Thermo Fisher Scientific) according to the manufacturer’s instructions. The isolated RNAs were reverse transcribed and amplified via qPCR as described above. MALAT1 were used as an endogenous nuclear control, while β-actin were used as an endogenous cytoplasmic control.

RNA interference

SiRNA against FOXD2-AS1 was designed and synthesized by Genomeditech Shanghai. The Smart Silence Kit was designed and synthesized by Guangzhou RiboBio Co., Ltd. (Guangzhou, China). These were transfected into the cells with Lipofectamine 3000 reagent (Invitrogen) according to the manufacturer’s instructions. The sequences are provided in Supplementary Table 2.

Cell counting kit-8 (CCK-8) analysis

Cells transfected for 24 h with Smart Silencer/SiRNA were seeded into 96- well plates at a density of 1000 cells/well in triplicate. 10 μL of CCK-8 reagent (Dojindo, Kumamoto, Japan) was added to 100 μL of culture medium. The cells were subsequently incubated for 2 h at 37° C, and the optical density was measured at 450 nm and 600 nm using a microplate reader (Multiskan™ Sky Spectrophotometer, Thermo Scientific, USA).

Colony formation assay

Cells transfected for 24 h with Smart Silencer/SiRNA were seeded into 6-well plates at a density of 1000 cells per well and incubated for 10–14 days to form cell colonies. The colonies were fixed with 4% paraformaldehyde (Sangon Biotech, Shanghai, China) for 15 min and stained with Giemsa stain (Solarbio, Beijing, China) for 30 min; those with more than 50 cells were counted under a dissecting microscope.

Wound healing assay

The Cal27 cells were seeded into 6-well plates and transfected when cultures reached 50% confluence. The cell layers were then scratched using a 10-μL plastic pipette tip to produce wounds, which were photographed at 0 and 48 h under an inverted phase-contrast microscope. Three random fields were marked and measured. All assays were carried out in triplicate.

Transwell migration assay

Cell migration assays were performed using 24-well Transwell chambers with 8-μm porosity polycarbonate filters and Transwell insert chambers (Corning, USA). A total of 200 μL of cell suspension in serum-free medium was added into each upper chamber, while 500 μL of DMEM supplemented with 10% FBS was added to the lower chambers as a chemoattractant. After incubating for 24 to 36 h, the migrated cells were fixed with 4% paraformaldehyde (Sangon Biotech) for 15 min and stained with Giemsa stain (Solarbio) for 30 min. After the cells on the upper surface of the filter were removed, at least 3 randomly selected microscopic fields of the fixed cells per filter were imaged using an inverted phase-contrast microscope. The cells were counted and cell counts were averaged.

In vivo assays

SPF nude mice (5 weeks) were purchased from the Shanghai Laboratory Animal Center (Shanghai, China) and housed under SPF conditions at the animal care facility of the Experimental Animal Center of Shanghai University of Traditional Chinese Medicine. A total of 2×106 Cal27 cells in 200 μL serum-free DMEM were subcutaneously injected into the flanks of the nude mice. The initial tumor xenograft volumes were similar. When the tumor size reached 100 mm³, the mice were randomized into 2 treatment groups: (a) ASO-NC group (5 nM/3 days); (b) ASO-FOXD2-AS1 group (5 nM/3 days). The ASO used was chemically modified by Ribobio. Its sequence was the first one of the Smart Silence Kit. Tumor size and animal weight were monitored at the time of ASO injection, and tumor volume was calculated. Three days after the 7^th injection, the mice were sacrificed and tumor tissues were collected and weighed.

Statistical analysis

Statistical analysis was performed using R(3.6.1). A comparison between the FOXD2-AS1 expression in tumor and normal tissues was performed using Wilcoxon rank-sum tests. The samples were divided into high- and low-expression groups (median FOXD2-AS1 expression level as cut-off value). Next, the ROC curve was made to test the performance of FOXD2-AS1 as a diagnostic marker using the survivalROC R package. The relationships between clinicopathologic features and FOXD2-AS1 expression were analyzed using the Wilcoxon rank-sum test and logistic regression. Clinicopathologic characteristics associated with OS and DSS were analyzed with Cox regression and the Kaplan-Meier method. The multivariate Cox analysis was used to identify the influence of FOXD2-AS1 expression on survival, along with other clinical features. A nomogram was constructed based on the results of the multivariate analysis using the rms package in R. A risk score for each patient was calculated as the sum of the scores for each parameter, which were obtained by multiplying the value or level of each parameter and its coefficient. According to the median risk scores, patients were divided into low- and high-risk groups. All hypothetical tests were two-sided, and P values<0.05 were considered significant.