1. Introduction
Wilms tumor (WT) is the most common genitourinary renal tumor that typically occurs in early childhood [Reference Breslow, Olshan, Beckwith and Green1], which is usually found in one or both kidneys and may metastasize to other important organs [Reference Vujanić and Sandstedt2]. With new medicines and innovative immune therapies performed in WT [Reference Hingorani, Zhang and Zhang3], the survival rate of WT has remarkably increased [Reference Dome, Graf and Geller4]. 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 [Reference Hol, Kuiper and van Dijk5] 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 [Reference Zhu, Jia and Wu6], 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 [Reference Little, Hanks and King-Underwood7–Reference Andrade, Cardoso and Ferman9].
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) [Reference Salmena, Poliseno, Tay, Kats and Pandolfi10]. Since a miRNA modulates the expression of several target messenger RNAs (mRNAs), each mRNA may be regulated by multiple miRNAs [Reference Thomas, Lieberman and Lal11]. 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 [Reference Zan and Li12, Reference Fan and Liu13]. The key RNAs are regarded as cancer biomarkers and contribute to tumor progression [Reference Cheng, Geng and Wang14–Reference Yan, Lu and Mao16].
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 [Reference Stroup, Yeu and Budhipramono17], G401 [Reference Garvin, Re, Tarnowski, Hazen-Martin and Sens18, Reference Pritchard-Jones and Perotti19], and SK-NEP-1 [Reference Pritchard-Jones and Perotti19, Reference Smith, Morton, Phelps, Girtman, Neale and Houghton20], 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.
TABLE 1: The main characteristics of 125 Wilms tumor patients.
FHWT, favorable histology; DAWT, diffuse anaplasia.
FIGURE 1: Flow diagram of the study.
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 P < 0.05 as cutoff criteria. Similarly, differentially expressed miRNAs were screened with |log2 FC| > 0.5 and P < 0.05 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” [Reference Yu, Wang, Han and He21]. Disease ontology (DO) analysis was conducted using the R package “DOSE” [Reference Yu, Wang, Yan and He22]. P < 0.05 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) [Reference Szklarczyk, Gable and Lyon23]. The PPI network was reconstructed via Cytoscape software [Reference Kohl, Wiese and Warscheid24] (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/) [Reference Jeggari, Marks and Larsson25]. 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 [Reference Agarwal, Bell, Nam and Bartel26], miRTarBase [Reference Chou, Chang and Shrestha27], and miRDB [Reference Wong and Wang28]). 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) [Reference Shannon, Markiel and Ozier29].
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 (P < 0.05) 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)).
FIGURE 2: Identification of WT-specific mRNAs. Differentially expressed mRNA in TARGET (a) and GSE66405 (b) with the aforementioned cutoff thresholds is depicted in the volcano map. The red/green points represent the upregulation/downregulation of mRNAs, and the black points represent RNAs with no significant difference. The 116 upregulated (c) and 334 downregulated (d) WT-specific mRNAs convergence in both databases is presented in Venn diagrams
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.
FIGURE 3: Functional enrichment analysis for WT-specific mRNAs. (a) The top 10 predominant GO.BP, GO.CC, and GO.MF terms are shown in the dot plot. (b) All the 14 significantly enriched KEGG pathways are shown in the dot plot. (c) The DO enrichment analysis showed the top 20 WT-specific mRNAs associated with diseases.
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 [Reference Wagner and Wagner30]. 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).
FIGURE 4: Identification and expression validation of hub genes. (a) The PPI network of WT-specific mRNAs shown by Cytoscape. (b) The degree distribution of the PPI network of the top 30 WT-specific mRNAs. (c) Expression validation of the hub genes performed with the Yusenko renal dataset. Mean ± SEM. (d) Expression validation of the hub genes performed with the Cutcliffe renal dataset. Mean ± SEM.
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)).
FIGURE 5: Identification of WT-specific miRNAs. Differentially expressed miRNA in TARGET (a) and GSE57370 (b) with the aforementioned cutoff thresholds was depicted in the volcano map. The red/green points represent the upregulation/downregulation of mRNAs, and the black points represent RNAs with no significant difference. the 5 upregulated (c) and 14 downregulated (d) WT-specific miRNAs convergence in both databases is presented in Venn diagrams.
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).
FIGURE 6: Construction of the ceRNA network.
TABLE 2: The DERNAs enrolled in the ceRNA network.
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 (P < 0.05, Figure 7).
FIGURE 7: Kaplan–Meier survival curves for the prognosis-associated DERNAs. Kaplan–Meier survival curves for lncRNA (a) AL445228.2, (b) LMO7-AS1, (c) DLEU2, (d) MEG3, (e) RMST, (f) HNF1A-AS1 and (g) PVT1) and miRNAs (h) hsa-mir-429 and (i) hsa-mir-200a.
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 [Reference Little, Hanks and King-Underwood7–Reference Andrade, Cardoso and Ferman9]. 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 [Reference Huang, Qian, Zhu, Huang, Luo and Qing31–Reference Wang, Zhou and Wang33]. miRNA profile changes in WT have been used as predictors of chemo responsiveness in WT blastema [Reference Watson, Bryan and Williams34]. LINC00473 mediates the pathogenesis of WT by antagonizing the tumor suppressor hsa-mir-195 [Reference Zhu, Fu and Zhang35]. 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 [Reference Apelt, Hubertus and Mayr36]. 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 [Reference Xu, Gong and Zi37]. Protective lncRNA MEG3 was regarded as a novel tumor suppressor by inhibiting tumor cell proliferation in many cancers [Reference Zhang, Yao, Wang, Yu, Liu and Lin38–Reference Guo, Qian, Yan, Li and Huang40]. Protective lncRNA RMST (rhabdomyosarcoma 2-associated transcript) could enhance cell apoptosis in triple-negative breast cancer (TNBC) [Reference Wang, Liu and Wu41]. HNF1A-AS1 was considered to function as a regulator of cell proliferation and migration in esophageal adenocarcinoma [Reference Yang, Song and Cheng42] and lung cancer [Reference Wu, Liu, Shi, Yao, Yang and Song43]. 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 [Reference Yang, Zang, Zhong, Li, Zhao and Feng44] and multidrug resistance in gastric cancer [Reference Zhang, Bu, Liu, Zhang and Li45].
It is well-known that miRNA-regulated pathways are indispensable in studies of tumorigenesis [Reference Niu, Sun, Guo, Niu and Liu46, Reference Song and Meltzer47]. miRNAs are involved in the occurrence, development, incursions, and metastasis of tumors [Reference Huang, Gumireddy and Schrier48, Reference Evans-Knowell, LaRue and Findlay49]. 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 [Reference Humphries and Yang50]. hsa-mir-200a was reported to inhibit cell growth, migration, and invasion in meningioma [Reference Senol, Schaaij-Visser and Erkan51] and nasopharyngeal carcinoma [Reference Williams, Renthal, Condon, Gerard and Mendelson52] and determines prognosis in colorectal cancer patients [Reference Pichler, Ress and Winter53] and ovarian tumorigenesis [Reference Mateescu, Batista and Cardon54]. As the has-mir-200 family is important for maintaining the epithelial phenotype [Reference Lin, Wang, Fewell, Cameron, Yin and Flemington55], 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).
Open access
