FOXP4‐AS1 is a favorable prognostic-related enhancer RNA in ovarian cancer

Abstract Ovarian cancer (OV) is the main cause of deaths worldwide in female reproductive system malignancies. Enhancer RNAs (eRNAs) are derived from the transcription of enhancers and has attracted increasing attention in cancers recently. However, the biological functions and clinical significance of eRNAs in OV have not been well described presently. We used an integrated data analysis to identify prognostic-related eRNAs in OV. Tissue-specific enhancer-derived RNAs and their regulating genes were considered as putative eRNA–target pairs using the computational pipeline PreSTIGE. Gene expression profiles and clinical data of OV and 32 other cancer types were obtained from the UCSC Xena platform. Altogether, 71 eRNAs candidates showed significant correlation with overall survival (OS) of OV samples (Kaplan–Meier log-rank test, P<0.05). Among which, 23 were determined to be correlated with their potential target genes (Spearman’s r > 0.3, P<0.001). It was found that among the 23 prognostic-related eRNAs, the expression of forkhead box P4 antisense RNA 1 (FOXP4-AS1) had the highest positive correlation with its predicted target gene FOXP4 (Spearman’s r = 0.61). Moreover, the results were further validated by RT-qPCR analysis in an independent OV cohort. Our results suggested the eRNA FOXP4-AS1 expression index may be a favorable independent prognostic biomarker candidate in OV.


Introduction
Ovarian cancer (OV) is the fifth malignant reproductive tumor in women leading to mortality [1]. The Surveillance, Epidemiology and End Results (SEER) Program reported United States would have 21750 new diagnosed cases and 13940 new deaths in women would be seen in 2020 [2]. During the past 30 years, the survival in OV has barely improved though chemotherapy drugs and surgical approaches had obtained impressive advances, with 5-year survival remaining approx. 30% for patients in advanced stages [3]. Thus, it remains an urgent issue to explore underlying molecular mechanism for OV progression.
Enhancer RNA (eRNA) is considered to be a type of long noncoding RNAs (lncRNAs) transcribed from putative enhancer regions [4]. However, the functions of eRNAs remain enigmatic presently. It was suggested that enhancer transcription is the noisy byproduct of transcription machinery [5]. However, more and more studies have demonstrated the importance of eRNAs in transcriptional machinery to mediate the target genes transcription [6,7]. Moreover, some studies showed that knockdown of the eRNA was associated with the down-regulation of target genes [8][9][10]. In human cancers, the enhancer activation and production of eRNAs were highly associated with the dysregulation of tumor oncogenes, tumor suppressor genes and other stimuli [11,12]. Current models suggested that eRNAs interact with RNA polymerase II (RNA pol II), mediators and transcription factors to promote promoter-enhancer looping and the consequent up-regulation of the corresponding target genes [13]. Research by Li et al. [14] pointed out that widely estrogen-induced eRNA transcription was associated with up-regulated corresponding genes in hormone-dependent tumor breast cancer cells. Zhao et al. [15] found PSA eRNA can cis and trans regulate expression of a subclass of genes involved in androgen action and cancer progression in castration-resistant prostate cancer cells. Some eRNAs are related to clinical features of cancer or even patients survival [16]. Zhao et al. [17] identified strong correlation of eRNA expression and smoking history (SCRIBe), grade (APH1Ae), stage (CELF2e), subtype (EN1e) and survival (TAOK1e and NET1e). However, the biological functions and clinical significance of eRNAs in OV have not been well described presently.
Above all, our study aims to explore certain prognostic-related eRNAs in OV patients based on bioinformatics analysis. A total of 71 eRNAs were significantly correlated with overall survival (OS) of OV patients from The Cancer Genome Atlas (TCGA). Among these 71 eRNAs, 23 were determined to be correlated with their potential target genes. The functional lncRNA forkhead box P4 antisense RNA 1 (FOXP4-AS1) located within the tissue-specific enhancer and its nearby gene forkhead box P4 (FOXP4) showed the highest correlation using PreSTIGE algorithm and UCSC Xena browser (http://xenabrowser.net). Furthermore, we also validated the effect of FOXP4-AS1 on clinical outcome as an important prognostic-related eRNA in an independent OV cohort.

Acquisition of eRNAs data
PreSTIGE was able to predict the enhancers through identifying protein-coding genes with increased tissue-specific expression. PreSTIGE was a method based on the assumption that these genes are targets of tissue-specific enhancers within the specified domain size [18]. In this current study, a slightly modified version of PreSTIGE was employed to address the association of lncRNA with enhancer function, in which CTCF domains are excluded with the domain size expanded to 200 kb across the transcription start sites of the protein-coding genes. In addition, all enhancers that were overlapping with their predicted targets were screened from the enhancer datasets used in the subsequent analysis [19]. A previous study showed that AP001056.1 is a key immune-related eRNA in squamous cell carcinoma of the head and neck, which was associated with clinical outcomes through the above-mentioned method [20].
Mapping between gene symbol and Ensembl transcript ID was obtained by Ensembl BioMart. UCSC Xena browser [21] was also3 used to collect clinical data of cancer and putative eRNAs' levels. Candidate key eRNAs in OV were set as eRNAs with significant correlations with target genes levels (P<0.001, r > 0.3) and OS (P<0.05).

Functional enrichment analysis
Gene Ontology (GO) functional analysis was performed using cluster Profiler package of R software, so as for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of eRNA-related coding genes based on co-expression screening.

Survival and hazard analyses
In the present study, FOXP4-AS1 expressions were divided into FOXP4-AS1 high and FOXP4-AS1 low groups according to median expression value. Kaplan-Meier with log-rank test, univariate and multivariate Cox regression analyses were used to compare prognostic differences between different FOXP4-AS1 expression groups. Independent factors for survival were determined by Cox regression model. FOXP4-AS1 levels have also been investigated on its association with clinical characteristics, during which comparison between clinical variables in two groups was performed using Wilcoxon rank-sum test.

Sampling of tumor specimens
Tissue samples were collected from 42 patients with histologically confirmed OV between May 2017 and May 2018. All patients received cytoreductive surgery debulking followed platinum-based chemotherapy and received no treatment prior to surgery. Tumor tissue specimens were stored in RNA later solution (Ambion, Carlsbad, CA, U.S.A.) immediately after surgical removal, incubated at 4 • C for at least 24 h and subsequently stored at −20 • C until needed.
The study was approved by the Ethics Committee of Affiliated Xing Tai People Hospital of Hebei Medial University Hospital (2020 [015]). All patients provided written informed consent.

Quantification of FOXP4-AS1 and FOXP4 transcript levels in tissue samples
Total RNA was extracted using TRIzol reagent (Generay Biotech, Shanghai, Co., Ltd., China), according to manufacturer's protocol. Revert-Aid First Strand cDNA Synthesis Kit (Thermo Scientific, U.S.A.) was used to synthesize cDNA

Statistical analysis
All statistical analyses were performed using R software (Version 3.6.3). Comparisons between two groups were performed using Wilcoxon rank-sum test and comparisons between three or more groups were done by using Kruskal-Wallis test. Survival analysis was done by Kaplan-Meier analysis. The independent factors were determined by univariate and multivariate Cox regression analyses. P-values <0.05 were set as statistically significant.

Prognostic-related eRNAs screening in OV
Taking advantage of the PreSTIGE algorithm, 2695 lncRNA transcripts labeled by active tissue-specific enhancers derived ENCODE (Encyclopedia of DNA Elements database) and 2303 predicted target genes were identified. The clinical data of TCGA OV cohort came from UCSC Xena browser, as shown in Table 1. The conversion of transcript ID into gene ID were carried out by Ensembl BioMart. After that, 2695 transcripts were mapped to their corresponding 1288 genes. And eventually using Kaplan-Meier log-rank test, 71 gene were identified to be associated with OS significantly from these 1288 genes (Supplementary Table S1; Kaplan-Meier log-rank test, P<0.05). Spearman's rank correlation results showed that significant correlations were observed in just 23 of these 71 eRNAs when comparing them with predicted target gene mRNA levels, as shown in Table 2 (coefficient r > 0.3, P<0.001).

eRNA FOXP4-AS1 showed the highest positive correlation with its target gene FOXP4 in OV
Among the 23 prognostic-related eRNAs, it was found that the expression of FOXP4-AS1 has the highest positive correlation with its predicted target gene FOXP4 ( Figure 1A; Spearman's r = 0.61, P<0.001). Therefore, FOXP4-AS1 was confirmed as the research object from the 23 eRNAs in the present study. Another 185 co-expression genes also     Table S3; r > 0.3, P<0.001). To analyze pathways and biological functions, KEGG pathway and GO term analyses were performed on FOXP4-AS1 co-expression genes ( Figure 1B,C). The significantly enriched GO terms included transmembrane receptor protein serine/threonine kinase signaling pathway, RNA polymerase activity and ribosomal subunit. KEGG pathway analysis presented spliceosome, breast cancer, transcriptional misregulation in cancer and Wnt signaling pathway were also enriched according to KEGG pathway analysis.

High expression of FOXP4-AS1 correlates with the favorable prognosis of OV patients
The patients were divided into low and high groups based on the median value of FOXP4-AS1 expression in TCGA cohort. The FOXP4-AS1 high group showed longer OS when compared with FOXP4-AS1 low group (P=0.005; Figure  2A). Multivariable analysis of survival showed the FOXP4-AS1 expression levels were significantly associated with OS when tumor residual size, grade, stage and age were included (P=0.023, HR = 0.741, 95% CI = 0.572-0.959, Table 3). Moreover, lower FOXP4-AS1 was more associated with higher grade and lymphatic invasion compared with higher FOXP4-AS1 ( Figure 2B,C). These results showed that for OV patients FOXP4-AS1 may be a favorable independent prognostic biomarker candidate.

Validation of the role of FOXP4-AS1 in an independent OV cohort
Archived information of 42 OV patients was obtained from the Affiliated Xing Tai People Hospital of Hebei Medial University. All patients received platinum-based chemotherapy following primary debulking surgery and were followed up for 3 years. The median age was 55.2 years (ages ranged from 36 to 71). Among them, 42 patients were  Figure 3A). The expressions of FOXP4-AS1 and FOXP4 in stage III-IV were both lower than those in stages I-II OV (P=0.023, P=0.009, Figure 3B). The 42 patients were divided into low and high groups based on the median value of FOXP4-AS1 expression. Kaplan-Meier analysis indicated that the progression-free survival (PFS) and OS of patients with FOXP4-AS1 low group were shorter compared with FOXP4-AS1 high group ( Figure 3C,D).

Pan-cancer analysis of FOXP4-AS1
To evaluate tissue-specificity of FOXP4-AS1, we further investigated its role in pan-cancer using UCSC Xena browser (http://xenabrowser.net). As summarized in Table 4, the significant correlations between FOXP4-AS1 and FOXP4 existed in the 27 tumor types. Moreover, the prognostic effect of FOXP4-AS1 in pan-cancer was also explored. Besides OV, the significant impact of FOXP4-AS1 on survival was also observed in Adrenocortical carcinoma (ACC), Esophageal carcinoma (ESCA), Kidney renal clear cell carcinoma (KIRC), Lower Grade Glioma (LGG), Liver hepatocellular carcinoma (LIHC) and Mesothelioma (MESO) (P<0.05, Figure 4).  Discussion eRNAs have attracted particular interest due to their potential roles in mediating enhancer functions and gene transcription associated with cancers [22]. The aim of the present study is to explore the eRNAs as prognostic related biomarker in OV. The results revealed that: (1) for the first time, FOXP4-AS1 was determined to be an eRNA with top positive correlation with its target gene FOXP4 among the prognosis-related eRNAs in OV.
(2) FOXP4-AS1 may be an independent biomarker in OV patients' prognosis. (3) The prognostic significance of FOXP4-AS1 as well as the positive correlation with FOXP4 were further validated in an independent OV cohort. (4) The impact of FOXP4-AS1 on survival was presented in pan-cancer analysis. To date, several models have been developed to show how eRNAs regulate the transcription of targeted genes. However, the mechanisms are quite different [23]. In addition to promoter-enhancer looping and histone modifications, eRNAs promote targets transcription through enhancing the binding of RNA Pol II or leasing negative regulator [22]. As shown in Figure 5A, eRNAs can independently recruit RNA Pol II to both enhancer and promoter loci. Once recruited, eRNAs interact with the NELF and P-TEFb complexes to regulate pause-release of RNA Pol II and promote enhancer and gene transcription [24]. However, compared with sense eRNAs, the mechanism of antisense eRNAs in regulation of target genes was kind of different. Pan et al. [25] found that antisense eRNA mediated mRNA through binding to DNMT1. Promoter-gene-ending interaction increases the specific antisense eRNA targeting gene-ending region. In the present study, FOXP4-AS1 was found to be top correlated with its predicted target FOXP4 in OV. Furthermore, the correlation was validated in an independent 42 OV patients. Similarly, Li et al. [26] reported that FOXP4-AS1 and FOXP4 were up-regulated in esophageal carcinoma samples and were positively correlated with each other based on bioinformatics prediction. Therefore, we speculated that the regulation of antisense eRNA in the predicted targets may be a possible mechanism of how FOXP4-AS1 contributes to FOXP4 transcription ( Figure 5B). Our study presented the putative regulatory roles of eRNA FOXP4-AS1 on target gene FOXP4 in OV, nevertheless, further experiments are still necessary to confirm the underlying mechanism.
FOXP4-AS1 is an emerging cancer-related biomarker, whose roles in different types of cancers were inconsistent. As previous research reported, FOXP4-AS1 was confirmed to be an oncogene in colorectal cancer associated with a dismal prognosis [27]. The consistent results were also shown in prostate cancer [28], pancreatic ductal adenocarcinoma [29], and osteosarcoma [30]. More recently, Binang et al. pointed out that increased FOXP4-AS1 expression was associated with beneficial outcomes considering gastric cancer disease-free survival based on public databases [31]. In this study, high FOXP4-AS1 expression was related to better OS based on TCGA OV patients. Multivariable analysis of survival showed the FOXP4-AS1 expression levels were significantly associated with OS when tumor residual size, grade, stage, and age were included. Additionally, FOXP4-AS1 expression showed a correlation with some cancer-related clinical features including clinical stage, grades and lymphatic invasion. In the validation stage, the FOXP4-AS1 high group is closely associated with better survival than the FOXP4-AS1 low group in a single OV cohort. In summary, FOXP4-AS1 may play a role of tumor suppressor in OV. However, strong correlations between the high expression of FOXP4-AS1 and poor prognosis were found in pan-cancer survival analysis including ACC, ESCA, LIHC, brain LGG, KIRC, MESO. Compared with these types of cancer, the different prognosis values of FOXP4-AS1 in OV may be attributed to the cancer type-specific eRNAs. Approximately 59% of all identified eRNAs (5332 out of 9108) were confirmed to be cancer type-specific from a previous report [17], from which we know cancer type and unique pattern were closely associated for eRNA expression. Thus, FOXP4-AS1 expression pattern may be obviously diverse in the different cancer cells. Our results supplied more information on the potential prognosis-related role of FOXP4-AS1 in different cancer types. The converse impact of FOXP4-AS1 on survival in OV made it a noteworthy eRNA for this tumor type.
FOXP4, one member of the forkhead box P family [32], was shown to be most positively correlated with FOXP4-AS1 expression among the 186 potential target genes depending our criteria in this study. Moreover, positive correlations of the two genes were also observed in the pan-cancer analysis. FOXP4 was observed in different cancer types as a tumor suppressor or oncogene previously [33]. FOXP4 positively regulated by eRNA FOXP4-AS1 may be one of the mechanisms that involved in the biological processes of OV. In addition, the other potential target genes and transcripts form the functional clusters may be majorly involved in cancer progression in OV. It cannot be ignored the influence through sponging miRNAs, lncRNAs regulate mRNA expression of these targets as ceRNAs on post-transcriptional level. The functional role of eRNA FOXP4-AS1 in OV remains need to determine in vivo and in vitro experiments.

Conclusion
In this current study, we explored the prognosis-related eRNAs in OV patients as well as related target genes. Study results suggested FOXP4-AS1 may be a favorable independent prognostic biomarker candidate in OV with top positive correlation with its target gene FOXP4. However, the role of FOXP4-AS1 in prognosis was inconsistent in pan-cancer analysis. Moreover, the target genes of FOXP4-AS1 are mainly involved in cancers, implying the importance of eRNA FOXP4-AS1 in the biological progresses. These results provided new ideas in exploring the regulatory network of eRNAs in the process of OV evolution and prognosis, and thus to reveal the panorama of the disease. It should be a valuable field for further investigation of cancer therapy.

Data Availability
All datasets generated for the present study are included in the Supplementary material.