Downregulation of enhancer RNA EMX2OS is associated with poor prognosis in kidney renal clear cell carcinoma

Enhancer RNAs are a subclass of long non-coding RNAs transcribed from enhancer regions that play an important role in the transcriptional regulation of genes. However, their role in kidney renal clear cell carcinoma (KIRC) is largely unknown. Herein, we identified the key enhancer RNAs in KIRC via an integrated data analysis method. Gene expression profiles and clinical data of KIRC and 32 other cancer types were acquired using the University of California Santa Cruz Xena platform. Reported enhancer RNAs and genes regulated by them were selected as putative enhancer RNA-target pairs. Kaplan-Meier survival and correlation analyses were performed to identify the key enhancer RNAs. Finally, EMX2OS was identified as the enhancer RNA most associated with survival, with EMX2 as its target. EMX2OS downregulation was significantly associated with higher histological grade, advanced stage, and poorer prognosis. The results were validated in pan-cancer data from The Cancer Genome Atlas and RT-qPCR analysis of 12 pairs of KIRC and normal real-world samples. Functional enrichment analysis indicated that several metabolism-associated signaling pathways were enriched. This study demonstrated that EMX2OS is a key metabolism-associated enhancer RNA in KIRC with a favorable impact on survival and may be a novel therapeutic target in KIRC.


INTRODUCTION
Kidney renal clear cell carcinoma (KIRC) is the most common histological type of renal cell carcinoma and is associated with the highest metastasis and mortality rates [1]. The number of new cases of KIRC was 73,820 in the United States, and it had caused approximately 14,770 deaths based on cancer statistics data in 2019 [2]. Contrast-enhanced, triplephase helical computed tomography (CT) and surgery are the most important diagnostic and treatment methods, respectively [3]. Although the treatment for KIRC has improved significantly, its mortality rate is still increasing [2]. Up to 50% of patients who undergo nephrectomy will progress to distant metastasis, and the 5-year survival rate of patients with metastasis is below 12% [4]. The TNM stage remains the most important prognostic predictor for KIRC, but the survival rate can differ greatly for patients even with the same stage because of tumor heterogeneity [5]. Therefore, it is essential to identify effective prognostic biomarkers to identify high-risk patients with poor survival.

AGING
Enhancer RNAs (eRNAs) are a subclass of long noncoding RNAs transcribed from enhancer regions. eRNAs are known to promote the formation of enhancer-promoter loops and play diverse and emerging roles in the regulation of gene expression and cell fate [6][7][8]. In addition, increasing evidence suggests that eRNAs are closely associated with carcinogenesis [9][10][11][12][13]. However, few studies have investigated the role of eRNAs in KIRC. In this study, we explored and identified the eRNAs associated with survival in patients with KIRC. We found that eRNA EMX2OS, located in the tissue-specific enhancer of a tumor suppressor gene EMX2 [14][15][16], was positively correlated with favorable clinicopathological characteristics and was significantly related to overall survival (OS) in patients with KIRC.

Putative prognostic eRNAs in KIRC
After matching the patients' clinical information and gene expression data, 535 patients with KIRC were enrolled in the study. Among them, 531 patients had information regarding survival. The patients' clinicopathological features are summarized in Table 1. A total of 2695 eRNAs, which were annotated using the Encyclopedia of DNA Elements database and derived from tissue-specific enhancers, as well as 2303 predicted target genes, have been identified previously using the PreSTIGE algorithm [17]. We used these eRNA-target pairs to identify potential key eRNAs in KIRC. Finally, 16 eRNA-target pairs were identified with certain conditions (Kaplan-Meier log rank of p < 0.001 and correlation coefficient of r > 0.6 and p < 0.001) ( Table 2).

EMX2OS is a key eRNA in KIRC
As shown in Table 2, EMX2OS was the putative eRNA most associated with survival, and it showed an obviously positive correlation with its predicted target EMX2. Thus, we performed further analysis for EMX2OS. Based on the expression data of 535 KIRC samples and 72 normal kidney samples, we found that EMX2OS expression was significantly lower in tumor tissues than in normal tissues (p < 0.001) ( Figure 1A). According to the median expression level of EMX2OS, 531 patients with KIRC with survival information were divided into the high-EMX2OS and low-EMX2OS expression groups. Kaplan-Meier survival analysis demonstrated that EMX2OS downregulation was significantly associated with poorer OS (hazard ratio [HR] = 0.64, p < 0.001, Figure 1B). Lower EMX2 expression was also related to poorer OS (HR = 0.68, p < 0.001, Figure 1C). Furthermore, there was a strong co-expression relationship between EMX2OS and its predicted target EMX2 (correlation coefficient r = 0.81; p < 0.001) ( Figure 1D). We next investigated the relationship between EMX2OS expression and the clinicopathological characteristics of KIRC. We found that the EMX2OS expression level was significantly related to several clinicopathological features of KIRC, including gender (p = 0.001), cancer status (p < 0.001), race (p = 0.021), histological grade (p < 0.001), tumor size (T stage, p < 0.001), lymph node stage (N stage, p = 0.008), distant metastasis (M stage, p = 0.001), and American Joint Committee on Cancer (AJCC) stage (p < 0.001) (Table 1). Moreover, EMX2 expression showed a similar relationship with the clinicopathological characteristics of KIRC (Table 1). Further analysis showed that EMX2OS expression was negatively correlated with histological grade (p < 0.001), T stage (p < 0.001), and AJCC stage (p < 0.001) ( Figure 1E-1G). To validate the results, we explored the differential expression and prognostic value of EMX2OS and its correlation with EMX2 in pan-cancer data (32 other cancer types) from The Cancer Genome Atlas (TCGA) as an internal validation. Interestingly, EMX2OS expression was significantly lower in many other types of cancers, including breast invasive carcinoma, colon adenocarcinoma, kidney chromophobe, and liver hepatocellular carcinoma, than in normal tissues ( Figure 1H). EMX2OS also played a prognostic role in adrenocortical carcinoma, cervical squamous cell carcinoma and endocervical adenocarcinoma, stomach adenocarcinoma, and uveal melanoma (Table 3). Of interest, EMX2OS and EMX2 showed significant correlations in all cancers, except cholangiocarcinoma (Table 3). These results suggest that EMX2OS functions as a tumor suppressor in KIRC.

Independent prognostic value of EMX2OS in KIRC
As Kaplan-Meier survival analysis showed that patients with different EMX2OS expressions had significant different OS, we further explored whether EMX2OS had an independent prognostic value in KIRC using univariate and multivariate Cox regression analyses. Because information on the M stage, N stage, T stage, cancer status, and race were missing or were unevenly distributed in most cases, these data were excluded. Finally, 520 patients with information on age, gender, histological grade, AJCC stage, and EMX2OS expression were included in the Cox regression analysis. We found that there was a significantly prognostic difference between the high-EMX2OS and low-EMX2OS expression groups in both univariate (HR, 0.654; 95%CI, 0.579 -0.739; p < 0.001) and multivariate (HR, 0.783; 95%CI, 0.684 -0.896; p < 0.001) Cox regression analysis (Figure 2A, 2B). EMX2OS expression, age, histological grade, and AGING  AJCC stage had independent prognostic values in KIRC. Additionally, the stratified analysis revealed that the patients in the low-EMX2OS expression group had significantly poorer OS than those in the high-EMX2OS expression group regarding age ≤ 60 years (p = 0.001), age > 60 years (p < 0.001), male sex (p = 0.002), female sex (p < 0.001), AJCC stage I/II (p = 0.044), AJCC stage III/IV (p = 0.007), and histological grade 3/4 (p < 0.001) but not histological grade 1/2 (p = 0.235) ( Figure  3A-3D).

Determination of EMX2OS and EMX2 levels using RT-qPCR
Quantitative reverse transcription polymerase chain reaction (RT-qPCR) was performed to measure the relative expression of EMX2OS and EMX2 in 12 patients with KIRC. Compared with the matched tumorfree samples, EMX2OS and EMX2 expression levels were downregulated in 91.7% (11 of 12) and 83.3% (10 of 12) of KIRC samples, respectively ( Figure 4A).
However, compared with normal samples, only EMX2OS was significantly downregulated in KIRC samples (p = 0.013), although it was nearly significantly different for EMX2 (p = 0.054), possibly because of the small sample size ( Figure 4B). Additionally, there was an obviously significant positive correlation between EMX2OS and EMX2 in both normal (r = 0.85; p < 0.001) and tumor samples (r = 0.89; p < 0.001) ( Figure  4C).

Functional enrichment analysis
To further investigate the function of EMX2OS, we identified 2124 significantly co-expressed genes with EMX2OS in KIRC (r > 0.4; p < 0.001) (Supplementary Table 1). Gene Ontology (GO) functional enrichment analysis showed that several catabolic processes, including small molecule, carboxylic acid, and fatty acid catabolic processes, were significantly enriched ( Figure 5A). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that the AGING metabolism of several important energy substances, including fatty acid, pyruvate, and tryptophan, was enriched. Moreover, FOXO and PPAR signaling pathways were also enriched ( Figure 5B).

DISCUSSION
The treatment for KIRC remains a challenge across the world because of its substantial morbidity and mortality   [18]. Thus, identifying novel markers for prognosis and therapy is needed to improve the survival of patients with KIRC. Recently, with rapid advanced in highthroughput sequencing technology and bioinformatics, more novel biomarkers have been identified. One of the most surprising discoveries is eRNAs. Accumulation evidence has shown that the dysregulation of eRNAs is closely associated with various human diseases and that eRNAs are a new therapeutic target [19][20][21][22].
eRNAs are transcribed from putative enhancer regions distinguished by high levels of histone H3 lysine 27 acetylation and histone H3 lysine 4 monomethylation and low levels of histone H3 lysine 4 trimethylation. PreSTIGE is a recently developed algorithm to predict tissue-specific enhancers and their targets based on the H3K4me1 marker and tissue-specific expression of mRNAs [23]. Using this enhancer prediction approach, 2695 eRNAs and their predicted targets were identified in a previous study [17]. As the role of eRNAs had been rarely studied in KIRC, we used these 2695 eRNAtarget pairs as candidates to identify key eRNAs in KIRC. Using Kaplan-Meier survival and correlation analyses, EMX2OS was found to be the most survival-  AGING associated eRNA with high correlation with its predicted target EMX2.
EMX2OS is a long non-coding RNA, an antisense transcript from EMX2 [24]. EMX2 encodes a homeobox protein and is a vital gene that promotes the formation of the urogenital and central nervous systems during embryonic development [25,26]. In cancer, EMX2 usually acts as a tumor suppressor [15,16].
Although rarely reported, we speculate that EMX2OS also serves as a tumor suppressor because of its strong co-expression relationship with EMX2. In this study, we found that EMX2OS was significantly downregulated in KIRC tissues compared with that in normal kidney tissues. EMX2OS downregulation was significantly associated with unfavorable clinicopathological features such as higher histological grade, larger tumor size, lymph node metastasis, distant  AGING metastasis, and advanced AJCC stage. Kaplan-Meier survival analysis showed that downregulated EMX2OS expression was significantly associated with poorer OS. Furthermore, multivariate Cox regression analysis confirmed that EMX2OS played an independent prognostic role in KIRC. Additionally, patients with lower EMX2OS expression had significantly poorer OS than patients with higher EMX2OS expression in stratified analysis according to age, gender, and AJCC stage. However, a significant difference in survival was found in patients with histological grade 3/4 but not in those with histological grade 1/2. A possible reason for this difference is that the tumors with low histological grades were relatively less malignant.
To validate the aforementioned results, we used data of 32 other types of cancer from TCGA as internal validation and RT-qPCR analysis of 12 pairs of KIRC and normal real-world samples as external validation. Pan-cancer analysis showed that significant expression and survival differences of EMX2OS existed in several types of cancers. The RT-qPCR results further confirmed the downregulation of EMX2OS in KIRC. All results were consistent and suggested that EMX2OS serves as a tumor suppressor in KIRC.
GO and KEGG functional enrichment analyses provided some clues on how EMX2OS influences patients' survival. These analyses showed that the metabolism of several energy-involved substances, such as fatty acid and pyruvate, were enriched, which suggests that EMX2OS influences the energy metabolism of the tumor. In the KEGG pathway analysis, FOXO and PPAR signaling pathways were also enriched. The FOXO signaling pathway plays an important role in the regulation of metabolic homeostasis and suppression of tumor growth [27][28][29]. PPARs are nuclear receptors that regulate cellular and whole-body energy homeostasis during carbohydrate and lipid metabolism, cell growth, and cancer development [30]. Alternatively, energy metabolism reprogramming is an emerging hallmark of cancers [31]. Thus, we speculated that EMX2OS may regulate energy metabolism by enriching the FOXO and PPAR signaling pathways to influence the prognosis of KIRC.
Although this study discovered the potential prognostic value of EMX2OS in KIRC for the first time, several limitations need to be considered. First, our data were mainly sourced from TCGA, and most patients were white. Second, the tissue specimens used for RT-qPCR analyses were limited. Therefore, additional data and samples from different races are necessary to confirm the results of this study. Finally, basic experiments are essential to detect the molecular mechanism of EMX2OS in KIRC.
In conclusion, this study demonstrated that EMX2OS is a key survival-associated eRNA in KIRC. With a potential role in energy metabolism, EMX2OS may be a novel therapeutic target in patients with KIRC.

Data extraction and identification of prognostic eRNAs in KIRC
The gene expression profiles and clinical data of KIRC and 32 other types of cancers from TCGA were acquired using the University of California Santa Cruz Xena database (https://xena.ucsc.edu/). We matched the patients' clinical information and gene expression data. Thus, patients with clinical information and gene expression data were enrolled for further analysis. Next, we obtained a list of eRNAs transcribed from active tissue-specific enhancers and their target predicted using the Predicting Specific Tissue Interactions of Genes and Enhancers (PresSTIGE) method [17,23]. The association between the level of putative eRNAs and OS of patients with KIRC was investigated using the R packages "survival" and "survminer" (The R Foundation for Statistical Computing, Vienna, Austria). Co-expression analysis was also performed to evaluate the correlation between the level of eRNAs and their predicted targets. Putative eRNAs were considered candidates if they satisfied the following two criteria: significant association with OS (Kaplan-Meier log rank of p < 0.001) and co-expression with the predicted target (r > 0.6 and p < 0.001). Then, we selected the most significant survival-associated eRNA for further analysis. The differential expression between tissues with KIRC and normal tissues was explored. The independent prognostic role of eRNAs in KIRC was investigated. The correlation between the expression level of the selected eRNA and clinicopathological characteristics in KIRC was also assessed. Finally, we validated the results using the pan-cancer data from TCGA (32 other types of cancers) and RT-qPCR results of 12 pairs of KIRC and normal real-world samples as internal and external validations, respectively.

RT-qPCR
RT-qPCR was used to further validate the RNAsequencing data obtained from TCGA. We collected KIRC tissues and paired adjacent normal tissues from 12 patients who underwent nephrectomies or partial nephrectomies at Meizhou People's Hospital between 2019 and 2020. Informed consent was obtained from all patients. We extracted total RNA using the TRIzol™ reagent (Waltham, Massachusetts, USA). The PrimeScript RT reagent kit (Takara Bio, Inc., Dalian, China) was used to synthesize complementary DNAs following the manufacturer's protocols. The SYBR Green PCR kit (Takara Bio, Inc., Dalian, China) was used to conduct quantitative real-time PCR using the ABI 7500 fluorescent quantitative PCR system (Applied Biosystems Inc., Foster City, CA, USA). We used glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as the internal control. The primer sequences were as follows: EMX2OS (forward: GTGACTTG CACAAGGACACAA; reverse: CCTGTCTGGCCAT TCCTCT), EMX2 (forward: CGGCACTCAGCT ACGCTAAC; reverse: CAAGTCCGGGTTGGAGT AGAC), and GAPDH (forward: ATGACATCAA GAAGGTGGTG; reverse: CATACCAGGAAATGA GCTTG). The expressions of EMX2OS and EMX2 were measured using 2−ΔΔCt method.

Functional enrichment analysis
To determine the underlying molecular mechanisms, we explored the co-expressed genes with the selected eRNA (correlation coefficient r > 0.4, p < 0.001). Then, GO and KEGG analyses were performed.

Statistics
All data processing and statistical analysis were performed using R (version 3.6.1; The R Foundation for Statistical Computing, Vienna, Austria), Strawberry Perl (version 5.30.1.1; http://strawberryperl.com/), and Statistical Package for Social Sciences (version 25.0; IBM, Armonk, New York, USA). Analysis of variance or t-test was used to compare the gene expression level among different subgroups. Spearman's rank correlation coefficient was used to evaluate the correlation strength. A result was considered statistically significant when the p value was < 0.05, except when stated otherwise.

AUTHOR CONTRIBUTIONS
HMJ conceived and designed the study, performed the RT-qPCR, and wrote the manuscript. NHC designed the study and revised the manuscript. HBC analyzed the data and revised the manuscript. PW analyzed the data. SDS generated the tables and figures. The authors read and approved the final manuscript.