NF-κB–Dependent Lymphoid Enhancer Co-option Promotes Renal Carcinoma Metastasis

Cancers develop the ability to spread, or metastasize, from their site of formation to distant organs through a multistep selection process (1). A variety of mechanisms are known to increase the probability of successful completion of one or more metastasis steps (2, 3). However, the mechanisms by which cancer clones acquire metastatic traits have remained poorly understood. For example, genetic analyses of both model systems and patient samples have thus far failed to identify mutational signatures of metastasis (4–7). On the other hand, metastatic progression in various experimental and clinical contexts is associated with specific gene expression programs, indicating the existence of underlying gene regulatory mechanisms that provide cancer cells with metastatic capabilities (5). How these aberrant gene expression programs arise is unknown.

Worldwide, renal cell carcinomas (RCC) resulted in an estimated 338,000 new cancer diagnoses and 144,000 cancer deaths in 2012 (8). The most common form of RCC, accounting for 70% to 75% of the cases, is clear cell renal cell carcinoma (ccRCC; ref. 9). Genetically, ccRCCs are characterized by inactivation of the von Hippel–Lindau tumor suppressor (VHL), which is lost in up to 90% of ccRCCs (10, 11). VHL loss leads to normoxic stabilization of the hypoxia-inducible factors HIF1A and HIF2A (also known as EPAS1), which in turn drive the expression of multiple target genes involved in ccRCC development (12). Of the two HIF isoforms, HIF2A is more important for ccRCC progression (12), whereas HIF1A may have tumor-suppressive functions (13). In addition to supporting tumor initiation and early progression, HIF2A is also critical for ccRCC metastasis. For example, HIF2A drives the expression of the chemokine receptor CXCR4, which promotes metastatic progression in ccRCC and other tumor types (14–19). Furthermore, VHL loss can also lead to the activation of other canonical oncogenic pathways, such as NF-κB (20, 21), which is associated with ccRCC progression and resistance to drug-induced apoptosis (22–25). The VHL–HIF2A pathway has thus been unequivocally linked to ccRCC initiation and progression at multiple levels. Yet, the molecular mechanisms that shape HIF2A transcriptomic output during ccRCC evolution remain incompletely understood.

Transcriptional enhancers are distal DNA elements that play an important role in gene regulation (26). Enhancer activation landscapes are highly tissue-specific, implying that enhancers and enhancer clusters are critical determinants of cell fate and identity (27, 28). Abnormal enhancer activation has also been linked to various disease states, including cancer (29). Specifically, altered enhancer activity can modulate cancer risk and activate oncogenes (30–38), suggesting that aberrant enhancer usage may also be critical for cancer progression and metastasis. Recently, genome-wide chromatin profiling has also linked VHL loss and the consequent HIF2A stabilization to the maintenance of enhancer activation landscapes in ccRCC (39). However, whether and how specific enhancers contribute to metastatic cancer progression remains largely unknown.

To further understand the mechanisms of metastasis gene activation, we set out to test whether perturbed enhancer function underlies metastatic carcinoma progression. Combining high-throughput enhancer profiling by chromatin immunoprecipitation and sequencing (ChIP-seq) and chromatin conformation capture by Hi-C in metastasis models of ccRCC, a malignancy with significant metastatic potential (40), we identified enhancers associated with metastatic progression. Using CRISPR interference (CRISPRi) and CRISPR/Cas9 genome editing, we demonstrate the functional importance of distal enhancer elements for metastatic colonization in vivo. We further show how a specific metastasis-promoting enhancer cluster is activated by the coregulatory action of two signaling pathways, the ccRCC core driver HIF2A and the NF-κB pathway that acts via a tissue-restricted lymphoid enhancer element. Finally, cross-species analysis revealed that the metastatic ccRCC enhancer complement consists of multiple developmental enhancer programs co-opted from different normal cell lineages. We propose a model whereby metastatic transcriptional programs in carcinoma cells arise through cross-lineage co-option of developmentally distinct physiologic enhancer activation states.

To identify metastasis-associated enhancers (MAE), we performed ChIP-seq using an antibody for acetylated histone H3 lysine 27 (H3K27ac), a chromatin modification associated with active gene regulatory elements (41), in four different VHL-mutant models of ccRCC: the nonmetastatic human cell lines 786-O and OS-RC-2 and the corresponding metastatic subpopulations 786-M1A (M1A) and OS-LM1 (LM1), respectively (Fig. 1A). In immunocompromised mice, the M1A and LM1 cells form tumors that histologically resemble human ccRCCs and they display a substantially increased metastatic phenotype when compared with their parental populations (19). At the molecular level, these systems harbor clinically relevant mutations in well-characterized ccRCC genes (42), and their gene expression profiles resemble those seen in patients with ccRCC with poor clinical outcome (19). Recent evidence suggests that their response to targeted therapy is also similar to that seen in patients (43, 44). Approximately 25,000 enhancers were detected in each sample, of which a shared set of 314 H3K27ac-enriched and 161 H3K27ac-depleted regions showed conserved changes across the replicates and the two model systems (Fig. 1B and C; Supplementary Fig. S1A–S1C; Supplementary Tables S1 and S2). The majority of elements were located in introns and intergenic regions >10 kb from known transcription start sites (Fig. 1D; Supplementary Fig. S1D), supporting their putative role as enhancers. Enrichment of a second enhancer-associated chromatin mark, monomethylated histone H3 lysine 4 (H3K4me1), was also detected in the same enhancer regions (Fig. 1E; Supplementary Fig. S1E and S1F), suggesting that MAEs were not poised enhancers but rather represented de novo activated regions (41).

MAEs often appeared relatively close to each other (e.g., Fig. 1C). Indeed, the distribution of distances between both enriched and depleted MAEs followed a bimodal pattern with an overrepresentation of MAEs being <100 kb apart from each other, when compared with an expected genome-wide distribution of similarly sized groups of H3K27ac-enriched elements (Supplementary Fig. S2A and S2B). The clustering of MAEs suggested that metastatic programming could be linked to superenhancers, i.e., genomic regions with particularly high enrichment of H3K27ac that have been linked to the control of cell identity (27). We found no evidence for the de novo formation of superenhancers in the metastatic cells; however, MAEs tended to be constituent enhancers within super enhancer loci (Supplementary Fig. S2C).

To test the association between enhancer activity and gene expression, we performed RNA-sequencing (RNA-seq) mRNA profiling. In general, the expression of the most proximal gene correlated best with MAE activity (Fig. 1F). A subtler correlation was seen on the second most proximal gene, and no general effect was detected for genes further away (Fig. 1F). MAE clusters, defined as two or more MAEs proximal to the same gene, were more strongly associated with gene expression than single MAEs (Fig. 1G; Supplementary Fig. S2D), but superenhancer-associated MAEs were not more significantly correlated with gene expression when compared with other MAEs (Supplementary Fig. S2E). Pathway analysis performed on the MAE-associated gene set highlighted functional classes associated with cancer and tumor progression (Supplementary Fig. S3; Supplementary Tables S3 and S4). Importantly, the enhancer activation landscapes in our model systems strongly correlated with those present in a set of human ccRCC samples, confirming that MAE activation patterns were relevant in vivo (Supplementary Fig. S4A and S4B). To further test the clinical relevance of the MAEs identified in our experimental systems, we analyzed H3K27ac patterns and corresponding RNA-seq transcriptomic profiles in a patient cohort of 10 matched cases of normal kidney and ccRCC, 9 of which harbored confirmed VHL mutations (39). In general, H3K27ac-enriched MAEs showed stronger H3K27ac signal in ccRCCs when compared with normal kidney (Supplementary Fig. S4C). All 10 ccRCCs showed activation of a subset of H3K27ac-enriched MAEs when compared with matched normal kidney, but the precise MAE activation landscape varied among the cases (Supplementary Fig. S4D). Importantly, comparison of the genes near active and inactive MAEs, respectively, revealed a strong correlation between gene expression and proximal MAE activity (Fig. 1H), closely mirroring the data from model systems (Fig. 1F). Finally, the expression pattern of genes correlating with proximal MAE activation in our model systems was associated with poor patient outcome in the large clinical ccRCC data set of The Cancer Genome Atlas (TCGA; Fig. 1I; ref. 11). Thus, metastatic ccRCC populations show clustered enhancer activation patterns that correlate with clinically relevant gene expression traits.

To identify the functionally relevant regions within MAEs, we performed ChIP-seq using an antibody for the enhancer-associated transcriptional cofactor p300 (Fig. 2A; ref. 45). Of the 314 H3K27ac-enriched MAEs, 205 contained at least one p300 peak. We chose the two most statistically significant p300-containing MAEs, MAE-1 and MAE-2, for functional analysis by lentivirally delivered CRISPRi (46) in M1A cells (Fig. 2B and C). Additionally, we focused on MAE-126, as it belonged to the same MAE cluster as MAE-1 (Fig. 2D). ChIP-qPCR and ChIP-seq confirmed specific CRISPRi-mediated reduction in H3K27ac and p300 (Fig. 2E and F; Supplementary Fig. S5A–S5I). We then inoculated the CRISPRi cells into the tail vein of immunocompromised mice. Targeting MAE-1 or MAE-126 strongly reduced the lung colonization potential of M1A cells; MAE-2 inhibition showed a nonsignificant trend (Fig. 3A and B). Histologic analysis revealed that MAE inhibition had a significant negative effect on the number of lung metastatic colonies, indicating a role for MAE-1 and MAE-126 in the early steps of metastatic colonization at the distant site (Fig. 3C and D; Supplementary Fig. S6A). The functional relevance of MAE-1 was confirmed in an additional metastatic clone (Fig. 3E–H; Supplementary Fig. S6B and S6C). MAE-1 and MAE-126 inhibition did not affect cell proliferation in vitro (Fig. 3I), and targeting MAE-1 had only a small and statistically nonsignificant effect on subcutaneous tumor growth (Supplementary Fig. S6D), indicating phenotypic specificity for MAE-1/MAE-126 as a mediator of metastatic colonization. ChIP-qPCR confirmed efficient enhancer targeting in vivo (Supplementary Fig. S6E). Orthotopic implantation of cancer cells gave concordant results: MAE-1 inhibition did not affect growth in the kidney, but metastatic colonization in the lungs was reduced (Supplementary Fig. S6F–S6I). Clinical ccRCC specimens also demonstrated evidence of MAE-1/MAE-126 activity (Supplementary Fig. S6J). Collectively, these experiments demonstrate that specific MAEs are functional mediators of metastatic colonization.

The human genome is organized into topologically associating domains (TAD) that are thought to limit the range over which enhancers can act (26). To identify putative target genes for MAE-1 and MAE-126, we determined the TAD structure at the MAE-1/MAE-126 locus using Hi-C (47). We did not observe differences in the TAD structure between 786-O and M1A cells in the immediate MAE-1/MAE-126 region, indicating that perturbation of the higher order chromatin structure was unlikely to cause MAE-1 or MAE-126 activation in the metastatic cells (Fig. 4A). However, other genomic loci, such as that immediately distal to the MAE-1/MAE-126, did show evidence for altered TAD structure in the metastatic cells (Fig. 4A; Supplementary Fig. S7A). Although in this particular case the TAD change was not associated with altered gene expression within this locus, it suggests that TAD alterations could be relevant for metastasis gene activation in other loci. Both MAE-1 and MAE-126 resided in the same 600-kb TAD that also included two protein-coding genes, DARS and CXCR4, of which CXCR4 was located closer to the MAE-1/MAE-126 cluster (Fig. 4B). Only CXCR4 expression was reduced in the CRISPRi cells targeting MAE-1 and MAE-126 (Fig. 4C and D), supporting the observation that the most proximal genes were the most likely MAE targets (Fig. 1F and H). Loss of CXCR4 expression was also evident at the protein level (Fig. 4E). The DARS locus, on the other hand, harbored strong H3K27ac signal already in the parental cells, possibly explaining the insensitivity of DARS expression to MAE-1/MAE-126 perturbation (Supplementary Fig. S7B). In a genome-wide analysis, CXCR4 was the most significantly downregulated gene in the MAE-1/MAE-126 CRISPRi cells (Supplementary Fig. S7C; Supplementary Table S5). We also detected increased interaction between MAE-1, MAE-126, and the CXCR4 promoter in the M1A cells when compared with the nonmetastatic parental cells using 3C (Supplementary Fig. S7D).

To evaluate the specificity of the CXCR4-regulatory landscape in ccRCC, we used H3K27ac ChIP-seq to identify putative enhancers in the CXCR4 locus in an experimental system of triple-negative breast cancer metastasis derived from MDA-MB-231 cells (17, 48). As previously reported, the bone metastatic 1833-BoM cells in this series expressed high levels of CXCR4 mRNA (Supplementary Fig. S8A). We found no evidence of MAE-1 activity in these cells; however, there was a prominent H3K27ac peak, MAE-Br, upstream of CXCR4 in a region that was not active in ccRCC cells (Supplementary Fig. S8B). Comparative CRISPRi-mediated targeting of MAE-Br, MAE-1, and MAE-126 in 1833-BoM and M1A cells, respectively, demonstrated that MAE-Br was active only in the breast cancer cells, whereas MAE-1 and MAE-126 supported CXCR4 expression only in M1A cells (Supplementary Fig. S8C–S8F). These experiments confirmed the specificity of CRISPRi-mediated enhancer targeting and suggested that metastatic clones arising from different cancer types may use alternative enhancer configurations to activate shared metastasis genes.

Previously, the Polycomb Repressive Complex 2 (PRC2) has been shown to repress CXCR4 expression in nonmetastatic ccRCC clones (19). Subsequently, however, it has been reported that PRC2 does not play an instructive role in gene silencing, but rather accumulates at inactive promoters to strengthen a silenced state (49, 50). We therefore tested whether MAE-1 activation in nonmetastatic ccRCC could evict PRC2 from the CXCR4 promoter and induce gene expression. In line with the hypothesis, recruitment of the dCas9-VPR activator (51) to the MAE-1 locus was sufficient to decrease SUZ12 binding and histone H3 lysine 27 trimethylation (H3K27me3) at the CXCR4 promoter, with associated induction of CXCR4 expression (Fig. 4F–I). These results, together with previous observations (19), support a model whereby the PRC2 complex contributes to CXCR4 silencing and consequent suppression of metastatic potential when sufficient signals for gene activation are not present.

We then tested whether CXCR4 was the functional mediator of the prometastatic activity of the MAE-1/MAE-126 cluster in ccRCC. Reexpression of CXCR4 completely rescued the lung colonization phenotype in MAE-1–targeted M1A cells (Fig. 5A–D; Supplementary Fig. S9A and S9B). The MAE-1/MAE-126 cluster thus supports metastatic colonization by activating the expression of the chemokine receptor CXCR4, a gene linked to metastatic progression through clinical association analysis, as well as experimental approaches in multiple cancer types, including ccRCC (16–19).

The transcription factor HIF2A, a critical oncogenic driver in ccRCC (12), is necessary but not sufficient for CXCR4 expression in ccRCC (Supplementary Fig. S10A–S10E; ref. 19). We therefore tested whether HIF2A contributed to the activation of MAE-1/MAE-126. ChIP-seq revealed that HIF2A and its dimerization partner HIF1B strongly bound MAE-126, at a region that contained an HIF2A binding motif, in both nonmetastatic and metastatic cells, whereas no significant binding was observed elsewhere in the CXCR4 locus, including MAE-1 (Supplementary Fig. S11A–S11D). Moreover, mutating the HIF2A binding site decreased MAE-126 reporter activity (Fig. 6A), and CRISPR/Cas9-based mutational targeting of this site in the endogenous MAE-126 locus abolished CXCR4 mRNA and protein expression in M1A cells (Supplementary Fig. S11E and S11F). The contribution of HIF2A to CXCR4 expression was thus mediated via MAE-126, but HIF2A signaling alone was insufficient for MAE-1 activation.

We then explored the mechanisms by which MAE-1 activity is controlled. Transient reporter assays showed increased MAE-1 activation in metastatic cells when compared with the parental lines (Supplementary Fig. S12A and S12B), suggesting that the MAE-1 DNA sequence determined its activity at least partially independently of the chromatin context. Sequence analysis identified an evolutionarily conserved binding motif for the NF-κB transcriptional effector p65/RELA centered at the peak of the p300 bound MAE-1 region (Supplementary Fig. S12C), suggesting that the NF-κB pathway was responsible for MAE-1 activity. In line with this possibility, CRISPR/Cas9-based mutational targeting of this site abolished CXCR4 mRNA and protein expression (Fig. 6B–D). Mutating the putative p65/REL binding site decreased MAE-1 reporter activity (Fig. 6E). We observed a close correlation between NF-κB activity and CXCR4 mRNA levels as well as poor patient outcome in the large TCGA ccRCC data set (Fig. 6F; Supplementary Fig. S12D). The metastatic cells showed higher level of NF-κB activity at the transcriptome level when compared with the nonmetastatic cells (Supplementary Fig. S12E), and ChIP-seq analyses demonstrated p65 binding at the MAE-1 locus in M1A cells (Fig. 6G; Supplementary Fig. S12F). Importantly, mutating the MAE-1 locus decreased the metastatic colonization potential of M1A cells (Fig. 6H–J; Supplementary Fig. S12G and S12H), providing direct genetic evidence for the role of MAE-1 in metastatic colonization. In agreement, metastatic ccRCC cells displayed significantly higher activation of NF-κB (Supplementary Fig. S13A–S13D), and inhibition of NF-κB pathway activity through expression of a super repressor mutant (Supplementary Fig. S13E–S13G) resulted in downregulation of CXCR4 in the M1A cells (Supplementary Fig. S13H). Conversely, expression of a constitutively active IKK-2 mutant, a kinase that can activate the NF-κB pathway but also other pathways, resulted in increased NF-κB activity, MAE-1 activation, and an increase in CXCR4 expression in 786-O cells (Supplementary Fig. S14A–S14F). Taken together, the expression of the ccRCC metastasis gene CXCR4 is controlled by the combined action of two independent enhancers, MAE-1 and MAE-126. MAE-126 is activated by the canonical VHL–HIF2A pathway already at the premetastatic stage, and MAE-1 is subsequently induced and regulated by the NF-κB pathway, with possible contribution from other MAE-1/MAE-126–binding factors.

The highly conserved p65/RELA motif in MAE-1 and its sensitivity to mutation suggest that MAE-1, and possibly other MAEs, have a role in normal gene regulation and are not simply de novo cancer enhancers. Consistent with this idea, H3K27ac ChIP-seq data from multiple tissue types in the mouse ENCODE data set (52) revealed that the majority of MAEs belong to tissue-restricted enhancer modules representing multiple developmentally distinct cell linages (Supplementary Fig. S15A). The orthologous mouse MAE-1 locus is active in the spleen and thymus, whereas no general activation is seen in embryonic, postnatal, or adult mouse kidney (Fig. 7A). Furthermore, interrogation of the relationship between ccRCC MAEs and 2,328,936 normal human enhancers identified by the Roadmap Epigenomics project (53) revealed that 94% of MAEs colocalized with normal enhancers belonging to several different physiologic regulatory modules, such as the prominent lymphoid regulatory module that harbors MAE-1 (Fig. 7B; Supplementary Fig. S15B). Similar analysis of the enhancers activated in the metastatic breast cancer clones 1833-BoM and 4175-LM2 gave concordant results, with the CXCR4-associated MAE-Br being part of the lymphoid enhancer cluster (Supplementary Fig. S16A–S16D). Pathway analysis of the MAE-associated genes also identified several functional categories involved in immune cell biology (Supplementary Fig. S3). Finally, p65 was found to bind MAE-1 in several human lymphoblastoid cell lines which are derived from the B-cell lineage, further supporting the role of the NF-κB pathway in the activation of MAE-1 (Fig. 7C; ref. 54). Thus, MAE-1 exemplifies a functionally conserved and tissue-restricted lymphoid enhancer co-opted by aggressive ccRCC cells during metastatic progression.

Here, we show that metastatic cancer clones activate a distinct set of enhancers and that these drive the expression of metastasis-associated genes. We provide functional evidence for the importance of specific distal enhancers as mediators of metastatic colonization and demonstrate how the activity of these enhancers is supported in metastatic cells. Furthermore, we find that the MAE landscape in ccRCC consists of multiple superimposed enhancer modules aberrantly co-opted from diverse physiologic cell lineages, providing clues to the mechanisms of metastasis gene activation more generally (Fig. 7D).

Despite significant efforts, metastasis driver mutations have thus far not been identified (4, 6, 7). Yet, as demonstrated by extensive transcriptional profiling of human cancer data sets and experimental systems, metastatic progression is associated with recurrent gene expression patterns (55). How these patterns emerge has remained largely unknown. Our results support a model whereby metastatic transcriptional programs arise from de novo combinatorial activation of multiple distinct and developmentally distant transcriptional modules. These results are in line with recent reports that have demonstrated extensive chromatin remodeling in association with metastatic progression in other tumor types as well (56–58). The complex patterns of MAEs could explain why individual mutational events do not seem to be sufficient for the establishment of a metastatic state. Also, our results and those by Morrow and colleagues (57) suggest that metastatic gene expression programs in human systems are driven by multiple transcriptional effectors, not individual factors as seen in mice (56, 58). This discrepancy may reflect the differences in the kinetics of metastatic progression in mice and humans, respectively.

A specific example of cross-lineage enhancer co-option as a promoter of metastasis is provided by our discovery of the prometastatic enhancer cluster that is supported by the ccRCC core driver HIF2A and an NF-κB–dependent lineage-restricted lymphoid enhancer. Both of the two main constituent enhancers of this cluster, MAE-1 and MAE-126, are needed for its metastasis-promoting function. Given the tissue-specific nature of VHL loss–driven tumorigenesis, it is likely that other factors in addition to HIF2A and NF-κB also contribute to MAE-1/MAE-126 function. The identification of NF-κB as a regulator of a prometastatic enhancer is in line with previous observations that have linked high NF-κB activity to cancer progression and poor patient survival in ccRCC (23, 25). However, the level of NF-κB activity that directly follows VHL inactivation (20, 21) is insufficient for the activation of the prometastatic program. Rather, our model suggests that further enhancement of NF-κB is needed for the activation of MAE-1 and possibly other prometastatic target loci. The detailed mechanisms of such NF-κB hyperactivation remain to be investigated.

Renal cancers can invade the venous vasculature as tumor thrombi. Yet, even in cases in which vascular involvement has been macroscopically confirmed, a significant fraction never forms distant metastasis (59, 60). Thus, in line with evidence from other experimental and clinical cancer contexts (3), early invasion and intravasation do not seem to be the major rate-limiting steps of renal cancer progression. We show that the MAE-1/MAE-126 cluster mediates metastatic colonization by regulating the expression of CXCR4, a gene encoding a G-protein coupled chemokine receptor that supports metastatic progression in several different cancer types by increased survival and chemotactic migration (16–19). Whether the MAE-1/MAE-126 cluster also contributes to earlier steps of metastasis remains unclear, however. In line with previous data (19), we find that PRC2, a critical multiprotein complex involved in gene silencing (61), is linked to the CXCR4 repressed nonmetastatic cell state. That MAE-1 activation is sufficient for PRC2 eviction from CXCR4 promoter suggests that PRC2 contributes to the maintenance of CXCR4 repression when sufficient activating signals are not present, but it does not determine CXCR4 silencing. This model is in line with data from recent mechanistic analyses of PRC2 function (49, 50), and it suggests that reduced PRC2 function could increase transcriptional variability, which in turn could favor tumor progression (19, 62). Additionally, our comparative functional dissection of the mechanisms underlying CXCR4 expression suggests that alternative regulatory modules may drive its expression in different cancer types. Thus, depending on the tissue of origin, different enhancer configurations may converge on shared metastatic mediators. Analysis of the regulatory basis of metastasis and the mechanisms that maintain transcriptional stability in nonmetastatic clones is therefore needed for a comprehensive understanding of the origins of metastatic programming.

In genetically engineered mouse models, the transition from normal homeostasis to neoplastic growth is characterized by aberrantly activated tissue regeneration programs that in conjunction with perturbed cell identity networks support the establishment of oncogenic enhancers (63). Our findings of cross-lineage enhancer co-option at the metastatic stage complement this model and suggest the possibility that tumor development from benign to metastatic is driven by progressive lineage infidelity that facilitates the activation of increasingly aggressive cancer phenotypes. The biological implications of such broad perturbations in cellular identity and possible resultant conflicts in coexisting cell fate programs (64) remain to be explored.

The recent findings that modulation of enhancer activity can alter cancer risk (31–33, 38), that cancer-specific enhancer elements activate oncogenes (30), and that cancer-associated enhancer clusters are susceptible to pharmacologic inhibition (65, 66) have raised interest in understanding the mechanisms by which gene-regulatory elements contribute to cancer progression and metastasis (56–58). Our work demonstrates that in renal cancer, metastatic transcriptional programs are supported by a combination of multiple superimposed physiologic enhancer modules co-opted by metastatic cancer clones. Such cross-lineage enhancer co-option as a source of metastatic traits may be applicable to other cancer types as well. The complex regulatory origins of metastatic traits suggested by this model could have implications for therapeutic strategies that aim to reduce metastasis in patients.

The renal cancer cell lines 786-O TGL (786-O), OS-RC-2 TGL (OS-RC-2), and their metastatic derivatives 786-M1A (M1A), 786-M2B (M2B), and OS-LM1 (LM1), respectively, have been previously described (19), and they were obtained from J. Massagué (Memorial Sloan Kettering Cancer Center, New York) in 2014. The breast cancer cell line MDA-MB-231 and its bone and lung metastatic derivatives 1833-BoM and 4175-LM2, respectively, have been described (17, 48) and were obtained from J. Massagué in 2014. The identity of the renal cancer cell lines has been confirmed by Sanger sequencing–based detection of known unique homozygous VHL mutations in April 2017. The breast cancer cell lines have been exome sequenced as previously reported and they carry the expected mutational complement (42). All cancer cell lines have been confirmed negative for Mycoplasma by biannual tests using the MycoAlert Mycoplasma Detection Kit (Lonza, LT07-318). At the time of the study, none of the cell lines were listed in the ICLAC database of commonly misidentified cell lines. For lentivirus production, we used HEK293T cells. Renal cancer cells were cultured in either RPMI-1640 supplemented with 10% FBS, penicillin (100 U/mL), and streptomycin (μg/mL) or DMEM/F12 supplemented with B27 (Invitrogen), streptomycin (μg/ml), and EGF/FGF (Peprotech 10 ng/mL). Other cells were grown in DMEM, supplemented with 10% FBS, penicillin (100 U/mL), and streptomycin (μg/mL). Cell lines were used in experiments within 10 passages from thawing.

For ChIP-seq, subconfluent cells were cross-linked for 10 minutes at room temperature in 1% formaldehyde in growth media, followed by 5 minutes of quenching with 0.125 mol/L glycine. Cells were washed twice with PBS, the supernatant was aspirated, and the cell pellet was frozen in liquid nitrogen and stored at −80°C. For tumor ChIP, samples were homogenized in PBS using a Precellys instrument (Bertin) before cross-linking. Protein A/G magnetic beads (100 μL; Thermo Scientific, 26162) were blocked with 0.5% BSA in PBS, followed by incubation with antibody at 4°C for a minimum of 4 hours. The following antibodies were used: H3K27ac (Abcam, ab4729), H3K4me1 (Abcam, ab8895), p300 (Santa Cruz Biotechnology, sc-585), SUZ12 (Abcam, ab12073), H3K27me3 (Millipore, 17-622), HIF1B (Novus, NB100-110), HIF2A (Abcam, ab199), and p65 (Abcam, ab16502). Control experiments were done using rabbit polyclonal IgG (Abcam, ab27478). Cross-linked cells were resuspended and sonicated in lysis buffer (20 mmol/L Tris–HCl pH 8.0, 150 mmol/L NaCl, 2 mmol/L EDTA pH 8.0, 0.1% SDS, and 1% Triton X-100). Sonication was performed in a Bioruptor (Diagenode) for 14 cycles (30″ on/30″ off) at max output to obtain fragments of 100 to 500 bp. Sonicated lysates were cleared and incubated overnight at 4°C with antibody-bound magnetic beads. Beads were sequentially washed three times with low-salt buffer (50 mmol/L HEPES pH7.5, 140 mmol/L NaCl, 1% Triton) and once with high-salt buffer (50 mmol/L HEPES pH7.5, 500 mmol/L NaCl, 1% Triton). DNA was eluted in elution buffer (50 mmol/L NaHCO₃, 1% SDS) and cross-links were reversed for 3 hours (65°C, 1,000 rpm shaking). DNA was purified using the QuickClean II PCR Extraction Kit (Genescript L00419-100) according to the manufacturer's recommendations and eluted with 100 μL H₂O. ChIP-qPCR primers are listed in Supplementary Table S6.

Purified ChIP DNA was used to prepare Illumina multiplexed sequencing libraries with the KAPA Hyper Prep Kit (KR0961) Illumina platforms sample preparation protocol (v1.14). After adapter ligation, libraries were size-selected to 150 to 350 bp using Agencourt AMPure XP beads (Beckman Coulter, A63880) according to the protocol recommendations. Size-selected libraries were amplified for 15 cycles using the KAPA HiFi HotStart ReadyMix. PCR libraries were then pooled in equimolar concentrations and sequenced.

RNA was extracted in triplicate from subconfluent cells using the RNeasy Mini Kit (Qiagen 74104) following the manufacturer's recommendations. RNA quality was assessed with the RNA Nano 6000 Kit (Agilent 5067-1511) and Agilent Bionalyzer 2100 instrument. Purified RNA (1 μg) was used for the preparation of Illumina multiplexed sequencing libraries following the KAPA Stranded mRNA-seq Kit (KR0960) Illumina platforms sample preparation protocol (v3.15). After adapter ligation, libraries were cleaned using Agencourt AMPure XP beads (Beckman Coulter, A63880) according to the manufacturer's recommendations. Purified adapter-ligated DNA was amplified for 15 cycles. PCR libraries were then pooled in equimolar concentrations and sequenced.

All ChIP-seq and RNA-seq was performed on Illumina HiSeq 2500 or HiSeq 4000 systems using 50 bp single-end runs. For ENCODE data, fastq files were downloaded and processed locally. Raw ChIP-seq sequencing reads were aligned to hg19 or mm10 using bowtie2 (67), and the resulting sam files were converted into sorted bam files using samtools (68). Genome-wide correlation between replicates was calculated using 1,000 bp windows; regions of no signal were excluded. For cell lines, H3K27ac peaks were called using MACS2 (69), with options –broad and -p 1e-8 and the corresponding input samples as controls. Overlap between peaks was identified using the intersect function in bedops (70). For the identification of H3K27ac regions specifically enriched/depleted in the metastatic ccRCC populations, replicate ChIP-seq experiments for each cell line were combined and peaks were called using MACS2 with the corresponding parental/metastatic sample as a control. Overlaying the 786 and OS systems resulted in a shared set of 416 enriched and 223 depleted H3K27ac regions in the metastatic clones, respectively. Promoter-associated H3K27ac peaks within this set were filtered out by removing regions closer than 1,000 bp to H3K4me3 peaks as determined previously (19), resulting in a final set of 314 enriched and 161 depleted MAEs. The genomic distribution of these regions relative to known transcription start sites was determined using Homer (71). MEME suite modules Tomtom (72) and FIMO (73) were used for motif identification within sequences of interest. For peak visualization, we used the R packages BSgenome.Hsapiens.UCSC.hg19, GenomicAlignments, and GenomicRanges (74). Briefly, reads were extended to the modal length of 250 bp and duplicates were removed, after which genome coverage was calculated for the region of interest. For final ranking of MAEs, we subtracted the signal from the corresponding input sample and used DESeq2 (75) to calculate a fold change and P value for H3K27ac enrichment. Overlap with Roadmap Epigenomics enhancers was determined using bedops. A minimum of 20% overlap with 314 enriched MAEs was needed for an element to be assigned to an MAE. The activity of an MAE in the 111 Roadmap epigenetic states was determined as the sum of active Roadmap enhancers overlapping with each of the MAEs. The correlation between cell lines and clinical samples was assessed by averaging cell line and clinical sample H3K27ac data across 38,625 H3K27ac-enriced regions identified in the ccRCC cell lines.

RNA-seq sequencing reads were aligned to hg19 using RSEM (76) and bowtie2 with default settings. For determining the association between the expression of transcripts that correlated with MAE activity and clinical outcome, we used ccRCC RNA-seq data from the TCGA project (11). Briefly, we calculated sums of z-scores (MAE activation score) using all genes that were first or second most proximal to an MAE and showed a corresponding average gene expression fold change of ±1.5 in the metastatic M1A and LM1 cells when compared with the corresponding parental populations. Pathway analysis was performed using Ingenuity Pathway Analysis, and B-H multiple testing corrected P value of <0.001 was considered statistically significant. The P value for clinical association was determined using a Cox proportional hazards model, where the MAE activation score was treated as a continuous variable. The NF-κB score was calculated as a sum of z-scores of a previously compiled set of ccRCC-associated NF-κB target genes (25).

The sequencing data that support the findings of this study have been deposited to the Gene Expression Omnibus under the accession code GSE98015. The ENCODE data were downloaded from https://www.encodeproject.org. The Roadmap Epigenomics data were downloaded from http://www.roadmapepigenomics.org.

Lentiviruses were produced by transfecting HEK293T with standard packaging vectors using Lipofectamine 2000 (Thermo Scientific, 11668027). Viral supernatants were harvested 48 hours following transfection and filtered through a 0.45-mm PVDF syringe filter. Virus supernatant was then used for infection in the presence of polybrene (Millipore, 8 μg/mL). For CRISPRi enhancer inactivation cells were first stably transduced with a lentiviral construct expressing KRAB-dCas9-mCherry (Addgene 60954; ref. 77). For CRISPR activation (CRISPRa)–mediated enhancer activation, cells were transduced with a lentiviral construct expressing dCas9-VPR-mCherry (A. Obenauf, Research Institute of Molecular Pathology, Vienna). Single-guide RNA (sgRNA) pairs with separate U6 promoters targeting the p300 enhancer peak were cloned in tandem into the pKLV-U6gRNA plasmid (Addgene 50946; ref. 78) and lentivirally transduced into cells expressing either KRAB-dCas9-mCherry or dCas9-VPR-mCherry, for CRISPRi or CRISPRa, respectively. Double-positive cells for the mCherry and BFP markers were selected by fluorescence-activated cell sorting (FACS). Control cells expressed a nontargeting sgRNA pair. CRISPRi-mediated reduction in H3K27ac and p300 was tested by ChIP-qPCR. Genome-wide assessment of CRISPRi was performed by ChIP-seq. Briefly, enrichment in a combined set of H3K27ac regions was determined and the difference between the CRISPRi and control sample was assessed using DESeq single sample mode (79). For CRISPR/Cas9 mutagenesis, cells were stably transduced with a modified version of the Cas9-expressing plasmid Lenti-Cas9 Blast (Addgene 52962; ref. 80) in which the blasticidin resistance gene was replaced with the fluorescent marker mCherry. The sgRNAs targeting the HIF binding site of MAE-126, the evolutionary conserved region in MAE-1 or HIF2A locus, were cloned into pKLV-U6gRNA and stably transduced into the cells. Double-positive cells for mCherry and BFP were then FACS selected. sgRNAs were designed as previously described (81). A full list of the sgRNAs used in this study is available in Supplementary Table S6.

Anonymized ccRCC RNA-seq data and associated clinical annotations were downloaded from the TCGA data portal (https://tcga-data.nci.nih.gov/). Four tumor samples for ChIP-seq analyses were collected from nephrectomy specimens. The study was conducted in accordance with the Declaration of Helsinki. All patients were recruited for tissue donation by providing written consent and after study approval from an Institution Review Board (NRES Committee East of England–Cambridge Central; REC 03/018). Raw ccRCC and corresponding normal kidney RNA-seq and ChIP-seq (H3K27ac and inputs) data were downloaded from the National Center for Biotechnology Information (NCBI) Short Read Archive (Bioproject ID PRJNA396588 and PRJNA340239; ref. 39). Illumina adapters and poor-quality sequence were removed, and the data were processed as described above for cell lines. H3K27ac signal in H3K27ac-enriched MAEs was normalized to the corresponding input control for all samples, after which each tumor was compared with the corresponding normal sample. A fold change of >3 was used as a cutoff for calling an MAE active. Gene expression fold change over the corresponding normal tissue was calculated for genes in the proximity of all MAEs in all samples, after which the “MAE on” and “MAE off” groups were compared.

All animal experiments were performed in accordance with protocols approved by the Home Office (UK) and the University of Cambridge ethics committee (PPL 70/7990). For experimental lung metastasis assays, we inoculated 3 × 10⁵ tumor cells, resuspended in 100 μL of PBS, into the lateral tail vein of 5- to 7-week-old male NOD/SCID or NSG mice, obtained from Charles River Laboratories. Tumor growth was monitored by IVIS bioluminescence imaging (PerkinElmer). At the experimental endpoint, lungs were harvested for immunohistochemistry. For the orthotropic tumor assay, a 1-cm incision was made above the left kidney of 6-week-old male NOD/SCID mice (Charles River Laboratories). The kidney was exposed, and 1 × 10⁶ tumor cells, resuspended in 25 μL of Matrigel (BD), were injected under the renal capsule. Tumor growth was monitored by IVIS bioluminescence imaging (PerkinElmer). At the experimental endpoint, lungs were harvested for IVIS imaging. For subcutaneous tumor growth assays, 5 × 10⁵ cells, resuspended in Matrigel (BD), were injected in the flanks of male NSG mice (Charles River Laboratories). Tumor growth was monitored by caliper measurement. Tumor volume was calculated using the equation V = (L × W²) × 0.5, where L is the length and W is the width of the tumor. After 6 weeks, animals were euthanized and tumors were harvested.

Lungs and subcutaneous tumors were collected and fixed overnight in 4% paraformaldehyde, washed, embedded in paraffin, and sectioned. Hematoxylin and eosin staining was performed by standard methods. Human Vimentin (Cell Signaling Technology; cat. 5741, 1:100) staining was performed in a Bond-Max instrument (Leica) using Bond Polymer Refine Detection reagents (Leica) according to the manufacturer's protocol (IHC Protocol F).

In-nucleus Hi-C was performed on 20 million cells per cell line as previously described (82), with only minor modifications: DNA extraction and all purifications were performed using Agencourt Ampure XP beads (Beckman Coulter), and ligations were performed in a volume of 1 mL. Final libraries were amplified using 6 cycles of PCR and 2 samples pooled per sequencing lane. The 150 bp paired-end sequencing was performed on the Illumina HiSeq 4000 instrument.

For Hi-C data analysis, reads were aligned, filtered, and summarized into frequency matrices with HiC-Pro (83), using the default settings and by setting the minimum and maximum fragment sizes to 100 and 600, respectively. The interaction matrices for each experiment were normalized with ICE (Iterative Correction and Eigenvector Decomposition) using the implementation in the iced python package included in HiC-Pro, with 100 iterations. TADs were calculated with TADbit (84) from interaction matrices based on 40-kb genomic bins. HiCPlotter (85) was used for creating figures incorporating the normalized matrices, genes, and enhancer information tracks, as well as the TAD distribution for each experiment. For plotting purposes, the normalized matrices corresponding to replicates of the same experiment were summed together.

Samples for 3C were cross-linked, lysed, and HindIII digested similarly to Hi-C samples (82). After digestion, ligation was performed using T4 DNA Ligase (NEB) at 16°C overnight. DNA was purified using Agencourt Ampure XP beads (Beckman Coulter) and quantified (Quant-iT Pico Green, Life Technologies). qPCRs were performed, in triplicate, using JumpStart Taq ReadyMix for High-Throughput qPCR (Sigma-Aldrich) with 400 ng 3C library and 0.3 μmol/L each primer. Reactions (25 μL) were performed according to the manufacturer's protocol but using a Tm of 62°C and 50 cycles. Primer efficiencies were determined for each primer pair by standard curve generation on serial dilutions of a digested and ligated sample using the above conditions. The ratio of interactions for each sample was calculated using: Ratio metastatic/parental = ([Efficiency of test primers − Ct_{(test metastatic)} − Ct_{(test parental)}]/Efficiency of control primers − Ct_{(control metastatic)} − Ct_{(control parental)}]).

Cell line and tumor total RNA was extracted using RNAzol (Sigma) according to the manufacturer's protocol. RNA (1 μg) was used to generate cDNA with the High-Capacity cDNA Reverse Transcription Kit (Thermo Scientific). For qPCR, we used the StepOnePlus instrument (Thermo Scientific) with predesigned TaqMan gene expression assays (Thermo Scientific): CXCR4 (Hs00607978_s1), DARS (Hs00154683_m1), IL6 (Hs00174131_m1), IL8 (Hs00174103_m1), and TBP (Hs00427620_m1). TBP served as housekeeping control, and data were analyzed using the double delta Ct method. For immunoblotting, subconfluent cells were pelleted, washed with cold PBS, and lysed with RIPA buffer. Nuclear and cytoplasmic extracts were obtained using the NE-PER Nuclear and Cytoplasmic extraction reagents (Thermo Scientific) according to the recommendations of the manufacturer. Proteins were separated by SDS-PAGE, transferred to PVDF membranes, and blotted with antibodies against CXCR4 (Abcam, ab124824, 1:500), HIF2α (Novus, NB100-122, 1:500), VHL (Cell Signaling Technology, 2738, 1:500), phospho-p65 (Cell Signaling Technology, 3033, 1:1,000), p65 (Santa Cruz Biotechnology, sc8008, 1:1,000), TBP (Abcam, ab63766, 1:5000), B-tubulin (Abcam, ab6046, 1:10,000), and β-actin (Sigma, A5441, 1:5,000), β-tubulin (Sigma, T5201, 1:5,000). Secondary antibodies were HRP-conjugated (Dako, 1:5,000).

We cloned the CXCR4 cDNA into a modified pLVX-Puro plasmid (Clonetech) where the puromycin resistance gene was replaced with a hygromycin resistance gene. The empty vector was used as a control. Lentivirus were generated in HEK293T cells as described above. Hygromycin selection started the following day (Invivogen, 900 μg/mL). CXCR4 expression was determined by qPCR. For NF-kB inactivation in metastatic ccRCC cells, we cloned the IKBalpha-mut super repressor cDNA (Addgene 15291; ref. 86) into the pLVX-Puro vector (Clonetech). Similarly, for NF-kB activation in parental ccRCC cells, the IKK-2 S177E S181E mutant (Addgene 11105; ref. 87) was cloned into the pLVX-Puro vector. For VHL reintroduction in metastatic cells, HA-VHL (Addgene 19234; ref. 88) was cloned into the pLVX-Puro vector. Empty vector was used as a control. Lentiviruses were generated in HEK293T cells as described above. Puromycin selection started the following day (Invivogen, 4 μg/mL). NF-kB activity was determined by transient reporter assay experiments and expression of downstream target genes (IL6 and IL8).

Enhancer activity was determined by luciferase reporter assays. The region comprising the p300 peak of MAE-126 and MAE-1 (approximately 900 bp) was cloned into the pNL1.1[Nluc] vector (Promega, N1001) containing a NanoLuc Luciferase downstream of a minimal promoter using the following primers: 5′-ACGTGCTAGCGTTGGAGGGGAGATCTAGGC-3′ and 5′-ACGTCTCGAGGCCTGTCCATCCACTACCTA-3′ for MAE-126, and 5′-ACGTGCTAGCGTGATGACCTTTTCCACTTGCT-3′ and 5′-ACGTCTCGAGGTCCAGATTCCGACCTTCC-3′ for MAE-1. Genomic DNA from 786-O cells was used as template. Parental and metastatic ccRCC cells were cotransfected with 2 μg of enhancer-NanoLuc reporter plasmid and 200 ng of pRL-TK Renilla luciferase control vector (Promega E2241) as a normalization control using an Amaxa Nucleofector (program T-20, reagent V, Lonza). After 48 hours of incubation NanoLuc luciferase activity was measured using the NANO-Glo Luciferase assay (Promega, N1110) and Renilla luciferase activity was measured with the Renilla luciferase assay system (E2810). Mutation of the HIF2A binding site in pNL1.1[Nluc]-MAE-126 was generated using the QuickChange II site-directed mutagenesis kit (Agilent) with the following oligos: 5′-CTCAGCCCTGGCATTTCATCTTTTAGGCGCCGTTCTTAAAGGG-3′ and 5′-CCCTTTAAGAACGGCGCCTAAAAGAT GAAATGCCAGGGCTGAG-3′. Mutation of the p65 binding site in pNL1.1[Nluc]-MAE-1 was generated using the QuickChange II site-directed mutagenesis kit (Agilent) with the following oligos: 5′-GCAGAGCCTTGGGGATTTAAAGGGTGGGGTGATGAATC-3′ and 5′-GATTCATCA CCCCACCCTTTAAATCCCCAAGGCTCTGC-3′. To determine NF-kB activity, five tandem NF-kB response elements were cloned into pNL1.1[Nluc] vector using the following oligos: 5′-CTAGCTTGGGGACTTTCCGCTGGGGACTTTCCGCTGGGGACTTTCCGCTGGGGACTTTCCGCTGGGGACTTTCCGCAC-3′ and 5′-TCGAGTGCGGAA AGTCCCCAGCGGAAAGTCCCCAGCGGAAAGTCCCCAGCGGAAAGTCCCCAGCGGAAAGTCCCCAAG-3′. Cells were cotransfected with 50 ng of pNL1.1[Nluc]-NF-kB(5×) or respective empty-vector control and 10 ng of the Renilla luciferase plasmid (pRL-TK). NanoLuc and Renilla luciferase activities were measured as described above.

For in vitro proliferation assays, 1 × 10⁴ cells were seeded in triplicate in 12-well plates (one plate per day up to 5 days). Each day, cells were resuspended and counted using a Vi-Cell XR instrument (Beckman Culture).

Statistical analyses were conducted in R and GraphPad Prism. Sample sizes are denoted in figure legends. P values lower than 0.05 were considered statistically significant. No statistical method was used to predetermine sample size. For animal bioluminescence data, one-tailed Wilcoxon rank-sum test was used and for experiments with ≥2 tests a Bonferroni correction was applied. No animals were excluded from the analyses. For histologic tumor count analyses, the experimental groups were blinded; for other experiments, the experimental groups were not randomized or blinded. For parametric tests, data normality was assessed using the Shapiro–Wilk test, variance between the groups was not assumed to be equal. Log transformation was applied as indicated in the figure legends. For tumor count data, a one-tailed Student t test with the Welch approximation of degrees of freedom was used. Bonferroni correction was applied for experiments with ≥2 tests. For reporter assays, a two-tailed Student t test with the Welch approximation of degrees of freedom was used. For testing the association between MAE activation and mRNA expression of the nearby genes, a one-sample two-sided Wilcoxon rank-sum test was used. For survival analysis, P value was derived using a Cox proportional hazards model in which sums of Z-scores were used as a continuous variable. For visualization of survival data, Kaplan–Meier curves were used. For qRT-PCR, three independent experiments are shown, each of which is the average of three technical replicates. Box plots represent median and 25th and 75th percentiles, the whiskers extend to data extremes or 1.5× interquartile range as indicated in the figure legends.

A.Y. Warren participated in a Roche Scientific advisory board on the potential for PD-L1 testing in RCC (one meeting only). G.D. Stewart reports receiving a commercial research grant from Pfizer, has received honoraria from the speakers bureau of Pfizer, and is a consultant/advisory board member for Pfizer and EUSA Pharma. No potential conflicts of interest were disclosed by the other authors.

Conception and design: P. Rodrigues, S. Vanharanta

Development of methodology: P. Rodrigues, K. Saeb-Parsy

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): P. Rodrigues, S.A. Patel, L. Harewood, E. Vojtasova, S.E. Syafruddin, M.N. Zaini, E.K. Richardson, J. Burge, A.Y. Warren, G.D. Stewart, K. Saeb-Parsy, S. Vanharanta

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): P. Rodrigues, I. Olan, S.A. Samarajiwa, S. Vanharanta

Writing, review, and/or revision of the manuscript: P. Rodrigues, A.Y. Warren, G.D. Stewart, K. Saeb-Parsy, S. Vanharanta

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): P. Rodrigues, G.D. Stewart

Study supervision: S. Vanharanta

Other (sample acquisition): A.Y. Warren

We thank A. Speed for technical assistance, P. Jones for discussion, A. Obenauf for discussion and the CRISPR activator construct, and V. Gnanapragasam for administration of the DIAMOND study, which allowed access to clinical RCC samples. Infrastructure for the Cambridge Urological Bio-repository was funded by the Cambridge Biomedical Research Campus and CRUK Cambridge Centre. The Human Research Tissue Bank is supported by the NIHR Cambridge Biomedical Research Centre. This work was supported by the Medical Research Council (MC_UU_12022/7).