Identification of Core Genes Involved in the Metastasis of Clear Cell Renal Cell Carcinoma

Introduction

Renal cell carcinoma (RCC) is one of the most common malignant tumors in the world. It is estimated that approximately 74,000 new cases and 15,000 deaths of RCC may occur in the US in 2020, among which clear cell RCC (ccRCC) is the most common histological type of RCC, accounting for 85–90% of all RCC cases.^1,2 The five-year survival rate of patients with distant metastasis drops to 12%, despite that the likelihood of survival in early and localized ccRCC patients is comparatively higher.³ The established studies in the past 10 years have been used in clinical diagnosis and treatment of renal tumors.^4,5 For example, it has been demonstrated that VHL (Von Hippel–Lindau) deletion directly led to mRNA alteration and facilitated RCC development. Thus, a specific drug targeting pVHL-HIF-VEGF has been utilized in clinical practice and its efficacy was shown to be better than that of cytokine therapy.⁶ However, the clinical prognosis of advanced and metastatic ccRCC remained unsatisfactory due to chemotherapy and radiotherapy resistance. Therefore, it is of great necessity to develop specific markers that can better predict the prognosis of ccRCC, provide new insights for the development of related drugs, reduce the incidence of advanced ccRCC, and ultimately improve the survival rate.

In recent years, microarrays based on mathematical tools can be utilized to extract useful biological information from a mass of data, identifying and analyzing the association and molecular mechanisms between key genes and central signaling pathways in tumorigenesis. In addition, data mining from high-throughput platforms can be conducted to identify significant genetic or epigenetic changes in carcinogenesis, which in turn help to select promising biomarkers for cancer diagnosis and prognosis. Mounting studies were conducted to explore the molecular mechanisms involved in RCC development through mining public datasets,^7–10 whereas few have been validated by experiments.

In this study, we downloaded gene expression profiles GSE53757 and GSE84546 from the Gene Expression Omnibus (GEO) containing paired advanced-stage ccRCC and normal samples to identify the differentially expressed genes (DEGs). In addition, gene ontology (GO) analysis was used to search for pathways, and a protein–protein interaction (PPI) network was constructed by the STRING database. The hub genes were identified by Cytoscape software. Finally, in vitro experiments were performed to verify the biological effects of typical hub genes. Collectively, the identified genes associated with the metastasis of ccRCC may shed light on novel therapeutic targets and pave the way for further exploration of the underlying mechanisms in ccRCC metastasis.

Methods

Download and Conversion of Microarray Data

Two gene expression profiles (GSE53757 and GSE84546)^11,12 were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). GSE53757 dataset contained 15 metastatic ccRCC (mccRCC) samples and 15 paired normal samples and GSE84546 dataset contained three mccRCC samples and three paired normal samples (Table 1). GEO-2R online tools (http://www.ncbi.nlm.nih.gov/geo/geo-2r/) were performed to normalize and log 2 transform the matrix data of each GEO dataset and screen out DEGs in each microarray. The value of | log FC (fold change) | ≥2 with P<0.05 was considered of statistical significance. Hierarchical clustering analysis was performed to classify the data into two groups with similar expression patterns. The Venn diagram package of R software was used to make a Venn diagram and screen out DEGs for further study.

Table 1 Statistics of the Two Microarray Databases Derived from the GEO Database

Enrichment Analysis of Functional Pathways

Cytoscape software was utilized to construct unified conceptual framework through integrating biomolecular interaction network, high-throughput expression data and other molecular states.¹³ In this study, GO term enrichment¹⁴ and the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment pathway¹⁵ were conducted using the plug-in ClueGO¹⁶ in Cytoscape software to analyze DEGs. The enrichment analysis of GO terms includes biological process (BP), cell composition (CC) and molecular function (MC).

Protein–Protein Interaction (PPI) Network Construction and Hub Gene Analysis

STRING and Cytoscape software were used to construct PPI network and screen key genes of PPI network.¹⁷ In this study, STRING database (http://string-db.org) was used to conduct PPI network analysis of DEGs in ccRCC and normal renal samples. Key genes in PPI network were screened by CytoHubba and MCC (Maximum Clique Centrality) topology analysis algorithm was selected as the node sorting method.

Correlation Analysis of Hub Genes Expression and Prognosis in ccRCC

In order to verify the mRNA expression level of hub genes in ccRCC, we employed HPA database (http://www.proteinatlas.org/) to obtain IHC staining. The total survival analyses for hub genes were conducted through online cBioPortal platform (http://www.cbioportal.org). UALCAN database (http://ualcan.path.uab.edu/index.html) was used to identify survival rate for ccRCC patients. Genes significantly associated with overall survival were considered potential biomarkers of prognosis.

Cell Lines and Culture

The 786-O RCC cell line was purchased from the American Type Culture Collection (ATCC), and cultured in RPMI-1640 (Gibco, Thermo Fisher Scientific, USA) supplemented with 10% FBS and antibiotics. The cell line was maintained at 37°C in a humidified chamber supplemented with 5% CO₂.

Oligonucleotide and Plasmid Transfection

Short interfering RNA (siRNA) against KIF20A, CCNB2 and CDCA8 (si-KIF20A, si-CCNB2 and si-CDCA8), and negative control siRNA with nonspecific sequences (si-NC) were purchased from Santa Cruz Biotechnology (Santa Cruz, CA, USA). For cell transfection, an equal dose (40 pmol) of si-KIF20A, si-CCNB2, si-CDCA8, or si-NC was transiently transfected into 786-O cells, which were seeded in six-well plates. The cells were employed in in-vitro loss-of function assays 48 hafter transfection.

CCK-8 (Cell Counting Kit-8) Assay

Cells were seeded in 96-well plates with a density of 2×10⁴ cells/well. Approximately 20 µl CCK-8 reagent was added to each well, incubated at 37°C for 10 min, and then placed on a shaker for five minutes. The OD value of each well was detected. Each assay was performed in triplicate and repeated three times.

Transwell Assay

After transfection by si-KIF20A, si-CCNB2, and si-CDCA8 for 48 h, 786-O cells were diluted to a final concentration of 5.0×10⁶/mL for the invasion assay. The transfected cells were placed on the filter membrane of the upper chamber of transwell and incubated at 37°C for one hour. In the lower chamber, 600 µl of complete cell culture medium was added, and 100 µl of cell suspension was added in the upper chamber. Afterward, 786-O cells were fixed with methanol, stained with 0.2% crystal violet, counted, and photographed under a microscope.

RNA Extraction and Quantitative Real-Time PCR Analysis

Total cellular RNA isolation was carried out following the manufacturer’s protocol. The synthesis of cDNA was performed using the PrimeScript RT reagent (Takara, Kusatsu, Japan) after measuring the concentration and purity of the extracted RNA. The fluorescence quantitative PCR detection of mRNAs was conducted using SYBR Premix Ex Taq I. PCR amplification was carried out with reverse-transcribed cDNA as a template and U6 snRNA as an internal control. The experimental results were analyzed quantitatively by using the standard 2^–ΔΔCt method.

Western Blot Analysis

Total protein was extracted from 786-O cells after lysis in RIPA buffer, and the concentration was detected by the BCA method. The buffer solution was added according to the protein concentration, and SDS-polyacrylamide gels were used in electrophoresis after adding samples. Next, the protein was transferred onto polyvinylidene fluoride (PVDF) membranes, which were then blocked using a 5% skimmed milk powder solution. The membranes were then placed in a shaker at 37°C for one hour, incubated with primary antibodies overnight at 4°C, washed by TBST buffer three times, and added with secondary antibody (1:2000) for one hour. ECL liquid was used for luminescence. X-ray film exposure was carried out in the darkroom. The primary antibodies used in this experiment were as follows: anti-KIF20A (ab70791, Abcam, USA), anti-CCNB2 (ab185622, Abcam, USA), and anti-CDCA8 antibody (ab67126, Abcam, USA). GAPDH was used as an internal control.

Statistical Analysis

SPSS version 24.0 software was used for all statistical analyses. The log-rank method was applied in Kaplan–Meier analysis to determine the survival difference of hub genes. Student’s-t test and one‐way ANOVA analysis were used to compare the significance of two groups. The qPCR and Western blot results were analyzed via Wilcoxon signed-rank tests. All experiments above were repeated three times. A P-value of <0.05 was considered statistically significant.

Results

Identification of DEGs

After the standardization of GSE53757 and GSE84546 profiles, Venn analysis was performed to obtain the DEG profiles intersection (2791 in GSE84546 and 923 in GSE53757). The result showed that a total of 312 genes were markedly differentially expressed in two groups. Limma R package was used to identify 312 DEGs, among which 148 genes were upregulated and 164 downregulated. The DEGs were analyzed by hierarchical cluster analysis (Figure 1).

Figure 1 Identification of DEGs. Venn diagram of DEGs (2791 in GSE84546 and 923 in GSE53757) was used to obtain the DEG profiles intersection. Three hundred and twelve DEGs were differentially expressed in 2 groups.

Enrichment Analysis

To better understand the targeted DEGs, GO function analysis was conducted by using Cytoscape plug-in ClueGo. The results presented that the DEGs were mostly enriched in the biological processes (BP) such as blood vessel development (25.0%), tube morphogenesis (23.33%) and regulation of cell adhesion (13.33%) (P<0.01) (Figure 2). Furthermore, KEGG pathway analysis demonstrated that DEGs were mostly aggregated in the intestinal immune network for IgA production (50.0%) and bile secretion (50.0%) (P<0.05) (Figure 3, Table 2).

Table 2 Significantly Enriched GO Terms and KEGG Pathways of DEGs

Figure 2 Enrichment analysis of functional pathways. (A) GO analysis of the DEGs of metastatic ccRCC. (B) Pie chart of different GO terms of the DEGs.

Figure 3 Enrichment analysis of functional pathways. Enriched key genes involved in KEGG analysis.

PPI Network Analysis and Screening for Hub Genes

PPI analysis of DEGs was predicted by online STRING software. The processed data were imported into Cytoscape to generate PPI network (Interaction score setting >0.7, strong correlation). The results revealed 301 nodes and 215 edges among the DEGs, as shown in Figure 4A. MCC topology analysis algorithm of CytoHubba tool was conducted to screen the top 10 hub genes in PPI network, including downregulated C3, CXCR4, CCL4, ACKR3, KIF20A, CCNB2 and CDCA8. It upregulated CCL28, S1PR5 and CCL20, with the connectivity scores of the top 10 central genes listed in Table 3 (Figure 4B, Table 3).

Table 3 The Top 10 Hub Genes in PPI Network with Highest Degree

Figure 4 PPI network analysis and screening for hub genes. (A) PPI network of DEGs. (B) The PPI network of top 10 hub genes.

Analysis of Hub Genes

The overall survival rate analysis of identified hub genes was conducted by using the Kaplan–Meier survival curve with a log-rank test. In the samples of the ccRCC patients at an advanced stage, CCL4, KIF20A, CCNB2, and CDCA8 were related to worse overall survival, whereas CCL20 and CCL28 were accompanied by a favorable prognosis (Figure 5). Next, UALCAN online tools were employed to validate the associations between these four genes and the poor prognosis of RCC patients. KIF20A, CCNB2, and CDCA8 were concordantly related to unfavorable survival rates. (Figure 6). Consistently, analysis of TCGA data from the UALCAN database illustrated that KIF20A, CCNB2, and CDCA8 were overexpressed in RCC tissues compared to normal ones, and the expressions of the screened hub genes were positively related to clinical stage and lymph node metastasis. This result suggested that KIF20A, CCNB2, and CDCA8 might play a crucial role in the development and progression of ccRCC (Figure 7). After analyzing the HPA database, it was also found that KIF20A, CCNB2, and CDCA8 were upregulated in RCC tissues (Figure 8). Afterward, the association between KIF20A, CCNB2, and CDCA8 was investigated. By using the LinkedOmics database, it was discovered that CDCA8 was positively correlated with CCNB2 and KIF20A, and CCNB2 was positively correlated with KIF20A in RCC (P<0.001) (Figure 9A–C).

Figure 5 Analysis for hub genes. The overall survival rate of identified hub genes was analyzed using the Kaplan–Meier survival curve.

Figure 6 Analysis for hub genes. The overall survival analyses of hub genes were performed by cBioPortal online platform using data of 531 renal cell carcinoma patients from TCGA database.

Figure 7 Analysis for hub genes. (A) Association of the expression of CCNB2, CDCA8, and KIF20A  with normal and tumor tissues. (B) Association of the expression of CCNB2, CDCA8, and KIF20A with different clinical stages. (C) Association of the expression of CCNB2, CDCA8, and KIF20A with different nodal metastasis events. Abbreviations: N0, renal cell carcinoma without nodal metastasis; N1, renal cell carcinoma with nodal metastasis.

Figure 8 Analysis for hub genes. The protein expression levels of CCNB2, CDCA8 and KIF20A in normal and RCC tissues from HPA.

Figure 9 Analysis for hub genes. (A) The correction between CDCA8 and CCNB2 in ccRCC was analyzed by LinkedOmics. (B) The correction between CDCA8 and KIF20A in ccRCC was analyzed by LinkedOmics. (C) The correction between CCNB2 and KIF20A in ccRCC was analyzed by LinkedOmics.

Verification of Hub Gene Expressions in RCC Tissues

Considering that KIF20A, CCNB2, and CDCA8 were closely associated with ccRCC metastasis, the effects of these genes on the biological activities of RCC cells were investigated. Firstly, we validated three hub genes in normal and RRCC cells were treated with the inhibitors of CCNB2, CDCA8, and KIF20A, respectively. The knockdown efficiency after siRNA transfection was detected via qPCR (Figure 10A, C and E). Western blot analysis was also conducted to measure the expression difference of hub genes after the knockdown (Figure 10B, D and F). The results demonstrated that the expression of these genes was markedly decreased at mRNA and protein levels. Expectedly, the knockdown of CCNB2, CDCA8, and KIF20A suppressed RCC cell growth and invasion ability compared with that of si-NC-transfected cells (Figure 11A–F).

Figure 10 Analysis for hub genes. (A, C, and E) The knockdown efficiency after si-CCNB2, si-CDCA8, and si-KIF20A transfection was detected via qPCR. (B, D, and F) The knockdown efficiency after si-CCNB2, si-CDCA8, and si-KIF20A transfection was detected via Western blot. *P<0.05.

Figure 11 Verification of hub genes expression in RCC tissues. (A) CCK-8 assay was performed to determine cell viability of CCNB2 on 786-O cells. (B) Invasion assay was performed to detect invasion ability of CCNB2 on 786-O cells. (C) CCK-8 assay was performed to determine cell viability of CDCA8 on 786-O cells. (D) Invasion assay was performed to detect invasion ability of CDCA8 on 786-O cells. (E) CCK-8 assay was performed to determine cell viability of KIF20A on 786-O cells. (F) Invasion assay was performed to detect invasion ability of KIF20A on 786-O cells. Each bar represents the mean ±SD of three independent experiments. *P<0.05.

Discussion

Despite that molecularly targeted therapies, representatively, drugs targeting vascular endothelial growth factor (VEGF) pathway and mammalian/mechanistic target of rapamycin (mTOR) pathway, have played a promising role in treating metastatic RCC, the median survival time is only 13 months, and the five-year survival rate is less than 10% due to the paucity of effective molecular targets. Even worse, the vast majority of the patients with metastatic RCC have eventually developed drug resistance, leading to unfavorable survival time. It is, therefore, important to further explore the underlying mechanisms or molecular markers for the diagnosis and treatment of ccRCC.

In this study, GO enrichment analysis demonstrated that DEGs are mainly enriched in the BPs such as blood vessel development, tube morphogenesis and regulation of cell adhesion. KEGG pathway analysis revealed that DEGs are mostly aggregated in intestinal immune network regulating IgA production and bile secretion. The top 10 hub genes were identified by PPI network analysis. By searching through HPA and UALCAN databases, we found that CCNB2, CDCA8, and KIF20A were significantly overexpressed in ccRCC tissues and closely associated with the poor prognosis of ccRCC patients. Survival analyses performed by the Kaplan–Meier method demonstrated that CCNB2, CDCA8, and KIF20A production were related to the survival rate of ccRCC patients. The correction between CCNB2, CDCA8 and KIF20A in RCC was also positively correlated with each other, illustrating that they might serve pivotal roles in the progression and metastasis of ccRCC.

CCNB2 (cyclin B2) is a member of cyclin family proteins, which are important components in the cell cycle regulation mechanism. It has been fully demonstrated to trigger G2/M transition by activating CDK1 kinase.¹⁸ Accumulating evidence proved that CCNB2 was overexpressed in various human cancers, such as lung cancer, colorectal adenocarcinoma, and ovarian carcinoma.^19–21 Moreover, it has been found that advanced clinical stage and worse metastasis status were associated with the elevated expression of circulating CCNB2.²² As a biomarker, it has been reported that overexpression of CCNB2 could act as an independent predictor of poor prognosis in patients with lung adenocarcinoma and breast cancer. However, little is known about the prognostic role of CCNB2 in metastatic ccRCC. Integrating the results from survival analysis and the differential expression in different stages, we hypothesized that CCNB2 might act as a potential indicator and therapy target in mccRCC treatment.

CDCA8 is a member of the chromosomal passenger complex (CPC), which dominates cell division.²³ It is commonly acknowledged that dysregulation of cell cycle and uncontrolled cell proliferation are hallmarks of cancer, and loss of CDCA8 leads to defective cell proliferation and early fetal death. Established studies have illustrated that CDCA8 overexpression is essential for tumor development and is associated with poor prognosis in patients with cutaneous melanoma and breast cancer.^24,25 In this study, we observed that CDCA8 overexpression had the same effect in ccRCC patients, especially in stage IV patients with lymph node metastasis. Our study found that CDCA8 was upregulated in ccRCC and might act as a prognostic biomarker.

Previous studies have shown that KIF20A was essential for cell cycle mitosis, and its aberrant accumulation promoted cell proliferation.^26,27 In addition, abnormal expression of KIF20A participates in the carcinogenesis of multiple malignancies, including cell differentiation, invasion, and distant metastasis, which ultimately lead to drug resistance^28–30. It is worth noting that KIF20A can not only be used as a novel biomarker for predicting prognosis but also a new potential therapeutic target in various malignant tumors through different mechanisms. Representatively, Asahara et al conducted a single-center phase I/II clinical trial utilizing a cancer vaccine reagent KIF20A-66 to treat advanced pancreatic cancer, which significantly prolonged the overall survival rate of patients.³¹ Therefore, further study on the carcinogenic mechanism of KIF20A may contribute to the development of molecularly targeted drugs. Consistently, our results echoed with previous research showing that KIF20A might be a key molecule affecting the development of ccRCC and as a result impeding the prognosis.

ccRcc is the most common subtype of RCC, and its tumor progression and prognosis are associated with complex gene interactions. Therefore, exploring potential molecular mechanisms or markers of ccRCC is of great necessity for the diagnosis and treatment of ccRCC. However, there are a lack of molecular markers related to the pathological stage of ccRCC. Notably, few effective treatments for metastatic ccRCC have been developed, despite that the response rate to VEGFR inhibitors is relatively favorable (≈50%). Further exploration is required considering that there is a lack of available drugs to treat metastatic RCC. In this study, we identified DEGs in metastatic ccRCC samples by using bioinformatics analysis, and screened out three hub genes CCNB2, CDCA8, and KIF20A associated with mccRCC. These genes can be used as personalized therapeutic targets and diagnostic biomarkers for metastatic ccRCC. Our study may contribute to a better understanding of ccRCC and have important clinical significance in improving risk stratification, therapeutic strategy, and prognosis prediction of ccRCC patients.