Thyroid hormone dependent transcriptional programming by TRβ requires SWI/SNF chromatin remodelers.

Transcriptional regulation in response to thyroid hormone (3,5,3'-triiodo-l-thyronine, T3) is a dynamic and cell-type specific process that maintains cellular homeostasis and identity in all tissues. However, our understanding of the mechanisms of thyroid hormone receptor (TR) actions at the molecular level are actively being refined. We used an integrated genomics approach to profile and characterize the cistrome of TRβ, map changes in chromatin accessibility, and capture the transcriptomic changes in response to T3 in normal human thyroid cells. There are significant shifts in TRβ genomic occupancy in response to T3, which are associated with differential chromatin accessibility, and differential recruitment of SWI/SNF chromatin remodelers. We further demonstrate selective recruitment of BAF and PBAF SWI/SNF complexes to TRβ binding sites, revealing novel differential functions in regulating chromatin accessibility and gene expression. Our findings highlight three distinct modes of TRβ interaction with chromatin and coordination of coregulator activity.

INTRODUCTION

Thyroid hormones and their cognate thyroid hormone receptors (TRs) play key roles in maintaining cell identity, metabolism, and homeostasis in all tissues. Thyroid hormone receptors (TRs) function primarily as transcription factors that regulate broad networks of target genes in response to fluctuations in thyroid hormone (3,5,3′-triiodo-l-thyronine, T3) levels. Two distinct genes encoding TRα (THRA) and TRβ (THRB) are expressed in tissue-specific patterns. TRβ is predominant in the liver, kidney, and thyroid, while TRα is prevalent in the heart, bone, and brain (1). TRs are known master regulators of chromatin structure that dynamically regulate transcriptional activation and repression. The classic bimodal switch model (reviewed in (2,3)), in which TRs are constitutively bound to chromatin and T3 binding promotes dissociation of corepressors and recruitment of coactivators, was the predominant model used to describe transcriptional regulation by TRs for many years. However, more recent genome-wide studies of TRβ binding have shown significant T3-dependent recruitment of TRβ to chromatin and de novo chromatin remodeling, reflecting a more nuanced and complex model of transcriptional regulation (4–6). Moreover, a large proportion of TRβ binding occurs in distal regulatory elements and thereby directs histone acetylation of enhancer regions and higher order chromatin structure (7). Recruitment of cofactors by TRβ has been recently described as a T3-dependent coregulator shift, rather than a complete loss of corepressors and gain of co-activators upon ligand binding (8). Our data supports a multi-modal regulation model for TRβ interaction with chromatin that integrates each of these concepts.

Functional organization of the genome is dynamically regulated by numerous epigenetic mechanisms that ensure a rapid integration of different signals. ATP-dependent chromatin remodeling complexes are known to act as coregulators for nuclear receptors and have been implicated in chromatin remodeling specifically by TRβ. Members of the SWI/SNF complex were found to be associated with NCOR, a well-known co-repressor for TRs and other non-steroidal hormone receptors (9). This provided indirect evidence that SWI/SNF chromatin remodeling might be required for target gene repression by TRs. We described an interaction between the SWI/SNF core subunit BRG1 and TRβ for repression the oncogene RUNX2 (10). Recruitment of the SWI/SNF complex to T3-activated promoters has been suggested to be dependent upon interactions between SRC and p300, rather than a direct interaction between BRG1 and TR (11). Direct interaction between BAF57, a key SWI/SNF subunit, and TRβ have been observed at T3-responsive promoters (12). Given the limited studies of SWI/SNF participation in T3-regulated gene expression, a specific mechanism for TR interactions with SWI/SNF complexes has yet to be defined.

In this study, we mapped the binding sites of endogenous TRβ in normal human thyroid cells and integrated this data with changes in chromatin accessibility to classify three distinct modes of TRβ chromatin interaction and remodeling. Of note, a majority of TRβ binding sites were found in distal regulatory elements. In order to identify the protein-protein interactions involved in T3-dependent transcriptional regulation within these modes, we performed a proximity labeling assay followed by mass spectrometry. Some interactions between TRβ and its binding partners are gained or lost with T3, but notably the majority do not exhibit a T3-dependent switch. Many novel interactions were revealed in this screen, along with several well-characterized TRβ binding partners. Differential enrichment of interactions between TRβ and members of BAF and PBAF SWI/SNF complex subclasses were identified. We further demonstrated that TRβ differentially recruits SWI/SNF complexes to its binding sites. Based on our comprehensive genomic and proteomic analyses, we propose a new model whereby selective recruitment of BAF and PBAF SWI/SNF complexes to TRβ binding sites regulates chromatin accessibility and gene expression.

MATERIALS AND METHODS

Cell culture and hormone treatments

The immortalized normal human thyroid epithelial cell line Nthy-ORI (Sigma) was routinely cultured in RPMI 1640 growth media with l-glutamine (300 mg/l) (Sigma), sodium pyruvate and nonessential amino acids (1%) (Cellgro/Mediatech), supplemented with 10% fetal bovine serum (Peak Serum) and penicillin-streptomycin (200 IU/l) (Cellgro/Mediatech) at 37°C, 5% CO2, and 100% humidity. For T3 treatments, cells were hormone-starved for 24 h in growth media substituted with phenol-red free RPMI 1640 and charcoal-stripped fetal bovine serum (Sigma) prior to the addition of 10 nM T3 or the equivalent volume of 1 N NaOH vehicle for the indicated time course. Nthy-ORI cell line was authenticated by the Vermont Integrative Genomics Resource at the University of Vermont with short tandem repeat profiles using the Promega GenePrint10 platform.

CUT&RUN

Sample collection and library preparation

CUT&RUN was performed as previously described (13). Briefly, Nthy-ORI cells were harvested after 6 h of treatment with T3 or vehicle control, washed, and bound to activated Concanavalin A coated magnetic beads (Epicypher 21-1401). Cells were then permeabilized with Wash buffer (20 mM HEPES pH 7.5, 150 mM NaCl, 0.5 mM spermidine 0.05% digitonin). Permeabilized cells were then incubated with the indicated antibody (Supplemental Table S3) at 4°C with constant agitation overnight. Cells were washed twice more before incubation with recombinant p-AG MNase (Epicypher 15-1016) at 4°C for 2 h. Liberated DNA was purified, and libraries were prepared using the NEB Ultra FS II DNA Library Kit (NEB E6177) and amplified with 14 cycles of PCR. Amplified libraries were then purified with AMPure beads (Agencourt), quantified via Qubit (Life Technologies), and quality was assessed using the BioAnalyzer (Agilent) High-Sensitivity DNA kit. CUT&RUN libraries were pooled and sequenced on the Illumina HiSeq 1500 platform with 100 bp paired-end reads.

Data analysis

Quality scores across sequenced reads were assessed using FASTQC. Illumina adapters were removed using Trim-Galore v.0.6.5. Paired-end reads were mapped to GRCh38 using Bowtie2 (14) v.2.3.4.1 with –very-sensitive-local –dovetail –no-unal –no-mixed –no-discordant parameters. DeepTools (15) was used to filter bam files and remove sequences longer than 120 bp. Peaks were then called using MACS2 (16) v.2.1.1 with –keep-dup = all -f BAMPE options. Consensus peak sets for downstream analysis were derived using IDR (17) with two replicates (Supplemental Figure S1) per target and a cut-off of 0.1. DeepTools (15) v.2.5.1 bamCompare function was used to generate fold enrichment bigWig tracks for each sample compared to the IgG control.

ATAC-seq

ATAC-Seq was performed as previously described (18) using 50 000 Nthy-ORI cells with two biological replicates per condition. Libraries were generated using custom Nextera barcoded primers (18) and were amplified by PCR for a total of 10 cycles. Amplified libraries were then purified with AMPure beads (Agencourt), quantified using a Qubit (Life Technologies), and the quality was assessed using the BioAnalyzer (Agilent) High-Sensitivity DNA kit. ATAC-seq libraries were then pooled in equimolar ratios and sequenced on the Illumina HiSeq 1500 platform with 100 bp paired-end reads.

Quality scores across sequenced reads were assessed using FASTQC. Nextera adapters were removed using Trim-Galore v.0.6.5. Paired-end reads were mapped to GRCh38 using Bowtie2 (14) v.2.3.4.1 with –very sensitive –no-discordant –no-mixed parameters indicated. Reads mapped to chrM were excluded. Peaks were called using MACS2 (16) v.2.1.1 (–nomodel –keep-dup all -p 0.01). Consensus peak sets for downstream analysis were derived using IDR (17) with two replicates (Supplemental Figure S2) per target and a cut-off of 0.05. DiffBind (19) was used to identify regions of differential accessibility. DeepTools (15) v.2.5.1 was used to generate bigWig tracks with the –normalizeUsing RPKM option indicated.

RNA-seq

Nthy-ORI cells were treated 10 nM T3 or vehicle for 6 or 24 h prior to sample collection. Total RNA was extracted and purified using RNeasy Plus Kit (Qiagen) according to manufacturer's protocol. This was repeated to collect a total of three biological replicates per condition. Purity of the total RNA samples was assessed via BioAnalyzer (Agilent) and samples with an RNA integrity score >8 were used for library construction. rRNA was depleted from 1 μg of total RNA with the RiboErase kit (KAPA Biosystems). Strand-specific Illumina cDNA libraries were prepared using the KAPA Stranded RNA-Seq library preparation kit with 10 cycles of PCR (KAPA Biosystems). Library quality was assessed by BioAnalyzer (Agilent) to ensure an average library size of 300 bp and the absence of excess adaptors in each sample. RNA-Seq libraries were pooled in equimolar ratios and sequenced on the Illumina HiSeq 1500 platform yielding an average of 25 million 50 bp single-end reads.

Quality scores across sequenced reads were assessed using FASTQC. All samples were high quality. For alignment and transcript assembly, the sequencing reads were mapped to GRCh38 using STAR (20). Sorted reads were counted using HTSeq (21) and differential expression analysis was performed using SARtools (22) and DESeq2 (23). Genes with a P-value of <0.05 and a log2 fold change greater than 0.5 or less than –0.5 were considered differentially expressed (Supplemental Table S2). Gene Ontology pathway enrichment analysis was performed using Ingenuity Pathway Analysis software.

Proximity labeling by miniTurboID

Cloning of 3xHA-miniTurbo-TRβ

3xHA-miniTurbo-NLS_pCDNA3 vector was a gift from Dr Alice Ting (Addgene plasmid # 107172; http://n2t.net/addgene:107172; RRID: Addgene_107172) (24). 3xHA-miniTurbo-NLS_pCDNA3 was linearized by PCR and THRB cDNA was inserted via HiFi Assembly Cloning (NEB E5520S) to create the 3X-HA-miniTurbo-TRβ vector. Site directed mutagenesis (NEB E00554) was used to create mutant 3X-HA-miniTurbo-TRβ-GS125 and 3X-HA-miniTurbo-TRβ-PV vectors. Successful insertion of the THRB cDNA and site-directed mutagenesis were confirmed by Sanger sequencing. Expression of fusion constructs from the cloned vectors was confirmed by Western blot (Supplemental Figure S3A).

Transfection and biotin labeling

Nthy-ORI cells were grown as a monolayer in DMEM-F12 (Cellgro/Mediatech) supplemented with 10% fetal bovine serum and penicillin-streptomycin (200 IU/l) in 15cm cell culture dishes. Cells were transfected at ∼80% confluency with 20 μg of plasmid DNA using 25 μl Lipofectamine 3000 for 24 h. BioID samples were simultaneously labeled using 250 μM biotin and treated with 10 nM T3 or vehicle for 15 min. Labeling was stopped by placing cells on ice and washing three times with ice-cold PBS. Cells were detached from the plate and collected by centrifugation. The cell pellet was subjected to nuclear protein extraction using the NE-PER Protein Extraction Kit (ThermoFisher 78833) with the addition of Protease Inhibitor Cocktail (Thermo Scientific 781410) per the manufacturer's instructions.

Sample preparation and mass spectrometry

To enrich biotinylated proteins, 1 mg of nuclear extract was incubated for 30 min rotating at room temperature with 100 ul of streptavidin-coated magnetic beads (Invitrogen 65001). The supernatant was removed, and the beads were washed three times with high salt RIPA buffer (100 mM Tris pH 9.0, 500 mM LiCl, 150 mM NaCl, 1% Igepal/NP-40, 1% deoxycholic acid). Washed beads were then resuspended in Laemmli sample buffer and boiled for 10 min to denature and release the biotinylated proteins from the beads. The eluents were loaded onto 10% Tris–glycine gels (Invitrogen XP00100BOX), and separated by SDS-PAGE. Gels were then silver stained (Pierce 24600) prior to band excision for mass spectrometry (Supplemental Figure S4C). LC–MS was performed using an LTQ-Orbitrap instrument (ThermoFisher).

Data acquired by mass spectrometry was quantified using MaxQuant label-free quantification (LFQ) workflow, and LFQ values were used to calculate differential enrichment of identified proteins between experimental conditions using the DEP Bioconductor package (25) (Supplemental Figure S4D, Table S1). Proteins with a P-value of <0.05 and a log2 fold change greater than 1 or less than –1 were considered differentially enriched. Proteins found to be enriched in the empty vector control group were excluded from wildtype and mutant TRβ datasets and were not used for downstream analysis.

RESULTS

TRβ binds to and remodels chromatin in three distinct modes

We profiled the genomic binding patterns of TRβ in the human normal thyroid epithelial cell line Nthy-ORI using CUT&RUN (13). This approach allows high-resolution detection of genomic binding sites of low abundance transcription factors compared with ChIP-seq. Replicate CUT&RUN peaks were used to generate high confidence peak sets for TRβ with and without 10nM T3 for 6 h (Supplemental Figure S1A). We identified a total of 8200 TRβ binding sites that we classified into three distinct groups based on response to T3 treatment (Figure 1A). The largest group is the liganded TRβ binding sites (n = 6768), which have significantly increased TRβ enrichment when T3 is present. Unliganded TRβ binding sites (n = 959) are enriched with TRβ only in the absence of T3, and ligand-independent binding sites (n = 473) are detected both in the presence and absence of T3. The enrichment profile of TRβ binding at each of the three groups of sites was visualized in a heat map (Figure 1B), and in representative genome browser shots (Figure 1C). Accordingly, we observed a nearly complete loss of TRβ enrichment at unliganded binding sites with T3 treatment, while the enrichment remained stable at ligand-independent binding sites with T3 treatment. Liganded binding sites exhibited a significant gain in TRβ enrichment upon T3 treatment. The gain in signal at liganded binding sites is consistent with the dynamic assisted loading model (26,27), the concept whereby nuclear receptors have sparse transient interactions with chromatin and ligand binding increases the residency time to stabilize interactions with response elements.

Figure 1.

CUT&RUN reveals three distinct modes of TRβ binding in thyroid cells. (A) Venn diagram illustrates peak overlap analysis and defines three modes of TRβ occupancy. (B) Heatmap demonstrating regions classified as unliganded, ligand-independent, or liganded. (C) Representative genome browser shots of regions containing each of the three modes of TRβ binding. (D) Distribution of TRβ binding sites annotated to proximal promoters (<1000 bp from TSS), distal regulatory elements (1000 bp–10 kb from TSS), and intragenic regions. (E) Distance from TRβ binding sites to the nearest TSS. (F) Fraction of TRβ peaks containing a full-length DR4 thyroid hormone response element (TRE) or a TRE half-site.

Annotation of TRβ binding sites revealed differences in the distribution of TRβ binding relative to transcriptional start sites (TSS) (Figure 1D). 30% of ligand-independent sites mapped to proximal promoter regions (within 500 bp of TSS), while only 18% of liganded and 9% of unliganded TRβ binding sites mapped to a proximal promoter. Unliganded binding sites were more likely to map to distal regulatory elements (28%), compared to liganded binding sites (19%) or ligand-independent binding sites (10%). A large proportion of TRβ binding sites, regardless of binding mode, were annotated to intergenic regions. Calculating the distance to the nearest TSS of each of these three types of binding sites revealed that ligand-independent binding sites are found closer to transcriptional start sites than others (Figure 1E). Together these results suggest that there is significant TRβ functional activity directed towards more indirect mechanisms for regulation of gene expression such as regulation of histone modifications and chromatin accessibility across broad regions, and maintenance of higher order chromatin structure. Transcription factor motif analysis revealed that thyroid hormone response elements (TRE) were strongly enriched in each of the three groups of TRβ binding sites (Supplemental Figure S2). Ligand-independent binding sites had the highest frequency of the full-length direct-repeat palindromic TRE (DR4), while liganded sites have the highest frequency of TRE half-sites (Figure 1F). Binding sites that contain full-length TREs are likely to be bound by TRβ homodimers or TRβ/retinoid X receptor (RXR) heterodimers, while sites that contain half-sites are more likely to be bound by TRβ monomers based on prior studies of TRβ interaction with response elements (3,6,28–30).

To map changes in chromatin accessibility in response to T3 in Nthy-ORI cells, we performed ATAC-seq after 6 h (early) and 24 h (late) after treatment with T3 (Supplemental Figure S3). Differential accessibility analysis revealed a significant increase in the chromatin accessibility with early and late T3 treatments (7,754 and 21 678 opened regions, respectively; FDR < 0.05). In contrast, we found that only 107 regions were closed by T3 treatment, all of which were closed with early T3 treatment (Figure 2A). We next examined the chromatin accessibility associated with TRβ binding. Unliganded TRβ binding sites, which are lost upon T3 treatment, had relatively little T3-induced chromatin accessibility with a modest decrease in accessibility with T3 treatment (Figure 2B). However, chromatin accessibility increased at liganded binding sites with early T3 treatment and further increased with late T3 treatment (Figure 2C). Ligand-independent binding sites also had significant increases in chromatin accessibility after 6 and 24 h of T3 treatment (Figure 2D). The bimodal distribution of ATAC-seq tag density observed near liganded binding sites (Figure 2D) that is not observed near unliganded or ligand-independent binding sites may be due to an enrichment of bidirectional promoters within this subset of regions (Supplemental Figure S3E). To estimate the size of differentially accessible regions of chromatin following T3 treatment, we compared the average peak width of differentially accessible ATAC-seq peaks within 5 kb of each of the groups of TRβ binding sites (Figure 2E). Unliganded TRβ binding sites had the smallest average peak width (190 bp), consistent with the modest changes in accessibility near those sites. The average width of the differentially accessible peaks was greatest for ligand-independent TRβ binding sites (355 bp) followed by liganded sites (253 bp). This suggests that remodeling around liganded binding TRβ binding sites is more focused while remodeling around ligand-independent TRβ binding sites stretches across broader regions.

Figure 2.

T3 induces changes in chromatin accessibility. (A) Scatter plot highlights differentially accessible ATAC-seq peaks after 6 and 24 h of T3 treatment. T3 induced accessibility at 7754 peaks after 6 h (light pink) and an additional 21 678 peaks after 24 h (dark pink). 107 peaks were closed by T3 (blue). Differentially accessible peaks are defined as having a log2FC ≥ 1 or ≤–1 and an FDR ≤ 0.05. ATAC-seq tag density is plotted near unliganded (B), ligand-independent (C), and liganded (D) TRβ CUT&RUN peaks. (E) Distribution of peak width of differentially accessible ATAC-seq peaks within 5kb of TRβ binding sites. Statistical significance was determined by one-way ANOVA followed by multiple comparisons; * indicates P < 0.001.

T3 induces changes in the TRβ interactome

To detect protein–protein interactions between TRβ and potential coregulators, including transient cofactors that are exchanged in a T3-dependent fashion, we used a live cell proximity-dependent biotin labeling assay (Figure 3A). Nthy-ORI cells were transfected with vectors that express a TRβ-miniTurboID fusion construct that can rapidly biotinylate proximal proteins, which were isolated via biotin-affinity purification and subsequently identified by mass spectrometry (Figure 3A, Supplemental Figure S4). We detected a total of 1328 high-confidence proteins that interact with TRβ either in the presence and absence of T3. Several well-characterized TRβ binding partners were identified in our dataset including RXR, HDAC1, HDAC2, SMRT and THRAP (Supplemental Figure S4B, Table S1). Differential enrichment analysis revealed that 75 interactions were gained in the presence of T3, and 70 were lost in the presence of T3 (Figure 3B, Supplemental Table S1). The largest group are the unchanged interactions (1183 proteins). Pathway analysis revealed an enrichment of distinctive biological processes associated with the different groups of TRβ interaction partners that change with T3 treatment (Figure 3D). Interacting proteins that were lost with T3 treatment were enriched for transcription repressor activity, demethylase activity, and translation factor activity, while gained proteins were enriched for modification-dependent protein binding, C-terminal protein binding, and RNA polymerase binding activity. Interacting proteins that remained stable were classified as transcriptional coactivator, nucleosome remodeling, ATPase activity and histone binding proteins.

Figure 3.

The TRβ interactome is altered by T3 binding. (A) Schematic of miniTurboID proximity labeling assay used to identify TRβ-binding partners. (B) Volcano plot depicts differentially enriched TRβ binding partners (DEPs) defined as gained with T3 treatment (75), lost with T3 treatment (70), or stable (1,183). Differentially enriched proteins are defined as having a log2FC ≥ 1 or ≤ –1 and a P-value ≤0.05. (C) Venn diagrams show DEPs that interact with wildtype TRβ compared with TRβPV (LBD mutant), and TRβGS125 (DBD mutant). (D) GO Molecular Function analysis of gained, lost, and stable TRβ binding partners. Color corresponds to P-value and size corresponds to the number of interacting proteins within that pathway divided by the total (enrichment score). (E) Multi-subunit chromatin remodeling complexes found to interact with TRβ, and each member of the complex identified in our proteomic dataset. (F) Complex interaction scores, defined as the sum of the signal intensity of each subunit in the complex divided by the total subunits identified, of chromatin remodeling complexes.

To examine the impact of ligand-dependent interactions and DNA-binding dependent interactions, we compared the interaction profiles of a ligand-binding domain mutant (TRβPV) (31) and a DNA-binding domain mutant (TRβGS125) (32) with that of wildtype TRβ. The T3-induced differential interactions with wildtype TRβ were altered by TRβPV and TRβGS125 mutations (Figure 3C). The TRβPV mutant retained many of the wildtype interactions, but did not exhibit differential binding upon T3 treatment. This suggests that the conformational changes that result from ligand binding are needed to stabilize a small number of protein-protein interactions with TRβ. However, this relatively small list of binding partners may be critical for facilitating T3-dependent changes in chromatin accessibility and gene expression. The loss of 826 interactions with TRβGS125 compared to wildtype TRβ indicates that many of these are dependent upon direct TRβ interaction with DNA for proper docking, consistent with the concept of specific TRβ recruitment of coregulator complexes to targeted sites. Our analysis of the TRβ interactome is consistent with a co-regulatory shift model (8), where, rather than an all-or-nothing switch, T3 alters the ratio of corepressor to coactivator binding partners.

We focused our analysis on interacting proteins that were likely to participate in chromatin remodeling. Several multi-subunit chromatin remodeling complexes were identified in our proximity labeling assay such as SWI/SNF, Mi-2/NURD, the NCOR and SIN3 co-repressors and the CBP/p300 and TRAP/DRIP/Mediator co-activators (Figure 3E). Each of these has been previously linked to TRβ gene regulation (33–37). We identified several members of each of these complexes as in our proteomic dataset. An interaction score was calculated for TRβ-associated chromatin remodeling complexes by dividing the sum of the signal intensity of each subunit in the complex by the total subunits identified to offset inflation of the interaction score based on the number of subunits in the complex (Figure 3F). The Mi-2/NURD, SIN3 and TRAP/DRIP/Mediator complexes did not have significant changes in their interaction scores. There was an increase in the CBP coactivator complex interaction score and decrease in the NCOR corepressor complex interaction score with T3, consistent with previous studies of their interaction with TRβ (33). The SWI/SNF complex also significantly increased with T3. TRβPV and TRβGS125 mutations both inhibited these T3-dependent interactions.

TRβ differentially recruits SWI/SNF complexes

Upon closer examination of the SWI/SNF complex subunits that were identified as TRβ interacting proteins, it became apparent that TRβ interacts with two different subclasses of SWI/SNF complexes: canonical BAF and polybromo-associated BAF (PBAF). Each of these contains one of the mutually exclusive ATPase catalytic subunits, (BRG1 or BRM) which are necessary for nucleosome displacement, array of accessory components, and a few subclass-specific subunits (38) (Figure 4A). Our proteomic analysis identified the BRG1 core subunit, many of the accessory components, and three subclass-specific subunits. Enrichment of the BAF-specific subunit ARID1A (BAF250) was decreased with T3 treatment, while the PBAF-specific subunits ARID2 (BAF200) and PBRM1(BAF180) were increased in the presence of T3 (Figure 4B). BAF57 enrichment was also increased in the presence of T3. BRG1 was not differentially enriched.

Figure 4.

BAF and PBAF complexes differentially interact with TRβ. (A) Diagram of BAF and PBAF complex subunits. Core subunits are colored in grey, BAF specific subunits are colored in blue, and PBAF specific subunits are colored in pink. Labels of subunits identified as TRβ binding partners by miniTurboID proximity-labeling assay are bolded. (B) Signal intensity of differentially enriched SWI/SNF subunits identified in miniTurboID proximity-labeling assay in the presence and absence of T3. Significance (*P < 0.05) was determined by paired t-test. (C) Clustered heatmap of genome-wide SWI/SNF subunit signal shows patters of differential occupancy. Plotted signal for BRG1, ARID1A, and PBRM1 is RPKM normalized and log scaled. Median signal for each factor within each cluster is plotted to the right. Bar plots show the direct overlap between peaks within each cluster with TRβ binding sites. Box plots demonstrate differential CUT&RUN tag density of BRG1 (D), ARID1A (E) and PBRM1 (F) subunits at TRβ binding sites. Significance (*P < 0.05) was determined by paired t-test.

Based on differential enrichment observed in our proximity ligation data, we performed CUT&RUN to determine T3-induced changes in genomic binding of ARID1A as a representative of BAF complexes, PBRM1 as a representative of PBAF complexes, as well as BRG1 (Supplemental Figure S1B–D). CUT&RUN peaks for all three factors were clustered to reveal genome-wide patters of BAF and PBAF complex localization (Figure 4C). Cluster 1 represents the group of peaks with BRG1 binding, but low ARID1A and PBRM1 signal. These BRG1-enriched peaks overlap most with unliganded TRβ binding sites (20%). Cluster 2 peaks have BRG1 and ARID1A binding, indicating the presence of BAF complexes, and overlap most frequently with liganded TRβ binding sites (4%). Peaks within cluster 3 have BRG1 and PBRM1 binding, denoting PBAF complex localization, and overlap with the largest fraction of unliganded TRβ binding sites (7%). PBRM1 signal increases upon T3 treatment at these sites. Cluster 4 sites show enrichment of both BAF and PBAF complexes. These sites overlap with a large fraction of all three types of TRβ binding sites, but most frequently with ligand-independent binding sites (78%).

We then examined the tag density of each of these factors at TRβ binding sites. BRG1 CUT&RUN tag density was unchanged with T3 at unliganded sites, increased at ligand-independent sites, and decreased at liganded sites (Figure 4D). This indicates that unliganded TRβ may recruit BRG1 to prime its liganded binding sites a manner similar to a mechanism described for the glucocorticoid receptor (39). ARID1A tag density was low with a slight increase with T3 at unliganded sites, increased with T3 at ligand-independent sites, and stable at liganded sites (Figure 4E). PBRM1 tag density was low at unliganded sites, high with a slight increase with T3 at ligand-independent sites, and stable at liganded sites (Figure 4F). Of note, ligand-independent sites had an increase in all three SWI/SNF components with T3, indicating that both BAF and PBAF complexes are recruited.

To further examine differential recruitment of BAF and PBAF complexes, TRβ peaks were annotated based on whether they occurred at a proximal (<1 kb from nearest TSS) or at a distal (1–10 kb from nearest TSS) regulatory site (Figure 5A). The group of proximal TRβ binding sites (n = 1330) is made up of liganded binding sites (61%), ligand-independent binding sites (32%), and a small number of unliganded binding sites (7%) (Figure 5B). Nearly all the ligand-independent binding sites fall within this group (425/473). The distal sites (n = 2801) are comprised of mostly liganded binding sites (86%), some unliganded binding sites (13%), and a small fraction of ligand-independent binding sites (0.4%) (Figure 5B). CUT&RUN and ATAC-seq signal at proximal and distal TRβ binding sites was then visualized in heatmaps (Figure 5C) and as line plots (Figure 5D). This analysis revealed distinct patterns of BAF and PBAF recruitment to TRβ binding sites. As expected, BRG1 was recruited to both proximal and distal sites in the presence and absence of T3. ARID1A was recruited in a T3-dependent manner to both TRβ-bound proximal and distal sites. PBRM1, however, was recruited preferentially to proximal sites and not to distal sites. The T3-induced increase in chromatin accessibility is greater near proximal sites than distal binding sites. Given that both BAF and PBAF complexes are recruited to promoters, both complexes may be required to facilitate the large changes in chromatin accessibility observed.

Figure 5.

BAF and PBAF complexes are differentially recruited to proximal and distal TRβ binding sites. (A) Venn diagram illustrates distribution of TRβ peaks classified as proximal (less than 1 kb from nearest TSS) and distal (1–10 kb from nearest TSS) binding sites. Binding sites >10 kb were excluded. (B) Stacked bar blot shows the proportion of proximal and distal binding sites that correspond to each of the three modes of TRβ binding. The proximal group is made up of ligand-independent (34%) and liganded (62%) sites. The distal group is made up of unliganded (11%) and liganded (88%) sites. (C) Heatmap shows RPKM normalized CUT&RUN signal of TRβ, BRG1, ARID1A, PBRM1 and ATAC-seq signal. (D) CUT&RUN tag density of BRG1, ARID1A (BAF-specific) and PBRM1 (PBAF-specific) near proximal TRβ binding sites classified as proximal or distal. BRG1 and ARID1A are recruited to proximal and distal sites, PBRM1 is preferentially recruited to proximal sites. ATAC-seq tag density is increased upon T3 treatment for 6 and 24 h near both proximal and distal TRβ binding sites.

TRβ chromatin interactions are correlated with target gene expression

To determine the effect of T3 treatment on gene expression, Nthy-ORI cells were treated with 10nM T3 or vehicle control for 6 and 24 h and global transcriptomic analysis by RNA-seq was performed. Differential gene expression analysis was performed comparing each treatment to the corresponding control. We determined that 366 and 368 genes were up and downregulated, respectively, by T3 at 6 h (Figure 6A, Supplemental Table S2). The effect was increased at the 24-h period where 480 and 667 up and downregulated genes, respectively. Distinct groups of pathways were altered by T3 at early and late time points (Supplemental Table S4). Enriched biological function pathways within early and late DEGs were compared using Ingenuity Pathway Analysis software (Figure 6C). Highly enriched upregulated pathways included cellular homeostasis, survival, viability and cell cycle progression. Transcription and protein synthesis switched from down- to upregulated between early and late T3 treatment, while carbohydrate metabolism switched from up- to downregulated. Apoptosis, cell transformation, and ER stress response were downregulated. BETA (40) was used to predict whether T3-regulated genes are likely to be direct transcriptional targets of TRβ based on proximity of a TRβ CUT&RUN peak to the TSS (Figure 6B). 63% (230/366) of early upregulated genes and 51% (245/480) of late upregulated genes were predicted to be direct targets of TRβ. This suggests that increases in chromatin accessibility that occur in the time between early and late T3 timepoints (Figure 2A) may allow TRβ to access additional binding sites for induction of a subsequent set of target genes. Conversely, 48% (177/368) of early and 18% (120/667) of late downregulated genes were predicted to be direct targets of TRβ, which indicates that direct TRβ-induced gene repression occurs more immediately and the additional downregulation we observed may be secondary effects downstream of the early transcriptomic effects of T3-treatment.

Figure 6.

TRβ chromatin interactions are correlated with target gene expression. (A) T3-induced and -repressed genes determined by RNA-seq after 6 and 24 h of T3 treatment in two independent biological replicates. Differentially expressed genes (DEGs) are defined as having a log2FC ≥ 0.5 or ≤ –0.5 and an P-value ≤0.05. (B) Activating/repressive function prediction of TRβ direct targets among DEGs after 6 and 24 h of T3 treatment. The red and the blue lines represent upregulated and downregulated genes, respectively. The dashed line indicates the non-differentially expressed genes as background. Genes are cumulated by rank based on the regulatory potential score from high to low. (C) IPA analysis reveals differentially enriched biological functions after early and late T3 treatment. (D) Proportion of DEGs with each of the three types of TRβ binding sites within 10 kb of the TSS. (E) Average log2(Fold Change) of DEGs with 0, 1, 2, 3 or >4 TRβ binding sites within 10 kb of the TSS. (F) Proportion of DEGs with ARID1A (yellow) or PBRM (red) binding exclusively, both ARID1A and PBRM1 (blue), or neither within 10 kb of the promoter. (G) Genome browser shots highlight differential recruitment of SWI/SNF complexes to a representative TRβ target gene. Arrows are pointing to a liganded TRβ binding site upstream KLF9 promoter, and two ligand-independent TRβ binding sites within the proximal promoter.

TRβ peaks associated direct targets predicted by BETA (Figure 6B) were then annotated based on the mode of TRβ binding they belonged to (Figure 6D). Ligand-independent peaks make up the largest proportion of binding sites that are predicted to have direct regulatory potential, regardless of timepoint or direction, because they are most likely to be located near the TSS. Liganded peaks were also predicted to have regulatory potential, particularly for upregulated DEGs. Unliganded peaks were the least likely to directly regulate a DEG. However, many DEGs have multiple TRβ peaks nearby, which was not considered in the BETA prediction. We therefore grouped the DEGs based on the number of TRβ peaks that were contained within 10 kb of their TSS (0, 1, 2, 3 or >4). This analysis illustrated that the more TRβ peaks are located near the TSS of a target gene, the greater the average fold change is in response to T3 treatment at both 6 and 24 h (Figure 6E). This implies that access to multiple TRβ sites that exhibit different functional modes of binding may be needed to facilitate the larger fold changes in gene expression that we observed.

We examined SWI/SNF complex recruitment to promoters of T3-induced differentially expressed genes (Figure 6F). A majority of early upregulated genes have a SWI/SNF complex binding site within 10 kb of the TSS; most have both ARIDIA and PBRM1 peaks, while a minority had exclusively one or the other. There was a similar trend in the late upregulated genes. In contrast, most downregulated genes, both early and late, did not have SWI/SNF binding near their promoter. However, among those that did, there was a preference for ARID1A recruitment over PBRM1. Visualization of representative binding sites in genome browser shots (Figure 6G) further demonstrated differential recruitment. Combined, these results suggest that SWI/SNF complexes participate in T3-induced gene regulation near promoters of target genes, particularly in the context of upregulation. For upregulation of gene expression, both BAF and PBAF complexes are likely to be recruited to remodel the proximal promoter region, while BAF complexes may be preferred for remodeling of downregulated promoters.

DISCUSSION

Although it is now appreciated that the classic model for TRβ interaction with chromatin is oversimplified, a consensus has yet to be reached on a more complex and nuanced model. Based on the data presented here, we propose a multi-modal regulation model where TRβ has at least three distinct modes of binding and associated chromatin remodeling (Figure 7). In agreement with previous genome-wide studies (4–6,8), we observed significant shifts in TRβ occupancy in the presence and absence of T3. At unliganded sites, TRβ binds to chromatin with limited effects on chromatin accessibility; T3 treatment results in a loss of binding. Liganded binding sites represent the vast majority of TRβ binding, and they are associated with significant changes in localized chromatin accessibility (Figures 1A, B and 2C). Intriguingly, liganded binding sites show a clear increase in enrichment upon the addition of T3, however most sites have some, albeit low, enrichment before T3 is added (Figure 1B). This suggests that TRβ interactions with this type of site are transient and are stabilized by ligand binding and recruitment of coregulators, consistent with a dynamic assisted loading model (27,41). Since they are numerous and broadly distributed across the genome, there are likely multiple functions of liganded binding sites that further studies may clarify. While some of these binding sites occur near proximal promoters and facilitate direct induction or repression of gene expression, many occur in distal regulatory elements and may regulate enhancers or coordinate higher order chromatin structure. Ligand-independent binding sites are defined by enrichment both in the presence and absence of T3, and substantial induction of chromatin accessibility (Figures 1A, B and 2D). While they have a clear functional importance for regulation of gene expression, these binding sites represent a small minority of the overall cistrome of TRβ. Notably, ligand-independent binding sites have many of the characteristics originally described in the bimodal switch model such as their proximity to transcriptional start sites, high frequency of full-length DR4 TREs, and high regulatory potential.

Figure 7.

Model for TRβ chromatin interaction and recruitment of SWI/SNF complexes. (A) BRG1 is recruited to unliganded TRβ binding sites within inaccessible chromatin for modest changes in chromatin accessibility. (B) BAF complexes are specifically recruited to liganded TRβ binding sites to facilitate T3-induced changes in chromatin accessibility at distal regulatory regions. (C) BAF and PBAF complexes are both recruited to ligand-independent TRβ binding sites near the TSS of target genes to facilitate broad changes in chromatin accessibility.

TRβ binding sites are also characterized by differential recruitment of SWI/SNF chromatin remodelers. BRG1 alone is recruited to unliganded and liganded sites in the absence of T3 (Figure 4C) suggesting that that it may be recruited to prime binding sites for TRβ occupancy, a mechanism which has been clearly delineated for the glucocorticoid receptor (39). BAF and PBAF complexes, each with unique subunits which dictate their precise function, are recruited to TRβ binding sites within promoters and likely contribute directly to the changes in chromatin accessibility and recruitment of other transcription factors to alter target gene expression. Both BAF and PBAF complexes have been implicated previously in hormone-dependent gene regulation (42,43), however the distinct functional role of each when they are recruited by a nuclear receptor to the same location remains unclear. BAF complexes are preferentially recruited by TRβ to binding sites that occur in distal regulatory regions. This might be a mechanism by which TRβ primes these binding sites or organizes higher order chromatin structure to promote persistent changes in gene expression. These multifaceted interactions with a single chromatin remodeling complex illustrate the importance of the proposed coregulator shift model (8) over a complete coregulator switch. Importantly, while the recruitment of SWI/SNF complexes tracks with much of the T3-induced chromatin accessibility we observed, it does not directly overlap with all the differentially accessible regions. The additional changes that do not directly overlap with a SWI/SNF binding site are likely regulated by other chromatin remodelers or are a result of regulation within a region of higher order chromatin organization.

As a deeper understanding of the variety of ways in which TRs regulate gene expression and coordinate a network of cofactors is developed, it is important that the models we use reflect their multidimensional function. We suggest a model for multi-modal regulation by TRβ that has at least three distinct modes defined by their T3-dependent occupancy, changes in accessibility, and differential recruitment of chromatin remodelers. Further investigation of other nuclear complexes that contribute to TRβ regulation of gene expression, association with binding sites, and chromatin remodeling, will enhance our understanding of the nature of TRβ interaction with the epigenome. Collectively, this study provides a next-generation model for TRβ interactions with chromatin, and lays a foundation for further studies of TRβ regulation of gene expression and recruitment of key cofactors in both normal cells and in disease models.

DATA AVAILABILITY

All raw and processed next generation sequencing data associated with this manuscript is available for download from the NCBI Gene Expression Omnibus repository under accession code GSE168954. The proteomics data associated with this manuscript is available for download from the PRIDE repository under accession code PXD030403.

Click here for additional data file.