Genetic alteration of histone lysine methyltransferases and their significance in renal cell carcinoma

Samples with genomic and clinical data

The DNA copy number, gene expression, mutation, clinicopathological data and overall survival datasets of 882 RCC samples used in this research were obtained from The Cancer Genome Atlas (TCGA) at https://genome-cancer.ucsc.edu, including 528 ccRCC, 288 pRCC, and 66 chRCC samples. A total of 50 human HMTs were analyzed. In addition, 78 samples of Tokyo university downloaded from cBioportal were added to persuasively compare the tumor stage and overall survival between 62 SET-domain mutated patients and 16 non-SET-domain mutated patients. The copy number of HMTs is generated by the copy number algorithm genomic identification of significant targets in cancer algorithm analysis: “−2” represents a homozygous deletion, “−1” indicates a heterozygous deletion, “0” represents diploid, “1” indicates a low-level gain, and “2” signifies a high-level amplification. For mRNA expression data, the relative expression and gene expression profiles in the RCC samples were analyzed.

Statistical analysis

Statistical analyses were performed using the R software (R Core Team, 2013), Graphpad Prism (version 7.01; GraphPad, La Jolla, CA, USA), and SPSS (version 18.0; SPSS, Chicago, IL, USA). The correlations between copy number alteration (CNA) and mutation status of 50 HMTs and 882 phenotypes of specimens were analyzed using Chi-square test. The CNA and mRNA expression of 46 HMTs from 882 sequenced RCC specimens were analyzed using Spearman, Kendall, and Pearson correlation tests. KMT2A, KMT2C, KMT2D, and KMT2E were excluded for lack of data. Heatmap of HMTs expression profiles in different subtypes of RCC was conducted by R statistical software. The Student’s t-test was used in calculating differences in mRNA expression levels of each HMT between ccRCC and other subtypes. Kaplan–Meier survival curve was conducted to analysis the impact of CNA or gene expression of different HMTs on survival. Multivariate survival analysis was performed to investigate prognosis factors by Cox regression. The mRNA relative expression levels were analyzed by one-way analysis of variance. P < 0.05 was considered statistically significant.

Gene set enrichment analysis

Gene set enrichment analysis (GSEA) was performed by GSEA software (Version 2.2.2), which was downloaded from the Broad Institute (http://www.broad.mit.edu/gsea). Enrichment map was generated for visualization of the GSEA results. False discovery rate value, normalized enrichment score, and adjusted P-value were calculated to identify the Hallmarks enriched in each phenotype.

Cell culture

293T, HK-2 and ACHN cells were maintained in Dulbecco’s modified Eagle’s medium supplemented with 10% fetal bovine serum in an atmosphere at 37 °C with 5% CO₂. 786-O and OSRC-2 cells were cultured in RPMI-1640 supplemented with 10% fetal bovine serum in an atmosphere at 37 °C with 5% CO₂.

Quantitative real-time PCR

Total RNAs were extracted by MagZol (Invitrogen, Carlsbad, CA, USA) and cDNAs were synthetized by using SYBR Premix Ex TaqTM (TaKaRa, Kusatsu, Japan). Real-time PCR was performed by SYBR Green Realtime PCR Master Mix (TOYOBO, Osaka, Japan). Amplification conditions were as follows: 95 °C for 15 s, 60 °C for 30 s, 72 °C for 30 s for 40 cycles in a 20 μl reaction mix containing 2× SYBR Green. Primers for the reaction are provided in Table S1.

Plasmids

The shRNA plasmids for SETD2 and EZH2 knockdown were constructed from pSicoR (#11579; Addgene, Watertown, MA, USA) with target sequences of shSETD2: TAGTACACCAAGACTCCAG, and shEZH2: CCAACACAAGTCATCCCATTA. All plasmids were verified by sequencing.

Cell viability, cell proliferation, cell migration, and invasion assays

Cell viability was assessed at 0, 24, 48, 72, and 96 h upon treatments by the 3-(4,5-dimethylthiazol-2-yl)-5-(3-carboxymethoxyphenyl)-2-(4-sulfophenyl)-2H-tetrazolium, inner salt (MTS) method (Promega, Madison, WI, USA, #0000253755) according to the manufacturer’s instructions. The MTS have six replications. Cell proliferation was estimated using the cell-light EdU Apollo 568 in vitro kit (Ribobio, GuangZhou, China. #C10310-1) according to the manufacturer’s instructions. Migration and invasion assays were performed using uncoated and Matrigel-coated Transwell inserts according to manufacturer’s instructions. All experiments were performed in triplicate.

Genetic alterations of HMTs in renal cell carcinoma

Copy number alteration and somatic mutations are crucial mechanisms for oncogenesis or inactivating tumor suppressor genes in the occurrence and development of cancer (Albertson et al., 2003). We hypothesized that genetic alterations of HMTs play significant roles in the development and progression of RCC (Pires-Luís et al., 2015). To systematically identify HMTs’ potential of being biomarkers of diagnosis and prognosis of RCC, we first analyzed copy numbers and mutations of 882 RCC samples from the TCGA database via Cancer Browser (Yao et al., 2016). In Cancer Browser CNA datasets of RCC, as previously described, CNA was counted as “−2,” “−1,” “0,” “1,” and “2,” and the average CNA rate of 50 HMTs was −0.0066. We found a different pattern of altered copy number and mutation of HMTs in RCC. Strikingly, as showed in Table 2, seven HMTs (NSD1, PRDM6, MECOM, KMT2C, EZH2, PRDM14, and KMT2E) exhibited high-level amplification in more than 0.5% of RCC samples, and two of these seven HMTs (NSD1 and PRDM6) had a much higher amplification rate in over 10% of ccRCC samples. Three HMT genes, SETD2, SETD5, and SETMAR, showed homozygous deletion in more than 10% ccRCC samples. Intriguingly, NSD1 and PRDM6, the highest frequency of the top two HMTs with high-level amplification, were located in chromosome 5q; whereas SETD2, SETD5, and SETMAR, the highest frequency of the top three HMTs with homozygous deletion, were located in chromosome 3p. Additionally, KMT2C and SETD2 exhibited somatic mutations in more than 6% of ccRCC samples. Several studies revealed that SETD2 was frequently mutated in ccRCC and SETD2-mutated signal pathway played a vital role in the process of oncogenesis (Li et al., 2016a; Wang et al., 2015; Li et al., 2016b). Analogously, in our study SETD2 showed the highest mutation rate (11.51%) among 50 HMTs in ccRCC (Table 2).

Furthermore, HMTs showed different frequencies of CNA and mutation in different subtypes of RCC. Among the seven most frequently amplified HMTs, the frequencies of NSD1 and PRDM6 of high-level amplification were markedly higher in 528 ccRCC with more than 13.45% of tumor samples exhibiting amplificated when compared with that of 288 pRCC and 66 chRCC subtypes, which was less than 1.04% of tumor samples (Fig. 1A). Also, SETD2 and SETD5 exhibited the highest frequency of homozygous deletion in ccRCC. However, neither of them exhibited homozygous deletion in pRCC and chRCC subtypes (Fig. 1B). Additionally, in 213 ccRCC, 113 pRCC, and 65 chRCC, of the most commonly mutated HMTs, SETD2, and KMT2C were most frequently mutated in ccRCC and pRCC subtypes, whereas SETD2 was mutated in less than 1.54% of tumor samples and KMT2C did not exhibit mutation in chRCC (Fig. 1C). The above data indicates that the subtype of ccRCC has a higher CNA and somatic mutation frequency in several HMTs, including amplification of NSD1 and PRDM6, homozygous deletion of SETD2, SETD5, and SETMAR, and mutation of KMT2C and SETD2.

Gene expression and CNA profiling of HMTs in renal cell carcinoma

The correlation between gene expression and copy number has been associated with the detection of human oncogenes. Therefore, we analyzed the correlation between copy number and gene expression level of 46 HMTs from 882 sequenced RCC samples which were divided into ccRCC and non-ccRCC (KIRP + KICH) groups. Four HMTs (KMT2A, KMT2C, KMT2D, and KMT2E) were excluded for lack of data. The rank correlation coefficients in the three statistical tests were similar for HMTs (Table 3). As shown in Table 3, HMTs were ranked based on the Spearman correlation coefficient, and the correlations between CNA and mRNA expression of 46 HMT genes were positive, and three of them (NSD1, WHSC1L1, SETDB1) showed a Spearman correlation coefficient greater than 0.5. (P < 0.0001). WHSC1L1 exhibited the highest correlation coefficient by Pearson (r = 0.669), Spearman (r = 0.689) and Kendall (r = 0.596). Meanwhile, to investigate the degree of discretization of the data sets, the standard deviations of mRNA expression between the ccRCC and non-ccRCC subtypes was performed. As shown in the Table 3, PRDM6 in ccRCC group had a higher degree of discretization than that of in non-ccRCC group.

Expression levels of the 46 HMTs was compared in the heatmap. Compared with non-ccRCC, mRNA levels of 10 HMTs (PRDM1, PRDM8, MECOM, PRDM16, SETD7, PRDM5, ASH1L, NSD1, SUV39H2, and SETDB1) were significantly higher (P < 0.001) and 10 HMTs (PRDM12, SUV420H2, SETMAR, SETD1A, SETD2, SETD4, PRDM4, SETD8, DOT1L, and SETD1B) were significantly lower (P < 0.001) in ccRCC (Fig. 2; Table S2).

SETD2 and KMT2C mutations in clear cell renal cell carcinoma

As described previously, we found that SETD2 and KMT2C are most frequently mutated HMTs in ccRCC, at rates of 11.51% and 6.10% (Table 2). A systematic analysis of these mutation profiles was performed in ccRCC samples. Results showed that a total of 78 SETD2 mutations were found whereas 14 mutations were excluded because their data were untested. A total of 64 mutations were valid, including 18 nonsense mutations, 21 missense mutations, five frameshift insertions, eight frameshift deletions, six splices, and three other mutations. In addition, 24 KMT2C gene mutations were identified, including 10 missense mutations, five nonsense mutations, three frameshift deletions, two splice, and four other mutations (Fig. 3A). A mutation map was performed to display the distribution of SETD2 and KMT2C mutations (Fig. 3B). By systematic analysis of the mutation distribution, we found that SETD2 mutations were more likely to occur at SET domain area. Taking account of the crucial function of SET domain in SETD2, we predicted that mutations at the SET domain might result in the loss of methyltransferases features of SETD2 and poor prognosis of ccRCC patients. Therefore, we performed a Kaplan–Meier plots to investigate the clinical features between SETD2 SET domain mutation group and non-SET domain mutation group; results indicated that SETD2 SET domain mutation group was featured with advanced tumor stage and poor prognosis in ccRCC (Fig. 4A).

Identification of SET-Domain associated biological process by GSEA

To identify SET-Domain associated biological process and function loss of SETD2 on a generalized level, GSEA was performed by using high throughput RNA-sequencing data of the TCGA cohort. The mutation state of SET domain was used as the phenotype label. Among all the predefined Hallmarks gene sets, DNA repair, E2F targets, G2M checkpoint, and mitotic spindle were found to be significantly associated with SET domain mutation in the TCGA cohort (Fig. 4B).

HMTs CNA and expression and clear cell renal cell carcinoma patient survival

To explore the clinical association of genetic alterations of HMTs in RCC, we investigated the association between CNA, mRNA expression, and overall patient survival in 452 ccRCC samples, 76 ccRCC samples were excluded because their detailed survival data is not available in TCGA database. First of all, samples were divided into the following three groups for each HMT: amp/gain (high-level amplification/low-level gain), diploid, and deletion. For six HMTs (EZH2, NSD1, PRDM6, SETD2, SETD5, and SETMAR), copy number amp/gain and deletion were significantly related to poorer survival in RCC patients (P < 0.05). Deletions of KMT2C and PRDM6 were related to shorter survival, however, only amp/gain of EZH2 and PRDM14 was more likely related to poorer survival. More importantly, deletion of NSD1 was significantly related to poorer survival; amp/gain of NSD1 was significantly related to longer survival, compared with diploid or deletion groups (Fig. 5A; Fig. S1).

To analyze the relationship between HMTs expression and overall survival of ccRCC patients, they were divided into low (n = 226) and high (n = 226) expression groups based on the mRNA expression of each HMT. High mRNA levels of EZH2, PRDM6, SETD5, and SETMAR were significantly associated with shorter survival in ccRCC patients, whereas only high NSD1 expression was correlated with longer survival in ccRCC (P < 0.05) (Fig. 5A; Fig. S1).

A multivariate analysis by Cox model (n = 428) was performed to investigate the capability to predict poor prognosis of each HMT compared with standard prognostic markers, such as age at diagnosis, gender and tumor stage (stage I–stage IV). Results indicated that amp/gain of ASH1L had a hazard radio (HR), a ratio of death probabilities, of 1.533 compared with non-amp/gain of ASH1L in RCC patients. In addition, deletion of PRDM8 was significantly associated with shorter survival (HR = 1.516, P < 0.05) in RCC patients. High mRNA levels of SETD1A (HR = 1.275) and PRDM9 (HR = 1.258) was significantly associated with shorter survival in RCC patients (P < 0.05). However, higher expression of PRDM8 (HR = 0.739), EHMT1 (HR = 0.723), ASH1L (HR = 0.711), WHSC1L1 (HR = 0.674), and SUV420H1 (HR = 0.616) was negatively associated with shorter survival in RCC patients (P < 0.05). Also, amp/gain of PRDM6 (HR = 0.703), PRDM9 (HR = 0.662), PRDM7 (HR = 0.640), SETD1A (HR = 0.638), NSD1 (HR = 0.598), and DOT1L (HR = 0.592) was negative correlated with shorter survival in RCC patients (P < 0.05). (Fig. 5B; Table S3).

RT-PCR analysis of mRNA expression of important HMTs

Quantitative RT-PCR was performed to measure the expression level of eight important HMTs (NSD1, PRDM6, MECOM, EZH2, ASH1L, SETD1A, WHSC1L1, and SETD2). HK-2, a renal tubular epithelial cell line, was used as the control group. Relative expression of the eight HMTs in kidney cancer cell lines compared with HK-2 cell line was shown in Fig. 6. It showed that mRNA levels of EZH2 were more than fourfold higher in kidney cancer cell lines. In contrast, mRNA levels of WHSC1L1, ASH1L, and NSD1 were more than onefold lower in kidney cancer cell lines. Notably, for SETD2 and MECOM, mRNA levels were more than threefold and sixfold lower in kidney cancer cell lines, respectively. However, mRNA level of PRDM6 did not show significant difference between RCC and control cell lines. These results indicated that there was a correlated change between CNAs and mRNA expression.

Comprehensive identification of important HMTs in RCC

According to previous data of these CNA, mutation, mRNA expression, and clinical outcome, Table 4 listed the integrative score of important HMTs. Every category counted as “+” when an HMT met the criteria. Especially, SETD2 and KMT2C counted as “++” in CNA/Mutations category because they met the two criteria. As shown in Table 4, six HMTs (ASH1L, PRDM6, NSD1, EZH2, WHSC1L1, SETD2) had a score of more than 2, suggesting that these six HMTs may play important roles in RCC oncogenesis.

SETD2 inhibits cell proliferation in vitro, EZH2 promotes cell proliferation, migration and invasion in vitro

We performed MTS assays, Transwell migration/invasion assays and EdU assays in 786-O cell line to detect the function on cell proliferation and migration of SETD2 and EZH2. Assays showed a higher proliferative ability in 786-O cell line with SETD2 knockdown. However, the migration and invasion ability showed no significant differences with or without SETD2 knockdown. Meanwhile EZH2 knockdown reduced the ability of cell proliferation, migration, and invasion significantly. Results mentioned above (Figs. 7 and 8) were consistent with our preceding findings.