Prognostic value of epithelial-mesenchymal transition markers in clear cell renal cell carcinoma

Differential expression of EMT related genes in ccRCC patients

As shown in Figure 1A, Oncomine datasets demonstrated CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 in 20 types of cancers between tumor and normal tissues. Transcriptional expressions of VIM were significantly elevated in cancer tissue compared with normal tissues, while expression of CDH1, a key epithelial marker, was significantly decreased in cancer tissue in multiple datasets.

Figure 1. Analysis of the six EMT related genes in Oncomine database and TCGA database. (A) The Oncomine database was queried for the expression of CDH1, CDH2, SNAI1, SNAI2, VIM, and TWIST1 in the available datasets based on the following criteria: 1) “Cancer Type”; 2) “Gene: CDH1, CDH2, SNAI1, SNAI2, VIM, or TWIST1”; 3) “Data Type: mRNA”; 4) “Analysis Type: Cancer vs Normal Analysis”, and 5) Threshold Setting Condition (p<0.001, fold change >2, gene rank = top 10%). The 'red cells' represents gene overexpression and the 'blue cells' represent gene underexpression. The color intensity equals the percentile, i.e. Top 1%, 5%, or 10% significantly over- or underexpressed (see the legend below the grid). We found that CDH1 and SNAI2 was underexpressed in the kidney cancer vs normal datasets, while VIM was highly overexpressed. (B–G) Differential mRNA expression of six EMT related genes in clear cell renal cell carcinoma (ccRCC) tumor samples and adjacent normal tissues from TCGA. Epithelial marker CDH1 mRNA expression was significantly lower in tumor samples compared with adjacent normal tissues (B); Most mesenchymal markers (CDH2, SNAI1, VIM, and TWIST1) mRNA expression was elevated in tumor samples compared with adjacent normal tissues (C–G).

We further compared the mRNA expression of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 between ccRCC tumor samples and adjacent normal tissues respectively based on RNA-sequence data from TCGA database. Consistent with Oncomine data, CDH1 expression was lower in ccRCC primary tumors in comparison with adjacent normal tissues (Figure 1B). However, expression of mesenchymal markers, CDH2, SNAI1, VIM, TWIST1, were significantly higher in ccRCC primary tumors (Figure 1C–1G).

Neighbor genes and hierarchical partitioning of hub genes from TCGA

As was shown in Supplementary Figure 1, network of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 and their 49 frequently genetic altered neighbor genes was integrated and constructed using cBioPortal.

Demographic and clinical characteristics of ccRCC patients in TCGA and FUSCC cohorts

The TCGA cohort comprised 337 (65.31%) male patients and 179 (34.69%) female patients. The median age of the 516 ccRCC patients was 60.5 years, with a range from 26 to 90 years. Information of TNM stage, AJCC stage, ISUP grade, laterality was shown in Table 1. The median follow-up time was 40.6 months and 172 (33.33%) patients died during follow-up. Besides, 114 (20.09%) patients developed progression or recurrence after surgery.

The FUSCC cohort consisted of 248 (67.57%) male patients and 119 (32.43%) female patients. The median age of the 367 ccRCC patients was 56 years, with a range from 21 to 86 years. The detailed clinical data are shown in Table 1. During follow-up (median: 60 months), 135 (36.78%) patients died and 196 (53.41%) patients developed progression or recurrence after surgery.

Prognostic value of expression of EMT related genes in ccRCC patients

We first examined the prognostic value of mRNA expression of these six EMT key markers in the TCGA cohort. A Kaplan-Meier plot indicated that patients with lower expression of epithelial marker, CDH1, had worse PFS (Figure 2A) as well as OS (Figure 2G). Besides, patients with higher expression of mesenchymal markers, VIM, or TWIST1, had both worse PFS (Figure 2E, 2F) and worse OS (Figure 2K–2J). Higher SNAI1 expression was also associated with moderately worse PFS (Figure 2C) and OS (Figure 2I).

Figure 2. Kaplan Meier survival plot of ccRCC patients in TCGA database according to high and low mRNA expression of six EMT related genes.CDH1 mRNA expression was associated with both worse progression-free survival (p=0.015) and worse overall survival (p=0.003) of ccRCC patients (A, G); CDH2 mRNA expression was not an indicator of either progression-free survival (p=0.593) or overall survival (p=0.075) of ccRCC patients (B, H); Higher SNAI1 mRNA expression was moderately associated with both worse progression-free survival (p=0.054) and worse overall survival (p=0.010) of ccRCC patients (C, I); SNAI2 mRNA expression was not an indicator of either progression-free survival (p=0.105) or overall survival (p=0.242) of ccRCC patients (D, J); Higher VIM mRNA expression was associated with both worse progression-free survival (p<0.001) and worse overall survival (p=0.005) of ccRCC patients (E, K); Higher TWIST1 mRNA expression was associated with both worse progression-free survival (p<0.001) and worse overall survival (p<0.001) of ccRCC patients (F, L).

To confirm the prognostic value of these EMT related genes, we validated it in the FUSCC ccRCC cohort. As shown in Figure 3, patients with lower expression of CDH1 would have worse PFS as well as OS (Figure 3A and 3G). However, elevated expression of mesenchymal markers, SNAI1, SNAI2, VIM, and TWIST1, were associated with both worse PFS and worse OS (Figure 3C–3F, 3I–3L).

Figure 3. Kaplan Meier survival plot of ccRCC patients in FUSCC cohort according to high and low mRNA expression of six EMT related genes. Lower CDH1 mRNA expression was associated with both worse progression-free survival (p=0.016) and worse overall survival (p<0.001) of ccRCC patients (A, G); CDH2 mRNA expression was not an indicator of either progression-free survival (p=0.288) or overall survival (p=0.202) of ccRCC patients (B, H); Higher SNAI1 mRNA expression was associated with both worse progression-free survival (p<0.001) and worse overall survival (p<0.001) of ccRCC patients (C, I); Higher SNAI2 mRNA expression was associated with both worse progression-free survival (p=0.005) and worse overall survival (p<0.001) of ccRCC patients (D, J); Higher VIM mRNA expression was associated with both worse progression-free survival (p<0.001) and worse overall survival (p<0.001) of ccRCC patients (E, K); Higher TWIST1 mRNA expression was associated with both worse progression-free survival (p<0.001) and worse overall survival (p<0.001) of ccRCC patients (F, L).

Integrated prognostic and diagnostic model

After integrating all the significant clinicopathological parameters and gene expression profiles in the Cox regression models (Table 2), we generated the formula: = -0.708×CDH1 expression (ref. Low) + 1.360×SNAI1 expression (ref. Low) + 1.905×VIM expression (ref. Low) + 2.179×TWIST1 expression (ref. Low) + 1.274×T stage (ref. T1-T2) + 1.919×M stage (ref. M0) + 2.021×AJCC stage (ref. I-II) + 2.013×ISUP grade (ref. 1-2) for PFS, and another formula: = -0.564×CDH1 expression (ref. Low) + 1.532×SNAI1 expression (ref. Low) + 1.804×VIM expression (ref. Low) + 1.714×TWIST1 expression (ref. Low) + 1.226×T stage (ref. T1-T2) + 1.778×M stage (ref. M0) + 2.515×AJCC stage (ref. I-II) + 1.954×ISUP grade (ref. 1-2) for OS. The Kaplan–Meier method was used to determine the significant survival outcomes (PFS: p<0.0001; OS: p<0.0001), shown in Figure 4A, 4B. Meanwhile, ROC curves were generated to validate the ability of the logistic model to predict prognosis. The AUC index for the integrated model were 0.886 (p<0.001) for PFS (Figure 4C) and 0.814 for OS (p<0.001) (Figure 4D). To further verify the practical value of the model, ROC curves were constructed to perform external validation using clinicopathological parameters and mRNA expression profiles from TCGA cohort. The AUC index for the integrated model were 0.720 (p<0.001) for PFS (Figure 4E) and 0.684 for OS (p<0.001) (Figure 4F).

Figure 4. Construction and internal validation of integrated prognostic and diagnostic model. All significant clinicopathologic parameters and gene expression profiles was integrated in the Cox regression models, which indicated this formula: = -0.708×CDH1 expression (ref. Low) + 1.360×SNAI1 expression (ref. Low) + 1.905×VIM expression (ref. Low) + 2.179×TWIST1 expression (ref. Low) + 1.274×T stage (ref. T1-T2) + 1.919×M stage (ref. M0) + 2.021×AJCC stage (ref. I-II) + 2.013×ISUP grade (ref. 1-2) for PFS (A), and another formula: = -0.564×CDH1 expression (ref. Low) + 1.532×SNAI1 expression (ref. Low) + 1.804×VIM expression (ref. Low) + 1.714×TWIST1 expression (ref. Low) + 1.226×T stage (ref. T1-T2) + 1.778×M stage (ref. M0) + 2.515×AJCC stage (ref. I-II) + 1.954×ISUP grade (ref. 1-2) for OS (B). The Kaplan–Meier method was used to determine the significant survival outcomes (PFS: p<0.0001; OS: p<0.0001). ROC curves were generated to validate the ability of the logistic model to predict prognosis. The AUC index for the integrated model were 0.886 for PFS (p<0.001) (C) and 0.814 for OS (p<0.001) (D). (E, F) External validation of integrated model using TCGA cohorts. ROC curves were constructed to perform external validation using clinicopathological parameters and mRNA expression profiles from TCGA cohort. The AUC index for the integrated model were 0.720 for PFS (p<0.001) (E) and 0.684 for OS (p<0.001) (F).

EMT markers and ccRCC microenvironment

As seen in Figure 5, EMT makers, including CDH1 (Figure 5A), CDH2 (Figure 5B), SNAI1 (Figure 5C), SNAI2 (Figure 5D) and TWIST1 (Figure 5F), showed significant correlation with stromal process in ccRCC environments (p<0.001). In addition, CDH1 showed negative association with stromal score (r²=-0.191), while stromal score positively correlated CDH2 (r²=0.337), SNAI1 (r²=0.199), SNAI2 (r²=0.201) and TWIST1 (r²=0.305) mRNA expression in ccRCC patients from TCGA cohort.

Figure 5. CDH1, CDH2, SNAI1, SNAI2 and TWIST1 significantly involved in stromal process of ccRCC tumor environment. EMT makers, including CDH1 (A), CDH2 (B), SNAI1 (C), SNAI2 (D), TWIST1 (F), but not VIM (E), showed significant correlation with stromal process in ccRCC environments (p<0.001). In addition, CDH1 showed negative association with stromal score (r²=-0.191), while stromal score positively correlated CDH2 (r²=0.337), SNAI1 (r²=0.199), SNAI2 (r²=0.201) and TWIST1 (r²=0.305) mRNA expression in ccRCC patients from TCGA cohort.

Module analysis and functional annotation

In this study, PPI network, activation and indirect relation was predicted in Figure 6A and PPI network derived from active interaction sources was detailed illustrated with required interaction score equal 0.400 in Figure 6B. Function annotations of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 was enriched in hemophilic cell adhesion and cell-cell adhesion of GO: BP, adherens junction, anchoring junction and cell junction of GO: CC, RPTP-like protein binding, phosphatase binding, protein phosphatase binding and enzyme binding of GO: MF. Participating upstream or downstream signaling pathways enrichment include adherens junction, cell adhesion molecules (CAMs) of KEGG pathways. The details were illustrated in Figure 6C. Hierarchical partitioning using transcriptional expression profiles of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 from FUSCC cohort was performed in Figure 6D. Hierarchical partitioning using transcriptional expression profiles of six hub genes from TCGA cohort was performed in a heat map in Figure 6E.

Figure 6. Module analysis and functional annotations of the six EMT related gene in silico. Protein-protein interaction (PPI), activation and indirect relation were predicted and displayed in association with sig EMT related genes (A). PPI network derived from active interaction sources was detailed illustrated with required interaction score equal 0.400 (B). GO and KEGG functional annotations analysis of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 was enriched in hemophilic cell adhesion and cell-cell adhesion of biologic process, adherens junction, anchoring junction and cell junction of cellular component, RPTP-like protein binding, phosphatase binding, protein phosphatase binding and enzyme binding of molecular function. Participating upstream or downstream signaling pathways enrichment include adherens junction, cell adhesion molecules (CAMs) of KEGG pathways (C). Hierarchical partitioning using transcriptional expression profiles of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 from FUSCC cohort (D). Hierarchical partitioning using transcriptional expression profiles of six hub genes from TCGA cohort was performed in the heat map (E).

Significant genes and pathways obtained by GSEA

A total of 100 significant genes were obtained by GSEA with positive and negative correlation. Importantly, GSEA was used to perform hallmark analysis for CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1. Results suggested that the most involved significant pathways included apical junction, epithelial mesenchymal transition, estrogen response, hypoxia, kras signaling pathway up, inflammatory response, myogenesis, TNF-alpha signaling via NF-κB, etc. The details are shown in Supplementary Figure 2A–2F.

Oncomine database

In this study, transcriptional expression profiles of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 in 20 different common neoplasms were obtained from Oncomine database using Oncomine online database (http://www.oncomine.com) [39]. Difference of transcriptional expression was compared by Student’s t-test. Cut-off of p value and fold change were as following: p value=0.01, Fold Change=1.5, gene rank=10%, Data type: mRNA.

Patients from TCGA cohort and transcriptional expression profile

A total of 533 ccRCC patients with available RNA-sequence data from TCGA database were consecutively recruited in analyses [40]. The gene expression profile was measured experimentally using the Illumina HiSeq 2000 RNA Sequencing platform by the University of North Carolina TCGA genome characterization center. Level 3 data was downloaded from the Cancer Genomics Browser of the University of California Santa Cruz (https://genome-cancer.ucsc.edu/proj/site/hg Heatmap/). Differential transcriptional expression profiles of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 was measured between ccRCC tumor samples and adjacent normal tissues, respectively.

Statistical analysis of TCGA cohort

In this study, significant co-regulated network of genomic profiles of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 and their frequently altered neighbor genes was constructed using cBioPortal, an open-access online resource accessing the TCGA genomics data [41]. Hierarchical partitioning was performed using transcriptional expression profiles of six hub genes in a heat map. Color gradients suggest high (red) or low (blue) expression level.

Among the 533 ccRCC patients, 516 patients had complete clinicopathologic information. Phenotype and expression profiles of hub genes in the 516 ccRCC patients from TCGA were analyzed and displayed. Survival comparison between distinct mRNA expressions groups of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 were analyzed in ccRCC patients. The primary end point for patients was progression-free survival (PFS), and overall survival (OS) was the secondary end point, which was evaluated from the date of first therapy to the date of death or last follow-up. X-tile software was utilized to take the cut-off value of mRNA expression of six hub genes, in concordance of which overall participants were divided to two groups, respectively [42]. The follow-up duration was estimated using the Kaplan-Meier method with 95% confidence intervals (95% CI) and log-rank test in separate curves. Univariate and multivariate analysis were performed with Cox logistic regression models to find independent variables, including age at diagnosis, age at surgery, T stage (ref. T1-T2), N stage (ref. N0), M stage (ref. M0), AJCC (American Joint Committee on Cancer) stage (ref. I-II), ISUP (The International Society of Urological Pathology) grade (ref. 1-2), CDH1 expression (ref. Low), CDH2 expression (ref. Low), SNAI1 expression (ref. Low), SNAI2 expression (ref. Low), VIM expression (ref. Low) and TWIST1 expression (ref. Low). All hypothetical tests were two-sided and p-values less than 0.05 were considered significant in all tests.

Patients and variables from FUSCC cohort

To further validate prognostic implications of six hub genes in real world, a total of 367 ccRCC patients, who have undergone radical nephrectomy in the Department of Urology of Fudan University Shanghai Cancer Center (FUSCC) (Shanghai, China) from Aug. 2008 to Sept. 2016, were consecutively recruited in analyses, with electronic medical records or pathology reports available. Clinical and pathological parameters, specifically age at surgery, clinical manifestation, tumor laterality, TNM stage, ISUP grade classification were summarized. Tissue samples, including ccRCC and normal tissues, were collected during surgery and available from FUSCC tissue bank. A central review of pathology was performed by an experienced pathologist. Clinicopathological characteristics were obtained from electronic records. Patients were regularly followed up by telephone, mail, or in the clinic once every 3 months. All the study designs and test procedures were performed in accordance with the Helsinki Declaration II. The Ethics approval and consent to participate of the current study was approved and consented by the ethics committee of FUSCC.

Real-time quantitative PCR analysis

Transcription level of six hub genes was measured using RT-PCR analysis in 367 ccRCC patients from FUSCC cohort. Total RNA sequence was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA) from 367 paired tumor samples. Total RNA reverse-transcribed reaction was performed using the SuperScript First-Strand cDNA Synthesis System (Invitrogen, Carlsbad, CA). ABI Prism 7900 Sequence Detector (Applied Biosystems) was utilized to realize Real-time PCR reactions. Forward and reverse PCR primers of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 were listed in Supplementary Table 1. According to SYBR Green PCR master mix (Applied Biosystems) manufacturer protocols, a total of 10μL reaction mixture was prepared for each test.

Statistical analysis of FUSCC cohort

Kaplan-Meier method with 95% confidence intervals (95%CI) and log-rank test was applied to analyze different survival outcomes in separate curves in 367 ccRCC patients from FUSCC cohort. Univariate analyses were performed with Cox logistic regression models to find independent variables, including pT stage, pN stage, M stage, AJCC stage, and ISUP grade.

Integrated score was identified as sum of the weight of each significant hub gene. X-tile software was utilized to take the cut-off value. All hypothetical tests were two-sided and P-values less than 0.05 were considered significant in all tests. The receiver operating characteristic curve (ROC) was constructed by predicting the probability of a diagnosis being of high or low integrated score of significant hub gene expression. Area under curve (AUC) analysis was performed to determine the diagnostic ability. In addition, validation of diagnosis and prognostic value of integrated score was constructed using mRNA data from TCGA.

EMT markers and ccRCC microenvironment

ESTIMATE algorithm, available from "estimate" R package, was obtained and utilized to measure to measure stromal components in ccRCC tumor microenvironment. Transcriptional expression levels of EMT markers and stromal scores were matched to sample identification names in the TCGA database. Pearson correlation coefficient analysis was used to reflect the degree of linear correlation between two random variables. The value of r is between -1 and 1. When the value is 1 or -1, the two random variables are completely negatively or positively correlated.

Module analysis and functional annotations

In the present study, protein-protein interaction (PPI), activation, indirect relation and activation expression modular was predicted using Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) (version 10.0) online database was used to predict PPI network of DEGs and analyzing the functional interactions between proteins [43]. An interaction with a combined score >0.4 was considered statistically significant. Significant co-regulated network of six hub genes was constructed using cBioPortal. Subsequently, the gene ontology (GO): BP (biological process), GO: CC (cellular component), GO: MF (molecular function) and KEGG pathways analyses for a total of 48 genes in this module were performed using DAVID [44], and then visualized in bubble chart. Hierarchical partitioning using transcriptional expression profiles of CDH1, CDH2, SNAI1, SNAI2, VIM, TWIST1 from FUSCC cohort was performed in a heat map. Color gradients suggest high (red) or low (blue) expression level.

Data processing of gene set enrichment analysis (GSEA)

TCGA database were implemented with GSEA method using the Category version 2.10.1 package. For each separate analysis, Student’s-t-test statistical score was performed in consistent pathways and the mean of the differential expression genes was calculated. A permutation test with 1000 times was used to identify the significantly changed pathways. The adjusted P values (adj. P) using Benjamini and Hochberg (BH) false discovery rate (FDR) method by default were applied to correct the occurrence of false positive results [45]. The significant related genes were defined with an adj. P less than 0.01 and FDR less than 0.25. Statistical analysis and graphical plotting were conducted using R software (Version 3.3.2).

Summary

We preliminarily tested six key EMT markers and validated their prognostic value in both TCGA cohort and USCC cohort. Among these six EMT markers, CDH1, SNAI1, VIM, and TWIST1 were found to be independent PFS and OS of ccRCC patients.