Identification and in vitro validation of prognostic lncRNA signature in head and neck squamous cell carcinoma

ABSTRACT Long non-coding RNAs (lncRNAs) are promising cancer prognostic markers. However, the clinical significance of lncRNA signatures in evaluating overall survival (OS) outcomes of head and neck squamous cell carcinoma (HNSCC) has not been explored. This study aimed to assess the significance of lncRNA in HNSCC and to develop a lncRNA signature related to OS in HNSCC. LncRNA expression matrices were retrieved from the Cancer Genome Atlas (TCGA) data. Least Absolute Shrinkage and Selection of the Operator (LASSO), univariate and multivariate Cox regression were used for establishing a prognostic model. In vitro experiments were carried out to demonstrate the biological role of lncRNA. A prognosis model based on 7 DElncRNAs was finally established.The patients were then divided into high-risk and low-risk groups. Relative to the low-risk group, overall survival times for patients in the high-risk group were significantly low (P=2.466e−07). Risk score remained an independent prognostic factor in univariate (HR=1.329, 95%CI=1.239−1.425, p < 0.001) and multivariate (HR=1.279, 95%CI=1.184−1.382, p < 0.001) Cox regression analyses. The area under the curve (AUC) of the signature was as high as 0.78. Expressions of FOXD2-AS1 in tumor tissues were elevated, and significantly correlated with OS (P=0.008). FOXD2-AS1 silencing then significantly reduced HNSCC cell proliferation, invasion, and migration. In conclusion, a lncRNA signature was established for HNSCC prognostic prediction and FOXD2-AS1 was identified as an HNSCC oncogene.


Introduction
Head and neck squamous cell carcinomas (HNSCC) are a cluster of cancers that occur in the mouth, oropharynx, and laryngeal regions [1]. HNSCC is the sixth most common malignancy in the world, with over 600,000 new cases diagnosed each year and 380,000 deaths worldwide [2,3]. The current treatments for HNSCC include surgery, chemotherapy, radiotherapy, and targeted therapy. Despite advancements in therapeutic technology, HNSCC is a particularly aggressive solid tumor. Local recurrence and distal metastasis are still the two major challenges of HNSCC treatments, and the current 5-year overall survival rate of HNSCC patients is about 50%, while that of advanced patients is less than 50% [4,5]. However, there is still a lack of effective clinical prognosis prediction methods for HNSCC patients [6]. Therefore, there is an urgent need to find effective markers to predict the prognosis of HNSCC, which is necessary for guiding personalized anti-cancer treatment strategies.
Long non-coding RNAs (lncRNAs), which cannot encode proteins, are RNA molecules whose transcriptional length is more than 200 nucleotides [7]. LncRNAs have been shown to regulate gene expression in a variety of ways at various levels, including epigenetic, transcriptional, post-transcriptional, translational, and post-translational, by interacting with mRNA, DNA, proteins, and miRNAs [8][9][10]. LncRNAs play vital roles in the pathogenesis of several cancers, including HNSCC [11]. Zheng et al found that long noncoding RNA01134 can interact with microRNA-4784 and down-regulate the expression of structure-specific recognition protein 1, thereby promoting the progression of hepatocellular carcinoma [12]. In HNSCC, Long intergenic noncoding RNA HOTAIR has been reported to be overexpressed and regulate PTEN methylation, which is involved in tumor development [13]. In addition to solid tumors, lncRNA is also considered to be related to the occurrence and development of hematological tumors. In acute myeloid leukemia, Long noncoding RNA SATB1-AS1 participates in the tolerance process of chemotherapy [14]. The expression of lncRNA varies greatly in different diseases and tissues and different disease stages. Therefore, compared with other tumor characteristics, lncRNA is a potential marker for diagnosis and prognosis [15,16]. In the past, although studies have reported that Long noncoding RNA is related to the pathogenesis and prognosis of HNSCC [6,13], there is a lack of systematic analysis of the prognosis of HNSCC and lncRNA at the omics level. In this study, we hypothesized that there is a lncRNA expression signature that is significantly related to the prognosis of HNSCC, which is promising for clinical application. We developed a lncRNA-based signature based on a large number of HNSCC patient-related data in the database and comprehensively analyzed genomic data through bioinformation methods to predict the prognosis of HNSCC, and verified it in experiments in vitro.

Data acquisition and processing
RNA-Seq and clinical data for HNSCC patients were retrieved from the TCGA database (https:// portal.gdc.cancer.gov/). These data were extracted and integrated by downloading Perl (https://www. perl.org/get.html). RNA data were normalized using the R software (v.4.0.2 http://www.r-pro ject.org). In total, 498 HNSCC patients were randomly assigned to the training and validation groups in a 6:4 ratio and analyzed using the R software package 'caret.' Differences between RNA expression patterns were analyzed using the R software 'edgeR' package to determine DElncRNAs. Threshold was set at | log 2 foldchange (FC) >2.0, using adjusted P < 0.01.

Establishment of DElncRNA-associated prognostic model
Screening of DElncRNAs with significant associations with overall survival (OS) was carried out with Univariate Cox regression analysis [17]. Validation of OS-related DElncRNAs was performed by LASSO regression using the R software 'glmnet' package. Then, Multivariate COX regression analyses were conducted on candidate lncRNAs obtained from LASSO regression, and a model was finally established. Risk scores for every patient were determined by the linear combination of lncRNA expression levels (x) and regression coefficient (β) from the multivariate COX regression analysis. Risk score = β 1 x 1 + β 2 x 2 + · ···· +β n x n .

Confirmation of lncRNA signature
The participants were divided into high-and lowrisk groups based on their median risk scores, and Kaplan Meier survival curves were drawn [18]. The univariate Cox proportional risk regression model was used to analyze data. The prognostic values were determined using an independent receiver operating characteristic (ROC) curve. Specificity, the sensitivity of risk score in predicting survival was compared to determine prognostic value. In addition, the associations between lncRNA risk score prediction and other clinical factors were verified by multivariate Cox regression. The model's predictive accuracy was validated using the test group (n = 199).

Extraction of total RNA and qRT-PCR analysis
The RNA extraction from cells was carried out with the Mirneasy Mini Kit extraction reagent box (QIAGEN, Germany), after which cDNA synthesis was carried out using the PrimeScript™ RT reagent Kit (TaKaRa, Japan). The qRT-PCR reaction was undertaken on SYBR® Premix Ex Taq™ (TaKaRa, Japan) using various primers: FOXD2-AS1, forward: TGGACCTAGCTGCAGCTCCA, reverse: AGTT GAAGGTGCACACACTG, and GAPDH, forward: CCCATCACCATCTTCCAGGAG, reverse: GTTG TCATGGATGACCTTGGC.

Fluorescence in situ hybridization (FISH)
The FOXD2-AS1-FISH probe was designed and synthesized by Gemma Gene (Shanghai, China). About 1 × 10 4 cells were inoculated overnight on a 48-well plate, fixed with 4% paraformaldehyde for 15 minutes at room temperature, aspirated the paraformaldehyde, and permeated with 0.1% Buffer. A (main component triton x 100) for 15 minutes, and incubate at 37°C for 30 minutes. Added the FISH probe to the hybridization mixture and incubated at 37°C overnight (12-18 hours). Then, the cells were washed with 20x SSC buffer. DAPI was stained to observe cell nuclei. Finally, the cells were observed under a Leica fluorescence microscope (DMi8, Leica, Germany).

Wound healing and Transwell migration assays
The wound healing and Transwell migration assays were carried out following the methods reported in a previously published study [19]. In short, Cal-27 and FADU cells were seeded in 6-well plates containing DMEM without FBS and scraped using 200-µl pipette tips. Cell migration and imaging were assessed at 0and 24-hours following damage using inverted microscopy (OLYMPUSIX73, Olympus, Tokyo, Japan). Invasiveness of cells was detected using Transwell migration chamber (Corningcostar; Aperture 8 μm) coated with Matrigel (Sigma). Transfected cells were added onto the top compartment with a resuspended medium without FBS, and the bottom compartment was filled with 30% FBS. Uninvaded cells were removed from the top of the membrane after 48 hours, and invaded cells were then fixed at the bottom of the membrane and stained. Cells were counted from 5 random fields (100x magnification) and photographs were taken under a light microscope.

Results
To identify the prognostic lncRNA signature in head and neck squamous cell carcinoma, the lncRNA expression was analyzed in HNSCC tumor and non-tumor patient samples, and specific lncRNAs correlated with overall survival outcomes of HNSCC patients were involved. Then, the prognostic lncRNA signature including 7 lncRNAs was identified by multiple regression followed by its validation in the lncRNA model. Finally, we chose FOXD2-AS1 in the 7 lncRNAs and validated the cancer-promoting function in vitro.

Establishment of the lncRNA signature
In this study, 72 DElncRNAs were found to have a significant correlation with overall survival outcomes (P < 0.01) (Supplementary Table 1). DElncRNAs were further verified by Lasso regression and 10 DElncRNAs were obtained for the model to achieve the best performance ( Figure S2A, B). The 10 DElncRNAs were then analyzed using multiple regression. Findings showed that 7 HNSCC prognosis-related DElncRNAs including HOXC-AS2, FOXD2-AS1, LINC00460, AC022031.2, AP002957.1, AP002478.1, and AC011369.1 were obtained. The Forest plot for the relationship between each lncRNA and OS is shown in Figure 2a. Each lncRNA coefficient is shown in Supplementary Tables 2. The expression of lncRNAs was considerably upregulated in HNSCC, according to TCGA transcriptome data ( Figure 2b).

Validation of the lncRNA signature
The risk scores for each patient were computed using the expression levels and risk coefficients of seven prognostic lncRNAs in the model. Patients' allocation into low-and high-risk groups was done using the median risk score. Expression levels of 7 lncRNAs were distributed along with the training group as shown in Figure S3. The risk scores and the survival time distribution for the training group are shown in Figure 3a and b, respectively. Kaplan-Meier curves of low-as well as high-risk groups, are shown in Figure 3c. The high-risk patients were found to have shorter survival times (P = 2.466e−07) compared with low-risk patients. In addition, the prognostic model based on lncRNA biomarkers had an AUC of 0.78, which was higher than other indicators including age, gender, and tumor stage (Figure 3d). Furthermore, univariate analysis revealed that the risk score had a significant association with OS (HR = 1.329, 95%CI = 1.239 − 1.425, P < 0.001) (Figure 3e). In the multivariate analysis, the risk score was shown to be an independent prognostic factor in HNSCC (HR = 1.279, 95% CI = 1.184 − 1.382, P< 0.001) (figure 3f). Therefore, the lncRNA signature had better HNSCC prognostic abilities, when compared with the other indicators.

Validation and verification of the lncRNA model
Distributions of lncRNA expressions, risk scores, as well as survival times in the validation group are shown in Figure S4and

Expression levels of FOXD2-AS1 in HNSCC tissues are significantly elevated and correlated with poor prognostic outcomes
The expression of the LncRNA FOXD2-AS1 is abnormal in a range of tumors, and it is associated with patient survival and prognosis [13]. Therefore, the role of FOXD2-AS1 in HNCSS needs to be investigated further. The current study evaluated the transcription level of FOXD2-AS1 in HNSCC by evaluating the TCGARNA-seq data. In tumor tissues, FOXD2-AS1 expressions were significantly elevated, compared to normal tissues (Figure 5a). In addition, we analyzed FOXD2-AS1 expressions in tumors and adjacent normal tissues by TCGA data (Figure 5b). Patients with low levels of FOXD2-AS1 expression had considerably longer survival times than those with elevated levels of expression ( Figure 5c). Furthermore, FOXD2-AS1 was found to be an independent predictor in HNSCC patients (Figure 5d,e).

FOXD2-AS1 functional analysis and its relationship with HNSCC clinicopathological characteristics
GSEA analysis was undertaken in the current study to determine the biological effects of FOXD2-AS1. High scores of FOXD2-AS1 were shown to be significantly enriched in pathways such as cell cycle, base excision repair, DNA replication, and the p53 signaling pathway (figure 5f). The relationship between lncRNA expression levels and clinico-pathological features was determined using TCGA data. Expressions of FOXD2-AS1 were significantly correlated with grade (P = 0.006), T stage (P = 0.027), N stage (P = 0.002), and TNM stage (P = 0.004) (Figure 6a-d).

FOXD2-AS1 was amplified in HNSCC, and suppressing FOXD2-AS1 inhibited cell proliferation and migration in HNSCC cells
FOXD2-AS1 expression levels were detected in HNSCC cell lines (CAL27, FADU, and TSCCA) and normal HOK cell lines. In the HNSCC cell lines, FOXD2-AS1 was highly expressed when compared to the normal HOK cell lines, whereas expression levels of CAL27 and FADU were higher compared with that of TSCCA. CAL27 and FADU cell lines were, therefore, selected for further experiments (Figure 7a). The current study evaluated cell transfection efficiency using qRT-PCR and established that FOXD2-AS1 expression levels were significantly reduced after shRNA transfection ( Figure 7b). Meanwhile, the expression and localization of FOXD2-AS1 have been verified by FISH experiments in cell lines including HOK, TSCCA, CAL-27, FADU cell lines. The results showed that FOXD2-AS1 expression levels were higher in HNSCC cell lines (TSCCA, CAL-27, FADU) compared with normal HOK cell lines. Moreover, FOXD2-AS1 was mainly located in the nucleus of these cell lines, as shown in Figure 7c. CCK-8 was used to establish the effects of FOXD2-AS1 knockdown on cell proliferation. When FOXD2-AS1 was silenced, CAL27 and FADU cell proliferation was significantly reduced compared to sh-NC cell proliferation (Figure 8a,b). Furthermore, wound healing tests showed that FOXD2-AS1 silencing significantly inhibited wound healing in the two cell lines (Figure 8c). The invasion and migration rates of FADU cells and CAL27 transfected with shRNA were significantly suppressed as compared with those transfected with sh-NC (Figure 8d,e).

Discussion
According to previous studies, HNSCC patients are prone to recurrence and metastasis, indicating the urgent need for aggressive research into novel biomarkers and biomarker-based therapies to improve survival [20]. In recent years, lncRNAs have received a lot of attention as emerging regulators with a variety of biological functions [21]. lncRNAs are transcripts with over 200 nucleotides and are involved in several cellular properties, including cell growth, pluripotency, and cell differentiation [22]. Several studies have suggested that dysregulated lncRNAs could influence cancer progression [23,24]. Many lncRNAs are potential biomarkers and targets for cancer diagnosis and treatment [25]. Several prognostic models based on lncRNA have been developed for tumors [17,18,26]. Previous studies have also reported the association between Long noncoding RNA and the pathogenesis and prognosis of HNSCC [6,13], but there is a lack of systematic analysis of Long noncoding RNA at the omics level in HNSCC. The current study used bioinformatics and univariate Cox analysis to evaluate OSrelated DElncRNAs using HNSCC transcriptome data from the TCGA database. A total of 72 OSrelated DElncRNAs were identified. Lasso regression analysis was used to modify the list of candidate DElncRNAs. Multivariate Cox regression analysis showed prognostic signatures including 7 DElncRNAs. Kaplan-Meier, Cox regression analyses as well as time-dependent ROC curves were then used to verify the prognostic value of the lncRNA signature, which were considered independent predictors of HNSCC prognosis. Further validation in the internal validation set was undertaken in the current study.
Previous studies have reported that FOXD2-AS1 is associated with tumorigenesis and development. According to studies, FOXD2-AS1 is upregulated in a variety of cancers, including gastric cancer, bladder cancer, nasopharyngeal cancer, lung cancer, esophageal cancer, liver cancer, thyroid cancer, colorectal cancer, and skin cancer, and promotes cancer cell proliferation, migration, and invasion, and is negatively correlated with prognosis [27]. The mechanism by which FOXD2-AS1 promotes tumor progression has been reported differently in different tumors. FOXD2-AS1 can regulate microRNA; Ni et al. [28] reported that overexpression of FOXD2-AS1 promotes glioma growth via modulating the miR-185-5P/HMGA2 axis and PI3K/AKT signaling pathways. In rectal cancer, up-regulated FOXD2-AS1 has been reported to regulate miR-25-3p/Sema4c Axis to promote the invasion and metastasis of colorectal cancer cells [29]; in addition, FOXD2-AS1 is abnormally elevated in breast cancer tissues and promotes FOXD2-AS1/miR-150-5p/PFN2 axis, accelerating the growth of transplanted tumors in vivo [30]. Conversely, knock-down of FOXD2-AS1 can inhibit tumor progression. Knockdown of FOXD2-AS1 can inhibit the progression of gallbladder cancer by regulating the methylation of MLH1 [31]. FOXD2-AS1 is highly expressed in proliferative hemangioma tissues, and its knockdown inhibits cell proliferation, migration, invasion, and colony formation through the miR-324-3p/PDRG1 pathway [32]. FOXD2-AS1 depletion can also regulate the expression of CDX1 and inhibit the proliferation of cervical cancer cells [33]. Therefore, FOXD2-AS1 induces cancerpromoting activities, which are mostly thought to be due to its interaction with microRNAs, thereby altering gene expression or epigenetic modification. The mechanism may be different in different cancers. But the mechanism also may be similar, and FOXD2-AS1 changes the entire gene expression network. In this study, a single-gene analysis of FOXD2-AS1 was undertaken to evaluate the role of FOXD2-AS1 in HNSCC. We compared FOXD2-AS1 expression levels in the tumor and normal or adjacent tissues in the TCGA-HNSCC cohort. Both unpaired and paired tests indicated increased expression levels of FOXD2-AS1 in HNSCC. FOXD2-AS1 expression levels were significantly correlated with clinicopathological characteristics, including grade, T stage, TNM stage, and N stage. Patients with elevated FOXD2-AS1 expression levels had worse prognostic outcomes compared with those with low FOXD2-AS1 expression levels. We found that FOXD2-AS1 independently predicted tumor prognosis. Furthermore, GSEA analysis of FOXD2-AS1 expression data in TCGA-HNSCC was performed. The findings revealed that FOXD2-AS1 influences the occurrence and progression of HNSCC through 9 pathways. Some of these biological effects of related pathways on tumors have been reported in the literature, such as cell cycle regulation. Several post-translational modifications control cell cycle progression, and cyclin-dependent kinases (CDKs) are key cell cycle regulators [34]. Zhang et al. [35] reported that lncRNA CASC11 enhanced gastric cancer cell proliferation, migration as well as invasion by regulating the expression of cell cycle-related protein CDK1. Studies by Qin et al. [36] have shown that microRNA-99a-5p inhibits breast cancer progression as well as cell cycle pathway by down-regulating cell cyclerelated protein CDC25A. Furthermore, inhibiting CDK4/6 increases the radiosensitivity of HPVnegative HNSCC [37]. In addition, the role of FOXD2-AS1 in HNSCC was also related to the p53 signaling pathway, and the p53 signaling pathway plays a part in the regulation of expression of several genes and is closely related to the development of human tumors, whose activation initiates cell cycle arrest as well as DNA repair, and inhibits proliferation of OSCC cells [38]. Overexpression of hsa-miR-128a-mediated p53 signaling pathway increases the role of Pembrolizumab in Laryngeal Squamous Cell Carcinoma Cell lines [39]. The findings of the current study confirmed that FOXD2-AS1 is elevated in HNSCC tissues and is a potential predictor of HNSCC in patients, thereby promoting the proliferation of HNSCC cell lines. The current study's findings also revealed that FOXD2-AS1 has an impact on HNSCC biology via these pathways. However, this potential mechanism should be explored in future research.
There were certain limitations to the current study, which have been noted. Firstly, the study lacked a significant number of clinical specimens to detect differential expression of lncRNAs, and only a database was used for paired analysis. Second, no further in vivo or in vitro research has been done to confirm the prognostic value of the proposed lncRNA signatures. Instead, bioinformatics was used to infer from online data sets. Finally, lncRNAs lacked regulatory mechanisms. The current study intends to refine these findings in the future.

Conclusion
In conclusion, we identified novel lncRNA signatures including HOXC-AS2, FOXD2-AS1, LINC00460, AC022031.2, AP002957.1, AP002478.1, and AC011369.1 that could predict HNSCC prognosis independently, and we furtherly validated the cancerpromoting function of FOXD2-AS1 in HNSCC in vitro, which was similar to that in other kinds of cancers. Our findings may contribute to the creation of prognostic prediction systems for HNSCC in the clinic. However, the current study's findings need to be validated, and the functions of another six lncRNAs in HNSCC advancement have yet to be investigated.