Abstract Clear cell renal cell carcinoma (ccRCC) is a malignant tumor originating from the renal tubular epithelium. Although the microRNAs (miRNAs) transcriptome of ccRCC has been extensively studied, the role of miRNAs editing in ccRCC is largely unknown. By analyzing small RNA sequencing profiles of renal tissues of 154 ccRCC patients and 22 normal controls, we identified 1025 miRNA editing sites from 246 pre-miRNAs. There were 122 editing events with significantly different editing levels in ccRCC compared to normal samples, which include two A-to-I editing events in the seed regions of hsa-mir-376a-3p and hsa-mir-376c-3p, respectively, and one C-to-U editing event in the seed region of hsa-mir-29c-3p. After comparing the targets of the original and edited miRNAs, we found that hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u had many new targets, respectively. Many of these new targets were deregulated in ccRCC, which might be related to the different editing levels of hsa-mir-376a-3p, hsa-mir-376c-3p, hsa-mir-29c-3p in ccRCC compared to normal controls. Our study sheds new light on miRNA editing events and their potential biological functions in ccRCC. Subject terms: Renal cell carcinoma, Renal cell carcinoma, Computational biology and bioinformatics Introduction Renal cell carcinoma (RCC), which is formed by malignant proliferation of tubular epithelial cells, is the most common malignant tumor of kidney^[42]1. According to the histological classification, there are many subtypes of RCC, among which clear cell renal cell carcinoma (ccRCC) is the most prevalent subtype, accounting for about 75% of all RCCs^[43]1. Although some ccRCC cases can be surgically resected, the metastatic rate of ccRCC is high, and about 30% of patients who present with metastases at first screening are not candidates for surgery^[44]2,[45]3. Therefore, early detection and targeted therapy of ccRCC is the most effective way to reduce the number of deaths from this disease^[46]4. MicroRNAs (miRNAs) are a class of single-stranded non-coding RNAs with approximately 22 nucleotides that can undergo extensive post-transcriptional modifications^[47]5,[48]6. As important regulatory molecules in biological processes, miRNAs are involved in almost all cellular pathways and pathological processes, including cancer initiation, progression and metastasis^[49]7. Mature miRNAs pair with their target mRNAs through the seed regions (the first eight nucleotides from the 5 [MATH: ′ :MATH] ends of miRNAs) to induce mRNA degradation or translational repression^[50]8. The regulation of gene expression by miRNAs is diverse and complex. A miRNA can bind to many different targets, and a target can be regulated by multiple miRNAs^[51]9. Therefore, even a single nucleotide change in a miRNA, especially if it occurs within the seed region, leads to severe changes of its targets and also affects the expression of its targets. In ccRCC, miRNAs are reported as either oncomiRNAs or tumor suppressors^[52]10. For example, hsa-miR-21-5p was overexpressed in ccRCC tissues compared to normal controls and correlated with the downregulation of PPAR- [MATH: α :MATH] in ccRCC^[53]11. Ji et al. reported that hsa-miR-155-5p was upregulated in ccRCC, and inhibition of hsa-miR-155-5p could significantly suppress the proliferation, colony formation, migration and invasion and induce G1 phase arrest and apoptosis^[54]12. In addition, hsa-miR-193a-3p and hsa-miR-224-3p as oncogenic miRNAs^[55]13, activated the PI3K/Akt pathway and targeted glycosylation-related enzymes to mediate cell proliferation, migration and invasion in ccRCC. hsa-miR-30a-5p^[56]14 and hsa-miR-30c-5p^[57]15 were downregulated in ccRCC, and were associated with tumor aggressiveness^[58]16. Cui et al. found that hsa-miR-99a-5p was downregulated in ccRCC and correlated with overall survival^[59]17. Moreover, hsa-miR-187-3p presents lower levels in tumor tissue and plasma of patients with ccRCC, and lower levels of hsa-miR-187-3p were associated with higher tumor grade and stage^[60]18. In summary, a large number of miRNAs play crucial roles in the initiation and progression of ccRCC. RNA editing is a prevalent and conserved post-transcriptional mechanism that plays a crucial role in the diversification of gene expression and transcriptome complexity across various organisms. This process involves specific proteins that catalyze chemical modifications to RNA molecules during their generation processes, leading to the alterations, deletions, and/or insertions of bases^[61]19,[62]20. RNA is edited in a variety of ways. The most widely reported editing type is adenosine to inosine (A-to-I) editing, catalyzed by adenosine deaminase acting on RNA (ADAR) enzymes^[63]19. This is followed by editing of cytosine (C) to uracil (U), which is catalyzed by the apolipoprotein B mRNA editing catalytic polypeptide-like (APOBEC) proteins^[64]19,[65]20. Furthermore, some RNAs are modified on their 3 [MATH: ′ :MATH] ends to add additional nucleotides, such as A or U^[66]21–[67]23. It has been found that a large number of RNA editing events also exist in miRNAs^[68]24–[69]37 and the A-to-I miRNA editing events in cancer have attracted some attention^[70]38–[71]40. For example, ADAR2 is responsible for editing hsa-miR-21-3p/5p, hsa-miR-221-5p and hsa-miR-222-3p, which exhibit high levels and are known to promote cancer progression in glioblastoma^[72]41. In a study on miRNA editing in pan-cancer, the edited hsa-miR-200b-3p promotes cell invasion and migration by targeting LIFR, a metastasis suppressor^[73]39. A-to-I edited hsa-miR-589-3p targets ADAM12 and original hsa-miR-589-3p targets PCDH9, a tumor suppressor related to glioma progression^[74]42. A-to-I editing level of hsa-miR-589-3p is reduced in glioblastoma, which promotes the proliferation and invasion of brain cancer cells^[75]42. In addition, miRNA editing also plays an important role in neurological diseases. For example, A-to-I editing of hsa-miR-497-5p is enhanced in Parkinson’s disease (PD), which potentially promotes progressive neurodegeneration of PD patients^[76]37. In summary, miRNA editing is associated with the initiation and progression of numerous cancers. However, miRNA editing events in ccRCC are largely unknown, which has become a major obstacle and blind spot for the study of ccRCC. In this study, to comprehensively characterize miRNA editing and/or modification sites in ccRCC, we systematically analyzed small RNA sequencing profiles of 154 ccRCC patients and 22 healthy controls. We identified 1025 miRNA mutation and/or editing (M/E) sites with significant editing levels in these samples. In addition, 122 M/E sites showed significantly different editing levels in ccRCC patients compared to control samples. Among them, we focused on 3 editing events that occurred in the seed regions and predicted the target genes of their corresponding edited miRNAs. We explored the potential functions of new target genes of edited miRNAs and examined expression patterns of these new targets in ccRCC using existing gene expression data. Our results suggested that A-to-I edited hsa-mir-376a-1-3p and hsa-mir-376c-3p targeted BCAT1 and MTHFD2, respectively, and C-to-U edited hsa-mir-29c-3p targeted RASSF8. Presumably, because the A-to-I editing levels of hsa-mir-376a-1-3p and hsa-mir-376c-3p were reduced in ccRCC, BCAT1 and MTHFD2 were upregulated in ccRCC. Similarly, potentially due to increased C-to-U editing level of hsa-mir-29c-3p, RASSF8 was significantly downregulated in ccRCC. These results suggest that miRNA editing contributes to initiation and progression of ccRCC. Results Summary of mutation and editing sites identified in miRNAs The MiRME pipeline^[77]43 was used to identify miRNA M/E sites from the collected 176 sRNA-seq profiles. In total, 1025 M/E sites with significant editing levels were identified (as listed in Supplementary Table [78]S2.1). As mentioned above, these M/E sites were divided into 9 categories (Fig. [79]1A and Supplementary Table [80]S2.2). We found that there were more 3 [MATH: ′ :MATH] -addition sites than 5 [MATH: ′ :MATH] editing sites (Fig. [81]1A and Supplementary Tables [82]S2.3, [83]S2.4). In addition, M/E sites in the central regions of mature miRNAs include 9 (0.88%) A-to-I, 30 (2.93%) C-to-U, 5 (2.54%) Other and 2 (0.20%) SNP sites. Figure 1. [84]Figure 1 [85]Open in a new tab The number of significant M/E sites in miRNAs and their categories in the analyzed sRNA libraries. (A) The categories of significant M/E sites in miRNAs. (B) The editing sites and numbers in the central region of mature miRNAs. (C) The distribution of the numbers of pre-miRNAs with different numbers of [MATH: 3′ :MATH] -, [MATH: 5′ :MATH] -editing and Central editing sites. We further investigated the miRNA editing events that occured in the central regions of mature miRNAs (Supplementary Table [86]S2.5), i.e., A-to-I, C-to-U and Other. We counted the base changes of these three editing types (Fig. [87]1B), and found 9 A-to-I, 30 C-to-U, and 5 Other editing events. One insertion and two deletion events were also found (Fig. [88]1B). Furthermore, we also investigated the distributions of these events within different regions of mature miRNAs, including the 5 [MATH: ′ :MATH] ends, 3 [MATH: ′ :MATH] ends, and central regions. Our results indicated that 3 [MATH: ′ :MATH] -addition was the most common type of RNA editing events of miRNAs. As shown in Fig. [89]1C, the numbers of editing events at the 3 [MATH: ′ :MATH] ends were significantly higher than those at the 5 [MATH: ′ :MATH] ends and the central regions. For the vast majority of miRNAs, only 1 or 2 editing events happened at the 5 [MATH: ′ :MATH] ends and the central regions. Nine miRNA A-to-I editing sites in ccRCC Totally, 9 A-to-I editing sites with significant editing levels were identified from the 176 samples (Fig. [90]2A). The editing levels of hsa-mir-376a-1 and hsa-mir-376a-2 were almost identical, probably due to their high sequence similarity. We then counted the nucleotides around these 9 A-to-I editing sites, and found that the 5 [MATH: ′ :MATH] and 3 [MATH: ′ :MATH] ends of these A-to-I editing sites had clear U and G base preference, respectively (Fig. [91]2B), which was consistent with results in literature^[92]27,[93]29,[94]43. As examples, the detailed information of two A-to-I editing sites in hsa-mir-411 and hsa-mir-7977 in two ccRCC samples (SRR11873730 and ERR4367258, respectively) were presented in Fig. [95]2C and D, respectively. As shown in Fig. [96]2E and F, a large number of sequencing reads in these samples supported the two A-to-I editing sites. Figure 2. [97]Figure 2 [98]Open in a new tab The details of 9 identified A-to-I editing sites in miRNAs. (A) The editing levels of the 9 A-to-I editing sites in the 176 selected sRNA profiles. (B) The percentages of nucleotides around of 9 A-to-I sites. (C) The MiRME map of hsa-mir-411 in one of the ccRCC samples selected (SRR11873730). (D) The MiMRE map of hsa-mir-7977 in one of the ccRCC samples selected (ERR4367258). (E) The details of hsa-mir-411_20_A_g in SRR11878730. (F) The details of hsa-mir-7977_6_A_g in ERR4367258. In Part (C) to (D), The three panels from top to bottom show the total number of reads (Tags Per Ten Million sequencing reads, tptm), the number of M/E reads with variations (tptm), and the p-values of M/E events, respectively. In Part (E) to (F), the blue nucleotides are original and edited. And the three numbers on the right indicate the frequency of the read, the length of the read, and the weight of the read at this locus, respectively. Thirty miRNA C-to-U editing sites in ccRCC We carefully studied the 30 C-to-U editing sites and found that their editing levels were generally low in the 176 sample (Fig. [99]3A). Furthermore, nucleotides on the 5 [MATH: ′ :MATH] sides of these C-to-U editing sites were biased toward C, while those on the 3 [MATH: ′ :MATH] sides were biased toward A (Fig. [100]3B). Examples of two C-to-U editing sites were shown in Fig. [101]3C and D, respectively. In addition, the detailed editing sites and editing levels of these two miRNAs were shown in Fig. [102]3E and F, respectively. These two C-to-U editing sites were supported by thousands of sequencing reads in the samples used, indicating good reliabilities of these two editing sites. Figure 3. [103]Figure 3 [104]Open in a new tab The details of 30 identified C-to-U editing sites in miRNAs. (A) The editing levels of the 30 C-to-U sites in the 176 selected sRNA-seq profiles. (B) The percentages of nucleotides around the 30 C-to-U sites. (C) The MiRME map of hsa-mir-126 in one of the ccRCC samples selected (ERR4367258). (D) The MiRME map of hsa-mir-451a in one of the ccRCC samples selected (ERR4367208). (E) The details of hsa-mir-126_58_C_u in ERR4367258. (F) The details of hsa-mir-451a_21_C_u in ERR4367208. In Part (C) to (D), The three panels from top to bottom show the total number of reads (Tags Per Ten Million sequencing reads, tptm), the number of M/E reads with variations (tptm), and the p-values of M/E events, respectively. In Part (E) to (F), the blue nucleotides are original and edited. And the three numbers on the right indicate the frequency of the read, the length of the read, and the weight of the read at this locus, respectively. Two miRNA SNP sites in ccRCC The identified M/E sites were compared with SNPs in miRNAs previously reported^[105]44 and NCBI dbSNP (v151) according to the criteria mentioned in Materials and Methods. Two SNPs (rs11614913 and rs2155248) were shown in Supplementary Fig. [106]S1 and Supplementary Table [107]S2.6, and the editing levels of these two SNPs in one normal and one ccRCC sample, respectively, were shown in Supplementary Fig. [108]S1A. As shown in Supplementary Fig. [109]S1B and [110]S1C, these SNPs were supported by many sequencing reads in SRR11873716 and SRR11873719, respectively. Additionally, the editing levels of these two sites reached 100% (Supplementary Fig. [111]S1D and [112]S1E), consistent to the fact that these sites were SNPs. 122 M/E sites with significantly different editing levels in ccRCC We obtained 122 miRNA M/E sites with significantly different editing levels in ccRCC samples (n = 154) compared to normal controls (n = 22) by using Mann-Whitney U-tests (corrected [MATH: p<0.05 :MATH] ) (Supplementary Table [113]S3). The 122 miRNA M/E sites were divided into 7 categories. 112 of these 122 sites were 3 [MATH: ′ :MATH] -addition events (3 [MATH: ′ :MATH] -A, 3 [MATH: ′ :MATH] -U, and 3 [MATH: ′ :MATH] -Other). In addition, there were 5 A-to-I editing sites, 2 C-to-U editing sites, and the remaining two were 5 [MATH: ′ :MATH] end sites and one Pseudo site, i.e., potential false positive site (Fig. [114]4A). Most of the 122 M/E sites had increased editing levels (77 M/E sites, 63.1%) and 45 M/E sites (36.9%) had decreased editing levels in ccRCC samples when compared to normal controls (Fig. [115]4B). The categories of M/E sites with increased and decreased editing levels in ccRCC samples were shown in Fig. [116]4C. In M/E sites with increased editing levels in ccRCC, the 3 [MATH: ′ :MATH] -A editing events accounted for the largest portion, while the 3 [MATH: ′ :MATH] -U editing events accounted for the largest part in those with decreased editing levels in ccRCC compared to normal controls. Figure 4. [117]Figure 4 [118]Open in a new tab The 122 miRNA M/E sites that have significantly different editing levels in ccRCC samples when compared to normal samples. (A) The categories of the 122 M/E sites. (B) The change of editing levels of 122 M/E sites in ccRCC samples compared to normal samples. (C) The percentages of different categories of 122 M/E sites with significantly different editing levels in ccRCC samples. (D) Comparison of editing levels of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-6503_59g in ccRCC and normal samples. (E) Comparison of editing levels of hsa-mir-29c_59u in ccRCC and normal samples. (F) Comparison of expression levels of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-6503_59g in ccRCC and normal samples. (G) Comparison of expression levels of hsa-mir-29c_59u in ccRCC and normal samples. (H) The MiRME map of hsa-mir-376a-1 in one of the ccRCC samples (SRR11873721). (I) The MiRME map of hsa-mir-376c in one of the normal samples (SRR11873716). (J) The MiRME map of hsa-mir-29c in one of the normal samples (ERR4367209). (K) The details of hsa-mir-376a-1_49_A_g in SRR11873721. (L) The details of hsa-mir-376c_48_A_g in SRR11873716. (M) The details of hsa-mir-29c_59_C_u in ERR4367209. In Part (D) to (G), “*”: [MATH: p< :MATH] 0.05; “**”: [MATH: p< :MATH] 0.01; “***”: [MATH: p< :MATH] 0.001 and “NS”: not significant, i.e., [MATH: p> :MATH] 0.05. In ccRCC, the editing levels of hsa-mir-376a-1_49_A_g and hsa-mir-376c_48_A_g were significantly decreased (Fig. [119]4D), while those of hsa-mir-6503_59_A_g and hsa-mir-29c_59_u were significantly increased (Fig. [120]4D and E, respectively). We also investigated the levels of these edited miRNAs and found that hsa-mir-376a-1_49g and hsa-mir-376c_48g were significantly downregulated in ccRCC samples, while the levels of hsa-mir-6503_59g and hsa-mir-29c_59u were significantly upregulated in ccRCC (Fig. [121]4F and G, respectively). Figure [122]4H–M illustrated the details of the three edited miRNAs, hsa-mir-376a-1_49_A_g, hsa-mir-376c_48_A_g and hsa-mir-29c_59_u in one ccRCC sample (SRR11873721), one normal sample (SRR11873716) and one ccRCC sample (ERR4367209), respectively. The expression levels of ADAR and APOBEC genes in ccRCC As shown in Supplementary Fig. [123]S2 and Supplementary Table [124]S7, the expression levels of the ADAR genes that catalyzed A-to-I editing of miRNAs were examined. In the TCGA data, ADAR1 was downregulated in ccRCC compared to normal controls (Supplementary Fig. [125]S2A) which was consistent with the decreased editing levels of hsa-mir-376a-1_49_A_g, hsa-mir-376a-2_55_A_g, and hsa-mir-376c_48_A_g. However, in the other two cohorts of gene expression profiles, ADAR1 was upregulated, suggesting that more researches were needed to reveal the role of ADAR1 in the editing of miRNAs in ccRCC. We found that ADAR2 was significantly upregulated in ccRCC samples compared to normal controls in the three selected cohorts of gene expression datasets (Supplementary Fig. [126]S2B). Meanwhile, the editing level of hsa-mir-6503_59_A_g was significantly increased in ccRCC samples compared to normal controls (Fig. [127]4D), suggesting that ADAR2 mediated the editing of hsa-mir-6503 in ccRCC and contributed to its increased editing levels in ccRCC compared to normal controls. ADAR3 was downregulated in one cohort of gene expression profiles, and had no significant change in the other two datasets (Supplementary Fig. [128]S2C). The expression of the APOBEC genes that mediated C-to-U editing of miRNAs was also examined. As shown in Supplementary Fig. [129]S3 and Supplementary Table [130]S8, the expression levels of most APOBEC family members showed a significantly upward trend in ccRCC tumors, including AICDA, APOBEC3A, APOBEC3B, APOBEC3C, APOBEC3D, APOBEC3F, APOBEC3G and APOBEC3H genes (Supplementary Fig. [131]S3A, D–J, respectively). In our results, the editing level of hsa-mir-29c_59_C_u was significantly increased in ccRCC samples (Fig. [132]4E). More researches are needed to clarify which APOBEC gene is the most promising enzyme that mediates the C-to-U editing of hsa-mir-29c_59_C_u and other miRNA C-to-U editing sites in ccRCC. The expression levels of TENT genes in ccRCC Terminal nucleotidyltransferases (TENTs) modify RNA at post-transcriptional level, which regulates the stability and activity of RNA^[133]45. TENT2 and TENT4B proteins catalyze 3 [MATH: ′ :MATH] end adenylation of miRNAs^[134]21,[135]46–[136]50. TUT4 and TUT7 proteins are terminal uridine transferases that belong to the TENT3 subfamily, which have the ability to catalyze the 3 [MATH: ′ :MATH] uridylation of miRNAs^[137]23,[138]51–[139]53. The levels of these TENT genes in ccRCC and normal controls were examined (Supplementary Table [140]S9). In ccRCC, the level of TENT2 was significantly higher (Supplementary Fig. [141]S4A–C), while TENT4B was significantly downregulated (Supplementary Fig. [142]S4A and C). There are more 3 [MATH: ′ :MATH] -A sites with increased editing levels in ccRCC than 3 [MATH: ′ :MATH] -A sites with decreased editing levels in ccRCC (Fig. [143]4C), suggesting that TENT2 might be the major enzyme for the 3 [MATH: ′ :MATH] -A editing of miRNAs in ccRCC, while the downregulation of TENT4B could be responsible for the decreased 3 [MATH: ′ :MATH] -A editing levels of some miRNAs in ccRCC. The levels of TUT4 had no significantly different expression in ccRCC compared to normal controls (Supplementary Fig. [144]S4A to C), but the expression of TUT7 was significantly upregulated in ccRCC (Supplementary Fig. [145]S4A–C). These results suggest that some 3 [MATH: ′ :MATH] -U editing sites of miRNAs with increased editing levels in ccRCC (Fig. [146]4C, left column) are mainly catalyzed by TUT7. Identification of targets for A-to-I and C-to-U edited miRNAs Three M/E sites (hsa-mir-376a-1_49_A_g and hsa-mir-376c_48_A_g and hsa-mir-29c_59_C_u) in seed regions of mature miRNAs were selected for subsequent analysis, because these sites had significantly different editing levels in ccRCC and the levels of their corresponding edited miRNAs were also significantly different and in the same directions of changes when compared their editing levels in ccRCC with normal controls. To decrease false positives and increase the reliabilities of predicted miRNA targets, we used 11 PAR-CLIP sequencing profiles. The principle of PAR-CLIP (Photoactivatable Ribonucleoside-Enhanced Crosslinking and Immunoprecipitation) technology is used to detect the crosslinking of RNA-binding proteins (RBPs) to RNA molecules, followed by immunoprecipitation and sequencing to identify the binding sites^[147]54. Through PAR-CLIP experiments, the interaction between RISC (including AGO proteins and miRNAs) and mRNA can be identified^[148]54. The sequences of the miRNA binding sites were obtained through PAR-CLIP experiments (targeting one of the AGO proteins) and then were sequenced^[149]54. Next, the MiCPAR algorithm^[150]55 was used to identify the targets of original and edited miRNAs by analyzing 11 PAR-CLIP profiles of AGO proteins (see Methods for details). We found that the original and edited hsa-mir-376a-1-3p had 102 common target genes, and 484 new target genes of hsa-mir-376a-1_49g were found (Fig. [151]5A and Supplementary Tables [152]S4.1–[153]S4.3). The original and edited hsa-mir-376c-3p had 107 common target genes, and hsa-mir-376c_48g had 568 new target genes (Fig. [154]5B and Supplementary Tables [155]S5.1–[156]S5.3). The original and edited hsa-mir-29c-3p had 107 identical target genes, and hsa-mir-29c_59u had 480 new target genes (Fig. [157]5C and Supplementary Tables [158]S6.1–[159]S6.3). Next, we compared the new target genes of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u to deregulated genes in ccRCC. We found that in the three gene expression datasets selected ([160]GSE151419, [161]GSE126964 and TCGA), 88 and 82 target genes of hsa-mir-376a-1_49g and hsa-mir-376c_48g were commonly significantly upregulated in ccRCC (Fig. [162]5D,E and Supplementary Tables [163]S4.4, [164]S5.4), respectively. Additionally, 51 target genes of hsa-mir-29c_59u were commonly significantly downregulated in ccRCC (Fig. [165]5F and Supplementary Table [166]S6.4). Figure 5. [167]Figure 5 [168]Open in a new tab Comparisons of targets of edited miRNAs and deregulated genes in clear cell renal cell carcinoma. Venn diagrams were prepared with Venny (https://bioinfogp.cnb.csic.es/tools/venny/). (A–C) The venn diagrams depict the targets of the original and edited hsa-mir-376a-1, hsa-mir-376c and hsa-mir-29c, respectively. (D) The overlap between the new target genes of hsa-mir-376a-1_49g and genes upregulated in ccRCC in three gene expression datasets. (E) The overlap between the new target genes of hsa-mir-376c_48g and genes upregulated in ccRCC in three gene expression datasets. (F) The overlap between the new target genes of hsa-mir-29c_59u and genes downregulated in ccRCC in three gene expression datasets. (G)–(I) The numbers of PAR-CLIP reads for the 88, 82 and 51 overlapped targets of hsa-mir-376a-1_49g, hsa-mir-376c_48g, and hsa-mir-29c_59u in Part (D) to (F), respectively. The numbers of PAR-CLIP reads for the 88, 82 and 51 overlapped target genes in Fig. [169]5D–F were shown in Fig. [170]5G–I, respectively. Then, among these target genes of hsa-mir-376a-1_49g, hsa-mir-376c_48g, and hsa-mir-29c_59u, genes with more than 10 PAR-CLIP reads with T-to-C variations were retained for further analysis (Supplementary Tables [171]S4.5, [172]S5.5, [173]S6.5, respectively). Among the retained target genes of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u, BCAT1, MTHFD2 and RASSF8 were further analyzed after examining their functional relevance in ccRCC by reading their literature, respectively. Figure [174]6A–C showed the distribution of PAR-CLIP reads for BCAT1, MTHFD2 and RASSF8, respectively. These results indicated that PAR-CLIP reads were significantly accumulated at the complementary sites of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u, respectively. The identified miRNA complementary sites and their P [MATH: s :MATH] values on BCAT1, MTHFD2 and RASSF8 were shown in Fig. [175]6D–F, respectively. The complementary sites of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u were located in the 3 [MATH: ′ :MATH] UTR, CDS and 3 [MATH: ′ :MATH] UTR of BCAT1, MTHFD2 and RASSF8, respectively (as shown in Fig. [176]6G to [177]6I, respectively). BCAT1 and MTHFD2 were upregulated in ccRCC, while RASSF8 was downregulated in ccRCC (Supplementary Fig. [178]S5 and Supplementary Tables [179]S10.1–[180]S10.3). Figure [181]6J–L showed the details of complementary sites of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u on BCAT1 ([182]NM_001178093.1), MTHFD2 ([183]NM_006636.3) and RASSF8 ([184]NM_001164746.1), respectively, as well as the PAR-CLIP reads at these sites. Figure 6. [185]Figure 6 [186]Open in a new tab The details of the selected target genes of edited hsa-mir-376a-3p, hsa-mir-376c-3p and hsa-mir-29c-3p. (A–C) The distributions of PAR-CLIP reads on BCAT1 ([187]NM_001178093.1), MTHFD2 ([188]NM_006636.3) and RASSF8 ([189]NM_001164746.1), respectively. The peaks pointed by blue arrows were the complementary sites of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u. (D–F) The identified miRNA sites and their P [MATH: s :MATH] values on BCAT1 ([190]NM_001178093.1), MTHFD2 ([191]NM_006636.3) and RASSF8 ([192]NM_001164746.1), respectively. (G–I) The loci of the complementary sites of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u on BCAT1 ([193]NM_001178093.1), MTHFD2 ([194]NM_006636.3) and RASSF8 ([195]NM_001164746.1), respectively. (J–L) The details of complementary sites of hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u on BCAT1 ([196]NM_001178093.1), MTHFD2 ([197]NM_006636.3) and RASSF8 ([198]NM_001164746.1), respectively, as well as PAR-CLIP reads on these sites. In Part (J) to (L), blue nucleotides indicate T-to-C on mRNA and PAR-CLIP sequencing reads. The orange nucleotides in the seed regions of hsa-mir-376a-1_49g, hsa-mir-376c_48g, and hsa-mir-29c_59u indicate the Guanine and Uracil introduced by the RNA editing. And the three numbers on the right indicate the number of raw sequencing read, the length of the read, and the weight of the read at this locus, respectively. In Parts (A), (D) and (G), the locus pointed by blue arrow is the complementary site of hsa-mir-376a-1_49g on BCAT1 ([199]NM_001178093.1) in Part (J), and the same relations apply for Parts (B), (E), (H) and (K), and Parts (C), (F), (I) and (L), respectively. Functional analysis of A-to-I and C-to-U edited miRNAs We performed Gene Ontology (GO) and KEGG pathway enrichment analysis to investigate the functions of the new target genes that were commonly deregulated in the three gene expression datasets selected for hsa-mir-376a-1_49g, hsa-mir-376c_48g, and hsa-mir-29c_59u, respectively (Supplementary Fig. [200]S6, Supplementary Tables [201]S11.1–[202]S11.3 and [203]S12.1–[204]S12.3). The target genes of hsa-mir-376a-1_49g were enriched in many cancer-related pathways, including “Transcriptional misregulation in cancer”, “Small cell lung cancer”, “Prostate cancer”, “Pathways in cancer”, “p53 signaling pathway”, “Non-small cell lung cancer”, “Melanoma”, “Glioma”, “Endometrial cance”, “Colorectal cancer”, “Chronic myeloid leukemia” and “Breast cancer” (Supplementary Fig. [205]S6A). The target genes of hsa-mir-376c_48g were mainly enriched in the “Pathways in cancer” and the “Rap1 signaling pathway” (Supplementary Fig. [206]S6C). Additionally, “Metabolic pathways” was the most enriched pathway of the 51 target genes of hsa-mir-29c_59u (Supplementary Fig. [207]S6E). These results suggest that hsa-mir-376a-1_49g and hsa-mir-376c_48g are directly related to ccRCC by targeting genes in many cancer-related pathways, while hsa-mir-29c_59u seems to be involved in the deregulation of genes in metabolic pathways. Discussion We identified 1025 M/E sites with significant editing levels by analyzing 176 sRNA-seq profiles of ccRCC patients (n = 154) and healthy controls (n = 22). The M/E sites at the 3 [MATH: ′ :MATH] end accounted for the majority, which was consistent with previous work^[208]37,[209]40,[210]43,[211]56–[212]59. We identified 122 editing events with significantly different editing levels in ccRCC samples compared to healthy control samples. The editing levels of 63.1% and 36.9% of these 122 miRNA editing events were increased and decreased in ccRCC patients, respectively, which may be related to the pathogenesis of ccRCC. Subsequently, we focused on A-to-I and C-to-U miRNA editing events that occurred in the seed regions and were dysregulated in ccRCC. The editing levels of hsa-mir-376a-1_49g and hsa-mir-376c_48g were significantly lower in ccRCC than those in normal samples. In addition, their levels also showed the same significant downregulation in ccRCC samples compared to normal controls. The editing level of hsa-mir-29c_59_C_u and level of hsa-mir-29c_59u showed significant increases in ccRCC compared to normal controls. We compared the targets of original hsa-miR-376a-3p, hsa-miR-376c-3p and hsa-miR-29c-3p with the targets of edited miRNAs, hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u. We found that hsa-mir-376a-1_49g, hsa-mir-376c_48g and hsa-mir-29c_59u could target 484, 568 and 480 additional novel targets, respectively. C-to-U editing of miRNA was not widely reported. Interestingly, we detected 30 C-to-U editing events, much more than A-to-I editing events identified, in ccRCC. The more prevalence of C-to-U miRNA editing in ccRCC may be related to the high levels of APOBEC gene family in ccRCC (Supplementary Fig. [213]S4). For examples, APOBEC3C, APOBEC3F and APOBEC3G were highly expressed in ccRCC samples compared to normal controls (as shown in Supplementary Figs. [214]S4F, H and I, respectively). The fifth position of hsa-miR-100-5p was specifically edited in Treg^[215]60. C-to-U edited hsa-miR-100-5p repressed a new target SMAD2, which reduced the fraction of Treg in peripheral blood mononuclear cells^[216]60. The C-to-U editing of hsa-mir-29c_59_C_u in ccRCC might represent another case of functional C-to-U editing sites in miRNAs, because its editing levels significantly increased in ccRCC and the C-to-U edited hsa-mir-29c-3p potentially repressed RASSF8 in ccRCC. As 11 of these 30 C-to-U editing sites were in seed regions (Fig. [217]3 and Supplementary Table [218]S2.5), more of these C-to-U editing might be functional in tissue- or time-specific manners. In our study, BCAT1 and MTHFD2 were upregulated in ccRCC and identified as target genes of hsa-mir-376a-1_49g and hsa-mir-376c_48g in ccRCC, respectively. BCAT1 is overexpressed in many cancers and has been proposed as a marker for predicting cancer prognosis^[219]61. For example, overexpression of BCAT1 promotes tumor growth in renal clear cell carcinoma^[220]62, lung cancer^[221]63, glioma^[222]64, ovarian cancer^[223]65, breast cancer^[224]66, hepatocellular carcinoma^[225]67,[226]68, myeloid leukemia^[227]69, gastric cancer^[228]70,[229]71, endometrial cancer^[230]72, non-small cell lung cancer^[231]73 and pan-cancer^[232]74. Moreover, BCAT1 and its metabolites participate in the metabolism of cancer cells through different mechanisms. For example, the PI3K/AKT/mTOR signaling pathway is activated by BCAT1 to promote the proliferation and angiogenesis of gastric cancer cells in vitro^[233]70. At the same time, in the study of lung adenocarcinoma, BCAT1 promotes lung adenocarcinoma progression by enhancing mitochondrial function and NF- [MATH: κ :MATH] B pathway^[234]75. Recent reports have shown that MTHFD2 is also highly expressed in many cancers, including renal clear cell carcinoma^[235]76,[236]77, pancreatic cancer,^[237]78 breast cancer^[238]79, myeloid leukemia^[239]80, non-small cell lung cancer^[240]81,[241]82, hepatocellular carcinoma^[242]83, colorectal cancer^[243]84, breast cancer^[244]85,[245]86, ovarian cancer^[246]87 , bladder cancer^[247]88, glioblastoma^[248]89, nasopharyngeal carcinoma^[249]90, oral squamous cell carcinoma^[250]91, and pan-cancer^[251]92. In addition, studies have reported the role of MTHFD2 in tumor immune evasion. MTHFD2 is induced by IFN- [MATH: γ :MATH] and promotes basal and IFN- [MATH: γ :MATH] -induced PD-L1 expression. Mechanistically, MTHFD2 drives folate circulation to maintain intracellular UDP-GlcNAc and cMYC O-GlcNAcylation, which promotes PD-L1 transcription^[252]93. Huang et al. found that MTHFD2 promoted breast cancer cell proliferation through AKT signaling pathway, and high level of MTHFD2 reduced the survival rate of breast cancer patients^[253]79. In summary, the upregulation of BCAT1 and MTHFD2 in ccRCC may be due to the lower levels of hsa-mir-376a_49g and hsa-mir-376c_48g in ccRCC, respectively, which may contribute to the initiation and/or progression of ccRCC. Reduced A-to-I editing level of miRNAs were also noticed previously^[254]38,[255]39. The editing levels of hsa-mir-376a-1_49_A_g and/or hsa-mir-376c_48_A_g were also reduced in many other cancers, such as glioblastoma^[256]94 and lung cancer^[257]95. The increased expression of BCAT1 and MFTHD2 might be due to the reduced editing levels of hsa-mir-376a-1_49_A_g and hsa-mir-376c_48_A_g in many cancers. Furthermore, it has been shown that the higher the editing level of hsa-mir-376c_48_A_g site in ccRCC, the better the prognosis and the longer the survival time of patients^[258]38. Therefore, reduced editing levels of hsa-mir-376a-1_49_A_g and hsa-mir-376c_48_A_g might represent a general mechanism in the initiation and/or progression of different cancers. RASSF8 has been reported to be a tumor suppressor gene with lower expression in gastric cancer^[259]96,[260]97, colorectal cancer^[261]98,[262]99, cervical cancer^[263]100, ovarian cancer^[264]101, melanoma^[265]102, esophageal squamous cell carcinoma^[266]103, lung cancer^[267]104,[268]105 and osteosarcoma^[269]106. In addition, RASSF8 is a potential therapeutic target for the prevention of many cancers^[270]97,[271]100,[272]103,[273]105. In this study, RASSF8, the target gene of hsa-mir-29c_59u, was downregulated in ccRCC. Therefore, RASSF8 may also play similar tumor suppressor function in ccRCC, and the higher level of hsa-mir-29c_59u in ccRCC may lead to the decreased expression of RASSF8 in ccRCC, which consequently promotes the initiation and/or progression of ccRCC. As shown in Fig. [274]5G–I, some genes with many PAR-CLIP reads were identified as targets of hsa-mir-376a-1_49g and hsa-mir-376c_48g, and hsa-mir-29c_59u, respectively. Although there currently are limited evidences about their functions in ccRCC, our results suggest that these genes represent promising directions for revealing the mechanisms of ccRCC in the future. In summary, we presented a comprehensive view of miRNA editing and/or modification events in ccRCC, which promoted our understanding of miRNA editing and/or modification events in ccRCC. More researches are needed in the future to clarify the functional roles of different miRNA editing patterns in the pathogenesis of ccRCC. Our results also provided new clues for the clinical treatment of ccRCC, such as to promote the editing levels of hsa-mir-376a_49_A_g and hsa-mir-376c_48_A_g and to repress the editing level of hsa-mir-29c_59_C_u. Methods Small RNA-Seq data used We collected 176 small RNA-Seq (sRNA-seq) profiles of ccRCC patients and healthy controls, including 154 samples of cancer tissue and 22 normal kidney tissue samples. These 176 sRNA-seq profiles were obtained from previous studies^[275]107,[276]108, and their accession numbers were shown in Supplementary Table [277]S1.1. The qualities of these sRNA-seq sequence data were evaluated with FastQC (v0.11.9)^[278]109. Genome and annotation of miRNAs used The human genome sequence used was GRCh38, and downloaded from the UCSC Genome Browser^[279]110. The pri-miRNA sequences in GFF3 format and mature miRNA annotation files were derived from miRBase (v21)^[280]111. Identification of mutation and editing sites in miRNAs Totally, 176 sRNA-seq profiles were analyzed using the MiRME pipeline with default settings^[281]43. Briefly, the raw reads were checked to retain the qualified reads with the sequencing scores greater than 30 for the first 25 nucleotides. Then, reads of at least 18 nt were retained after removal of the 3 [MATH: ′ :MATH] adapters. BLASTN^[282]112 was used to compare the retained reads with the pre-miRNAs and the reads mapped to pre-miRNAs were retrieved. These reads, which were mapped to pre-miRNAs, were then aligned with the genome using Bowtie (v1.0.0)^[283]113. Then, the cross-mapping correction algorithm^[284]28 was used to examine the alignment of reads to the genome to calculate the weights or percentages of reads at different genomic loci. Results from different samples were then combined using separate programs in the MiRME package^[285]43,[286]55. According to the location of M/E sites in miRNAs and mutations status of dbSNP, the M/E sites were divided into nine different editing types, i.e., A-to-I, C-to-U, 3 [MATH: ′ :MATH] -A, 3 [MATH: ′ :MATH] -U, 3 [MATH: ′ :MATH] -Other, 5 [MATH: ′ :MATH] -editing, Other, SNP and Pseudo^[287]43. The identified M/E site was named with the pre-miRNA name, M/E site position in pre-miRNA, original nucleotide in upper case, the edited/mutated nucleotide in lower case. For example, hsa-mir-376c_48_A_g means an A-to-I editing at the 48th nucleotide of hsa-mir-376c. And edited miRNA was named by the pre-miRNA name, the M/E site position in pre-miRNA, and the edited/mutated nucleotide in lower case. For example, hsa-mir-376c_48g is the A-to-I edited miR-376c-3p. The criteria for defining M/E sites with significant editing levels were as follows: (i) the relative levels of editing were at least 5%; (ii) editing events were supported by at least 10 reads; (iii) the threshold for sequencing reads score was 30; (iv) multiple test corrected p-values (with the Benjamini and Hochberg method^[288]114) were smaller than 0.05. In order to remove M/E sites due to random sequencing errors, 1025 M/E sites that had significant editing levels in at least 10% of the 176 samples (18 samples) used in this study were retained for further analysis. The identified M/E sites were compared with known human miRNA editing sites, including the DARNED database^[289]115, the RADAR database^[290]116 and relevant literature^[291]27,[292]29,[293]34,[294]43,[295]95,[296]117,[297]118. Finally, A-to-I, C-to-U, and Other predicted M/E sites were manually checked. Identification of SNPs sites The identified M/E sites were compared with the dbSNP (v151) database^[298]119 and previously reported SNPs in miRNAs^[299]44. An M/E site was considered as SNP if it satisfied the following criteria: (i) the genomic position of M/E site and SNP was identical, (ii) the original and edited nucleotides had the same nucleotides as the SNP’s allele, (iii) the M/E site had editing level of 100% in at least one of the 176 samples selected, and (iv) the M/E site occurred in the center of the miRNA. Identification of M/E sites with significantly different editing levels in ccRCC The Mann-Whitney U-tests were used to analyze the differences between the editing levels of 1025 miRNA M/E sites in ccRCC (n = 154) and normal samples (n = 22). The p-values were corrected with the Benjamini-Hochberg correction method^[300]114. The M/E sites with corrected p-values smaller than 0.05 were defined as having significantly different editing levels in ccRCC. DESeq2^[301]120 was used to analyze different expression levels (TPTM) of edited miRNAs between ccRCC and normal control samples. The miRNAs with corrected p-values smaller than 0.05 were defined as having different expression levels in ccRCC and normal control samples. Identification of targets for original and edited miRNAs Among the editing sites with significant differences in editing levels between ccRCC tumors and normal tissues, a total of five editing events occurred in seed regions of mature miRNAs. We chose two A-to-I editing sites (hsa-mir-376a-1_49_A_g and hsa-mir-376c_48_A_g) and one C-to-U editing site (hsa-mir-29c_59_C_u), because their editing levels and the levels of their corresponding edited miRNAs had significant differences in ccRCC tumor tissues compared to normal controls. Then, we predicted the targets of original and edited miRNAs using the MiCPAR algorithm^[302]55. As shown in Supplementary Table [303]S1.2, 11 PAR-CLIP sequences were used in the identification of miRNA targets and were downloaded from the NCBI SRA database. Seven PAR-CLIP profiles of HEK293 cell lines stably expressing FLAG/HA-tagged AGO1, AGO2, AGO3, and AGO4 proteins were reported by^[304]54. Another study included four PAR-CLIP profiles derived from HEK293 cell lines stably expressing HIS/FLAG/HA-tagged AGO1 and AGO2 proteins^[305]121. In order to obtain reads of acceptable quality, the raw reads of these 11 PAR-CLIP profiles were filtered to keep reads with sequencing scores of at least 30 for their first 25 nucleotides from 5 [MATH: ′ :MATH] ends. The qualified reads were combined and the targets of miRNA were identified by the MiCPAR algorithm. The targets with at least one PAR-CLIP read with T-to-C variation were retained for further analysis. Analyzing expression of targets of edited miRNAs In order to understand the expression of targets of edited miRNAs in ccRCC tumor tissues, we identified genes that were significantly differentially expressed in ccRCC samples compared to normal samples. We collected three batches of gene expression profiles of ccRCC and controls (Supplementary Table [306]S1.3), of which two batches were obtained from^[307]107,[308]122. Another set of gene expression profiles was obtained from TCGA (https://portal.gdc.cancer.gov/). The edgeR (v3.34.1) package^[309]123 was used for differential analysis of target genes. The glmFit and glmLRT functions in edgeR were used to build generalized linear models and to perform likelihood ratio tests, respectively. Genes with corrected p-values smaller than 0.05 were defined as having significantly different expression levels in ccRCC. Functional analysis of new target genes of edited miRNAs We kept targets with at least one PAR-CLIP read with T-to-C variation. Then, the targets of original and edited miRNAs were compared to obtain the new targets of edited miRNAs. For edited miRNAs with increased editing levels in ccRCC, we compared their targets with genes that were downregulated in ccRCC because miRNAs normally repressed their target genes, and vice versa. We also analyzed the numbers of PAR-CLIP reads for the targets of edited miRNAs. We manually examined the functional relevance of some selected targets by reading the literature of the genes in the NCBI Gene database. The targets with more PAR-CLIP reads and functional relevance in ccRCC were preferred in our analysis. KOBAS3.0^[310]124 was used to perform enrichment analysis of GO terms and KEGG pathways for new target genes of edited miRNAs. Significantly enriched GO terms and KEGG pathways were filtered with multiple test corrected p-values smaller than 0.05. Supplementary Information [311]Supplementary Information 1.^ (2.9MB, docx) [312]Supplementary Information 2.^ (32.2KB, xlsx) [313]Supplementary Information 3.^ (20MB, xlsx) [314]Supplementary Information 4.^ (146.3KB, xlsx) [315]Supplementary Information 5.^ (1.2MB, xlsx) [316]Supplementary Information 6.^ (1.3MB, xlsx) [317]Supplementary Information 7.^ (1.2MB, xlsx) [318]Supplementary Information 8.^ (36.5KB, xlsx) [319]Supplementary Information 9.^ (65.7KB, xlsx) [320]Supplementary Information 10.^ (41.7KB, xlsx) [321]Supplementary Information 11.^ (54.8KB, xlsx) [322]Supplementary Information 12.^ (225.6KB, xlsx) Author contributions Y.L.: data curation, formal analysis, investigation, validation, visualisation, writing - original draft; S.G.: data curation, formal analysis, investigation, validation, visualisation; W.X.: data curation, formal analysis, investigation, validation; H.Y.: data curation, formal analysis; W.L.: investigation, validation; N.Z.: formal analysis, visualisation; J.Y.: formal analysis; G.Z.: validation; C.M.: validation; Y.Z.: conceptualisation, formal analysis, funding acquisition, methodology, project administration, resources, software, supervision, writing - review & editing. Funding The research was supported in part by a grand (No. 31460295) of the National Natural Science Foundation of China and an Open Research Fund (No. SKLGE-2107) of State Key Laboratory of Genetic Engineering, Fudan University, China to YZ. Data availibility The 176 sRNA-seq sequencing profiles and 11 PAR-CLIP sequencing profiles were obtained from NCBI SRA database, and 3 gene expression profiles were obtained from NCBI GEO and TCGA database with accession numbers shown in Table S1. Competing interests The authors declare no competing interests. Footnotes Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Supplementary Information The online version contains supplementary material available at 10.1038/s41598-023-42302-y. References