Integrated Analysis of the lncRNA-Associated ceRNA Network in Wilms Tumor via TARGET and GEO Databases

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 (P < 0.05). 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.


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 [7][8][9]. e 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 largescale regulatory network. e dysregulation of RNAs involved in the ceRNA network has been previously studied in many cancers [12,13]. e key RNAs are regarded as cancer biomarkers and contribute to tumor progression [14][15][16].
Wilms tumor is relatively rare in absolute numbers (1 in 10,000 children). e 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.
us, 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.

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. e 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. e main characteristics of patients enrolled in the study are shown in Table 1, and the study flow diagram is shown in Figure 1.

Identification of DERNAs between WT Tissues and Healthy Tissues.
e R 3.6.0 software (https://www.rproject.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 |log 2 fold change (FC)| >1 and P < 0.05 as cutoff criteria. Similarly, differentially expressed miRNAs were screened with |log 2 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 |log 2 FC| > 2 and a false discovery rate (FDR) < 0.01 as the cutoff criteria.

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]. P < 0.05 was set as the statistical significance threshold criterion for enrichment analysis.

Construction of the Protein-Protein Interaction Network.
mRNAs convergence in TARGET and GSE66405 were determined as WT-specific mRNAs. e 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]. e 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. e mRNAs in the top cluster selected were taken as hub genes.

Expression Validation of Hub Genes in Oncomine
Database.
e 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

Survival
Analysis. WT patients were divided into a specific RNA high-expressing group and a specific RNA lowexpressing 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. e survival-associated DERNAs (P < 0.05) were identified as prognosis-associated RNAs.

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. e 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. ese 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.

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)). e 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. ese 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.

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. e distribution of all differentially expressed miRNAs is depicted in the volcano maps shown in

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. en, 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. e five miRNAs recognized above were mapped into the TargetScan, miRTarBase, and miRDB databases to identify their target mRNAs. e 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). e rhombuses represent DElncRNAs, the squares represent DEmiRNAs, and the rounds represent DEmRNAs. Red represents upregulation and green represents downregulation.

Discussion
Several genes, including CTNNB1, WTX, WT1, and TP53, have been reported to be involved in the tumorigenesis and progression of WT [7][8][9]. 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 [31][32][33]. 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]. e 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]. erefore, 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. en, 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."

Genetics Research
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." ese 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.
e 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. erefore, 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.
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 WTprotective 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 WTpatients. 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.

Conclusions
In the present study, we successfully constructed a lncRNAmediated ceRNA network with WT-specific lncRNAs and miRNAs/mRNAs convergence in both GEO and TARGET databases. e 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.

Conflicts of Interest
e authors declare that there are no conflicts of interest.