Research Paper Volume 11, Issue 15 pp 5666—5688

Overexpressed pseudogenes, DUXAP8 and DUXAP9, promote growth of renal cell carcinoma and serve as unfavorable prognostic biomarkers

Jing Chen 1, 2, , Weiyang Lou 3, , Bisha Ding 3, , Xian Wang 1, ,

  • 1 Department of Medical Oncology, Sir Run Run Shaw Hospital, Medical School of Zhejiang University, Zhejiang Province, Hangzhou 313100, China
  • 2 First Affiliated Hospital of Jiaxing University, Zhejiang Province, Jiaxing 314000, China
  • 3 Department of Surgery, First Affiliated Hospital, College of Medicine, Zhejiang University, Hangzhou 313100, China

received: May 11, 2019 ; accepted: July 31, 2019 ; published: August 13, 2019 ;

https://doi.org/10.18632/aging.102152
How to Cite

Copyright © 2019 Chen et al. This is an open-access article distributed under the terms of the Creative Commons Attribution (CC BY) 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Growing studies have reported that pseudogenes play key roles in multiple human cancers. However, expression and roles of pseudogenes in renal cell carcinoma remains absent.

Results: 31 upregulated and 16 downregulated pseudogenes were screened. Higher expression of DUXAP8 and DUXAP9 indicated poorer prognosis of kidney cancer. 33 and 5 miRNAs were predicted to potentially binding to DUXAP8 and DUXAP9, respectively. miR-29c-3p was identified as the most potential binding miRNAs of DUXAP8 and DUXAP9 based on expression, survival and correlation analyses. 254 target genes of miR-29c-3p were forecast. 47 hub genes with node degree >= 10 were identified. Subsequent analysis for the top 10 hub genes demonstrated that COL1A1 and COL1A2 may be two functional targets of DUXAP8 and DUXAP9. Expression of DUXAP8, DUXAP9, COL1A1 and COL1A2 were significantly increased in cancer samples compared to normal controls while miR-29c-3p expression was decreased. Luciferase reporter assay revealed that miR-29c-3p could directly bind to DUXAP8, DUXAP9, COL1A1 and COL1A2. Functional experiments showed that DUXAP8 and DUXAP9 enhanced but miR-29c-3p weakened growth of renal cell carcinoma.

Conclusions: In conclusion, upregulated DUXAP8 and DUXAP9 promote growth of renal cell carcinoma and serve as two promising prognostic biomarkers.

Methods: Dysregulated pseudogenes were obtained by dreamBase and GEPIA. The binding miRNAs of pseudogene and targets of miRNA were predicted using starBase and miRNet. Kaplan-Meier plotter was utilized to perform survival analysis, and Enrichr database was introduced to conduct functional enrichment analysis. Hub genes were identified through STRING and Cytoscape. qRT-PCR, luciferase reporter assay, cell counting assay and colony formation assay were performed to validate in silico analytic results.

Introduction

Kidney malignancy is the twelfth most common cancer all over the whole world [1]. Renal cell carcinoma is the most common type of primary kidney malignancy and accounts for 90%–95% of all kidney cancer cases [2]. Among all subtypes of renal cell carcinoma (RCC), clear cell renal cell carcinoma (ccRCC) is the most common one, approximately accounting for 75% of all RCC cases [3]. In addition to ccRCC, there are also other subtypes, such as chromophobe renal cell carcinoma (chRCC) and papillary renal cell carcinoma (pRCC). Over the past few decades, great improvements have been achieved in the diagnosis, therapy and prognosis of RCC. However, the patients’ prognosis remains unsatisfactory, partially owing to the high rate of cancer recurrence after undergoing surgical resection [4]. Thus, it is necessary to unlock the molecular mechanisms of RCC onset and progression and develop novel effective diagnostic, therapeutic and prognostic biomarkers for RCC.

Pseudogenes, genomic DNA sequences similar to normal genes, are regarded as defunct relatives of functional genes [5, 6]. Pseudogenes are generally divided into two main types: processed (or retrotransposed) pseudogene and non-processed (or duplicated) pseudogenes [79]. After the first discovery of pseudogene structure in 5S DNA of Xenopus laevis by C. Jacq et al. in 1977 [10], increasing studies have demonstrated that dysregulation of pseudogenes plays key roles in development of human disorders and these deregulated pseudogenes may act as promising therapeutic targets for diseases [11, 12]. Recently, pseudogenes have emerged as critical regulators in the pathogenesis of cancer, including lung cancer [13], pancreatic cancer [14], hepatocellular carcinoma [9], brain glioma [15], colon cancer [16], gastric cancer [17, 18], breast cancer [19, 20] as well as kidney cancer [21].

Pseudogenes could serve as competing endogenous RNAs (ceRNAs) of their parental genes and other unrelated genes through competitively binding to shared miRNAs, thereby exerting a variety of biological functions [22]. For example, pseudogene OCT4-pg4 regulates OCT4 expression by competing for miR-145 in hepatocellular carcinoma [23]; pseudogene TUSC2P promotes TUSC2 function by binding multiple microRNAs [24]; a FTH1 gene-pseudogene-microRNA network modulates development of prostate cancer [25]. However, to the best of our knowledge, pseudogene-miRNA-gene network in renal cell carcinoma has not yet been fully elucidated. Exploration of functional pseudogenes associated with renal cell carcinoma and their corresponding ceRNA mechanism may provide novel insights for diagnosis, therapy and prognosis of renal cell carcinoma in the future.

In this study, we first screened differentially expressed pseudogenes in RCC by analyzing dreamBase and GEPIA. Next, we further found those prognosis- related pseudogenes among all the differentially expressed pseudogenes. Using starBase and miRNet databases, a pseudogene-miRNA-gene network was subsequently constructed to explore the underlying molecular mechanisms of key pseudogenes in RCC. Finally, a series of experiments, including qRT-PCR, luciferase reporter assay, cell counting assay and colony formation assay, were conducted to validate analytic results.

Results

Screening dysregulated pseudogenes in human ccRCC

As mentioned above, ccRCC is the most common type of RCC. Therefore, differentially expressed pseudogenes in ccRCC was first obtained to explore the functional pseudogenes and excavate their potential mechanisms of action in RCC using human pseudogene database, namely dreamBase. Based on the cut-off criteria, a total of 399 dysregulated pseudogenes, containing 335 upregulated and 64 downregulated pseudogenes, were finally identified in human ccRCC (Supplementary Table 1). In order to confirm the preliminary results, we determined expression levels of these dysregulated pseudogenes in human ccRCC using another online database GEPIA. 31 out of 335 upregulated pseudogenes and 16 out of 64 downregulated pseudogenes were found to be consistent with the previous analytic results from dreamBase database as listed in Supplementary Table 2. The 47 differentially expressed pseudogenes were defined as potential pseudogenes and were selected for further analysis.

Identification of DUXAP8 and DUXAP9 as two key pseudogenes in RCC

Next, we determined the prognostic values of 47 potential pseudogenes in human ccRCC based on TCGA data using GEPIA database. Two prognostic indexes, consisting of overall survival (OS) and disease free survival (RFS), were included in this part. Among the 47 potential pseudogenes, only high expression of three upregulated pseudogenes (AC007326.9, DUXAP8 and DUXAP9) and four downregulated pseudogenes (NUDT4P2, RP11-255H23.2, AF186192.5 and SLC2A3P1) indicated poorer and better OS and RFS in patients with ccRCC, respectively (Supplementary Table 3). The expression levels of the 7 pseudogenes between cancer and normal samples were displayed in Figure 1A1-A7. Moreover, expression analysis among various major stages of the 7 pseudogenes demonstrated that all of them possessed statistical significance (Figure 1B1B7). Figure 1C1C7 and Figure 1D1D7 revealed OS and RFS of the 7 pseudogenes, respectively. To discover the most potential functional pseudogenes in RCC, we performed survival analysis (OS and RFS) of the 7 pseudogenes in chRCC and pRCC. As listed in Table 1 and Figure 2, only chRCC and pRCC patients with higher expression of DUXAP8 and DUXAP9 predicted poorer prognosis. All these findings suggest that DUXAP8 and DUXAP9 may serve as two potential prognostic biomarkers in human RCC.

Expression and prognostic values of 7 potential pseudogenes in clear cell renal cell carcinoma (ccRCC) determined by GEPIA database. (A1) Expression of AC007326.9 in ccRCC compared with normal controls; (A2) expression of DUXAP8 in ccRCC compared with normal controls; (A3) expression of DUXAP9 in ccRCC compared with normal controls; (A4) expression of NUDT4P2 in ccRCC compared with normal controls; (A5) expression of RP11-255H23.2 in ccRCC compared with normal controls; (A6) expression of AF186192.5 in ccRCC compared with normal controls; (A7) expression of SLC2A3P1 in ccRCC compared with normal controls; (B1) expression of AC007326.9 among major stages in ccRCC; (B2) expression of DUXAP8 among major stages in ccRCC; (B3) expression of DUXAP9 among major stages in ccRCC; (B4) expression of NUDT4P2 among major stages in ccRCC; (B5) expression of RP11-255H23.2 among major stages in ccRCC; (B6) expression of AF186192.5 among major stages in ccRCC; (B7) expression of SLC2A3P1 among major stages in ccRCC; (C1) prognostic role (overall survival) of AC007326.9 in ccRCC; (C2) prognostic role (overall survival) of DUXAP8 in ccRCC; (C3) prognostic role (overall survival) of DUXAP9 in ccRCC; (C4) prognostic role (overall survival) of NUDT4P2 in ccRCC; (C5) prognostic role (overall survival) of RP11-255H23.2 in ccRCC; (C6) prognostic role (overall survival) of AF186192.5 in KIRC; (C7) prognostic role (overall survival) of SLC2A3P1 in ccRCC; (D1) prognostic role (disease free survival) of AC007326.9 in ccRCC; (D2) prognostic role (disease free survival) of DUXAP8 in ccRCC; (D3) prognostic role (disease free survival) of DUXAP9 in ccRCC; (D4) prognostic role (disease free survival) of NUDT4P2 in ccRCC; (D5) prognostic role (disease free survival) of RP11-255H23.2 in ccRCC; (D6) prognostic role (disease free survival) of AF186192.5 in ccRCC; (D7) prognostic role (disease free survival) of SLC2A3P1 in ccRCC. “*” represents P-value less than 0.05. Three horizontal lines in the box plot represent minimum, median and maximum, respectively.

Figure 1. Expression and prognostic values of 7 potential pseudogenes in clear cell renal cell carcinoma (ccRCC) determined by GEPIA database. (A1) Expression of AC007326.9 in ccRCC compared with normal controls; (A2) expression of DUXAP8 in ccRCC compared with normal controls; (A3) expression of DUXAP9 in ccRCC compared with normal controls; (A4) expression of NUDT4P2 in ccRCC compared with normal controls; (A5) expression of RP11-255H23.2 in ccRCC compared with normal controls; (A6) expression of AF186192.5 in ccRCC compared with normal controls; (A7) expression of SLC2A3P1 in ccRCC compared with normal controls; (B1) expression of AC007326.9 among major stages in ccRCC; (B2) expression of DUXAP8 among major stages in ccRCC; (B3) expression of DUXAP9 among major stages in ccRCC; (B4) expression of NUDT4P2 among major stages in ccRCC; (B5) expression of RP11-255H23.2 among major stages in ccRCC; (B6) expression of AF186192.5 among major stages in ccRCC; (B7) expression of SLC2A3P1 among major stages in ccRCC; (C1) prognostic role (overall survival) of AC007326.9 in ccRCC; (C2) prognostic role (overall survival) of DUXAP8 in ccRCC; (C3) prognostic role (overall survival) of DUXAP9 in ccRCC; (C4) prognostic role (overall survival) of NUDT4P2 in ccRCC; (C5) prognostic role (overall survival) of RP11-255H23.2 in ccRCC; (C6) prognostic role (overall survival) of AF186192.5 in KIRC; (C7) prognostic role (overall survival) of SLC2A3P1 in ccRCC; (D1) prognostic role (disease free survival) of AC007326.9 in ccRCC; (D2) prognostic role (disease free survival) of DUXAP8 in ccRCC; (D3) prognostic role (disease free survival) of DUXAP9 in ccRCC; (D4) prognostic role (disease free survival) of NUDT4P2 in ccRCC; (D5) prognostic role (disease free survival) of RP11-255H23.2 in ccRCC; (D6) prognostic role (disease free survival) of AF186192.5 in ccRCC; (D7) prognostic role (disease free survival) of SLC2A3P1 in ccRCC. “*” represents P-value less than 0.05. Three horizontal lines in the box plot represent minimum, median and maximum, respectively.

Table 1. Prognostic roles of 7 potential pseudogenes in chRCC and pRCC determined by GEPIA database.

Pseudogene namechRCC apRCC bEffecte
OScRFSdOScRFSd
AC007326.9NoNoNoNoNo
DUXAP8PoorPoorPoorNoPoor
DUXAP9PoorPoorPoorPoorPoor
NUDT4P2NoNoNoNoNo
RP11-255H23.2NoNoNoNoNo
AF186192.5NoNoGoodGoodNo
SLC2A3P1NoNoNoNoNo
achRCC, chromophobe renal cell carcinoma;
bpRCC, papillary renal cell carcinoma;
cOS, Overall Survival;
dRFS, Disease Free Survival;
eEffect, total effect of individual pseudogene in chRCC and pRCC.
Survival analysis of DUXAP8 and DUXAP9 in chRCC (chromophobe renal cell carcinoma) and pRCC (papillary renal cell carcinoma). (A) Prognostic value (overall survival) of DUXAP8 in chRCC; (B) prognostic value (disease free survival) of DUXAP8 in chRCC; (C) prognostic value (overall survival) of DUXAP8 in pRCC; (D) prognostic value (disease free survival) of DUXAP8 in pRCC; (E) prognostic value (overall survival) of DUXAP9 in chRCC; (F) prognostic value (disease free survival) of DUXAP9 in chRCC; (G) prognostic value (overall survival) of DUXAP9 in pRCC; (H) prognostic value (disease free survival) of DUXAP9 in pRCCs.

Figure 2. Survival analysis of DUXAP8 and DUXAP9 in chRCC (chromophobe renal cell carcinoma) and pRCC (papillary renal cell carcinoma). (A) Prognostic value (overall survival) of DUXAP8 in chRCC; (B) prognostic value (disease free survival) of DUXAP8 in chRCC; (C) prognostic value (overall survival) of DUXAP8 in pRCC; (D) prognostic value (disease free survival) of DUXAP8 in pRCC; (E) prognostic value (overall survival) of DUXAP9 in chRCC; (F) prognostic value (disease free survival) of DUXAP9 in chRCC; (G) prognostic value (overall survival) of DUXAP9 in pRCC; (H) prognostic value (disease free survival) of DUXAP9 in pRCCs.

Identification of miR-29c-3p as a potential binding miRNA of DUXAP8 and DUXAP9

Increasing evidences have suggested that ceRNA hypothesis is an important regulatory mechanism for ncRNA, including lncRNA, circRNA as well as pseudogenes [22, 26]. Therefore, to uncover underlying mechanisms and functions of DUXAP8 and DUXAP9, we first predicted potential binding miRNAs by starBase database. Finally, 33 and 5 miRNAs were predicted to potentially bind to DUXAP8 and DUXAP9, respectively. For better visualization, DUXAP8/DUXAP9-miRNA network was established as presented in Figure 3. According to the ceRNA mechanism, there should be a negative correlation between miRNA expression and DRXAP8/9 expression. Thus, we determined the expression correlation using starBase database. As shown in Table 2, among these pseudogene-miRNA pairs, only four interactions (DUXAP8-miR-29c-3p, DUXAP8-miR-92b-3p, DUXAP8-miR-500a-3p and DUXAP9-miR-29c-3p) demonstrated significantly negative relationship in ccRCC. In the next step, we determined the expression and prognostic values of the three miRNAs (miR-29c-3p, miR-92b-3p and miR-500a-3p) in ccRCC. The results from expression analysis revealed that miR-29c-3p was significantly downregulated in ccRCC samples when compared with normal controls (Figure 4A). Figure 4B told us a marked upregulation of miR-92b-3p expression in cancer samples. We also found that miR-500a-3p expression in cancer tissues was significantly lower than that in normal tissues (Figure 4C). Survival analysis for the three miRNAs demonstrated that ccRCC patients with higher expression of miR-29c-3p and miR-92b-3p had better and worse prognosis as presented in Figure 4D and Figure 4E, respectively. Figure 4F showed that there was no significant survival difference between miR-500a-3p low expression group and miR-500a-3p high expression group. These findings demonstrated that miR-29c-3p expression was conversely correlated with DUXAP8 and DUXAP9 expression, and miR-29c-3p was significantly downregulated in ccRCC and the expression downregulation indicates a poor prognosis of ccRCC patients. As shown in Supplementary Figure 1A and Supplementary Figure 1B, miR-29c-3p was also significantly downregulated in chRCC and pRCC. Furthermore, we also found that pRCC patients with high expression of miR-29c-3p had a favorable prognosis (Supplementary Figure 1C). Taken together, miR-29c-3p may be a potential binding miRNA of DUXAP8 and DUXAP9 in renal cell carcinoma.

Establishment of the potential DUXAP8/DUXAP9-miRNA regulatory network in renal cell carcinoma.

Figure 3. Establishment of the potential DUXAP8/DUXAP9-miRNA regulatory network in renal cell carcinoma.

Table 2. Correlation of potential pseudogene and miRNA in clear cell renal cell carcinoma (ccRCC) determined by starBase version 3.

PseudogenemiRNARP-value
DUXAP8hsa-miR-29b-3p0.0874.69E-02
DUXAP8hsa-miR-29c-3p−0.1582.97E-04
DUXAP8hsa-miR-29a-3p0.061.75E-01
DUXAP8hsa-miR-382-5p0.2393.87E-08
DUXAP8hsa-miR-488-3p−0.0855.33E-02
DUXAP8hsa-miR-4735-3p01.00E+00
DUXAP8hsa-miR-18a-5p0.1155.09E-02
DUXAP8hsa-miR-18b-5p−0.0029.71E-01
DUXAP8hsa-miR-5581-3p0.0088.99E-01
DUXAP8hsa-miR-378b−0.0365.37E-01
DUXAP8hsa-miR-378e0.0256.67E-01
DUXAP8hsa-miR-378d0.0931.13E-01
DUXAP8hsa-miR-378f−0.0029.72E-01
DUXAP8hsa-miR-378i0.0069.22E-01
DUXAP8hsa-miR-378a-3p0.1155.09E-02
DUXAP8hsa-miR-422a01.00E+00
DUXAP8hsa-miR-378h−0.0573.34E-01
DUXAP8hsa-miR-378c0.0951.06E-01
DUXAP8hsa-miR-33b-5p0.0523.78E-01
DUXAP8hsa-miR-33a-5p0.0841.52E-01
DUXAP8hsa-miR-490-3p0.1568.06E-03
DUXAP8hsa-miR-32-5p0.1018.68E-02
DUXAP8hsa-miR-363-3p−0.0138.20E-01
DUXAP8hsa-miR-367-3p01.00E+00
DUXAP8hsa-miR-25-3p0.0682.49E-01
DUXAP8hsa-miR-92a-3p0.238.07E-05
DUXAP8hsa-miR-92b-3p0.116.07E-02
DUXAP8hsa-miR-500a-3p0.1451.37E-02
DUXAP8hsa-miR-361109.95E-01
DUXAP8hsa-miR-1914-3p0.018.72E-01
DUXAP8hsa-miR-5194−0.0167.82E-01
DUXAP8hsa-miR-6738-5p01.00E+00
DUXAP8hsa-miR-193a-5p0.1645.10E-03
DUXAP9hsa-miR-488-3p−0.0393.73E-01
DUXAP9hsa-miR-29a-3p0.0661.36E-01
DUXAP9hsa-miR-29b-3p0.0711.06E-01
DUXAP9hsa-miR-29c-3p0.1421.24E-03
DUXAP9hsa-miR-382-5p0.2471.34E-08
Expression and survival analysis of miRNAs in clear cell renal cell carcinoma (ccRCC or KIRC). (A) Box-whisker plot showed expression of hsa-miR-29c-3p in ccRCC compared with normal controls (three horizontal lines in the box plot represent minimum, median and maximum); (B) box-whisker plot showed expression of hsa-miR-92b-3p in ccRCC compared with normal controls (three horizontal lines in the box plot represent minimum, median and maximum); (C) box-whisker plot showed expression of hsa-miR-500a-3p in ccRCC compared with normal controls (three horizontal lines in the box plot represent minimum, median and maximum); (D) prognostic value of hsa-miR-29c-3p in ccRCC; (E) prognostic value of hsa-miR-92b-3p in ccRCC; (F) prognostic value of hsa-miR-500a-3p in ccRCC.

Figure 4. Expression and survival analysis of miRNAs in clear cell renal cell carcinoma (ccRCC or KIRC). (A) Box-whisker plot showed expression of hsa-miR-29c-3p in ccRCC compared with normal controls (three horizontal lines in the box plot represent minimum, median and maximum); (B) box-whisker plot showed expression of hsa-miR-92b-3p in ccRCC compared with normal controls (three horizontal lines in the box plot represent minimum, median and maximum); (C) box-whisker plot showed expression of hsa-miR-500a-3p in ccRCC compared with normal controls (three horizontal lines in the box plot represent minimum, median and maximum); (D) prognostic value of hsa-miR-29c-3p in ccRCC; (E) prognostic value of hsa-miR-92b-3p in ccRCC; (F) prognostic value of hsa-miR-500a-3p in ccRCC.

KEGG pathway enrichment, GO functional annotation, and PPI network analysis for miR-29c-3p target genes

For better understanding the potential functions of pseudogene DUXAP8 and DUXAP9, the target genes of miR-29c-3p were first predicted using a comprehensive tool, miRNet. As shown in Supplementary Table 4, a total of 254 potential target genes were obtained. Next, we mapped them into Enrichr database for KEGG pathway enrichment analysis and GO functional annotation. Three GO items (biological process, BP; cellular component, CC; molecular function, MF) were included in GO functional annotation. KEGG pathway analysis for these target genes demonstrated that they were significantly enriched in several human cancers and cancer-associated pathways, such as small cell lung cancer, renal cell carcinoma, endometrial cancer, colorectal cancer, melanoma, TNF signaling pathway and ECM-receptor interaction (Figure 5A). A variety of significant enriched GO terms were discovered and the top 10 enriched GO terms were shown in Figure 5B5D, including extracellular matrix organization (GO:0030198), negative regulation of cellular senescence (GO:2000773) and positive regulation of substrate adhesion-dependent cell spreading (GO:1900026) in the BP category, endoplasmic reticulum lumen (GO:0005788), chromatin (GO:0000785) and chromatin silencing complex (GO:0005677) in the CC category, and platelet-derived growth factor binding (GO:0048407), transcriptional repressor activity, RNA polymerase II activating transcription factor binding (GO:0098811) and RNA polymerase II regulatory region sequence-specific DNA binding (GO:0000977) in the MF category. It has been widely acknowledged that genes exert their biological functions by interacting with each other instead of acting alone. Thus, we established the PPI network of aforementioned potential target genes of miR-29c-3p using STRING database. The result showed that many gene-gene interactions among these target genes were existed (Data were not shown). To obtain the hub genes, we re-entered the gene-gene interactions into Cytoscape software (version 3.6.0). The node degree of each target gene was calculated by CytoHubba plugin of Cytoscape software. For better visualization, the hub genes with node degree >= 10 were screened (Table 3) and a sub-PPI network were re-constructed as shown in Figure 5E.

Table 3. The hub genes in the PPI network (node degree >= 10).

Gene nameNode degree
GAPDH46
VEGFA46
PTEN38
JUN35
MMP231
SIRT128
ITGB128
CDC4225
COL1A125
COL1A225
MDM224
CASP824
FOS24
PDGFRB23
COL3A122
COL4A122
MYCN21
COL6A221
COL4A221
KLF420
COL7A120
SPARC20
COL5A220
CDK620
SP119
CCNA219
FBN118
LAMC118
MCL118
AKT218
SERPINH117
VHL17
AKT316
DNMT3A16
DNMT3B16
CCND216
LOX15
DICER115
COL21A114
KDM6B14
COL10A114
COL15A113
LEPRE112
FGA11
REL10
ITGA610
FGG10
KEGG pathway enrichment, GO functional annotation and PPI network analysis for target genes of hsa-miR-29c-3p. (A) The top 10 enriched KEGG pathway items; (B) the top 10 enriched biological process (BP) items; (C) the top 10 enriched cellular component (CC) items; (D) the top 10 enriched molecular function (MF) items; (E) the top 47 hub genes (node degree >= 10) in PPT network of target genes.

Figure 5. KEGG pathway enrichment, GO functional annotation and PPI network analysis for target genes of hsa-miR-29c-3p. (A) The top 10 enriched KEGG pathway items; (B) the top 10 enriched biological process (BP) items; (C) the top 10 enriched cellular component (CC) items; (D) the top 10 enriched molecular function (MF) items; (E) the top 47 hub genes (node degree >= 10) in PPT network of target genes.

COL1A1 and COL1A2 are identified as functional target genes of DUXAP8 and DUXAP9

Next, we explored the potential functional targets of DUXAP8 and DUXAP9. Genes with higher node degree usually have stronger functions. Therefore, the top 10 hub genes (GAPDH, VEGFA, PTEN, JUN, MMP2, SIRT1, ITGB1, CDC42, COL1A1 and COL1A2) were selected for subsequent analysis. The prognostic values of the top 10 hub genes in human ccRCC were first assessed by Kaplan-Meier plotter. As shown in Table 4, high expression of GAPDH, COL1A1 and COL1A2 indicated a poor prognosis in patients with ccRCC. Four hub genes (PTEN, SIRT1, ITGB1 and CDC42) were favorable prognostic biomarkers in ccRCC. For VEGFA, JUN and MMP2, no prognostic values of them in ccRCC were discovered. Subsequently, we analyzed the expression relationship between DUXAP8/DUXAP9 and the top 10 hub genes (Table 5). The analytic result showed that six hub genes (GAPDH, PTEN, MMP2, ITGB1, COL1A1 and COL1A2) and five hub genes (PTEN, MMP2, ITGB1, COL1A1 and COL1A2) were positively correlated with DUXAP8 and DUXAP9, respectively. By combination of ceRNA hypothesis and the results from survival analysis and correlation analysis, we suggested that COL1A1 and COL1A2 may be two most potential functional targets of DUXAP8 and DUXAP9. Then, we evaluated the expression levels of the two target genes in ccRCC using two online databases, GEPIA and UALCAN. As shown in Figure 6A and Figure 6B, both COL1A1 and COL1A2 were significantly upregulated in cancer samples compared with normal samples. A similar analytic result was acquired when we determined COL1A1 and COL1A2 expression using UALCAN database (Figure 6C and Figure 6D). Furthermore, COL1A1 and COL1A2 expression among various major stages showed statistical difference as presented in Figure 6E and Figure 6F, respectively. Generally, they were markedly upregulated in advanced stages of ccRCC. Finally, the expression correlation of miR-29c-3p with COL1A1/COL1A2 was assessed. Figure 6G and Figure 6H demonstrated that miR-29c-3p was negatively associated with COL1A1 and COL1A2 expression in ccRCC.

Table 4. Prognostic values of top ten hub genes in clear cell renal cell carcinoma (ccRCC) determined by Kaplan-Meier plotter.

Gene namesP-valueHRsa95%CIbPrognostic outcomes
GAPDH3.00E-041.771.29–2.42Poor
VEGFA0.160.790.571.10Good
PTEN9.00E-050.560.410.75Good
JUN0.05870.750.551.01Good
MMP20.281.190.871.64Poor
SIRT15.00E-060.470.340.65Good
ITGB12.40E-040.570.430.78Good
CDC425.60E-080.450.330.60Good
COL1A10.000421.731.27–2.35Poor
COL1A20.00981.551.11-2.17Poor
aHRs: hazard ratios;
b95%CI: 95% confidence interval.

Table 5. Correlation analysis between DUXAP8/DUXAP9 with top 10 hub genes expression in kidney renal clear cell carcinoma using GEPIA database.

GeneDUXAP8DUXAP9
RP-valueRP-value
GAPDH0.13000.00320.09700.0270
VEGFA0.02000.65000.09600.0290
PTEN0.15000.00060.16000.0002
JUN−0.03900.37000.02900.5200
MMP20.41000.00000.36000.0000
SIRT10.02600.55000.00790.8600
ITGB10.24000.00000.25000.0000
CDC420.03600.41000.06300.1500
COL1A10.36000.00000.30000.0000
COL1A20.43000.00000.36000.0000
COL1A1 and COL1A2 were identified as two potential target genes of DUXAP8 and DUXAP9. (A) Box-whisker plot represented expression of COL1A1 in clear cell renal cell carcinoma (ccRCC or KIRC) compared with normal controls determined by GEPIA database; (B) Box-whisker plot represented expression of COL1A2 in ccRCC compared with normal controls determined by GEPIA database; (C) expression of COL1A1 in ccRCC compared with normal controls determined by UALCAN database; (D) expression of COL1A2 in ccRCC compared with normal controls determined by UALCAN database; (E) expression of COL1A1 among major stages in ccRCC determined by GEPIA database; (F) expression of COL1A2 among major stages in ccRCC determined by GEPIA database; (G) correlation analysis between has-miR-29c-3p and COL1A1 in ccRCC; (H) correlation analysis between has-miR-29c-3p and COL1A2 in ccRCC. “*” represents P-value less than 0.05. Three horizontal lines in the box plot represent minimum, median and maximum, respectively.

Figure 6. COL1A1 and COL1A2 were identified as two potential target genes of DUXAP8 and DUXAP9. (A) Box-whisker plot represented expression of COL1A1 in clear cell renal cell carcinoma (ccRCC or KIRC) compared with normal controls determined by GEPIA database; (B) Box-whisker plot represented expression of COL1A2 in ccRCC compared with normal controls determined by GEPIA database; (C) expression of COL1A1 in ccRCC compared with normal controls determined by UALCAN database; (D) expression of COL1A2 in ccRCC compared with normal controls determined by UALCAN database; (E) expression of COL1A1 among major stages in ccRCC determined by GEPIA database; (F) expression of COL1A2 among major stages in ccRCC determined by GEPIA database; (G) correlation analysis between has-miR-29c-3p and COL1A1 in ccRCC; (H) correlation analysis between has-miR-29c-3p and COL1A2 in ccRCC. “*” represents P-value less than 0.05. Three horizontal lines in the box plot represent minimum, median and maximum, respectively.

DUXAP8/DUXAP9-miR-29c-3p-COL1A1/COL1A2 axis controls growth of renal cell carcinoma

These bioinformatic analytic findings indicate that pseudogene DUXAP8 and DUXAP9 may fuel onset and development of renal cell carcinoma by targeting COL1A1 and COL1A2 through competitively binding miR-29c-3p. In this part, we conducted corresponding experimental validation. Firstly, we determined the expression levels of five constituents (DUXAP8, DUXAP9, miR-29c-3p, COL1A1 and COL1A2) in 20 pairs of clinical cancer tissues and normal tissues of renal cell carcinoma. As shown in Figure 7A, 7B, two pseudogenes DUXAP8 and DUXAP9 were significantly upregulated in cancer samples compared with normal samples. miR-29c-3p expression level (Figure 7C) was markedly lower but COL1A1 (Figure 7D) and COL1A2 (Figure 7E) expression levels were higher in cancer samples than that in normal samples. Next, luciferase reporter assay was introduced to affirm the direct bind of miR-29c-3p to DUXAP8 (Figure 7F), DUXAP9 (Figure 7G), COL1A1 (Figure 7H) and COL1A2 (Figure 7I). Subsequently, we employed siRNAs of DUXAP8/ DUXAP9 to knockdown DUXAP8/DUXAP9 and used mimic of miR-29c-3p to upregulate miR-29c-3p. The knockdown and overexpression effects were presented in Figure 7J and Figure 7K, respectively. Finally, we detected proliferative roles of DUXAP8, DUXAP9 and miR-29c-3p in renal cell carcinoma cell lines. As shown in Figure 7L, knockdown of DUXAP8 or DUXAP9 could significantly suppress cell growth and overexpression of miR-29c-3p also inhibited cell growth. Moreover, we found that cells with knockdown of DUXA8/DUXAP9 and overexpression of miR-29c-3p grew the slowest. Identical with cell counting assay, colony formation assay demonstrated downregulation of DUXAP8/ DUXAP9 or upregulation of miR-29c-3p decreased colony number (Figure 7M). Taken together, amplification of pseudogenes DUXAP8 and DUXAP9 may lead to increase expression of COL1A1 and COL1A2 by competitively binding to miR-29c-3p, resulting in uncontrolled cell proliferation in renal cell carcinoma (Figure 8).

Pseudogenes DUXAP8 and DUXAP9 promote tumor growth via suppression of miR-29c-3p-COL1A1/COL1A2 axis in renal cell carcinoma. (A) The expression of pseudogene DUXAP8 in renal cell carcinoma samples and normal controls; (B) the expression of pseudogene DUXAP9 in renal cell carcinoma samples and normal controls; (C) the expression of miR-29c-3p in renal cell carcinoma samples and normal controls; (D) the expression of COL1A1 in renal cell carcinoma samples and normal controls; (E) the expression of COL1A2 in renal cell carcinoma samples and normal controls; (F) miR-29c-3p suppressed Renilla luciferase activity of the reporters containing the wild-type but not mutant DUXAP8 in ACHN cell line; (G) miR-29c-3p suppressed Renilla luciferase activity of the reporters containing the wild-type but not mutant DUXAP9 in A498 cell line; (H) miR-29c-3p suppressed Renilla luciferase activity of the reporters containing the wild-type but not mutant COL1A1 in ACHN cell line; (I) miR-29c-3p suppressed Renilla luciferase activity of the reporters containing the wild-type but not mutant COL1A2 in A498 cell line; (J) the knockdown effect of si-DUXAP8 and si-DUXAP9 in ACHN and A498 cell lines; (K) the overexpression effect of miR-29c-3p mimic in ACHN and A498 cell lines; (L) knockdown of pseudogene DUXAP8/DUXAP9 or/and overexpression of miR-29c-3p inhibited cell growth in vitro; (M) knockdown of pseudogene DUXAP8/DUXAP9 or/and overexpression of miR-29c-3p suppressed the colony formation of cell lines in renal cell carcinoma. “*” represents P-value less than 0.05. Abbreviations: MUT, mutant; WT, wild type.

Figure 7. Pseudogenes DUXAP8 and DUXAP9 promote tumor growth via suppression of miR-29c-3p-COL1A1/COL1A2 axis in renal cell carcinoma. (A) The expression of pseudogene DUXAP8 in renal cell carcinoma samples and normal controls; (B) the expression of pseudogene DUXAP9 in renal cell carcinoma samples and normal controls; (C) the expression of miR-29c-3p in renal cell carcinoma samples and normal controls; (D) the expression of COL1A1 in renal cell carcinoma samples and normal controls; (E) the expression of COL1A2 in renal cell carcinoma samples and normal controls; (F) miR-29c-3p suppressed Renilla luciferase activity of the reporters containing the wild-type but not mutant DUXAP8 in ACHN cell line; (G) miR-29c-3p suppressed Renilla luciferase activity of the reporters containing the wild-type but not mutant DUXAP9 in A498 cell line; (H) miR-29c-3p suppressed Renilla luciferase activity of the reporters containing the wild-type but not mutant COL1A1 in ACHN cell line; (I) miR-29c-3p suppressed Renilla luciferase activity of the reporters containing the wild-type but not mutant COL1A2 in A498 cell line; (J) the knockdown effect of si-DUXAP8 and si-DUXAP9 in ACHN and A498 cell lines; (K) the overexpression effect of miR-29c-3p mimic in ACHN and A498 cell lines; (L) knockdown of pseudogene DUXAP8/DUXAP9 or/and overexpression of miR-29c-3p inhibited cell growth in vitro; (M) knockdown of pseudogene DUXAP8/DUXAP9 or/and overexpression of miR-29c-3p suppressed the colony formation of cell lines in renal cell carcinoma. “*” represents P-value less than 0.05. Abbreviations: MUT, mutant; WT, wild type.

Model of DUXAP8/DUXAP9-miR-29c-3p-COL1A1/COL1A2 axis in renal cell carcinoma.

Figure 8. Model of DUXAP8/DUXAP9-miR-29c-3p-COL1A1/COL1A2 axis in renal cell carcinoma.

Discussion

For a long time, pseudogenes are considered as non-functional or ‘junk’ DNA. Growing studies have suggested that pseudogenes act as ceRNAs of functional genes by competitively binding to shared miRNAs, thereby regulating a variety of physiological and pathological processes, including carcinogenesis [2325]. However, pseudogenes and their underlying ceRNA regulatory mechanisms in RCC remain largely unknown. To the best of our knowledge, a pseudogene-miRNA-target gene network in renal cell carcinoma is still not been constructed.

In this study, we first screened several differentially expressed pseudogenes in ccRCC using the dreamBase database, after which expression levels of these pseudogenes were further validated using another database (GEPIA). Finally, a total of 47 pseudogenes, consisting of 31 upregulated and 16 downregulated pseudogenes in ccRCC were identified. Subsequently, survival analysis and differential expression analysis among various major stages demonstrated that only high expression of pseudogene DUXAP8 and DUXAP9 indicated poor OS and RFS in ccRCC, chRCC and pRCC. DUXAP8 has been reported to function as an oncogene in glioma [27], pancreatic cancer [14], bladder cancer [28, 29], gastric cancer [30], non-small-cell lung cancer [31] as well as renal cell carcinoma [32]. DUXAP9 has also been found to promote non-small-cell lung cancer progression by directly binding with Cbl-b and thus augmenting EGFR signaling [33]. All these reports together with the previous analytic results suggest that pseudogene DUXAP8 and DUXAP9 may serve as two promising therapeutic targets ad prognostic biomarkers for RCC.

ceRNA hypothesis is an important action mechanism of ncRNAs, including pseudogenes [22]. Therefore, we predicted potential miRNAs that can theoretically bind to pseudogene DUXAP8 and DUXAP9 using starBase database. 33 and 5 potential miRNAs were predicted to bind to DUXAP8 and DUXAP9, respectively. By combination of correlation, expression and survival analysis for these miRNAs, miR-29c-3p was identified as the most potential binding miRNA of DUXAP8 and DUXAP9. It has been well documented that miR-29c-3p is an anti-tumor miRNA in multiple human cancers, such as hepatocellular carcinoma [34], breast cancer [35], gastric cancer [36], pancreatic cancer [37] as well as renal cell carcinoma [38]. Next, we predicted the downstream target genes of miR-29c-3p by miRNet, and 254 target genes were obtained. To further understand and explore the roles of DUXAP8 and DUXAP9 in RCC, pathway enrichment analysis for these target genes were performed. The results demonstrated that these target genes were significantly enriched in renal cell carcinoma and some cancer-associated pathways, including extracellular matrix organization.

It has been widely acknowledged that genes with more node degrees usually possess more functions. The top 10 hub genes among the PPI network of target genes of miR-29c-3p were selected for further analysis. Based on ceRNA mechanism, pseudogenes’ expression levels and prognostic values should be positively correlated with their functional targets. After performing survival analysis and correlation analysis, only two hub genes (COL1A1 and COL1A2) were met this requirement. The expression levels of COL1A1 and COL1A2 in ccRCC were obviously higher than that in normal controls, further suggesting that COL1A1 and COL1A2 may be two potential functional targets of DUXAP8 and DUXAP9. Previous studies have reported that COL1A1 and COL1A2 could enhance cancer progression. For example, Liu et al. suggested that COL1A1 promoted metastasis of breast cancer [39]; Zhang et al. found that COL1A1 facilitated metastasis by regulating the WNT/PCP pathway in colorectal cancer [40]; Yu et al. confirmed that COL1A2 inhibited colorectal cancer cell proliferation, migration and invasion [41]. The group of Li J also supposed that COL1A1 and COL1A2 might predict poor clinical outcomes in gastric cancer patients [42]. Moreover, Yamada et al. reported that high expression of ten target genes of miR-29 family, including COL1A1, predicted poor patient prognosis in renal cell carcinoma patients [43]. Subsequently, correlation analysis revealed that miR-29c-3p expression was inversely linked to COL1A1 and COL1A2 expression. All these findings suggest that pseudogene DUXAP8/DUXAP9-miR-29c-3p-COL1A1/COL1A2 may be a critical signaling pathway in onset and progression of renal cell carcinoma. qPCR showed that DUXAP8, DUXAP9, COL1A1 and COL1A2 were significantly upregulated in renal cell carcinoma cancer samples compared with normal controls whereas miR-29c-3p expression in cancer tissues was decreased. miR-29c-3p was confirmed to directly bind to DUXAP8, DUXAP9, COL1A1 and COL1A2 by luciferase reporter assay. Moreover, functional experiments demonstrated that knockdown of pseudogenes DUXAP8/DUXAP9 reduced proliferative ability of renal cell carcinoma. All these findings were in accordance with our previous analytic results, further partially suggesting accuracy of these bioinformatic analyses

Although in silico analyses and preliminary experimental validation have shown potential oncogenic roles of DUXAP8 and DUXAP9 and its underlying ceRNA regulatory mechanism in renal cell carcinoma, much more experiments and big clinical trials need to be launched to further confirm these preliminary data in the future.

Materials and Methods

Screening for differentially expressed pseudogene

RNA-seq expression data of pseudogenes in human ccRCC were directly downloaded from the dreamBase database (http://rna.sysu.edu.cn/dreamBase/index.php). DreamBase database is an integrated platform for analyzing regulatory features of pseudogenes from multi-dimensional high-throughput sequencing data [44]. |log2FC| > 2.0 was set as the cut-off criterion for differentially expressed pseudogenes. FC = Tumor/ Normal. Subsequently, they were validated by pseudogene names and expression levels using dreamBase database again. These screened pseudogenes were considered as candidate pseudogenes.

Confirmation pseudogene expression

A newly developed interactive tool, named Gene Expression Profiling Interactive Analysis (GEPIA) (http://gepia.cancer-pku.cn/detail.php), provides customizable functions such as tumor/normal differential expression analysis, profiling according to cancer types or pathological stages, patient survival analysis, similar gene detection, correlation analysis and dimensionality reduction analysis [45].

In this study, GEPIA was utilized to further determine expression levels of candidate pseudogenes in human KIRC. P-value < 0.05 and |log2FC| > 1.0 were set as the thresholds for identifying potential pseudogenes. Expression differences among various major stages of these potential pseudogenes were also studied using this database. The target name was entered in the “stage plot” page, with selection of “Yes” for “use major stage” and “log scale”. After click the “Plot”, the stage plot and corresponding F-value and P-value were automatically generated. P-value < 0.05 was considered as statistically significant.

Kaplan-Meier survival analysis

The prognostic values (OS, overall survival; RFS, disease free survival) of potential pseudogenes in patients with ccRCC, chRCC and pRCC were assessed via GEPIA database (http://gepia.cancer-pku.cn/detail.php) [45]. Longrank P-value < 0.05 was considered as statistically significant. The pseudogenes with significant OS and RFS in ccRCC, chRCC and pRCC were re-defined as key pseudogenes and were chosen for subsequent investigation.

Prediction of miRNA binding to pseudogene

The potential miRNA-pseudogene interactions were predicted by starBase (version 3) (http://starbase.sysu.edu.cn/), which is an open-source platform for studying the miRNA-ncRNA, miRNA-mRNA, ncRNA-RNA, RNA-RNA, RBP-ncRNA and RBP-mRNA interactions from CLIP-seq, degradome-seq and RNA-RNA interactome data [46, 47]. The predicted miRNAs could be directly exported as excel, after which a visual pseudogene-miRNA network was established using Cytoscape software (version 3.6.0).

Correlation analysis between pseudogene/gene and miRNA

StarBase (version 3) (http://starbase.sysu.edu.cn/) was also utilized to perform correlation analysis between pseudogene and miRNA in ccRCC [46, 47]. The pseudogene-miRNA interactions with R < -0.1 AND P-value < 0.05 was considered as potential items. Accordingly, these miRNAs involved in the potential interactions were defined as potential miRNAs. This database was also employed to analyze the expression association between miRNA and gene. Similarly, miRNA-gene interactions with R < -0.1 and P-value < 0.05 were regarded as significant.

Expression analysis of potential miRNA

The expression levels of potential miRNAs in ccRCC, chRCC and pRCC were also analyzed by starBase (version 3) (http://starbase.sysu.edu.cn/) [46, 47]. P < 0.05 was considered as statistically significant.

Survival analysis of potential miRNA

The prognostic values of miRNAs in ccRCC and pRCC were determined through Kaplan-Meier plotter database (http://kmplot.com/analysis/), which is capable to access the effect of 11,000 miRNAs from 20 different cancer types, including clear cell carcinoma [48, 49]. Simply, the miRNAs were first typed into this database. Then, Kaplan Meier survival plots were directly created, and logrank P-value, hazard ratio (HR) and 95% confidence interval (CI) were displayed on the webpage. Logrank P-value < 0.05 was considered as statistically significant.

Prediction of potential target gene of miRNA

miRNet is an easy-to-use tool with comprehensive support for statistical analysis and functional interpretation of data generated in miRNAs studies. In this database, the miRNA target gene data were collected from several well-annotated databases, including miRTarBase v7.0, TarBase v7.0 and miRecords [50, 51]. Therefore, miRNet was selected to predict potential downstream target genes of potential miRNAs.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis

Enrichment analysis for predicted target genes was conducted using Enrichr database (http://amp.pharm.mssm.edu/Enrichr/), which is a comprehensive gene set enrichment analysis web server [52, 53]. KEGG pathway analysis and three categories of GO functional annotation, including biological process (BP), cellular component (CC) and molecular function (MF), were included in this study. The bar charts of BP, CC, MF and KEGG were directly downloaded from the webpage.

Establishment and analysis of protein-protein interaction (PPI) network

STRING database (http://string-db.org) is a database of known and predicted PPIs, in which interactions are derived from five main sources: genomic context predictions, high-throughput lab experiments, (conserved) co-expression, automated text mining and previous knowledge in databases [54, 55]. PPI network of predicted target genes was constructed using STRING database. Firstly, these predicted targets were mapped into this database. Then, gene pairs with combined score >= 0.4 were exported. Subsequently, these gene pairs were re-loaded into Cytoscape software (version 3.6.0), after which the hub genes among the PPI network were identified using Cytohubba plugin of Cytoscape software (version 3.6.0) as previously described [56].

Correlation analysis between pseudogene and hub gene

To determine the most potential target of pseudogene, the correlation analysis between pseudogene and hub gene in ccRCC was performed via GEPIA database (http://gepia.cancer-pku.cn/detail.php) [45]. Pseudogene -hub gene pairs with R > 0.1 and P-value < 0.05 were considered as potential pairs and were selected for further analysis.

Survival analysis of hub gene

The prognostic values of hub genes in ccRCC were determined using Kaplan-Meier plotter database (http://kmplot.com/analysis/) as mentioned above. Logrank P-value < 0.05 was considered as statistically significant.

Expression analysis of hub gene

Two online databases, named GEPIA and UALCAN, were introduced to detect the expression levels of hub genes in ccRCC [45, 57]. The expression differences among different major stages of hub genes in clear cell carcinoma were also evaluated by GEPIA database. P-value < 0.05 was regarded as statistically significant.

Cell culture and clinical samples

The human renal cell carcinoma cell lines ACHN and A498 were purchased from ATCC. The cells were maintained in minimal essential medium (MEM) supplemented with 10% FBS, penicillin (100 U/ml) and streptomycin (100ug/ml). 20 pairs of clear cell renal cell carcinoma tumor tissues and adjacent normal kidney tissues were collected from the First Affiliated Hospital of Zhejiang University. All procedures performed in this study involving human participants were conducted in accordance with the ethical standards of the First Affiliated Hospital of Zhejiang University and written informed consent from every participant was obtained.

Primer, RNA oligoribonucleotide, miRNA mimic and vectors

ALL primers, RNA oligoribonucleotides and miRNA mimic used in this study were designed and purchased from RiboBio Co. Ltd (Guangzhou, China).

Analysis of gene or miRNA expression

The expression level of DUXAP8, DUXAP9, miR- 29c-3p, COL1A1 and COL1A2 was analyzed by quantitative real-time PCR as we previously described [58].

Cell transfection

RNA oligios and miRNA mimic were transfected into cells with a final concentration of 50 nM using LipofectamineTM 3000 (Invitrogen) according to manufacturer’s instruction.

Luciferase reporter assay

To determine if miIR-29c-3p can directly bind to DUXAP8, DUXAP9, COL1A1 and COL1A2, luciferase reporter assay was performed as we described previously [59].

Cell counting assay

Cell counting assay was employed to evaluate cell growth. Cells were first transfected as mentioned above. At 24 h post-transfection, 4 x 104 cells were re-seeded into 12-well plates. After culturing for 3 days, cell numbers were counted.

Colony formation assay

Twenty-four hours after transfection, 500 of transfected cells (ACHN and A498) were re-seeded into 6 well-plates and cultured for two weeks. Two weeks later, colonies were fixed in methanol for 30 minutes and then stained with a 0.1% crystal violet solution for 15 minutes. Finally, visible colonies were counted.

Statistical analysis

Most of the statistical analysis has been done by the bioinformatic tools mentioned above. Experimental results were shown as mean±SD. Paired Student’s t-test was used to evaluate expression differences of genes or miRNA between paratumor group and tumor group. Differences between two groups were estimated by unpaired Student’s t-test. A two-tailed value of P < 0.05 was considered as statistically significant.

Conclusions

To sum up, this study for the first time comprehensively and systematically explored the function and ceRNA regulatory mechanism of pseudogenes in renal cell carcinoma. Our current results indicated that DUXAP8 and DUXAP9 may act as two key oncogenes in renal cell carcinoma through upregulation of COL1A1 and COL1A2 by competitively binding miR-29c-3p. Moreover, our results suggested that DUXAP8 and DUXAP9 and other components in this pathway possessed significant prognostic values in patients with kidney cancer. All these findings are worthy of more experimental validation in the future.

Author Contributions

JC and XW conceived this project. JC wrote the first manuscript. JC, WYL and BSD performed analysis. JC, WYL and XW polished the manuscript. All authors read and approved the final manuscript submitted for publication.

Conflicts of Interest

The authors have no conflicts of interest to declare.

Funding

This project was supported by the phosphorus project of the first hospital of Jiaxing (2018QMX008).

References

  • 1. Ferlay J, Soerjomataram I, Dikshit R, Eser S, Mathers C, Rebelo M, Parkin DM, Forman D, Bray F. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. Int J Cancer. 2015; 136:E359–86. https://doi.org/10.1002/ijc.29210 [PubMed]
  • 2. Kuthi L, Jenei A, Hajdu A, Németh I, Varga Z, Bajory Z, Pajor L, Iványi B. Prognostic Factors for Renal Cell Carcinoma Subtypes Diagnosed According to the 2016 WHO Renal Tumor Classification: a Study Involving 928 Patients. Pathol Oncol Res. 2017; 23:689–98. https://doi.org/10.1007/s12253-016-0179-x [PubMed]
  • 3. Yuan J, Dong R, Liu F, Zhan L, Liu Y, Wei J, Wang N. The miR-183/182/96 cluster functions as a potential carcinogenic factor and prognostic factor in kidney renal clear cell carcinoma. Exp Ther Med. 2019; 17:2457–64. https://doi.org/10.3892/etm.2019.7221 [PubMed]
  • 4. Porta C, Cosmai L, Leibovich BC, Powles T, Gallieni M. The adjuvant treatment of kidney cancer: a multidisciplinary outlook. Nat Rev Nephrol. 2019; 15:423–433. https://doi.org/10.1038/s41581-019-0131-x [PubMed]
  • 5. Sasidharan R, Gerstein M. Genomics: protein fossils live on as RNA. Nature. 2008; 453:729–31. https://doi.org/10.1038/453729a [PubMed]
  • 6. Vanin EF. Processed pseudogenes. Characteristics and evolution. Biochim Biophys Acta. 1984; 782:231–41. https://doi.org/10.1016/0167-4781(84)90057-5 [PubMed]
  • 7. Graça T, Ku PS, Silva MG, Turse JE, Hammac GK, Brown WC, Palmer GH, Brayton KA. Segmental Variation in a Duplicated msp2 Pseudogene Generates Anaplasma marginale Antigenic Variants. Infect Immun. 2019; 87:e00727–18. https://doi.org/10.1128/IAI.00727-18 [PubMed]
  • 8. Molineris I, Sales G, Bianchi F, Di Cunto F, Caselle M. A new approach for the identification of processed pseudogenes. J Comput Biol. 2010; 17:755–765. https://doi.org/10.1089/cmb.2009.0027 [PubMed]
  • 9. Song H, Yang J, Zhang Y, Zhou J, Li Y, Hao X. Integrated analysis of pseudogene RP11-564D11.3 expression and its potential roles in hepatocellular carcinoma. Epigenomics. 2019; 11:267–80. https://doi.org/10.2217/epi-2018-0152 [PubMed]
  • 10. Jacq C, Miller JR, Brownlee GG. A pseudogene structure in 5S DNA of Xenopus laevis. Cell. 1977; 12:109–20. https://doi.org/10.1016/0092-8674(77)90189-1 [PubMed]
  • 11. Chen Z, Huang Z, Chen LX. The Olfactory Receptor Pseudo-pseudogene: A Potential Therapeutic Target in Human Diseases. Biomed Environ Sci. 2018; 31:168–70. https://doi.org/10.3967/bes2018.022 [PubMed]
  • 12. Nuerzhati Y, Dong R, Song Z, Zheng S. Role of the long non-coding RNA-Annexin A2 pseudogene 3/Annexin A2 signaling pathway in biliary atresia-associated hepatic injury. Int J Mol Med. 2019; 43:739–48. https://doi.org/10.3892/ijmm.2018.4023 [PubMed]
  • 13. Shang J, Wang Z, Chen W, Yang Z, Zheng L, Wang S, Li S. Pseudogene CHIAP2 inhibits proliferation and invasion of lung adenocarcinoma cells by means of the WNT pathway. J Cell Physiol. 2019; 234:13735–46. https://doi.org/10.1002/jcp.28053 [PubMed]
  • 14. Lian Y, Yang J, Lian Y, Xiao C, Hu X, Xu H. DUXAP8, a pseudogene derived lncRNA, promotes growth of pancreatic carcinoma cells by epigenetically silencing CDKN1A and KLF2. Cancer Commun. 2018; 38:64. https://doi.org/10.1186/s40880-018-0333-9 [PubMed]
  • 15. Zhang Y, Li Y, Wang J, Lei P. Long non-coding RNA ferritin heavy polypeptide 1 pseudogene 3 controls glioma cell proliferation and apoptosis via regulation of the microRNA-224-5p/tumor protein D52 axis. Mol Med Rep. 2018; 18:4239–46. https://doi.org/10.3892/mmr.2018.9491 [PubMed]
  • 16. Lynn H, Sun X, Ayshiev D, Siegler JH, Rizzo AN, Karnes JH. Single nucleotide polymorphisms in the MYLKP1 pseudogene are associated with increased colon cancer risk in African Americans. PLoS One. 2018; 13:e0200916. https://doi.org/10.1371/journal.pone.0200916 [PubMed]
  • 17. Guo X, Deng L, Deng K, Wang H, Shan T, Zhou H, Liang Z, Xia J, Li C. Pseudogene PTENP1 Suppresses Gastric Cancer Progression by Modulating PTEN. Anticancer Agents Med Chem. 2016; 16:456–64. https://doi.org/10.2174/1871520615666150507121407 [PubMed]
  • 18. Ma H, Ma T, Chen M, Zou Z, Zhang Z. The pseudogene-derived long non-coding RNA SFTA1P suppresses cell proliferation, migration, and invasion in gastric cancer. Biosci Rep. 2018; 38:BSR20171193. https://doi.org/10.1042/BSR20171193 [PubMed]
  • 19. Tan L, Mai D, Zhang B, Jiang X, Zhang J, Bai R, Ye Y, Li M, Pan L, Su J, Zheng Y, Liu Z, Zuo Z, et al. PIWI-interacting RNA-36712 restrains breast cancer progression and chemoresistance by interaction with SEPW1 pseudogene SEPW1P RNA. Mol Cancer. 2019; 18:9. https://doi.org/10.1186/s12943-019-0940-3 [PubMed]
  • 20. Yndestad S, Austreid E, Skaftnesmo KO, Lønning PE, Eikesdal HP. Divergent Activity of the Pseudogene PTENP1 in ER-Positive and Negative Breast Cancer. Mol Cancer Res. 2018; 16:78–89. https://doi.org/10.1158/1541-7786.MCR-17-0207 [PubMed]
  • 21. Yu G, Yao W, Gumireddy K, Li A, Wang J, Xiao W, Chen K, Xiao H, Li H, Tang K, Ye Z, Huang Q, Xu H. Pseudogene PTENP1 functions as a competing endogenous RNA to suppress clear-cell renal cell carcinoma progression. Mol Cancer Ther. 2014; 13:3086–97. https://doi.org/10.1158/1535-7163.MCT-14-0245 [PubMed]
  • 22. Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011; 146:353–58. https://doi.org/10.1016/j.cell.2011.07.014 [PubMed]
  • 23. Wang L, Guo ZY, Zhang R, Xin B, Chen R, Zhao J, Wang T, Wen WH, Jia LT, Yao LB, Yang AG. Pseudogene OCT4-pg4 functions as a natural micro RNA sponge to regulate OCT4 expression by competing for miR-145 in hepatocellular carcinoma. Carcinogenesis. 2013; 34:1773–81. https://doi.org/10.1093/carcin/bgt139 [PubMed]
  • 24. Rutnam ZJ, Du WW, Yang W, Yang X, Yang BB. The pseudogene TUSC2P promotes TUSC2 function by binding multiple microRNAs. Nat Commun. 2014; 5:2914. https://doi.org/10.1038/ncomms3914 [PubMed]
  • 25. Chan JJ, Kwok ZH, Chew XH, Zhang B, Liu C, Soong TW, Yang H, Tay Y. A FTH1 gene:pseudogene:microRNA network regulates tumorigenesis in prostate cancer. Nucleic Acids Res. 2018; 46:1998–2011. https://doi.org/10.1093/nar/gkx1248 [PubMed]
  • 26. Karreth FA, Pandolfi PP. ceRNA cross-talk in cancer: when ce-bling rivalries go awry. Cancer Discov. 2013; 3:1113–21. https://doi.org/10.1158/2159-8290.CD-13-0202 [PubMed]
  • 27. Zhao X, Hao S, Wang M, Xing D, Wang C. Knockdown of pseudogene DUXAP8 expression in glioma suppresses tumor cell proliferation. Oncol Lett. 2019; 17:3511–16. https://doi.org/10.3892/ol.2019.9994 [PubMed]
  • 28. Jiang B, Hailong S, Yuan J, Zhao H, Xia W, Zha Z, Bin W, Liu Z. Identification of oncogenic long noncoding RNA SNHG12 and DUXAP8 in human bladder cancer through a comprehensive profiling analysis. Biomed Pharmacother. 2018; 108:500–507. https://doi.org/10.1016/j.biopha.2018.09.025 [PubMed]
  • 29. Lin MG, Hong YK, Zhang Y, Lin BB, He XJ. Mechanism of lncRNA DUXAP8 in promoting proliferation of bladder cancer cells by regulating PTEN. Eur Rev Med Pharmacol Sci. 2018; 22:3370–77. https://doi.org/10.26355/eurrev_201806_15158 [PubMed]
  • 30. Ma HW, Xie M, Sun M, Chen TY, Jin RR, Ma TS, Chen QN, Zhang EB, He XZ, De W, Zhang ZH. The pseudogene derived long noncoding RNA DUXAP8 promotes gastric cancer cell proliferation and migration via epigenetically silencing PLEKHO1 expression. Oncotarget. 2016; 8:52211–24. https://doi.org/10.18632/oncotarget.11075 [PubMed]
  • 31. Sun M, Nie FQ, Zang C, Wang Y, Hou J, Wei C, Li W, He X, Lu KH. The Pseudogene DUXAP8 Promotes Non-small-cell Lung Cancer Cell Proliferation and Invasion by Epigenetically Silencing EGR1 and RHOB. Mol Ther. 2017; 25:739–751. https://doi.org/10.1016/j.ymthe.2016.12.018 [PubMed]
  • 32. Huang T, Wang X, Yang X, Ji J, Wang Q, Yue X, Dong Z. Long Non-Coding RNA DUXAP8 Enhances Renal Cell Carcinoma Progression via Downregulating miR-126. Med Sci Monit. 2018; 24:7340–47. https://doi.org/10.12659/MSM.910054 [PubMed]
  • 33. Zhu T, An S, Choy MT, Zhou J, Wu S, Liu S, Liu B, Yao Z, Zhu X, Wu J, He Z. LncRNA DUXAP9-206 directly binds with Cbl-b to augment EGFR signaling and promotes non-small cell lung cancer progression. J Cell Mol Med. 2019; 23:1852–1864. https://doi.org/10.1111/jcmm.14085 [PubMed]
  • 34. Wang CM, Wang Y, Fan CG, Xu FF, Sun WS, Liu YG, Jia JH. miR-29c targets TNFAIP3, inhibits cell proliferation and induces apoptosis in hepatitis B virus-related hepatocellular carcinoma. Biochem Biophys Res Commun. 2011; 411:586–92. https://doi.org/10.1016/j.bbrc.2011.06.191 [PubMed]
  • 35. Cittelly DM, Finlay-Schultz J, Howe EN, Spoelstra NS, Axlund SD, Hendricks P, Jacobsen BM, Sartorius CA, Richer JK. Progestin suppression of miR-29 potentiates dedifferentiation of breast cancer cells via KLF4. Oncogene. 2013; 32:2555–64. https://doi.org/10.1038/onc.2012.275 [PubMed]
  • 36. Wang L, Yu T, Li W, Li M, Zuo Q, Zou Q, Xiao B. The miR-29c-KIAA1199 axis regulates gastric cancer migration by binding with WBP11 and PTP4A3. Oncogene. 2019; 38:3134–50. https://doi.org/10.1038/s41388-018-0642-0 [PubMed]
  • 37. Lu Y, Hu J, Sun W, Li S, Deng S, Li M. MiR-29c inhibits cell growth, invasion, and migration of pancreatic cancer by targeting ITGB1. Onco Targets Ther. 2015; 9:99–109. https://doi.org/10.2147/OTT.S92758 [PubMed]
  • 38. Qu Y, Xiao H, Xiao W, Xiong Z, Hu W, Gao Y, Ru Z, Wang C, Bao L, Wang K, Ruan H, Song Z, Chen K, et al. Upregulation of MIAT Regulates LOXL2 Expression by Competitively Binding MiR-29c in Clear Cell Renal Cell Carcinoma. Cell Physiol Biochem. 2018; 48:1075–1087. https://doi.org/10.1159/000491974 [PubMed]
  • 39. Liu J, Shen JX, Wu HT, Li XL, Wen XF, Du CW, Zhang GJ. Collagen 1A1 (COL1A1) promotes metastasis of breast cancer and is a potential therapeutic target. Discov Med. 2018; 25:211–23. [PubMed]
  • 40. Zhang Z, Wang Y, Zhang J, Zhong J, Yang R. COL1A1 promotes metastasis in colorectal cancer by regulating the WNT/PCP pathway. Mol Med Rep. 2018; 17:5037–42. https://doi.org/10.3892/mmr.2018.8533 [PubMed]
  • 41. Yu Y, Liu D, Liu Z, Li S, Ge Y, Sun W, Liu B. The inhibitory effects of COL1A2 on colorectal cancer cell proliferation, migration, and invasion. J Cancer. 2018; 9:2953–62. https://doi.org/10.7150/jca.25542 [PubMed]
  • 42. Li J, Ding Y, Li A. Identification of COL1A1 and COL1A2 as candidate prognostic factors in gastric cancer. World J Surg Oncol. 2016; 14:297. https://doi.org/10.1186/s12957-016-1056-5 [PubMed]
  • 43. Yamada Y, Sugawara S, Arai T, Kojima S, Kato M, Okato A, Yamazaki K, Naya Y, Ichikawa T, Seki N. Molecular pathogenesis of renal cell carcinoma: Impact of the anti-tumor miR-29 family on gene regulation. Int J Urol. 2018; 25:953–965. https://doi.org/https://doi.org/10.1111/iju.13783 [PubMed]
  • 44. Zheng LL, Zhou KR, Liu S, Zhang DY, Wang ZL, Chen ZR, Yang JH, Qu LH. dreamBase: DNA modification, RNA regulation and protein binding of expressed pseudogenes in human health and disease. Nucleic Acids Res. 2018; 46:D85–91. https://doi.org/10.1093/nar/gkx972 [PubMed]
  • 45. Tang Z, Li C, Kang B, Gao G, Li C, Zhang Z. GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 2017; 45:W98–102. https://doi.org/10.1093/nar/gkx247 [PubMed]
  • 46. Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 2014; 42:D92–97. https://doi.org/10.1093/nar/gkt1248 [PubMed]
  • 47. Yang JH, Li JH, Shao P, Zhou H, Chen YQ, Qu LH. starBase: a database for exploring microRNA-mRNA interaction maps from Argonaute CLIP-Seq and Degradome-Seq data. Nucleic Acids Res. 2011; 39:D202–09. https://doi.org/10.1093/nar/gkq1056 [PubMed]
  • 48. Lou W, Chen J, Ding B, Chen D, Zheng H, Jiang D, Xu L, Bao C, Cao G, Fan W. Identification of invasion-metastasis-associated microRNAs in hepatocellular carcinoma based on bioinformatic analysis and experimental validation. J Transl Med. 2018; 16:266. https://doi.org/10.1186/s12967-018-1639-8 [PubMed]
  • 49. Wang W, Lou W, Ding B, Yang B, Lu H, Kong Q, Fan W. A novel mRNA-miRNA-lncRNA competing endogenous RNA triple sub-network associated with prognosis of pancreatic cancer. Aging (Albany NY). 2019; 11:2610–27. https://doi.org/10.18632/aging.101933 [PubMed]
  • 50. Fan Y, Siklenka K, Arora SK, Ribeiro P, Kimmins S, Xia J. miRNet - dissecting miRNA-target interactions and functional associations through network-based visual analysis. Nucleic Acids Res. 2016; 44:W135–41. https://doi.org/10.1093/nar/gkw288 [PubMed]
  • 51. Fan Y, Xia J. miRNet-Functional Analysis and Visual Exploration of miRNA-Target Interactions in a Network Context. Methods Mol Biol. 2018; 1819:215–33. https://doi.org/10.1007/978-1-4939-8618-7_10 [PubMed]
  • 52. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, Koplev S, Jenkins SL, Jagodnik KM, Lachmann A, McDermott MG, Monteiro CD, Gundersen GW, Ma’ayan A. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016; 44:W90–7. https://doi.org/10.1093/nar/gkw377 [PubMed]
  • 53. Lou W, Ding B, Fan W. High Expression of Pseudogene PTTG3P Indicates a Poor Prognosis in Human Breast Cancer. Mol Ther Oncolytics. 2019; 14:15–26. https://doi.org/10.1016/j.omto.2019.03.006 [PubMed]
  • 54. Lou W, Ding B, Xu L, Fan W. Construction of Potential Glioblastoma Multiforme-Related miRNA-mRNA Regulatory Network. Front Mol Neurosci. 2019; 12:66. https://doi.org/10.3389/fnmol.2019.00066 [PubMed]
  • 55. Szklarczyk D, Franceschini A, Kuhn M, Simonovic M, Roth A, Minguez P, Doerks T, Stark M, Muller J, Bork P, Jensen LJ, von Mering C. The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res. 2011; 39:D561–68. https://doi.org/10.1093/nar/gkq973 [PubMed]
  • 56. Lou W, Liu J, Ding B, Xu L, Fan W. Identification of chemoresistance-associated miRNAs in breast cancer. Cancer Manag Res. 2018; 10:4747–57. https://doi.org/10.2147/CMAR.S172722 [PubMed]
  • 57. Chandrashekar DS, Bashel B, Balasubramanya SA, Creighton CJ, Ponce-Rodriguez I, Chakravarthi BV, Varambally S. UALCAN: A Portal for Facilitating Tumor Subgroup Gene Expression and Survival Analyses. Neoplasia. 2017; 19:649–58. https://doi.org/10.1016/j.neo.2017.05.002 [PubMed]
  • 58. Lou W, Liu J, Ding B, Chen D, Xu L, Ding J, Jiang D, Zhou L, Zheng S, Fan W. Identification of potential miRNA-mRNA regulatory network contributing to pathogenesis of HBV-related HCC. J Transl Med. 2019; 17:7. https://doi.org/10.1186/s12967-018-1761-7 [PubMed]
  • 59. Chen D, Si W, Shen J, Du C, Lou W, Bao C, Zheng H, Pan J, Zhong G, Xu L, Fu P, Fan W. miR-27b-3p inhibits proliferation and potentially reverses multi-chemoresistance by targeting CBLB/GRB2 in breast cancer cells. Cell Death Dis. 2018; 9:188. https://doi.org/10.1038/s41419-017-0211-4 [PubMed]