Oncogenic Role of ZFAS1 lncRNA in Head and Neck Squamous Cell Carcinomas

Background: Head and neck squamous cell carcinoma (HNSCC) is a heterogeneous disease with high mortality. The identification of specific HNSCC biomarkers will increase treatment efficacy and limit the toxicity of current therapeutic strategies. Long non-coding RNAs (lncRNAs) are promising biomarkers. Accordingly, here we investigate the biological role of ZFAS1 and its potential as a biomarker in HNSCC. Methods: The expression level of ZFAS1 in HNSCC cell lines was analyzed using qRT-PCR. Based on the HNSCC TCGA data, the ZFAS1 expression profile, clinicopathological features, and expression of correlated genes were analyzed in patient tissue samples. The selected genes were classified according to their biological function using the PANTHER tool. The interaction between lncRNA:miRNA and miRNA:mRNA was tested using available online tools. All statistical analyses were accomplished using GraphPad Prism 5. Results: The expression of ZFAS1 was up-regulated in the metastatic FaDu cell line relative to the less aggressive SCC-25 and SCC-040 and dysplastic DOK cell lines. The TCGA data indicated an up-regulation of ZFAS1 in HNSCCs compared to normal tissue samples. The ZFAS1 levels typically differed depending on the cancer stage and T-stage. Patients with a lower expression of ZFAS1 presented a slightly longer disease-free survival and overall survival. The analysis of genes associated with ZFAS1, as well its targets, indicate that they are linked with crucial cellular processes. In the group of patients with low expression of ZFAS1, we detected the up-regulation of suppressors and down-regulation of genes associated with epithelial-to-mesenchymal transition (EMT) process, metastases, and cancer-initiating cells. Moreover, the negative correlation between ZFAS1 and its host gene, ZNFX1, was observed. The analysis of interactions indicated that ZFAS1 has a binding sequence for miR-150-5p. The expression of ZFAS1 and miR-150-5p is negatively correlated in HNSCC patients. miR-150-5p can regulate the 3′UTR of EIF4E mRNA. In the group of patients with high expression of ZFAS1 and low expression of miR-150-5p, we detected an up-regulation of EIF4E. Conclusions: In HNSCC, ZFAS1 displays oncogenic properties, regulates important processes associated with EMT, cancer-initiating cells, and metastases, and might affect patients’ clinical outcomes. ZFAS1 likely regulates the cell phenotype through miR-150-5p and its downstream targets. Following further validation, ZFAS1 might prove a new and valuable biomarker.

Mikroorganismen und Zellkulturen GmbH, Leibniz Institut, Braunschweig, Germany), respectively. The FaDu cell line was cultured as described previously [32]. All cell lines were cultured with penicillin-streptomycin antibiotic (Merck Millipore, Burlington, MA, USA), and mycoplasma detection tests were performed routinely using the VenorGeM Mycoplasma PCR Detection Kit (Minerva Biolabs, Berlin, Germany).
The spheres forming capacity ability was checked by soft agar assay using low melting temperature SeaPlaque Agarose (Lonza, Basel, Switzerland). The wells of the culture plates were coated with bottom agar (1%), next the single cells (5000 cells/mL) were suspended in 0.3% agarose with optimal culture media, and 1 mL of this mixture onto bottom agar was placed. Cells were incubated under standard conditions and were supplemented with fresh media every 3 days. After 2 weeks, the spheres were measured using a microscope with cellSens Entry software (Olympus, IX70 Fluorescence Microscope, Olympus, Tokyo, Japan).
Total RNA from the cell lines was isolated using a High Pure miRNA isolation kit (Roche, Basel, Switzerland), according to the isolation protocol for total RNA from tissue and cell line samples. Quality and quantity of RNA samples were analyzed using a NanoDrop spectrophotometer (Thermo Scientific, Waltham, MA, USA). cDNA synthesis reactions were performed using 1 µg of RNA and EvoScript Universal cDNA Master (Roche) according to manufacturer s instruction. ZFAS1 (F: 5 -AAGCCACGTGCAGACATC TA-3 and R: 5 -CTACTTCCAACACCCGCATT-3 ) [33] and reference B2M (F: 5 -TTCTGGCCTGGAG GCTATC-3 and R: 5 -TCAGGAAATTTGACTTTCCATTC-3 ) genes were quantified using LightCycler 480 SYBR Green I Master buffer (Roche) and LightCycler 96 (Roche) according to manufacturer's instruction. All data were shown as 2 −∆Ct values and normalized to the B2M. Gene quantification was carried out using three independent cDNA replicates for each of the cell lines.

TCGA Data
The TCGA expression data of lncRNA ZFAS1, expression of selected genes, and clinical data were downloaded from cBioPortal (Head and Neck Squamous Cell Carcinoma, TCGA, Provisional, 530 samples data set) [34], from the UALCAN databases (http://ualcan.path.uab.edu) [35], and from StarBase v3.0 (http://starbase.sysu.edu.cn) [36] for 520 cancers and 44 normal tissue samples. All data is available online, and access is unrestricted and does not require patients consent or other permissions. The use of the data does not violate the rights of any person or any institution.

Gene Analysis
Genes positively and negatively correlated with ZFAS1 (Pearson correlation >+0.3 or <−0.3, respectively) were analyzed using the PANTHER Classification System, classifying them into specific biological processes and cellular pathways [37].
The panel of genes connected with the EMT process and migration, as well as influence on cancer-initiating cells, was created based on previous reports [38][39][40][41][42][43] and analyzed in the ZFAS1 lowand high-expressing groups of patients.

Statistical Analysis
All statistical analyses were performed using GraphPad Prism 5 (GraphPad, San Diego, CA, USA). The Shapiro-Wilk normality test, t-test, and Mann-Whitney U test were used for ZFAS1 and ZNFX1 level (depending on clinical parameters) and gene expressions (depending on ZFAS1 subgroups). The expression level of ZFAS1 and ZNFX1 (depending on the cancer location) was checked using one-way ANOVA obtained using Dunn's multiple comparisons test. All qRT-PCR and TCGA data are presented as mean with SEM. For DSF and OS analyses, the Log-Rank (Mantel-Cox) and Gehan-Breslow-Wilcoxon tests were used, and Hazard Ratio (Mantel-Haenszel; HR) and 95% Confidence Interval (CI) of ratio were calculated. In all analyses, p < 0.05 was used to determine statistical significance.

Availability of Data and Materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request. Raw data are available on the cBioPortal, UALCAN and StarBase v3.0 databases.

ZFAS1 is Up-Regulated in HNSCC Cell Lines and Cancer Samples of HNSCC Patients
The analyzed HNSCC cell lines were characterized by different morphology and tumorigenic potential. The FaDu cells were spindly, more fibroblast-like compared to DOK, SCC-25, and SCC-040, which are scale-like, cube-shaped, epithelial cells ( Figure 1A). Moreover, the FaDu cells were more aggressive and had a higher sphere forming ability (number and size of spheres) compared to the SCC-25 and SCC-040 cell lines (mean sphere diameter: 70.2 µm vs. 35.2 µm vs. 57.2 µm, respectively) and to DOK cell line, which did not form spheres ( Figure 1B).
According to the database (cBioportal and UALCAN), the expression of ZFAS1 was significantly up-regulated in cancer samples of HNSCC patients compared to normal tissue (median expression of 226.109 vs. 175.467 transcripts per million; p = 2.24 × 10 −14 ) (Figure 2A).
HNSCC patients were divided into three main localization groups: oral cavity (n = 314), pharynx (n = 90) and larynx (n = 116), according to the National Institute of Health (NIH) classification, and expression levels of ZFAS1 were analyzed. No differences between tumors from the oral cavity, pharynx, and larynx localizations were observed (p = 0.7093), Figure 2B. HNSCC patients were divided into three main localization groups: oral cavity (n = 314), pharynx (n = 90) and larynx (n = 116), according to the National Institute of Health (NIH) classification, and expression levels of ZFAS1 were analyzed. No differences between tumors from the oral cavity, pharynx, and larynx localizations were observed (p = 0.7093), Figure 2B. (A) Expression in normal (n = 44) and cancer (n = 520) tissues; (B) Expression depending on HNSCC localization (n = 520); Graphs from UALCAN database, modified; Un-paired T-test; the graphs show mean of value presented as transcripts per million; and box and whiskers with 5-95 percentile, oneway ANOVA obtained using Dunn's multiple comparisons tests; ns-no significant, **** p < 0.0001.  HNSCC patients were divided into three main localization groups: oral cavity (n = 314), pharynx (n = 90) and larynx (n = 116), according to the National Institute of Health (NIH) classification, and expression levels of ZFAS1 were analyzed. No differences between tumors from the oral cavity, pharynx, and larynx localizations were observed (p = 0.7093), Figure 2B. and cancer (n = 520) tissues; (B) Expression depending on HNSCC localization (n = 520); Graphs from UALCAN database, modified; Un-paired T-test; the graphs show mean of value presented as transcripts per million; and box and whiskers with 5-95 percentile, oneway ANOVA obtained using Dunn's multiple comparisons tests; ns-no significant, **** p < 0.0001. and cancer (n = 520) tissues; (B) Expression depending on HNSCC localization (n = 520); Graphs from UALCAN database, modified; Un-paired T-test; the graphs show mean of value presented as transcripts per million; and box and whiskers with 5-95 percentile, one-way ANOVA obtained using Dunn's multiple comparisons tests; ns-no significant, **** p < 0.0001.

ZFAS1 Levels Differ Depending on Clinicopathological Parameters
The expression levels of ZFAS1 were analyzed depending on the group division based on available clinicopathological parameters in all HNSCC samples. The significant differences between expression levels of ZFAS1 were observed in patients with various cancer stage (p = 0.0091) and T-stage (p = 0.0169). Other analyzed parameters did not differ between the studied groups (Table 1).

Association of ZFAS1 Expression and DFS and OS in the Studied Patients
HNSCC samples were divided into low, medium, and high ZFAS1 expression groups using the <25, 25-75 and >75 percentile of ZFAS1 expression as a cutoff, respectively. We observed a slightly longer DFS of low ZFAS1 expression patients compared to the high expression group (p = 0.0598; HR = 0.6554; 95% CI = 0.4029-1.066). We also detected a slight longer OS in the low ZFAS1 expression group compared to the high group (p = 0.0356; HR = 0.6922; 95% CI = 0.4623-1.037) ( Figure 3).

ZFAS1 is Involved in Important Cellular Processes
Next, genes positively and negatively correlated with ZFAS1 expression were analyzed. Fourhundred-forty-one genes were positively, and 112 genes were negatively correlated with the studied lncRNA (Pearson correlation >+0.3 or <−0.3, respectively). The classification analysis revealed that the genes positively correlated with ZFAS1 genes are associated with the regulation of multiple cellular

ZFAS1 is Involved in Important Cellular Processes
Next, genes positively and negatively correlated with ZFAS1 expression were analyzed. Four-hundred-forty-one genes were positively, and 112 genes were negatively correlated with the studied lncRNA (Pearson correlation >+0.3 or <−0.3, respectively). The classification analysis revealed that the genes positively correlated with ZFAS1 genes are associated with the regulation of multiple cellular processes and pathways, such as cell cycle, cell adhesion, signal transduction, death, response to stimulus, apoptosis signaling pathway, FAS signaling pathway, integrin signaling pathway, and mRNA splicing. The genes negatively correlated with ZFAS1 are associated with processes such as cell adhesion, signal transduction, cell differentiation, death, response to stimulus, angiogenesis, oxidative stress response, and various pathways (apoptosis, cadherin and integrin signaling pathways, EGFR, endothelial, FAS, FGF, insulin/IGF, TGF-beta, VEGF, interleukin, JAK/STAT, PDGF, PI3K, p53, p38, Ras, Toll receptor, and Wnt signaling pathways) ( Table 2).

lncRNA ZFAS1 is Negatively Correlated with ZNFX1 mRNA in HNSCC
A previous report has indicated that lncRNA ZFAS1 shares the same transcription start sites with ZNFX1 (Zinc Finger NFX1-Type Containing 1) gene, and that expression of ZFAS1 and ZNFX1 are positively correlated [10]. Surprisingly, using the StarBase v3.0 database, the negative correlation between ZFAS1 and ZNFX1 in HNSCC patients was observed (r = −0.308, p = 1.75 × 10 −12 ) ( Figure 4A).
Next, the expression levels of ZNFX1 were checked depending on cancer localization. No differences between tumors from the pharynx (−0.2584 ± 0.08905) and larynx (−0.2309 ± 0.07557) localizations were observed (p > 0.9999), but in the case of oral cavity significantly up-regulation of ZNFX1 compared to pharynx or larynx was observed (p < 0.0001), Figure 4D.
HNSCC patients were divided into low, medium, and high ZNFX1 expression groups and DSF as well as OS were analyzed. No differences between groups of patients in the case of DFS and OS were observed (p > 0.05) ( Figure 4E). The expression levels of ZNFX1 were also analyzed depending on the group division based on available clinicopathological parameters in all HNSCC samples. The significant differences between expression levels of ZNFX1 were observed in the case of gender (p = 0.0004), cancer stage (p < 0.0001) and T-stage (p = 0.0240), cancer grade (p = 0.0158), perineural invasion (p = 0.0022) or HPV status (p = 0.0086). Other analyzed parameters did not differ between the studied groups (Table 3).  Next, the expression levels of ZNFX1 were checked depending on cancer localization. No differences between tumors from the pharynx (−0.2584 ± 0.08905) and larynx (−0.2309 ± 0.07557) localizations were observed (p > 0.9999), but in the case of oral cavity significantly up-regulation of ZNFX1 compared to pharynx or larynx was observed (p < 0.0001), Figure 4D.
HNSCC patients were divided into low, medium, and high ZNFX1 expression groups and DSF as well as OS were analyzed. No differences between groups of patients in the case of DFS and OS were observed (p > 0.05) ( Figure 4E). The expression levels of ZNFX1 were also analyzed depending on the group division based on available clinicopathological parameters in all HNSCC samples. The significant differences between expression levels of ZNFX1 were observed in the case of gender (p = 0.0004), cancer stage (p < 0.0001) and T-stage (p = 0.0240), cancer grade (p = 0.0158), perineural invasion (p = 0.0022) or HPV status (p = 0.0086). Other analyzed parameters did not differ between the studied groups ( Table 3).

Role of ZFAS1 in the EMT Process, Cancer-Initiating Cells Maintenance, and Metastasis Process in HNSCC
ZFAS1 is described as a modulator of the EMT process, cancer-initiating cell maintenance, and metastasis in many cancers [47], so its role in HNSCC was also checked.

Discussion
The major finding of the study is a delineation of the biological role of lncRNA ZFAS1 and its potential utility as a biomarker in HNSCC. We report the up-regulation of ZFAS1 in HNSCC cell lines and cancer tissue samples derived from patients. Moreover, compared to SCC-25 and SCC-040 or DOK cell lines, higher levels of ZFAS1 are observed in the FaDu cell line, which is highly tumorigenic and possesses fibroblast-like features. Interestingly, the ZFAS1 expression level did not differ in various HNSCC localizations.
Higher expression of ZFAS1 was found in HNSCC patients with more advanced disease. Moreover, patients with a lower level of ZFAS1 displayed slight longer DFS and OS compared to the high-expressing group. Gao et al. presented a similar observation, where higher ZFAS1 expression was significantly correlated with advanced tumor stage and worse OS in glioma patients [11]. In the case of skin melanoma, higher ZFAS1 expression was associated with higher clinical stage, primary tumor thickness, and with the presence of lymph node metastases. Also, it served as a predictive marker of DFS and OS [25]. Shi et al. based on the retrospective analysis of 398 lymph node-negative esophageal squamous cell carcinoma patients, reported an association of higher ZFAS1 expression with less differentiated cancers [28].
The analysis of genes positively and negatively correlated with ZFAS1 in HNSCC indicated their association with some important cellular processes. Genes positively correlated with ZFAS1 were associated with cell cycle, cell adhesion, signal transduction, death, response to stimulus, apoptosis signaling pathway, FAS signaling pathway, integrin signaling pathway, and mRNA splicing. The genes negatively correlated with ZFAS1 were associated with processes such as: cell adhesion, signal transduction, cell differentiation, death, response to stimulus, angiogenesis, oxidative stress response, and multiple pathways (apoptosis, cadherin and integrin signaling pathways, EGFR, endothelial, FAS, FGF, insulin/IGF, TGF-beta, VEGF, interleukin, JAK/STAT, PDGF, PI3K, p53, p38, Ras, Toll receptor, and Wnt signaling).
Askarian-Amiri et al. described that the lncRNA ZFAS1 and ZNFX1 (Zinc Finger NFX1-Type Containing 1) genes share the same transcription start sites, and that expression of ZFAS1 and ZNFX1 are positively correlated [10]. Surprisingly, our analysis did not confirm the above observation in HNSCC, and ZNFX1 was negatively correlated with ZFAS1. However, ZNFX1 is up-regulated in cancer compared to normal samples and its expression depends on cancer localization. Our analysis also indicated, that expression of ZNFX1 depends on clinicopathological parameters and is up-regulated in the case of: female patients, lower cancer stage, T-stage and cancer grade, it is associated with cancer invasion to the space surrounding the nerves, and it is higher in HPV negative patients. Moreover, no difference between ZNFX1 level and patients' survival (DFS neither OS) was observed. Unfortunately, there is lack of reports indicated the role of ZNFX1 in HNSCC or other cancers.
Previous studies have indicated the role of ZFAS1 in the regulation of EMT process, migration, and influence on cancer-initiating cells in different cancer types [11,12,14,15]. In our study, we also analyzed the panel of target genes studied in previous reports and associated with these processes [38][39][40][41][42][43]. We found genes up-regulated in the ZFAS1 low expression group of patients compared to the ZFAS1 high group. These displayed a suppressor function for EMT processes, metastases, and cancer initiating cells maintenance, and down-regulation of the genes supporting these processes. These data support the hypothesis that ZFAS1 is an oncogene and its high expression is associated with the more aggressive phenotype of HNSCC. It has been proposed that ZFAS1 is a key activator of the EMT process in glioma, colorectal cancer, and gastric cancer [12,15]. However, the authors analyzed only a limited number of markers associated with the EMT process. Our analysis was based on multiple marker genes, which sometimes display an opposite function to ZFAS1 in these processes. Examples include patients with high level of ZFAS1 with low expression of NOTCH1, one of the important elements of the pathway described in the context of EMT and cancer-initiating cells [38]. Gao et al. showed that ZFAS1 affects the NOTCH signaling pathway. Knockdown of ZFAS1 caused down-regulation of the HES-1 (HES family bHLH transcription factor 1) and NICD (Notch intracellular domain), which are NOTCH signal-related proteins, but the mechanism of NOTCH signaling regulation by ZFAS1 remains unknown [11]. However, in HNSCC, a high level of NOTCH1 was associated with better survival [48], which supports our findings, where low ZFAS1 expressing patients displayed higher NOTCH1 level and better survival.
Moreover, we observed a high expression level of EGFR and CD274 (PD-L1) in the group of ZFAS1 low-expressing patients. EGFR and PD-L1 are well-known targets for immunotherapy in HNSCC patients [4]. Accordingly, the patients with low expression of ZFAS1 might benefit from anti-EGFR (e.g., cetuximab) and anti-PDL1 (e.g., atezolizumab) therapy.
The direct regulation mechanism of mRNAs by lncRNA ZFAS1 remains unknown. However, some previous reports have indicated that ZFAS1 can act as a molecular sponge and reduce the abundance of miRNAs, such as miR-9, miR-150, miR-484 or miR-200b/c, and reduce their activity in the cell and have an indirect influence on mRNAs [47]. We analyzed this possible mechanism and indicated that, indeed, a sequence of ZFAS1 possesses the binding site for miR-150-5p. Moreover, the negative correlation between ZFAS1 and miR-150-5p was observed in HNSCC patients. Next, we checked if, in the group of genes associated with ZFAS1, any targets for miR-150-5p are present. We found nine potential mRNAs targets, CPEB4, GAB1, ARHGEF10L, KALRN, DSP, IL13RA1, EIF4E, UQCR11, and MET. Only UQCR11 and EIF4E were significantly down-regulated in the group of patients with high level of miR-150-5p and low ZFAS1, which supports our assumption of direct regulation by miR-150-5p. There is no association between the UQCR11 (ubiquinol-cytochrome c reductase, complex III sub-unit XI) gene and cancer. However, the second gene, eukaryotic translation initiation factor 4E (EIF4E), is activated in cancers [49] and is required for translation of some mRNAs involved in proliferation and survival [50], as well as in EMT process and cancer invasion [51,52]. The phosphorylation of EIF4E is very frequently observed in HNSCC [49]. Moreover, EIF4E is up-regulated in surgical margins of HNSCC patients with local recurrence and could serve as a prognostic biomarker [53]. DeFatta et al. indicated that the FaDu cell line displays a high EIF4E protein level and that this is similar to the level observed in patients. Knock-down of EIF4E results in suppression of the tumorigenic and angiogenic properties of the FaDu cell line manifested by loss of capacity to grow in soft agar, reduced expression of angiogenic factors (FGF-2 and VGF), and loss of tumor growth in nude mice [54]. The FaDu cell line has the highest ZFAS1 level among HNSCC cell lines, and our results suggest that ZFAS1 reduced the level of suppressor miR-150-5p and maintained a high level of EIF4E. It seems likely that the oncogenic EIF4E, in turn, up-regulates expression of some genes associated with EMT metastasis and could result in poor patient outcome ( Figure 6). However, the above hypothesis needs to be further verified by in vitro and in vivo analysis of ZFAS1 function in HNSCC.

Conflicts of Interest:
The authors declare that there is no conflict of interest regarding the publication of this paper. References:

Conflicts of Interest:
The authors declare that there is no conflict of interest regarding the publication of this paper.