Identification of a five‐mRNA signature as a novel potential prognostic biomarker in pediatric Wilms tumor

Abstract Background The aim of this study was to generate a prognostic model to predict survival outcome in pediatric Wilms tumor (WT). Methods The data including mRNA expression and clinical information of pediatric WT patients were downloaded from the Therapeutically Available Research to Generate Effective Treatments (TARGET) database. The differentially expressed genes were identified and a prognostic signature of pediatric WT was generated according to the results of univariate and multivariate Cox analysis. Receiver operating characteristic (ROC) curve was used to evaluate the five‐mRNA signature in pediatric Wilms tumor patients. Bootstrap test with 500 times was used to perform the internal validation. Results We identified 6,964 differentially expressed mRNAs associated with pediatric WT, including 3,190 downregulated mRNAs and 3,774 up‐regulated mRNAs. Univariate and multivariate Cox analysis identified five mRNAs (SPRY1, SPIN4, MAP7D3, C10orf71, and SPAG11A) to establish a predictive model. The risk score formula is as follows: Risk score = 0.3036*SPIN4 + 0.8576*MAP7D3 −0.1548*C10orf71 −0.7335*SPRY1 −0.2654*SPAG11A. The pediatric WT patients were divided into low‐risk group and high‐risk group based on the median risk score (value = 1.1503). The receiver operating characteristic (ROC) curve analysis revealed good performance of the 5‐mRNA prognostic model (the area under the curve [AUC] was 0.821). Bootstrap test (Bootstrap resampling times = 500) was used to perform the internal validation and revealed that the AUC was 0.822. REACTOME, KEGG, and BIOCARTA pathway analyses demonstrated that these survival‐related genes were mainly enriched in ErbB2 and ErbB3 signaling pathways, and calcium signaling pathway. Conclusion The five‐mRNA signature can predict the prognosis of patients with pediatric WT. It has significant implication in the understanding of therapeutic targets for pediatric WT patients. However, further study is needed to validate this five‐mRNA signature and uncover more novel diagnostic or prognostic mRNAs candidates in pediatric WT patients.

Besides stage and histology, a variety of clinical and biological factors was used to define treatment, including age, tumor size and volume, the loss of heterozygosity at chromosomes 1p and 16q, and response to chemotherapy (Dome et al., 2013;Dome, Perlman, & Graf, 2014). The treatment studies of children with WT have been evaluated by two different clinical groups, including COG Renal Tumor Committee (COG RTC) (D'Angio et al., 1989) and SIOP (Graf, Tournade, & de Kraker, 2000). The standard approach to WT treatment in the COG RTC group was immediate surgery, while the first step in treatment in SIOP was preoperative chemotherapy. Postoperative chemotherapy was used in both groups. The long-term survival outcomes of pediatric WT patients has improved gradually in the recent years, however, the subsequent chronic health conditions, including renal failure, cardiac toxicity, and subsequent malignancies should not be ignored (Aldrink et al., 2018;Gratias et al., 2016;Wong et al., 2016).
High risk groups compose 25% of patients with WT, including those with unfavorable histological, bilateral disease, and recurrence disease (Dome et al., 2015). As we all know, the heterogeneity among individuals often makes conventional prognostic systems. For instance, the risk stratification of TNM staging system is not sufficient. Besides, it is also insufficient to provide an accurate estimation of survival outcome. Thus, it is urgent to generate an accurate prognostic model to predict the survival outcomes in pediatric WT patients. Prognostic model plays a crucial role in the management of tumors, such as prostate specific antigen, alpha fetoproteinca, and carcinoembryonic antigen. Although a meta-analysis (Cone et al., 2016) reported that a large number of tumor biomarkers have been used to predict the prognostic outcomes in pediatric WT, there has been no prior study which has focused on an mRNA signature to predict the prognosis of WT patients.
The present study aimed to conduct an integrated study to develop a five-mRNA signature for the prognostic predication of WT patients by analyzing pediatric WT patients from Therapeutically Available Research to Generate Effective Treatments (TARGET) database.

| Acquisition of TARGET pediatric WT data
The RNA-seq data (level 3) and corresponding clinical information of pediatric WT in TARGET database were downloaded from Genomic Data Commons Data Portal (portal.gdc.cancer.gov/). We identified 136 cases investigated in this study, including 6 normal samples and 130 WT samples. No further normalization was needed for the expression data downloaded from TARGET database which have already been normalized. The data with no expression were deleted previously. The level 3 RNA-seq data between normal tissues and WT tissues were analyzed by edgeR package based on R language for differential expression analysis. Genes with absolute log 2 fold change > 1 and p < .05 were regarded as differentially expressed mRNAs. Since the data come from the TARGET database, no further approval was required from the Ethics Committee.

| Survival analysis
Clinical data were combined with those pediatric patients with WT in TARGET database to identify the prognostic differential expressed mRNAs signature. The survival curves of those samples with differential expressed mRNAs were plotted by using the "survival" package in R. The primary endpoint was overall survival. Univariate Cox analysis and multivariate Cox analysis were performed in this study. All identified differential expressed mRNAs were performed by univariate Cox analysis. The hazard ratio and P value of all differential expressed mRNAs were calculated. Receiver operating characteristic (ROC) curve has been used to prove the sensitivity and specificity of the calculated riskscore in predicting the overall survival of pediatric WT patients. The area under the curve (AUC) was generated and bootstrap was used to estimate 95%CI with the AUC.

| Statistical analysis
Kaplan-Meier survival analyses were used to determine the overall survival of pediatric Wilms patients who were classified as high expression and low expression group based on the median expression level of each differentially expressed mRNA. Log-rank test with the R package "survival" was used to determine the difference in the survival of pediatric patients. p < .05 was considered as statistically significant.

| Survival analysis by Kaplan-Meier method among differentially expressed mRNAs in pediatric WT patients
We identified 6,964 differentially expressed mRNAs, including 3,190 downregulated mRNAs and 3,774 upregulated mRNAs. Survival analyses among each deferentially expressed mRNAs were performed by Kaplan-Meier method subsequently. The high expression and low ex-

| Survival analysis by univariate cox analysis and multivariate cox analysis among differentially expressed mRNAs in pediatric WT patients
Univariate Cox analysis for all differentially expressed mRNAs was assessed to determine the survival-related mRNAs (Table S1). The primary endpoint for survival analysis was overall survival. The significant level cutoff threshold was set as 0.001 (p < .001) to identify the candidate mRNAs (Table 1). Multivariate Cox analysis was then performed by using these candidate mRNAs identified by univariate Cox analysis. Finally, five mRNAs (SPRY1, SPIN4, MAP7D3, C10orf71, and SPAG11A) were identified (Table 1). The results of multivariate Cox analysis also revealed the independent prognostic value of these 5 hub mRNAs. Two were associated with high risk of death in pediatric WT (SPIN4 and MAP7D3). SPIN4 and MAP7D3 were associated with a poor overall survival of pediatric WT patients. Specifically, the risk of death in patient with high expression of SPIN4 was 1.355 times higher than patient with low expression of SPIN4. Plus, the risk of death in patient with higher expression of MAP7D3 was 2.358 times higher than patient with low expression of MAP7D3. Also, three were associated with low risk of death in pediatric WT (SPRY1, C10orf71, and SPAG11A). the risk of death in patients with low expression of SPRY1 was 2.083 times higher than patient with high expression of SPRY1. Plus, the risk of death in patients with low expression of C10orf71 was 1.167 times higher than patients with high expression of C10orf71. Also, the risk of death in patients with low expression of SPAG11A was 1.304 times higher than patients with high expression of SPAG11A.

| The development of the 5-mRNA prognostic model
For each patient, a risk score analysis was conducted among the five mRNAs to determine the risk score (Table 1)  of five-mRNA signature was demonstrated in Figure 2c. The pediatric WT patients were divided into low-risk group and high-risk group based on the median risk score (value = 1.1503). Survival analysis between high-risk group and low-risk group was performed by using the log-rank test (Figure 2d). The result revealed that low-risk group was related to a better prognosis (p < .001  (Figure 2f), which was consistent with primary results of AUC (0.821). The results demonstrated that the 5-mRNA prognostic model had a promising sensitivity and specificity in predicting the survival outcomes of pediatric WT patients.

| REACTOME, KEGG, and BIOCARTA pathway analyses among survivalrelated mRNAs
We then included 466 survival-related mRNAs screened by univariate Cox analysis (p < .05) into pathway analyses. A total of 47 pathway ways were enriched in this study, including 16 pathways enriched by KEGG database, 29 pathways enriched by REACTOME database, and 2 pathways enriched by BIOCARTA database. The top five enriched pathways ranked as the P value were demonstrated in Table 2. The results demonstrated that these survival-related genes were mainly enriched in ErbB2 and ErbB3 signaling pathways, and calcium signaling pathway.

| DISCUSSION
Evidence has proved that mRNAs play crucial roles in the tumorigenesis and progression of pediatric WT (Apelt et al., 2016;Martins, Pinto, Domingues, & Cavaco, 2018;Zhu et al., 2018). Although several previous studies have identified several mRNAs with prognostic value in pediatric WT, they were not focused on the correlations between mRNA signature model and prognosis of pediatric WT (Gadd et al., 2017;Ludwig et al., 2016;Wari et al., 2017). Moreover, with the development of detection technology, the single mRNA expression pattern was no longer sufficient for accurate predication of prognosis of pediatric WT.
To the best of our knowledge, it is the first time to screen out the DEGs between pediatric WT and paired tissues from TARGET database. A novel five-mRNA signature (SPRY1, SPIN4, MAP7D3, C10orf71, and SPAG11A) was then generated. ROC curve proved that this five-mRNA signature revealed a high sensitivity and specificity in predicating the survival outcomes of pediatric WT patients. The predictive value of the five-mRNA signature was validated in TARGET dataset of 136 pediatric WT patients. Based on these five prognostic mRNAs, we established a five-mRNA prognostic model which can classify pediatric WT patients into low-risk and high-risk groups with different survival outcomes.
Wilms tumor are most common types of childhood kidney cancers. It has been reported that for children younger than 15 years with Wilms tumor, the 5-year survival rate has increased over the same time from 74% to 88% (Smith, Altekruse, Adamson, Reaman, & Seibel, 2014). The 5-year survival rate for Wilms tumor with favorable histology has been consistently above 90% since the 1980s (Smith, Altekruse, Adamson, Reaman, & Seibel, 2014). The results of this manuscript demonstrated that among these five mRNAs, SPRY1 and SPIN4 were associated with high risk of development of pediatric WT, and MAP7D3, C10orf71, and SPAG11A were associated with low risk of development of pediatric WT.
In mammals, SPRY1 was reported to be consisted of four members and was inhibitor of receptor tyrosine kinase signaling (Rozen et al., 2009). In mice, SPRY1 plays an important role during kidney morphogenesis by antagonizing GDNF signaling (Basson et al., 2005). SPRY1 also plays an important role in the early steps of glomerulus formation and represents a physiologically associated target gene of WT1 during the development of kidney (Gross et al., 2003). SPRY1 was reported to be associated with many kinds of tumors, such as breast cancer (He et al., 2016), colorectal cancer (Zhang et al., 2016), and human epithelial ovarian cancer (Masoumi-Moghaddam, Amini, Wei, Robertson, & Morris, 2015). The protein encoded by MAP7D3 belongs to the MAP7 family.
There is little known about the role of MAP7 with respect to cancer progression (Blum et al., 2008). Many important cellular processes attributed to microtubules involvement, including cell division, motility, and changes in cell shape (Bhat & Setaluri, 2007). Yan et al. (Yan et al., 2013) showed that miR-16 targeting MAP7 played an important role in regulating proliferation in cancer cells. Also, Blum et al. (Blum et al., 2008) demonstrated that the expression ratio of MAP7/B2M can be regarded as a prognostic factor for survival in patients with colon cancer. Peng et al.  demonstrated that SPAG11A was involved in the biological process of papillary thyroid cancer. However, the SPIN4 and C10orf71 have not been reported associated with the development and progression of cancer. To the best of our knowledge, given the potential molecular mechanism of the five mRNAs signature, no reports of the function and mechanism of these five mRNAs, SPRY1, SPIN4, MAP7D3, C10orf71, and SPAG11A, have been published concerning WT. The development of pediatric WT is a multi-step process. A large number of genetic alterations were involved in this multi-step biological process (Morrison, Viney, Saleem, & Ladomery, 2008). For the sake of elucidating the effects and functions of these survival-related mRNAs screened by univariate Cox analyses, we used REACTOME, KEGG, and BIOCARTA pathway databases to perform pathway analyses. The results demonstrated that these survival-related genes were mainly enriched in ErbB2 and ErbB3 signaling pathways and calcium signaling pathway.
Both BIOCARTA and REACTOME pathway databases revealed that these survival-related genes were mainly enriched in ErbB2 and ErbB3 signaling pathways. ErbB2 and ErbB3 belong to the family of human epidermal growth factor receptors consisting of EGFR (ErbB1), ErbB2, ErbB3, and ErbB4 (Vermeulen, Segers, & De Keulenaer, 2016). ErbB2 amplification plays a critical role in tumor growth. Amplified ErbB2 can bind to ErbB3 to form an oncogenic ErbB2/ErbB3 complex (Holbro et al., 2003). ErbB3 interacts with the regulatory p85 subunit of PIK3 in this complex to activate the PI3K/Akt pathway and intense cell growth and proliferation. (Schoeberl et al., 2009) Therefore, ErbB3 plays an important role in oncogenic ErbB2 signaling pathway. Rotter et al. (Rotter, Block, Busch, Thanner, & Hofler, 1992) reported that the expression of ErbB2 was downregulated in the renal cell carcinoma when compared with normal kidney tissue. To the best of our knowledge, the molecular mechanisms behind the alteration of ErbB2 in renal cell carcinoma compared with normal kidney was still unknown. Plus, the expression of ErbB3 has not been thoroughly studied in renal cell carcinoma. KEGG pathway database revealed that these survival-related mRNAs were mainly enriched in calcium signaling pathway. Previous studies (Cole & Kohn, 1994;Soboloff, Zhang, Minden, & Berger, 2002;Sukumaran, Sun, Vyas, & Singh, 2015) have been reported that inhibition of calcium influx can cause either growth arrest or cell death in numbers of cancer cells. However, the role of calcium signaling pathway in the development and progression of WT has not been elucidated yet. Xu et al. (Xu, Chen, Ye, Zhong, & Chen, 2015) reported that calcium signaling pathway has been involved in inducing the apoptosis in non-small cell lung cancer cells, for the overload of calcium has been reported to play a crucial role in the initiation and regulation of apoptosis.
There are some limitations in this study. The predictive value of the five-mRNA signature was not validated in another independent dataset because it is very difficult for us to obtain tumor specimens, especially pediatric tumor samples.

| CONCLUSION
In conclusion, the five-mRNA signature can predict the prognosis of patients with pediatric WT. It has significant implication in the understanding of therapeutic targets for pediatric WT patients. However, further study is needed to validate this five-mRNA signature and uncover more novel diagnostic or prognostic mRNA candidates in pediatric WT patients.