VULCAN integrates ChIP-seq with patient-derived co-expression networks to identify GRHL2 as a key co-regulator of ERa at enhancers in breast cancer

Abstract

Background: VirtUaL ChIP-seq Analysis through Networks (VULCAN) infers regulatory interactions of transcription factors by overlaying networks generated from publicly available tumor expression data onto ChIP-seq data. We apply our method to dissect the regulation of estrogen receptor-alpha activation in breast cancer to identify potential co-regulators of the estrogen receptor's transcriptional response.

Results: VULCAN analysis of estrogen receptor activation in breast cancer highlights the key components of the estrogen receptor complex alongside a novel interaction with GRHL2. We demonstrate that GRHL2 is recruited to a subset of estrogen receptor binding sites and regulates transcriptional output, as evidenced by changes in estrogen receptor-associated eRNA expression and stronger estrogen receptor binding at active enhancers after GRHL2 knockdown.

Conclusions: Our findings provide new insight into the role of GRHL2 in regulating eRNA transcription as part of estrogen receptor signaling. These results demonstrate VULCAN, available from Bioconductor, as a powerful predictive tool.

Keywords: Breast cancer; ChIP-seq; Dynamics; ER; GRHL2; H3K27ac; Master regulator; Network analysis; P300; VULCAN.

Fig. 1

An overview of VULCAN. (1) ChIP-seq analysis from multiple conditions is undertaken to generate cistrome data at multiple time points (or conditions). Binding events are then compared using differential binding analysis to establish log-fold change values for individual binding events between each time point. (2) Network generation was undertaken with ARACNe-AP by inferring all pairwise TF-target co-expression from patient datasets (e.g., TCGA breast and METABRIC datasets). (3) All the targets of each specific TF in the network, i.e., the individual regulons, are tested against the established changes in ER binding through the msVIPER algorithm [15] to identify proteins that interact with the target transcriptional factor and final prediction is given for potential interacting cofactors

Fig. 2

Dynamic behavior during early activation of ER. ChIP-qPCR of the TFF1 gene (a) at three time points shows increased binding of ER at 45 min after MCF7 cells are stimulated by estradiol. The previously reported maximum is followed by a decrease in the TFF1 promoter occupancy at 90 min. p values are generated by one-tailed t test. The maximal point at 90 min was identified as an outlier (> median + 2 × IQR); however, the removal did not alter the significance of results. (b) Differential binding analysis of ChIP-seq data at three time points to monitor the activation of ER. The ER exhibits a strong increase in binding at 45 min vs 0 min (c), and the majority of sites still display binding at 90 min

Fig. 3

ER occupancy after estradiol treatment in terms of TF network activity. (a) Global TF network behavior as predicted by VULCAN in our ChIP-seq dataset, highlighting the ESR1 TF at time 0 and 45/90 min after estradiol treatment. (b) Global TF activity after estradiol treatment in MCF7 cells, inferred using the METABRIC network, highlighting TFs significantly upregulated at 45 min and 90 min. (c) Global TF activity after estradiol treatment in MCF7 cells, inferred using the METABRIC network, highlighting TFs significantly downregulated at 45 min and 90 min. (d) Global TF activity after estradiol treatment in MCF7 cells, inferred using the METABRIC network, highlighting TFs significantly upregulated at 45 min but not at 90 min. (e) Global TF activity after estradiol treatment in MCF7 cells, inferred using the METABRIC network, highlighting TFs significantly upregulated at 90 min but not at 45 min. (f) Most enriched motif in peaks upregulated at both 45 and 90 min after estradiol treatment, as predicted by HOMER

Fig. 4

Global TF activity after estradiol treatment using different network models. XY scatter showing the TF activity as calculated by VULCAN for our differential ChIP-seq analysis of ER binding at 45 min (a) and at 90 min (b) after stimulation with 100 nM E2. Comparison of the results calculated using the METABRIC (y-axis) and TCGA (x-axis) networks shows consistent results know ER interactors including PGR, RARA, GATA3, and GRHL2. GRHL2 activity is notably enriched against. The regulon of ER is also consistently enriched in both networks. Pearson’s correlation coefficient (PCC) shown along with the significance

Fig. 5

Inferring TF co-occupancy in public datasets with VULCAN. (a) VULCAN activity scores for a few TFs derived from the ER-targeted ChIP-seq breast cancer patient-derived xenograft (PDX) dataset GSE110824. The behavior of ESR1, FOXA1, and GATA3 is correlated, while FOXC1 shows an inversely correlated pattern (blue line). Interestingly, the sample with the lowest Allred score (V0980 U) has the lowest activity and the other luminal markers. (b) VULCAN activity scores for FOXA1 in ChIP-seq experiments targeting the androgen receptor (AR) in LNCaP-1F5 prostate-derived cells (dataset GSE39880). The bar plots show the relative VULCAN normalized enrichment score calculated on absolute peak intensities after treating cells with dihydrotestosterone (DHT) and partial AR modulators cyproterone acetate (CPA) and mifepristone (RU486). FOXA1 network binding is higher in the presence of the strong AR recruiter DHT. This shows an increased FOXA1/AR promoter co-occupancy in DHT-treated cells, in agreement with the conclusions of the study that originated the dataset. Two replicates for each treatment were produced and are reported in matching colors

Fig. 6

GRHL2 differential ChIP-seq between 0 and 45 min. (a) Activation of the ER with estro-2-diol results in a genome-wide increase in GRHL2 binding. (b) VULCAN analysis of the same data shows a significant enrichment for ESR1 sites in both the context of the METABRIC and TGCA networks. The regulon for FOXA1 is also not enriched. Inspection of known FOXA1/GRHL2 sites (e.g., RARa promoter) shows GRHL2 already bound. (c) Overlap of GRHL2 binding with public datasets shows that E2-responsive GRHL2 sites show considerable overlap with ER, FOXA1, and P300 sites; H3K4Me1 and H3K4Me3 show little enrichment. (d) Analysis of P300 binding showed a greater overlap of GRHL2 ER-responsive sites in the presence of E2 than in control conditions. (e) Overlap with ER ChIA-PET sites showed enrichment for GRHL2 sites at ER enhancers. (f) Analysis of Gro-SEQ data (GSE43836) at GRHL2 sites. Blue lines are control samples, pink lines are samples after stimulation with E2. In general, GRHL2 sites (left) show no change in the levels of transcription on the addition of E2; however, E2-responsive GRHL2 sites (right) show a robust increase in transcription on the activation of the ER. (g) Motif analysis of differentially bound sites gave the top two results as GRHL2 and ER

Fig. 7

Estrogen time course and Co-IP of GRHL2. Analysis by western blot of the GRHL2 showed no changes in the levels of GRHL2 at 45 min, 90 min, or 24 h after stimulation with estradiol in either MCF7 or T47D. Co-IP of ER (bait, red, Santa Cruz:sc-8002) identified GRHL2 (green, Atlas: HPA004820) as an interactor in estrogenic conditions (M = marker, I = input, FT = flow through, IP = immunoprecipitation). siRNA knockdown of GRHL2 in MCF7 (right) resulted in a loss of the ~ 75-kDa band

Fig. 8

Effect of GRHL2 knockdown after 24 h on eRNA at E2-responsive binding sites and overexpression of GRHL2 Δ425–437. (a) Overexpression of GRHL2 in MCF7 resulted in a reduction of eRNA transcribed from the GREB1, TFF1, and XBP1 enhancers. The effect was significant at TFF1 and XBP1 enhancers (p < 0.05, paired t test). (b) Overexpression of GRHL2 Δ425–437 (delta) compared to empty vector (EV) and GRHL2 wild type (OE) at 24 h. In all three cell lines at all three loci, overexpression of the wild type (WT) led to a reduction in the mean eRNA production at GREB1, TFF1, and XPB1. This effect was significant in six out of nine experiments (p < 0.05, t test, one-tailed, paired). Overexpression of GRHL2 Δ425–437 had a reduced effect that led to a significant reduction in only two out of nine experiments (p < 0.05, t test, one-tailed, paired). Importantly, in four out of nine experiments, WT overexpression had significantly less eRNA production than GRHL2 Δ425–437, suggesting the P300 inhibition domain plays a role in the regulation of eRNA production

Fig. 9

Changes in H3K27ac on knockdown of GRHL2. (a) The effect of silencing GRHL2 on H3K27ac at 48 h in MCF7 and T47D cell lines was monitored by ChIP-seq. Analysis of sites proximal to TFF1, XBP1, and GREB1 showed significant changes in acetylation at all three sites in MCF7. Significant changes were only found at GREB1 in T47D (top right). While XBP1 and GREB1 show an increase in histone acetylation on silencing GRHL2, TFF1 showed the reverse effect. (b) Genome-wide, the effects of silencing GRHL2 led to a significant redistribution of H3K27ac in both the MCF7 and T47D cell lines, with both showing an increase and decrease in the histone mark dependent on site. c From left to right. Coverage as calculated by Homer. H3K27ac was found at GRHL2 sites in both MCF7 and T47D cells, in particular at the E2-responsive sites. The same mark was also found at P300 sites as expected. Analysis of ER binding at H3K27ac sites showed an enrichment for ER binding at the H3K27ac sites that were most responsive to knockdown of GRHL2 in MCF7 cells

Fig. 10

Overview of the role of GRHL2 in ER activation. On activation of the ER by the ligand E2, the protein is released from a complex containing HSPs and translocates to the nucleus. The holo-ER dimer forms a core complex at estrogen response elements (ERE) with FOXA1 (pioneer factor) and GATA3. ER further recruits P300 and GRHL2. GRHL2 has an inhibitory effect on P300 (a transcriptional activator interacting with TFIID, TFIIB, and RNAPII), thereby reducing the level of eRNA transcription at enhancer sites. Overexpression of GRHL2 further suppresses transcription, while knockdown of GRHL2 reverses the process