VHL Deficiency Drives Enhancer Activation of Oncogenes in Clear Cell Renal Cell Carcinoma

Xiaosai Yao1,2, Jing Tan3,4, Kevin Junliang Lim5, Joanna Koh1, Wen Fong Ooi1, Zhimei Li3, Dachuan Huang3, Manjie Xing1,5,6, Yang Sun Chan1, James Zhengzhong Qu1, Su Ting Tay5, Giovani Wijaya3, Yue Ning Lam1, Jing Han Hong5, Ai Ping Lee-Lim1, Peiyong Guan3, Michelle Shu Wen Ng2, Cassandra Zhengxuan He1, Joyce Suling Lin1, Tannistha Nandi1, Aditi Qamra1,7, Chang Xu5,8, Swe Swe Myint3, James O. J. Davies9, Jian Yuan Goh1, Gary Loh1, Bryan C. Tan10, Steven G. Rozen5, Qiang Yu1, Iain Bee Huat Tan1,11, Christopher Wai Sam Cheng12, Shang Li5, Kenneth Tou En Chang13, Puay Hoon Tan14, David Lawrence Silver9, Alexander Lezhava15, Gertrud Steger16, Jim R. Hughes9, Bin Tean Teh2,3,5,8,17, and Patrick Tan1,5,8,17

iNtRODUctiON Clear cell renal cell carcinoma (ccRCC) is the most common subtype of kidney cancer, with 338,000 new cases in 2012 worldwide (1). With most ccRCCs being radiochemoresistant, patients with metastatic ccRCC exhibit dismal 8% five-year overall survival (2). Although targeted therapies inhibiting angiogenesis and mTOR pathways can lead to initial tumor control, most patients develop resistance in less than a year (3,4). A better understanding of ccRCC molecular dependencies and vulnerabilities is thus needed to develop new therapeutic strategies for patients who fail standard-of-care treatment.
Loss of the von Hippel-Lindau (VHL) tumor suppressor is a defining feature of ccRCC (5,6). When partnered with inactivation of additional tumor suppressors (Pbrm1, Bap1, Trp53, Rb1, and/or Cdkn2a) and/or activation of oncogenes Delineating the global ccRCC cis-regulatory landscape may also identify master regulators involved in tissue-specific disease processes. Compared with promoters that are largely cell-type invariant, distal enhancers integrate multiple lineage-and context-dependent signals, catering to the specialized needs of diverse cell types and diseases (32,33). In cancer, such master regulators are frequently located near "superenhancers" or "stretch-enhancers" marked by long stretches of H3K27ac signals (34,35). For example, subtype-specific genomic alterations such as EGFRvIII in glioblastoma (36) and EWS-FLI in Ewing sarcoma (37) induce de novo enhancers, causing reactivation of developmental master regulators required for self-renewal and lineage specification (36). Although VHL inactivation has been shown to modulate protein levels of different histone modifiers (e.g., KDM5C/ JARID1C, ref. 38 41), the impact of these protein alterations at specific epigenomic loci remains unclear. Moreover, previous studies profiling histone modifications in ccRCCs have also been limited by small sample sizes (2 cases; ref. 42), reliance on in vitro systems, and the lack of long-range interactome data and functional enhancer testing to accurately assign cognate enhancer targets.
In this study, we establish the most comprehensive collection of ccRCC histone profiles to date, annotating the precise genomic locations of altered promoters, enhancers, and superenhancers in ccRCC. Using isogenic cell lines with or without wild-type VHL, we further demonstrate that besides its well-defined role in oxygen sensing, VHL also safeguards the chromatin landscape; its loss induces tumor-specific enhancer gains around ccRCC hallmark genes such as angiogenic and metabolic targets through the stabilization of HIF2α/HIF1β (ARNT) heterodimers and recruitment of p300 histone acetyltransferase (EP300). One important target of epigenetic activation is ZNF395, a master regulator of ccRCC tumorigenesis. Taken collectively, our results reveal an epigenetic framework by which the major ccRCC-specific driver mutation, VHL, induces de novo enhancers, contributing to oncogenic transcription.
Specific histone modifications can distinguish different categories of functional regulatory elements-H3K4me3 is generally associated with promoters, H3K4me1 with enhancers, and H3K27ac with active elements (33,43). Integrating signals from three histone marks and GENCODE v19 annotated transcription start sites (TSS), we defined active promoters as H3K27ac + /H3K4me3 + /±2.0 kb TSS regions, and distal enhancers as H3K27ac + /H3K4me1 + regions not overlapping with promoters. Focusing on epigenomic events specific to somatic cancer cells, we derived cell lines from five primary tumors and, combined with the commercial lines, excluded peaks not found in any of the cell lines to reduce confounding effects from stromal cells. On average, we observed 80% overlap of ChIP-seq peaks between primary tumors and matched lines ( Supplementary Fig. S1B). Using these criteria, we identified 17,497 putative promoters and 66,448 putative enhancers (Fig. 1A), numbers comparable with previous studies in other tumor types (43)(44)(45). The numbers of defined promoters and enhancers reached saturation after 4 and 16 samples, respectively, suggesting that a sample size of 20 (10 tumor/ normal pairs) is sufficiently powered to discover the majority of cis-regulatory elements in ccRCC (Supplementary Fig.  S1C and S1D). Principal components analysis (PCA) using the first two components of global H3K27ac intensities at promoters or enhancers (representing 83% and 64% of total variance, respectively; Supplementary Fig. S1E and S1F) successfully separated normal and tumor samples, indicating that genome-wide pervasive alterations in cis-regulatory elements are a salient feature of ccRCC (Fig. 1B).
We performed differential analysis to identify altered promoters and enhancers. To define gained or lost regions, we applied a fold difference of H3K27ac RPKM ≥ 2, an absolute difference ≥ 0.5, and for greater stringency no alterations in the reverse direction in the remaining tumor/normal pairs (see Methods and Supplementary Fig. S1G for distribution of altered elements by number of patients). At the threshold of ≥5/10 patients, 80% of the altered regions achieved statistical significance (q-value < 0.1, paired t test, with Benjamini-Hochberg correction; Supplementary Fig. S1H), and at this same threshold, the increase in the fraction of samples meeting statistical significance reached a saddle point ( Supplementary Fig. S1I). Applying these criteria, we obtained a high-confidence and comprehensive set of 4,719 gained promoters, 592 lost promoters, 4,906 gained enhancers, and 5,654 lost enhancers (  Table S4). Representative regions are presented in Supplementary Fig. S2.
Supporting the reliability of these data, gained promoters and enhancers exhibited increased chromatin accessibility measured by higher FAIRE-seq signals (46) in tumor tissues than normal tissues, respectively (P < 0.0001) and also decreased November 2017 CANCER DISCOVERY | OF4 H3K27ac + /H3K4me3 + /±2 kb TSS   1D). Interestingly, we also noted elevated expression of long noncoding RNAs (48) adjacent to gained promoters and enhancers in tumor tissues compared with normal tissues (P < 0.0001, respectively). Lastly, we confirmed that many of our cis-regulatory elements involved regions previously implicated in ccRCC; for example, we observed gains of H3K27ac signals and enrichment of H3K4me1 at a distal enhancer of CCND1 overlapping with an RCC susceptibility locus (rs7105934; refs. 21 and 49; Fig. 1E). Our ability to rediscover this important enhancer in our unbiased profiling supports our data reliability.

Tumor-Specific Enhancers Are Associated with Hallmarks of ccRCC
To identify genes modulated by the tumor-specific regulatory elements, we assigned enhancers using three approaches. The first approach utilized predefined linear proximity rules involving a set of highly confident genes (GREAT algorithm;  ref. 50; Supplementary Table S5). MSigDB pathway analysis using GREAT-assigned genes revealed that gained enhancers exhibit a highly significant RCC-specific signature compared with gained promoters (enhancer q-value = 3.2 × 10 −26 ; promoter q-value = 1.5 × 10 −1 , binomial FDR; Fig. 2A). Although gained promoters were involved in general cancer processes (e.g., cell cycle, transcription, and RNA metabolism; see Supplementary Table S6 for a complete list of promoter pathways), gained enhancers were enriched in disease-specific features of ccRCC, including HIF1α network activity, proangiogenic pathways (platelet activation and PDGFRβ signaling), and SLC-mediated transmembrane transport ( Fig. 2A  expression of the nearest ccRCC lncRNA are compared between normal and tumor tissues. ***, P < 0.001, two-sided t test. E, Shown are tracks of histone ChIP-seq (H3K27ac, H3K4me1, and H3K4me3) and RNA sequencing (RNA-seq) at the CCND1 locus in a tumor-normal pair of patient 40911432. The histone ChIP-seq profiles of normal adult kidney tissue from the Epigenome Roadmap are displayed above the normal tissue generated by Nano-ChIPseq for comparison. A cell line was derived from the tumor tissue, and its histone profile is displayed below its matching primary tissue. This enhancer is known to interact with the CCND1 promoter from a previous study (21) and is situated close to an RCC susceptibility SNP rs7105934 (49).
We also used a second method of enhancer-gene assignment based on correlations between H3K27ac signals and expression of genes within the same topologic associated domain (TAD; ref. 34). Using a q-value of <0.05 based on Spearman correlation, we assigned 2,311 gained enhancers to 2,186 protein-coding targets (Supplementary Table S9). Reassuringly, H3K27ac signals of many gained enhancers were highly correlated with gene expression of their putative target genes. For example, H3K27ac levels of a VEGFA enhancer exhibited high correlation with VEGFA gene expression (r = 0.83, Spearman correlation), whereas H3K27ac signals of an SLC2A1 enhancer were highly correlated with SLC2A1 gene expression (r = 0.72, Spearman correlation; Fig. 2B; Supplementary Fig. S3D). Similar to the GREAT approach, the TAD correlation approach also highlighted hypoxia (Krieg_Hypoxia_not_via_KDM3A, FDR q-value = 7 × 10 −120 ) and metabolism (Chen_Metabolic_Syndrome_ Network, FDR q-value = 2 × 10 −91 ) as highly enriched pathways (Supplementary Table S10).
Third, to independent validate the GREAT and TAD approaches in the specific context of ccRCC, we experimentally explored the interactome of ccRCC tumor-specific enhancers by performing Capture-C assays (56    both high-resolution (down to single Kb resolution) and high-throughput interrogation of user-defined regions (a usual working range of 10-500 regions). We designed probes against a subset of 56 gained enhancers and examined their interactions with protein-coding genes in 786-O cells. Each gene-enhancer pair revealed by Capture-C was further filtered by correlations between gene expression and H3K27ac levels (q-value < 0.05). The 56 gained enhancers were paired with 36 protein-coding genes ( Supplementary Table S11)-of these, 58% were predicted by GREAT, and 80% by gene correlations within TADs. The median distance of interactions detected by Capture-C was 16 kb, and 83% of the interactions fell within a 100-kb window ( Supplementary Fig. S3E). As a visual example, Capture-C confirmed interactions between VEGFA enhancer and the VEGFA TSS, spanning a distance of ∼100 kb (Fig. 2B), and interactions between the SLC2A1 enhancer and its promoter (Fig. 2C). Taken collectively, these findings highlight the disease-specific nature of enhancer elements (33) and an important role for enhancer malfunction in modulating ccRCC pathology.

Tumor Superenhancers Identify ZNF395 as a Master Regulator of ccRCC Tumorigenesis
The importance of enhancers in ccRCC led us to examine the landscape of "superenhancers" or "stretch-enhancers"dense clusters of enhancers located near master regulators of Research.
Because current therapeutic targets in kidney cancer are limited to angiogenesis and mTOR pathways (3), we sought to examine these less-understood genes uncovered    Fig. 3E). Furthermore, among 12 types of cancer profiled by TCGA, ZNF395 was exclusively overexpressed in ccRCC tumors, consistent with the proposed lineage-and diseasespecific nature of superenhancers (Fig. 3F).
No study to date has functionally tested the tumorigenic requirement of ZNF395 in ccRCC or any other cancer type. We validated ZNF395's tumor-promoting effect using individual shRNA clones (Supplementary Fig. S4E and S4F). Two independent ZNF395 shRNA clones drastically decreased in vitro colony formation (Fig. 3G) and cell viability (  S4G). In vivo, tumor formation studies in mouse xenograft models revealed marked tumor suppression by ZNF395 depletion (Fig. 3J). Knockdown of ZNF395 led to elimination of A-498 tumors up to day 74, when tumors in the control group began to exceed the size limits imposed by institutional animal protocols. Similarly, ZNF395 depletion significantly slowed in vivo tumor growth of 786-O cells (Fig. 3J). Taken together, we showed the indispensable role ZNF395 plays in ccRCC tumorigenesis.
Even though gained enhancers were expected to show only depletion after VHL restoration, changes in H3K27ac levels were bidirectional (Fig. 4A). However, only gained enhancers with H3K27ac depletion were uniquely active in VHL-mutated ccRCC cell lines (786-O, A-498, and 12364284) compared with VHL-wild-type ccRCC cells (86049102L), normal kidney cell lines (PCS-400, HK2, and HKC-8), and 31 other cell lines of various cancer types (Fig. 4B). The lack of H3K27ac signals in normal kidney cell lines argues against tissue lineage as the dominant contributor to the high H3K27ac ChIP-seq signals seen in ccRCC cell lines. On the other hand, gained enhancers with H3K27ac enrichment after VHL restoration showed high activity across multiple cancer types, suggesting that these enhancers are not unique to ccRCC (Fig. 4B).
Furthermore, only gained enhancers showing H3K27ac depletion after VHL restoration were significantly associated with a concomitant downregulation of gene expression of their putative targets in both 786-O and 12364284 cells, whereas enhancers gained in primary ccRCCs and further H3K27acenriched after VHL restoration did not lead to significant gene upregulation on a global level ( Fig. 4C; Supplementary  Fig. S6D). These results suggest that the former enhancers (H3K27ac depletion) are likely to represent ccRCC-and VHLspecific epigenomic alterations, whereas the latter enhancers (H3K27ac enrichment) are likely to represent signify generic, compensatory mechanisms in response to VHL restoration.
Combining data from multiple lines, a total of 1,564 enhancers were depleted by VHL restoration in ≥1 cell line, representing almost a third (32%) of all gained enhancers identified in primary ccRCC tumors (Supplementary Table    Fig. S7C). Collectively, pathway analysis of enhancers depleted in ≥2 cell lines highlighted direct p53 effectors, integrin-linked kinase signaling, and HIF1α transcription factor networks as the top five pathways (Supplementary Table S14), covering genes such as EGFR (Fig. 4D), CCND1 (Fig. 4E), ITGB3 (Fig.  4F), VEGFA (Fig. 4G), SLC2A1 (Supplementary Fig. S7D), and HK2 ( Supplementary Fig. S7E). These results support a major role for VHL loss in ccRCC enhancer malfunction, even in the presence of other driver mutations.
We also examined whether other histone marks were concomitantly altered with H3K27ac marks. We found a surprisingly high degree of correlation between H3K27ac and H3K4me1 in response to VHL restoration in both 786-O  Supplementary Fig. S7I), H3K27me3 levels remained low at gained enhancers even after VHL restoration (Fig. 4H). These findings suggest that VHL restoration may result in a loss of enhancer identity by codepletion of H3K27ac and H3K4me1, but not a formal transition to a poised enhancer state that would have retained H3K4me1 but acquired H3K27me3.

HIF2`-HIF1a Heterodimers Are Enriched at VHL-Responsive Enhancers
We sought to investigate which transcription factors might mediate VHL-dependent chromatin remodeling at gained enhancers. Beginning with the primary ccRCC dataset, we looked for enrichment of trans-regulators in gained enhancers over lost enhancers. Using HOMER (65), we found that the top enriched motifs were the AP1 family, ETS family, and NF-κB-p65-Rel and HIF1α/2α motifs (Fig. 5A, full list of motifs in Supplementary Table S15). For subsequent in vitro validation, we chose c-Jun as a representative AP1 family member because of its activation in ccRCC (66) and ETS1 as an ETS family representative because of its known interaction with HIF2α (67), but acknowledge that other family AP1 and ETS family members may play a role in ccRCC. Immunoblotting of c-Jun, Research.
on October 26, 2017. © 2017 American Association for Cancer cancerdiscovery.aacrjournals.org Downloaded from ETS1, and NF-κB-p65 showed variable protein expression in both normal and tumor cell lines, but expression of HIF1α and HIF2α restricted to tumor cells only (Fig. 5B). HIF2α was expressed in a higher proportion of ccRCC cell lines than HIF1α (Fig. 5B). We further examined gene expression of these transcription factors in the TCGA cohort and found that ETS1, RELA (subunit of NF-κB-p65), and HIF2a were significantly overexpressed in tumors compared with normal tissues, with a range of tumor-association expression patterns similar to variations in ccRCC lines (Supplementary Fig. S8A).
To further investigate chromatin occupancy of these factors, we generated ChIP-seq binding profiles of c-Jun, ETS1, and NF-κB cells and reexamined HIF2α, HIF1α, and HIF1β binding profiles from the previous literature (21,30), all performed in 786-O cells. Of note, because 786-O cells have lost endogenous HIF1α expression through genomic deletion, the HIF1α ChIP-seq was performed on 786-O cells genetically manipulated to reexpress HIF1α protein (30). ChIP-seq results showed that all six transcription factors exhibited increased occupancy at gained enhancers compared with lost enhancers, validating the HOMER predictions (Fig. 5C).
To determine which of these transcription factors might be directly dependent on VHL, we then compared their protein expression in VHL-mutated isogenic cell lines with and without wild-type-VHL restoration. As shown in Fig. 5D, VHL restoration consistently downregulated HIF2α expression in both 786-O and 12364284 cell lines, but protein levels of other factors displayed contrasting trends between the two cell lines, implying that among the six factors examined, HIF2α protein expression was the most VHL dependent. Indeed, supporting an important role for HIF2α in VHLdependent enhancer remodeling, only HIF2α and HIF1β were significantly enriched at enhancers showing VHLdependent H3K27ac depletion (Fig. 5E). Moreover, among all known motifs in the HOMER database, HIF2α was the most enriched motif at VHL-responsive enhancers exhibiting H3K27ac depletion (P = 1 × 10 −11 ; Supplementary Table S16).
In contrast, HIF1α was not enriched at enhancers showing H3K27ac depletion (Fig. 5E). Despite sharing many binding sites with HIF2α, HIF1α predominantly localized to promoter-proximal regions, whereas HIF2α frequently occupied introns and intergenic regions in 786-O cells ( Supplementary  Fig. S8B), consistent with a promoter-centric occupancy of HIF1α and an enhancer-centric occupancy of HIF2α (Fig. 5F). Gained enhancers displayed a HIF2α occupancy twice that of tumor-specific promoters (P < 1 × 10 −16 , proportions test) in 786-O cells, suggesting that HIF2α may play a greater role in regulating enhancers than promoters.
To extend these HIF1α and HIF2α occupancy-pattern findings to a system that expresses endogenous levels of both factors, we then performed HIF1α and HIF2α ChIP-seq in 40911432 ccRCC cells, which abundantly coexpress both HIFα subunits (Fig. 5B). Similar to 786-O, in 40911432 cells, HIF1α showed a preferential occupancy at promoter-proximal regions, whereas a large proportion of HIF2α was found in distal regions (introns and distal intergenic regions; Supplementary Fig. S8C). A higher proportion of HIF1α binding sites overlapped with gained promoters than HIF2α (68% of HIF1α vs. 41% of HIF2α, P = 0.002, proportions test; Fig. 5G). Conversely, a higher proportion of HIF2α binding sites overlapped with gained enhancers than HIF1α (29% of HIF1α vs. 51% of HIF2α, P < 2.2 × 10 −16 , proportions test). HIF2α's preferential occupancy at enhancers was further substantiated by its higher enrichment at enhancers showing H3K27ac depletion after VHL restoration than HIF1α (Fig. 5H). Specific examples of VHL-responsive enhancers bound exclusively by HIF2α but not HIF1α included an enhancer near UBR4 (Fig. 5I) and a superenhancer near CMIP (Fig. 5J). Therefore, even in HIF1α/HIF2α coexpressing ccRCC cells, these results suggest that HIF2α plays a greater role in VHL-mediated enhancer remodeling than HIF1α.
We sought to establish a causal link between HIF2α-bound enhancers and control of gene expression. We performed CRISPR-mediated genomic depletion of the ZNF395 enhancer region with the highest HIF2α peak (Fig. 6G). All four clones with the homozygous deleted ZNF395 enhancer consistently downregulated their ZNF395 expression compared with clones with the intact enhancer (P < 0.05), providing evidence that ZNF395 expression is epigenetically controlled by this HIF2α-HIF1β-bound enhancer (Fig. 6G). Taken together, these results indicate that that HIF2α is likely an important mediator of VHL-driven enhancer remodeling.

VHL Restoration Reduces p300 Recruitment but Preserves Promoter-Enhancer Interactions
Finally, we sought to investigate why VHL restoration caused a decrease in H3K27ac levels. Previous pulldown assays have reported that both HIF2α and HIF1β can interact with histone acetyltransferase p300 (68)(69)(70)(71). Indeed, p300 frequently marks enhancers (43) and is thought to be recruited     Distance from HIF2α ChIP-seq summits by tissue-specific transcription factors (72). However, chromatin profiles of p300 have not been previously established in kidney cancer cell lines, so the contribution of p300 in shaping enhancers in ccRCC remains unclear. Therefore, we performed p300 ChIP-seq in 786-O cells and confirmed its enrichment at gained enhancers over lost enhancers (Fig. 7A). Comparing p300 ChIP-seq with HIF2α ChIP-seq yielded a surprisingly high degree of overlap between HIF2α and p300 (96%), even more than that of HIF2α and HIF1β (89%; Fig.  7B and C). In contrast, other transcription factors such as c-Jun, ETS1, and NF-κB did not exhibit such a high degree of overlap (≤60%; Fig. 7B).
We compared p300 binding at tumor enhancers with and without VHL. Despite increased p300 protein levels in 786-O cells after VHL restoration (Fig. 7D), binding of p300 decreased across all four enhancers examined (Fig. 7E). HIF2a depletion by siRNA knockdown also decreased p300 recruitment (Fig. 7F), suggesting that loss of HIF2α may interfere with p300 recruitment.
We investigated whether VHL restoration and the subsequent loss of p300 binding disrupted promoter-enhancer interactions. We performed Capture-C of enhancer regions in paired 786-O cell lines with and without VHL restoration. Surprisingly, Capture-C interactions showed a relatively high Research.
on October 26, 2017. © 2017 American Association for Cancer cancerdiscovery.aacrjournals.org Downloaded from correlation between VHL-deficient and VHL-restored 786-O cells at VHL-responsive regions (r = 0.74, Pearson correlation), even higher than correlations observed at non-VHL-responsive regions (r = 0.57, Pearson correlation; Fig. 7G). As a visual example, interactions between the VEGFA promoter and enhancer were intact even after VHL restoration (Fig. 7H), indicating that loss of enhancer activity is likely insufficient to dissociate promoter-enhancer interactions. Furthermore, many of these promoter-enhancers were lineage specific; for example, the interaction between SLC2A1 enhancer with its promoter was not detected in KATOIII, a gastric cancer cell line (Supplementary Fig. S9). Therefore, promoter-enhancer interactions often preexist in kidney cells, frequently in a tissue-specific manner.

DiscUssiON
Understanding epigenomic alterations and their genetic origin can identify new disease mechanisms (34), vulnerabilities (73,74), and therapeutic strategies (75)(76)(77). Through comprehensive profiling of histone modifications in primary normal-tumor pairs and cell lines, we generated a compendium of ccRCC-associated promoters and enhancers. Our study demonstrates that the most frequent ccRCC mutational event-VHL inactivation-leads to genome-wide enhancer and superenhancer remodeling, which directly imparts ccRCC hallmarks including angiogenesis and metabolic reprogramming. ZNF395, epigenetically controlled by a VHL-responsive superenhancer, emerged as a crucial regulator of ccRCC tumorigenesis.
Our work has three main advances. First, to our knowledge, this is the most comprehensive atlas of histone profiles in ccRCC and will likely provide an invaluable resource to the ccRCC field. Using high-resolution multiplexed interactome data (Capture-C, ref. 56) and H3K27ac-expression correlation, we minimized ambiguity in enhancer assignment and further confirmed the dependency of enhancers on VHL/HIF status by reporter assays. Second, using isogenic cell lines, we show that VHL loss contributes significantly to enhancer remodeling. Even though another mutation in ccRCC, SETD2, can mediate widespread increases in chromatin accessibility (46) and DNA hypomethylation (78), its relatively low mutation frequency at ∼10% in all ccRCC tumors (78) cannot explain epigenetic abnormalities in the vast majority of SETD2-wild-type tumors. Lastly, an examination of somatically altered superenhancers enabled us to identify a master regulator crucial to the pathogenesis of ccRCC, ZNF395. Even though ZNF395 overexpression in ccRCC has been previously reported (79)(80)(81) and its proximity to a superenhancer was independently noted (42), our study is the first to pinpoint the specific VHL-dependent enhancer required for ZNF395 expression and to show ZNF395's indispensable functional role for ccRCC tumorigenesis in vitro and in vivo.
Our data suggest that mechanistically, loss of VHL stabilizes HIF2α occupancy at tumor-specific gained enhancers, which in turn recruits histone acetyltransferase p300 (28,82) to maintain H3K27 acetylation, upregulating expression of ccRCC-specific genes such as ZNF395 (Fig. 7I). Restoration of wild-type VHL resulted in codepletion of H3K27ac and H3K4me1 marks and thus abrogation of active enhancer identity at tumor-associated enhancers. Compared with the promoter-centric occupancy of HIF1α, HIF2α is predominantly found at enhancers, pointing toward a key difference between HIF1α and HIF2α. We also found that HIF2a siRNA knockdown specifically attenuates the activity of HIF2αbound enhancers/superenhancers. Interestingly, the majority of promoter-enhancer interactions remained largely unaltered by VHL status, suggesting that these promoter-enhancer couplings are largely stable and preformed. This is consistent with a recent report demonstrating that promoter-enhancer interactions remain largely unchanged between normoxia and hypoxia (29). Our study demonstrating VHL's impact on chromatin remodeling also suggests that other cancer genes with high tumor-type-specific mutational penetrance, such as BRAF in melanoma (83) and APC in colon cancer (84), may also act to modify cellular epigenomes to effect broad yet disease-specific cellular changes, despite these genes not being classic chromatin modifiers.
Besides VHL, mutations in other genes such as PBRM1, SETD2, ARID1A, SMARCA4, JARID1C, and KDM6A/UTX have been reported in ccRCC (47,85), and these are likely to augment VHL's core transcriptional effects (86), contributing to heterogeneity in disease phenotypes (87) and progression patterns (88). Using the example of 786-O cells, at least two other mutations in these cells may have a direct impact on chromatin-MLL3 (p.A3902G) and gain-of-function TP53 mutations (p.R248W). MLL3, a histone 3 lysine 4 methyltransferase, is directly responsible for formation of the H3K4me1 enhancer mark (89,90) and plays a critical role for enhancer regulation (91). Gain-of-function TP53 mutants also bind aberrantly to chromatin, especially near methyltransferases MLL1 and MLL2, potentially contributing to tumor growth via chromatin deregulation (92). Besides mutations, structural variants are also known to alter enhancers via enhancer hijacking (93) or copy-number gains (94) in other cancers. Given the multitude of driver and bystander mutations in ccRCC, it is thus unlikely that VHL alone can account for all epigenomic changes observed in this tumor type. Nevertheless, by integrating data across multiple ccRCC cell lines, our data suggest that VHL inactivation is likely to account for almost a third (32%) of all gained enhancer regions, supporting its role as a dominant driver of epigenetic abnormalities in ccRCC despite the presence of other genetic changes.
Our epigenetic maps contain a wellspring of both wellvalidated and uncharacterized targets that may contribute to ccRCC tumorigenesis. We found extensive enhancer gains around well-characterized hypoxia-related targets (ref.  (98), and adipogenesis (PLIN2,refs. 51,99). New targets revealed in this study include SMPDL3A, which could be another important ccRCC-specific oncogene given its role in lipid and cholesterol metabolism (61,100). Genes associated with lost superenhancers, which could only be identified with normal-tumor pairs, implicated potential tumor suppressors (EHF, MAL, GCOM1, and HOXB9) that warrant further investigation.
One notable finding from this epigenomic study is the tumorigenic requirement of ZNF395 in ccRCC. ZNF395 is also known as HDBP2 (101)  PPARγ2 to promote adipogenesis (103). ZNF395 has been shown to bind to the promoters of Huntington gene (101) and interferon-induced genes, and to cause upregulation of cancer-related genes (MACC1, PEG10, CALCOCO1, and MEF2C; ref. 104) and proangiogenic chemokines including IL6 and IL8 under hypoxia (105). It remains to be elucidated in future studies the precise mechanism contributing to ZNF395's tumorigenic role.
The enhancer landscapes profiled in this study have implications beyond ccRCC. The poorly perfused tumor core makes hypoxia a feature of virtually all solid tumors (106). MCF7 cells under hypoxia (but not normoxia) share similar H3K27ac profiles as 786-O (29). Although ZNF395 is highly expressed in ccRCC, its low basal expression can be upregulated upon hypoxia in other cancer types, including glioblastoma and skin cancer (104,105). Targeting ZNF395 or its downstream effectors in future studies may be therapeutically relevant to both ccRCC and other hypoxic solid malignancies. Direct targeting of ZNF395 using a peptidebased cancer vaccine is undergoing phase I trials in patients with sarcoma (107)(108)(109), opening up the possibility of using immunotherapy to target the extracellular fragments of nuclear master regulators. Our study suggests that initiating a similar trial in ccRCC may be worthwhile. Moreover, given the recent progress in targeting transcription factors using various modalities, including small molecules and stapled peptides (110)(111)(112), inhibitors of ZNF395 may provide an important therapeutic inroad for ccRCC treatment.

Patient Information
Fresh-frozen normal-tumor tissues were obtained from nephrectomy cases under approvals from institutional research ethics review committees and patient consent under institutional review board protocol 2010/735/B. Normal tissues were harvested from sites distant from the tumor. Refer to Supplementary Table S1 for detailed patient information.

Cell Lines
Commercial cell lines (786-O, A-498, HK2, and PCS-400) were purchased from ATCC. Cell lines were maintained in RPMI (Invitrogen) with 10% FBS with the exception of primary renal proximal tubule epithelial cells, PCS-400, which were maintained in Renal Epithelial Cell Basal Medium (ATCC). Cell line authentication was performed by short tandem repeat (STR) analysis (Cancer Science Institute of Singapore) in 2015 against publicly available STR profiles. Mycoplasma testing was performed using the MycoSensor PCR assay kit (Stratagene).

Establishment of Tumor-Derived Cell Lines from Primary Tumors
Tumor cells were disassociated from primary tumors by collagenase, seeded, and maintained in RPMI with 10% FBS. At 80% to 90% confluency, the cells were passaged at a 1:3 ratio. Cultured cells were considered to be successfully immortalized after 60 passages. Correct pairing of tumor tissues and cell lines was achieved by comparing the percentage identity of single nucleotide polymorphisms (SNP) based on targeted sequencing. All tumor-cell line pairs showed identities of > 90% whereas shuffling of pairing showed identities < 80%. Tumors and cell lines from 12364284 and 40911432 showed the same VHL mutations, but 86049102 tissue (named 86049102T) is VHL-mutant, whereas the cognate 86049102 cell line (named 86049102L) is VHL-wild-type. DNA/well of a 6-well plate using Lipofectamine 3000 (LifeTechnologies). A medium change was performed 10 to 16 hours after transfection. The supernatant from PlatA cells containing retroviruses was harvested 48 hours later and added to ccRCC cells, which were then selected with puromycin for 3 days after transduction.

Histone Nano-ChIP-seq
Nano-ChIP-seq was performed as previously described (113) with slight modifications. Fresh-frozen cancer and normal tissues were dissected using a razor blade to obtain ∼5 mg of tissue. The tissues were fixed in 1% formaldehyde for 10 minutes at room temperature. Fixation was stopped by addition of glycine to a final concentration of 125 nmol/L. Tissue pieces were washed 3 times with TBSE buffer. Pulverized tissues were lysed in 100 μL lysis buffer and sonicated for 16 cycles (30s on, 30s off) using a Bioruptor (Diagenode). The following antibodies were used: H3K27ac (ab4729, Abcam), H3K4me3 (07-473, Millipore), H3K4me1 (ab8895, Abcam), and H3K27me3 (07-449, Millipore). The total volume of immunoprecipitation was 1 mL and the amount of antibody used was 2 μg. The input DNA was precleared with protein G Dynabeads (Life Technologies) for 1 hour at 4°C and then incubated with antibodies conjugated protein G beads overnight at 4°C. The beads were washed 3 times with cold wash buffer. After recovery of ChIP and input DNA, whole-genome amplification was performed using the WGA4 kit (Sigma-Aldrich) and BpmI-WGA primers. Amplified DNA was digested with BpmI [New England Biolabs (NEB)]. After that, 30 ng of the amplified DNA was used with the NEBNext ChIP-seq library prep reagent set (NEB). ChIP-seq in cell lines was performed using the same Nano-ChIP-seq protocol described above but with 1 × 10 6 cells. Each library was sequenced to an average depth of 20 to 30 million raw reads on HiSeq2500 using 101-bp single end reads.

Histone ChIP-seq Analysis
Sequencing tags were mapped against the human reference genome (hg19) using Burrows-Wheeler Aligner (BWA-mem; ref. 114; version 0.7.10). Reads were trimmed 10 bp from the front and the back to produce 81 bp. Only reads with mapQ >10 and with duplicates removed by rmdup were used for subsequent analysis. Significant peaks were called using CCAT (P < 0.05; ref. 115). The strength and quality of immunoprecipitation were assessed using CHANCE (116).

Transcription Factor ChIP-seq
For each transcription factor, 3 × 10 7 cells were cross-linked with 1% formaldehyde for 10 minutes at room temperature and stopped by adding glycine to a final concentration of 125 nmol/L. Chromatin was extracted and sonicated to ∼500 bp (Vibra cell, SONICS). The following antibodies were used for chromatin immunoprecipitation: c-Jun (sc-1694, Santa Cruz Biotechnology), NF-κB p65 (sc-372, Santa Cruz Biotechnology), ETS1 (sc-350, Santa Cruz Biotechnology), HIF1α (610959, BD Biosciences), HIF2α (NB100-122, Novus Bio), and p300 (sc-585, Santa Cruz Biotechnology). The total volume of immunoprecipitation was 1.5 mL and the amount of antibody used was 15 μg. Input DNAs were precleared with protein G Dynabeads (LifeTechnologies) for 2 hours at 4°C and then incubated with antibody-conjugated protein G beads overnight at 4°C. The beads were washed 6 times with wash buffer at room temperature. At least 10 ng of the DNA was used with the NEBNext ChIP-seq library prep reagent set (NEB). Each library was sequenced to an average depth of 30 to 50 million reads on a HiSeq2500 using 101-bp single-end reads.

Capture-C
Capture-C was performed as previously described (56). Briefly, 1 × 10 7 cells were cross-linked by 2% formaldehyde, followed by lysis, homogenization, DpnII digestion, ligation, and de-cross-linking. DNA was sonicated using a Covaris to 150 to 200 bp to produce DNA suitable for oligo capture. A total of 3 μg of sheared DNA was used for sequencing library preparation (NEB). Enhancer sequences were double captured by hybridization to customized biotinylated oligos (IDT) and enriched with Dynabeads (Life Technologies). Captured DNA was sequenced to an average depth of 2 million reads per probe on the HiSeq Illumina platform using 150-bp paired-end reads.

Capture-C Analysis and Gene Assignment
Preprocessing of raw reads was performed to remove adaptor sequences (trim_galore, http://www.bioinformatics.babraham.ac.uk/ projects/trim_galore/), and overlapping reads were merged using FLASH (117). In order to achieve short read mapping to the hg19 reference genome, the resulting preprocessed reads were then in silico digested with DpnII and aligned using Bowtie (using p1, m2, best, and strata settings). Aligned reads were processed using Capture-C analyzer (118) to (i) remove PCR duplicates; (ii) classify subfragments as "capture" if they were contained within the capture fragment, "proximity exclusion" if they were within 1 Kb on either side of the capture fragment, or "reporter" if they were outside of the "capture" and "proximity exclusion" regions; and (iii) normalize read counts per 100,000 interactions in bigwig format. We additionally used the r3Cseq package (119) on the capture and reporter fragments to identify significant interactions of the viewpoint against a scaled background (q-value < 0.05). Gene assignment is defined by the overlap of significant Capture-C peaks with genes with start and end defined by GENCODE v19. Interactions were plotted using Epigenome Gateway v40.0.

Identification of Differentially Enriched Regions
Significant H3K27ac peaks called by CCAT were merged across all normal-tumor samples. The same was performed with H3K4me1 and H3K4me3 ChIP-seq data. TSS were based on GENCODE v19. Promoters were defined as regions of overlap between H3K27ac and H3K4me3 and also overlapping with ±2.0 Kb around the TSS. Enhancers were defined as regions of overlap between H3K27ac and H3K4me1 but not overlapping with promoters. To minimize stromal contamination, we performed further filtering using cell line data, where enhancers and promoters not overlapping with H3K27ac peaks in any of the cell lines were discarded. Wiggle files of window size 50 bp were generated using MEDIPs (120) from bam files. The inputsubtracted signal for each promoter or enhancer region was computed using bigWigAverageOverBed to yield reads per kilobase per million (RPKM). The RPKM of H3K27ac, H3K4me1, and H3K4me3 ChIP-seq from promoters and enhancers were corrected for batch effects using Combat. Tumor-specific regions were defined as regions that have a fold difference of ≥2 and a difference of 0.5 RPKM from patientmatched normal tissue. Normal regions were defined as regions that have a fold difference of ≤0.5, and a difference of −0.5 RPKM from the corresponding regions in patient-matched tumors. Recurrently gained regions were defined as gain in ≥5/10 patients and no loss in any patients. Recurrently lost regions were defined as loss in ≥5/10 patients and no gain in any patients. Statistical testing for each cisregulatory region was performed using paired t tests with Benjamini-Hochberg correction. The differential regions were visualized using NGSplot (121).

Identification of Superenhancer Regions
Superenhancer regions were identified using ROSE (ref. 35; with promoter excluded), using H3K27ac peak regions merged from all patients (both normal and tumor tissue). Wiggle files of window size 50 bp were generated using MEDIPs (120) from bam files. The input-subtracted signal for each superenhancer was computed using bigWigAverageOverBed (sum of reads over covered bases). The superenhancer regions were ranked by the average difference of normal-tumor H3K27ac ChIP-seq signals. Gained superenhancers were defined as regions that have average differential H3K27ac ChIPseq signals >0. Lost superenhancers were defined as regions that have average differential H3K27ac ChIP-seq signals <0.
Additional methods can be found in Supplementary Methods.