Potential Five-MicroRNA Signature Model for the Prediction of Prognosis in Patients with Wilms Tumor.

BACKGROUND Wilms tumor (WT) is the most common type of pediatric renal malignancy, and is associated with poor prognosis. The aim of the present study was to identify microRNA (miRNA) signatures which might predict prognosis and categorize WTs into high- and low-risk subgroups. MATERIAL AND METHODS The miRNA expression profiles of WT patients and normal samples were obtained from the Therapeutically Applicable Research to Generate Effective Treatment database. Differentially expressed miRNAs between WT patients and normal samples were identified using the EdgeR package. Subsequently, correlations between differentially expressed miRNAs and the prognosis of overall survival were analyzed. Enrichment analyses for the targeted mRNAs were conducted via the Database for Annotation, Visualization, and Integration Discovery. RESULTS A total of 154 miRNAs were identified as differentially expressed in WT. Of those, 18 miRNAs were associated with overall survival (P<0.05). A prognostic signature of 5 differentially expressed miRNAs (i.e., has-mir-149, has-mir-7112, has-mir-940, has-mir-1248, and has-mir-490) was constructed to classify the patients into high- and low-risk subgroups. The targeted mRNAs of these prognostic miRNAs were primarily enriched in Gene Ontology terms (i.e., protein autophosphorylation, protein dephosphorylation, and stress-activated MAPK cascade) and the Kyoto Encyclopedia of Genes and Genomes signaling pathways (i.e., MAPK, AMPK, and PI3K-Akt). CONCLUSIONS The 5-miRNA signature model might be useful in determining the prognosis of WT patients. As a promising prediction tool, this prognosis signature might serve as a potential biomarker for WT patients.

Background

Wilms tumor (WT) is a rare common type of pediatric renal malignancy, affecting approximately one in 10 000 children (aged <15 years) [1]. WT is associated with a poor prognosis. However, the rate of 5-year overall survival (OS) is constantly improving in response to the advancement of disease-associated therapies [2]. Prognostic tumor biomarkers play a vital role in the diagnosis and treatment of tumors. The identification of low-risk WT patients might lead to the administration of treatment with less morbidity. In addition, clinicians might apply more aggressive treatments in WT patients identified as at high risk of poor prognosis according to the results of the prognostic biomarker test. Therefore, novel biomarkers or prognostic models for the prediction of survival in WT are urgently warranted. This research regarding prognostic biomarkers might be an important step in improving the prognosis of patients with WT.

MicroRNAs (miRNAs) are single-stranded (18–25 nucleotide-long) noncoding RNAs which target mRNAs to control gene expression [3]. MiRNAs restrain protein synthesis through de-adenylation and repression of translation, and the subsequent degradation of mRNA targets [4]. The expression of >30% of human genes is regulated by miRNAs, indicating that miRNAs might play a pivotal role in numerous biological functions [5]. Aberrant expression of miRNAs is realized as an important step in the development of cancer [6]. Notably, the expression of miRNAs has been shown to be stable even in samples obtained from formalin-fixed paraffin-embedded tissues [4]. Such dysregulated miRNAs might be promising therapeutic targets. Previous studies reported that numerous key miRNAs (e.g., miR-483-5p [7] and miR-195 [8]) are closely related to the pathogenesis of WT. However, these studies were focused on a limited number of investigated miRNAs from a small number of patients or different miRNA-chip platforms. It is likely that multiple miRNAs (signature model) might provide more statistically predictive results. Thus, there is an urgent need to apply larger patient cohort data to the investigation of more specific prognostic classifiers in WT patients.

The Therapeutically Applicable Research to Generate Effective Treatment (TARGET) database has publicized large amounts of miRNA sequencing data, which might provide insight into the molecular mechanism involved in the progression and tumorigenesis of WT. The aim of the present study was to identify patterns of differentially expressed miRNAs (DEMs) between normal and WT samples via the use of advanced RNA sequencing analysis. Additionally, the prognostic value of these DEMs was evaluated. Furthermore, a 5-miRNA signature model was constructed to effectively predict the prognosis of patients with WT.

Material and Methods

Data collection and preprocessing

All miRNA expression data were downloaded from the publicly available TARGET database. The miRNA expression data and corresponding clinical information of 138 samples, including 6 normal samples and 132 WT samples, were analyzed in the present study. Considering that the data were obtained from the TARGET database, approval of the study by an ethics committee was not required. The EdgeR package in the R software (version 3.2.2) was used to analyze DEMs between normal and WT samples. The cutoff criteria were defined as |logFC(fold change)| >2 and a false discovery rate of < 0.05.

Risk stratification and receiver operating characteristic (ROC) curve

DEMs with a P<0.05 in the univariate Cox regression analysis were considered as candidate variables for inclusion in the multivariate Cox regression analysis. The risk scores of 127 patients with WT were calculated according to the predictive 5-miRNA signature model. The classification cutoff value between the low- and high-risk subgroups was defined as the median risk score. Differences in OS time between the low- and high-risk WT subgroups were compared via the log-rank test and Kaplan-Meier analysis. The receiver operating characteristic (ROC) curve and area under curve (AUC) values were used to evaluate the specificity and sensitivity of this 5-miRNA prognosis signature model. All the analyses were conducted using the R software (version: 3.3.2).

Survival analysis

Prognostic DEM signatures were identified via the log-rank test and Kaplan-Meier analysis. The survival curves were constructed using the “survival” package in the R software (version: 3.3.2). The OS was defined as the survival endpoint in the present survival analysis. A P<0.05 denoted statistical significance.

Prediction of target genes and functional enrichment analysis

The online tools MiRTarBase [9] (), MiRDB [10] (), and TargetScan [11] () were used to predict target mRNAs for all prognostic miRNAs. To minimize false positives and screen most putative targets, target mRNAs, which were predicted by all 3 online programs, were selected as input for the functional enrichment analyses. Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEMs were performed via the Database for Annotation, Visualization, and Integrated Discovery in an attempt to elucidate the mechanism of tumorigenesis in WT.

Results

DEMs in WT

The miRNA expression profiles between normal and WT samples acquired from the TARGET database were analyzed. According to the cutoff threshold of |logFC(fold change)| >2 and a false discovery rate of <0.05, 105 upregulated and 49 downregulated DEMs were identified in WT. The heatmap plot and volcano plot of the related DEMs between normal and WT samples are presented in Figure 1. The top 10 upregulated and downregulated DEMs are shown in Table 1.

Figure 1

Heatmap analysis (A) and volcano plot (B) of differentially expressed micoRNAs.

Table 1

The top 10 upregulated and downregulated differentially expressed microRNAs.

Genesymbol	LogFC	FDR	Change
hsa-mir-934	−8.458976144	8.18E-109	Down
hsa-mir-203a	−5.559616525	2.43E-29	Down
hsa-mir-29c	−4.184218278	8.62E-46	Down
hsa-mir-506	−4.068436968	7.11E-10	Down
hsa-mir-29a	−4.000404843	1.54E-26	Down
hsa-mir-30a	−3.996893892	7.42E-27	Down
hsa-mir-514a-1	−3.966692055	6.42E-12	Down
hsa-mir-514a-2	−3.896476522	2.33E-13	Down
hsa-mir-514a-3	−3.863700282	5.84E-12	Down
hsa-mir-203b	−3.503204575	1.32E-06	Down
hsa-mir-383	9.298659023	4.06E-06	Up
hsa-mir-1269a	7.525185748	0.000415743	Up
hsa-mir-1269b	6.968194063	0.046671257	Up
hsa-mir-767	6.775915413	0.006257893	Up
hsa-mir-548f-1	6.766938092	1.24E-07	Up
hsa-mir-2115	6.512683834	4.88E-05	Up
hsa-mir-301b	6.276573984	1.89E-10	Up
hsa-mir-105-1	6.260382994	0.019280135	Up
hsa-mir-483	6.210785374	1.38E-11	Up
hsa-mir-105-2	5.578850103	0.025007756	Up

FC – fold change; FDR – false discovery rate.

Establishment of a five-miRNA signature model was associated with OS of WT patients

Univariate and multivariate Cox regression analyses were used to evaluate associations between the expression levels of DEMs and the OS of WT patients. Five DEMs (i.e., has-mir-149, has-mir-7112, has-mir-940, has-mir-1248, and has-mir-490) were selected from the stepwise multivariate Cox regression analysis to establish a predictive signature model. In the present study, this model was defined as follows: survival risk score=(0.263391572×expression value of has-mir-149) + (−0.229979121288481×expression value of has-mir-7112) + (0.427344897140062×expression value of has-mir-940) + (0.227391121990691×expression value of has-mir-1248) + (−0.147764549439855 x expression value of has-mir-490).

Risk stratification and ROC curve analysis

The 5-miRNA expression-based scores of 127 WT patients were calculated. According to the median risk score, these WT patients were assigned to a high-risk group (n=63) or a low-risk group (n=64) (Figure 2). In the current study, the AUC value of the ROC curve was applied to evaluate the prognostic power of the 5-miRNA signature model (Figure 3). Higher AUC values were indicative of greater predictive ability by the model. In addition, AUC values >0.70 denoted excellent model performance. In the signature model of the present study, the AUC value obtained from the ROC curve analysis was 0.767. This finding indicated good specificity and sensitivity of the 5-miRNA signature model in predicting the survival risk of WT patients.

Figure 2

MicroRNA predictor-score analysis of Wilms tumor patients. (A) Heatmap of miRNA expression profiles of Wilms tumor patients. (B) MicroRNA predictor-score distribution. (C) Patients’ survival status.

Figure 3

Kaplan-Meier and (receiver operating characteristic) ROC curves for the 5-miRNA signature model. (A) The Kaplan-Meier curves for entire Wilms tumor cohort divided by the optimum cutoff point. (B) The ROC curve for predicting overall survival in the Wilms tumor cohort.

A Kaplan-Meier univariate survival analysis was performed by combining information of gene expression profiles and the clinical information of WT samples. Among the 154 DEMs examined in this study, 18 miRNAs (i.e., hsa-mir-29b-2, hsa-mir-99a, hsa-mir-135a-1, hsa-mir-135a-2, hsa-mir-149, hsa-mir-181b-2, hsa-mir-200a, hsa-mir-483, hsa-mir-490, hsa-mir-509-2, hsa-mir-636, hsa-mir-660, hsa-mir-934, hsa-mir-940, hsa-mir-1248, hsa-mir-1269a, hsa-mir-3651, and hsa-mir-5002) were closely related to OS in WT patients (Figure 4). Low expression of hsa-mir-29b-2, hsa-mir-149, hsa-mir-483, hsa-mir-636, hsa-mir-934, hsa-mir-940, hsa-mir-1248, and hsa-mir-3651, was related to a high OS rate in WT patients. Notably, high expression of hsa-mir-99a, hsa-mir-135a-1, hsa-mir-135a-2, hsa-mir-181b-2, hsa-mir-200a, hsa-mir-490, hsa-mir-509-2, hsa-mir-660, hsa-mir-1269a, and hsa-mir-5002 was related to high survival rates in WT patients.

Figure 4

(A–R) Prognostic significance of differentially expressed microRNAs.

Prediction of target mRNA and functional enrichment analysis

The online analysis tools TargetScan, miRDB, and miRTarBase were used to predict target mRNAs of prognostic DEMs. The overlapping target genes predicted by all 3 online programs were identified as target mRNAs. The results of the GO analysis showed that target mRNAs were mainly enriched in protein autophosphorylation, protein dephosphorylation, the stress-activated MAPK cascade, regulation of binding, and regulation of organ growth (Figure 5A, Table 2). According to the KEGG analysis, 6 statistically significant pathways were identified: MAPK signaling pathway, signaling pathways regulating the pluripotency of stem cells, AMPK signaling pathway, endocrine resistance, PI3K-Akt signaling pathway, and prostate cancer pathway (Figure 5B).

Figure 5

Functional Gene Ontology (A) and Kyoto Encyclopedia of Genes (B) enrichment analysis of targeted mRNAs.

Table 2

The Enriched Gene Ontology Terms (top 10) and Kyoto Encyclopedia of Genes and Genomes of Target mRNAs.

Category	Term	ID	FDR	Targeted mRNA
GO	Focal adhesion	GO: 0005925	0.015653541	MYH9/ARHGAP31/ITGA3/CALR/GIT1/ARHGAP22/DLC1/GNA13/KLF11/EPHA2/STX16/ALCAM/YWHAE/ADD1/FGFR3
GO	Cell-substrate adherens junction	GO: 0005924	0.015653541	MYH9/ARHGAP31/ITGA3/CALR/GIT1/ARHGAP22/DLC1/GNA13/KLF11/EPHA2/STX16/ALCAM/YWHAE/ADD1/FGFR3
GO	Cell-substrate junction	GO: 0030055	0.015653541	MYH9/ARHGAP31/ITGA3/CALR/GIT1/ARHGAP22/DLC1/GNA13/KLF11/EPHA2/STX16/ALCAM/YWHAE/ADD1/FGFR3
GO	Protein autophosphorylation	GO: 0046777	0.036590366	AKT1/MINK1/FGFR1/MARK2/STK11/RAP2C/EPHA7/MAP3K9/TAOK1/EPHA4/FGFR3/MTOR
GO	Protein dephosphorylation	GO: 0006470	0.036590366	CTDSP2/DUSP7/STK11/PTPRK/DLC1/UBASH3B/PTPRD/PPP1R15B/ARPP19/DUSP1/YWHAE/CTDSPL/MTMR3/MTOR
GO	Stress-activated MAPK cascade	GO: 0051403	0.036590366	AKT1/MINK1/WNT7B/ZFP36L1/FOXM1/ZEB2/ZNF675/MAPKAPK2/MAP3K9/TAOK1/EPHA4/DUSP1/HMGB1
GO	Regulation of binding	GO: 0051098	0.036590366	AKT1/MARK2/CARM1/DISC1/UBASH3B/ZNF675/FOXC1/IGF1/EPHA4/SYAP1/PYGO2/RAPGEF2/HMGB1/STPG1/ADD1
GO	Regulation of organ growth	GO: 0046620	0.036590366	AKT1/FGFR1/CARM1/YAP1/FOXC1/IGF1/RBPJ/MTOR
GO	Ras protein signal transduction	GO: 0007265	0.037327922	RAB11B/PLEKHG2/IQSEC3/ITGA3/NGFR/RAB7A/FOXM1/RAP2C/GRB2/OGT/DLC1/GNA13/RAB8B/RAB5B/IGF1/RAPGEF2
GO	Positive regulation of organ growth	GO: 0046622	0.037327922	AKT1/FGFR1/YAP1/IGF1/RBPJ/MTOR
KEGG	MAPK signaling pathway	hsa04010	0.000605175	AKT1/MKNK2/FGFR1/NGFR/CACNB1/DUSP7/GRB2/EPHA2/MAPKAPK2/IGF1/TAOK1/FOS/DUSP1/RAPGEF2/FGFR3
KEGG	Prostate cancer	hsa05215	0.018238418	AKT1/FGFR1/GRB2/CCNE2/ZEB1/IGF1/MTOR
KEGG	Signaling pathways regulating pluripotency of stem cells	hsa04550	0.020225637	AKT1/ZFHX3/FGFR1/WNT7B/GRB2/IGF1/SKIL/FGFR3
KEGG	AMPK signaling pathway	hsa04152	0.033440098	AKT1/RAB11B/STK11/SREBF1/SCD5/IGF1/MTOR
KEGG	Endocrine resistance	hsa01522	0.044793852	AKT1/CARM1/GRB2/IGF1/FOS/MTOR
KEGG	PI3K-Akt signaling pathway	hsa04151	0.044793852	AKT1/ITGA3/FGFR1/NGFR/STK11/GRB2/CCNE2/EPHA2/IGF1/YWHAE/FGFR3/MTOR

GO – Gene Ontology, KEGG – Kyoto Encyclopedia of Genes; FDR – false discovery rate.

Discussion

WT is a type of pediatric renal malignancy. Poor prognosis remains the main cause of cancer-specific death in patients with WT [2]. The current treatment strategies against WT rely on attempts to investigate biomarkers reflecting prognosis and predict the risk of patients. MiRNAs are promising biomarkers for tumor diagnosis and prognosis. In recent years, miRNAs have been reported as prognostic biomarkers for the diagnosis and prognosis of WT patients. WT is a multifactorial and highly complex tumor disease, and numerous dysregulated miRNAs are involved in tumorigenesis and progression. Therefore, a model involving the combination of multiple miRNAs might provide more detailed and accurate information regarding the diagnosis and prognosis of WT versus a single miRNA.

MiRNAs target mRNAs to control gene expression [3]. In the present study, a total of 105 upregulated and 49 downregulated DEMs were identified. Of the 154 DEMs detected, 18 DEMs were associated with OS in WT patients (P<0.05). Several DEMs detected via our analysis have been previously investigated in WT. For example, miR-483 was shown to cooperate with IGF2 and serve as an autonomous oncogene in WT [12]. The oncogenic mechanism of miR-483 was partially demonstrated by research findings indicating that miR-483 might modulate the proapoptotic protein (BBC3/PUMA) to protect cells from apoptosis [12]. Additionally, information regarding the regulatory role of hsa-mir-29b-2, hsa-mir-99a, hsa-mir-135a-1, hsa-mir-135a-2, hsa-mir-200a, hsa-mir-181b-2, hsa-mir-509-2, hsa-mir-636, hsa-mir-660, hsa-mir-934, hsa-mir-1248, hsa-mir-1269a, hsa-mir-3651, hsa-mir-149, hsa-mir-940, hsa-mir-490, and hsa-mir-5002 in WT is currently limited. A previous study showed that these DEMs might be related to oncogenesis. For example, miR-29b-2 was reported to improve prognosis in pancreatic ductal adenocarcinoma by targeting Cbl-b and promoting the expression of p53 [13]. In addition, miR-99a might downregulate 2 novel oncogenic proteins (i.e., E2F2 and EMR2), and repress stemness in lung cancer [14]. Moreover, miR-135a-1 might target EGFR to inhibit cell growth and migration in prostate cancer [15]. Furthermore, miR-200a was shown to target FOXA1 and act as a tumor suppressor on the survival, proliferation, and invasion of glioma cells [16]. Previous studies also reported that miR-200a is associated with the development of cancer (i.e., esophageal, breast, and endometrial), by targeting specific genes such as CRMP1, EPHA2, and PTEN [17-19]. Notably, miR-509-2 was differentially expressed in early-stage cervical cancer and might serve as a diagnostic biomarker in these patients [20]. Also, miR-636 acted as a tumor suppressor miRNA, which might affect the tumorigenesis of hepatocellular carcinoma through the downregulation of Ras [21]. Suppression of ANT2 by shRNA might restore the expression of miR-636 to play an anticancer role in hepatocellular carcinoma [21]. Of note, mir-660 might reduce the expression of MDM2 and stabilize the protein levels of p53 to act as a tumor suppressor miRNA in lung tumors [22]. Furthermore, miR-934 was significantly upregulated in estrogen receptor-negative breast cancer, and its expression was related to that of VGLL1 [23]. Suppression of miR-1248 might induce AGTR1-mediated cancer cell death to improve the outcome of chemotherapy against osteosarcoma [24]. MiR-1269a promotes tumor metastasis and forms a positive feedback loop with TGF-β signaling in colorectal cancer [25]. The role of these DEMs in the progression of WT is currently not well understood. Therefore, additional targeted molecular studies are required to evaluate their involvement in this process. In addition, further studies are warranted to illuminate the molecular and biological mechanisms of these DEMs.

GO and KEGG analyses of 18 prognostic DEMs (targeted mRNAs) were conducted to further investigate the related cellular mechanisms in WT. The results of the GO analysis showed that that target mRNAs were mainly enriched in protein autophosphorylation, protein dephosphorylation, the stress-activated MAPK cascade, regulation of binding, and regulation of organ growth. The results of the KEGG analysis indicated that target mRNAs were mainly enriched in the MAPK signaling pathway, signaling pathways regulating the pluripotency of stem cells, AMPK signaling pathway, endocrine resistance, PI3K-Akt signaling pathway, and prostate cancer pathway. In recent years, several biological processes and signaling pathways identified through our analysis have been investigated in WT. The stress-activated MAPK signaling pathway might serve as a tumor suppressor and common pro-oncogenic signal in the process of tumorigenesis [26,27]. Notably, a new candidate drug induced apoptosis of WT cells through the MAPK pathway [28]. Berberine was reported to activate AMPK signaling for the regulation of tumor progression and metastasis in WT [29]. Recent evidence indicated that the PI3K-AKT-p53 signaling pathway participated in the tumorigenesis of WT. This finding might represent a potential treatment strategy in the future [30]. The aforementioned evidence was partially consistent with the results obtained through the GO and KEGG analyses in the present study.

An increasing number of studies suggest that the expression of miRNAs in WT plays an important role in diagnosis and prognosis. Ludwig et al. identified circulating miRNAs (i.e., miR-130b-3p, miR-100-5p, and miR-143-3p) with a high diagnostic potential in WT [31]. Apart from circulating miRNAs, miRNAs in WT tissues have also been studied. For example, miR-483-5p was shown to be predictive of an adverse prognosis in WT [32]. Koller et al. investigated the expression levels of miR-204 and found that its expression was significantly decreased in 23 WT samples versus normal renal parenchyma [33]. An upregulated expression of miR-21 was reported to have a possible correlation with aggressive progression and poor prognosis in WT [34]. In the present study, we demonstrated that a miRNA signature model – instead of a single miRNA – constructed from the TARGET database might be a suitable prognostic biomarker. A 5-miRNA model was established, and this signature model was used to assign patients to different risk subgroups according to the characteristics of their tumors. Similar multi-miRNA prognostic models have been established in numerous types of tumors. For example, a 4-miRNA model (i.e., miR-4758, miR-876, miR-142, and miR-190b) was identified to serve as an effective independent prognostic factor for the prediction of OS in patients with endometrial cancer [35]. In addition, an 8-miRNA model was reported to effectively predict OS in WT and assist clinicians in screening WT patients at high risk of unfavorable OS.

However, there were several limitations in the present study. Firstly, the miRNA expression data used for the signature model were obtained from a single database (TARGET). The performance of the 5-miRNAs signature model might be more accurately determined following additional validation with independent external miRNA expression data. The Gene Expression Omnibus and the Cancer Genome Atlas datasets were investigated for further validation. However, currently, there are no databases containing full survival information and miRNA expression data for further analysis. Secondly, this study applied bioinformatics data to construct the 5-miRNA signature model. Further experimental research validating the functions of the 5-miRNA signature model, and clinical studies investigating its predictive value are warranted.

Conclusions

We identified DEMs and constructed a 5-miRNA signature model, which might act as a potential predictor of prognosis in patients with WT. Significantly altered miRNAs might serve as prognostic biomarkers and therapeutic targets for WT. Further validation and functional studies are required to elucidate the molecular mechanisms of these miRNAs in WT.