Bioinformatical analysis of the key differentially expressed genes for screening potential biomarkers in Wilms tumor

Wilms tumor (WT) is the most common pediatric renal malignant tumor in the world. Overall, the prognosis of Wilms tumor is very good. However, the prognosis of patients with anaplastic tumor histology or disease relapse is still poor, and their recurrence rate, metastasis rate and mortality are significantly increased compared with others. Currently, the combination of histopathological examination and molecular biology is essential to predict prognosis and guide the treatment. However, the molecular mechanism has not been well studied. Genetic profiling may be helpful in some way. Hence, we sought to identify novel promising biomarkers of WT by integrating bioinformatics analysis and to identify genes associated with the pathogenesis of WT. In the presented study, the NCBI Gene Expression Omnibus was used to download two datasets of gene expression profiles related to WT patients for the purpose of detecting overlapped differentially expressed genes (DEGs). The DEGs were then uploaded to DAVID database for enrichment analysis. In addition, the functional interactions between proteins were evaluated by simulating the protein–protein interaction (PPI) network of DEGs. The impact of selected hub genes on survival in WT patients was analyzed by using the online tool R2: Genomics Analysis and Visualization Platform. The correlation between gene expression and the degree of immune infiltration was assessed by the Estimation of Stromal and Immune cells in Malignant Tumor tissues using the Expression (ESTIMATE) algorithm and the single sample GSEA. Top 12 genes were identified for further study after constructing a PPI network and screening hub gene modules. Kinesin family member 2C (KIF2C) was identified as the most significant gene predicting the overall survival of WT patients. The expression of KIF2C in WT was further verified by quantitative real-time polymerase chain reaction and immunohistochemistry. Furthermore, we found that KIF2C was significantly correlated with immune cell infiltration in WT. Our present study demonstrated that altered expression of KIF2C may be involved in WT and serve as a potential prognostic biomarker for WT patients.

qRT-PCR experiment.qRT-PCR was used to measure KIF2C mRNA expression.Total RNA was extracted from 3 randomly selected WT tissues and 3 adjacent WT tissues with Tissue Total RNA kit (5100050, Simgen, China) according to the protocol.Subsequently, reverse transcription was then performed on the extracted RNA to produce cDNA.Real-time PCR was performed using a StepOne Plus Real-Time PCR System, using an Applied Biosystems Prime Script tm RT Master Mix, containing 20 ng cDNA and 10 μL of each primer.Cycling conditions were 1 cycle at 95 °C for 30 s, 40 cycles at 95 °C for 5 s, and 60 °C for 34 s.PCR amplification specificity was determined by melting curve analysis for each reaction.The expression level of KIF2C was normalized to GAPDH and calculated using 2 −∆∆Ct method.The primer sequences for KIF2C were forward (GAA GAG AGC TCA GGA GTA TG) and reverse (TTC TGT GCT CTT CGA TAG GA).The GAPDH primer sequences were forward (GGG GCT CTC CAG AAC ATC ) and reverse (TGA CAC GTT GGC AGTGG).
Immunohistochemical staining.Seventeen samples for IHC were collected from WT tissues removed from post-operative patients in the Department of Surgical Oncology.Three of the samples contained both tumor tissue and adjacent normal tissue.Paraffin-embedded histological slices were cut into 5 μm slices after 10% formalin fixation.Sections were stained with hematoxylin and eosin (H&E) or deparaffinized with dimethylbenzene and dehydrated through 100, 95, 85, and 75% ethanol.The antigen retrieval procedure was performed in 0.01 mol/L sodium citrate buffer by heating in a pressure cooker for 10 min and cooling naturally.Endogenous peroxidase was blocked by incubation with 3% H 2 O 2 for 15 min at room temperature.The sections were subsequently washed in triethanolamine buffered saline with Tween (TBST) solution, blocked with bovine serum albumin (BSA) for 20 min and incubated with anti-KIF2C (28372-1-AP, 1:200, Proteintech, China) overnight at 4 ℃.Following three TBST solution washes, the sections were incubated with HRP-conjugated secondary antibody for 20 min at room temperature.All slides were counterstained with diaminobenzidine (DAB) solution (Sangon Biotech, Shanghai, China).Finally, sections were counterstained with hematoxylin, dehydrated, cleared and mounted.The expression level of KIF2C was determined by immunohistochemical score (IHS).IHS were determined based on staining intensity (SI) and percentage of immunoreactive cells (PR).Staining intensity (SI): 0: no neoplastic cells showed membranous-like staining; 1: ≤ 10% of neoplastic cells showed incomplete, weak circumferential membranous-like staining; 2: ≤ 10% of neoplastic cells showed complete, weak circumferential membranous-like staining.3: > 10% of neoplastic cells showed incomplete circumferential membranous-like staining, 4: > 10% of neoplastic cells showed complete circumferential membranous-like staining; Neoplastic cells with different staining intensities were assigned different scores, and the sum was calculated.

Immune infiltration analysis. Estimation of Stromal and Immune cells in Malignant
Tumor Tissues using the Expression date (ESTIMATE) algorithm was used to calculate ESTIMATE score, Immune score and Stromal score by the ESTIMATE 1.0.13R package.Based on these scores, we assessed the difference between tumor microenvironment exhibited by high and low KIF2C expression groups.Furthermore, the microenvironment cell population counter (MCP-counter) 1.2.0 R package was used to assess the composition of immune cells in the high and low expression groups of KIF2C.A further analysis of the enrichment scores of 28 immune-related gene sets was carried out by using ssGSEA in the high and low KIF2C expression groups.

Statistical analysis.
All data were analyzed using GraphPad Prism 8.0 statistical software and R programming language (version 4.0.2).Differences between groups were compared using t-test.Statistical significance was determined by a P value of < 0.05.

Ethics approval and consent to participate.
Informed consent was obtained for tissue samples used for qRT-PCR and IHC.This study has been reviewed and approved by the Ethics Committee of the Children's Hospital of Zhejiang University School of Medicine (2020-IRB-049-A1).All methods were carried out in accordance with relevant guidelines and regulations.

Consent for publication.
Written informed consent was obtained from all legally authorized open-source data for the publication of identifiable data.

Result
Screening of differentially expressed genes.A total of 1598 and 965 DEGs compared between WT and normal kidney tissue were identified from the GSE66405 and GSE11024 datasets, respectively.Based on the two sets of sample data, the differential expression of multiple genes across the two microarrays is shown as volcano plots (Fig. 1A).In total, 506 genes between the two datasets were shown to be identical based on bioinformatics analysis (Fig. 1B).Compared with normal samples, 142 genes were up-regulated and 364 genes were down-regulated in WT samples (Table 1).Heat maps of the top 50 genes up-regulated and down-regulated expression are shown in Fig. 1C.

GO and KEGG pathway enrichment analysis.
To further understand the diversity of functions performed by DEGs, we conducted GO functions and KEGG pathway enrichment analyses for DEGs on the DAVID platform.By extending analysis of the enriched biological process (BP), cellular components (CC) and molecular function (MF), we gained a deeper understanding of the biological functions of overlapping DEGs.According to P values, the ggplot2 R package was used to construct the bubble map of the top 10 biological processes (Fig. 2).The up-regulated genes were mainly related to the cell division, sister chromatid cohesion, mitotic nuclear division in enriched GO terms for BP; meanwhile nucleus, spindle, and chromosome centromeric region for CC; and protein binding and microtubule binding for MF.The down-regulated genes were mostly related to excretion, long transmembrane transport in enriched GO terms for BP; extracellular exosome for CC; and transporter activity and drug transmembrane transporter activity for MF (Fig. 2).Furthermore, in the KEGG pathway enrichment analysis, there were four main pathways that were associated with up-regulated DEGs, including the cell cycle, HTLV-I infection, oocyte meiosis, and focal adhesion, while almost only one metabolic pathway was enriched in the down-regulated genes (Fig. 3).

PPI network construction and Hub gene selection.
The PPI network was constructed by selecting 458 nodes and 2638 edges consisting of 142 up-regulated genes and 364 down-regulated genes (Fig. 4A).The top 12 hub genes including KIF2C, centromere protein E (CENPE), kinesin family member 20A (KIF20A), thymidylate synthetase (TYMS), topoisomerase II alpha (TOP2A), kinesin family member 15 (KIF15), ubiquitin-conjugating enzyme E2C (UBE2C), cyclin B2 (CCNB2), TPX2 Microtubule Nucleation Factor (TPX2), enhancer of zeste 2 polycomb repressive complex 2 subunit (EZH2), cell division cycle 20 (CDC20), and centromere protein F (CENPF) were identified by the cytoHubba plug-in (Fig. 4B).All the 12 hub genes identified in the PPI network were up-regulated in WT compared to normal kidney samples.The MCODE plug-in was used to determine the most significant module from the PPI network based on the degree of importance, which included 48 nodes and 1048 edges (Fig. 4C).An analysis of gene enrichment revealed that genes in this module are mainly involved in the cell cycle, oocyte meiosis, and HTLV-1 infection (Fig. 4D).

Kaplan-Meier survival plotter.
The prognostic information related to the 12 hub genes can be found on the platform R2: Genomics Analysis and Visualization.High expression level of KIF2C (p = 3.2e-03), KIF15 (p = 7.6e-03), KIF20A (p = 9.0e-03), CENPF (p = 0.031) and CENPE (p = 0.035) was found to be associated with worse OS of WT patient (Fig. 5).High expression level of TOP2A (p = 0.017) was related to better performance at patients' overall survival.Among these genes, KIF2C played the most important role in prognosis.Thus, KIF2C was chosen for further analysis.www.nature.com/scientificreports/qRT-PCR to verify the expression of KIF2C.The expression of KIF2C transcriptional mRNA level was verified in tissues using qRT-PCR in 3 available WT patient samples.Compared with 3 adjacent WT normal tissues, there was a significant increase in the expression level of KIF2C in WT tissues (p = 0.0065) (Fig. 6).

Immunohistochemistry to verify the expression of KIF2C.
Immunohistochemical staining demonstrated the expression of KIF2C in blastemal, epithelial and stromal cells of WT tissues as well as in normal kidney cells.Compared with normal kidney tissues, the positive immunohistochemistry staining of WT tissues was mainly manifested as membranous-like staining, and most of neoplastic cells showed incomplete, weak circumferential membranous-like staining (Fig. 7A).The expression socre of KIF2C in WT tissues was significantly higher (p = 0.037) than that of adjacent normal tissues (Fig. 7B).

The correlaion between KIF2C expression and immune infiltration.
In order to assess the relationship between KIF2C expression levels and tumor microenvironment, we first compared the differences in ESTI-MATE, Immune and Stromal score among the high and low KIF2C expression groups.As compared to the low expression group, the high KIF2C expression group's ESTIMATE score, immune score and stromal score were all significantly lower (Fig. 8A).Furthermore, the MCP-counter algorithm and ssGSEA analyses were performed identify the difference in infiltrating immune cells in the tumor microenvironment of WT patients.The results of MCP-counter showed a decreased number of tumor-infiltrating immune cells, including monocytic lineage, myeloid dendritic cells and T cells in the high KIF2C expression group compared with the low expression group (Fig. 8B).Consistently, ssGSEA results also showed that the vast majority of infiltrating immune cells were fewer in the high-expression group than low-expression groups, except that Type 2T helper cell were more in the highexpression group (Fig. 8C).

Discussion
Despite advances in the surgery and treatment of WT, the prognosis of patients with anaplastic tumor histology or disease relapse is still poor.Therefore, it is important to find reliable tumor biomarkers and investigate possible molecular mechanisms to understand the pathogenesis and abnormal biological behavior of WT, as well as to identify novel therapeutic targets of WT.Although numerous studies have demonstrated that genes contribute significantly to WT development and incidence, most of them were limited to cohort studies 39,40 .In this study, we identified genes with significant differences in two GEO databases (GSE66405 and GSE11024 datasets), and finally found 506 common altered DEGs through bioinformatic method, including 142 up-regulated genes and 364 down-regulated genes.The DAVID database was used to analyze GO terms and KEGG pathways for these 506 DEGs.To screen out the hub genes closely related to WT, a PPI network was constructed.The top 12 hub genes with higher degree of connectivity were filtered from the PPI network.The gene hub genes were primarily associated with the cell cycle, meiosis of the oocyte, and HTLV-1 infection according to functional and KEGG pathway enrichment analysis.We plotted survival curves for each of these hub genes using the Kaplan-Meier method and found that KIF2C had the most predictive value for WT survival.Finally, we verified the expression of KIF2C in WT tissues by qRT-PCR and IHC respectively.It is worth noting that the IHC results showed  www.nature.com/scientificreports/an expression of KI2C primarily in circumferential membranous-like areas.Most neoplastic cells showed weak incomplete membranous-like staining, and some neoplastic cells showed complete membranous-like staining.It is consistent with the results of the bioinformatic analysis we conducted.According to our analysis of GSE66405 and GSE11024, the log fold change values for KIF2C are 2.5 and 1.907768, respectively.This suggests that the expression of KIF2C in WT is different but not very strong compared with normal kidney tissues.Therefore, we used a scoring standard similar to that used for HER2 in breast cancer when assessing IHC results 41,42 .Our results show that the expression score of KIF2C in WT tissues was significantly higher (p = 0.037) than that of adjacent normal tissues.The expression scores of blastemal cells appear to be higher, which may be associated with a poor prognosis 21 .
A major finding of our study concerned the KIF2C gene.Kinesins are the microtubule-associated motors that produce mechanical work associated with adenosine-triphosphate (ATP) hydrolysis 43 .Kinesin family-13 is the critical regulator of microtubule dynamics during mitosis 44 .KIF2C, a mitotic centromere-associated kinesin, is the most representative member of the Kinesin family-13, which is located in the cytoplasm throughout the cell cycle 45 .According to previous studies, KIF2C regulates mitosis and cell cycle by attaching microtubules to kinetochores, assembling spindles, and regulating chromosome congression [20][21][22]46 . In ddition, KIF2C participates in cytoskeleton remodeling during tumor metastasis 47 .It has been reported that KIF2C can stimulate the proliferation and migration of gastric cancer cells as well as non-small cell lung cancer cells 25,48 .It would therefore be worthwhile to further research the role of KIF2C in WT.
The function of KIF2C is critical for the assembly of the spindle to repair microtubules and chromosomal abnormalities.Therefore, abnormal expression of KIF2C could theoretically contribute to the development of cancer by causing mitotic anomalies, chromosomal instability, and uncontrolled transcription 49 .A number of studies have demonstrated abnormal expression of KIF2C in cancers such as breast cancer, oral tongue cancer, hepatocellular carcinoma, and non-small cell lung cancer. 22,26,50,51.In lung adenocarcinoma, Bai et al. demonstrated that KIF2C expression is closely associated with the relapse of the cancer and the stage of the tumor. 49.Nakamura et al. identified that KIF2C gene expression in gastric cancer tissues (65 cases) was significantly higher than the expression in non-malignant tissues, and the higher expression of KIF2C was associated with observed a significantly increased expression of KIF2C in malignant compared with autologous healthy colorectal tissues 53 .However, the KIF2C expression level in WT remains unclear.According to our study, bioinformatic analysis indicates that KIF2C expression is significantly higher in WT tissues.The results of qRT-PCR and IHC staining validated our bioinformatic analysis at mRNA and protein levels, respectively.These results, combined with our findings, suggest that KIF2C may act as a novel biomarker that may be associated with prognosis in patients with WT.
Recent studies have explored how the KIFs family regulates chromosome and spindle movements during mitosis and meiosis 54,55 .However, the mechanism of KIF2C dysregulation in WT remains largely unknown.It www.nature.com/scientificreports/has been demonstrated that KIF2C seems to be related to a poor prognosis by regulating cell cycle signaling pathway.An et al. illustrated that KIF2C promotes endometrial cancer cell proliferation, migration, and invasion in vitro.In addition, they found that cells were prevented from entering the G2 phase when KIF2C expression was reduced, suggesting that KIF2C knockdown may induce cell cycle arrest by prolonging the G1 phases 27 .Li et al. found that down-regulation of KIF2C effectively inhibited the growth of liver cancer cells by decreasing the proliferation and increasing the G1 arrest 51 .Based on Gan et al. 's findings, KIF2C knockdown cells showed significantly reduced viability and colony formation 26 .It has therefore been suggested that KIF2C may contribute to the proliferation of malignancies.Additionally, the results of GO and KEGG analyses also revealed that some genes closely related to cell mitosis may play important roles in WT.The hub genes identified in this article were also found to be related to the regulation of cell cycle.The findings mentioned above further suggest that KIF2C may be involved in during the cell cycle in WT.
In mitosis, KIF2C regulates microtubule dynamics and ensures that chromosomes attach correctly to spindles 28 .Considering KIF2C's role in regulating microtubule dynamics, it has the potential to be developed as a drug target.Microtubule-targeting agents disrupt microtubule dynamics, resulting in prolonged mitotic arrest  and eventually cell death 56 .Hedrick et al. reported that using anti-microtubule drugs such as paclitaxel and vinblastine, the knockdown of KIF2C in HeLa cells significantly exacerbated the morphological defects of the microtubule cytoskeleton 57 .It has been hypothesized that KIF2C silencing synergizes and enhances the effects of paclitaxel on cancer cells by interfering with spindle separation, which prolongs mitotic time, resulting in mitosis and apoptotic disorders 58,59 .It is believed that blocking mitotic exit could be a promising anticancer strategy, particularly for treating cancers such as breast cancer that are resistant to chemotherapy 60,61 .Studies by Wei et al. demonstrated the role of KIF2C in the positive regulation of mTORC1 signaling in hepatocellular carcinoma, indicating that an mTOR inhibitor could be a potential therapeutic intervention in hepatocellular carcinoma 62 .
The infiltrating immune cells are a crucial aspect of the tumor microenvironment and have a profound impact on patient prognosis and the efficacy of immunotherapy 63 .It has been demonstrated by Tu et al. that increased expression of KIF2C correlates with altered levels of immune cell infiltration in gliomas 64 .In hepatocellular carcinoma tissues, Li et al. observed that all T cell infiltration level decreased significantly compared to non-tumour tissues, indicating that low levels of T cell infiltration may represent a poor clinical prognosis 27 .However, the relationship between KIF2C expression and immune infiltration in WT is unknown.In this study, our analysis suggests that KIF2C is significantly negatively correlated with most immune cells in WT.First, we found that the proportion of infiltrating immune cells, including monocytic lineage, myeloid dendritic cells and T cells, was reduced in the high KIF2C expression group compared with the low expression group.Secondly, the results of ssGSEA analysis showed that the vast majority of infiltrating immune cells were fewer in the high-expression group than low-expression groups.Those findings are in agreement with those of the published studies.Mardanpour et al.'s research has shown that patients with high CD8 + tumor-infiltrating lymphocyte levels are associated with favorable clinical outcomes 65 .Tian et al. found that riskscore and immune score were negatively correlated in WT patients.Their research suggested that impaired anti-tumor immunity may be responsible for the poor prognosis of WT patients at high risk.Research into the mechanisms that disrupt the anti-tumor immune response will contribute to the development of more effective treatments 66 .
Our study has some limitations.Firstly, the reason why multiple datasets were selected for integrated analysis in this research is to obtain DEGs with higher specificity.However, due to experimental variations between the different datasets, the analysis based on public datasets may lack some reliability.The experiment type of two www.nature.com/scientificreports/datasets selected are both expression profiling by array.We used an online tool named GEO2R to process and normalize gene expression, so that the two datasets can be effectively compared under the same conditions.Although we have validated in PCR and IHC with a small number of tissue samples, there is still a lack of relevant cell experiments, large clinical sample experiments and reliable corresponding clinical data.It is important to note that although our research identified some mechanisms that may contribute to WT genesis and prognosis, they must be further investigated in cell experiments.In brief, we utilized bioinformatics methods to integrate two WT datasets in the GEO database to identify hub genes involved in WT development.The hub genes were all found to be up-regulated and might serve as new biomarkers for predicting the prognosis of WT.Survival analysis was used to predict the prognosis of these hub genes.KIF2C was the gene most closely associated with WT survival.We found that KIF2C was up-regulated in WT patients.There is a considerable correlation between the high expression of KIF2C and the poor OS of WT patients.The expression of KIF2C in WT was further verified by qRT-PCR and IHC.The expression of KIF2C may be associated with regulating mitosis and cell cycle.Further investigation is needed to explore these mechanisms.The biological behavior, molecular mechanism and the role of KIF2C in WT immune cell infiltration need to be further examined and characterized to obtain more accurate correlation results.Further analysis will be carried out to determine which specific immune cells play a significant role in the WT tumor microenvironment and which functions these infiltrating immune cells carry out in order to better understand the relationship between high expression of KIF2C and poor prognosis in WT patients.This study may provide new insights into the diagnosis and effective therapeutic strategies for WT in the future.

Figure 1 .
Figure 1.The identification of DEGs compared between WT and normal kidney tissue in gene expression profiling datasets.(A) Volcano plots of differential gene expression.Red indicates up-regulated DEGs; green indicates down-regulated DEGs.(B) Venn plot for DEGs that were shared between the datasets.(C) Hierarchical cluster analysis (heat maps) of the top 50 up-regulated and down-regulated common DEGs of the two datasets.DEGs differentially expressed genes.

Figure 2 .
Figure 2. Top 10 functional analyses of the overlapping up-regulated and down-regulated DEGs according to biological process, cellular components and molecular function.DEGs differentially expressed genes.

Figure 3 .
Figure 3. Top 10 KEGG pathway enrichment analyses of the overlapping up-regulated and down-regulated DEGs.KEGG Kyoto Encyclopedia of Genes and Genomes, DEGs differentially expressed genes.

Figure 4 .
Figure 4. Construction of protein-protein interaction networks of DEGs by STRING.(A) Differential gene protein interaction network.Red indicates up-regulated DEGs; Green indicates down-regulated DEGs.(B) Top 12 hub genes with a higher degree of connectivity.(C) Top module screened from the PPI network by MCODE plug-in.(D) KEGG pathway enrichment analysis of the top module.DEGs differentially expressed genes, STRING Search Tool for Retrieval of Interacting Genes, KEGG Kyoto Encyclopedia of Genes and Genomes, MCODE Molecular Complex Detection.

Figure 6 .
Figure 6.The KIF2C mRNA level was up-regulated in WT tissues by quantitative real-time polymerase chain reaction.

Figure 7 .
Figure 7. Upregulation of KIF2C expression in WT. (A) Representative immunohistochemical staining results showed the protein level expression of KIF2C in in blastemal, epithelial and stromal cells of WT tissues as well as in normal kidney tissue cells.(B) Expression scores between WT and normal tissues.↑, neoplastic cells showed incomplete, weak circumferential membranous-like staining.⇡, neoplastic cells showed complete circumferential membranous-like staining.WT Wilms tumor, HE hematoxylin-eosin staining, KIF2C Kinesin family member 2C.

Figure 8 .
Figure 8. Correlation between KIF2C expression and immune cell infiltration in Wilms tumor.(A) ESTIMATE scores, Immune scores and Stromal scores between high KIF2C group and low KIF2C group, respectively.(B, C) The abundance of various cells in the tumor microenvironment analyzed by MCP-counter algorithm and ssGSEA.*P < 0.05, **P < 0.01, ***P < 0.001.

Table 1 .
List of consistent DEGs in two datasets, including 142 up-regulated and 364 down-regulated DEGs.DEGs, differentially expressed genes.