Comprehensive molecular characterization of clear cell renal cell carcinoma.

Genetic changes underlying clear cell renal cell carcinoma (ccRCC) include alterations in genes controlling cellular oxygen sensing (for example, VHL) and the maintenance of chromatin states (for example, PBRM1). We surveyed more than 400 tumours using different genomic platforms and identified 19 significantly mutated genes. The PI(3)K/AKT pathway was recurrently mutated, suggesting this pathway as a potential therapeutic target. Widespread DNA hypomethylation was associated with mutation of the H3K36 methyltransferase SETD2, and integrative analysis suggested that mutations involving the SWI/SNF chromatin remodelling complex (PBRM1, ARID1A, SMARCA4) could have far-reaching effects on other pathways. Aggressive cancers demonstrated evidence of a metabolic shift, involving downregulation of genes involved in the TCA cycle, decreased AMPK and PTEN protein levels, upregulation of the pentose phosphate pathway and the glutamine transporter genes, increased acetyl-CoA carboxylase protein, and altered promoter methylation of miR-21 (also known as MIR21) and GRB10. Remodelling cellular metabolism thus constitutes a recurrent pattern in ccRCC that correlates with tumour stage and severity and offers new views on the opportunities for disease treatment.

Kidney cancers, or renal cell carcinomas (RCC), are a common group of chemotherapy resistant diseases that can be distinguished by histopathological features and underlying gene mutations.[1] Inherited predisposition to RCC has been shown to arise from genes involved in regulating cellular metabolism, making RCC a model for the role of an oncologic-metabolic shift, commonly referred to as the “Warburg effect”, leading to malignancy.[2] The most common type of RCC, clear cell renal cell carcinoma (ccRCC), is closely associated with VHL gene mutations that lead to stabilization of hypoxia inducible factors (HIF-1α and HIF-2α) in both sporadic and familial forms. PBRM1, a subunit of the PBAF SWI/SNF chromatin remodeling complex, as well as histone deubiquitinase BAP1 and histone methyltransferase SETD2, were recently found to be altered in 41%, 15% and 12% of ccRCCs, respectively[3-5], implicating major roles for epigenetic regulation of additional functional pathways participating in the development and progression of the disease. Oncogenic metabolism and epigenetic reprogramming have thus emerged as central features of ccRCC.

In the present study, clinical and pathological features, genomic alterations, DNA methylation profiles, and RNA and proteomic signatures were evaluated in ccRCC. We accrued more than 500 primary nephrectomy specimens from patients with histologically confirmed ccRCC that conformed to the requirements for genomic study defined by the Cancer Genome Atlas (TCGA)[6], together with matching ‘normal’ genomic material. Samples were restricted to those that contained at least 60% tumor nuclei, (median 85%) by pathological review (clinical data summary provided in Table S1). A data freeze representing 446 samples was generated from at least one analytical platform (‘Extended’ data set) and data from all platforms were available for 372 samples for coordinated, integrative analyses (‘Core’ data set)(Data File S1, Table S2). No substantial batch effects in the data that might confound analyses were detected (Figures S1–S20).

Somatic Alterations

The global pattern of somatic alterations, determined from analysis of 417 samples, is shown in Figure 1A. DNA hybridizations showed that recurrent arm-level and focal somatic copy number alterations (SCNAs) occurred at a fewer sites than is generally observed in other cancers (p<0.0004; Figures S21–S22, Table S3). However, SCNAs that were observed more commonly involved entire chromosomes or chromosome arms, rather than focal events (17% vs 0.4%, Figure 1b). Notably, the most frequent arm-level events involved loss of chromosome 3p3 (91% of samples), encompassing all of the four most commonly mutated genes (VHL, PBRM1, BAP1 and SETD2).

Figure One

Somatic alterations in ccRCC

(a) Top histogram, mutation events per sample; left histogram, samples affected per alteration. Upper heat map, distribution of fusion transcripts and VHL methylation across samples (n=385 samples, with overlapping exome/CNA/RNA-seq/Methylation data); middle heatmap, mutation events; bottom heatmap, copy number gains (red) and losses (blue). Lower chart, mutation spectrum by indicated categories. (b) Left panel, frequency of arm level copy number alterations versus focal copy number alterations. Right panel, comparison of the average numbers of arm level and focal copy number changes in ccRCC, colon cancer (CRC), glioblastoma (GBM), breast cancer (BRCA) and ovarian cancer (OVCA). (c) Circos plot of fusion transcripts identified in 416 samples of ccRCC, with recurrent fusions highlighted.

The data also suggested lower and more variable tumor cellularity[7] in the accrued samples, compared to conventional pathological review (median 54% ± 14%). This may reflect stromal or endothelial cell contributions, or tumor cell heterogeneity. A recent study of multiple samples from single tumors has demonstrated significant regional genomic heterogeneity, but with shared mutations in frequently mutated genes and convergent evolution of other common gene level events.[8] The mutation frequencies of key genes (VHL, PBRM1, etc.), as well as copy number gains and losses found here, were, however, consistent with previous reports. Tumor purity was therefore not determined to be a limitation in the current study.

Arm level losses on chromosome 14q, associated with loss of HIF1A, which has been predicted to drive more aggressive disease[9], were also frequent (45% of samples). Gains of 5q were observed (67% of samples) and additional focal amplifications refined the region of interest to 60 genes in 5q35, which was particularly informative as little has been known about the importance of this region in ccRCC since the 5q gain was initially described. Focal amplification also implicated the protein kinase C member, PRKCI[5] and the MDS1 and EVI1 complex locus MECOM at 3p26, the p53 regulator MDM4 at 1q32, MYC at 8q24 and JAK2 on 9p24. Focally deleted regions included the tumor suppressor genes CDKN2A at 9p21 and PTEN at 10q23, putative tumor suppressor genes NEGR1 at 1p31, QKI at 6q26, and CADM2 at 3p12 and the genes that are frequently deleted in cancer, PTPRD at 9p23 and NRXN3 at 14q24.[10]

Whole exome sequencing (WES) of tumors from 417 patients identified 36,353 putative somatic mutations, including 16,821 missense mutations, 6,383 silent mutations and 2,999 indels, with an average of 1.1 ± 0.5 non-silent mutations per megabase (Figures S23–S25). Mutations from 50 genes with high apparent somatic mutation frequencies (Table S4) were independently validated using alternative sequencing instrumentation (Figure S26). In tumors from 22 patients, whole genome sequencing was also used to validate and calibrate the WES data and confirmed 83% of the WES mutation-calls (Tables S5–S6). In line with results of previous studies (Tables S7–S8), the validated mutation data identified nineteen significantly mutated genes (SMGs) (q< 0.1), with VHL, PBRM1, SETD2, KDM5C, PTEN, BAP1, MTOR and TP53 representing the eight most extreme members (q<0.00001) (Figure 1a). Eleven additional SMGs were of considerably lower significance (q < 0.1 – 0.5) but included known cancer genes. Among all SMGs, only mutation of BAP1 correlated with poor survival outcome (Figure S27). Approximately 20% of cases had none of the 19 recorded SMGs, although many contained rare mutations in other known oncogenes or tumor suppressors, involving survival associations, illustrating the genetic complexity of ccRCC[8] (Figures S28–S30 and Table S9).

Eighty-four putative RNA fusions were identified in 416 ccRCC samples.[11] Eleven of thirteen predicted events (Figure 1c) were validated using targeted methods, consistent with an 85% true-positive rate (Table S10 and Figures S31–S35). A recurrent SFPQ-TFE3 fusion (previously linked to non-clear cell translocation-associated RCC[12]) was found in five samples, all of which were VHL wildtype, indicating either that these tumors are a clear cell variant or that translocation-associated renal tumors may be histologically indistinguishable from conventional ccRCC. Furthermore, the TFE3 protein as well as an X(p11) rearrangement was found in three of those samples, where there were available slides.

DNA Methylation Profiles

We observed epigenetic silencing of VHL in about 7% of ccRCC tumors, which was mutually exclusive with mutation of VHL (Figure 1a), reflecting the central role of this locus in ccRCC.[13] An additional 289 genes displayed evidence of epigenetic silencing in at least 5% of tumors. The top-ranked gene by inverse correlation between gene expression and DNA methylation was UQCRH, hypermethylated in 36% of the tumors. UQCRH has been previously suggested to be a tumor suppressor[14], but not linked to ccRCC. Interestingly, increasing promoter hypermethylation frequency correlated with higher stage and grade (Figure 2a, b).

Figure Two

DNA methylation and ccRCC

(a–b) Overall promoter DNA hypermethylation frequency in the tumor increases with rising stage (a) and grade (b). The promoter DNA hypermethylation frequency is calculated as the percentage of CpG loci hypermethylated among 15,101 loci which are unmethylated in the normal kidney tissue and normal white blood cells (boxplots, median with 95% confidence interval). (c) Volcano plots showing a comparison of DNA methylation for SETD2 mutant versus non-mutant tumors (n=224, Human Methylation 450 platform). Unshaded area: CpG loci with Benjamini-Hochberg FDR=0.001 and difference in mean beta value >0.1 (n=2,557). (d) Heatmap showing CpG loci with SETD2 mutation-associated DNA methylation (from part c); blue to red indicates low to high DNA methylation. The loci are split into those hypomethylated (top panel; n=1,251) or hypermethylated (bottom panel; n=1,306) in SETD2 mutants. Top color bars indicate SETD2 mRNA expression (red: high, green: low) and SETD2 mutation status. Gray-scale row-side color bar on left-hand side represents the relative number of overlapping reads, based on H3K36me3 ChIP-seq experiment in normal adult kidney (http://nihroadmap.nih.gov/epigenomics/); black, high read count. DNA methylation patterns include 14 normal kidney samples. Among the tumors without SETD2 mutations, six (arrowhead) have both the signature pattern of SETD2 mutation and low SETD2 mRNA expression.

We also evaluated the global consequences of mutation in specific epigenetic modifiers. Mutations in SETD2, a non-redundant H3K36 methyltransferase, were associated with increased loss of DNA methylation at non-promoter regions (Figures 2c, 2d). This discovery is consistent with the emerging view that H3K36 trimethylation may be involved in the maintenance of a heterochromatic state[15], whereby DNA methyltransferase 3A (DNMT3A) binds H3K36me3 and methylates nearby DNA.[16] Thus, reductions of H3K36me3 through SETD2 inactivation could lead indirectly to regional loss of DNA methylation.

RNA Expression

Unsupervised clustering methods identified four stable subsets in both mRNA (m1–m4) and miRNA (mi1–mi4) expression datasets (Figures 3a, S36–S39). Supervised clustering revealed the similarity of these new mRNA classes to the previously reported ccA and ccB expression subtypes[17], with cluster m1 corresponding to ccA and ccB divided between m2 and m3 (Table S11). Cluster m4 probably accounts for the roughly 15% of tumors previously unclassified in the ccA/ccB classification scheme. Similarly, the survival advantage previously observed for ccA cases was again identified for m1 tumors (Figure 3b).

Figure Three

mRNA and miRNA patterns reflect molecular subtypes of ccRCC

(a) By unsupervised analyses, tumors separated into four sample groups (i.e. “clusters”), based on either differentially expressed mRNA patterns (left panel, showing 500 representative genes: m1–4) or differentially expressed miRNA patterns (right panel, showing 26 representative miRNAs: mi1–4). (b) Significant differences in patient survival were identified among either the mRNA-based clusters (left panel) or the miRNA-based clusters (right panel). (c) Numbers of samples overlapping between the two sets of clusters, with significant concordance observed between m1 and mi3 and between m3 and mi2; Red, significant overlap (P<1E-5, chi-squared test). (d) mRNA-miRNA correlations, for predicted targeting interactions. Rows indicate miRNAs from part a (indicated by cluster specific color bar); columns, mRNAs (5000 differentially regulated genes selected for average RPKM>10 and at least one predicted miRNA interaction); mRNA-miRNA entries with no predicted targeting show as white. To the right of the correlation matrix, t-statistics (Spearman’s rank) indicate group target enrichment.

The m1-subtype was characterized by gene sets associated with chromatin remodeling processes and a higher frequency of PBRM1 mutations (39% in m1 vs 27% in others, p= 0.027). Deletion of CDKN2A (53% vs 26%; p< 0.0001) and mutations in PTEN (11% vs 1%; p< 0.0001) were more frequent in m3 tumors (Figure S5). The m4 group displayed higher frequencies of BAP1 mutations (17% vs 7%; p=0.002) and base excision repair; however, this group also harbored more mTOR mutations (12% vs 4%; p=0.01) and ribosomal gene sets.

Survival differences evident in miRNA-based subtypes (Figures S40–S44) correlated with the mRNA data (Figures 3b–d). For example, miR-21, previously shown to demonstrate strong regulatory interactions in ccRCC[18] and with established roles in metabolism[15,19,20] correlated strongly with worse outcome, and DNA promoter methylation levels inversely correlated with expression of miR-21, miR-10b, and miR-30a (Table S12–S14). miRNA interactions thus represent a significant component of the epigenetic regulation observed in ccRCC.

Integrative data analyses

We used a combination of approaches for integrative pathway analysis. The HotNet[21] algorithm employs a heat diffusion model, to find subnetworks distinguished by both the frequency of mutation in genes (nodes in the network) and the topology of interactions between genes (edges in the network). In ccRCC, Hotnet identified twenty-five subnetworks of genes within a genome-scale protein-protein interaction network (Table S15 and Figure S45). The largest and most frequently mutated network contained VHL and interacting partners. The second most frequently mutated subnetwork included PBRM1, ARID1A and SMARCA4, key genes in the PBAF SWI/SNF chromatin remodeling complex.

We also inferred activities for known pathways, by using the PARADIGM algorithm to incorporate mutation, copy, and mRNA expression data, with pathway information catalogued in public databases. This method identified a highly significant subnetwork of 2,398 known regulatory interactions, connecting 1,218 molecular features (645 distinct proteins) (Figures S46–S49, Tables S16–S17). Several “active” transcriptional “hubs” were identified, by searching for transcription factors with targets that were inferred to be active in the PARADIGM network. The active hubs found included HIF/ARNT, the transcription factor program activated by VHL mutation, as well as MYC/Max, SP1, FOXM1, JUN, and FOS. Together, these hubs and several less well-studied transcription factors, interlink much of the transcriptional program promoting glycolytic shift, de-differentiation, and growth promotion in ccRCC.

We next searched for causal regulatory interactions connecting ccRCC somatic mutations to these transcriptional hubs, using a bi-directional extension to HotNet (‘TieDIE’) and identified a chromatin-specific sub-network (Figures 4a, S50–S52). TieDIE defines a set of transcriptional targets, whose state in the tumor cells is hypothesized to be influenced by one or more of the significantly mutated genes. The chromatin modification pathway intersects a wide variety of processes, including the regulation of hormone receptors (e.g. ESR1), RAS signaling via the SRC homolog (SHC1), immune-related signaling (e.g. NFκB1 and IL6)[26], transcriptional output (e.g. HIF1A, JUN, FOS, and SP1), BRCA1 function (via BAP1) and Beta-catenin (CTNNB1) and TGF-beta (TGFBR2) signaling via interactions with a SMARC-PBRM1-ARID1A complex. The complexity of these interactions reflects the potential for highly pleiotrophic effects following primary events in chromatin modification genes.

Figure Four

Genomically-altered pathways in ccRCC

(a) Alterations in chromatin remodeling genes were predicted to impact a large network of genes and pathways (larger implicated network in supplemental). Each gene is depicted as a multi-ring circle with various levels of data, plotted such that each ‘spoke’ in the ring represents a single patient sample (same sample ordering for all genes). ‘PARADIGM’ ring, bioinformatically inferred levels of gene activity (red, higher activity); ‘Expression’, mRNA levels relative to normal (red, high); ‘Mutation’, somatic event; center, correlation of gene expression or activity to mutation events in chromatin-related genes (red, positive). Protein-protein relationships inferred using public resources. (b) For the PI3K/Akt/mTOR pathway (altered in ~28% of tumors), the MEMo algorithm identified a pattern of mutually exclusive gene alterations (somatic mutations, copy alterations, and aberrant mRNA expression) targeting multiple components, including 2 genes from the recurrent amplicon on 5q35.3. The alteration frequency and inferred alteration type (blue for inactivation, and red for activation) is shown for each gene in the pathway diagram.

The mutations in the chromatin regulators PBRM1, BAP1, and SETD2 were differentially associated with altered expression patterns of large numbers of genes when compared to samples bearing a background of VHL mutation (Table S18–S21, Figure S53). Each chromatin regulator had a distinct set of downstream effects, reflecting diverse roles for chromatin remodeling in the transcriptome.

Additionally, an unsupervised pathway analysis using the MEMo algorithm[22], identified mutually exclusive patterns of alterations targeting multiple components of the PI3K/Akt/mTOR pathway in 28% of the tumors (Figures 4b, Table S22). Interestingly, the altered gene module included two genes from the broad amplicon on 5q35.3: GNB2L1/RACK1 and SQSTM1/p62. Both these genes have previously been associated with activation of PI3K signaling[23,24]. Furthermore, mRNA expression levels of these two genes were correlated with both DNA copy number level increases and alteration status of the PI3K pathway (Figures S54–S55). The mutual exclusivity module also includes frequent over-expression of EGFR, which correlates with increased phosphorylation of the receptor (Figure S56), and which has been previously associated with Lapatinib response in ccRCC[25].

Correlations with survival

Where unsupervised analyses had indicated that common molecular patterns were associated with patient survival, we sought to further define molecular prognostic signatures at the levels of mRNA, miRNA, DNA methylation, and protein. Data were divided into ‘discovery’ (N=193) and ‘validation’ (N=253) sets and platform-specific signatures were defined using Cox analyses[26]. Kaplan-Meier analysis for each signature showed statistically significant associations with survival in the validation subset (Figures 5a and S57). Multivariate Cox analyses, incorporating established clinical variables, showed that the mRNA, miRNA, and protein signatures provided additional prognostic power (Table S23). In addition, these signatures could provide molecular clues as to the drivers of aggressive cancers.

Figure Five

Molecular correlates of patient survival involve metabolic pathways

(a) Sample profiles were separated into discovery and validation subsets, with the top survival correlates within the discovery subset being defined for each of the four platforms examined (mRNA, microRNA, protein, DNA methylation). Kaplan-Meier plots show results of applying the four prognostic signatures to the validation subset, comparing survival for patients with predicted higher risk (red, top third of signature scores), lower risk (blue, bottom third), or intermediate risk (gray, middle third); successful predictions were observed in each case. (b) When viewed in the context of metabolism, the molecular survival correlates highlight a widespread metabolic shift, with tumors altering their usage of key pathways and metabolites (red and blue shading representing the correlation of increased gene expression with worse or better survival respectively, univariate Cox based on extended cohort). Worse survival correlates with up-regulation of pentose phosphate pathway genes (G6PH, PGLS, TALDO and TKT), fatty acid synthesis genes (ACC and FASN), and PI3K pathway enhancing genes (miR-21). Better survival correlates with up-regulation of AMPK complex genes, multiple Krebs cycle genes, and PI3K pathway inhibitors (PTEN, TSC2). Additionally, specific promoter methylation events, including hypermethylation of PI3K pathway repressor GRB10, associate with outcome. (c) Heat map of selected key features from the metabolic shift schematic (b) demonstrating coordinate expression by stage at DNA methylation, RNA, and protein levels (data from validation subset).

Top protein correlates of worse survival included reduced AMP-activated kinase (AMPK) and increased acetyl-CoA carboxylase (ACC) (Figure S58). Together, down-regulation of AMPK and up-regulation of ACC activity contribute to a metabolic shift towards increased fatty acid synthesis.[27] A metabolic shift to an altered use of key metabolites and pathways was also apparent when considering the full set of genes involved in the core metabolic processes, including a shift towards a “Warburg effect”-like state (Figure 5b). Poor prognosis correlated with down-regulation of AMPK complex and the Krebs cycle genes, and with up-regulation of genes involved in the pentose phosphate pathway (G6PD, PGLS, TALDO, TKT) and fatty acid synthesis (FASN, ACC).

Examination of potential genetic or epigenetic drivers of a glycolytic shift led us to identify methylation events involving miR-21 and GRB10, with decreased promoter methylation of each gene (thereby higher expression) being associated with worse or better outcome, respectively (Figure 5B, Figure S59, Table S24). Both genes regulate the PI3K pathway: miR-21 is inducible by high glucose levels and down-regulates PTEN[20]; while the tumor suppressor Grb10 negatively regulates PI3K and insulin signaling.[28] Promoter methylation of miR-21 and GRB10 were coordinated with their mRNA expression patterns, as well as with the mRNA expression of other key genes and protein expression in the metabolic pathways (Figures 5C and S60). In addition to the PI3K pathway (Figures 5B and S61), molecular survival correlations involved several pro-metastatic matrix metalloproteinases (Figure S62).

Discussion

Our study sampled a single site of the primary tumor, in a disease with a potentially high level of tumor heterogeneity[8]. The extent to which convergent evolutionary events are a common theme in ccRCC remains to be determined, but may indicate that critical genes will be represented across the tumor landscape for an individual mass. In general, the large sample size appeared to overcome the intrinsic challenges of studying a genetically complex disease, revealing rare variants at rates similar to what has been described previously[29]. The samples, taken from primary tumor specimens, were reflective of patients fit for either definitive or cytoreductive nephrectomy, while future work could explore the genomic landscape of metastatic lesions.

Pathway and integrated analyses highlighted the importance of the well-known VHL/HIF pathway, the newly emerging chromatin remodeling/histone methylation pathway, and the PI3K/AKT pathway. The observation of chromatin modifier genes being frequently mutated in ccRCC strongly supports the model of nucleosome dynamics providing a key function in renal tumorigenesis. Although the mechanistic details remain to be defined as to how such modulation promotes tumor formation, the data presented here revealed alterations in DNA methylation associated with SETD2 mutations. As an epigenetic process that can potently alter many transcriptional outputs, these mutational events have the potential to alter the landscape of the tumor genome via altered expression of global sets of genes and genetic elements. Molecular correlates of patient survival further implicated PI3K/AKT as having a role in tumor progression, involving specific DNA methylation events. The PI3K/AKT pathway presents a strong therapeutic target in ccRCC, supporting the potential value of MTOR and/or related pathway inhibitor drugs for this cancer.[30,31]

Cross-platform molecular analyses indicated a correlation between worsened prognosis in patients with ccRCC and a metabolic shift involving increased dependence on the pentose phosphate shunt, decreased AMPK, decreased Krebs cycle activity, increased glutamine transport and fatty acid production. These findings are consistent with the isotopomer spectral analysis of a pair of VHL −/− clear cell kidney cancer cell lines, both of which were notably derived from patients with aggressive, metastatic disease, which revealed a dependence on reductive glutamine metabolism for lipid biosynthesis.[32] The metabolic shift identified in poor prognosis ccRCC remarkably mirrors the Warburg metabolic phenotype (increased glycolysis, decreased AMPK, glutamine dependent lipogenesis) identified in type 2 papillary kidney cancer characterized by mutation of the Krebs cycle enzyme, fumarate hydratase.[32] Further studies to dissect out the role of the commonly mutated chromosome 3 chromatin remodeling genes, PBRM1, SETD2 and BAP1, in ccRCC tumorigenesis and their potential role in the metabolic remodeling associated with progression of this disease will hopefully provide the foundation for the development of effective forms of therapy for this disease.

Methods Summary

Specimens were obtained from patients, with appropriate consent from institutional review boards. Using a co-isolation protocol, DNA and RNA were purified. In total, 446 patients were assayed on at least one molecular profiling platform, which platforms included: (1) RNA sequencing; (2) DNA methylation arrays; (3) miRNA sequencing; (4) Affymetrix SNP arrays; (5) exome sequencing; and (6) reverse phase protein arrays. As described above and in the Supplemental Methods, both single platform analyses and integrated cross-platform analyses were performed.