Identification of potential novel differentially-expressed genes and their role in invasion and migration in renal cell carcinoma

Identification of DEGs in renal cell carcinoma

In this study, 39 ccRCC samples and 39 paracancerous normal kidney tissue samples were analyzed. Using the limma software package, 4224, 3879 and 549 DEGs were extracted from GSE66272, GSE100666 and GSE105261, respectively. Figure 1A–1C show the DEGs in two sets of sample data from each of the three microarrays. Then use Venn diagram software to identify the DEGs common to these three datasets. A total of 310 common DEGs were detected in the ccRCC organization, including 177 genes were down-regulated (logFC <1) and 133 genes were up-regulated (logFC> 1) (Figure 1D).

Figure 1. Differential expression of genes in the two sets of samples, Venn diagram, PPI network and the most significant DEG module. (A) GSE66272 data, (B) GSE100666 data, and (C) GSE105261 data. Red points represent upregulated genes screened on the basis of fold change > 1.0 and a corrected P-value of < 0.05. Green points represent downregulation of gene expression screened on the basis of fold change < 1.0 and a corrected P-value of < 0.05. Black points represent genes with no significant difference. (D) DEGs were selected with a |fold change| >1 and P-value <0.05 among the mRNA expression profiling sets GSE66272, GSE100666 and GSE105261. The 3 datasets exhibited an overlap of 310 genes. (E) The PPI network of DEGs was constructed using Cytoscape. (F) The most significant module was obtained from a PPI network with 21 nodes and 74 edges. Upregulated genes are marked in light red; downregulated genes are marked in light blue. Abbreviations: FC: fold change; GEO: Gene Expression Omnibus; DEGs; differentially expressed genes; PPI: protein–protein interaction.

Gene ontology of DEGs and KEGG pathway analysis in renal cell carcinoma

In order to analyze the biological classification of DEGs, DAVID was used for functional and pathway enrichment analysis. GO analysis showed that BP changes in DEGs were significantly enriched in extracellular matrix organization, excretion, response to drugs, angiogenesis and response to hypoxia. The changes of CC were mainly concentrated in exosome, extracellular region, apical plasma membrane, extracellular space, and basolateral plasma membrane. The changes of MF mainly focus on protein homodimerization activity, ATPase binding, the same protein binding, extracellular matrix binding and extracellular matrix structure composition. Analysis of the KEGG pathway suggested that the DEGs were significantly enriched in carbon metabolism, biosynthesis of antibiotics, aldosterone-regulated sodium reabsorption, collecting duct acid secretion and phagosome (Table 1).

Modular analysis via DEGs protein-protein interaction (PPI) network

To determine the modular structure, STRING online database (available online: http://string-db.org) and the Cytoscape software were used the combined 310 DEGs to construct the PPI network of DEGs (Figure 1E) with the most important module obtained by using Cytoscape (Figure 1F). The functional analyses of genes involved in this module were analyzed using DAVID. The results of GO analysis showed that the changes of BPs of important module genes were significantly rich in collagen catabolic process, extracellular matrix organization, response to peptide hormone, cellular response to amino acid stimulus and platelet degranulation (Table 2). The CC changes of important module genes are mostly concentrated in extracellular region, collagen trimer, endoplasmic reticulum lumen, extracellular matrix, and collagen type IV trimer (Table 2). The MF changes of important module genes are mainly concentrated on extracellular matrix structural constituent, platelet-derived growth factor binding and protein binding (Table 2). KEGG pathway analysis showed that important modular genes were mostly enriched in receptor interaction, protein digestion and absorption, PI3K-Akt signaling pathway, focal adhesion and amoebiasis (Table 2).

Hub gene selection and survival analysis based on TCGA

30 hub genes were generated through filtering according to the criterion of degree of connectivity >10 (each node had more than 10 interactions). The 20 genes with the most significant degree of connectivity were VCNA, CAV1, EPCAM, EGF, GPC3, COL1A1, TIMP1, CCL5, CSF1R, DCN, VEGFA, KNG1, ITGB2, IGFBP3, ALB, CD163, ITGA5, CASR, AQP2, SLC12A1 (Figure 2A). Table 3 shows the names, abbreviations and functions of these genes. The cBioPortal online platform was used to analyze the network of hub gene and their co-expressed genes (Figure 2B). The bioprocess analysis of central genes and the enrichment analysis of the Kyoto Encyclopedia of Genes and Genomes (KEGG) are shown in Figures 2C, 2D. Hierarchical clustering indicated that analysis of hub genes essentially allowed kidney cancer samples to be distinguished from non-cancer samples (Figure 2E).

Figure 2. Interaction network and analysis of hub genes. (A) The 20 most important hub genes were screened using the Cytoscape software plugin cytoHubba. (B) Hub genes and their co-expressed genes were analyzed using the cBioPortal. Nodes with a bold black outline represent hub genes. Nodes with thin black outlines represent co-expressed genes. (C) Biological processes functional annotation analysis of hub genes was performed using ClueGO and CluePedia. Different colors of nodes refer to the functional annotation of ontologies. Corrected P value <0.01 was considered statistically significant. (D) KEGG functional annotation analysis of hub genes was performed by ClueGO and CluePedia. Different colors of nodes refer to the functional annotation of ontologies. Corrected P value <0.01 was considered statistically significant. (E) Hierarchical clustering heatmap of the 20 most important hub genes was constructed from a TCGA cohort. Red indicates that the relative expression of genes was upregulated, green indicates downregulation, and black indicates that no significant change in gene expression was observed; gray indicates that signal strength was not high enough to be detected. Abbreviation: TCGA: the cancer genome atlas program; KEGG: Kyoto Encyclopedia of Genes and Genomes.

In order to determine whether the pivotal genes in Clear cell renal cell carcinoma had clinical relevance, correlation analysis was performed with clinical correlation to kidney cancer outcomes according to the report of the Cancer Genome Atlas (TCGA) on the renal cancer data set. Using data from GEPIA, it was noted that ccRCC patients exhibiting genomic alteration in TIMP1 demonstrated a reduction in overall and disease-free survival (P=6.8x10^-7 for overall survival and P=3.7 x10^-5 for disease-free survival) (Figure 3A, 3B). Therefore, it was considered that the hub gene TIMP1 may play a key role in the progression of clear cell renal cell carcinoma.

Figure 3. Univariate survival analysis of hub genes was performed using Kaplan-Meier curves. TIMP1 showed significant differences in both OS (A) and DFS (B) in ccRCC samples (Logrank P < 0.05). (C) Transcriptional levels of TIMP1 expression were found expressed in 533 ccRCC tissues compared with 72 normal tissues (p<0.0001). Abbreviation: DFS: disease-free survival; OS: overall survival; ccRCC: clear cell renal cell carcinoma.

Analysis of the most significant hub gene

Based on RNA sequence data from the TCGA database, TIMP1 mRNA expression was compared between kidney tumor samples and adjacent normal tissues. Transcriptional level data demonstrated that TIMP1 expression was highly in 533 ccRCC tissues compared with 72 normal tissues (Figure 3C). Consistent with this finding, Oncomine analysis of cancerous and normal tissues showed that TIMP1 was greatly overexpressed in ccRCC samples from different data sets (Figure 4A). As shown in Figure 4B, TIMP1 mRNA expression in ccRCC samples was significantly correlated with late clinical stage, and the highest TIMP1 mRNA expression was observed in stages 3 and 4. Similarly, the relationship between TIMP1 mRNA expression and different pathological grade was determined, which suggested that TIMP1 mRNA expression was significantly correlated with pathological grade (Figure 4C). In general, increased expression of TIMP1 mRNA in ccRCC patients was significantly correlated with advanced clinicopathologic parameters.

Figure 4. Transcriptional expression of TIMP1 significantly correlated with advanced clinicopathological parameters and poor survival outcomes in ccRCC patients. (A) Oncomine analysis of cancer vs. normal tissue of TIMP1. Heat maps of TIMP1 gene expression in clinical hepatocellular carcinoma samples vs. normal tissues. 1. Clear Cell Renal Cell Carcinoma vs. Normal, Higgins, Am J Pathol, 2003 [49]. 2. Clear Cell Renal Cell Carcinoma vs. Normal, Yusenko, BMC Cancer, 2009 [50]. 3. Clear Cell Renal Cell Carcinoma vs. Normal, Jones J, Clin Cancer Res, 2005 [51]. 4. Clear Cell Renal Cell Carcinoma vs. Normal, Gumz ML, Clin Cancer Res, 2007 [52]. (B) Transcriptional expression of TIMP1 was significantly correlated with AJCC stage, patients who were in a more advanced stage tended to express higher mRNA expression of TIMP1. (C) Transcriptional expression of TIMP1 was significantly correlated with ISUP grade. Patients in a more advanced grade tended to exhibit elevated TIMP1mRNA expression. Abbreviations: ccRCC: clear cell renal cell carcinoma.

GSEA was used for hallmark analysis of TIMP1. The results suggest that the most significantly involved pathways included epithelial-mesenchymal-transition, inflammatory-response and IL6/JAK/STAT3 signaling, as displayed in Figures 5A–5C. In addition, the expression profiles of the 100 most significant genes are displayed in a heat map (Figure 5D).

Figure 5. Significantly related genes and hallmarks pathways in ccRCC obtained by GSEA. GSEA was used to perform hallmark analysis for TIMP1 (A–C). The most significant pathways included epithelial-mesenchymal-transition, inflammatory-response, and IL6/JAK/STAT3 signaling. (D) Transcriptional expression profiles of the 100 most significant genes expressed as a heat map.

TIMP1 expression correlates with interstitial phenotype in ccRCC cells

GSEA enrichment analysis indicated that TIMP1 was involved in the EMT process. In order to understand the relationship between TIMP1 and the mesenchymal phenotype in ccRCC cells, TIMP1 and EMT represent markers (such as E-cadherin and N-) in 2 different ccRCC cell lines (A498 and Caki-1) and HK2 cells (normal renal cell lines) were examined by qRT-PCR and Western blot (Figure 6A–6D). As consistently shown by qRT-PCR and Western blot analysis, TIMP1 was abundantly expressed in A498 and Caki-1 mesenchymal cells, which exhibit higher N-cadherin and lower E-cadherin expression than HK2 cells, indicating that that TIMP1 level may be related to the mesenchymal phenotype of ccRCC cell line.

Figure 6. TIMP1 expression was associated with mesenchymal phenotype in ccRCC cell lines. Relative mRNA expression of TIMP1 (A), E-cadherin (B) and N-cadherin (C) in 2 ccRCC cell lines and a normal cell line presented separately as histograms. (D) TIMP1, E-cadherin, and N-cadherin protein levels as determined by Western blot analysis. β-actin was used as the internal control. (* p<0.05, ** p<0.01, *** p<0.001 by t-test).

To demonstrate the correlation between mesenchymal genes and TIMP1 expression, 200 about mesenchymal genes in GSEA Molecular Signatures Database HALLMARK EPITHELIAL MESENCHYMAL TRANSITION gene set (For more details, please refer to the follow website link: https://www.gsea-msigdb.org/gsea/msigdb/cards/HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION.html) were first retrieved. Next, Pearson and Spearman correlation test were carried out to calculate the correlation and P value between mesenchymal genes and TIMP1 expression within 3 datasets (GSE66727, GSE100666, GSE105261). The results were shown in Supplementary Table 1. Totally, expression values of 190/200 genes were found in three datasets and all genes showed significant correlation with TIMP1.

Knockdown of TIMP1 induces MET of A498 and Caki-1 cells

To study the effect of TIMP1 knockdown on EMT in ccRCC cells, lentiviruses with TIMP1 shRNA or NC were transfected into A498 and Caki-1 cells. Stably passaged A498-KD, A498-NC, Caki-1-KD and Caki-1-NC cells were thus obtained. qRT-PCR demonstrated a dramatic decrease in TIMP1 mRNA expression in A498 and Caki-1 cells following lentivirus transfection (Figure 7A, 7B). Western blot analysis confirmed that TIMP1 protein was inhibited in both A498 and Caki-1 cells (Figure 7C). Through TIMP1 knockdown, Western blot analysis and qRT-PCR confirmed the enhanced expression of epithelial markers (e.g. E-cadherin) and interstitial markers (e.g. N- cadherin) expression is reduced in A498-KD and Caki-1-KD cells (Figure 7C–7G). These results suggest that the knockdown of TIMP1 in A498 and Caki-1 cells induced transition to the epithelial phenotype.

Figure 7. Knockdown of TIMP1 induced MET in A498 and Caki-1 cells. qRT-PCR (A, B) analysis of TIMP1 expression in A498 and Caki-1 cells transfected with TIMP1 shRNA (KD) and normal controls (NC). (C) Western blot analysis of TIMP1 expression in A498 and Caki-1 cells transfected with TIMP1 shRNA (KD) from a normal control (NC). Expression levels of protein (C) and mRNA (D–G) of the EMT markers are shown. (* p<0.05, ** p<0.01, *** p<0.001 by t-test).

Knocking down TIMP1 inhibits migration and invasion of A498 and Caki-1 cells

The cell migration and invasion characteristics are considered to be important consequences of EMT [9, 10], so the effect of TIMP1 knockdown on the migration and invasiveness of A498 and Caki-1 cells was studied using transwell migration and invasion assays. The migration rate and invasion rate of A498-KD and Caki-1-KD cells decreased compared with A498-NC and Caki-1-NC cells, respectively (Figure 8A–8D).

Figure 8. TIMP1 knockdown inhibited the migration and invasion of A498 and Caki-1 cells. (A) Migration and Matrigel invasion assays of A498 cells transfected with TIMP1 shRNA or NC were evaluated. (C) Migration and Matrigel invasion assays of Caki-1 cells transfected with TIMP1 shRNA or NC were evaluated. (B, D) Migrated and invaded cells were counted in 3 random 100× fields (* p<0.05, ** p<0.01, *** p<0.001 using a t-test).

Gene expression profiling

We downloaded the publicly available GSE66272, GSE100666, and GSE105261 gene expression datasets from NCBI-GEO. These studies compared gene expression profiles in kidney tumor and normal kidney tissue samples. GSE66272 data were generated with the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array), and were based on 27 tumor and 27 normal tissue samples. GSE100666 data were generated with the GPL16951 platform (Phalanx Human OneArray Ver. 6 Release 1), and were based on 3 ccRCC and 3 normal tissue samples. GSE105261 data were generated with the GPL10558 platform (Illumina HumanHT-12 V4.0 expression BeadChip), and were based on analyses of 9 tumor and 9 paracancerous normal kidney tissue samples.

Standardization and elucidation of DEGs

After downloading matrix files for each GEO dataset, we utilized a Robust Multi-Array Average (RMA) algorithm in order to adjust data for background signals, and to conduct quantile normalization and final oligonucleotides-per-transcript summarization with a median polish algorithm. The R Impute package was used to impute missing values based on a k-nearest neighbor (KNN) approach. Bayes methods were used to adjust for batch effects with the sva R package, as published previously, with R v3.6.1 being used for all analyses [37].

DEGs were identified via comparing ccRCC and normal kidney tissue gene expression profiles based on adjusted P-values, fold change values, and Benjamini and Hochberg FDR values. We removed any probe sets that did not correspond to a gene symbol, while genes with multiple probe sets were averaged. DEGs were considered to be significant if they met the criteria: |log2FC| >1.00, P <0.05. An online Venn diagram program was then used to identify shared DEGs among these three datasets, with data being uploaded in a TXT format. Genes were up- and down-regulated if logFC values were > 0 and < 0, respectively.

Functional enrichment analysis

The DAVID (http://david.ncifcrf.gov) (version 6.7) database was used to conduct functional annotation and pathway enrichment analyses for identified DEGs [38]. The KEGG database was used for pathway enrichment analysis of these DEGs in an effort to gain high-level insights into their potential functional relevance [39]. GO enrichment analyses were conducted to assess the enrichment of these DEGs in specific biological processes, molecular functions, and cell components (BPs, MFs, and CCs, respectively) [40]. Herein, DAVID was utilized for these enrichment analyses.

PPI network analysis

DEGs were grouped into a PPI network with the STRING database (http://string-db.org) [41], with a score of >0.4 being used as a significance cutoff for network construction. The network was then visualized and analyzed with the Cytoscape tool (v3.7.1) [42]. The Cytoscape Molecular Complex Detection (MCODE) plugin (v1.5.1) was further used for topological clustering within the network as a means of identifying significantly interconnected gene modules within the overall PPI network [43]. The following criteria were used for significant module identification: degrees of cut-off = 2, node score cut-off = 0.2, Max depth = 100, K-score = 2.

Hub gene selection and analysis

Hub genes were identified as genes that had a ≥10 degree of connectivity. We then utilized cBioPortal (http://www.cbioportal.org) [44] to identify genes that were co-expressed with these hub genes. We additionally utilized the ClueGO Cytoscape plugin, which allows for the visualization of non-redundant terms associated with gene clusters in networks that have been grouped based on functionality [45]. We utilized ClueGO v2.5.4 and CluePedia v1.5.4 to identify and visualize biological processes associated with these hub genes [46]. We additionally performed hub gene hierarchical clustering.

Survival analysis

Differences in survival between groups were evaluated via a Kaplan-Meier approach, with disease-free survival (DFS) as the primary endpoint. DFS was defined as the time between curative treatment and progression, death, or second-line treatment. As a secondary endpoint in these analyses, we assessed overall survival (OS) from diagnosis/treatment to death or most recent follow-up. After comparing survival between groups using Kaplan-Meier curves, 95% confidence intervals (95% CIs) were generated and curves were compared via log-rank tests. The sum of hub gene weights was used to determine overall scores.

Data processing of Gene set enrichment analysis (GSEA)

A GSEA approach was used to analyze TCGA data with the Category v2.10.1 package. In individual analyses, Student’s t-tests were used to generate scores to assess consistent changes in DEG expression within pathways of interest. A 1000x permutation test was employed to identify pathways that were significantly altered, and the potential for false-positive result detection was controlled by correcting P-values using the Benjamini and Hochberg false discovery rate approach [47]. Genes were found to be significantly related when adj. P was < 0.01 and FDR < 0.25. R v3.6.1 was used for these statistical analyses.

Cell culture

F-12, McCoy’s 5A, and DMEM were used to culture HK-2, Caki-1, and A498 cells, respectively. In all three cases, media was supplemented with 10% FBS, 2 mM L-glutamine, and penicillin/streptomycin. All cell culture was conducted in a humidified 37°C 5% CO2 incubator.

Lentiviral preparation

HEK 293T cells (American Type Culture Collection, USA) were co-transfected with lentiviral vectors and packaging plasmids using Lipofectamine 3000 (Invitrogen, Shanghai, China) based on provided directions. Following a 48 h incubation, lentivirus-containing supernatants were harvested.

Lentiviral transduction

A total of 3×105 ccRCC cells were added to 6 cm tissue culture dishes and were incubated overnight, after which 200 μl of an appropriate lentiviral supernatant and polybrene (8 μg/ml) were added to each well. Cells were then incubated at 37°C, and on day three the virus-containing media was removed. Cells were then rinsed using PBS, and puromycin-containing media (2 μg/ml) was added. After overnight incubation, a subset of cells were collected for western blotting in order to identify stably transduced cells for further experimental utilization. Infection of target cells with lentiviral particles.

qRT-PCR

Trizol (Invitrogen) was then used to extract total cellular RNA, after which SYBR® Premix Ex Taq (TaKaRa) was used for qRT-PCR reactions based on previously published protocols [48]. Primers used were: TIMP1 F: 5’- ACCACCTTATACCAGCGTTATGA -3’; TIMP1 R: 5’- GGTGTAGACGAACCGGATGTC -3’. Β-actin was used to normalize gene expression, with the 2-ΔΔCt approach being used to compare gene expression levels.

Western blotting

RIPA buffer was used to lyse harvested cells. Protein was then separated via SDS-PAGE, transferred to NC membranes (GE Health Care Life Science). Blots were blocked with 5% skim milk in TBST, followed by an overnight incubation with appropriate primary antibodies at 4°C. Antibodies were: anti-TIMP1 (Abcam, ab61224), and anti-Actin (Zen bioscience, 220529). Blots were then probed for 1 h with appropriate secondary antibodies (1:5000), and chemiluminescence was then used for protein visualization.

Transwell assay

A 24-well Transwell system (Greiner Bio-one, Switzerland) was employed for assays examining ccRCC cell invasive activity. The upper surface of these Transwell inserts was first coated using Matrigel (BD Bioscience, USA), after which 200 uL of serum-free media containing 2×105 cells was added to the upper chamber. Next, 500 uL of media containing 10% FBS was added to the lower compartment, and cells were incubated for 48 h at 37°C. A swab was then used to carefully remove cells remaining on the upper surface, while cells that had migrated to the lower surface were treated with 100% methanol for fixation followed by 0.1% crystal violet staining. A total of five fields of view per insert were then visualized at random via microscope (Olympus, Japan; 200x), with the number of cells in each field being quantified.

Ethics approval

The Ethics Committee of Qilu Hospital of Shandong University approved the study.