A Survival-Related Competitive Endogenous RNA Network of Prognostic lncRNAs, miRNAs, and mRNAs in Wilms Tumor

Wilms tumor (WT) commonly occurs in infants and children. We evaluated clinical factors and the expression of multiple RNAs in WT samples in the TARGET database. Eight long non-coding RNAs (lncRNAs; AC079310.1, MYCNOS, LINC00271, AL445228.3, Z84485.1, AC091180.5, AP002518.2, and AC007879.3), two microRNAs (miRNAs; hsa-mir-152 andhsa-mir-181a), and nine messenger RNAs (mRNAs; TCTEX1D4, RNF133, VRK1, CCNE1, HEY1, C10orf71, SPRY1, SPAG11A, and MAGEB18) were screened from differentially expressed RNAs and used to construct predictive survival models. These models showed good prognostic ability and were highly correlated with tumor stage and histological classification. Additionally, survival-related ceRNA network was constructed using 35 RNAs (15 lncRNAs, eight miRNAs, and 12 mRNAs). KEGG pathway analysis suggested the “Wnt signaling pathway” and “Cellular senescence” as the main pathways. In conclusion, we established a multinomial predictive survival model and a survival-related ceRNA network, which provide new potential biomarkers that may improve the prognosis and treatment of WT patients.


INTRODUCTION
Wilms tumor (WT), a renal malignancy originating from the metanephric blastema, is widespread in infants and children (1). WT accounts for 7% of all pediatric malignancies and occurs in one out of every 10,000 children (2). Fortunately, with the development of treatments, the survival rate of children with WT has increased by nearly 60% (3). The International Society of Pediatric Oncology stated that the combination of nephrectomy and chemotherapy significantly improved overall survival (OS) by more than 90% (4). However, 25% of children still have a poor prognosis based on the tumor stage (5). Therefore, clarifying the cellular process involved in WT development and providing prognostic biomarkers are important steps to improving the survival of patients.
In recent years, several studies have suggested non-coding RNAs as key molecules involved in tumorigenesis and tumor progression (6). microRNAs (miRNAs) are non-coding RNAs composed of 18-25 nucleotides that can negatively regulate gene expression (7). By contrast, long non-coding RNAs (lncRNAs) are more than 200 nucleotides in length and they adjust the biological behavior of tumors through competing endogenous RNAs (ceRNAs) (8). It is suggested that there is a complex regulatory network among lncRNAs, miRNAs, and messenger RNAs (mRNAs). lncRNAs competitively inhibit the function of miRNAs by acting as a sponge, thus indirectly disrupting mRNA expression and ultimately affecting gene expression (9). In general, the instability of a ceRNA network may induce tumorigenesis (10,11). Several studies have suggested that a ceRNA network can be used as a prognostic biomarker for WT, but did not actually describe the effect of a ceRNA network in WT (12,13). Therefore, the establishment of a ceRNA network related to the survival of WT patients is of great significance for judging the prognosis of patients. By understanding the role of various RNAs in tumorigenesis, we can find potential targets to improve the prognosis of WT.
In this study, we explored the ability of multiple RNAs to prognosticate WT as a whole and established RNA models that can be used to predict survival. In addition, we established a survival-related ceRNA network and tried to understand its molecular mechanism through functional enrichment. Finally, we provided new potential biological biomarkers to improve the prognosis and treatment of WT.

Data Collection and Processing
Clinical and RNA sequencing data of all patients were obtained from the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) database, which can be downloaded from The Cancer Genome Atlas (TCGA) portal (https://portal.gdc.cancer.gov/; Data Release 25.0; release time: July 22, 2020). This study met the requirements for using TCGA database and did not require the approval of an ethics committee. We selected only the sequencing data from primary solid tumors for analysis. The mRNA and lncRNA sequencing data included 125 primary WT samples and six normal samples. The miRNA sequencing data included 127 primary WT samples and six normal samples. The clinical data of the 128 patients included in this study are shown in Table 1.

Identification of Differentially Expressed RNAs
The edgeR package in the R 4.0.2 software was used to analyze differentially expressed lncRNAs (DElncRNAs), differentially expressed miRNAs (DEmiRNAs), and differentially expressed mRNAs (DEmRNAs) between the WT and normal samples. The cutoff value for differentially expressed RNAs (DERNAs) was | log2 fold-change (FC)| >1 and false discovery rate (FDR) <0.05. Visualization of DERNAs was performed using ggplot2 package in R software.
Survival Analysis dentification of survival-associated RNAs was performed via univariate Cox regression analyses in the R software. Significantly correlated survival-associated RNAs (P < 0.003) were chosen for multivariate Cox regression analysis and to establish predictive survival models. Prognosis index (PI) = (expression gene1 × b gene1 ) + (expression gene2 × b gene2 ) +… + (expression genen × b genen ). Patients were divided into two groups based on the median PI. Survival prognosis of the two groups was compared by Kaplan-Meier analysis. The receiver operating characteristic (ROC) curve for evaluating the predictive ability of the model was depicted through R software.

Protein-Protein Interaction Network Construction
The Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) database was used to obtain PPI data of the significant survival-associated mRNAs (P < 0.003). Establishment of the PPI network was performed using the Cytoscape software.

Functional Enrichment Analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of mRNAs in the ceRNA network were conducted using KOBAS 3.0 (http:// kobas.cbi.pku.edu.cn/kobas3/). Visualization of the enrichment analyses was conducted using the R software.

Statistical Analysis
The correlation of RNAs with clinical characteristics was analyzed by rank sum test. Differences between survival curves were

Identification of Survival-Associated RNAs in WT
The relationship between DERNAs and survival was assessed by univariate Cox regression analysis. RNAs with P <0.05 were selected as survival-associated RNAs, yielding a total of 696 survival-associated RNAs (199 lncRNAs, 17 miRNAs, and 480 mRNAs). The top 15 survival-associated RNAs are shown in Figure 1. The gene networks of strongly correlated survivalassociated mRNAs (P < 0.003) were constructed by STRING ( Figure 2). The hub genes, including XAB2, SNRPA, PRPF19, and TP53, are shown in the PPI network.

Establishment of Predictive Survival Models
RNAs with strong correlation were selected by univariate Cox regression analysis (P < 0.003), and then multivariate Cox regression analysis was used to analyze survival-associated RNAs with strong correlation. Finally, a total of eight lncRNAs (AC079310 .1 PI values were calculated for each patient and divided into two groups. We found that among the three groups of RNAs, the low-risk group had better survival as determined by Kaplan-Meier analysis ( Figures 3A-C). The ability of the models to predict 3-year survival was evaluated by drawing the ROC curve. The areas under the curves (AUCs) of the three groups were 0.818, 0.701, and 0.848, respectively ( Figures 3D-F). These results suggest that the three groups of modules have great potential for predicting the clinical prognosis of WT. Figure 4 shows the risk scores, survival status, and RNA expression profiles in each group.
The correlation of these RNAs with other clinical characteristics was assessed by rank sum test. Clinical characteristics included age (<5/≥5), gender (male/female), race (white/non-white), tumor stage (I-II/III-IV), and pathological classification (FHWT/DAWT). We found that the RNAs in the models were significantly correlated with tumor stage and histological classification ( Figure 5). This implies that these RNAs can be used as potential indicators to judge the degree of WT tumor development. Next, we identified the relationship between clinical characteristics and OS. Multivariate Cox regression analysis suggested that tumor stage and risk level directly affected tumor prognosis ( Table 2).

Construction of a Survival-Related ceRNA Network in WT
Based on the survival-associated RNAs, 39 pairs of lncRNA-miRNA and 13 pairs of miRNA-mRNAs were detected from the databases. Subsequently, a ceRNA network containing 15 lncRNAs, eight miRNAs, and 12 mRNAs was constructed ( Figure 6). Many of these RNAs have been extensively studied as cancer-related molecules, such as miR-181a, CCNE1, and WIF1. Next, we further studied the molecular function of mRNAs in the ceRNA network. A total of 99 functional enrichment terms (68 biological processes, 13 cellular components, and 18 molecular functions) from the GO analysis and 15 KEGG pathways were observed. The biological processes were mainly enriched in "cell surface receptor signaling pathway involved in cell-cell signaling"; the cellular components were mainly enriched in "nucleus"; and the molecular function was mainly enriched in "binding". The KEGG analysis suggested that the "Wnt signaling pathway" and "Cellular senescence" were the main pathways ( Figure 7).

DISCUSSION
WT is a common pediatric tumor and poor prognosis is the main factor affecting long-term survival of patients (14). However, the exact molecular mechanism of WT is still unclear. The advancement in sequencing technology and the proposed ceRNA hypothesis provide a new perspective for the study of tumorigenesis and tumor progression (9). In the present study, we screened out survival-associated RNAs in the TARGET database and studied molecular events associated with WT. In particular, we proposed a novel survival-related ceRNA network that provides a new dimension for predicting the prognosis of WT patients.   Among all the DERNAs, we screened a total of 696 survivalrelated RNAs, including 199 lncRNAs, 14 miRNAs, and 480 mRNAs. mRNAs are key to the realization of molecular functions. PPI network analysis was carried out to identify the hub genes. We found that XAB2, SNRPA, PRPF19, and TP53 had a higher degree of connection among all genes. It is wellknown that p53 is a tumor suppressor gene, and some studies have confirmed that TP53 gene mutation worsens the prognosis of WT (15). The other hub genes also affect the progress of other tumors (16,17), suggesting that these hub genes may also affect the development of WT.
Most recent studies on WT have focused on the search for a single RNA that may influence prognosis. Zong et al. (18) found that miR-30d expression in WT tissues was significantly lower than that in normal tissues. Further animal experiments showed that miR-30d mimic significantly inhibited tumor growth. In addition, overexpression of SOX4 could reverse the effect of mir-30d on WT cells. Bao et al. (19) found that the expression of mirna-203a was low in WT tumor tissues, and that its expression level was directly related to the prognosis of WT. Knockout of miR-203a significantly enhanced the invasiveness of WT cells. Luciferase analysis confirmed that miR-203a targeted JAG1 to exert biological functions. Wang et al. (20) found that overexpression of miR-613 inhibited the G0/G1 phase transition of WT cells and hindered the expansion ability of WT cells. There is also a subset of studies evaluating the relationship between single RNA and OS. Tang et al. (12) predicted 32 DERNAs to be possibly associated with OS in WT patients through Kaplan-Meier analysis, and Zhang et al. (13) predicted 61 DERNAs to be possibly associated with OS by Kaplan-Meier analysis.
However, the development of WT is influenced by multifactorial causes, indicating that we should analyze prognostic markers at a broader level. We established a multinomial predictive survival model based on survival-associated RNAs by multivariate Cox analysis, which could suggest prognostic survival of WT patients. experiment only selected the sequencing data of primary tumor tissues rather than the sequencing data of all tumors including metastatic tumors. Further, we divided the patients into high-risk and low-risk groups according to PI value, and Kaplan-Meier analysis showed that there was a significant difference in survival between the two groups. Visualizing the relationship among survival time, PI score, and grouping intuitively showed the poor survival of high-risk patients. These results indicate that these three RNA models have good specificity and sensitivity in predicting the survival of WT patients. Notably, we found that these predictive RNAs were significantly associated with tumor stage and histological classification, which means that PI may have great value in the prognosis of WT patients. Multivariate Cox regression analysis also suggested that risk level can predict the prognosis of WT patients.
According to the establishment of the predictive survival models, we found some biomarkers that may have important clinical significance. lncRNAs participate in multiple cellular activities and have been shown to regulate tumorigenesis and metastasis (22,23). miRNAs, as a bridge between lncRNAs and mRNAs, participate in tumor development. MYCNOS promotes the invasion and metastasis of neuroblastoma and rhabdomyosarcoma by regulating MYCN protein (24,25) and may participate in the development of WT (26). miR-152 targeting DNMT1 inhibits the development of endometrial cancer (27), glioblastoma (28), and lymphomas (29). miR-181a mediates the Wnt/b-catenin pathway to accelerate the progression of colorectal cancer (30), oral squamous cell carcinoma (31), and acute lymphoblastic leukemia (32). The expression of VRK1 directly affects the proliferation of breast cancer and liver cancer (33). CircAGFG1 promotes triplenegative breast cancer progression through the circAGFG1/ miR-195/CCNE1 pathway (34). Precancerous lesions in squamous cell carcinoma depend on upregulation of the NOTCH4-HEY1 pathway (35). The HGF-mediated c-Met/ FRA1/HEY1 cascade may be the key to inducing the transition from cirrhosis to hepatocellular carcinoma (36). miRNA-21 promotes proliferation of human glioma cells through the PI3K/AKT/SPRY1 pathway (37). The expression of MAGEB18 affects cell proliferation and apoptosis in melanoma (38). However, the effects of the aforementioned RNAs on WT progression have not been reported at present. In addition, it remains unknown whether AC079310.1, LINC00271, AL445228.3, Z84485.1, AC091180.5, AP002518.2, AC007879.3, TCTEX1D4, RNF133, CCNE1, C10orf71, and SPAG11A can influence the development of tumors. These RNAs provide a potential direction for future investigations on WT.
The ceRNA hypothesis explains the complicated RNA molecular regulation mechanisms via construction of a lncRNA-miRNA-mRNA network (9). Most studies on the molecular regulation mechanism of WT constructed a network based on DERNAs (12,13). We proposed a novel survivalrelated ceRNA network that provided a new dimension for predicting the prognosis of patients with WT. Many cancerrelated molecules are included in this network, such as miR-181a, CCNE1, and WIF1. Recent studies have determined that some RNAs in the ceRNA network affect the prognosis of WT patients. For example, miR-15a targeted by SNHG6 (39), cyclin CCNE1 regulated by tumor suppressor gene WWOX (40), and Wnt signaling pathway inhibited by WIF1 (41) all affect tumor prognosis. In view of the fact that mRNA is the executor of ceRNA network function, GO and KEGG analyses were conducted to gain insight into molecular mechanisms. The GO analysis suggested that the cell surface receptor signaling pathway is the main mechanism by which the ceRNA network affects the prognosis of WT patients. KEGG analysis suggested that the "Wnt signaling pathway" and "Cellular senescence" were the main enrichment pathways. Recently, the Wnt/b-catenin signaling pathway was demonstrated to regulate the occurrence and development of WT, and Wnt-targeted agents may have great potential in the treatment of WT (42, 43).
Our study had several limitations. First, due to the scarcity of patients, we only obtained RNA sequencing data from the TARGET database; thus, we could not perform multicenter validation. Second, all the data in this study are from the same pathological tissue and lacked certain accuracy for highly heterogeneous WT. In the future, we will use single-cell sequencing to detect the expression RNAs in multiple WT samples from the same patient to improve the accuracy of the results. Finally, some new RNAs that may be involved in WT progression need further study.
To conclude, we established a multinomial predictive survival model, which has a good application prospect in future clinical practice. It is expected to improve the long-term prognosis of WT patients by screening high-risk patients through sequencing results and strengthening the personalized treatment. Meanwhile, we describe a survival-related ceRNA network, which provides some new potential prognostic indicators and therapeutic targets for improving the prognosis of WT patients.