Comprehensive analysis of the long non-coding RNA-associated competitive endogenous RNA network reveals novel prognostic biomarkers in Wilms' tumor.

Wilms' tumor (WT) is one of the most common types of renal carcinoma in children. The aim of the present study was to construct a competitive endogenous RNA (ceRNA) regulation network and explore novel prognostic biomarkers for WT. The expression profiles were downloaded from The Cancer Genome Atlas database to identify differentially expressed RNAs (DERNAs). Based on the interactions between microRNAs (miRNAs) and mRNAs/long non-coding RNAs (lncRNAs), a ceRNA network was constructed. Functional enrichment analyses were subsequently conducted to explore the functions of the ceRNA-associated DEmRNAs. Survival analysis was performed to screen for prognosis-associated RNAs and the χ2 test was used to assess the associations between prognosis-associated RNA expression and histology classification/clinical staging. The present study identified 1,784 lncRNAs, 114 miRNAs and 3,337 mRNAs, which were abnormally expressed in WT compared with that in normal samples. By prediction, pairing and network analysis, a ceRNA network consisting of 38 DElncRNAs, 18 DEmiRNAs and 99 DEmRNAs was established. These DEmRNAs were significantly enriched in pathways associated with the occurrence and development of WT. By combining the expression data with survival analysis, seven prognosis-associated RNAs were identified (P<0.05). Of these seven RNAs, two (zinc finger and BTB domain containing 4; and deleted in lymphocytic leukemia 2) were significantly associated with clinical staging and histology classification. Lastly, the expression levels of the seven RNAs were verified in the Gene Expression Omnibus database. The present study revealed that 7 RNAs might be considered as novel prognostic biomarkers and potential treatment targets for therapy in WT. In addition, the ceRNA regulation network could provide novel strategies for further studies on lncRNAs and miRNAs in WT. Copyright: © Liu et al.

Introduction

Wilms' tumor (WT), also known as nephroblastoma, is a malignant form of renal cancer that originates from the metanephric blastema and is one of the most common types of pediatric cancer (1). The incidence rate of nephroblastoma makes it the fifth most common type of pediatric malignant neoplasm, and ~90% of childhood renal tumors are of the WT in 2001 (2,3). With advancements in diagnosis and therapy, the survival rate of children with nephroblastoma has greatly improved and the 5-year survival rate is >85% in European children (1978–1997) (4,5). However, some children with WT still have a much poorer prognosis, such as those with bilateral WT and focal anaplastic WT, and cannot be treated completely; thus having a poor clinical outcome due to tumor recurrence and metastasis (6,7). Therefore, there is an urgent requirement to identify potential biomarkers or therapeutic targets for WT.

Recent evidence suggests that microRNAs (miRNAs) and long non-coding RNAs (lncRNAs) could play an important role in the occurrence and development of WT (8–10). lncRNAs are non-coding RNA transcripts of >200 nucleotides in length (11). In human cells, lncRNAs are widely expressed and play a key regulatory role in cellular activity (12). miRNAs are a class of non-coding small RNA molecules that can inhibit the expression of target genes by interacting with miRNA response elements (13). Cao et al (14) demonstrated that Stat3 can promote cell proliferation through the upregulation of miRNA-370, which directly inhibits the expression of Wilms’ tumor gene on X chromosome (WTX), a tumor repressor of WT. Jiang and Li (15) found that p73 was the target mRNA of miR-1180 and that the inhibition of miR-1180 could promote apoptosis in the SK-NEP-1 cell line and inhibit tumor growth in mice. Low expression levels of miR-21 in SK-NEP-1 cells were found to promote cellular proliferation and migration by targeting PTEN, which is a tumor suppressor (16). Liu et al (17) demonstrated that the knockdown of miR-19b can inhibit the proliferation, invasion and migration of SK-NEP-1 WT cells. However, up to now, only 4 lncRNAs (8,14,18,19) and 10 miRNAs (15–17,20,21) have been intensively studied in WT, and the function and mechanism of most of these remain unknown.

In 2011, Salmena et al (22) proposed the competitive endogenous (ceRNA) regulatory network hypothesis. This hypothesis suggests that lncRNAs could not only directly participate in the expression of the target regulatory gene, but also may contain the seed sequence of the core of the miRNAs. The target gene can be further regulated by adsorbing the corresponding miRNA, thereby affecting the number and abundance of miRNAs, ultimately affecting gene expression. With the development of high-throughput sequencing technology, an lncRNA-miRNA-mRNA ceRNA regulatory network has been constructed for various disease types, such as gastric cancer, colorectal cancer, breast cancer, liver cancer and periodontitis (23–27). However, the construction of a ceRNA network based on high-throughput sequencing has not yet been generated for WT. The present study aimed to construct a ceRNA regulatory network by investigating the associations between lncRNA, miRNAs and mRNAs, and to identify candidate prognostic biomarkers based on this network.

The Cancer Genome Atlas (TCGA) project was constructed to understand the causes and pathogenesis of cancer from a molecular perspective, thereby further improving early diagnosis, treatment and, ultimately, cancer prevention (28). In the present study, RNAsequencing (seq) and miRNAseq datasets were downloaded from TCGA database. The data obtained were processed using the edgeR software, in order to identify differentially expressed (DE)RNAs. Subsequently, DERNAs were predicted and integrated to construct a lncRNA-miRNA-mRNA ceRNA network for WT. Survival analysis was further conducted to identify prognostic biomarkers in the ceRNA network. The χ2 test was used to assess the association between the expression of prognostic RNA and histology classification/clinical staging. The present study aids our further understanding of the molecular basis of WT, as well as the discovery of potential prognostic markers for diagnosis and treatment.

Materials and methods

Data download and pre-processing

The RNA-sequencing (RNA-seq) data and the miRNA sequencing data were downloaded from the TCGA data portal [https://tcga-data.nci.nih.gov/tcga/; Data Release v.15.0; release time: Feb 20, 2019; DbGaP (The database of Genotypes and Phenotypes) study accession, phs000218). The mRNA sequencing data included 120 WT malignant and 6 adjacent normal tissues, the miRNA sequencing data included 126 WT malignant tissues and 6 adjacent normal tissues. The GENCODE database is currently the main genome annotation database (29). The GENCODE (22nd edition) GTF file was used to annotate the expression profile of RNAseq files and quantify lncRNAs and mRNAs, RNA not included in the GENCODE database was excluded. Next, mRNA and lncRNA expression profiles were extracted from the RNAseq expression matrix. Thus, three expression profiles were obtained for mRNA, lncRNA and miRNA. The clinical follow-up datasets from 126 patients with WT, including clinical characteristics such as sex, age, ethnicity, pathological stage, histology classification and survival status, were also obtained from the TCGA database. The clinical characteristics of patients with WT are summarized in Table I.

Table I.

Clinical characteristics of all patients (n=128) in The Cancer Genome Atlas cohort.

Characteristics	Value
Age at diagnosis, days
Mean	1,688
Range	156-5,698
Sex, n (%)
Male	54 (42.2)
Female	74 (57.8)
Ethnicity, n (%)
White	95 (74.2)
Black or African American	19 (14.8)
Other (Asian, Native Hawaiian, American Indian)	5 (3.9)
Not reported	9 (7.0)
Clinical stage, n (%)
I	17 (13.3)
II	49 (38.3)
III	47 (36.7)
IV	14 (10.9)
IV/V	1 (0.8)
Histology classification of primary tumor, n (%)
FAWT	84 (65.6)
DAWT	44 (34.4)
Adverse event
Relapse	94 (73.4)
None	27 (21.1)
Progression	7 (5.5)
Overall survival time, days
Mean	1,728
Range	180-4,795
Vital status, n (%)
Alive	76 (59.4)
Deaths	52 (40.6)

FAWT, focal anaplastic Wilms tumors; DAWT, diffuse anaplasia in Wilms tumors.

Screening of DERNAs

Expression profiles from tumor and normal samples were processed, and all data representing unexpressed RNA were firstly filtered out. The remaining data were further analyzed using the edgeR software to obtain Differentially expressed lncRNA (DElncRNA), Differentially Expressed mRNA (DEmRNA) and Differentially Expressed miRNA (DEmiRNA). EdgeR (v.3.28.0) is an R package specifically designed to analyze DE genes (30). By following the steps outlined in the edgeR operating guidelines, all P-values were corrected for multiple tests using the false discovery rate (FDR). Subsequently, the DERNAs were screened out using a specific cut-off value [FDR<0.01 and log|fold change (FC)|>2]. Hierarchical Clustering method was used to cluster DERNAs and samples. A volcano plot and heat map plot were generated using the gplots package (https://cran.r-project.org/web/packages/gplots/index.html; v.3.0.12).

Construction of the ceRNA regulatory network

The miRcode database (http://www.mircode.org/mircode/) is a web search platform, which is dedicated to the prediction of target miRNAs by uploading relevant lncRNA and miRNA (31). By uploading the DElncRNAs, miRNAs that interact with DElncRNAs were screened out, and then overlapped with DEmiRNAs in order to derive common miRNAs; the common miRNAs were paired with the corresponding DElncRNAs to derive lncRNA-miRNA pairing files. The miRDB (http://www.mirdb.org/; v.6.0) (32) is an online database for miRNA target prediction using the MirTarget bioinformatics tool (v.4.0), which was developed by analyzing thousands of miRNA-target interactions from high-throughput sequencing experiments. The miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php; version 7.0) (33) is also a miRNA target gene database, which is validated in in vivo experiments. TargetScan (http://www.targetscan.org/vert_72/; v.7.2) predicts miRNA targets by determining mRNAs with conserved sequence complementarity to the seed (nucleotides 2–7) of the input miRNA (34). The target mRNAs of the miRNAs in the lncRNA-miRNA pairs were predicted using these three databases (miRDB, miRTarBase and TargetScan), and defined target mRNAs that were predicted by all three databases, as the final screened-out target genes. Subsequently, the intersection elements were obtained between the target mRNAs and DEmRNAs and these intersection mRNAs were selected to build a miRNA-mRNA pairing file. Cytoscape (v.3.3.2) is an open source bioinformatics software platform for visualizing and analyzing molecular interaction networks (35). The NetWorkAnalyzer toolkit (v.3.3.2) from Cytoscape was used to analyze the characteristics of the ceRNA topology network, including network connections, the path length, the closest centrality of nodes, and the degree of connectivity. Based on the node degree, the hub lncRNAs (degree >5) and its associated mRNAs and miRNAs in the network were identified. It was reported that lncRNAs can positively regulate mRNAs by competitive combination with miRNAs. Therefore, correlation analysis was performed for each candidate ceRNA pair, and only pairs with correlation coefficients >0.4 and P<0.05 were selected as the final ceRNA pairs. Finally, a ceRNA topological network was constructed.

Functional and enrichment analyses of mRNAs in the ceRNA network

Gene Ontology (GO) is a database established by the Gene Ontology Consortium and aims to annotate genes and gene products from different organisms using three ontologies, including cellular component (CC), molecular function (MF) and biological process (BP) (36). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database that systematically analyzes the metabolic pathways of gene products and their functions (37). The clusterProfiler (v.3.14.3) is a bioconductor package, which can perform statistical analyses and visualization of functional clustering of gene sets or gene clusters (38). In order to explore the function of mRNAs in the ceRNA network, a GO and KEGG pathway analysis was conducted.

Protein-protein interaction (PPI) network analysis

To determine the interactive associations between DEmRNAs in the ceRNA network, a PPI network analysis was performed with the online software Search Tool for the Retrieval of Interacting Genes/Proteins (STRING (version 11.0; http://string-db.org); confidence score >0.4) (39).

Survival analysis of the ceRNA module

Clinical information was downloaded from TCGA database, and the survival data of patients were extracted and combined with the expression matrix of RNAs in the ceRNA topological network. A Kaplan-Meier (KM) survival curve was analyzed for each node in the ceRNA topological network using the survival package (v.3.1–8) in R. The median value was used as the cut-off for the gene expression value and the log-rank test was performed to determine the differences between the high and the low expression groups (40). P<0.05 was considered to indicate a statistically significant difference.

Screening of important prognosis-associated DERNAs

Following expression and survival analyses, four types of prognosis-associated DERNAs were identified: i) High expression of DERNAs were associated with a poor prognosis in patients with WT, ii) low expression of DERNAs were associated with a poor prognosis in patients with WT, iii) high expression of DERNAs were associated with a good prognosis in patients with WT and iv) low expression of DERNAs were associated with a good prognosis in patients with WT. The abnormal expression of DERNAs, which were associated with a poor prognosis in patients with WT, were selected as the prognosis-associated DERNAs for further experimentation. χ2 tests were performed to analyze the associations between the important prognosis-associated DERNAs and clinical characteristics.

Validation of the expression of prognosis-related DERNAs

In order to further confirm the expression levels of the candidate prognosis-associated DERNAs, validation datasets were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). A total of 28 samples of WT tissues and 4 adjacent non-tumor tissues were included in the GSE66405 dataset. The expression profile of GSE50505 contained 28 samples of WT tissues and 6 adjacent non-tumor tissues. The miRNA expression profiles were also obtained from the GEO database, including GSE50505 (26 WT samples and 12 adjacent non-tumor tissues) and GSE57370 (62 WT samples and 4 normal controls (41,42). All validation datasets were downloaded from the Gene Expression Omnibus (GEO) database. Unpaired t-tests were used to compare the differences in expression levels. P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of Significantly DERNAs

The clinical data, RNAseq data and miRNAseq data regarding WT were downloaded from TCGA database. Subsequently, the edgeR package was used to analyze the original expression profiles between WT and normal tissues. Based on the cut-off conditions (FDR<0.01 and log|FC|>2), a total of 3,337 DEmRNAs were screened, including 1,577 upregulated mRNAs and 1,760 downregulated mRNAs (Table SI). For lncRNAs, 1,784 DElncRNAs were identified, including 833 that were upregulated and 951 that were downregulated (Table SII). The differential expression profiles of miRNAs were also compared in tumor and normal tissues; 114 DEmiRNAs were identified, of which 49 were upregulated and 65 were downregulated (Table SIII). The volcano plot and heat maps are shown in Fig. 1.

Figure 1.

DERNAs from Wilms' tumor. Heat map of (A) DEmRNAs, (B) DElncRNAs and (C) DEmiRNAs. Each column represents a sample, and each row represents a transcript. Volcano plots of (D) DEmRNAs, (E) DElncRNAs and (F) DEmiRNAs. Upregulated transcripts are shown in red and downregulated transcripts in green. DE, differentially expressed; lncRNAs, long non-coding RNA; miRNAs, microRNAs; FC, fold change; FDR, false discovery rate.

Construction of the ceRNA network

To further the understanding of the role of these DERNAs in WT, a ceRNA network was constructed to understand the interaction between them. Firstly, 23 DEmiRNAs that interacted with DElncRNAs were predicted, based on the results from the miRcode database. The mRNAs targeted by these 23 DEmiRNAs were retrieved from the miRTarBase, TargetScan and miRDB databases. Subsequently, these targeted mRNAs were compared with the 3,337 DEmRNAs obtained from the aforementioned differential analyses. The mRNAs that were not included in the 3,337 DEmRNAs were excluded, resulting in 133 DEmRNAs and 19 DEmiRNAs in the ceRNA networks. Following this, the 19 DEmiRNAs were compared with miRNAs in the lncRNA-miRNA pairing file, resulting in 189 DElncRNAs. Finally, 735 DElncRNA-DEmiRNA pairs and 188 DEmiRNA-DEmRNA interaction pairs were identified from 189 DElncRNAs, 19 DEmiRNAs and 133 DEmRNAs. In order to verify the reliability of the network, network analyses were performed in order to understand the characteristics of the ceRNA network. With an increase in node degree, the number of nodes decreased (Fig. S1A). The closeness centrality is a measure centrality which describes how fast information spreads from a given node to other reachable nodes in the network. A number of nodes displayed the numbers of closeness centrality ~0.4, which indicated these nodes were relatively centralized within the network. Only a few nodes have closeness centrality ~0.5, which indicated those nodes were relatively sparsely distributed (Fig. S1B). The distribution of shared neighbors is shown in Fig. S1C; most nodes in the network had few shared neighbors. Shortest paths is a measure of a network's overall navigability (43). Fig. S1D demonstrates the distribution of the shortest path and the path length was relatively shorter (<4), which means the network has better navigability. As the hub genes with higher degree in biological networks were more likely to be important, a hub lncRNA (degree >5) and its linked mRNAs and miRNAs in the ceRNA network were screened. According to the correlation analyses (Cor >0.4; P<0.05; Table SIV), candidate ceRNA pairs were chosen as the final ceRNA pairs. Finally, a ceRNA topological network was reconstructed, which included 218 DElncRNA-DEmiRNA pairs and 138 DEmiRNA-DEmRNA interaction pairs from 38 DElncRNAs, 18 DEmiRNAs and 99 DEmRNAs (Table SV). Subsequently, the Cytoscape software was used to visualize this information, and the constructed ceRNA regulatory network from WT is shown in Fig. 2.

Figure 2.

Visualization of the lncRNA-mRNA-miRNA ceRNA network in Wilms' tumor. Diamonds, lncRNA; rectangles, miRNAs; ovals, mRNAs; red, upregulated RNA; purple, downregulated RNA. lncRNAs, long non-coding RNA; miRNAs, microRNAs; ce, competitive endogenous.

Functional enrichment of the DEmRNAs

In order to investigate the function of the 99 DEmRNAs in the ceRNA network and the signaling pathways involved, functional and pathway enrichment analyses were performed using the clusterProfiler package. From GO analyses, 73 GO entries (FDR <0.01) were obtained. The results showed that these pathways were mainly enriched in the ‘positive regulation of smooth muscle cell proliferation’, ‘muscle cell proliferation’ and ‘regulation of transforming growth factor β production’. MF was mainly directed to ‘transcription factor activity’. The CC mainly included ‘RNA polymerase II transcription factor complex’, ‘nuclear chromatin’, ‘nuclear transcription factor complex’ and ‘transcription factor complex’ (Fig. 3A). From KEGG pathway analyses, it was found that these DEmRNAs were significantly enriched in the pathways of the ‘cell cycle’, ‘small cell lung cancer’, ‘miRNAs in cancer’, ‘human papilloma virus infection’ and ‘bladder cancer’ (Fig. 3B). To better understand the role of DEmRNAs, a PPI network was established using the STRING online software, including 96 nodes and 298 edges (Fig. 4A), of which the main hub nodes were CHEK1 (checkpoint kinase 1), CDC25A (cell division cycle 25A), SKP2 (S-phase kinase associated protein 2), STAT3 (signal transducer and activator of transcription 3) and E2F1 (E2F transcription factor 1) (Fig. 4B).

Figure 3.

Functional and signaling pathway analyses of differentially expressed mRNAs in the competitive endogenous RNA network. (A) GO enrichment results. (B) Top 20 terms in the Kyoto Encyclopedia of Genes and Genomes pathways analysis. GO, Gene Ontology.

Figure 4.

PPI network for differentially expressed mRNAs in the competitive endogenous RNA network. (A) Visualization of the PPI network. Filled nodes represent known or predicted 3D protein structures. Empty nodes represent proteins of unknown 3D structure. Edges represent PPIs. Pink, experimentally determined; green, gene neighborhood; light blue, manually curated; dark blue, gene co-occurrence; red, gene fusions; black, co-expression. (B) Top 30 hub genes in the PPI network. PPI, protein-protein interaction.

Identification of prognosis-associated genes in WT

KM analysis was performed to investigate the overall survival of patients with WT for the DERNAs (38 DElncRNAs, 18 DEmiRNAs and 99 DEmRNAs) in the ceRNA network. As a result, the expression levels of 2 DEmRNAs [zinc finger and BTB domain containing 4 (ZBTB4) and PHD finger protein 19 (PHF19)] (Figs. 5C and D), 5 DElncRNAs [maternally expressed 3 (MEG3), rhabdomyosarcoma 2 associated transcript (RMST), ZNF503 antisense RNA 1 (ZNF503-AS1), HNF1A antisense RNA 1 (HNF1A-AS1), deleted in lymphocytic leukemia 2 (DLEU2)] (Figs. S2A-C and 5E and F) and 4 DEmiRNAs (hsa-miR-132, hsa-miR-200a, hsa-miR-429 and hsa-miR-506) (Figs. 5G-I and S2D) were associated with the overall survival time of patients with WT (P<0.05). High expression of PHF19, DLEU2, ZNF503-AS1 and hsa-miR-506 was associated with low survival time. High expression of ZBTB4, MEG3, RMST, HNF1A-AS1, hsa-miR-132, hsa-miR-200a and hsa-miR-429 was associated with high survival time.

Figure 5.

Screening of the key RNAs in WT. Identification of key (A) upregulated competitive endogenous RNA-associated DERNAs and (B) downregulated DERNAs by combining expression and prognosis analyses. Expression and prognostic value of (C) ZBTB4, (D) PHF19, (E) DLEU2, (F) HNF1A-AS1, (G) hsa-miR-132, (H) hsa-miR-200a, and (I) hsa-miR-429 in *P<0.05, **P<0.01, ***P<0.001. WT. WT, Wilms' tumor; DE, differentially expressed; ZBTB4, zinc finger and BTB domain containing 4; PHF19, PHD finger protein 19; DLEU2, deleted in lymphocytic leukemia 2; HNF1A antisense RNA 1; miR, microRNA.

By combining the results of expression and survival analysis, DERNAs whose abnormal expression resulted in a trend for decreased survival times in WT patients were screened out

Among the 11 prognosis-related DERNAs, 2 RNAs (PHF19 and DLEU2) that were not only significantly upregulated in WT but also for which increased expression indicated a poor prognosis (Fig. 5A, D and E), 5 RNAs (ZBTB4, HNF1A-AS1, hsa-miR429, hsa-miR-132 and hsa-miR-200a) had low expression in tumor samples and the low expression of these indicated a poor prognosis (Fig. 5B, C and F-I). High expression levels of MEG3 and RMSET were associated with a good prognosis in patients with WT (Fig. S2A and B), while low expression levels of ZNF503-AS1 and hsa-miR-506 were associated with a good prognosis in patients with WT (Fig. S2C and D). A total of 7 RNAs (PHF19, DLEU2, ZBTB4, HNF1A-AS1, hsa-miR-429, hsa-miR-132 and hsa-miR-200a) which were associated with a poor prognosis in patients with WT were selected for further analysis. The association between the expression of these seven RNAs and histology classification and clinical staging were subsequently analyzed. For mRNAs, the results demonstrated that the expression of ZBTB4 was associated with clinical staging and histology classification (P<0.05; Fig. 6A and C). For lncRNAs, DLEU2 was also associated with both histology classification and clinical staging (P<0.05; Fig. 6B and D). For miRNAs, hsa-miR-132 was only significantly associated with histology classification (P<0.05; Fig. 6E).

Figure 6.

Analyses of prognosis-associated RNAs associated with histology classification and clinical staging in Wilms' tumor by χ2 test. (A and B) Clinical staging and (C-E) histology classification. ZBTB4, zinc finger and BTB domain containing 4; DLEU2, deleted in lymphocytic leukemia 2; HNF1A antisense RNA 1; miR, microRNA. FAWT, focal anaplastic Wilms' tumors; DAWT, diffuse anaplasia in Wilms' tumors.

GEO dataset verification

A total of seven aforementioned prognosis-associated RNAs (PHF19, DLEU2, ZBTB4, HNF1A-AS1, hsa-miR-429, hsa-miR-132 and hsa-miR-200a) were selected for validation. Consistent with the aforementioned results, downregulation of ZBTB4 (Fig. 7A and D) and upregulation of PHF19 (Fig. 7B and E) were confirmed in GSE66405 and GSE110696 datasets, respectively. The mean expression levels of HNF1A-AS1 were significantly lower in WT tissues compared with that in normal tissues in the GSE66405 dataset (Fig. 7C). DLEU2 was demonstrated to be significantly higher in WT tissues compared with that in normal tissues in the GSE110696 dataset (Fig. 7F). DLEU was not obtained from the GSE66405 dataset, as well as HNF1A-AS1 in the GSE110696 dataset; this may be as the two genes were not considered when designing probe sequences. For miRNA, in the GSE50505 dataset, the expression levels of hsa-miR-200a were significantly lower in WT tissues compared with that in normal tissue (Fig. 7G). In the GSE57370 dataset, hsa-miR-200a and hsa-miR-429 were significantly lower in WT tissues (Fig. 7H and I). No significant differences were found in hsa-miR-132 or hsa-miR-429 expression in the GSE50505 (Fig. S3A and B). hsa-miR-132 was not found to be DE in the GSE57370 dataset (Fig. S3C). This may be due to the small number of patients studied.

Figure 7.

Validation of the differentially expressed RNAs associated with survival. Box plots of (A) ZBTB4, (B) PHF19 and (C) HNF1A-AS1 expression levels in the GSE66405 validation dataset. Box plots of (D) ZBTB4, (E) PHF19, and (F) DLEU2 expression levels in the GSE110696 validation dataset. (G) Box plots of hsa-miR-200a miRNA expression levels in the GSE50505 validation dataset. Box plots of (H) hsa-miR-429 and (I) hsa-miR-200a miR expression levels in the GSE57370 validation dataset. *P<0.05, **P<0.01, ***P<0.001. miR/miRNA, microRNA; ZBTB4, zinc finger and BTB domain containing 4; PHF19, PHD finger protein 19; DLEU2, deleted in lymphocytic leukemia 2; HNF1A antisense RNA 1.

Flow diagram representing construction and analysis of the ceRNA regulatory network

A flow diagram representing the construction of the lncRNA-miRNA-mRNA regulatory network and the screening of prognostic RNAs in WT is shown in Fig. 8.

Figure 8.

Flow diagram of the construction and analysis of the competitive endogenous RNA regulatory network in Wilms' tumor. Pale yellow denotes RNA, which was associated with clinical staging and histology classification. Light orange represents RNA, which was associated with histology classification. TCGA, The Cancer Genome Atlas; lncRNA, long non-coding RNA; miRNA/miR, microRNA; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction; KM, Kaplan-Meier; ZBTB4, zinc finger and BTB domain containing 4; DLEU2, deleted in lymphocytic leukemia 2; DE, differentially expressed; cor, correlation.

Discussion

WT is the most common renal tumor in children and its carcinogenesis and progression is driven by multiple interacting mechanisms (44). The ceRNA hypothesis provides important clues and directions for the study of tumor pathogenesis, and provides a novel theoretical basis for the diagnosis and treatment of tumor (45). Studies have shown that lncRNA, miRNA and mRNA in the ceRNA network are in a certain equilibrium state, and that disease occurs once the balance is broken (46,47). The ceRNA hypothesis integrates the relationship between lncRNAs, mRNAs and miRNAs and can better explain the interaction among a variety of types of RNAs (22). In order to lay a useful foundation for studying the regulatory function of ceRNA in WT, a ceRNA network was constructed at the transcriptome-wide level and screened for prognosis-associated biomarkers for diagnosis and treatment purposes.

In the present study, a ceRNA regulatory network for WT was successfully constructed. The network was composed of 38 DElncRNAs, 18 DEmiRNAs and 99 DEmRNAs. Among these RNAs, combined with expression and the survival analysis, 7 DERNAs (PHF19, DLEU2, ZBTB4, HNF1A-AS1, hsa-miR-429, hsa-miR-132 and hsa-miR-200a) were significantly associated with prognosis. These RNAs may be potential biomarkers for predicting prognosis in WT. The expression of DLEU2 was associated with the histology classification and clinical staging, indicating that DLEU2 is an important lncRNA and the over expression of DLEU2 may play a role in the pathogenesis and progression of WT. Similar studies have found that the inactivation of DLEU2 can promote cell proliferation and tumor progression through functional loss of miR-15a/miR-16-1 (48–50). The present study demonstrated that DLEU2 was associated with two key DEmiRNAs (hsa-miR-21 and hsa-miR-506) and competed to regulate the mRNA expression in WT (Table SV). Therefore, an in-depth study of the mechanism of action of DLEU2 and its associated miRNAs maybe beneficial, as a potential treatment target for WT.

Previous studies have revealed that miRNAs have been extensively investigated in cancer (51) and may play important biological functions by regulating target genes in WT. For example, miR-100-5p and miR-130b-3p could be potential biomarkers for WT and the expression levels of these miRNAs in serum was stable over time with different serum storage conditions (52). The present study demonstrated that the expression level of hsa-miR-132 was lower in WT tissues compared with that in normal tissue and could be a significant indicator for poor prognosis in patients with WT. In addition, it was found that the expression level of hsa-miR-132 was associated with histology classification. hsa-miR-132, as a member of the miR-212/132 family, which is highly conserved in vertebrates, has been reported as a tumor-associated miRNA in a variety of cancer types, including liver, colorectal, pancreatic and ovarian cancer (53–56). From the aforementioned studies, it has been suggested that miR-132 may have a potential involvement in the occurrence and development of WT.

In addition, among the 99 target mRNAs in the ceRNA regulatory network, functional enrichment results indicated that these genes are involved in pathways associated with cancer development. Survival analyses demonstrated that two target mRNAs (ZBTB4 and PHF19) were significantly associated with the prognosis of patients with WT. Patients with low expression of ZBTB4 have poor prognoses, whereas low expression of PHF19 contributes to prolongation of survival time in patients with WT. It was also found that the expression level of ZBTB4 was associated with histology classification and clinical staging (P<0.05), which indicated that ZBTB4 may have an important role in the tumorigenesis of WT. ZBTB4, is a mammalian DNA-binding protein, and contains C2H2 zinc fingers and a POZ/BTB domain, and functions as a transcriptional repressor protein (57). Loss of ZBTB4 expression has been observed in several types of cancer, including breast cancer, prostate cancer and Ewing sarcoma, and was associated with shorter relapse-free patient survival (58–60). Roussel-Gervais et al (61) reported that ZBTB4 had the capacity for methyl-CpG-binding and had a conserved role in the preservation of genomic stability in 8 tumor types: breast-invasive carcinoma, prostate adenocarcinoma, uterine endometrioid carcinoma, kidney renal cell carcinoma, cervical squamous cell carcinoma, lung adenocarcinoma, stomach adenocarcinoma and head and neck squamous cell cancer. Furthermore, the loss of ZBTB4 induced transcriptional alterations indicative of aneuploidy and mitotic checkpoint deregulation (61). A recent study showed that high PHF19 expression was associated with a shorter survival time in patients with ovarian carcinoma, and that the silencing of PHF19 could reduce cell proliferation (62). Thus, further study on the mechanism of action of these two mRNAs is required, which may prove beneficial as possible targets for the treatment of WT.

The present study successfully identified numerous DElncRNAs, DEmiRNAs and DEmRNAs in WT. Moreover, a WT-specific lncRNA-associated ceRNA regulatory network was constructed, providing a potential target for the diagnosis, treatment and prognosis evaluation of WT. Importantly, novel candidate biomarkers consisting of two DElncRNAs (HNF1A-AS1 and DELU2), three DEmiRNAs (hsa-miR-429, hsa-miR-132 and hsa-miR-200a) and two DEmRNAs (ZBTB4 and PHF19) were identified. These biomarkers were significantly associated with overall survival time in patients with WT. Notably, ZBTB4 and DLEU2 were significantly associated both with clinical staging and histology classification. ZBTB4 and DLEU2 may therefore be considered as promising targets for therapy in WT. Moreover, further studies are required in order to validate the role of newly discovered genes in the mechanism of WT progression.