Microarray datasets GSE85258 and GSE105288 in the Gene Expression Omnibus (GEO) database were used to analyze the transcriptomic signature of metastatic ccRCC. The GSE85258 dataset contained 15 primary tumors and 16 metastatic tumor samples. The GSE105288 dataset contained nine primary tumors and 26 metastatic tumor samples. By GEO2R tool analysis, genes with fold change (FC) > 1.5 and p < 0.05 in metastatic ccRCC compared with primary tumors were defined as metastasis-related genes. Next, we obtained gene expression profiles of 530 ccRCCs from the TCGA data portal and analyzed them by the limma package. The screening criteria for differentially expressed genes (DEGs) were p < 0.05 and FC > 2. The Ggplot2 and pheatmap packages were used for volcano plots and heatmaps generation.

The R package “glmnet” performed LASSO regression analysis to identify key prognostic genes. One standard error of the optimal penalty parameter λ value after cross-validation was used to build the simplest gene feature model. The risk score was determined by calculating the sum of the product of each gene expression level and its corresponding coefficient according to the risk score formula. Next, patients were stratified into two groups based on the median risk score. Differences in survival time between the two groups were estimated by the Kaplan-Meier method. Time-related receiver operating characteristic (ROC) analysis was used to verify the prediction accuracy of the model.

CIBERSORT analysis tool was used to infers the proportion of infiltrating immune cells based on cell-type-specific gene expression profiles from the mixed cell data. Inferred results were considered accurate at p < 0.05.

ClusterProfiler package performed Gene Ontology (GO) and the Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis (Yu et al., 2012). The PPI network was retrieved by the STRING database and reconstructed by Cytoscape software (Shannon et al., 2003; Szklarczyk et al., 2015). Cytohubba was used to analysis and select top 10 hub genes.

The C2 kegg gene set was used for GSEA enrichment analysis and FDR values <0.01 were considered significantly enriched.

The ccRCC tissues and normal samples of 10 ccRCC patients involved in this study were obtained from Peking University Third Hospital after approval by the Ethics Committee of Peking University Third Hospital. The ccRCC cell lines A498, caki-1, caki-2, 786-O and OS-RC2 and the immortalized human kidney cell line HK-2 were purchased from the American Type Culture Collection. 1640 medium (Biological Industries) containing 10% fetal bovine serum (FBS) (Biological Industries) was used to culture cells.

Cells were transfected with specific siRNAs targeting KIF23 (si-KIF23). si-KIF23 were purchased from GenePharma Biotech (Shanghai, China). The siRNA sequence is:

si-KIF23-1 sense: 5′-GGU​CCC​AAA​CGA​ACC​UUA​ATT-3′, antisense: 5′-UUA​AGG​UUC​GUU​UGG​GAC​CTT-3′;

Lipofectamine 2000 transfection reagent (Invitrogen, Carlsbad, CA, United States) was used to transfect cells following the manufacturer’s instructions.

2 × 10³ ccRCC cells were cultured in a 96-well plate, and cell proliferation changes were assessed by a Cell Counting Kit 8 (CCK-8) (Beyotime, Shanghai, China) at 24, 48, and 72 h after si-KIF23 transfection. Ethynyl deoxyuridine (EdU) kit (Ribobio, Guangzhou, China) was used to perform EdU assay following the manufacturer’s instructions.

100 μl of serum-free medium with 2 × 10⁴ ccRCC cells was added into the upper chamber of a transwell migration chamber (8-μm pore size, Costar, New York). The lower chamber was filled with medium containing 20% FBS. After 24 h, cells were fixed with 4% PFA and stained with crystal violet. Cells were counted under a microscope in five random fields after removal of cells on the upper surface. For invasion assay, upper chamber was coated with Matrigel (354480, Corning) before seeded cells.

The cells with different treatments were seeded in six-well plates, starved for 24 h, and wounded with a pipette tip. Then replace with fresh medium. After 24 h, the wound closure process was observed and photographed.

RIPA buffer (Beyotime) was used for cell and tissue lysis. 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) was used for protein separation. Then, proteins were transferred onto polyvinylidene difluoride (PVDF) membrane (Millipore, Billerica, MA, United States) by electroblotting for 90 min. After blocked in 5% non-fat dry milk for 2 h, the PVDF membranes incubated with primary antibodies at 4°C overnight, followed by the secondary antibody for at least 1 h. The blots were then visualized using a Western Blotting Luminol Reagent. The primary antibodies of KIF23 (1:1000, Proteintech, 28587-1-AP), β-catenin (1:1000, Cell Signaling Technology, 8480), GSK-3β (1:1000, Cell Signaling Technology, 12456), p-GSK-3β (1:1000, Cell Signaling Technology, 5558), c-Myc (1:1000, Cell Signaling Technology, 18583) and GAPDH (1:2000, Proteintech, 10494-1-AP). Cytoplasmic protein extraction kit (Beyotime) was used to extract nuclear proteins following the manufacturer’s manual.

ccRCC tissue sections were deparaffinized and the antigen was retrieved, followed by blocking with 5% bovine serum albumin. Next, ccRCC tissue sections were incubated with KIF23 primary antibody (1:100, Proteintech, 28587-1-AP) at 4°C overnight, a peroxidase-conjugated secondary antibody was used to detect antigen-antibody complexes. Thereafter, a color reaction was conducted using a 3,3′-Diaminobenzidine (DAB) substrate kit (ZsBio). For immunofluorescence staining, ccRCC cells were seeded on coverslips and fixed in 4% PFA. Then, cells were washed and treated with 0.25% Triton X-100 for 15 min. After blocked by 5% donkey serum at room temperature for 30 min, cells were incubated with β-catenin (1:100, Cell Signaling Technology, 8480) at 4°C overnight. Alexa Fluor 594 conjugated secondary antibody (1:500; Yeasen) was used to detect antigen-antibody complexes. 4′,6-Diamidino-2-phenylindole (DAPI) was used for the staining of nucleus.

R (version 3.6) and GraphPad Prism software were used to execute statistical analysis. Differences between groups were statistically compared using the student’s t test or two-way ANOVA. Each experiment was performed in triplicate, and values were expressed as mean ± standard deviation. A p-value <0.05 was regarded statistically significant. Survival curves were obtained by the Kaplan-Meier method and compared by log-rank tests. And levels of significance were presented at *p < 0.05, **p < 0.01, ***p < 0.001 and ****p < 0.0001 respectively.

DEGs between metastatic and primary tumors were analyzed by the GEO2R tool with FC > 1.5 and p < 0.05 as screening criteria. 693 genes (473 up and 220 down) were identified as DEGs in GSE85258.497 genes (307 up and 190 down) were identified as DEGs in GSE105288 (Figures 1A–D). Venn diagrams were used to integrate GSE85258 and GSE105288 datasets, and 74 DEGs were identified as metastasis-related genes (Figure 2A). To validate these metastasis-related genes, the ccRCC dataset from TCGA was downloaded and analyzed. The expression levels of these DEGs were visualized using a heatmap (Figure 2B). Univariate Cox proportional hazards regression analysis revealed 54 genes significantly associated with OS in ccRCC patients (Figure 2C). Forty seven genes with a hazard ratio (HR) greater than one were identified as “high risk factors”, and seven genes with HR < 1 were considered protective factors.

It is difficult for a single metastasis-associated gene to accurately predict patient prognosis. Thus, development of a multi-gene combined prediction model is essential. To construct the prognostic model, the metastasis-associated genes were used to LASSO Cox regression analysis, and the regression coefficients were calculated as shown in Figure 3A. The best performing prognostic model was obtained when 12 genes were included (Figure 3B). The risk score for each patient is calculated from the expression levels of 12 genes: risk score = (−0.0072 * PDGFRL) + (0.0012 * C1R) + (0.1278 * TRIP13) + (0.0004 * TIMP1) + (0.0016 * C1S) + (−0.0072 * NAPSA) + (0.0011 * SERPINF1) + (−0.06808 * CGN) + (0.07435 * ASPM) + (0.02810 * CEP55) + (0.00039 * COL5A1) + (0.01375 * SCG5). Then, patients were divided into low-risk (n = 263) and high-risk (n = 262) groups based on their risk scores (Figure 3C). Patients in high-risk group have a worse prognosis (Figure 3D). The area under curve (AUC) value for metastasis-associated prognostic gene signature was 0.7 and 0.661 for 3- and 5-year survival rates, respectively (Figures 3E,F).

In this study, CIBERSORT algorithm was used to estimate the proportion of infiltrating immune cells in tumor tissue. Immune cell infiltration signatures of the low-risk and the high-risk groups are showed in Figures 4A,B. There were higher proportions of resting memory CD4 T cells, activated natural killer (NK) cells, monocytes, resting dendritic cells (DCs) and resting mast cells in the low-risk group. Plasma cells, regulatory T cells, and M0 macrophages were higher in the high-risk group (Figure 4C).

Since the two groups showed variations in gene expression profiles and prognosis, gene expression patterns between the groups were analyzed and 588 DEGs were identified (408 and 181 genes were up-regulated in the high- and low-risk groups, respectively, Figures 5A,B).

GO enrichment analysis showed that DEGs were involved in “extracellular matrix organization”, “extracellular structure organization” and “positive regulation of secretion by cell” in the “biological process” category. “Collagen-containing extracellular matrix” was enriched in the “cellular components” category and “extracellular matrix structural constituent” was significantly enriched in the “molecular function” category. KEGG pathway enrichment analysis results indicated that “cytokine-cytokine receptor interaction” and “NF-kappa B signaling pathway” are significantly enriched (Figures 5C,D). Next, we constructed PPI network, and cytoHubba was employed to select the top 10 hub nodes, including PLK1, DLGAP5, KIF23, BIRC5, AURKB, CDCA8, RRM2, BUB1B, UBE2C, and CDC20 (Figure 6). KIF23, a member of kinesin proteins, and its role in tumor development has been demonstrated in many tumors (Gao et al., 2020; Hu et al., 2020), yet its role in ccRCC remains elusive. Therefore, we explored the expression and function of KIF23 in ccRCC.

The TCGA database showed that KIF23 mRNA was significantly elevated in ccRCC tissues (Figure 7A). And patients with high expression of KIF23 showed worse OS rates (Figure 7B). Univariate analysis showed no significant difference in survival based on gender. KIF23 expression level, age, tumor stage and grade were associated with OS in ccRCC patients (Supplementary Table S1). Multivariate Cox risk regression analysis showed that KIF23 was an independent factor of OS (HR = 1.919, 95% confidence interval, 1.439–2.560, p < 0.001) (Figure 7C, Supplementary Table S2). Moreover, results from GSEA analysis showed that cell adhesion molecules (CAMs), DNA replication, cell cycle and Wnt signaling pathway were significantly enriched in the KIF23 high-expression group (Figure 7D). As anticipated, the protein level of KIF23 was also increased in ccRCC tumor tissues (Figure 8A). The protein expression level of KIF23 was detected in HK-2 cells and five human ccRCC cell lines. As showed in Figure 8B, KIF23 was increased in ccRCC cell lines compared with HK-2 cells. Moreover, the results of immunohistochemistry staining showed that KIF23 expression were higher in tumor specimens (Figure 8C).

We further explored the role of KIF23 on the proliferation, migration and invasion of ccRCC cells. Silencing of KIF23 expression using siRNA resulted in down-regulation of KIF23 protein in 786-O and A498 (Figure 9A). Compared with the negative control (NC) group, knockdown of KIF23 expression inhibited the proliferation of ccRCC cells as measured by CCK-8 and EdU assays (Figures 9B,C). The migration and invasion ability were further measured using wound-healing assay and transwell assays. The results showed that interfering with KIF23 expression reduced the migration and invasion ability of ccRCC cells (Figures 9D,E). To further elucidate the mechanisms underlying KIF23-promoted ccRCC proliferation and migration, GSEA analysis was used to predict the downstream pathways, and Wnt signaling pathways were enriched in the KIF23 high-expression group (Figure 7D). We found that KIF23 knockdown decreased the expression levels of molecules on the Wnt signaling pathway include p-GSK3-β, β-catenin and c-Myc (Figure 10A). Moreover, analysis of subcellular protein fraction showed that knockdown of KIF23 reduced nuclear translocation of β-catenin (Figure 10B). Immunofluorescence staining was used to assess β-catenin translocation, and it was found that knockdown of KIF23 expression reduced nuclear localization of β-catenin in 786-O and A498 cells (Figures 10C,D).