Role of Non-Coding Regulatory Elements in the Control of GR-Dependent Gene Expression

Abstract

The glucocorticoid receptor (GR, also known as NR3C1) coordinates molecular responses to stress. It is a potent transcription activator and repressor that influences hundreds of genes. Enhancers are non-coding DNA regions outside of the core promoters that increase transcriptional activity via long-distance interactions. Active GR binds to pre-existing enhancer sites and recruits further factors, including EP300, a known transcriptional coactivator. However, it is not known how the timing of GR-binding-induced enhancer remodeling relates to transcriptional changes. Here we analyze data from the ENCODE project that provides ChIP-Seq and RNA-Seq data at distinct time points after dexamethasone exposure of human A549 epithelial-like cell line. This study aimed to investigate the temporal interplay between GR binding, enhancer remodeling, and gene expression. By investigating a single distal GR-binding site for each differentially upregulated gene, we show that transcriptional changes follow GR binding, and that the largest enhancer remodeling coincides in time with the highest gene expression changes. A detailed analysis of the time course showed that for upregulated genes, enhancer activation persists after gene expression changes settle. Moreover, genes with the largest change in EP300 binding showed the highest expression dynamics before the peak of EP300 recruitment. Overall, our results show that enhancer remodeling may not directly be driving gene expression dynamics but rather be a consequence of expression activation.

Keywords: ChIP-Seq; ENCODE; EP300; GR; NR3C1; RNA-Seq; enhancer sequences; expression regulation; glucocorticoid receptor; histone modifications.

Figure 1

Profile of glucocorticoid receptor (GR)-induced gene expression changes. (A) Schematic representation of the experiment that was performed to generate the data analyzed in this publication. Data were downloaded from the ENCODE project database [6,21]. A549 epithelial-like human cell line was treated with 100 nM dexamethasone, and cells were fixed at distinct time points between 30 min and 12 h. RNA-Seq was performed at each time point. Raw RNA-Seq data were normalized and filtered for false discovery rate (FDR) < 0.0000001. (B) Abundance levels of transcript for each gene were scaled with z-score transformation and leveled to 0 at 0 min time point. Blue cluster—downregulated genes, red cluster—upregulated genes, gray cluster—not regulated genes. Table S1 contains a full list of regulated genes together with the p-value, FDR and cluster information. Table S2 lists the same information for selected random genes. Full raw data from the experiments are available from both ENCODE and GEO databases (see Materials and Methods).

Figure 2

Time-course of average shapes of enhancer peaks associated with NR3C1. (A) Schematic representation of data analysis process. For each of the selected genes, the enhancer region was scanned for the strongest NR3C1 peak (+/− 100 kb with the exclusion of +/− 2 kb). EP300, H3K27ac, and H3K4me1 peaks were extracted from the same locations and averaged across genes from each cluster. (B) Each peak was centered according to the amplitude of NR3C1, and +/− 2 kb surrounding this position are plotted for each time point (average signal across genes). Two-way ANOVA time x regulation cluster: EP300: F = 22.6, p = 2 × 10⁻⁶; H3K27ac: F = 17.35, p = 3.2 × 10⁻⁵; H3K4me1: F = 14.96, p = 1.1 × 10⁻⁴; NR3C1: F = 6.43, p = 0.011 (also see Figures S5 and S6, Statistics: Tables S3, S4 and S5). Stars represent the p-values of post hoc t-tests that compared amplitudes of peaks divided by the amplitude of random peaks in each cluster (upregulated: red stars, downregulated: blue stars) at each time point compared to the same value at 0 min. ** p < 0.01, *** p < 0.001.

Figure 3

Analysis of gene expression dynamics and enhancer remodeling in upregulated genes with NR3C1 binding to enhancer sites. (A) Schematic illustration of weighted maximum time point (MWT) analysis. MWTs were extracted from analyzed ChIP-Seq regulatory factors and gene expression dynamics. For each factor, an average time-weighted by peak amplitude was calculated. For example, an MWT of 1.5 h means that for this factor, the maximum amplitude was observed both as 1 h and 2 h. For expression, MWTs were calculated based on the difference in transcript counts to illustrate expression dynamics. The higher the MWT, the later maximum binding/dynamics occur. (B) Boxplots of MWTs for all upregulated genes with identified enhancer peaks for each of the regulatory factors measured. For comparison of weighted time points, one-way ANOVA (F = 178.5, p = 2 × 10⁻¹⁶) with post hoc pairwise tests was used. Gene expression and NR3C1 weighted time points were tested against all other weighted time points with Bonferroni correction (expression vs. EP300:0.303, all other test p-values < 1 × 10⁻⁴, (See Tables S6 and S7 for details). (C) Boxplot of z-scored amplitudes of ChIP-Seq peaks for NR3C1 and EP300 during the time course with fitted 3rd-degree polynomial lines. (D) Boxplot of gene expression counts represented as log2(count/(mean count at 0 min)). The fitted line does not represent changes in counts linearly as it is fitted to show expression dynamics. The scales of the boxplots and the fitted line are not directly comparable. Fitted line: A 3rd-degree polynomial fit to the derivative of gene expression scaled by a factor of 100 to equalize the scales). ** p < 0.01, *** p < 0.001.

Figure 4

Analysis of gene expression dynamics and enhancer remodeling in upregulated genes with the highest changes in EP300 binding to enhancer sites. (A) Overview of the analysis steps. First, NR3C1 peaks were located. Next, EP300 peak amplitudes were collected. A list of 100 genes closest to the peaks with maximum EP300 changes was identified, and the upregulated genes from this list (n = 85) were subjected to MWT analysis (see Figure 3A). (B) Boxplots of MWTs for genes with largest EP300-binding changes (colored bars) compared with the data for upregulated genes (gray bars, same data as Figure 3B, gray bars). To compare weighted time points, two-way ANOVA with post hoc pairwise tests was used (See Table S9 and S9 for details). ** p < 0.01, *** p < 0.001 (p < 0.05 is not marked for clarity).