Integrated epigenomic profiling reveals endogenous retrovirus reactivation in renal cell carcinoma

Abstract

Background: Transcriptional dysregulation drives cancer formation but the underlying mechanisms are still poorly understood. Renal cell carcinoma (RCC) is the most common malignant kidney tumor which canonically activates the hypoxia-inducible transcription factor (HIF) pathway. Despite intensive study, novel therapeutic strategies to target RCC have been difficult to develop. Since the RCC epigenome is relatively understudied, we sought to elucidate key mechanisms underpinning the tumor phenotype and its clinical behavior.

Methods: We performed genome-wide chromatin accessibility (DNase-seq) and transcriptome profiling (RNA-seq) on paired tumor/normal samples from 3 patients undergoing nephrectomy for removal of RCC. We incorporated publicly available data on HIF binding (ChIP-seq) in a RCC cell line. We performed integrated analyses of these high-resolution, genome-scale datasets together with larger transcriptomic data available through The Cancer Genome Atlas (TCGA).

Findings: Though HIF transcription factors play a cardinal role in RCC oncogenesis, we found that numerous transcription factors with a RCC-selective expression pattern also demonstrated evidence of HIF binding near their gene body. Examination of chromatin accessibility profiles revealed that some of these transcription factors influenced the tumor's regulatory landscape, notably the stem cell transcription factor POU5F1 (OCT4). Elevated POU5F1 transcript levels were correlated with advanced tumor stage and poorer overall survival in RCC patients. Unexpectedly, we discovered a HIF-pathway-responsive promoter embedded within a endogenous retroviral long terminal repeat (LTR) element at the transcriptional start site of the PSOR1C3 long non-coding RNA gene upstream of POU5F1. RNA transcripts are induced from this promoter and read through PSOR1C3 into POU5F1 producing a novel POU5F1 transcript isoform. Rather than being unique to the POU5F1 locus, we found that HIF binds to several other transcriptionally active LTR elements genome-wide correlating with broad gene expression changes in RCC.

Interpretation: Integrated transcriptomic and epigenomic analysis of matched tumor and normal tissues from even a small number of primary patient samples revealed remarkably convergent shared regulatory landscapes. Several transcription factors appear to act downstream of HIF including the potent stem cell transcription factor POU5F1. Dysregulated expression of POU5F1 is part of a larger pattern of gene expression changes in RCC that may be induced by HIF-dependent reactivation of dormant promoters embedded within endogenous retroviral LTRs.

Keywords: Cancer epigenetics; Cancer stem cell; Kidney cancer; Regulatory genomics; Renal cell carcinoma; Transcription factors.

Fig. 1

Overview of patient samples and data sets used for integrated analyses. (a) Primary culture of tumor and matched-normal tubule cells from three patients. Renal cell carcinoma tumor nephrectomies from three patients were used to derive primary cultures of proximal tubules and renal cell carcinoma. Microscope images scale bar = 100 μm. (b) Cytogenetic analysis of primary tumor cultures. Karyotype analysis of the carcinoma cultures revealed loss of the short of arm of chromosome 3 in all three patient samples. Sanger sequencing of the VHL gene in these same samples identified inactivating mutations in the remaining copy. (c) Example of integrated analysis at the STC2 gene locus. DNase-seq and RNA-seq datasets were also generated from the primary tubule and carcinoma cultures and compared to HIF ChIP-seq datasets from the 786-O renal cell carcinoma cell line and RNA-seq expression data from TCGA. STC2, a canonical HIF target gene, exhibited several differential DHS (blue arrows), some of which coincided with HIF binding determined by ChIP-seq. Compared to normal tubules, the STC2 transcript was strongly induced in the primary tumor cultures and in the TCGA tumor samples (depicted with 10% outlier trim for clarity).

Fig. 2

Shared regulatory landscapes in tubules and matched renal cell carcinomas from three patients. (a) Comparison of the shared regulatory landscape among patient samples. The three tubule samples shared a significant proportion of DHSs. Each tumor's landscape of DHSs incorporated a different fraction of DHSs from its tubule of origin and activates de novo DHSs. In the tumors, most of the tubule-derived DHSs were shared with tubule-derived DHSs from other patients. In contrast, a smaller fraction of RCC-derived de novo DHSs was shared among patient tumors. (b) Comparison of DNase-seq data by principal component analysis. While the tubule cultures from all three patients (brown spheres, in replicate) were tightly clustered, each tumor (red spheres, in replicate) and the 786-O cell line (blue spheres, in replicate) occupied a unique position in regulatory space. (c) Differential DHSs show highly concordant patterns of accessibility across patient samples. In pairwise comparisons, the shared differential DHSs were classified as concordantly upregulated in the tumor samples (gold), downregulated in the tumor samples (blue) or discordant in the two patient samples being compared (grey). The majority (>95%) of shared differential DHS showed concordant up- or downregulation.

Fig. 3

Concordant tumor regulatory landscapes expose key transcription factor drivers of RCC. (a) HIF-binding only accounts for a small proportion of the differentially accessible RCC regulatory landscape. ChIP-seq datasets for HIF1A and HIF2A showed substantial overlap with each other and most of these peaks coincided with a DHS in the tubule and/or RCC DNase-seq datasets. Most HIF peaks in DHSs mapped to non-changing/constitutive DHS in the tubules and RCC. (b) Differentially accessible HIF-bound DHSs show different patterns of accessibility across patient samples. Of the 904 HIF peaks that mapped to differentially accessible DHSs in at least one patient sample, most did not show significant change across the other patients' samples. (c) Transcription factors with changing expression located near HIF binding sites. The expression levels of 213 transcription factors changed by >1.5× in at least one patient sample and exhibited at least one HIF-bound DHS within 250 kb of their transcription start site (TSS). Many of these contained numerous HIF binding sites in proximity to their TSS, including transcription factors linked to renal cell carcinoma (ZEB2, BHLHE41) and POU5F1. (d) Selective expression of transcription factors in cancer. The transcription factors that were expressed (FPKM>1) and changing by at least 1.5-fold in any of the three patient samples (from panel C) were examined for differential expression in a wide range of tumors that have matched normal tissues available in the TCGA RNA-seq expression dataset (209 transcription factors are depicted; 4 factors are not mapped in the TCGA RNA-seq data). Transcription factors with RCC-selective increased expression are highlighted (e.g. HSF4, BHLHE41, ZEB2, POU5F1, etc.).

Fig. 4

Correlation of DNA binding motif enrichments with gene expression identifies enrichment for POU5F1 in RCC. (a) Transcription factor enrichment. Examination of differentially accessible or non-changing DHSs revealed different classes of transcription factors whose DNA binding recognition sequences were enriched in each category. The motif families containing transcription factors with genetic evidence linked to renal cell carcinoma susceptibility (i.e. MYC, BHLHE41, ZEB2 and HIF) and the stem cell related transcription factor POU5F1 (OCT4) are indicated. (b) Examination of RNA expression identifies candidate POU-family transcription factors driving motif enrichments in the DHS landscape. Since multiple transcription factors within the POU family share redundant DNA binding motifs, examination of transcription factor expression patterns may identify specific family members that are driving motif enrichment signatures. Examination of the differential gene expression patterns of these family members in RCC vs. tubules in our primary cultures and in the TCGA RNA-seq dataset revealed upregulation of POU5F1 in RCC. (c) Expression of POU5F1 in diverse somatic tumors. The mRNA expression levels of the stem cell related transcription factor POU5F1 (OCT4) in several non-germ cell tumors was compared to their matched normal tissue controls. The ends of the bar plots represent the 25th and 75th quartiles with whiskers representing 1.5× inter-quartile range (10% outlier trim applied for clarity).

Fig. 5

A novel human-specific promoter initiates long RNA transcripts through the PSORS1C3-POU5F1 locus in RCC. Overview of the POU5F1-PSORS1C3 genomic locus (hg19 chr6:31,125,253-31,156,354). RNA-seq tracks for the primary patient samples and the RCC cell line 786-O revealed a novel transcript originating from a DHS ~16 kb upstream of the known human embryonic stem cell (hESC) TSS. This transcript reads through the PSORS1C3 lncRNA gene into the POU5F1 gene body (green shaded box). ChIP-seq in 786-O cells revealed binding of HIF components (HIF1α, HIF2α, HIF1β) to this DHS with evidence of histone modification typical of active transcription across the entire transcript (H3K36Me3). This DHS was also associated with histone modifications characteristic of an active promoter in human renal epithelial cells (HRE), i.e. positioned nucleosomes marked by H3K4Me3 and depletion of the repressive H3K27Me3 mark. Examination of sequence conservation showed that this alternate promoter lies within a complex tandem long terminal repeat (LTR) element that is unique to humans (blue shaded box). CRG, Center for Genomic Regulation.

Fig. 6

The novel transcript for POU5F1 exhibits HIF dependence and POU5F1 expression levels correlate with patient survival. (a) Schematic mapping of POU5F1 5′-RACE PCR products. 5′-RACE performed on 786-O RNA captured a transcription start site originating in the -16 kb DHS (blue shaded box) therefore defining a novel isoform of POU5F1 that is produced following read through of the PSORS1C3 gene. Reverse primers in known POU5F1 exons (e.g. EX1 = reverse primer in exon 1) were used to amplify the 5′-ends of the cDNA molecule captured by 5′-RACE and sequence mapped to the genome. The exon 2 primer (*) failed to yield mappable sequence. The exon 5 primer yielded 3 different products (EX5-1, EX5-2, EX5-3). The location of RT-PCR primers to detect the canonical and novel POU5F1 transcript variants are also indicated. The inferred isoform structure is compared to EST KY781167. (b) Canonical and novel POU5F1 transcripts exhibit HIF-dependence. RT-PCR primers were used to quantify the canonical and novel POU5F1 transcripts in 786-O cells stably transduced with VHL (786-O + VHL) or empty vector (786-O + EV) cultured in normoxia or hypoxia (2% O₂) for 24 h. Expression levels (relative quantification, RQ) were calculated using the β-actin housekeeping gene (ACTB). Compared to the empty vector control, reintroduction of VHL protein into 786-O cells suppressed expression of both POU5F1 transcripts. Exposing 786-O + VHL cells to hypoxia induced both POU5F1 transcripts. N.B. expression scale differences between the canonical and novel transcripts. Error bars indicate standard deviations of three experimental replicates. *p < 0.05, **p < 0.005. (c) Immunohistochemistry of POU5F1 protein in renal cell carcinoma samples. POU5F1 (OCT4) immunohistochemistry was performed on RCC samples from 20 patients (5 from each of ISUP grades 1–4) and showed patchy nuclear positivity (arrows) in a single random sample from 4 patients. No nuclear staining was seen in any of the matched normal renal parenchyma from the same patients. (d) Overall survival as a function of POU5F1 expression in TCGA. Patients with POU5F1 expression data from TCGA (KIRC) were evenly divided into two groups split at the median expression level (233 RSEM normalized) and Kaplan-Meier curves for overall survival were plotted using the UCSC Xena browser tool.

Fig. 7

Activation of cryptic LTR-derived promoters in RCC. (a) Enrichment of HIF-bound DHS in LTR families. The ERV1 and ERVK families of LTR show significant enrichment for HIF-bound DHSs. Of these, the LTR2B subfamily shows the greatest number of HIF-bound DHSs (n = 34). (b) TCGA expression of selected genes putatively induced by HIF-bound cryptic promoters in LTRs. The ends of the bar plots represent the 25th and 75th quartiles with whiskers representing 1.5× inter-quartile range (10% outlier trim applied for clarity). All tumor-normal comparisons are significant (p < 1 × 10⁻¹⁰) by one-tailed t-test. (c) HIF-bound DHSs in LTRs show strand-specific induction of RNA transcripts. Since LTRs are intrinsically directional, enumeration of RNA-seq reads up to 1 kb on either the same or opposite strand of the LTR identifies elements with HIF-dependent promoter-like activity (increased transcripts in RCC samples compared to tubules with the same directional orientation as the LTR). The heatmap represents the ratio of the RCC/tubule read counts for each patient on the indicated strand. (d) A HIF-bound LTR is transcriptionally active and is associated with increased expression of the UBE2D2 gene. Similar to the alternate POU5F1 promoter, some of the HIF-bound LTRs that show promoter-like activity drive the expression of novel transcripts and increase the expression of nearby genes. Shown is the expression of UBE2D2 transcripts, which increases 1.76× in Patient 1's RCC compared to its matched tubule control. (e) HIF-bound LTR-induced genes exhibit HIF-dependence. RT-PCR primers were used to quantify the indicated transcripts in 786-O cells stably transduced with VHL (786-O + VHL) or empty vector (786-O + EV) cultured in normoxia or hypoxia (2% O₂) for 24 h. Expression levels (relative quantification, RQ) were calculated using the β-actin housekeeping gene (ACTB). N.B. expression scale differences between the canonical and novel transcripts. Error bars indicate standard deviations of three experimental replicates. *p < 0.05, ***p < 0.005, N.S. not significant (two-tailed t-test).