A Risk Signature of Nine Autophagy-Related Genes Associated With Immune Microenvironment and Clinical Prognosis in Wilms Tumor

: Background: Autophagy is a biological process to eliminate dysfunctional organelles, aggregates or even long-lived proteins. . Nevertheless, the potential function and prognostic values of autophagy in Wilms Tumor (WT) are complex and remain to be clarifed. Therefore, we proposed to systematically examine the roles of autophagy-associated genes (ARGs) in WT. Methods: Here, we obtained differentially expressed autophagy-related genes (ARGs) between healthy and Wilms tumor from Therapeutically Applicable Research To Generate Effective Treatments(TARGET) and The Cancer Genome Atlas (TCGA) database. The functionalities of the differentially expressed ARGs were analyzed using Gene Ontology. Then univariate COX regression analysis and multivariate COX regression analysis were performed to acquire nine autophagy genes related to WT patients’ survival. According to the risk score, the patients were divided into high-risk and low-risk groups. The Kaplan-Meier curve demonstrated that patients with a high-risk score tend to have a poor prognosis. Results: Eighteen DEARGs were identifed, and nine ARGs were fnally utilized to establish the FAGs based signature in the TCGA cohort. we found that patients in the high-risk group were associated with mutations in TP53. We further conducted CIBERSORT analysis, and found that the infiltration of Macrophage M1 was increased in the high-risk group. Finally, the expression levels of crucial ARGs were verifed by the experiment, which were consistent with our bioinformatics analysis. Conclusions: we emphasized the clinical significance of autophagy in WT, established a prediction system based on autophagy, and identified a promising therapeutic target of autophagy for WT.


Introduction
Wilms tumor is the second most common intra-abdominal cancer in childhood and the fifth most common malignant tumor in childhood [1]. The annual incidence of the Wilms tumor is 8.1 per million children [2]. The rate of children under ten years old accounted for 98 , and most occur in children at the age of 3. The vast majority of children with nephroblastoma are inadvertently found abdominal lumps, such as touching the lumps when bathing the child, changing clothes, or touching the abdomen of the child. The death rate has dropped from 90 percent to 30 percent, and modern medicine shows excellent progress due to the operative treatment combined with radiotherapy and chemotherapy. The treatment of nephroblastoma is mainly divided into two primary treatment options: The COG/CCCG method supports direct resection of tumor tissue, while the SIOP method advocates preoperative chemotherapy by stages. However, the prognosis of patients treated with either up-front nephrectomy or preoperative chemotherapy have been excellent [3]. The therapeutic schedule of Wilms tumor relies on the individual risk of recurrence. Before any treatment begins, the personal risk has been assessed according to patients' clinical situations, pathological grading and prognostic factors. The research has suggested that the level of p73 in blood could act as a reliable biomarker for prior treatment of WT [4]. The expression of B7-H1 is correlated with an increased likelihood of early surgical treatment failure [5]. However, these exiting prognostic assessment factors have not fully considered the patients' consideration. Therefore, the primary task is to find a more accurate and appropriate evaluation model of WT patients. Autophagy is a biological process that occurs in both pathological and physiological processes, including the formation of autophagosomes, wrapping the contents which need to be degraded, and finally achieving the balance between energy metabolism and regeneration in the body. However, there are still some ideas that suggest a link between autophagy and Wilms tumors. At present, no systematic studies have focused on the predictive signals of autophagy-related genes construction. Here we collected and downloaded 126 WT patient' RNA-seq and clinical information from TCGA database. We then established a prognostic model that could serve as an independent risk factor, providing an important clinical reference for assessing both preoperative physical status and postoperative recurrence rates.

Data collection
The expression profiles of 232 autophagy-related genes were extracted from the Website: HADb (Human Autophagy Database). The RNA-seq data information of Wilms tumor patients was downloaded from the TARGET(Therapeutically Applicable Research To Generate Effective Treatments) database contained 130 WT tissues and 6 normal tissues. The clinical and OS data were simultaneously downloaded from this website.

Differently expressed ARGs and enrichment analysis
We used the limma package based on R to differentiate the ARGs expression between the normal tissue and pathological tissue. The screening criteria of different expressed autophagy-related genes(DE-ARGs) based on the following criterion: llog2 fold change (FC)I> 2 and adjusted P-value<0.05. To better understand the filtered DE-ARGs' biological function and signalling pathways. We employed the Gene Ontology (GO) analysis to explore the primary function of these genes in WT. The"cluster Profile" package and "GOplot" package were all used for analyses. Construction and validation of the prognostic model based on the autophagy-related gene Univariate Cox regression analysis was used to screen out genes related to the overall survival rate of WT patients. Then multivariate Cox regression analysis was performed for genes significantly associated with overall survival (P <0.05). Meanwhile, all genes were divided into two groups according to HR value (hazardous genome HR >1, protective genome HR <1). After this, genes with no or weak correlation with OS-associated markers will be excluded because they could not be independent prognostic indicators. The following formula was used to establish the prognostic model : risk score=IX J * coef J, the coef J represented the ARG proportion in this model, and ×J meaned each autophagy-related gene relative expression. Then, 130 WT patients were divided into two groups according to this standard, and survival analysis was conducted by Kaplan-Meier (K-M) method. The predictive values were conducted through the ROC curve employing the package of "survivalROC" in R.

Acquisition of somatic mutation data and estimation of TMB
We downloaded the somatic mutant spectrum of all tumor samples from the public database. The information contained four main data formats, and what we needed was the Mutation Annotation Format (MAF) of somatic information in TCGA database. Later, the "maftools" R package was utilized to visualize the samples' multiple gene mutation analysis results. We calculated the TMB score of each sample according to the total number of variants divided by the whole length of exons. Then Wilcoxon rank-sum test was used to analyze the difference between the low-risk and the high-risk group. The differential abundance of immune infiltrating cells among the two risk groups To better understand the distribution of tumor-infiltrating immune cells in both risk and low-risk groups. We used the CIBERSORT(http://cibersort.stanford.edu/) to calculate the proportion of immune cells. The results were presented in the format of a box diagram. Then the Wilcoxon ranked-sum test was conducted to find the differential abundance in the two groups.

Cell Culture
The human Human renal epithelial cell line 293T and WT cell line SK-NEP-1 were bought from the ATCC. The cells were cultured in a humidified environment with 5% CO2 at 37°C. The SK-NEP-1 cell line were cultured with McCoy's 5A supplemented with 15% fetal serum while the 293T cultured with high glucose DMEM supplemented with 10% fetal serum. The culture medium and supplements were purchased from HyClone (Northbrook, IL, USA). RNA Extraction and qRT-PCR Total RNA was extracted by Trizol reagent (Invitrogen,15596018). RNAs were reverse transcribed utilizing Evo M-MLV RT Kit with gDNA Clean for qPCR II (Accurate biology, AG11711). qRT-PCR was performed using the TB Green Fast qPCR Mix (TaKaRa,RR430S). The primer sequences were used as follows:

Result: Differentially expressed ARGs in Wilms tumor and functional verification of these genes
We collected RNA-seq and clinical data contained 130 WT tissue samples and 6 normal tissue from TARGET Database. The expression of 234 ARGs collected from HADb has been analyzed using the criterion of an FDR <0.05 and llog2(FoldChange)l>2. There were 24 differentially expressed genes, including 15 down-regulated genes and 9 up-regulated genes ( Figure 1A). We visualized these genes using the volcano plot and box plot( Figure 1B, C). To investigate the biological processes in which these 25 different genes may be involved, we performed a GO enrichment analysis. The DAVID analysis presented the top 3 of GO enriched terms in the biological processes, which consist of regulation of cytokine-mediated signaling pathway, cytokine stimulus and the neuron death regulation. The cellular components included autophagosome, autophagosome membrane and aggresome. The molecular function included ubiquitin binding, CARD domain binding and ubiquitin-like protein binding (Figure 2A,B).

Construction of the prognostic signature according to the ARGs
The univariate Cox regression analysis was utilized to calculate that 18 ARGs were strongly correlated with the OS of the Wilms tumor patients( Figure 3A). After that, Multivariate Cox regression analysis was used to select 9 ARGs that could form a prognostic model ( Figure  3B). It can be seen that NFE2L2, MAPK8 and ERBB2 are protective factors for WT patients, while SH3GLB1,SESN2, ATG5, VAMP, MTMR14 and ATG10 are related to poor prognosis ( Table 1). The model was established following this formula: Risk score=(-1.1527×NFE2L2 expression)+(-0.9022 ×MAPK8 expression)+(-0.6219 × ERBB2 expression)+(0.5934× SH3GLB1 expression)+(0.6309× SesN2 expression)+(0.6941× ATGS7 expression)+(0.8902× MTMR14 expression)+(1.2844×ATG10 expression). We then calculated each patients' risk score and divided them into two groups based on the risk scores. After obtaining the Kaplan-Meier overall survival (OS) curve, the results showed that the prognosis of the high-risk group was worse than that of the low-risk group (P = 2.458E-06). ( Figure 4A). The ROC curves showed the accuracy of this predictive model( Figure 4B). The Kaplan-Meier curves of each autophagy-related gene that participated in constructing this predictive model were shown in the additional figure(Supplement 1) . Each sample of Wilms tumor patients downloaded from the TCGA was visualized( Figure 4D,E). These dots represented individual patients, and the mortality rate increased as the risk score increased.

The Prognostic Model act as an independent factor
We also used univariate Cox regression analysis to investigate whether age, sex, first event, stage, histological classification, and risk score could be prognostic factors in Wilms tumor patients. The risk score assessed by our model formula could provide clinical value to predict the prognosis, with the Hazard ratio of 1.267 (p < 0.001). Multivariate Cox regression analysis further demonstrated that the prognostic model was effective compared with other commonly used clinical indicators. The Hazard ratio was 1.140 (p=0.006) ( Table 2). The ROC curve showed that OS-related prediction signatures had better accuracy and better estimation performance than other clinical variables, and the area under the ROC curve (ROC) was 0.819 (Supplement 2).

The correlations between the risk group and clinicopathological parameter
Further studies were needed to clarify the relationship between prognostic models and clinicopathological parameters. Chi-square test showed that the risk score calculated by the formula was related to sex, stage and first onset, but not to age or histological classification. The heatmap shows the correlations between the ultimate prognostic gene and clinicopathologic features( Figure 5A). The boxplot illustrates in more detail the correlation between risk scores and clinical characteristics (Supplement 3).

TMB and mutation landscape of risk model
We obtained somatic mutation profiles of 39 WT patients from the TARGET database. Then we used the "maftols" package to present a visual mutation difference among the two groups. The waterfall map showed the mutation information of the unstable genes in each group. We found that the low-risk group altered in 17(68) of 25 samples compared with the high-risk group altered in 13(92.86) of 14 samples( Figure 6A, C). Moreover, the horizontal bar chart showed the characteristics of the top 10 mutations for different risk types. Notably, TP53 mutations were present in 36 of patients in the high-risk group and only 20 in the low-risk group ( Figure 6B, D). We then calculated the TMB(tumor mutant load) score for each WT sample. TMB of the low-risk group was lower than that of the high-risk group (TMB=0.2463), but there was no statistical significance (P=0.11).

Composition of tumor-infiltrating immune cells in a prognostic model
The histogram showed the detailed composition of tumor infiltrating immune cells based on the algorithm of CIBERSORT( Figure 7A). CIBERSORT is a tool for deconvolution of the expression matrix for human leukocyte subtypes based on Linear Support Vector Regression. The violin plot revealed different proportions of tumor-infiltrating immune cells, depending on the two groups' risk assessment. Wilcoxon rank-sum test showed significant differences in macrophage M1 abundance and neutrophil distribution between the two groups ( Figure  7B). Moreover, the risk score was positively correlated with macrophage M1 (R = 0.24) and negatively correlated with macrophage M2 (r = -0.18) (Figure 7C, D). The Expression of nine autophagy-related gene in 293T and SK-NEP-1 Cells It has already been confirmed that NFE2L2, MAPK8 and ERBB2 may protective factors; moreover,SESN2,ATG5,VAMP,MTMR14,SH3GLB1 and ATG10 are carcinogenic factors in WT(table1). Therefore, we analyzed the expression of these ARGs in the WT cell line SK-NEP-1 and normal Human renal epithelial cell line 293T. Our results indicated that NFE2L2 and ERBB27 were lowly-expressed in SK-NEP-1. The expression of MAPK8,SESN2 and MTMR14 were not statistically different between SK-NEP-1 and 293T, while the high expression of ATG5, VAMP,SH3GLB1 and ATG10 were in SK-NEP-1 (Figure 8).

Discussion :
Autophagy is a biological process to eliminate dysfunctional organelles, aggregates or even long-lived proteins. This catabolic process appears in various biological stages from yeast to man to maintain the cellular homeostasis [6]. Autophagy mainly consists of four steps: initiation, nucleation, maturation, and degradation [7]. At present, the relationship between autophagy and tumor remains unclear. The basal level of autophagy is regarded as tumor suppression element. Meanwhile, the lack of autophagy results in tumor development attributed to the inhibition of elimination damaged components or proteins under stressful conditions [8,9]. Nephroblastoma is the most common form of childhood kidney cancer and is hereditary heterogeneous. Although modern surgery has significantly improved this disease's survival rate, a small number of patients with "high-risk" Wilms tumors still have poor prognosis. Despite modern treatments, one in ten children will die from Wilms's tumor [10,11]. Therefore, it is urgent to find a precise diagnostic and therapeutic biomarkers to provide a reference value for clinical. It's worth noting that autophagy is significantly related to Wilms tumor. There has been reported that Beclin-1 and ATG4B expression are both considered as a biomarker correlated with favourable histology. The up-regulated expression of these two genes also predict a better survival condition [12]. While autophagy dysfunctions occur in Wilms tumors, strategies targeting autophagy may provide a new approach for the treatment of highly malignant cases [13]. Considering the possible role in the autophagy, we suspect that the construction of autophagy-related prognostic signature may have far-reaching influence in early diagnosis and treatment for WT patients. First of all, we gathered the expression information of 204 ARGs from the database of TCGA. We screened 24 differentially expressed ARGs from WT and normal tissues, and put these genes into GO analysis to better understand whether the biological function of WT is related to autophagy. We filtered out 18 OS-related ARGs based on the univariate Cox regression analysis. Then the multivariate Cox regression analysis was used to screen out the ultimate prognostic genes. The risk signature was established based on 9 ARGs. All patients were divided into two groups based on the calculation of risk score, and the risk score was also related to clinicopathological parameters such as gender, T stage and first event. Furthermore, Multivariate Cox regression proved that risk characteristics could be an independent factor for the prognosis of WT patients. WT is a typical genetic illness, and the underlying biological mechanisms have not been elucidated. Although the most common mutations include WT1, WTX, and CTNNB1, they appear in a third of the population of the WT patient [14]. TP53(tumor protein P53) is frequently mutated in human diseases and has a certain degree of correlation with the occurrence and development of patients with nephroblastoma. In nephroblastoma, the mutant gene encoding TP53 is associated with intercellular variation, but the relationship between TP53 and histologically benign nephroblastoma has not been studied. We then compared the mutated genes between the low-risk group and high-risk group. The waterfall plot showed significant differences in the distribution of gene mutations in the two subgroups. The mutation rate of TP53 in the high-risk group was significantly higher than that in the low-risk group, leading us to speculate that autophagy may be involved in somatic cell mutations, but this requires more future verification. TP53 mutations are found in the anaplastic areas of WT [ 15], and immunostaining for p53 is positive in 76 of anaplastic WT compared with 8 of favorable histology [16]. These findings were consistent with the results based on our ARG-related risk signature. Immune escape and immunosuppression play an indispensable role in tumor microenvironment, and immunotherapy has been a hot research topic in recent years. We discovered the risk score positively associated with macrophages M1 and negatively related to macrophages M2 for the first time. We all know that tumor-associated macrophages play a crucial role in the tumor microenvironment. Studies have shown that with the occurrence and progression of the disease, the density of macrophage M2 increases, while the density of M1 changes in the opposite direction [18]. This is also consistent with our risk model: risk score is positively correlated with M2 macrophage density, and both predict poor prognosis. Our risk signature comprises of nine genes: ATG5, VAMP7, MAPK8, SH3GLB1, ERBB2, SESN2, MTMR14, ATG10, NFE2L2. Specifically, ERBB2, MAPK8, and NFE2L2 were protective factors, while the remaining six genes were risk factors that may contribute to tumor progression. ATG5 is one of the most common autophagy-related genes which participated in the ATG 12-ATG5-ATG16L1 complex and combined with MAP1LC3/LC3 to the phagophore membrane [19]. VAMP7 is a longin family R-SNARE protein with several biological functions such as synthesis of autophagosome and cell adhesion, and more importantly, participated in the process of autophagosome-lysosome fusion [20]. MAPK8 induces autophagy under starvation, but the underlying mechanism by which MAPK8 participates in autophagy remains unclear [21]. SH3GLB1 is a membrane curvature-inducing protein that interacts with BECN1 though UVRAG and regulates the post-Golgi trafficking of membrane-integrated ATG9A for autophagy [22]. Among all these prognosis genes, only one gene has been reported and studied in WT patients: ERBB2. The authors confirm that ERBB2 is a tumor-associated antigen in WT expressed in almost all tumor stages [23]. This finding is consistent with our previous findings that ERBB2 is a protective gene that may serve as a therapeutic target in Wilms tumors. SESN2 could induce the mitophagy in macrophages to eliminate the damaged or incapacitated mitochondria [24]. The biological function of transcription factor NFE2L2 is to maintain the cellular homeostasis, and control the expression of cytoprotective gene [25]. However, there are still some deficiencies in our study. First, we lacked the verification of autophagy-related gene function in WT patients no matter in vivo or in vitro. Second, we only used the data from the TARGRT database, and other available databases needed to be validated. Due to the lack of clinical samples, we intend to collect tissue samples and expression profiles by ourselves and improve the work mentioned above in future work. In conclusion, We constructed prognostic risk markers associated with autophagy and verified its reliability in predicting the prognosis of WT patients for the first time. This model provides new targets for WT treatment, both for preoperative evaluation and prediction of postoperative recurrence. In addition, we verified that the high-risk score was associated with p53 mutations and immune infiltration, providing a new approach for the selection of personalized treatment strategies. Further experiments are needed to clarify the relationship between WT patients and autophagy-related genes.               The Expression of nine autophagy-related gene in 293T and SK-NEP-1 Cells

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.