Ep300 sequestration to functionally distinct glucocorticoid receptor binding loci underlie rapid gene activation and repression.

The rapid transcriptional response to the transcription factor, glucocorticoid receptor (GR), including gene activation or repression, is mediated by the spatial association of genes with multiple GR binding sites (GBSs) over large genomic distances. However, only a minority of the GBSs have independent GR-mediated activating capacity, and GBSs with independent repressive activity were rarely reported. To understand the positive and negative effects of GR we mapped the regulatory environment of its gene targets. We show that the chromatin interaction networks of GR-activated and repressed genes are spatially separated and vary in the features and configuration of their GBS and other non-GBS regulatory elements. The convergence of the KLF4 pathway in GR-activated domains and the STAT6 pathway in GR-repressed domains, impose opposite transcriptional effects to GR, independent of hormone application. Moreover, the ROR and Rev-erb transcription factors serve as positive and negative regulators, respectively, of GR-mediated gene activation. We found that the spatial crosstalk between GBSs and non-GBSs provides a physical platform for sequestering the Ep300 co-activator from non-GR regulatory loci in both GR-activated and -repressed gene compartments. While this allows rapid gene repression, Ep300 recruitment to GBSs is productive specifically in the activated compartments, thus providing the basis for gene induction.

INTRODUCTION

Cells continuously respond to hormones, such as glucocorticoids, by rapidly remodeling their transcriptional programs. Specifically, glucocorticoids signal through binding to the glucocorticoid receptor (GR, encoded by Nr3c1), a transcription factor of the nuclear receptor superfamily. The hormone-activated GR translocates to the nucleus, where it binds to chromatin target sites and directly reprograms the transcription of multiple genes. GR binding to DNA is highly cell-type-specific and relies on the interplay of GR with transcription factors (TFs) that define the locations of accessible chromatin (1–3). GR can bind these accessible sites together with other factors to induce variable transcription patterns (4).

Gene activation is thought to be primarily controlled by direct binding of GR to its recognition element (GRE) and subsequent recruitment of co-activators such as SRC-1, GRIP1 and CBP/EP300 (5). Transcriptional repression models, on the other hand, include both direct and indirect GR binding to DNA. However, genome-wide studies have shown that both indirect binding through tethering TFs, such as AP1 and NF-κB, and direct binding to ‘negative’ GR elements (nGRE) are equally common for up- and down-regulated genes, and thus cannot explain the opposite transcriptional responses (1,6–8).

Another mechanism suggested for transcriptional repression by GR involves competition with other transcription factors for a limited amount of co-activators (9). This squelching model is supported by the fact that at the genomic level, co-activators are recruited to binding loci of GR as well as several other nuclear receptors, while depleted from other regulatory elements (10–12). However, since co-activator redistribution is global, i.e. occurring near both repressed and activated genes, its depletion from regulatory loci near repressed genes alone, as suggested for the estrogen receptor (13), cannot explain gene repression, without integrating the dynamic events at other associated regulatory loci. Another study indicated an increase of Ep300 signal at contact domains associated with GR-upregulated and a decrease at contact domains of GR-downregulated genes (14). However, the contact domains included both GR binding sites (GBSs) and non-GBSs, thus not distinguishing the dynamics of Ep300 signal at each type of regulatory element and not reflecting the possible exchange of Ep300 between them. Also, in reporter assays, variations of CBP/EP300 cellular levels demonstrated an inconsistent effect on GR-mediated gene repression (15–17). This may indicate that the complexity of multiple regulatory elements in the chromatin context in reporter assays has been overlooked. Nevertheless, a formal evidence for Ep300 being a limiting factor is required to support that its sequestering drives rapid gene repression.

Mouse (and human) chromosomes are partitioned into topologically associating domains (TADs), large chromosomal units (ranging from tens of kb to 3 Mb) of high spatial connectivity that encompass multiple genes and regulatory elements (18,19). Chromosomal contacts between GBSs and their target genes within TADs show cell-type-specificity (20). In addition, TADs are exclusive to either gene activation or repression by progesterone and estrogen receptors (21), indicating that the spatial organization of genomic information is an important component of gene regulation. The majority (∼90%) of GBSs are distant from the transcription start sites (TSSs) of their nearest gene, and GR binding events greatly exceed the number of GR-responsive genes (2,22), suggesting that multiple GBSs and genes are associated by chromatin folding. Thus multiple non-GR regulatory loci (non-GBSs) residing in these TADs may affect the transcriptional response to glucocorticoids. Hence, the molecular makeup and the residing transcription factors of GBS and non-GBS loci may define the specific transcriptional response to GR, and must be studied in their entirety.

This study aimed to explore the mechanism of GR-mediated regulation by dissecting the relevant regulatory components and understanding their spatial crosstalk. To this end, we applied high-throughput and high-resolution gene-centric circular chromosome conformation capture (4C-seq) to determine the scope of regulatory sites associated with GR-responsive genes via chromosomal looping. GR-responsive genes clustered in 3D with multiple GBSs and other regulatory elements exclusive to either transcriptional activation or repression by GR. To explore features of the bi-polar transcriptional responses to GR we compiled catalogs of regulatory elements associated with GR-induced or -repressed genes and combined these data with multiple genomic datasets in resting and hormone-treated cells. Our results revealed distinct transcription factor composition and chromatin dynamics for GBS and non-GBS regulatory elements associated with GR-activated or repressed genes. The TFs include KLF4 and STAT6, which repress GR up-and down-regulated genes, respectively, largely independently of GR-mediated regulation. In addition, we identified the transcription factors Rev-erb and ROR as modular co-regulators of GR-activated genes. Notably, following dexamethasone treatment, the co-activator Ep300 was rapidly sequestered from non-GBS loci, associated with GR-repressed genes, but unexpectedly also from non-GBSs associated with GR-activated genes. However, the H3K27ac signal, deposited by Ep300, was decreased specifically at the repressing non-GBS, while increased specifically at activating GBSs. Accordingly, sequestering Ep300 from the ensemble of non-GBS’s together with its non-productive recruitment to associated GBS loci is a major mechanism for rapid GR-mediated transcriptional repression. Thus, the spatial convergence of multiple GBS and non-GBS elements that vary between GR activated and repressed genes underlies the mechanistic complexity of the bipolar transcriptional response to hormone-activated GR.

MATERIALS AND METHODS

Cell growth

Mouse mammary epithelial adenocarcinoma (3134) (23) and human pulmonary adenocarcinoma (A549) cells were maintained in Dulbecco's modified Eagle's medium (DMEM, Biological Industries) supplemented with 10% fetal bovine serum (FBS, Biological Industries), 2 mM l-glutamine, 0.5 mg/ml penicillin–streptomycin, and 1 mM sodium pyruvate (Biological Industries) in a humidified 37°C incubator with 5% CO2. Before hormone treatment (100 nM dexamethasone (Dex; Sigma) for 1 h, or carrier [EtOH] control), cells were cultured in media containing 10% charcoal-stripped serum (CSS) (Sigma, USA). HepG2 cells (ATCC HB-8065) were cultured in DMEM media supplemented with 10% FBS, 1% penicillin–streptomycin and 1% l-glutamine.

RNA extraction and quantification by qPCR

RNA was extracted and treated on-column with DNase I with Quick-RNA Miniprep Kit (Zymo Research). Total RNA concentration and purity were measured by Nanodrop. RNA integrity was verified on 1% agarose gel and by RNA ScreenTape Assay (Agilent). cDNA was produced with High Capacity cDNA Reverse Transcriptase Kit (Applied Biosystems) or with qScript cDNA Synthesis Kit (Quanta Bioscience) following the manufacturer's instructions. cDNA was quantified on a real-time PCR detection system (CFX Connect; Bio-Rad) using iTaq Universal SYBR Green Supermix (Bio-Rad). Cycling conditions consisted of initial denaturation (3 min) at 95°C, followed by 40 cycles of 10 s denaturation at 95°C and annealing/extension for 30 s at 60°C. Primer sequences are provided in Supplementary Table S1. qPCR reactions were initiated in a final volume of 10 μl, containing 300 nM each of forward and reverse primers and 25 ng of cDNA. Transcription levels were normalized to beta-actin and elongation factor 2 (eEF2) as internal controls. Experiments were repeated three times and are presented as fold change relative to the control. Calculations were performed using CFX Manager software (Bio-Rad).

Library preparation for RNA-seq

Messenger RNA (mRNA) was enriched from 1 μg of total RNA by Poly(A) mRNA Magnetic Isolation Module (New England Biolabs) according to the manufacturer's instructions. cDNA libraries were constructed using the NEBNext Ultra RNA Library Prep Kit (New England Biolabs) following the manufacturer's protocol. Library concentration was measured by DNA High Sensitivity Kit (Invitrogen) on a Qubit Fluorometer (Invitrogen). The High Sensitivity D1000 ScreenTape Assay combined with the 2200 TapeStation System (Agilent) was used to assess the quality of the libraries.

Silencing

For gene silencing by siRNA, 3134 cells were cultured for 24 h in Opti-MEM I Reduced Serum Media (Gibco) and transfected with 50 nM siRNAs using Lipofectamine 3000 (Invitrogen) according to the manufacturer's instructions. siRNA sequences are listed in Supplementary Table S1. DNA transfection was performed using jetPRIME (Polyplus, France). Following transfection, cells were cultured for 24 h in DMEM supplemented with 10% FBS, 2 mM l-glutamine, 0.5 mg/ml penicillin–streptomycin and 1 mM sodium pyruvate, and for an additional 24 h in CSS-based media. After 48 h post-transfection, cells were treated with 100 nM Dex or EtOH for 1 h, and RNA was collected.

Klf4 chromatin immunoprecipitation (ChIP) assay and qPCR

3134 cells were cultured in CSS for 24 h and treated with either vehicle (EtOH) or 100 nM Dex for 2 h. Cells were crosslinked for 10 min at 37°C in 1% formaldehyde (Sigma), followed by quenching with 125 mM glycine for 10 min. Crosslinked cells were re-suspended in RIPA buffer (10 mM Tris pH 8, 1 mM EDTA pH 8, 140 mM NaCl, 0.2% SDS, 0.1% DOC), supplemented with protease inhibitors and sonicated for 38 cycles of 30 s ON and 30 s OFF (Bioruptor sonicator, Diagenode). Cleared chromatin was incubated overnight with 10 μg α-KLF4 (sc-20691). Complexes were washed twice with RIPA buffer, twice with high-salt buffer (10 mM Tris pH 8, 1 mM EDTA pH 8, 500 mM NaCl, 1% Triton, 0.2% SDS, 0.1% DOC), twice with LiCl wash buffer (10 mM Tris pH 8, 1 mM EDTA pH 8, 0.25 M LiCl, 0.5% NP-40, 0.5% DOC), and once in TE buffer (10 mM Tris, 1 mM EDTA, pH8). Crosslinks were reversed with 0.2 mg/ml Proteinase K overnight at 65°C. Purified DNA served as a template for qPCR. Primers used for qPCR amplification are listed in Supplementary Table S1. Fold enrichment was calculated using CFX Manager software (normalizing for the input, and relative to negative control primers).

Ep300 chromatin immunoprecipitation (ChIP) assay and sequencing

3134 cells were cultured in CSS for 24 h and treated with either vehicle (EtOH) or 100 nM Dex for 1 h. Samples were sent to Active Motif (Carlsbad, CA) for ChIP-Seq. Active Motif prepared chromatin, performed ChIP reactions, generated libraries, sequenced the libraries, and performed basic data analysis. In brief, chromatin was isolated by adding lysis buffer, followed by disruption with a Dounce homogenizer. Lysates were sonicated, and the DNA was sheared to an average length of 300–500 bp with Active Motif's EpiShear probe sonicator (cat# 53051). Genomic DNA (Input) was prepared by treating aliquots of chromatin with RNase, proteinase K, and heat for de-crosslinking, followed by SPRI bead clean up (Beckman Coulter) and quantitation by Clariostar (BMG Labtech). An aliquot of chromatin (30 ug) was precleared with protein A agarose beads (Invitrogen). Genomic DNA regions of interest were isolated using 20 ug of antibody against p300 (Santa Cruz Cat# sc-585, Lot #: B0711). Complexes were washed, eluted from the beads with SDS buffer, and subjected to RNase and proteinase K treatment. Crosslinks were reversed by incubation overnight at 65°C, and ChIP DNA was purified by phenol–chloroform extraction and ethanol precipitation.

ChIP sequencing (Illumina)

Illumina sequencing libraries were prepared from the ChIP and input DNAs on an automated system (Apollo 342, Wafergen Biosystems/Takara). After a final PCR amplification step, the resulting DNA libraries were quantified and sequenced on Illumina's NextSeq 500 (75 nt reads, single end). Reads were aligned to the mouse genome (mm9) using the BWA algorithm (default settings). Duplicate reads were removed, and only uniquely mapped reads (mapping quality ≥ 25) were used for further analysis. Alignments were extended in silico at their 3’-ends to a length of 200 bp, the average genomic fragment length in the size-selected library, and then assigned to 32-nt bins along the genome. The resulting histograms (genomic ‘signal maps’) were stored in bigWig files. Peak locations were determined using the MACS algorithm (v2.1.0) with a cutoff of p-value = 1e−7. Peaks on the ENCODE blacklist of known false ChIP-Seq peaks were removed.

Ep300 overexpression

Ep300 sequence was amplified by PCR from pcDNA3.1-p300 vector (Addgene 23252) and cloned in frame to EGFP located in a pFUGW-H1 based plasmid. A549 cells were transfected with 1 microgram a Ep300-GFP or GFP vector using JetPrime reagent (Polyplus-transfection company) and cultured in media supplemented with 10% CSS for 48h. Following 1h incubation with 100nM Dex or ethanol, cells were fixed in 75% ethanol and re-suspended in 5 ml PBS with 0.5 U/μl RNase inhibitor (Takara, #2313B). GFP expressing cells were sorted by flow cytometry (BDFACSAria III, BD Bioscience). RNA was extracted from sorted cells using E.Z.N.A. Total RNA Kit (Omega Bio-tek, #6834-01). For qPCR analysis, RNA was reverse transcribed using qScript cDNA Synthesis Kit (Quanta Bioscience, #95047-100), or processed by - SMARTer Stranded Total RNA-Seq Kit v2 (Takara, 634411) for Illumina sequencing.

4C-seq

4C was performed as described (24,25). Proximity ligation junctions, reflecting in vivo spatial proximity, were generated with DpnII (New England Biolabs), followed by circularization with Csp6I (Thermo Scientific). Chromosomal contacts with the TSS-containing fragment were amplified with inverse PCR primers (Supplementary Table S1), and sequenced on the Illumina 2000 platform.

Computational analysis

4C-seq analysis

4C sequenced reads were sorted into FASTQ files for each viewpoint and condition according to the bait sequences using a custom Perl script. Data were analyzed using the 4Cseqpipe program (26). Briefly, the algorithm calculates medians of normalized coverage for running windows of linearly increasing size (2–49 kb, 1 kb steps), within an X kb window surrounding the bait. For each (2–49) window, the procedure generated a normalized signal for each 1 kb across X kb, giving rise to a table with 48 columns and X rows. Cell i,j in the table represents the normalized signal in the ith range (2–49 kb) in the jth position on the chromosome (from the X kb surrounding the bait; each position is 1 kb). These windows are presented in the figures as color-coded multiscale diagrams. In addition, a high-resolution contact intensity trend line, depicting the medians of 5 kb running windows, is presented together with the 20th and 80th percentiles. To define contact domains, we first assigned each 1 kb j position the maximum ith range (2–49 kb) with a normalized signal higher than the top quartile of all normalized signals. Contact domains were obtained by merging j positions assigned with the top ith range (49 kb), so that their normalized signal is higher across all ranges and closer than 15 kb.

Analysis of ChIP-seq and DNase I-seq data

GR, AP-1, RNA Pol-II, Ep300 ChIP-seq and DNase I-seq data from 3134 cells treated with Dex for 1 h from (1,2,27) (SRP004871, SRP007111, GSE61236) were analyzed as follows. Sequenced reads were aligned to the mouse genome (mm9) using BOWTIE (28), and peaks were called using MACS2 (29) with default parameters for ChIP, and - -nomodel; - -shift -100; - -extsize 200 parameters for DNase I hypersensitive sites (DHS). A p-value cutoff of 1e−10 was used for AP1, CTCF and Ep300; 1e−4 for GR; and 1e−3 for DHS.

Pol-II ChIP-seq analysis

A Pol-II density signal was assigned for each gene. Genes were considered to have detectable expression if the normalized density signal was >0.5 (median expression) under either of the conditions (10,565 expressed genes). The ratio (log2 fold change) between vehicle and Dex-treated cells was calculated, and z-score normalization was applied to the log2FC values. Genes were scored as GR-regulated if the absolute value of Dex-dependent z-score was >2 standard deviations.

RNA seq analysis

Reads were aligned to the mouse genome using STAR (30) and were counted on genes using HTseq (31). For each factor that was silenced, four combinations of enrichment analysis were analyzed: si-Neg.C, with and without Dex; target genes siRNA with and without Dex; si-Neg.C and gene-specific siRNA without Dex; and si-Neg.C and specific siRNA with Dex. Enrichment analysis between each combination was done using edgeR, which assigns an FDR value to each gene and calculates the log FC between the two conditions (32). Genes with FDR <0.05 and log FC >|0.5| were considered differentially expressed.

Motif analysis

Motif discovery analysis was performed by the HOMER suite, using the ‘findMotifsGenome.pl’ program (33) with the parameters ‘-size given -N 20000 -cpg - -noweight’ for both DHS and GBS. The regions of local maxima of read count within DHS peaks (sub-peaks) were calculated using a custom R script. 100–250 bp regions with the highest read count in the DHS peaks and 200 bp centered at the GR peak summits were submitted to HOMER. Motif analysis was done for DHS or GR peaks within the domains of GR up-regulated genes (up domains) or GR down-regulated genes (down domains) with the whole genome as background (HOMER default). The relative motif enrichment in up domains or down domains was calculated as the ratio between the proportions of the motifs in each group. The most highly enriched motifs are presented in Figure 2.

Figure 2.

Transcription factor configuration at regulatory elements associated with GR-activated genes. (A) Motifs enriched in DHSs and GBSs associated with GR-activated or repressed genes. (B) Recognition motifs of Klf4 at DHSs and GBSs (GR ChIP-seq) (data from hormone-treated cells (100 nM Dex, 1 h)), within 4C domain of Tsc22d3 (depicted in red). Genomic mm9 coordinates. Vertical red line indicates locus validated by KLF4 ChIP. (C) KlLF4 binding by ChIP-qPCR to regulatory elements with the KLF4 motif (indicated by the vertical red line in Figure 2B, Supplementary Figure S3); n = negative control loci. (D) Transcriptional response to GR activation (1 h Dex) measured by RT-qPCR in 3134 cells transfected with scrambled siRNA negative control (si-Neg.C) or Klf4 siRNA. Error bars indicate the SD of three biological repeats. * P < 0.05, **P < 0.01; Student's t-test. (E) Log2 of the fold change of genes activated by Dex in 3134 cells transfected withscrambled siRNA negative control (si-Neg.C) or Klf4 siRNA. (F) Transcriptional changes (log2 fold change Klf4 siRNA negative control siRNA) of GR upregulated genes (GR-up) in cells treated with vehicle (−Dex) or Dex. The genome-wide (GW) changes are shown as control. p-values are indicated; Wilcoxon test. (G) A model for gene regulation by GR and Klf4. Klf4 binding to GRE and non-GRE elements within a defined spatial domain reduce gene activation by GR.

Enrichment and statistical analysis

The relative enrichment of GR-responsive, up- and down-regulated genes on specific chromosomes was calculated as the ratio between the fraction of GR-regulated genes (from all GR regulated genes) on the specific chromosome and the percentage of genes (from the entire genome) on that chromosome. The significance of the relative enrichment/depletion of genes was calculated using upper/lower tail proportional test, using ‘prop.test’ functions in R.

The closest GBS to each TSS and its distance were calculated using Bedtools ‘closestBed’ script (34). Overlap between ChIP peaks and/or DHS and/or domains was identified using Bedtools ‘intersectbed’ script. Normalized ChIP signal was calculated as reads per million reads (from the whole dataset) per kb. Calculation of the overall GR ChIP signal in the 4C domains and flanking regions was performed by summing the ChIP signal in the 50 kb flanking both ends of all the 4C domains. Mann–Whitney–Wilcoxon Test statistics were applied to compare the differences of each group in up or down domains versus the entire genome using R ‘Wilcox.test’ function. The overlap between the GBS and GRE motifs (HOMER motif bed file with all putative GREs in the genome) or nGREs (CTCCGGAGA, CTCCNGGAGA and CTCCNNGGAGA) was calculated using Bedtools ‘intersectbed’ script.

RESULTS

GR-responsive genes localize within spatial domains specialized for transcriptional induction or repression by hormone-activated GR

To study the direct effects of GR on gene transcription, we screened for variation of activated (S5 phosphorylated) RNA polymerase-II occupancy 1 h after hormone induction, using available ChIP-seq profiles in mouse mammary 3134 cells (27). We retrieved 263 up-regulated genes (>two-fold Pol-II occupancy along gene body regions, z-score > 2) and 251 down-regulated genes (<2-fold, z-score ←2). The transcriptional response was validated by qRT-PCR (Supplementary Figure S1).

Next, we applied 4C-seq to define long-range associations of transcription start sites (TSS) of GR-responsive genes. 4C-seq measurements showed that GR-responsive genes are embedded in domains of high spatial connectivity, which drop sharply at the edges of the domain. The 4C domains indicate TAD structure and nested neighborhoods or loops (35,36). For example, the Zfp36I promoter loops across two Hi-C sub-domains to a GBS located ∼700 kb downstream (Figure 1A). 4C-seq measurements for TSS of GR-responsive genes (11 GR-activated and 7 GR-repressed) in vehicle and hormone-treated cells indicated that the borders of the domains remained stable after Dex treatment (Supplementary Figure S2). To attain high genomic coverage, we measured by 4C-seq the chromosomal contacts of 25 GR-activated and 15 GR-repressed genes, with a wide transcription level range spreading across 13 mouse chromosomes (Supplementary Table S2). We identified 12 additional GR activated genes within the previously identified GR activation domains (out of 167 genes in the domains; P < 10E−6, multinomial distribution), and 3 additional GR repressed genes positioned within the domains of GR-repressed genes (out of 57 genes; P < 0.01, multinomial distribution). For example, the Cxcl1 and Cxcl5 genes, both repressed by GR, share the same domain (Supplementary Figure S2). We found preferential co-localization of GR activated or repressed genes in TADs also in human A549 cells treated with Dex for 1 h (14). The ratio between additional GR up and down-regulated genes in 257 TADs with at least one GR-upregulated gene was 2.87 fold higher than the GW proportion (27 up/7 down, P < 0.001, multinomial distribution). The proportion of additional down vs up regulated genes in 183 TADs with at least one GR-downregulated gene was 1.6-fold higher than the GW proportion (12 down/10 up, P < 0.02, multinomial distribution). This co-localization of genes with similar responses within TADs is greater than expected by chance, indicating that GR-responsive genes are embedded within spatial domains specialized for transcriptional induction or repression by hormone-activated GR.

Figure 1.

Features of spatial domains of GR-regulated genes. (A) High-resolution chromosomal contact (4C-seq signal) profile of the GR-repressed gene Zfp36I (depicted in blue) in 3134 cells. 4C contact profiles with the viewpoint (highlighted by a vertical dashed line) are shown in a trend line (upper panel, black line, using 5 kb sliding window), and in a color-coded scale domainogram, which shows relative interactions (red indicates the strongest interactions) in a sliding window ranging from 2 to 50 kb. Hi-C domains from murine CH12-LX cells (71) are shown. GBS-GR binding sites were obtained from ChIP analysis. (B) Boxplots showing the distance between TSS and the nearest GBS of GR up-regulated (red), down-regulated (blue), and genome-wide (GW) transcribed genes (white). *P = 3.2e−34, ** P = 2.9e−76; Wilcoxon test. (C) GBS density (peaks per kb) in 4C domains of GR up-regulated (red), down-regulated (blue), and GW (random 200 kb domains, white). *P = 0.006, **P = 1e−9; Wilcoxon test. (D) GR ChIP signal at GBS associated with up-, down-regulated genes, and genome-wide (white). *P = 5.7e−6, **P = 8.6e−19; Wilcoxon test. (E) Profiles of aligned 4C domain border regions (left and right borders combined) are shown for GR binding (sum of ChIP signal in 50 kb from border). Red/blue bar and positive coordinates, inside 4C domains; white bar and negative coordinates, outside 4C domains. (F) GRE motif in GBS and DHS associated with promoters of up-and down-regulated genes. (G) Chromatin accessibility (DHS) signal in GBS and non-GBS regulatory sites within domains of GR up-regulated (red), down-regulated (blue) and genome (white), after 1 h Dex (+) or vehicle (−) treatment. *P < 0.005; Wilcoxon test.

GR-activated and repressed domains differ in chromatin structure and dynamics

Within the 4C domains, GR-responsive genes associate in 3D with multiple GBS and other non-GBS regulatory elements identified by their accessibility to DNase I digestion and Ep300 binding (Supplementary Figure S2). Since the transcriptional response within a given spatial domain is specific to either activation or repression by GR, the regulatory elements within the domain may serve as a toolbox to elicit specific transcriptional responses to GR signaling. To identify features discriminating between activating and repressing GBS’s and DNase I hypersensitive sites (DHS’s), we used the 4C domains and pooled catalogs of 231 GBSs and 910 DHSs associated with GR-activated genes, and 76 GBSs and 328 DHSs associated with GR-repressed genes. A marked difference was shown in GR distribution within domains of up-and down-regulated genes in mouse mammary cells: GBSs were closer to promoters of GR-activated genes compared to repressed genes (Median distance of 4, 38 and 51 kb from activated-, repressed- and genome-wide expressed genes, respectively. Figure 1B), as also shown in human cells (14,22,37,38).

Considering all GBSs within the defined 4C chromatin domain of the regulated gene, we found that both GBS density and ChIP signal were higher in domains of GR-activated relative to GR-repressed genes (Figure 1C, D). Notably, the GR signal was higher within domains of both classes compared to their flanking genomic regions (Figure 1E), indicating that the enrichment of GR binding within domains of GR-responsive genes is specific. Moreover, the higher GR ChIP signal corresponds to a higher proportion of GR binding sites with a canonical GR recognition motif (GRE) in domains of GR-activated genes (Figure 1F), suggesting that indirect GR binding is more pronounced at GBSs associated with GR-repressed genes. Regulatory elements of the two groups also differed in chromatin structure and dynamics. Following Dex treatment chromatin accessibility was increased in GBSs and other non-GBS regulatory elements associated with GR-activated genes and remained stable at these loci associated with GR-repressed genes (Figure 1G), indicating that hormone treatment enhances chromatin accessibility, especially at GR sites associated with activated genes.

Regulatory loci in GR-activated and repressed domains differ in their transcription factor composition

The distinct features of regulatory sites associated in 3D with GR-activated versus GR-repressed genes suggest that their molecular makeup is different. To identify possible TFs that define regulatory elements associated with GR-activated or repressed genes, we performed a motif discovery analysis in GBSs and DHSs using the genome as background. For motifs with a p-value lower than 0.01, we calculated the ratio of the motif in the GR-activated relative to the GR-repressed domains and vice versa. Analysis of both GBSs and DHSs revealed motifs of several candidate factors (Figure 2A), some of which were previously linked to GR biology in other cell types, though less is known about their co-regulatory activity at the chromatin level.

The KLF4 binding motif was enriched in DHSs associated with promoters of GR-induced genes (Figure 2A, B, Supplementary Figure S3A). KLF4 binding was confirmed by ChIP-qPCR in the domains of the Tsc22d3, Tgm2, Bcl2l1 and Pkp1 genes (Figure 2C, Supplementary Figure S3A). Transcriptional activation of these genes was enhanced in cells with siRNA-mediated knock down of KLF4 expression (Figure 2D, Supplementary Figure S3B–D). On a global scale, the major repressive effect of KLF4 on GR-activated genes was independent of Dex application (Figure 2E–G, Supplementary Figure S3E).

The binding motif of REV-ERB (RevDR2) was highly enriched in GBSs and DHSs associated with activated genes (Figure 2A). This motif is also the recognition sequence of RORs (RAR-related orphan receptors), termed ROR response element (RORE). Rev-erb and ROR belong to the nuclear receptor superfamily of transcription factors, with antagonistic activities: REV-ERB represses transcription by recruiting the NCoR/SMRT co-repressor, while RORs recruit the SRC and/or EP300 co-activators to ROREs (39). From the members of the Rev-erb and ROR families, we focused on RORβ, RORɣ and Rev-erbβ (Nr1f2, Nr1Ff3 and Nr1d2, respectively), which were the abundant transcripts in 3134 cells (Supplementary Figure S4A). To study the contribution of ROR and Rev-erb factors to gene regulation by GR, we down-regulated their expression by siRNA (Supplementary Figure S4B–G). Similar to Rev-erbα in liver cells (40–43), knocking down Rev-erbβ in 3134 cells up-regulated Bmal1 expression (Supplementary Figure S4G). In line with its expected repressive activity, knocking down Rev-erbβ led to higher Dex-mediated activation of Nuclear factor IL-3 (Nfil3) and Perilipin4 (Plin4), which are associated in 3D with regulatory loci having RORE motifs (Figure 3A, B). An extension of this analysis to the whole genome revealed similar trends (Figure 3C, Supplementary Figure S4I). In line with their predicted activating role, knocking down RORβ and RORɣ reduced the transcriptional induction of Nfil3 and Plin4 (Figure 3E). Genome-wide RNA-seq analysis showed that transcriptional activation by GR is indeed significantly reduced under RORβ and RORɣ silencing (Figure 3F, Supplementary Figure S4H), suggesting that RORβ and RORɣ collaborate with GR in the up-regulation of these genes. Rev-erbβ, RORβ and RORɣ knock-downs did not affect the basal expression of GR-activated genes prior to Dex application (Figure 3D, G, Supplementary Figure S4H, I), indicating that these factors regulate the rapid GR-mediated transcriptional response. Overall, our results suggest that Rev-erb and ROR factors are transcriptional co-regulators of GR-mediated gene activation (Figure 3H).

Figure 3.

Cooperation between ROR, Rev-erb and GR transcription factors. (A) Recognition motifs of ROR and Rev-erb at DHSs and GBSs (ChIP-seq) (data from hormone-treated cells (100 nM Dex, 1 h)) within 4C domains of Nfil3 and Plin4 (depicted in red). Genomic mm9 coordinates. (B, E) Transcriptional response to GR activation (1 h Dex) in 3134 cells transfected with scrambled siRNA (si-Neg.C) or Rev-erbβ / RORβ and RORɣ siRNAs. Error bars indicate the SD of three biological repeats. * P < 0.05, **P < 0.01; Student's t-test. (C, F) Log2 of the fold change of genes activated by Dex in 3134 cells transfected with scrambled siRNA (si-Neg.C) or Rev-erbβ / RORβ and RORɣ siRNAs. (D, G) Transcriptional changes (log2 fold change of specific knock-down /negative control siRNA (Neg.C)) of GR upregulated genes (GR-up) in Rev-erbβ (D), and RORβ and RORɣ (G) siRNA transfected cells treated with vehicle (−Dex) or Dex. p-values are indicated; Wilcoxon test. (H) A model for gene regulation by GR, ROR and Rev-erb. Both, Rev-erbβ and RORβ/ɣ, negative and positive GR co-regulators, can bind a composite GR binding site harbouring a GRE.

An analysis of motifs enriched in GBSs of the down- vs. up-regulated domains revealed nuclear factor kappa-B (NF-κB) and signal transducer and activator of transcription 6 (STAT6) (Figure 2A). Notably, the NF-κB motif was also enriched in non-GBS loci associated with GR-repressed genes (Figure 2A), suggesting that GR-mediated repression of NF-κB target genes extends its implicated activity either by binding directly or by tethering to GBS (4,5,44,45). Previous studies showed that GR physically interacts with members of the STAT family in both gene activation and repression (46,47). Notably, while the half-palindromes (TTC) in binding motifs of STAT factors are separated by three nucleotides (N3), STAT6 preferentially binds N4 sites (48), suggesting that its role is not redundant with the other STATs. Downregulation of STAT6 by siRNA led to elevated expression of Cxcl5, Cxcl1, Ptgs2 and Zfp36l1 prior to hormone application but did not affect the suppressive effect of GR (Supplementary Figure S5). Thus, in line with a previous report in HeLa cells (49), GR and STAT6 act independently as negative regulators of these genes. In agreement with previous studies (6,7,14), we found that the nGRE motif does not discriminate between activating and repressing GBSs since it is present at similar frequencies in both groups of GR sites (0.12% in up-regulated, 0.08% in down-regulated, and 0.08% genome-wide).

Collectively, GBS and other regulatory loci associated with GR-activated or repressed genes vary their TFs assignment. The KLF4 and STAT6 pathways that converge with GR in activated or repressed domains, respectively, contrast the transcriptional effects of GR independently of GR activation. In addition, Rev-erbβ, RORβ and RORɣ regulate GR-mediated gene activation.

GR binding sequesters Ep300 from active enhancers

To examine whether the bi-polar transcriptional responses to GR are related to variations in enhancer activity, we analyzed the loading of Ep300, which marks active enhancers (27,50). Analysis of all the regulatory sites showed that the Ep300 ChIP-seq signal increased at enhancers associated with GR-activated genes but was reduced at regulatory sites associated with repressed genes (Figure 4A red and blue boxes, respectively). Similar trends were observed for loci associated with GR responsive genes in human A549 lung cells (14). However, sorting the accessible loci by their GR occupancy revealed that the increased Ep300 signal in the activated domains was derived from the subset of GR-binding loci, and decreased at regulatory elements that are not bound by GR (non-GBS) (Figure 4B red boxes, Figure 4DPlin4). The Ep300 signal was decreased also from non-GBSs associated with GR-repressed genes, but did not vary at their associated GBSs (Figure 4B blue boxes, 4D Ccl2), indicating that the reduction of the Ep300 signal in the GR-repressed domains derived from the non-GBSs. Since Ep300 acetylates histone 3 at lysine 27 (H3K27ac) at regulatory elements, we analyzed the dynamics of H3K27ac in human A549 cells treated with Dex for 1 hour (51,52). Remarkably, the predominant effect in GR-activated domains was an increase in H3K27ac signal at GBSs in GR-activated domains, with a decrease in H3K27ac signal at accessible loci that are not bound by GR (non-GBS) in GR-repressed domains (Figure 4C). Thus, following GR activation, the Ep300 ChIP signal and its enzymatic activity were depleted from non-GBS loci associated with down-regulated genes, while recruited to GBSs associated with up-regulated genes (Figure 4E).

Figure 4.

Features of regulatory elements associated with GR-regulated genes. (A) Ep300 signal at DHS peaks within 4C domains. P-Wilcoxon test. (B) A similar analysis for Ep300 loci at GBS’s and accessible sites where GR does not bind (non-GBS). p-values of Wilcoxon test are indicated. (C) H3K27ac signal in GBS or non-GBS sites within domains of GR up-regulated (red), down-regulated (blue), and genome (white) after 1 h Dex (D) or EtOH (E) treatment in A549 cells (data from (51,52)). (D) Binding profiles of Dex-activated GR and Ep300 in Vehicle (EtOH) and Dex treated cells. Plin4 is GR-upregulated, and Ccl2 is down-regulated by GR. GBS and non-GBS are highlighted in grey and white rectangles, respectively. (E) A model for gene regulation by cross-talk between regulatory elements within the defined spatial domain. The activated domains contain more GBSs closer to the TSS relative to the repressed domains. In gene activation, GR binding to DNA increases the Ep300 signal at GBSs. Ep300 is also sequestered from weak enhancers (orange). In gene repression, GR is tethered by a transcription factor such as AP-1 (green) and sequesters Ep300 from strong, active enhancers (blue) within the defined spatial domain. (L, M) 3134 cells transfected with either siRNA for Ep300 knock-down or scrambled (si-Neg.C.) siRNA as control, and treated with Dex (+) or EtOH (−) or for 1 h. Log2 fold change of transcript levels in response to Dex (F) and siEp300 (G) are shown. p-values are indicated; Wilcoxon test. (H, I) A similar analysis for overexpression (OE) of Ep300 or GFP control in A549 cells.

Gene down-regulation by sequestering Ep300 from non-GBSs requires that Ep300 be a limiting factor. Knocking down Ep300 by siRNA indeed attenuated GR-mediated gene repression (Figure 4F, blue boxes), primarily by reducing the expression of its target genes prior to hormone application (Figure 4G, blue boxes, E). However, Ep300 knockdown had only a minor effect on GR-mediated gene induction (Figure 4F, red boxes). Increasing the amount of Ep300, on the other hand, was expected to negate GR-mediated gene repression by providing a sufficient amount of Ep300 to maintain high gene expression. Overexpressing the large (264 kDa) Ep300 gene in 3134 cells was impractical due to low transformation efficiency and poor stability of expression. Introducing Ep300 fused to GFP to A549 cells was more efficient (detected in 5–20% of the cells) and allowed enrichment of the transformed cells by FACS. Nevertheless, the instability of Ep300 expression required shortening the time-frame of the experiment. Therefore, infected cells were treated with Dex for 1 h and then fixed with ethanol and sorted by FACS based on GFP expression. RNA was extracted from the sorted cells and sequenced. RNA-seq analysis indicated that elevating Ep300 levels attenuated both GR-mediated up-and down-regulation (Figure 4H, blue and red boxes, respectively). In contrast to Ep300 depletion, which reduced the expression of GR-downregulated genes prior to Dex treatment, elevating Ep300 levels diminished the suppressive effects of GR by maintaining high gene expression following GR induction. Under Ep300 overexpression, gene activation by GR was also impaired. Ep300 overexpression elevated the expression of GR up-regulated genes prior to Dex treatment (Figure 4I), obviating the effect of GR in gene activation. Notably, GR-responsive genes were more sensitive to alterations in Ep300 levels than genes that are expressed (RPKM > 1) or expressed at the average level of GR responsive genes (RPKM∼300) (Figure 4F–l gray and white boxes), indicating some specificity between GR-responsive genes and the Ep300 co-activator.

Altogether, in line with Ep300 as a limiting factor, its reduction affected primarily GR-mediated gene repression rather than gene-activation. On the other hand, supplementing the cells with Ep300 probably countered the effect of its sequestering from non-GBSs, and thereby abolished GR suppressive effects of GR. These results suggest that gene repression by GR is predominantly passive and is achieved through sequestering Ep300 from other non-GR binding sites and enhancers (Figure 4E).

DISCUSSION

Despite our capacity to measure transcript levels with great accuracy and our ability to describe transcription factor binding loci and chromatin structure, the mechanism through which transcription factors rapidly activate some genes, while at the same time repressing others is unclear. In the current study, we addressed this question using the hormone-activated GR based on its rapid, specific, and direct transcriptional regulatory activity. Using 4C-seq to map regulatory elements exclusive to either GR-mediated gene activation or repression, we investigated how these elements dictate a specific transcriptional response to GR. To focus on the primary and direct activities of GR, we treated cells with Dex for 1 hour.

This study shows that sequestering Ep300 from active enhancers is a major and global mechanism for GR-mediated gene repression, as previously suggested (53,54). We provide evidence that Ep300 is a limiting factor by showing that the effects of Ep300 sequestering upon GR activation could be abrogated by overexpressing Ep300, while a reduced pool of available Ep300 resulted in repression of genes that are downregulated upon GR activation.

Understanding how previously reported genomic redistribution of co-activators is related to gene activation and repression required integrating the dynamics of the Ep300 signal at GBSs together with their 3D-associated non-GBS in the context of a specific transcriptional response. We found that Ep300 was depleted from non-GBS enhancers associated with both GR-repressed and activated genes, suggesting that this depletion alone is not sufficient to mediate gene repression (13). Importantly, a key discriminating feature of the different responses was the elevated Ep300 signal, specifically at GBSs associated with GR-activated genes. Thus, the opposite transcriptional responses can be explained by the sum of events at GBSs and non-GBSs that associate with GR responsive genes in 3D. In addition, prior to Dex application, the Ep300 ChIP signal was higher at non-GBS enhancers of the repressed genes, relative to the activated ones. Thus, the magnitude of Ep300 loss following Dex application was greater in non-GBS enhancers of the repressed genes. In the same line, super enhancers with high MED1 occupancy are more sensitive to repression following NF-κB activation (12). Importantly, the genomic effects of over-expressing and down-regulating Ep300 on the transcriptional response to GR stimulation support that Ep300 is present at a limiting amount. Recent studies showing that clustered enhancers can form phase-separated condensates (55,56) support the idea that the competition for Ep300 is limited to a local sub-nuclear volume represented by the structural domain (Figure 4E).

Integrating gene expression and regulatory information in 3D allowed us to focus on specific TFs involved in specific transcriptional responses to GR. We found that KLF4 and STAT6, which have been implicated in GR biology (49,57–59), regulate GR target genes largely independently of GR activity. Rev-erbβ and RORβ/ɣ were identified as negative and positive regulators, respectively, of gene activation by GR. Notably, Rev-erb and ROR are critical components of the circadian clock and regulate many physiological processes that are also regulated by glucocorticoids (39). Genomic binding of GR and REV-ERBα on regulatory elements is prevalent in the liver, and both factors associate in-vitro (60). Circadian (24 h) oscillations of ROR expression are synchronized with glucocorticoid release from the adrenal glands, which occurs 12 h out of phase with oscillatory REV-ERB expression. This finding reinforces the notion that RORs amplify GR-mediated gene activation concomitant with an increase in plasma glucocorticoids, while REV-ERB helps to reduce the expression of these genes alongside a decline in plasma glucocorticoids. This further suggests that the antagonistic functions of Rev-erb and ROR are physiologically relevant.

The different combinations of TFs in activating and repressing domains suggests that the exchange of Ep300 occurs between different types of regulatory elements in these two types of domains. For example, NF-κB and STAT6 are enriched in the repressed domains, while ROR and KLF4 are in the activated domains. In addition, GBS’s in the repressed domains are characterized by relatively low density, low frequency of the GRE motif, and low GR ChIP signal relative to GBSs in the activated domains, and are suggested to have a lesser effect on gene activation.

GR interacts dynamically with other transcription factors, co-activators and the chromatin template (61,62), and the GR ChIP signal reflects a cell population average of highly dynamic protein-DNA interactions (63). As Ep300 is recruited to chromatin through its interaction with the GR (45), the role of Ep300 in both the activating and repressing GBS may stem from different association modes between GR and chromatin. Whether GR binding at loci lacking GRE, which are more prevalent in repressed domains, is indirect via tethering to an intermediate DNA-bound transcription factor (6,64,65) or occurs directly to cryptic GREs (44), its residence time on chromatin is expected to be shorter in repressing vs. activating GBSs, and thus its associated Ep300, as well. Consequently, even if the affinity between GR and Ep300 is similar for activating and repressing GBS’s, on average, a shorter residence time of GR on GBS’s associated with repressed genes is likely to be less productive. In line with this notion, it was recently shown that the affinity of the GR to DNA is related to the resulting transcription activity (66) and that transient GR occupancy (on a time scale of 1h, measured by ChIP-seq) is linked to loss of Ep300 signal at GBSs (67). In addition, a higher GR-DNA exchange rate would provide more opportunities for GR to sequester Ep300 from proximal non-GBS regulatory loci, allowing efficient sequestering, while using fewer GBSs in the repressed compartments.

A recent study showed that although tethered GR binding loci are predominantly found near GR-repressed genes, direct GR binding to chromatin is required for GR-mediated gene repression (64). It is thus likely that the interplay between direct and indirect GR binding loci within spatial domains is required for sequestering Ep300 from active enhancers without activating GR-bound enhancers. Our data suggest that the spatial crosstalk between multiple regulatory elements is critical for GR-mediated gene regulation. TFs that participate in regulating the expression of GR target genes or the GR-mediated transcriptional response, can act through both GBSs and non-GBSs. For example, the GR-repressive effect could result from both sequestering of Ep300 from non-GBSs bound by NF-κB and interfering with GBSs bound by NF-κB.

Therefore, understanding the ensemble of features of regulatory elements and gene promoters in the context of the spatial domain can provide a quantitative understanding of GR-mediated gene regulation. This requires a systematic perturbation of different enhancer combinations in their endogenous 3D context. Genomic redistribution of co-activators has been reported in response to activation of different signaling pathways including hormones, inflammation and UV exposure (10–13,68–70), suggesting that co-factor sequestering is a general mechanism for rapid genetic reprogramming.

DATA AVAILABILITY

The data described in this publication have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO Series accession numbers GSE106209 and GSE193988.

Click here for additional data file.