Next Article in Journal
Mangiferin Inhibits Apoptosis in Doxorubicin-Induced Vascular Endothelial Cells via the Nrf2 Signaling Pathway
Next Article in Special Issue
Identification of Rhythmically Expressed LncRNAs in the Zebrafish Pineal Gland and Testis
Previous Article in Journal
Skeletal and Cardiac Muscle Disorders Caused by Mutations in Genes Encoding Intermediate Filament Proteins
Previous Article in Special Issue
Fish-Ing for Enhancers in the Heart
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
IJMS
Volume 22
Issue 8
10.3390/ijms22084258
ijms-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Malgorzata Borczyk
*,
Mateusz Zieba
*,
Michał Korostyński
and
Marcin Piechota
Laboratory of Pharmacogenomics, Department of Molecular Neuropharmacology, Maj Institute of Pharmacology, Polish Academy of Sciences, 31-343 Krakow, Poland
*
Authors to whom correspondence should be addressed.
Int. J. Mol. Sci. 2021, 22(8), 4258; https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22084258
Submission received: 15 March 2021 / Revised: 6 April 2021 / Accepted: 8 April 2021 / Published: 20 April 2021
(This article belongs to the Special Issue Functions of Non-coding DNA Regions)
Browse Figures
Review Reports Versions Notes

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.

1. Introduction

Glucocorticoid receptor (GR), encoded by the NR3C1 gene, is a ligand-dependent transcriptional regulator. After binding to glucocorticoids, it translocates to the nucleus and interacts with the genome to influence gene expression. GR stimulation activates and represses hundreds of genes [1]. GR signaling plays an important role in the regulation of neuronal and glial cells, including their responses to drugs [2,3]. ChiP-Seq experiments allow for high-throughput analysis of chromatin binding, and studies with this method have brought great progress in understanding GR binding to DNA. The specificity of GR binding depends on pre-existing chromatin landscape, with factors, such as accessibility and motifs promoting the binding [4,5,6]. GR frequently binds to enhancers-cis-regulatory elements distal to transcription start site (TSS) of both repressed and activated genes [6,7].
Although it is difficult to define the precise features of enhancers, they are often associated with specific histone modifications and the presence of particular proteins [8]. Histone methylations and acetylations mark active and dormant enhancers. Monomethylation of lysine 4 and acetylation of lysine 27 of histone H3 (H3K4me1 and H3K27ac) mark active enhancers, whereas enhancers lacking H3K27ac are poised [9]. Another marker of active enhancers is the Histone acetyltransferase EP300 (also known as p300) that functions as a transcription coactivator [9,10]. Apart from enhancers, moderate EP300 enrichment is found in core promoters [11]. H3K4me1, on the other hand, is found predominantly in distal elements [9]. Enhancer sequences are also highly conserved [12,13]. Once an enhancer is active, it works to activate transcription of its targets via chromatin loops. Millions of such loops exist in the human genome [14]. The important physiological role of enhancers was confirmed by genomic studies that found many disease-associated polymorphisms in enhancer sequences [15,16].
GR effects on gene expression have been thoroughly investigated in several studies and extensively reviewed [17,18,19,20]. However, it is not known how enhancer remodeling induced by NR3C1 connects to gene expression dynamics. The interplay between NR3C1 binding, enhancer remodeling, and gene expression changes remain not thoroughly investigated, particularly at time points in the range of hours. Right after glucocorticoid (GC) treatment, NR3C1 binds to chromatin, and maximal binding in vitro can be detected after 30 min of GC exposure [6]. GR binds predominantly to distal sites that are typically already predetermined and marked as active enhancers [4,6,7]. Binding is followed by further enhancer remodeling by redistribution of EP300. In a recent study, McDowell et al. showed that changes in transcription factor (TF) binding and histone modifications within enhancers that occur after dexamethasone (synthetic GC) treatment correlates with gene expression changes [6]. They also showed that NR3C1 binding and EP300 mutually reinforce each other. Still, an outstanding question remains: how does NR3C1 binding to enhancers, EP300 recruitment, and histone modifications relate to the timing of transcriptional changes?
Here we present results of an investigation of the interplay between gene expression, GR binding, and remodeling of the enhancer landscape. The dataset used for this research was obtained from the Reddy Laboratory through the ENCODE project and includes both RNA-Seq and ChIP-Seq data of dexamethasone-treated human epithelial-like cell line collected at different time points between 30 min and 12 h [6,21]. The analysis started with the identification of genes robustly regulated by NR3C1. For each gene, dynamics of a single NR3C1-binding site outside of the promoter regions were analyzed. Apart from NR3C1 binding, we focused on markers of enhancer remodeling: EP300, H3K4me1, and H3K27ac. Time-course analysis of upregulated genes showed that in selected enhancers majority of NR3C1 binding occurs within the first hours of GC treatment, which is then followed by an increase in H3K27ac, a marker of active enhancers. Maximum expression dynamics occur ~1 h after the maximal NR3C1 binding is reached and nearly coincide with the largest EP300 recruitment. Additionally, genes with the largest change in EP300 binding were investigated and showed the highest expression dynamics before peak EP300 recruitment. Overall, our results show that enhancer remodeling may not directly be driving gene expression dynamics but rather be a consequence of or co-occur with expression activation.

2. Results

2.1. GR Regulates Expression Bidirectionally

Gene expression changes were first evaluated to indicate differentially regulated transcripts. As GR regulates hundreds of genes, the goal of this part of the analysis was to select a smaller subset of genes with clear and robust expression changes after dexamethasone [6]. After applying a detection threshold of FDR < 0.0000001, we identified 374 upregulated and 371 downregulated genes (Figure 1 and Figure S1, Table S1). An extra cluster of random genes with similar gene expression levels and transcript lengths was used for comparison (Figures S2 and S3). Enrichr was used to identify published gene-sets with overlapping genes [22]. The list of significantly overlapping groups of genes includes 346 published gene-sets (adjusted p-value < 0.05). In the top part of the list, groups of steroid-activated genes, as well as transcripts regulated in hormone-dependent cancers, were identified [23,24,25,26]. What is more, similar to others, we observe that GR activation leads to bidirectional gene expression changes and that the number of up- and downregulated genes are roughly equal [1,6].

2.2. Gene Activation Is Associated with Enhancer Remodeling

We investigated the relationship between the observed gene expression changes and enhancers’ activity in the vicinity of those genes. Enhancer peaks were defined as ChIP-Seq peaks occurring within −/+ 100 kbp (kilo base-pair) from the TSS with the exclusion of the promoter region (−/+ 2 kb surrounding the TSS). If a gene had more than one peak, the one with the strongest NR3C1 signal was chosen. Two hundred thirty-six peaks of upregulated genes were identified at 1 h of dexamethasone treatment. Thus, 63.8% of the upregulated genes had NR3C1 signals. One hundred one peaks were identified near downregulated genes (in 26.9% of genes), and 183 out of 743 random genes had NR3C1 peaks (24%). The distribution of identified enhancer peaks showed that in the case of upregulated genes, the highest peaks tend to occur closer to the promoter region than for the other gene clusters (Figure S4). Still, these detected peaks are well outside the core promoters, defined as a region of a few hundred base pairs surrounding the TSS [27].
Next, we identified signals for EP300, H3K4me1, and H3K27ac in the same range as discovered peaks. Figure 2 shows the time course of averaged ChIP-Seq signal in identified peaks. For NR3C1, EP300, and H3K27ac, the signal was strongest in peaks near upregulated genes throughout the time course (two-way ANOVA, time:regulation, p < 0.02 for all four regulatory factors, Supplementary Table S3). EP300 binding near upregulated showed a notable decrease in signal at the 4 h timepoint. This could indicate overlapping of distinct patterns of EP300 binding, but individual time-courses did not confirm this observation (Figure S6). As downregulated genes did not show robust changes in the ChIP-Seq signal for any of the analyzed regulatory factors, they were not included in further analyses (Figure 2 and Figure S5, Tables S4 and S5). For all analyzed regulatory factors, similar dynamics were observed at promoter regions (Figures S7 and S8).

2.3. The Largest EP300 Recruitment Coincides with the Most Dynamic Gene Expression Changes

We interrogated the timing of gene expression dynamics concerning enhancer remodeling for upregulated genes. For gene expression, the time point of the largest dynamics was calculated as a weighted time average, where the first derivative of normalized transcript abundance was used as weights. For each of the analyzed ChIP-Seq tracks, maximum binding was interrogated by calculating a weighted time point with amplitudes at each time used as weights. The median weighted time point of maximum gene expression changes was 5 h and 27 min. Although increased and decreased gene expression is visible after 30 min of dexamethasone treatment, the changes do not reach their maximum dynamics for hours. The median weighted time point of maximum NR3C1 binding was 4 h 36 min. For H3K27ac and H3K4me1, the time points were 4 h 59 min and 4 h 48 min, respectively. Thus, gene expression changes reach maximum dynamics 51 min after maximum NR3C1 recruitment, and H3K27ac and H3K4me1 maximum amplitudes are observed within this interval. Most strikingly, the median of maximum EP300 binding nearly coincided in time with maximum changes in gene expression and occurred at 5 h 33 min. Thus, maximum expression dynamics occur together with the highest EP300 recruitment to peaks occupied by NR3C1 and are preceded with the highest NR3C1 binding and histone modifications (Figure 3).

2.4. Peaks with the Highest Changes in EP300 Recruitment Occur near Upregulated Genes and EP300 Is Recruited to Them after the Largest Expression Dynamics Are Reached

To further analyze the relationship between dexamethasone treatment-induced enhancer activity and gene expression, genes with the biggest change in EP300 binding in NR3C1 peaks were identified. First, a list of all identified NR3C1 peaks in enhancer regions was prepared. Then top 100 peaks with the largest difference in EP300 binding between any two time points starting from 30 min were selected. Then, for each of those peaks, the nearest gene was found. In the case two peaks were attributed to the same gene, the peak with the higher amplitude was retained. Expression dynamics of these genes were then analyzed, and 46 out of 100 genes were differentially expressed (FDR < 0.1), and 23 genes passed the previously applied 0.0000001 FDR threshold. Here, again genes were divided into up- and downregulated clusters (15 downregulated and 85 upregulated, Figure S9). We investigated the time points of the highest dynamics of enhancer remodeling of upregulated genes with the highest EP300-binding differences and compared them with upregulated genes from the gene expression-based analysis. Upregulated genes that showed the largest differences in EP300 binding in their strongest NR3C1 peak showed the largest expression dynamics 4 h 54 min and the highest EP300 binding at 5 h 38 min (Figure 4 and Figure S10, Tables S9 and S9). Thus, for the upregulated genes with the highest EP300 dynamics, the time that elapses between maximum expression changes and enhancer remodeling is even more uncoupled than for upregulated genes chosen based on their expression changes.

3. Discussion

Here we present results of an investigation of the time course of gene expression and chromatin-binding changes after dexamethasone treatment in a human epithelial-like cell line [21]. In the first part of the analysis, 374 genes upregulated after dexamethasone treatment were identified, while at the same FDR threshold, 371 genes were downregulated. These results are consistent with the original analysis of gene expression in the same dataset, although here, a more stringent statistical threshold was applied; thus, the number of genes in each cluster is lower than in the previous analysis [6]. For each gene, we identified a single ChIP-seq peak in the enhancer region based on the strongest NR3C1 amplitude. Each of the identified peaks was analyzed in terms of EP300 and NR3C1 binding and two histone modifications H3K27ac and H3K4me1. Next, enhancer-sited dynamics were interrogated in relation to gene expression dynamics and the time course of alterations.
Enhancer remodeling was analyzed for genes upregulated after dexamethasone as downregulated genes did not exhibit robust changes, neither in NR3C1 nor EP300 binding when compared to a randomly selected list of genes. This observation is in concert with literature [28]. The analysis of temporal order of changes in transcription factor binding and histone modifications revealed that there is no point in which all analyzed regulatory factors are orchestrated and acting as one machinery. It is rather a chain of subsequent events that could act through multiple mechanisms. Although strongly supported by the data, the above conclusions have limitations as they describe average dynamics of gene regulation and not individual time courses. We suggest that the binding of NR3C1 to its response elements activates at least two separate pathways. First, direct activation of expression, and the second, starting an orchestrated sequence of chromatin-associated events. The strongest binding of NR3C1 occurs just after dexamethasone treatment, and at this moment, transcription of target genes increases. However, at the following time points, binding of NR3C1 decreases, but the dynamics of changes in expression increase. The latter is associated with an increase in EP300 binding and H3K27 acetylation.
The association of elevated expression and increased binding of EP300 to enhancers made it the perfect marker of active enhancers, and EP300 is considered to be a coactivator of transcription [9,10]. As a result, most studies concentrate on the requirement of EP300 in the process of direct transcriptional regulation [6]. However, the results of this study suggest that active transcription is followed by EP300 recruitment. This observation adds another layer of meaning to the tight association of EP300 and enhancer-driven expression and shows that recruitment of EP300 might not be as crucial as previously thought. The enhancer-mediated transcriptional activity requires EP300-dependent H3K27 acetylation [29]. However, it might also be acetylated through other redundant mechanisms [29,30]. This may explain why mutations in EP300 lead to milder forms of Rubinstein–Taybi syndrome than mutations in CREBBP [31,32]. This discovery also points to the potential for a positive feedback loop. Activated enhancers may recruit EP300 and, in consequence, become more sensitive to further activations. More EP300 would promote further binding of TFs, including GR. Positive feedback loops were previously suggested as one of the cell differentiation mechanisms [33].
Close temporal association between EP300 binding and acetylation of histone H3 was previously shown [34]. However, for genes downregulated by dexamethasone, this association was not observed. While EP300 binding increases mildly even for downregulated genes, acetylation of H3K27 sites decreases. This suggests that for downregulated genes, there is a robust mechanism driving deacetylation that overcomes EP300-dependent acetylation. H3K27 acetylation is also, out of the four analyzed alterations, the one that directly precedes dynamics of changes in gene expression.
Primed enhancers are marked by H3K4me1 [35]. They are further activated by acetylation of H3K27 through, among others, EP300 acetyltransferase. However, the inverse mechanism was also postulated through which transcription factor recruitment can result in H3K4me1 priming [36,37]. Here, we observe a late increase in H3K4me1 after activation of NR3C1. This immediately suggests a possible mechanism of steroid-dependent cell differentiation in which repeatedly activated NR3C1 responsive sites would be primed by modifying the pattern of accessible enhancers for a particular cell type. This is in concert with the view of H3K4me1 as a fine-tuner of enhancer activity [38,39,40].

4. Methods

All of the code used to download, preprocess and analyze data are available in the project’s GitHub repository: https://github.com/ippas/ifpan-chipseq-timecourse (accessed on 10 March 2021).

4.1. Datasets

RNA-Seq (.tsv files) and ChIP-Seq (.bed and .bigwig files) data were downloaded from GEO database (https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/geo/query/acc.cgi?acc=GSE91222 (accessed on 10 March 2021) and https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/geo/query/acc.cgi?acc=GSE113464 (accessed on 10 March 2021); BioProject link: https://0-www-ncbi-nlm-nih-gov.brum.beds.ac.uk/sra?linkname=bioproject_sra_all&from_uid=356880 (accessed on 10 March 2021)) [21]. For ChIP-Seq, the following tracks were used in the analyses: NR3C1 (GR), EP300, H3K4me1, and H3K27ac.

4.2. RNA-Seq Analysis

RNA-Seq analysis was performed in R 3.4. First, quantile normalization was performed on raw counts. After normalization, one-way ANOVA with FDR correction was used to evaluate gene expression changes. Genes with FDR < 0.0000001 were included into the cluster analysis (n = 745 genes). Genes were classified as either up- or downregulated by dexamethasone. Additional random (control) gene-set of genes with similar gene-expression levels to regulated genes was selected. For dynamic expression analyses (Figure 3 and Figure 4), outliers in terms of gene expression counts were removed according to the boxplot.stats function in R. Weighted time point of maximum gene expression changes was calculated as a count-change-weighted average of time.

4.3. ChIP-Seq analysis

For each TSS, defined as the gene start in the Ensembl database, lists of surrounding (+/− 100,000 bp) peaks for NR3C1 were selected from *.bed files. Peaks common for samples within each experimental group were obtained using bedtools intersect. The peaks were further divided into core promoter regions (within +/− 2000 bp from TSS) and enhancer regions (excluding +/− 2000 bp from TSS). Within the enhancer region, a peak with the highest amplitude for each gene in 60 min time point after dexamethasone treatment was considered for further analysis. This time point showed the highest overall GR binding signal. For each analyzed NR3C1 peak, ChIP-Seq coverage for EP300, H3K4me1 and H3K27ac were analyzed within the coordinates of each NR3C1 peak using .bigwig data taking maximum coverage as summary value. For EP300 dynamics association with expression, a list of all identified NR3C1 peaks (60 min time point) in enhancer regions was prepared. Next, coverage of EP300 for those peaks was obtained using bigWigAverageOverBed from the kentUtils package. The top 100 peaks with the largest difference in EP300 binding between any two time points starting from 30 min were selected. Then, for each of those peaks, the nearest gene was identified. In case two peaks were attributed to the same gene, the peak with the higher amplitude was retained. For the temporal analysis, the time point of maximum binding was defined as an amplitude-weighted average of time. For visualization of ChIP-Seq coverage, .bigwig files were used and summarized into 10 bp buckets.

5. Conclussions

In summary, this analysis focuses on temporal enhancer and gene expression activation patterns that lead to novel observations. The highest NR3C1 recruitment is not directly associated with the highest dynamic in gene expression but rather precedes it. Correspondingly the highest dynamic in gene expression preceded the highest recruitment of EP300, which was unexpected considering the role of EP300 in enhancers activity. Furthermore, while gene expression dynamics declined, EP300 binding to chromatin stayed at high levels. In this orchestrated cascade of events, we observed an increase of H3K27 acetylation in NR3C1 responsive enhancers before maximum EP300 recruitment was reached and followed by increased H3K4me1 priming. These observations altogether show complicated temporal orchestration of various transcription mechanisms. The observations also point to possible positive feedback loops that might be involved in transcription regulation in case of response to repeated molecular stimuli. In such loops, initial NR3C1 binding would induce gene expression and chromatin remodeling as two distinct phenomena. The chromatin remodeling would then serve as a priming factor and lead to distinct responses to consecutive stimuli. This phenomenon is not well investigated. A recent study showed that GC exposure during neurogenesis leads to enhanced transcriptional responses to subsequent GC challenges [41]. Chromatin priming by NR3C1 could explain the role GC priming plays in regulating neuronal and glial responses to stress [42]. These priming loops could be involved in more complex processes, such as cell differentiation in which GR is widely known to be involved [43].

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/ijms22084258/s1, Figure S1: Heatmap of gene expression of selected differentially-expressed genes; Figure S2 Mean abundance (RPKM) histogram of analysed genes.; Figure S3 transcript length histogram for analysed genes.; Figure S4 Histogram of GR strongest enhancer peak amplitude distance from transcription start site (TSS).; Figure S5 Line-plot normalised to random genes for enhancer peaks; Figure S6 Heatmap of data presented in Figure 2; Figure S7 Line-plot for promotor peaks (axes +/- 10,000 from TSS). Here peaks were not centered, but simply average signal was plotted for the chosen interval.; Figure S8 Line-plot normalised to random genes for promotor peaks (axes +/- 10,000 from TSS); Figure S9 Heatmap of gene expression of genes with highest EP300 changes; Figure S9 Heatmap of gene expression of genes with highest EP300 changes; Figure S10 Heatmap of gene ChIP-seq peaks for genes with highest EP300 changes; Table S1 Results of ANOVA on gene expression with the full list of regulated genes together with the raw data, p-value, FDR and cluster information; Table S2 Results of ANOVA on gene expression with the full list of randomly selected genes together with the raw data, p-value, FDR and cluster information.; Table S3 Two-way ANOVA (time x gene regulation cluster) results based on amplitudes of normalised peaks (data presented in Figure S2.2); Table S4 Post-hoc pairwise t-tests of amplitudes of normalised peaks, comparisons between each timepoint ant 0 min; Table S5 Post-hoc pairwise t-tests of amplitudes of normalised peaks, comparisons between upregulated and downregulated genes; Table S6 one-way ANOVA results of weighted timepoints (data presented in Figure 3A); Table S7 Post-hoc results of weighted timepoints (data presented in Figure 3A); Table S8 two-way ANOVA results of weighted timepoints presented in Figure 4; Table S9 Post-hoc results of weighted timepoints presented in Figure 3A

Author Contributions

M.P. designed the analyses; M.Z. and M.P. conducted the data analyses, M.B. and M.Z. prepared figures; M.B. and M.P. drafted the manuscript; M.K. and M.Z. provided essential revisions to the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

The work was funded by statutory funds of the Maj Institute of Pharmacology and the following projects awarded by the National Science Center in Poland: APARIC 2015/16/W/ST5/00005 and OPUS 2017/27/B/NZ7/00204.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data used in this study are available via the ENCODE project.

Acknowledgments

The authors would like to acknowledge the ENCODE project, Professor Tim Reddy and his team for providing the underlying dataset.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Reddy, T.E.; Pauli, F.; Sprouse, R.O.; Neff, N.F.; Newberry, K.M.; Garabedian, M.J.; Myers, R.M. Genomic determination of the glucocorticoid response reveals unexpected mechanisms of gene regulation. Genome Res. 2009, 19, 2163–2171. [Google Scholar] [CrossRef] [Green Version]
  2. Piechota, M.; Korostynski, M.; Golda, S.; Ficek, J.; Jantas, D.; Barbara, Z.; Przewlocki, R. Transcriptional signatures of steroid hormones in the striatal neurons and astrocytes. BMC Neurosci. 2017, 18, 1–9. [Google Scholar] [CrossRef] [Green Version]
  3. Slezak, M.; Korostynski, M.; Gieryk, A.; Golda, S.; Dzbek, J.; Piechota, M.; Wlazlo, E.; Bilecki, W.; Przewlocki, R. Astrocytes are a neural target of morphine action via glucocorticoid receptor-dependent signaling. Glia 2013, 61, 623–635. [Google Scholar] [CrossRef]
  4. D’Ippolito, A.M.; McDowell, I.C.; Barrera, A.; Hong, L.K.; Leichter, S.M.; Bartelt, L.C.; Vockley, C.M.; Majoros, W.H.; Safi, A.; Song, L.; et al. Pre-established Chromatin Interactions Mediate the Genomic Response to Glucocorticoids. Cell Syst. 2018, 7, 146–160.e7. [Google Scholar] [CrossRef] [Green Version]
  5. John, S.; Sabo, P.J.; Thurman, R.E.; Sung, M.-H.; Biddie, S.C.; Johnson, T.A.; Hager, G.L.; Stamatoyannopoulos, J.A. Chromatin accessibility pre-determines glucocorticoid receptor binding patterns. Nat. Genet. 2011, 43, 264–268. [Google Scholar] [CrossRef]
  6. McDowell, I.C.; Barrera, A.; D’Ippolito, A.M.; Vockley, C.M.; Hong, L.K.; Leichter, S.M.; Bartelt, L.C.; Majoros, W.H.; Song, L.; Safi, A.; et al. Glucocorticoid receptor recruits to enhancers and drives activation by motif-directed binding. Genome Res. 2018, 28, 1272–1284. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Kuznetsova, T.; Wang, S.-Y.; Rao, N.A.; Mandoli, A.; Martens, J.H.A.; Rother, N.; Aartse, A.; Groh, L.; Janssen-Megens, E.M.; Li, G.; et al. Glucocorticoid receptor and nuclear factor kappa-b affect three-dimensional chromatin organization. Genome Biol. 2015, 16, 264. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Shlyueva, D.; Stampfel, G.; Stark, A. Transcriptional enhancers: From properties to genome-wide predictions. Nat. Rev. Genet. 2014, 15, 272–286. [Google Scholar] [CrossRef] [PubMed]
  9. Creyghton, M.P.; Cheng, A.W.; Welstead, G.G.; Kooistra, T.G.; Carey, B.W.; Steine, E.J.; Hanna, J.; Lodato, M.A.; Frampton, G.M.; Sharp, P.A.; et al. Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc. Natl. Acad. Sci. USA 2010, 107, 21931–21936. [Google Scholar] [CrossRef] [Green Version]
  10. Visel, A.; Blow, M.J.; Li, Z.; Zhang, T.; Akiyama, J.A.; Holt, A.; Plajzer-Frick, I.; Shoukry, M.; Wright, C.; Chen, F.; et al. ChIP-seq accurately predicts tissue-specific activity of enhancers. Nat. Cell Biol. 2009, 457, 854–858. [Google Scholar] [CrossRef] [Green Version]
  11. Heintzman, N.D.; Stuart, R.K.; Hon, G.; Fu, Y.; Ching, C.W.; Hawkins, R.D.; Barrera, L.O.; Van Calcar, S.; Qu, C.; Ching, K.A.; et al. Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat. Genet. 2007, 39, 311–318. [Google Scholar] [CrossRef] [PubMed]
  12. Chen, L.; Fish, A.E.; Capra, J.A. Prediction of gene regulatory enhancers across species reveals evolutionarily conserved sequence properties. PLoS Comput. Biol. 2018, 14, e1006484. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  13. Clément, Y.; Torbey, P.; Gilardi-Hebenstreit, P.; Crollius, H.R. Enhancer–gene maps in the human and zebrafish genomes using evolutionary linkage conservation. Nucleic Acids Res. 2020, 48, 2357–2371. [Google Scholar] [CrossRef] [PubMed]
  14. Jin, F.; Li, Y.; Dixon, J.R.; Selvaraj, S.; Ye, Z.; Lee, A.Y.; Yen, C.-A.; Schmitt, A.D.; Espinoza, C.A.; Ren, B. A high-resolution map of the three-dimensional chromatin interactome in human cells. Nature 2013, 503, 290–294. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Kvon, E.Z.; Zhu, Y.; Kelman, G.; Novak, C.S.; Plajzer-Frick, I.; Kato, M.; Garvin, T.H.; Pham, Q.; Harrington, A.N.; Hunter, R.D.; et al. Comprehensive in vivo interrogation reveals phenotypic impact of human enhancer variants. Cell 2020, 180, 1262–1271. [Google Scholar] [CrossRef]
  16. Zhang, G.; Shi, J.; Zhu, S.; Lan, Y.; Xu, L.; Yuan, H.; Liao, G.; Liu, X.; Zhang, Y.; Xiao, Y.; et al. DiseaseEnhancer: A resource of human disease-associated enhancer catalog. Nucleic Acids Res. 2018, 46, D78–D84. [Google Scholar] [CrossRef] [Green Version]
  17. Ramamoorthy, S.; Cidlowski, J.A. Corticosteroids. Rheum. Dis. Clin. N. Am. 2016, 42, 15–31. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  18. Ratman, D.; Berghe, W.V.; Dejager, L.; Libert, C.; Tavernier, J.; Beck, I.M.; De Bosscher, K. How glucocorticoid receptors modulate the activity of other transcription factors: A scope beyond tethering. Mol. Cell. Endocrinol. 2013, 380, 41–54. [Google Scholar] [CrossRef]
  19. Sacta, M.A.; Chinenov, Y.; Rogatsky, I. Glucocorticoid Signaling: An Update from a Genomic Perspective. Annu. Rev. Physiol. 2016, 78, 155–180. [Google Scholar] [CrossRef] [PubMed]
  20. Weikum, E.R.; Knuesel, M.T.; Ortlund, E.R.W.E.A.; Yamamoto, M.T.K.K.R. Glucocorticoid receptor control of transcription: Precision and plasticity via allostery. Nat. Rev. Mol. Cell Biol. 2017, 18, 159–174. [Google Scholar] [CrossRef]
  21. ENCODE Project Consortium. An integrated encyclopedia of DNA elements in the human genome. Nature 2012, 489, 57–74. [Google Scholar] [CrossRef] [PubMed]
  22. Chen, E.Y.; Tan, C.M.; Kou, Y.; Duan, Q.; Wang, Z.; Meirelles, G.V.; Clark, N.R.; Ma’Ayan, A. Enrichr: Interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinform. 2013, 14, 128. [Google Scholar] [CrossRef] [Green Version]
  23. Farmer, P.; Bonnefoi, H.; Becette, V.; Tubiana-Hulin, M.; Fumoleau, P.; Larsimont, D.; MacGrogan, G.; Bergh, J.; Cameron, D.A.; Goldstein, D.R.; et al. Identification of molecular apocrine breast tumours by microarray analysis. Oncogene 2005, 24, 4660–4671. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Hendriksen, P.J.; Dits, N.F.; Kokame, K.; Veldhoven, A.; Van Weerden, W.M.; Bangma, C.H.; Trapman, J.; Jenster, G. Evolution of the Androgen Receptor Pathway during Progression of Prostate Cancer. Cancer Res. 2006, 66, 5012–5020. [Google Scholar] [CrossRef] [Green Version]
  25. Mikkonen, L.; Pihlajamaa, P.; Sahu, B.; Zhang, F.-P.; Jänne, O.A. Androgen receptor and androgen-dependent gene expression in lung. Mol. Cell. Endocrinol. 2010, 317, 14–24. [Google Scholar] [CrossRef] [Green Version]
  26. Wirapati, P.; Sotiriou, C.; Kunkel, S.; Farmer, P.; Pradervand, S.; Haibe-Kains, B.; Desmedt, C.; Ignatiadis, M.; Sengstag, T.; Schütz, F.; et al. Meta-analysis of gene expression profiles in breast cancer: Toward a unified understanding of breast cancer subtyping and prognosis signatures. Breast Cancer Res. 2008, 10, 1–11. [Google Scholar] [CrossRef]
  27. Liu, R.; States, D.J. Consensus Promoter Identification in the Human Genome Utilizing Expressed Gene Markers and Gene Modeling. Genome Res. 2002, 12, 462–469. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  28. Vandevyver, S.; Dejager, L.; Tuckermann, J.; Libert, C. New Insights into the Anti-inflammatory Mechanisms of Glucocorticoids: An Emerging Role for Glucocorticoid-Receptor-Mediated Transactivation. Endocrinology 2013, 154, 993–1007. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Raisner, R.; Kharbanda, S.; Jin, L.; Jeng, E.; Chan, E.; Merchant, M.; Haverty, P.M.; Bainer, R.; Cheung, T.; Arnott, D.; et al. Enhancer Activity Requires CBP/P300 Bromodomain-Dependent Histone H3K27 Acetylation. Cell Rep. 2018, 24, 1722–1729. [Google Scholar] [CrossRef] [Green Version]
  30. De Almeida Nagata, D.E.; Chiang, E.Y.; Jhunjhunwala, S.; Caplazi, P.; Arumugam, V.; Modrusan, Z.; Chan, E.; Merchant, M.; Jin, L.; Arnott, D.; et al. Regulation of Tumor-Associated Myeloid Cell Activity by CBP/EP300 Bromodomain Modulation of H3K27 Acetylation. Cell Rep. 2019, 27, 269–281.e4. [Google Scholar] [CrossRef] [Green Version]
  31. Fergelot, P.; Van Belzen, M.; Van Gils, J.; Afenjar, A.; Armour, C.M.; Arveiler, B.; Beets, L.; Burglen, L.; Busa, T.; Collet, M.; et al. Phenotype and genotype in 52 patients with Rubinstein-Taybi syndrome caused byEP300mutations. Am. J. Med Genet. Part A 2016, 170, 3069–3082. [Google Scholar] [CrossRef]
  32. Zimmermann, N.; Acosta, A.M.B.F.; Kohlhase, J.; Bartsch, O. Confirmation of EP300 gene mutations as a rare cause of Rubinstein–Taybi syndrome. Eur. J. Hum. Genet. 2007, 15, 837–842. [Google Scholar] [CrossRef] [Green Version]
  33. Becskei, A.; Séraphin, B.; Serrano, L. Positive feedback in eukaryotic gene networks: Cell differentiation by graded to binary response conversion. EMBO J. 2001, 20, 2528–2535. [Google Scholar] [CrossRef] [Green Version]
  34. Zhang, B.; Day, D.S.; Ho, J.W.; Song, L.; Cao, J.; Christodoulou, D.; Seidman, J.G.; Crawford, G.E.; Park, P.J.; Pu, W.T. A dynamic H3K27ac signature identifies VEGFA-stimulated endothelial enhancers and requires EP300 activity. Genome Res. 2013, 23, 917–927. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Wang, C.; Lee, J.-E.; Lai, B.; Macfarlan, T.S.; Xu, S.; Zhuang, L.; Liu, C.; Peng, W.; Ge, K. Enhancer priming by H3K4 methyltransferase MLL4 controls cell fate transition. Proc. Natl. Acad. Sci. USA 2016, 113, 11871–11876. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Sérandour, A.A.; Avner, S.; Percevault, F.; DeMay, F.; Bizot, M.; Lucchetti-Miganeh, C.; Barloy-Hubler, F.; Brown, M.; Lupien, M.; Métivier, R.; et al. Epigenetic switch involved in activation of pioneer factor FOXA1-dependent enhancers. Genome Res. 2011, 21, 555–565. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  37. Sharifi-Zarchi, A.; Gerovska, D.; Adachi, K.; Totonchi, M.; Pezeshk, H.; Taft, R.J.; Schöler, H.R.; Chitsaz, H.; Sadeghi, M.; Baharvand, H.; et al. DNA methylation regulates discrimination of enhancers from promoters through a H3K4me1-H3K4me3 seesaw mechanism. BMC Genom. 2017, 18, 1–21. [Google Scholar] [CrossRef] [PubMed]
  38. Dorighi, K.M.; Swigut, T.; Henriques, T.; Bhanu, N.V.; Scruggs, B.S.; Nady, N.; Still, C.D., II; Garcia, B.A.; Adelman, K.; Wysocka, J. Mll3 and Mll4 facilitate enhancer RNA synthesis and transcription from promoters independently of H3K4 monomethylation. Mol. Cell 2017, 66, 568–576.e4. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Local, A.; Huang, H.; Albuquerque, C.P.; Singh, N.; Lee, A.Y.; Wang, W.; Wang, C.; Hsia, J.E.; Shiau, A.K.; Ge, K.; et al. Identification of H3K4me1-associated proteins at mammalian enhancers. Nat. Genet. 2018, 50, 73–82. [Google Scholar] [CrossRef] [PubMed]
  40. Rickels, R.; Herz, H.-M.; Sze, C.C.; Cao, K.; Morgan, M.A.; Collings, C.K.; Gause, M.; Takahashi, Y.; Wang, L.; Rendleman, E.J. Histone H3K4 monomethylation catalyzed by Trr and mammalian COMPASS-like proteins at enhancers is dispensable for development and viability. Nat. Genet. 2017, 49, 1647. [Google Scholar] [CrossRef] [Green Version]
  41. Provençal, N.; Arloth, J.; Cattaneo, A.; Anacker, C.; Cattane, N.; Wiechmann, T.; Röh, S.; Ködel, M.; Klengel, T.; Czamara, D.; et al. Glucocorticoid exposure during hippocampal neurogenesis primes future stress response by inducing changes in DNA methylation. Proc. Natl. Acad. Sci. USA 2020, 117, 23280–23285. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Sorrells, S.F.; Caso, J.R.; Munhoz, C.D.; Sapolsky, R.M. The Stressed CNS: When Glucocorticoids Aggravate Inflammation. Neuron 2009, 64, 33–39. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  43. Chapman, A.B.; Knight, D.M.; Ringold, G.M. Glucocorticoid regulation of adipocyte differentiation: Hormonal triggering of the developmental program and induction of a differentiation-dependent gene. J. Cell Biol. 1985, 101, 1227–1235. [Google Scholar] [CrossRef] [PubMed]
Ijms 22 04258 g001 550
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 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).
Ijms 22 04258 g001
Ijms 22 04258 g002 550
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−6; H3K27ac: F = 17.35, p = 3.2 × 10−5; H3K4me1: F = 14.96, p = 1.1 × 10−4; 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 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−6; H3K27ac: F = 17.35, p = 3.2 × 10−5; H3K4me1: F = 14.96, p = 1.1 × 10−4; 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.
Ijms 22 04258 g002
Ijms 22 04258 g003 550
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−16) 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−4, (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 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−16) 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−4, (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.
Ijms 22 04258 g003
Ijms 22 04258 g004 550
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).
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).
Ijms 22 04258 g004
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Borczyk, M.; Zieba, M.; Korostyński, M.; Piechota, M. Role of Non-Coding Regulatory Elements in the Control of GR-Dependent Gene Expression. Int. J. Mol. Sci. 2021, 22, 4258. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22084258

AMA Style

Borczyk M, Zieba M, Korostyński M, Piechota M. Role of Non-Coding Regulatory Elements in the Control of GR-Dependent Gene Expression. International Journal of Molecular Sciences. 2021; 22(8):4258. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22084258

Chicago/Turabian Style

Borczyk, Malgorzata, Mateusz Zieba, Michał Korostyński, and Marcin Piechota. 2021. "Role of Non-Coding Regulatory Elements in the Control of GR-Dependent Gene Expression" International Journal of Molecular Sciences 22, no. 8: 4258. https://0-doi-org.brum.beds.ac.uk/10.3390/ijms22084258

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop