Abstract

Wilms tumor (WT) is the most common genitourinary renal tumor that typically occurs in children under 15 and is thought to be linked to somatic and germline mutations. However, the specific functional role of competing endogenous RNAs (ceRNAs) and their potential implications in WT remain unclear. In this study, we developed an lncRNA-mediated (long noncoding RNA-mediated) ceRNA network via the R packages for WT with expression data obtained from the tumor alterations relevant for genomics-driven therapy (TARGET) database. Unsupervised hierarchical clustering analysis revealed that the WT specimens could be clearly distinguished from healthy specimens with respect to the expression of disordered RNAs. A total of 1,607 differentially expressed (DE) lncRNAs, 116 DE microRNAs (DEmiRNAs), and 3,262 DE messenger RNAs (DEmRNAs) were identified as WT-specific RNAs, and a lncRNA-miRNA-mRNA ceRNA network with 159 DElncRNAs, 18 DEmiRNAs, 131 DEmRNAs, and 792 interactions was constructed. According to the clinical survival data, 12 DElncRNAs, 5 DEmRNAs, and 2 DEmiRNAs were selected from the ceRNA network that could significantly impact the overall survival of WT patients (). Functional enrichment analysis showed that the biological processes and pathways of DEmRNAs, such as cell cycle and virus infection, may be associated with WT. The present study constructed a dysregulated lncRNA-mediated ceRNA network in WT and discovered that lncRNA-mediated ceRNAs may serve as important regulators in WT development and progression. Survival-associated RNAs may serve as new potential biomarkers, suggesting that the constructed ceRNA network in WT might be important for determining optimal therapeutic strategies.

1. Introduction

Wilms tumor (WT) is the most common genitourinary renal tumor that typically occurs in early childhood [1], which is usually found in one or both kidneys and may metastasize to other important organs [2]. With new medicines and innovative immune therapies performed in WT [3], the survival rate of WT has remarkably increased [4]. However, as a malignant tumor mainly occurs in young children, a deep understanding of the occurrence and development of WT still requires more attention.

With the development of sequencing technology, increased dysregulation and mutations of DNA or RNA in human carcinoma including epigenetic alterations [5] were revealed, and the bioinformatics analysis showed the important roles they played in many tumors. However, the knowledge of the genetic underpinnings of WT was largely limited to somatic and germline mutations [6], such as aberrations of tumor protein 53 (TP53), the WT gene on the X chromosome (WTX), WT gene 1 (WT1), catenin beta 1 (CTNNB1), and the imprinted 11p15 region [79].

The competing endogenous RNA (ceRNA) hypothesis was first presented in 2011, which postulated that RNAs communicate with each other through shared microRNA (miRNA) response elements (MREs) [10]. Since a miRNA modulates the expression of several target messenger RNAs (mRNAs), each mRNA may be regulated by multiple miRNAs [11]. mRNAs, long noncoding RNAs (lncRNAs), and other RNA transcripts, which can act as endogenous miRNA sponges, in the ceRNA network can form a large-scale regulatory network. The dysregulation of RNAs involved in the ceRNA network has been previously studied in many cancers [12, 13]. The key RNAs are regarded as cancer biomarkers and contribute to tumor progression [1416].

Wilms tumor is relatively rare in absolute numbers (1 in 10,000 children). The scarcity of WT specimen makes it difficult to perform experimental research on WT. Moreover, the commercially available, previously widely used purported WT cell lines, including WT-CLS1 [17], G401 [18, 19], and SK-NEP-1 [19, 20], turned out to be misclassified. Thus, the studies about WT are limited by the paucity of available cell/tissue specimens, and bioinformatics analysis becomes an effective and economical strategy for WT study.

To the best of our knowledge, this study is the first time to construct a lncRNA-mediated WT-specific ceRNA network with both GEO and TARGET databases. In the present study, we identified WT-specific miRNAs and mRNAs convergence in both GEO and TARGET databases and provide novel insight into a better understanding of lncRNA-mediated ceRNA regulation in the tumorigenesis and progression of WT. We hope this work helps elucidate the lncRNA-miRNA-mRNA crosstalk in WT and provides further insight into its molecular mechanisms.

2. Materials and Methods

2.1. Data Acquisition and Processing

All datasets used in the present study were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) and the TARGET database (https://ocg.cancer.gov/programs/target). In the GSE66405 dataset, mRNA expression of 28 WT and 4 normal kidney samples were determined using the Agilent-039494 SurePrint G3 Human GE v2 8 × 60 K (GPL17077) platform. In the GSE57370 dataset, miRNA expression of 62 WT and 4 normal kidney samples were determined with the Agilent-031181 Unrestricted_Human_miRNA_V16.0_Microarray (GPL16770) platform. The expression quantification of mRNAs, lncRNAs, and miRNAs was obtained from the TARGET database. 125 WT patients with corresponding clinical data and 6 healthy controls were enrolled. The main characteristics of patients enrolled in the study are shown in Table 1, and the study flow diagram is shown in Figure 1.

2.2. Identification of DERNAs between WT Tissues and Healthy Tissues

The R 3.6.0 software (https://www.r-project.org/) was used to identify the differentially expressed RNAs between WT and normal samples. In the GSE66405 dataset, differentially expressed mRNAs were identified using the “limma” package with |log2 fold change (FC)| >1 and as cutoff criteria. Similarly, differentially expressed miRNAs were screened with |log2 FC| > 0.5 and in GSE57370. Additionally, the differentially expressed RNAs in the TARGET dataset were determined using the “edgeR” package with |log2 FC| > 2 and a false discovery rate (FDR) < 0.01 as the cutoff criteria.

2.3. Functional Enrichment Analysis

Functional enrichment analysis was performed for the WT-specific mRNAs involved in the ceRNA network to understand the underlying functional implications of these mRNAs in Wilms tumor. Gene ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were conducted using the R package “clusterProfiler” [21]. Disease ontology (DO) analysis was conducted using the R package “DOSE” [22]. was set as the statistical significance threshold criterion for enrichment analysis.

2.4. Construction of the Protein-Protein Interaction Network

mRNAs convergence in TARGET and GSE66405 were determined as WT-specific mRNAs. The direct and indirect correlations between WT-specific mRNAs were assessed from the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; Version 11.0; http://string-db.org/) database (with the minimum required interaction score of 0.9) [23]. The PPI network was reconstructed via Cytoscape software [24] (Version 3.8.0). Functionally related clusters were further identified with the molecular complex detection algorithm (MCODE; Version: 1.5.1), based on topology to locate densely connected regions. The mRNAs in the top cluster selected were taken as hub genes.

2.5. Expression Validation of Hub Genes in Oncomine Database

The Cutcliffe renal dataset and Yusenko renal dataset in Oncomine (https://www.Oncomine.org) were used to verify the differential expression of hub genes. In the Cutcliffe renal dataset, mRNA expression of 18 WT and 3 fetal kidneys were determined using human genome U133 A array in which ASPM and NUF2 were not tested. In the Yusenko renal dataset, mRNA expression of 4 WT and 2 fetal kidneys was determined using human genome U133 Plus 2.0 array. The log2 median-centered intensity value was visualized with GraphPad Prism (Version 5.0).

2.6. Construction of the ceRNA Network

The interactions between WT-specific lncRNAs and WT-specific miRNAs were predicted using the miRcode database (http://www.mircode.org/) [25]. Then, WT-specific miRNAs were modified with the StarBase v2.0 database (http://starbase. sysu.edu.cn/), and their target mRNAs were retrieved according to three different databases (TargetScan [26], miRTarBase [27], and miRDB [28]). Only mRNAs predicted by all three databases were considered as candidate DEmRNAs for further study. Then, the lncRNA-miRNA-mRNA ceRNA network was constructed and visualized using the Cytoscape software (Version 3.8.0) [29].

2.7. Survival Analysis

WT patients were divided into a specific RNA high-expressing group and a specific RNA low-expressing group, with the median expression value as the cutoff value. Survival analysis between the specific RNA high-expressing and specific RNA low-expressing groups was evaluated using the Kaplan–Meier survival curve and log-rank test analysis. The survival-associated DERNAs () were identified as prognosis-associated RNAs.

3. Results

3.1. Identification of WT-Specific mRNAs

A total of differentially expressed 3,262 differentially expressed mRNAs (1,756 upregulated and 1,506 downregulated) were altered significantly in the TARGET-WT dataset. In addition, 724 (170 upregulated and 554 downregulated) differentially expressed mRNAs in GSE66405 were determined with the aforementioned cutoff thresholds. The distribution of all differentially expressed mRNAs is depicted in the volcano maps shown in Figures 2(a) and 2(b). The 116 upregulated and 334 downregulated mRNAs convergence in both databases were identified as WT-specific mRNAs and presented in Venn diagrams (Figures 2(c) and 2(d)).

3.2. Functional Enrichment Analysis

To outline the potential function of the WT-specific mRNAs, functional enrichment analysis was performed with “ClusterProfiler” and “DOSE.” Functional enrichment analysis revealed that a total of 394 GO terms, including 277 biological process terms (GO.BP), 31 cellular component terms (GO.CC), and 86 molecular function terms (GO.MF), were enriched. The top 10 predominant BP terms, CC terms, and MF terms in GO functional enrichment analysis are shown in Figure 3(a) which were concerned with organic acid carboxylic and transport. These results suggested that tumorigenesis and the development of WT may be related to the dysfunction of renal ion transport.

14 significantly enriched KEGG pathways were identified and are shown in Figure 3(b). KEGG results showed that the “PPAR signaling pathway” was the most concerned pathway in WT. PPARs belong to the family of ligand-activated nuclear receptors. Specific PPAR ligands have been proposed as potential therapies for a variety of diseases such as metabolic syndrome and cancer [30]. Moreover, DO analysis (Figure 3(c)) showed that WT was not only a urinary system disease but also concerned with “nephrocalcinosis” and calcium/mineral metabolism disease, which were verified in GO and KEGG analyses. Moreover, obesity and overnutrition were also involved, which may originate partly in lifestyle, in particular via markedly reduced levels of physical activity after diagnosis.

3.3. Protein-Protein Interaction Network and Hub Genes

To better understand the interplay among the WT-specific mRNAs, a PPI network with 83 nodes and 276 edges was reconstructed by Cytoscape (Figure 4(a)). The degree distribution of the nodes in this PPI network was analyzed, and the mRNAs with top 30° are shown in Figure 4(b). BUB1 mitotic checkpoint serine (BUB1) was identified to have the highest degree in the overall network (degree = 25).

Based on the topology to locate densely connected regions, some genes are notably concentrated. With the help of MCODE, a hub gene cluster with 22 nodes and 203 edges was selected. These 22 hub genes were all upregulated and with a high degree in the overall network (Figure 4(b)). NUF2 (NUF2) was identified to be the seed gene of the cluster (degree = 14). Most of the hub genes were predominantly involved in the cell cycle (BUB1, mitotic checkpoint serine/threonine kinase B (BUB1B), cyclin B2 (CCNB2), and pituitary tumor-transforming 1 (PTTG1)) and oocyte meiosis pathways (BUB1, CCNB2, and PTTG1). Expression validation of the hub genes was performed with the Yusenko renal dataset (Figure 4(c)) and Cutcliffe renal dataset (Figure 4(d)) in Oncomine. Except for CENPF (centromere protein F) was unexpectedly downregulated, all the other hub genes were upregulated in WT tissues in both Oncomine datasets as expected.

3.4. Identification of WT-Specific miRNAs and lncRNAs

A total of 116 differentially expressed miRNAs (including 77 upregulated and 39 downregulated) were identified as altered significantly in the TARGET-WT dataset and 80 (29 upregulated and 51 downregulated) differentially expressed miRNAs in GSE57370 were determined with aforementioned cutoff thresholds. The distribution of all differentially expressed miRNAs is depicted in the volcano maps shown in Figures 5(a) and 5(b). The 19 miRNAs convergence in both databases were identified as WT-specific miRNAs and are presented in Venn diagrams (Figures 5(c) and 5(d)).

3.5. ceRNA Network Construction in WT

1607 differentially expressed lncRNAs (including 851 upregulated and 756 downregulated) identified in TARGET database with the aforementioned cutoff thresholds were taken as WT-specific lncRNAs. Then, we screened the miRcode database with the WT-specific lncRNAs and obtained 689 interactions between 90 WT-specific lncRNAs and five WT-specific miRNAs for further analysis. The five miRNAs recognized above were mapped into the TargetScan, miRTarBase, and miRDB databases to identify their target mRNAs. The 150 mRNAs predicted by all three databases were considered target mRNAs candidates, and only 27 of them which were differentially expressed in WT tissues were identified as DEmRNAs for further analysis. Finally, 90 DElncRNAs, 5 DEmiRNAs, and 27 DEmRNAs enrolled the ceRNA network with 177 interactions (145 DElncRNA–DEmiRNA interactions and 32 DEmiRNA–DEmRNA interactions), constructed, and visualized with Cytoscape (Figure 6 and Table 2).

The rhombuses represent DElncRNAs, the squares represent DEmiRNAs, and the rounds represent DEmRNAs. Red represents upregulation and green represents downregulation.

3.6. Prognosis-Associated RNAs in the WT ceRNA Network

To identify the prognosis-associated RNAs, the Kaplan–Meier survival curve and log-rank test analysis of all the RNAs involved in the ceRNA network were performed. According to the results, 7 DElncRNAs and 2 DEmiRNAs were identified as prognosis-associated RNAs (, Figure 7).

Among these RNAs, the high expression of three risky lncRNAs including AL445228.2, LMO7-AS1 LMO7 antisense RNA (1), and DLEU2 (deleted in lymphocytic leukemia 2) was associated with shorter overall survival of WT patients (Figures 7(a) and 7(c)). The remaining four DElncRNAs (MEG3, maternally expressed gene 3; RMST, rhabdomyosarcoma 2-associated transcript; HNF1A-AS1, HNF1A antisense RNA 1; and PVT1, and plasmacytoma variant translocation 1) and two DEmiRNAs (hsa-mir-429 and hsa-mir-200a) appeared to be protective (Figures 7(d) and 7(i)).

In the figure, Kaplan–Meier survival curves for lncRNA (a) AL445228.2, (b) LMO7-AS1, (c) DLEU2, (d) MEG3, (e) RMST, (f) HNF1A-AS1, (g) PVT1 and miRNAs, (h) hsa-mir-429, and (i) hsa-mir-200a are shown.

4. Discussion

Several genes, including CTNNB1, WTX, WT1, and TP53, have been reported to be involved in the tumorigenesis and progression of WT [79]. However, the regulatory function of dysregulated RNAs (including lncRNAs, miRNAs, and mRNA) in WT remains elusive. Increasing evidence has revealed that dysregulated RNAs play important roles in many cancers [3133]. miRNA profile changes in WT have been used as predictors of chemo responsiveness in WT blastema [34]. LINC00473 mediates the pathogenesis of WT by antagonizing the tumor suppressor hsa-mir-195 [35]. The elevated expression of FOXM1 (forkhead box protein M1) has been reported as a new prognostic biomarker that is associated with histology and prognosis of WT [36]. Therefore, identification of key RNAs is vital for understanding the pathogenesis and abnormal biological behavior of WT which may help to identify novel therapeutic targets.

First, we analyzed the microarray data and RNA-sequencing data from the GEO and TARGET databases. Then, 450 WT-specific mRNAs convergence in both databases was identified for bioinformatics analysis. Identify their potential associated cellular signaling pathways and functions. GO analysis revealed that the mRNAs were involved in the “organic anion transport” and “carboxylic acid biosynthetic process” biological processes terms. Furthermore, the mRNAs were demonstrated to serve a role in transmembrane transporter activity at the molecular function level. KEGG pathway analysis revealed that 12 mRNA were enriched in the “PPAR signaling pathway,” 10 were enriched in the “mineral absorption,” and 9 were enriched in “protein digestion and absorption.”

According to the MCODE analysis, a hub gene cluster with all the genes upregulated was identified. In line with the results of functional enrichment analysis of all WT-specific mRNAs, the hub genes were mostly enriched in the GO terms concerned with the “cell cycle” too, such as “nuclear division,” “organelle fission,” “mitotic nuclear division,” and “chromosome segregation.” These results suggested that drugs targeting the cell cycle may be effective in the treatment of WT. Expression validation of the hub genes was performed with Oncomine. The validation study showed that CENPF was unexpectedly downregulated, which may be due to the paucity of specimens in the Cutcliffe renal dataset and Yusenko renal dataset. Moreover, other hub genes were still upregulated in both datasets. Therefore, it can be hypothesized that abnormal regulation of the hub genes may be generally existed in WT and contribute to the tumorigenesis and progression of WT.

Subsequently, to comprehensively understand how dysregulated RNAs work in WT, we constructed a ceRNA with 90 DElncRNAs, 5 DEmiRNAs, 27 DEmRNAs, and 177 interactions. Nine of these RNAs, including AL445228.2, LMO7-AS1, DLEU2, MEG3, RMST, HNF1A-AS1, PVT1, hsa-mir-429, and hsa-mir-200a, were considered to be prognostic biomarkers as they were significantly associated with overall survival in WT patients. In line with the GO enrichment results, most of these RNAs were associated with tumor progression via the cell cycle.

Risky marker DLEU2 was reported to promote cancer cell proliferation and invasion in pancreatic cancer [37]. Protective lncRNA MEG3 was regarded as a novel tumor suppressor by inhibiting tumor cell proliferation in many cancers [3840]. Protective lncRNA RMST (rhabdomyosarcoma 2-associated transcript) could enhance cell apoptosis in triple-negative breast cancer (TNBC) [41]. HNF1A-AS1 was considered to function as a regulator of cell proliferation and migration in esophageal adenocarcinoma [42] and lung cancer [43]. Still, there is WT-protective lncRNA which was a risky factor in other cancers. For example, PVT1 was reported to promote tumorigenesis in nonsmall cell lung cancer [44] and multidrug resistance in gastric cancer [45].

It is well-known that miRNA-regulated pathways are indispensable in studies of tumorigenesis [46, 47]. miRNAs are involved in the occurrence, development, incursions, and metastasis of tumors [48, 49]. In this study, among the five DEmiRNAs involved in the ceRNA network, two miRNAs (hsa-miR-200a and hsa-miR-429) were associated with the poor prognosis of WT patients. Both of them belong to the miR-200 family, which has been shown to be closely associated with carcinogenesis and progression, and potentially be important for the diagnosis and treatment of cancer [50]. hsa-mir-200a was reported to inhibit cell growth, migration, and invasion in meningioma [51] and nasopharyngeal carcinoma [52] and determines prognosis in colorectal cancer patients [53] and ovarian tumorigenesis [54]. As the has-mir-200 family is important for maintaining the epithelial phenotype [55], the dysregulation of hsa-mir-200a and hsa-mir-429 in WT patients may be associated with migration and invasion via the epithelial-mesenchymal transition.

5. Conclusions

In the present study, we successfully constructed a lncRNA-mediated ceRNA network with WT-specific lncRNAs and miRNAs/mRNAs convergence in both GEO and TARGET databases. The ceRNA network may provide novel insight into a better understanding of ceRNA regulation in WT. Furthermore, the ceRNA network and the key RNAs identified in this study may provide potential biomarkers for diagnosis and prognosis and improve clinical outcomes for children with WT.

Data Availability

The datasets generated and analyzed during the current study are available in the GEO database (https://www.ncbi.nlm.nih.gov/geo/), the Oncomine (https://www.Oncomine.org), and the TARGET database (https://ocg.cancer.gov/programs/target).

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Authors’ Contributions

Biao An and Yuan Hu contributed equally to this work.

Acknowledgments

This study was funded by grants from the National Natural Science Foundation of China (81902662 and 81821002), the Natural Science Foundation of Sichuan Province (2022NSFSC1531), Sichuan Science and Technology Program (2021YJ0011 and 2018JY0609), the Basic Research Project of Nursing Department of West China Second Hospital (HLBKJ202126), and Foundation of Development and Related Diseases of Women and Children Key Laboratory of Sichuan Province (2022003).