In the research, the gene expression profile of ccRCC and its corresponding normal samples were downloaded from International Cancer Genome Consortium (ICGC) and Cancer Genome Atlas (TCGA) database, including FPKM, count number expression profile data and clinical information. Then we converted FPKM to TPM based on the author’s R code. GENCODE’s gene annotation file (version 24) was used to differentiate between long noncoding RNAs and protein-coding genes (PCGs). We also transformed the Ensembl ID of the chosen genes into symbol form. All expression profiles have been processed before analysis. Also, 450k microarray chip data were downloaded from the TCGA-KIRC project. Besides, all the missing values in renal cancer data were completed with the KNN method. Based on the Gene Expression Omnibus (GEO) database, the GSE86091 dataset with para-cancerous and tumor samples was downloaded, including h3k4me1, h3k4me3, h3k27ac histone information (Yao et al., 2017).

Then, the lncRNAs and PCGs which expressed differentially were identified by R package deseq2 with p-value determined by the Benjamin Hochberg method. LncRNAs or PCGs with FDR < 0.05 and |logFC| > 1 were considered significant. Next, we selected the unique peak in RCC in accordance with the physical location of the histone modified peak. Only the peak with a p-value < 0.05 was retained as the difference peak. Furthermore, we combined with the GTF file of GENCODE to get differentially histone-modified genes as well as identity of gene enhancers from the human enhancer database FANTOM5. The promoter of the gene is defined as 2 kb upstream and 0.5 kb downstream of the transcription initiation site (TSS), which can be identified by the R package ChIPseeker and the human lincRNAsTranscripts database. Finally, Using R-packet CHAMP, the Bumphunter method is used to identify DMR and the regions with Bumphunter DMR. p-values < 0.01 are considered significant DMRs. In the end, we made a definition that lncRNAs meet the following conditions as epiRNAs: i) there was differential expression between tumor and normal samples; ii) promoter or enhancer has more than one histone modified or methylated region.

To figure out the genomic characteristics of epigenetic dysregulation and none dysregulation lncRNA/PCG, we made comparations among the number and length of exons, transcripts, genes of epi-lncRNA, non-epi-lncRNA, epi-PCG, non-epi-PCG.

To investigate the epigenetic features of lncRNA caused by histone modification, we analyzed the distribution characteristics of promoter and enhancer among epi-lncRNA modified through different histones in genome level.

To learn more about the relation between epigenetic disorders lncRNA and carcinoma, we made an analysis of epi-lncRNA enrichment and cancer lncRNAs with lnc2cancer v3.0. Besides, epi-lncRNA related to ccRCC was screened out for the sake of exploring the regulation of epi-lncRNA. In addition, considering that most of the genes that affect the process of the diseases tend to show expression disorder, we identified the epi-lncRNA in ccRCC through the TCGA TPM expression profile.

Pathway enrichment and quantification of specific pathways activity were completed by Gene Set Enrichment Analysis (GSEA) and Single Sample GSEA (ssGSEA) algorithm, including the Kyoto Encyclopedia of Genes and Genomes (KEGG) and epigenetically dysregulated lncRNAs score (Subramanian et al., 2005; Pang Y et al., 2022).

The Kaplan-Meier method was applied to estimate the overall survival (OS) rate and survival time, and based on the method of log-rank test, the differences among different sample subgroups were obtained. Finally, the potential prognostic markers of epi-lncRNA were identified. A risk model was built on basis of multivariate Cox regression analysis.

The CIBERSORT algorithm (version 1.03) and the LM22 signature matrix were used to analyze immune infiltration according to the expression profile (Chen et al., 2018).

The Genomics of Drug Sensitivity in Cancer Database (GDSC) was used to analyze drug susceptibility (Yang et al., 2013). The evaluation of patients taken immunotherapy was made by application of the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (Fu et al., 2020).

All analysis was performed with the programs of R software (version 4.2.1). Statistical significance was determined using two-tailed p-value and less than 0.05 was considered statistically significant.

This study was approved by the Institutional Review Board of the First Affiliated Hospital of Nanjing Medical University. All patients participated in the study voluntarily donated samples for clinical research and signed informed consents before surgery. We collected samples from patients with renal cell carcinoma from March 2024 to June 2024. All the samples were from patients took partial nephrectomy in the Department of Urology Department of the First Affiliated Hospital of Nanjing Medical University. Samples were collected under the guidance of one experienced pathologist, and immediately stored in a liquid nitrogen tank.

TRIZOL reagent (Invitrogen) was used to extract total RNA in both tissues and cells according to the manufacturer’s instruction. The cDNA was obtained from the purified RNA using a PrimeScript RT Reagent Kit (Takara). SYBR Premix Ex Taq II Kit (Takara) was used for qRT-PCR assays. Results were normalized to β-actin expression. All experiments were performed at least three times. All specific primers used in this study have been listed in Supplementary Table S1.

Genes that are differentially expressed were identified using DESeq2, including 5906 PCGs and 5370 lncRNAs. The combination of histone modifications and 450K methylation microarray data enabled us to identify 149 epiRNAs, 13426 non-epiRNAs, 2241 epi-PCGs, and 17200 non-epi-PCGs. Moreover, lncRNAs in ccRCC are significantly less frequently than in PCG (Figure 1A). A comparison of epi-lncRNA transcripts and exons with non-epi-lncRNA transcripts, epi-PCG transcripts, and non-epi-PCG transcripts was followed to represent the genomic characteristics of epigenetically dysregulated lncRNAs. Epi-lncRNAs appear with a long transcript length and large numbers of transcripts (Figures 1B, C) and epi-PCG showed similar structural characteristics as epi-lncRNAs. Simultaneously, epi-lncRNAs were presented with more exons and epi-PCG owns more exons and length (Figures 1D, E).

Systematic analysis of epigenetic disorders of epi-lncRNA in renal cell carcinoma revealed the landscape of epi-lncRNAs resulted from various histone modifications and different methylation regions (Figure 2A). Besides, in these lncRNAs, we found several abnormal histone modifications including H3K4me3, H3K4me1, and H3K27ac. Also, most of the lncRNAs with apparent disorders coexist with various histone modifications. Moreover, we also found that the abnormal histone modifications were mostly located in the promoter region (Figure 2B).

Next, we completed an analysis of the relation between the expression of epi-lncRNA and the pathway in ccRCC. The ssGSEA was applied to calculate each sample enrichment score in all the lncRNAs affected by histone modifications. We observed that the scores in cancer samples were remarkably higher than adjacent samples, suggesting that these epi-lncRNAs may own a great possibility of being carcinogenic (Figure 3A). Further, we analyzed the relationship between these enrichment scores and OS and found that high scores were associated with poor prognosis (Figure 3B). Additionally, we observed that there are 18 pathways in the KEGG pathway related to six sorts of epi-lncRNA scores closely, suggesting that there is a certain consistency among various epi-lncRNA-related pathways. Of these 18 pathways, the tumor-related pathway of BLADDER_CANCER and the immune-related pathways of PRIMARY_IMMUNODEFICIENCY, CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION, NATURAL_KILLER_CELL_MEDIATED_CYTOTOXICITY and AUTOIMMUNE_THYROID_DISEASE have Positive correlations to the epi-lncRNAs score. On the contrary, the metabolism-related pathways of SELENOAMINO_ACID_METABOLISM, CITRATE_CYCLE_TCA_CYCLE, PYRUVATE_METABOLISM and PROPANOATE_METABOLISM are negatively correlated with the epi-lncRNAs score (Figure 3C). These results reveal that epi-lncRNA plays an important role in tumor development and immunity.

We calculated the association of the enrichment scores of six epi-lncRNAs and m6a, m5c and m1A genes, indicating that they are significantly related. The enhancer and promoter of H3K4me3, H3K4me1 and H3K27ac not only have a similar correlation with these genes but also have their unique correlation (Figure 4A). The results indicate that there might exist unknown regulatory modes of epi-lncRNA resulted from enhancer and promoter histone modification, with a tight link between epi-lncRNA and RNA modification.

Using Lnc2Cancer V3.0, it can be discovered that 51 epi-lncRNAs were associated with normal human cancers (Figure 4B). To explore how candidate lncRNAs play the regulatory role in ccRCC, the list of ccRCC-related lncRNAs was collected from Lnc2Cancer v3.0. Of 51 epi-lncRNAs, 12 lncRNAs were identified that have direct relationships with ccRCC. We subsequently calculated the differences of 12 epi-lncRNAs between ccRCC samples and normal samples, respectively (Figure 4C). Of that, 11 epi-lncRNAs are presented with significant differences, except H19. Except that the expression of HOXA11-AS and LINC00461 in normal samples is higher than that in tumor samples, the other 9 lncRNAs are highly expressed in tumor samples. Among these epi-lncRNAs, HOTAIR, HOTTIP, XIST, H19, HCP5, and LINC00160 were only related to histone-modified promoters; LINC00461 is only related to DMR enhancers; PVT1 is associated with histone modified enhancer and promoter (Figure 5A). Then, we collected a clinical ccRCC cohrt (normal group: 8, disease group: 8) and performed PCR validation on the relevant lncRNAs, and found significant differences between the disease group and the normal group (Figure 5B).

Considering the potential prognostic value of epi-lncRNA is unclear yet, an analysis of the 12 epi-lncRNAs survival mentioned above based on the sample in TCGA-KIRC was made. Results indicated that 8 epi-lncRNAs are significantly associated with survival based on the univariate Cox regression analysis (Figure 6A). Kaplan-Meier survival curves illustrated that the epi-lncRNAs HOTAIR, HOTTIP, HOXA11-AS1, LINC00461, LUCAT1, PVT1, XIST, H19, LINC00160, MIAT, GAS6-AS1 were correlated with poor prognosis, while HCP5 was correlated with better prognosis (Figures 6B–M).

Multivariate survival analysis was used to establish the 8-EpiLncRNA signature on basis of the TCGA dataset, in which patients with high-risk score had higher mortality (Figure 7A). Receiver operating characteristic curve (ROC) analysis shows that the designed signature owned a remarkable prediction ability of patient OS (Figure 7B, 1-year AUC of 0.75, 3-year AUC of 0.67, and 5-year AUC of 0.72). Similarly, Kaplan-Meier survival curves indicated that the high-risk score patients have a worse OS performance tendency (Figure 7C). Then, we selected the RECA-EU with prognostic information from the ICGC database for model validation. In the ICGC cohort, we also observed satisfactory performance (Figure 7D). The ROC result of 1-, 3- and 5- years were 0.75, 0.68, and 0.68, respectively (Figure 7E). Kaplan-Meier survival curves indicated that the high-risk score ones also incline to have poor prognosis (Figure 7F). Univariate and multivariate analysis suggested that the new model we built can be considered as an independent risk factor of other clinical features (Supplementary Figures S1A, B).

Next, we investigated the fundamental distinction in Supplementary Figures S1A, B drug sensitivity between patients at high and low risk. Results indicated patients with high-risk score might respond better to axitinib and nilotinib, but not sorafenib and sunitinib (Figure 8). Moreover, A positive association was found between the TIDE score and the risk score. (Supplementary Figure S2A, cor = 0.274, p < 0.001). Meanwhile, it was worth noted that immunotherapy non-responders had a much higher score (Supplementary Figure S2B).

Based on the Hallmark gene collection, GSEA analysis revealed that the pathways of epithelial-mesenchymal transition, P53 pathway, UV response up, E2F targets, myogenesis, pancreas beta cells, and G2M checkpoint were triggered in patients at high risk (Supplementary Figure S3). Analysis of immune infiltration revealed that there’s a chance that the high-risk patients have more plasma cell infiltration, activated CD4⁺ memory T cells, M0 macrophages, follicular helper T cells, and regulatory T cells (Tregs), while a lower level of monocytes, resting mast cells and eosinophils (Supplementary Figure S4).