Discovering the interactome, functions, and clinical relevance of enhancer RNAs in kidney renal clear cell carcinoma

Abstract

Enhancer RNA (eRNA) has emerged as a key player in cancer biology, influencing various aspects of tumor development and progression. In this study, we investigated the role of eRNAs in kidney renal clear cell carcinoma (KIRC), the most common subtype of renal cell carcinoma. Leveraging high-throughput sequencing data and bioinformatics analysis, we identified differentially expressed eRNAs in KIRC and constructed eRNA-centric regulatory networks. Our findings revealed that up-regulated eRNAs in KIRC potentially regulate immune response and hypoxia pathways, while down-regulated eRNAs may impact ion transport, cell cycle, and metabolism. Furthermore, we developed a diagnostic prediction model based on eRNA expression profiles, demonstrating its effectiveness in KIRC diagnosis. Finally, we elucidated the regulatory mechanism of an eRNA (ENSR00000305834) on the expression of SLC15A2, a potential prognostic biomarker in KIRC, through bioinformatics analysis and in vitro validation experiments. In summary, Our study highlights the clinical significance of eRNAs in KIRC and underscores their potential as therapeutic targets.

Introduction

The progress made in high-throughput sequencing technology in the last decade has substantially enhanced our comprehension of the genome and transcriptome. Enhancer RNA (eRNA), once dismissed as "transcriptional noise" or a "by-product" [1, 2], has garnered increased attention, with a growing number of studies unveiling its significant role. Recent systematic investigations into eRNA expression in diverse human cancer samples have illuminated its potential expansive function in tumorigenesis, providing novel insights into the mechanisms underlying enhancer actions [3,4,5].

Renal cell carcinoma (RCC) stands out as a prevalent tumor in the urinary tract, with kidney renal clear cell carcinoma (KIRC) representing the most common subtype, encompassing approximately 80% to 90% of all cases [6,7,8]. Regrettably, KIRC also holds the unfortunate distinction of being the most lethal among kidney cancer subtypes [9]. Previous studies usually identified some KIRC biomarkers through differential expression analysis and combined them with prognostic analysis to demonstrate their clinical significance. As a pivotal class of DNA regulatory elements, enhancers conventionally exert control over the expression of target genes by establishing spatial chromatin loops with target promoters, resulting in either up- or down-regulation [10]. Barata and Rini's prior transcriptomic analyses demonstrated that noncoding RNAs, such as microRNAs and long noncoding RNAs, regulate oncogenic pathways, including aberrant hypoxia-inducible factor (HIF) signaling and metabolic reprogramming [11]. However, the role of enhancer RNAs (eRNAs), a subclass of noncoding RNAs critical for gene transcription regulation, remains largely unexplored in KIRC. While individual eRNA functions in KIRC have been studied, their broader role is unexamined [12, 13]. This study fills that gap by mapping the eRNA interactome and functional pathways, highlighting their clinical relevance in tumor progression and therapeutic response.

Solute carrier family 15 (H + /peptide transporter) member 2 (SLC15A2) functions as a proton-coupled peptide transporter expressed in the mammalian kidney, facilitating the uptake of small peptides, β-lactam antibiotics, and various peptide drugs from tubular filtrate [14]. SLC15A2 encodes the H + /peptide transporter 2 (PEPT2), a constituent of the solute carrier 15 family (SLC15) [15]. Research conducted by Huang, H., et al., demonstrated the differential expression of SLC15A2 in a variety of cancer samples [16]. Consequently, various types of cancer research have honed in on investigating the role of PEPT2 (SLC15A2) in disease progression. Both PEPT1 and PEPT2 (SLC15A2) exert an impact on tumor metabolism by regulating the cellular uptake of Gly-Sar and L-histidine in prostate cancer cell lines PC-3 and LNCaP [17]. The ability of PEPT2 (SLC15A2) to regulate transmembrane transport is underscored by Zimmermann et al., whose study on the transport of dipeptides through PEPT2 in the glioma cell line U373-MG highlights its potential as a therapeutic target for glioma [18]. In hepatocellular carcinoma (HCC), whole-genome sequencing analysis revealed that a single nucleotide variant (rs2257212) of SLC15A2 correlates with better progression-free survival and modulates the response to sorafenib treatment [19]. Within the kidney, PEPT2 (SLC15A2) assumes a primary role in the tubular reabsorption of peptide-bound amino acids [20]. In addition, PEPT2 (SLC15A2) has also been found to be a potential prognostic biomarker in kidney cancer [21, 22]. In the course of our study, we discovered that the expression of SLC15A2 in KIRC samples may be regulated by enhancer RNA (eRNA) ENSR00000305834. Furthermore, other eRNAs may also play a role in regulating the expression of target genes, offering a valuable reference for eRNA biomarkers applicable to KIRC prognostic research.

Methods

eRNA data source

The chromatin coordinates, survival relevance, and expression levels of eRNAs in KIRC were sourced from the eRic database [4]. Additionally, the expression levels of eRNAs in the normal kidney cortex were extracted from the HeRA database [23].

H3K27ac ChIP-seq data analysis

H3K27ac ChIP-seq data were obtained from the publicly available study (GSE75597) [24]. The data preprocessing steps, including alignment, reads filtration, duplicate removal, and peak calling, closely adhered to the procedures outlined in the original study [24]. Briefly, the raw data was first trimmed and low quality reads were filtered out using cutadapt [25] (-q 15). The remaining reads were mapped to human hg38 reference using BWA-0.6.2 [26] with default parameters. Next, the PCR duplicates were removed by picard MarkDuplicates (http://broadinstitute.github.io/picard/). Finally, the peaks were called by Homer [27] v4.3 makeTaqDirectory. The eRNAs with transcription start site located within the H3K27ac peaks were detected by Bedtools [28] and defined as the eRNAs derived from putative enhancers.

Differentially expression analysis on eRNAs

The normal and KIRC expression profiles of eRNAs were obtained from HeRA (n = 38) and eRic database (n = 530), respectively. To compare the expression levels of eRNAs in KIRC and the normal kidney cortex, normalization to Fragments Per Kilobase Million (FPKM) was initially applied. Subsequently, Student's t-test was employed to identify significant differences in eRNA expression between tumor and normal samples. The eRNAs with fold changes greater than 1.5 and p-values lower than 0.05 were defined as the up-regulated eRNAs in KIRC, and those with fold changes lower than 1/1.5 and p-values lower than 0.05 were defined as down-regulated genes. This rigorous criterion aimed to pinpoint eRNAs with substantial alterations in expression levels, contributing to a more robust analysis.

Co-expression network construction

The co-expression network based on 526 TCGA-KIRC samples was constructed using eRNA-IDO. Briefly, the correlations of gene expression levels were calculated using GCEN [29]. Significantly correlated pairs ofeRNAs and annotated genes, determined by absolute Pearson correlation coefficients exceeding 0.5 and false discovery rates (FDRs) below 0.05, were identified. Since eRNAs typically facilitate enhancer-mediated transcriptional activation, only positively correlated pairs were retained. For correlations between two annotated genes, an absolute Pearson correlation coefficient threshold of 0.8 was applied. Subsequently, all correlated pairs were amalgamated to construct the gene co-expression network.

eRNA-centric regulatory network construction

The eRNA-centric regulatory network was contructed using eRNA-IDO [30], in which a comprehensive collection of 11,356 TF ChIP-seq datasets were utilized to pinpoint TF binding sites on eRNA regions, 200 HiChIP datasets involving 108 biosamples were obtained from HiChIPdb [31] to unveil eRNA-associated Enhancer-Promoter (E-P) loops. TF binding sites overlapping the eRNA regions were defined as eRNA-associated TFs. E-P loops with anchors overlapping eRNA-expressing regions were deemed eRNA-associated loops. TF-eRNA and eRNA-loop relationships were then consolidated to form TF-eRNA-loop regulatory axes. This regulatory information was further integrated with protein–protein interactions from the STRING database [32] to construct the eRNA-centric regulatory network.

Module extraction

SPICi [33] in the unweighted mode with default parameters was used to extract the modules from the networks. The networks with at least one eRNA were reserved for further analyses.

Network visualization

Networks were visualized using Cytoscape v3.9.1 [34].

Functional enrichment analyses

The functional enrichment analyses on eRNA-associated protein-coding genes were conducted using the script from ncFANs v2.0 (https://github.com/zhangyw0713/FunctionEnrichment) [35].

Establishment and evaluation of the diagnosis prediction model

The expression profiles and clinial information of GTEx normal kidney cortex samples (n = 38) and TCGA KIRC samples (n = 530) obtained from HeRA and eRic database were used to build the diagnosis prediction model. To get a balanced dataset (positive:negative = 1:1), KIRC samples were downsampled. Next, half of the samples were extracted to create the training set, and the remaining samples served as the test set. The univariate logistic regression model with eRNA expression levels as the input variable was used to obtain the DE eRNAs significantly related to the oncogenesis, which was further used to construct the diagnosis prediction model using the multivariate logistic regression model. Finally, receiver operating characteristic (ROC) curves were drawn using the R package ‘pROC’ to estimate the performance of the diagnosis prediction model.

Cell lines and reagent

The 786-O human KIRC cell line was obtained from the Shanghai Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences, situated in Shanghai, China. These cells were cultured in Complete RPMI 1640 medium supplemented with 10% fetal bovine serum and 1% penicillin–streptomycin solution (obtained from Solarbio, China). Maintenance of all cell lines occurred in a controlled environment set at 37 °C with a 5% CO2 concentration. Routine assessments for mycoplasma contamination were performed using PCR techniques.

The overexpression of eRNA in the 786-O cell line

We employed an overlap analysis between the predicted genomic location of the enhancer RNA (eRNA) ENSR00000305834 obtained from the eRic database (chr3:122288200–122294200) and gro-seq data from human embryonic kidney cell lines sourced from the GEO database (accession ID: GSM2428721) to determine the sequence of the eRNA for subsequent validation in in vivo experiments.

To establish cells with stable expression of ENSR00000305834, we infected the 786-O cell line with lentiviral particles obtained from GenePharma (Shanghai, China). Following a 48-h incubation period, cells were subjected to selection for 7 days in culture medium containing 5 μg/ml puromycin to enhance expression levels. The specific primers utilized in this experiment are provided in the Supplementary Table 1.

Western blot

The cellular lysates were prepared by subjecting the cells to lysis in ice-cold RIPA lysis buffer supplemented with 1 mM phenylmethylsulfonyl fluoride (PMSF) to ensure protein integrity. Subsequently, the protein content in the lysates was quantified using a BCA protein assay kit to ensure equal loading. Equal amounts of total protein were then resolved by electrophoresis on a 10% SDS-PAGE gel and transferred onto PVDF membranes. Following transfer, the membranes were probed with specific antibodies against SLC15A2 and β-actin to enable target protein detection. Following thorough washing with TBST, the membranes were incubated with an HRP-conjugated secondary antibody, facilitating the visualization of protein bands using the ChemiScopeTouch imaging system. The intensity of the immunoblot bands was quantified using ImageJ 1.8.0 software, ensuring accurate measurement of protein expression levels.

Quantitative real-time PCR (qRT-PCR)

Total RNA was extracted using TRIzol reagent (Invitrogen, USA). cDNA was obtained by reverse transcription system. Subsequently, qPCR was performed using LightCycler 480 fluorescence quantitative system (Roche, Basel, Switzerland). Each experiment was repeated five times, and then the 2^−ΔΔCT method was used to calculate relative expression and display it using the “ggplot2” R package. Detailed information of primers can be found in Supplementary Table 1.

Colony formation assay

Cells from different experimental groups were seeded into 6-well plates at a density of 5000 cells per well. These cells were then allowed to grow undisturbed for a period of 10 to 14 days, with regular replenishment of the culture medium every two days to facilitate optimal colony formation. After the designated incubation period, the colonies were carefully washed with PBS to remove any non-adherent cells. Next, the colonies were fixed by incubating them with a 4% formaldehyde solution for 30 min. Following fixation, the colonies were stained with a 0.5% crystal violet solution prepared in 10% methanol for an additional 30 min to enable visualization. Subsequently, the plates were gently washed to remove excess stain, air-dried, and the colonies were then visualized.

Cell growth assay

After transfection with siRNA and negative control RNA, 786-O cells were plated into 6-well plates at a density of 10,000 cells per well. The cell number was then assessed on days 2, 4, and 6 post-seeding to evaluate the growth rate. This allowed for a comparative analysis of cell proliferation between the experimental groups over the specified time points.

Results

Identification of the actively transcribed eRNAs in human KIRC

We retrieved 1275 eRNAs associated with human KIRC from the eRic database. To bolster credibility, we utilized ChIP-seq data specific to H3K27ac modification (Fig. 1B) for filtering, resulting in the identification of 252 eRNAs transcribed from active enhancers (Fig. 1A). We also found that the overall expression level of these 252 active eRNAs was significantly higher than the remaining more than 1,000 inactive eRNAs (Fig. 1C).

Fig. 1

Differentially expressed eRNAs and annotated DEGs in KIRC. A Research technology roadmap; B Normalized signal of H3K27ac modification or not; C The expression level of active eRNAs and inactive eRNAs; D The pie plot of the prognostic relevance of these 252 active eRNAs; E The Kaplan–Meier survival curve of the ENSR00000022815 in KIRC

To delve into the clinical significance of these 252 active eRNAs in KIRC, the eRNA-associated survival data based on 530 TCGA-KIRC samples were obtained from eRic database. Strikingly, 45.85% of these active eRNAs exhibited significant associations with KIRC prognosis (Fig. 1D), highlighting a substantial role in the onset and progression of the disease. As an illustrative example, ENSR00000022815 emerged as a potential prognostic eRNA. The Kaplan–Meier survival curve depicted in Fig. 1E demonstrates that the Overall Survival (OS) of KIRC patients in the high-expression group significantly lagged behind those in the low-expression group (p = 0.0181). This compelling evidence underscores the clinical relevance and potential prognostic value of the identified active eRNAs in the context of KIRC.

Up-regulated eRNAs in KIRC potentially regulate immune response and hypoxia

The eRNA expression data from the GTEx and TCGA databases, comparing KIRC and normal kidney cortex, were utilized to identify differentially expressed eRNAs. Among these, 121 eRNAs were found to be up-regulated, while 55 were down-regulated in KIRC patients (Fig. 2A). To delve deeper into the regulatory roles of these eRNAs, co-expression analysis was conducted on the differentially expressed eRNAs and 3,573 annotated Differentially Expressed Genes (DEGs) between KIRC and normal samples, employing default standards from the oncoDB database. Notably, 75.46% of the differentially expressed eRNAs exhibited a positive correlation with the expression levels of their target genes (Fig. 2B), aligning with the findings of previous studies [36, 37]. Typically, eRNAs regulate Enhancer-Promoter (E-P) loops, thereby stimulating the expression of target genes. Subsequently, we only focus on the parts where eRNA and target genes are positively correlated.

Fig. 2

Regulatory network of the up-regulated eRNAs in KIRC. A The volcano plots of the differentially expressed eRNAs between KIRC samples and normal samples; B Correlation of the differentially expressed eRNA and the differentially expressed annotated target genes; C The co-expression network between up-regulated eRNAs and positively correlated targets; D The barplot of GO enrichment analysis results of the positively correlated targets; E, H The two most densely connected modules from the co-expression network; F, I The barplot of GO enrichment analysis results of these two modules; G, J The barplot of KEGG enrichment analysis results of these two modules

A co-expression network was constructed between up-regulated eRNAs and their positively correlated targets (Fig. 2C). All links are provided in Supplementary Table 2. Gene Ontology (GO) analysis results for these target genes revealed that up-regulated eRNAs in KIRC upregulated genes associated with immune cell activation, cell communication, and response to hypoxia pathways (Fig. 2D). The functions of these target genes reflect the main regulatory role of this part of eRNAs.

The genes with similar functions and potential regulatory relationships tend to be concentrically distributed. Therefore, we selected the modules with the greatest number of links and at least one eRNA to investigate the potential functions of the involved eRNAs. To further explore the functions of some key modules, we used SPICi to extract the two modules with the most edges. (Fig. 2E, H). The GO enrichment analysis results of module 1 centered on eRNA (ENSR00000045242) showed that this module mainly regulates hypoxia-related pathways (Fig. 2F). Many studies have shown that hypoxia can promote tumor cell proliferation and survival, increase tumor invasiveness, induce angiogenesis, reduce treatment effects, and inhibit immune responses [38,39,40]. The KEGG pathway enrichment results reflect its high correlation with tumors (Fig. 2G). The GO enrichment analysis results of module 2 clearly show the immune regulatory potential of this module (Fig. 2I). The KEGG enrichment analysis results of module 2 also show that it may regulate immune responses and cell signaling (Fig. 2J).

Down-regulated eRNAs in KIRC potentially regulate ion transport, cell cycle, and metabolism

For the down-regulated eRNAs, we conducted a co-expression network analysis between these eRNAs and positively correlated targets, revealing three distinct modules within the network (Fig. 3A). All links are provided in Supplementary Table 3. Subsequent Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the target genes unveiled that down-regulated eRNAs in KIRC upregulate the expression of genes associated with ion transport, cell cycle, and metabolism (Fig. 3B). Specifically, the GO analysis results of the module from subnetwork1 showed that this module was closely associated with transport (Fig. 3C, F). The GO analysis results of the module from subnetwork2 illustrated that the candidate targets in this module were associated with JUN kinase activity, cell cycle, and DNA integrity pathways (Fig. 3D, G). The module from subnetwork 3 is relatively simple (Fig. 3E). And the GO enrichment results showed that these genes were strongly related to metabolism-related pathways (Fig. 3H). In conclusion, we assumed that these down-regulated eRNAs in KIRC potentially regulate ion transport, cell cycle, and metabolism.

Fig. 3

Regulatory network of the down-regulated eRNAs in KIRC. A The co-expression network between down-regulated eRNAs and positively correlated targets; B The barplot of GO and KEGG enrichment analysis results of the positively correlated targets; C, D, E The three modules of subnetwork from the co-expression network; F, G, H The barplot of GO enrichment analysis results of these three modules

eRNA-centric network reveals the loop regulatory axis

To delve into the regulatory mechanism of eRNAs, we utilized data from the Cistrome database to establish a TF-eRNA-EP loop regulatory axis network. Figure 4A depicts the regulatory network centered around up-regulated eRNAs, while Fig. 4E illustrates the regulatory network in which down-regulated eRNAs participate. Leveraging 11,356 TF ChIP-seq datasets from the Cistrome database and 200 HiChIP datasets from HiChIPdb, we identified TF binding sites on eRNA regions and discovered H3K27ac eRNA-associated E-P loops. The resulting TF-eRNA and eRNA-EP loop relationships were merged to construct TF-eRNA-EP loop regulatory axes, further integrated with protein–protein interactions from the STRING database to form the eRNA-centric regulatory network. Figure 4B and C illustrated that up-regulated eRNA and TF/E-P loop have more regulatory pairs.

Fig. 4

The analysis of the TF-eRNA-EP loop regulatory axis network. A The TF-eRNA-EP loop regulatory network centered on up-regulated RNA; B The average number of TFs per up- and down-regulated eRNA connection; C The average number of E-P loops per up- and down-regulated eRNA connection; D Visualizing the eRNA (ENSR00000106541), its TF, the E-P loop that may be regulated by it, and the possible downstream target genes; E The TF-eRNA-EP loop regulatory network centered on down-regulated RNA; F The barplot of GO enrichment analysis results of up-regulated eRNA-centric regulatory network; G The barplot of GO enrichment analysis results of down-regulated eRNA-centric regulatory network; H The barplot of Msigdb hallmark enrichment analysis results of these two regulatory networks

Taking the example of an up-regulated eRNA (ENSR00000106541), TF binding sites were identified using the HEK293 cell line, strongly related to the kidney. Figure 4D shows that TFs such as ZNF692 and PRDM1 may affect the expression of ENSR00000106541, thereby affecting the strength of the E-P loop that mediates the formation of this eRNA and the promoter of CD70, thereby regulating the expression of the target gene (CD70). Interestingly, CD70 is closely related to the occurrence, development, and prognosis of KIRC [41]. In humans, CD70 is mainly expressed in activated lymphocytes [42]. Ruf et al. found that CD70 is highly expressed in RCC, of which KIRC expresses up to 78%, but normal renal cells do not express CD70 [43]. This may lead to immune function exhaustion and induce immune escape [44, 45]. In recent years, clinical studies of some novel chimeric antigen receptor (CAR) T cell therapies targeting CD70 in renal cancer patients have shown encouraging anti-tumor activity [46,47,48]. In summary, CD70 is expected to become a potential target with both efficacy and safety in the treatment of renal cancer. According to our research, the high expression level of CD70 in KIRC may be related to the up-regulated eRNA (ENSR00000106541). Therefore, this is a potential more upstream regulatory mechanism, which is very meaningful for our existing research on cancer regulatory genes.

Then, we did functional enrichment analyses on the genes involved in Fig. 4A and E. The results of the enrichment analyses showed that the pathways enriched by these genes are shown to be very related to DNA-binding transcription factors (Fig. 4F and G). Most of the genes used for enrichment are TFs, so this result is very reasonable. It was worth noting that the genes involved in the network constructed centered on the up-regulated RNAs were enriched in tumor necrosis factor receptor superfamily binding (Fig. 4F). Besides, the results of Msigdb hallmark enrichment analyses also show similar results. Among the genes enriched in the up-regulated-eRNA-centric network, the most significantly related one is TNFA signaling via NFKB, which is also significantly related to pathways such as glycolysis and hypoxia (Fig. 4H). These three pathways with the highest correlation significance were closely related to tumor progression.

The eRNA-based diagnosis prediction model

To further elucidate the clinical significance of these eRNAs in clear cell renal cell carcinoma (KIRC), we developed an accessible eRNA-related diagnostic model (Fig. 5A). Employing a univariate logistic regression model with eRNA expression levels as input variables, we identified differentially expressed (DE) eRNAs significantly associated with oncogenesis (Fig. 5B). Notably, 113 eRNAs exhibited a positive correlation with the occurrence of KIRC, suggesting that higher expression of these eRNAs is indicative of a higher likelihood of KIRC. Conversely, 62 eRNAs showed a negative correlation with KIRC occurrence, with the previously studied eRNA (ENSR00000305834) being among the negatively correlated factors. Subsequently, these eRNAs were utilized to construct a diagnostic prediction model employing a multivariate logistic regression model. The ROC curve illustrated that the diagnostic predictive model demonstrated high sensitivity and specificity, with an area under the curve (AUC) value of 0.834 (Fig. 5C). This reflects the effectiveness of this model.

Fig. 5

The KIRC-associated diagnosis model and eRNA-centric regulatory network. A Pattern diagram for model establishment; B The volcano plot of DE eRNAs significantly associated with tumorigenesis; C ROC curves of diagnosis eRNA factors (lambda. min = 29, lambda.1se = 4); D The TF-eRNA-EP loop regulatory network centered on up-regulated RNA, which only includes KIRC-annotated DEGs; E The barplot of the average number of edges of eRNA and overall in this network; F The most densely connected module from this up-regulated eRNA-centric regulatory network; G The barplot of GO enrichment analysis results of up-regulated eRNA-centric regulatory network; H The TF-eRNA-EP loop regulatory network centered on down-regulated RNA, which only include KIRC-annotated DEGs; I The barplot of the average number of edges of eRNA and overall in this network; J The most densely connected module from the this down-regulated eRNA-centric regulatory network; K The barplot of GO enrichment analysis results of down-regulated eRNA-centric regulatory network;

KIRC-specific eRNA-centric network

While the eRNA-centric networks in Fig. 4 (Fig. 4A and E) are constructed around differentially expressed eRNAs in KIRC samples, not all related TFs and E-P loops exhibit differential expression in KIRC. To create networks more closely related to KIRC, we integrated the network (Figs. 2C and 3A), including only differentially expressed eRNAs and positively correlated target genes.

The refined eRNA-centric networks (Fig. 5D and H) center around up-regulated and down-regulated eRNAs, respectively. Interestingly, in networks centered on up-regulated eRNAs, eRNAs were predominantly positioned at the network's edge, a pattern not observed in networks centered on down-regulated eRNAs. Upon analyzing the average number of edges of eRNAs, we found a lower average in the network centered on up-regulated eRNAs compared to the overall network (Fig. 5E), while the opposite was true for the network centered on down-regulated eRNAs (Fig. 5I). Further examination revealed that multiple up-regulated eRNAs often acted on a single E-P loop or TF in the network centered around up-regulated eRNAs (Fig. 5D), a scenario less common in the network centered around down-regulated eRNAs (Fig. 5H), indicating potential differences in the roles played by up-regulated and down-regulated eRNAs.

Utilizing SPICi, we extracted the two modules with the most edges (Fig. 5F and J) and conducted GO pathway enrichment analysis on the genes within these modules. The module centered around up-regulated eRNAs exhibited enrichment in DNA-binding transcription factors (Fig. 5G), while the downregulated eRNA-centric network module showed enrichment for genes associated with renal absorptive pathways (Fig. 5K). Despite the small size of this module (only three genes), it suggests a specific connection to the renal absorption pathway. These refined networks provide a nuanced understanding of the distinct roles played by up-regulated and down-regulated eRNAs in the context of KIRC.

The eRNA regulates SLC15A2 in KIRC

We use the module in Fig. 5J as an example for elaboration. Figure 6A illustrates that the expression of eRNA (ENSR00000305834) in cancer samples is significantly lower than in normal samples. Intriguingly, patients with high expression of ENSR00000305834 exhibit significantly better survival outcomes than those with low expression (Fig. 6B). The positive correlation between the expression of this eRNA and its target gene (SLC15A2) and transcription factor (TFAP2A) is evident (Fig. 6C and D). Additionally, the expression levels of both SLC15A2 and TFAP2A in cancer samples are significantly lower than in normal samples (Fig. 6E).

Fig. 6

eRNA (ENSR00000305834) regulatory mechanism. A The expression level of eRNA (ENSR00000305834); B The Kaplan–Meier survival curve of ENSR00000305834 in KIRC; C The contour map showed the correlation of the eRNA (ENSR00000305834) and SLC15A2; D The contour map showed the correlation of the eRNA (ENSR00000305834) and TFAP2A; E The expression level of SLC15A2 and TFAP2A; F Visualizing theENSR00000305834), its TF, and the E-P loop that may be regulated by it; G The mechanism diagram of eRNA (ENSR00000305834) regulating SLC15A2 expression. H Immunoblot detecting the expression of SLC15A2 in the ENSR00000305834 OE transduced 786-O cell line. I, J Quantitative real-time PCR (qRT-PCR) analysis of ENSR00000305834 and SLC15A2 levels in the TFAP2A OE transduced 786-O cell line (P values were determined by the paired student’s t-test). K Immunoblot detecting the expression of SLC15A2 in the ENSR00000305834 OE transduced 786-O cell line. L, M Colony formation assays and cell growth curves were performed to assess the proliferation of ENSR00000305834 OE transduced 786-O cell line (P values were determined by the one-way ANOVA, **P < 0.01, **** P < 0.0001)

Notably, the SLC15A family is believed to play a key role in the development and progression of many tumors due to their involvement in tumor metabolism [16]. The study of Huo, X. et al. elaborated on the specific distribution and physiological function of SLC15A2 in the kidney [21]. Multiple studies have also emphasized SLC15A2 as a potential prognostic biomarker and therapeutic target in kidney and liver cancers [16, 20,21,22].

Collectively, our findings suggest that the low expression of the TF (TFAP2A) of eRNA (ENSR00000305834) in KIRC samples reduces the eRNA's expression, diminishing the intensity of the E-P loop it mediates. This ultimately leads to decreased expression of the target gene (SLC15A2) (Fig. 6F and G).

Validation of eRNA regulation in vitro

To investigate the biological role of eRNA in tumor progression and its association with KIRC, we conducted a series of in vitro experiments. Initially, we performed overexpression experiments of TFAP2A, the TF of ENSR00000305834, in the human renal cancer cell line (786-O). Immunoblotting analysis depicted in Fig. 6H confirmed successful overexpression of TFAP2A. Subsequently, through PCR testing of eRNA (ENSR00000305834) and its downstream target gene (SLC15A2) in cell lines, we observed significant upregulation of both ENSR00000305834 and SLC15A2 in TFAP2A-overexpressing cells (Fig. 6I and J). These findings corroborate earlier bioinformatics analyses results. Then, to establish the regulatory influence of eRNA on downstream targets, we generated a 786-O cell line overexpressing ENSR00000305834. Western blot analysis demonstrated that heightened expression of ENSR00000305834 led to increased levels of the E-P loop target gene, SLC15A2 (Fig. 6K). Collectively, these results delineate the regulatory role of eRNA in the context of KIRC and underscore the reliability of this regulatory pathway.

Furthermore, colony formation assays revealed a significant reduction in the formation of tumor cell colonies upon overexpression of ENSR00000305834 (Fig. 6L), underscoring its pivotal role in modulating the clonal growth capability of kidney cancer cells. These findings were further supported by cell proliferation assays, which demonstrated a notable decrease in the growth rate of the renal cancer cell line following the overexpression of ENSR00000305834 (Fig. 6M). This highlights the significance of further exploring eRNA-based therapeutic strategies for managing KIRC and potentially other types of cancer.

Discussion

Over the last decade, numerous studies have linked enhancer RNA (eRNA) to disease occurrence and development [3]. However, the precise functions and mechanisms of eRNA in cancer progression remain unclear. The human transcriptome boasts an extensive number of eRNAs, with over 65,000 enhancers exhibiting eRNA transcripts, displaying specificity for tumor types and individual patients [3, 49,50,51]. This abundance suggests the potential for eRNAs as biomarkers and therapeutic targets.

In colorectal cancer, studies have revealed that eRNAs such as CCAT1 regulate chromatin looping, influencing oncogenic gene expression patterns [52]. Similarly, in prostate cancer, androgen receptor-regulated eRNAs play critical roles in driving tumorigenesis [37]. These findings demonstrate eRNAs' versatile involvement in transcriptional control across different cancers.

Clear cell renal cell carcinoma (KIRC), marked by high incidence and poor prognosis, served as a focal point for our exploration into the role of eRNA in renal cancer. Survival analysis revealed that nearly half of active eRNAs significantly influenced the prognosis of renal cancer patients. This highlighted the pivotal role of eRNA in kidney cancer development, particularly a substantial proportion of active eRNAs.

To delve deeper into the functions of these eRNAs, we constructed an eRNA regulatory network connecting them with target genes differentially expressed in cancer and normal samples. The findings indicated that up-regulated eRNAs in KIRC potentially regulate immune response and hypoxia, while down-regulated eRNAs potentially regulate ion transport, cell cycle, and metabolism.

To enhance the clinical relevance of our study, we developed a practical eRNA-based diagnostic model. What is gratifying is that the diagnostic model demonstrated notable effectiveness. This provides promising clues for KIRC's future diagnostic and therapeutic efforts.

Subsequent investigation into the regulatory mechanisms of eRNA in renal cancer involved the construction of eRNA-centric TF-eRNA-EP loop networks, revealing the TF-eRNA-EP loop regulatory axis. This model unveiled a potential regulatory mechanism of eRNA (ENSR00000305834) on an important biomarker (SLC15A2) in KIRC. The proposed model suggests that the expression of the TF of the eRNA influences the eRNA's expression level, impacting EP-loop formation and thereby regulating downstream target genes, exemplified by the regulation of SLC15A2 in KIRC.

Given the functional importance of SLC15A2 in renal physiology, our findings highlight its potential clinical implications. SLC15A2 encodes a solute carrier protein involved in cellular transport processes, which may be critical in tumor progression and metabolic reprogramming in KIRC. Targeting the regulatory axis of ENSR00000305834-SLC15A2 could represent a novel therapeutic strategy. Specifically, inhibiting the formation of TF-eRNA-EP loops may alter downstream oncogenic signaling, presenting opportunities for precision medicine. Furthermore, the association of SLC15A2 with cancer-specific pathways in our study underscores its potential as a diagnostic biomarker and a therapeutic target.

Although we have found how the expression of eRNA affects the occurrence and development of KIRC through the development of functional networks and diagnostic and prognostic models, as well as the important clinical significance of eRNA for KIRC. However, the limitations of this study must be acknowledged. This study offers foundational insights, yet further comprehensive research is warranted to achieve clinical applicability and identify specific, representative biomarkers for KIRC.