Anti-diabetic rosiglitazone remodels the adipocyte transcriptome by redistributing transcription to PPARγ-driven enhancers

Regulation of nascent gene transcription by rosi correlates with changes in mRNA levels

To assess the temporal and gene-specific relationship between nascent gene transcription and steady-state mRNA levels, we determined the adipocyte transcriptome using gene expression microarrays after 30 min, 1 h, 2 h, 6 h, 12 h, 24 h, 36 h, and 48 h of rosi treatment, with three biological replicates at each time point (Fig. 2). While the nascent transcription of many genes was regulated as early as 10 and 30 min after rosi, very few mRNA transcripts changed during this time period and for at least 2 h after rosi treatment. This was not surprising, as the time required to reach new steady-state levels is related to the rate of degradation rather than the rate of synthesis (Schimke and Doyle 1970). However, the correlation between nascent transcription and steady-state mRNA levels was high at later time points, with the greatest correlation noted between the transcription regulation at 3 h and the mRNA regulation measured 6 h after rosi. The lower correlation with later microarray time points suggests that steady-state regulation at these later times may be dependent on secondary transcriptional changes occurring later than 3 h of treatment. Importantly, direct transcriptional regulation by rosi in this model was highly similar to recently published regulation of in vivo adipose tissue gene expression in rosi-treated mice, suggesting that many of the rosi-dependent changes that we describe are physiologically relevant (Supplemental Fig. S1; Ohno et al. 2012).

Open in a separate window

Figure 2.

Adipocyte gene body transcription levels correlate with steady-state mRNA levels. Correlation between microarray mRNA levels and GRO-seq transcription levels is plotted at each time point. The Pearson correlation coefficient is given for each pair of time points.

Rosi regulates transcription of eRNAs that correlate with gene transcription

In addition to gene body transcription, GRO-seq revealed robust bidirectional transcripts at enhancers or eRNAs (Core et al. 2008; Kim et al. 2010), which were identified and quantified in an unbiased, genome-wide analysis (Fig. 3A). For example, bidirectional eRNAs were identified at enhancers upstream of the Fabp4 locus, and their transcription was observed to be up-regulated by rosi (Fig. 3B). Indeed, unbiased de novo calling of bidirectional intergenic transcripts confirmed that the transcription of many eRNAs was strongly and rapidly regulated by rosi, and down-regulated eRNAs greatly outnumbered up-regulated eRNAs (Fig. 3C). This was similar to the effect of rosi on gene body transcription, and, as has been observed in other systems (De Santa et al. 2010; Kim et al. 2010), the effect of rosi on eRNA transcription correlated strongly with transcription at the nearby gene bodies (Fig. 3D; Supplemental Fig. S2A). The correlation was confirmed at eRNA/gene body pairs by quantitative PCR (qPCR), which also demonstrated that eRNA induction often preceded gene induction (Fig. 3E; Supplemental Fig. S2B). The correlation was validated at repressed eRNA/gene pairs as well (Supplemental Fig. S2C). Although intragenic eRNAs were excluded from downstream analysis because of difficulties in identifying them reliably, we observed that many follow a similar pattern of correlation with the target gene (Supplemental Fig. S3).

Open in a separate window

Figure 3.

Adipocyte eRNA transcription is stimulated by rosi and correlates with gene body transcription. (A) Genome-wide average signal of intergenic bidirectional transcripts in untreated adipocytes from the plus and minus strands. (B) Two bidirectional eRNAs are transcribed at enhancers upstream of the Fabp4 TSS and up-regulated by rosi treatment. eRNA centers are indicated by arrows. (C) Heat map showing all rosi-regulated eRNAs found in an unbiased manner (N = 2251). (D) Correlation between rosi-regulated eRNAs and the regulation of the nearby gene for all pairs of matching time points (N = 462). For each gene, eRNAs within 100 kb of the TSS were included in the analysis. (E) Correlation between gene and eRNA rosi regulation as measured by RT-qPCR for two example genes, Fabp4 and Pdk4. N = 3. Error bars indicate SEM.

PPARγ directly mediates the induction of eRNA and gene transcription by rosi

De novo motif-finding analysis at sites of rosi-induced eRNAs revealed strong enrichment for a sequence that is highly similar to the canonical PPARγ/RXR-binding site (Fig. 4A; Supplemental Fig. S4). Indeed, up-regulated eRNAs were extremely likely to have PPARγ bound nearby (85%), much more so than unregulated or down-regulated eRNAs (Fig. 4B). In fact, down-regulated eRNAs were relatively devoid of PPARγ binding, as discussed below. Furthermore, strong PPARγ binding as measured by total normalized tag counts from chromatin immunoprecipitation (ChIP) coupled with deep sequencing (ChIP-seq) was enriched at up-regulated eRNAs, but not at down-regulated eRNAs, relative to unregulated eRNAs (Fig. 4C). PPARγ-binding sites that overlap with eRNAs had higher PPARγ occupancy than those that lack eRNAs, further suggesting that enhancers with eRNAs are more likely to be functional (Supplemental Fig. S5). Moreover, knockdown of PPARγ abrogated basal eRNA transcription at rosi-up-regulated sites (Fig. 4D; Supplemental Fig. S6), indicating that PPARγ binding was required for eRNA transcription at these sites.

Open in a separate window

Figure 4.

Rosi-up-regulated eRNAs are dependent on genomic binding of PPARγ. (A) Top hit from Homer de novo motif search at up-regulated eRNAs. The closest known motif, the DR1 motif, is shown for reference. The enrichment for this motif was 45% in the target sites and 13.5% in the background sites. (B) Percentage of up-regulated, down-regulated, and unregulated eRNA sites that overlap with a called PPARγ peak from ChIP-seq. There were 14,604 called PPARγ peaks in our data set. Among those, 5819 sites were extragenic, and 3642 sites had eRNAs. (C) Total PPARγ tag count in reads per million (RPM) within 1 kb of up-regulated, down-regulated, and unregulated eRNA sites. (*) P = 7.7 × 10⁻⁹⁵ versus unregulated; (**) P = 1.8 × 10⁻²⁰ versus unregulated. (D) RT-qPCR of eRNAs 24 h after siRNA knockdown of PPARγ. N = 3. Error bars indicate SEM. (E) Distance from TSS to the closest PPARγ sites for regulated genes. P < 10⁻¹⁰ for up-regulated versus unregulated sites by χ² test. (F) Number of PPARγ-binding sites within 100 kb of the TSS for regulated genes. P < 10⁻¹⁰ for up-regulated versus unregulated sites by χ² test.

In addition to regulating eRNAs, PPARγ binding was critical for the regulation of gene body transcription as well. Up-regulated genes were statistically significantly enriched for PPARγ sites closer to the transcription start site (TSS) compared with down-regulated or unregulated genes (Fig. 4E). Furthermore, up-regulated genes had significantly more PPARγ-binding sites within 100 kb of the TSS (Fig. 4F). Together, these data strongly support the model that rosi up-regulates gene expression by binding to PPARγ at enhancer sites proximal to the TSS, where it directly stimulates eRNA as well as gene body transcription.

Repression of eRNA transcription by rosi is not directly mediated by PPARγ

In contrast to the up-regulation of transcription by rosi, the PPARγ motif and strength of binding were not enriched at rosi-repressed eRNAs. Although 23.6% of down-regulated eRNAs had PPARγ bound (Fig. 4B), this percent was far less than at unregulated eRNAs, and the binding tended to be extremely weak (Fig. 4C). In addition, in contrast to its effect on up-regulated eRNAs, knockdown of PPARγ did not affect basal expression of rosi-down-regulated eRNAs (Fig. 4D). Moreover, unlike up-regulated genes, the number, strength, and proximity of PPARγ binding were not enriched near the TSS of down-regulated genes relative to unchanged genes (Fig. 4E,F). Together, these data strongly suggest that PPARγ does not mediate the repression of eRNA and gene body transcription by rosi directly; i.e., by binding in cis with the regulated gene body, as was the case at the majority of rosi-up-regulated eRNAs.

De novo motif analysis of sites with down-regulated eRNAs produced two highly enriched motifs that most resembled a C/EBP:AP-1 hybrid motif and the canonical AP-1 motif (Fig. 5A; Supplemental Fig. S7). ChIP-seq analysis for C/EBP and AP-1 factors that are abundant in adipocytes revealed significant enrichment at down-regulated eRNAs of both C/EBPα and FOSL2 (Fig. 5B), which has been shown to play a role in adipocyte gene expression (Wrann et al. 2012). A higher percentage of down-regulated eRNAs, compared with up-regulated and nonregulated, had C/EBPα and FOSL2 bound, and, interestingly, although C/EBPα has been shown to be enriched near genes repressed after rosi treatment (Vernochet et al. 2009; Haakonsson et al. 2013), the enrichment was higher for FOSL2 (Fig. 5C). Other factors tested, including C/EBPβ, ATF2, and JUND, also showed enrichment at down-regulated eRNAs (Supplemental Fig. S8), and the strength of this binding did not change upon rosi treatment (Supplemental Fig. S9). The finding that both C/EBP factors and all three AP-1 factors are enriched at down-regulated sites suggests that each of these factors, potentially along with other TFs but notably not PPARγ, is located at eRNAs in cis with gene bodies whose transcription is negatively regulated by rosi.

Open in a separate window

Figure 5.

Enrichment of C/EBP and AP-1 factors at down-regulated eRNAs. (A) Top two hits from Homer de novo motif search at down-regulated eRNAs. The closest known motif for each is shown for reference. (B) Total C/EBPα and FOSL2 tag counts in reads per million (RPM) within 1 kb of the center of regulated eRNAs from ChIP-seq. (*) P = 4.6 × 10⁻⁸ versus unregulated; (**) P = 8.9 × 10⁻⁴⁶ versus unregulated. (C) Percentage of regulated eRNAs that overlap with called C/EBPα or FOSL2 peaks.

Gene expression analysis

RNA was isolated from cells using TRIzol (Invitrogen) followed by the RNeasy minikit (Qiagen). RT–PCR was performed using 1 μg of RNA (Applied Biosystems) following the manufacturer’s instructions, and qPCR was performed using primers listed in Supplemental Table S1 using Power SYBR Green master mix (Applied Biosystems) on the PRISM 7500 and 7900HT instruments (Applied Biosystems). Analysis was performed using the standard curve method, and all genes were normalized to the housekeeping gene Arbp. For the microarray, RNA integrity was examined using an Agilent 2100 Bioanalyzer. RNA samples (150 ng) with RNA integrity number >7 were used for target amplification and labeling via the Ambion WT Expression kit (no. 4411974) and Affymetrix WT Terminal-Labeling kit (no. 900671) following the manufacturer’s protocol. Mouse Gene 1.1 ST array plates (no. 901418, Affymetrix) were used for microarray hybridization, wash, stain, and scan with GeneTitan hyb-wash-stain kits (no. 901622, Affymetrix) and a GeneTitan instrument. GeneTitan scanner data were collected with default parameters and further analyzed using the Partek genomics suite. Data were normalized using the default Robust Multichip Average (RMA) method. For each gene, an average was taken across replicates for the comparative analysis with GRO-seq gene transcriptional level.

GRO-seq library preparation

GRO-seq was performed as previously described (Core et al. 2008; Wang et al. 2011). Cells were washed twice with ice-cold PBS and then swelled in cold swelling buffer (10 mM Tris at pH 7.5, 2 mM MgCl₂, 3 mM CaCl₂) for 5 min on ice. Cells were scraped off and centrifuged at 400g for 10 min, resuspended in 10 mL of lysis buffer (swelling buffer with 10% glycerol, 1% Igepal), and incubated for 5 min on ice. Nuclei were washed twice with lysis buffer and then resuspended in freezing buffer (50 mM Tris at pH 8.3, 40% glycerol, 5 mM MgCl₂, 0.1 mM EDTA). Nuclei were counted and pelleted, and 5 × 10⁶ nuclei were resuspended in 100 μL of freezing buffer. For each library, run-on was performed on four tubes of 5 × 10⁶ nuclei.

For the run-on, cells were mixed with an equal volume of run-on buffer (10 mM Tris at pH 8.0; 5 mM MgCl₂; 1 mM DTT; 300 mM KCl; 20 U of SUPERase-In; 1% Sarkosyl; 500 μM ATP, GTP, and Br-UTP; 2 μM CTP) preheated to 30°C and incubated for 5 min at 30°C. Nuclear RNA was extracted with TRIzol (Invitrogen) and chloroform and precipitated with NaCl and ethanol overnight. The RNA pellet was resuspended in water, and the RNA was DNase-treated (Ambion) for 30 min. RNA was hydrolyzed using fragmentation reagents (Ambion) for 13 min at 70°C and purified through a Micro Bio-Spin p-30 column (Bio-Rad) according to the manufacturer’s instructions. RNA was incubated with 1.5 μL of T4 polynucleotide kinase (New England Biolabs) for 1 h at 37°C and then with an additional 1 μL for 1 h more. RNA was denatured for 5 min at 65°C.

Anti-BrU agarose beads (Santa Cruz Biotechnology) were rotated for 1 h in blocking buffer (0.5× SSPE, 1 mM EDTA, 0.05% Tween-20, 0.1% PVP, 1 mg/mL BSA). Run-on RNA was rotated with beads for 1 h, followed by 5-min washes twice in binding buffer (0.5× SSPE, 1 mM EDTA, 0.05% Tween-20), twice in low-salt buffer (0.2× SSPE, 1 mM EDTA, 0.05% Tween-20), once in high-salt buffer (0.5× SSPE, 1 mM EDTA, 0.05% Tween-20, 150 mM NaCl), and twice in TET buffer (TE at pH 7.4, 0.05% Tween-20). BrU-labeled RNA was eluted from the beads four times for 15 min with 100 μL of elution buffer preheated to 42°C. RNA was ethanol-precipitated overnight.

The RNA pellet was resuspended in water, denatured, and treated with poly(A)-polymerase (New England Biolabs) for 30 min at 37°C. cDNA synthesis was performed as described previously (Wang et al. 2011) using the oNTI223 primer (5′-pGATCGTCGGACTGTAGAACTCT;CAAGCAGAAGACGGCATACGATTTTTTTTTTTTTTTTTTTTVN-3′), where “p” represents 5′ phosphorylation, “;” represents the abasic dSpacer furan, and “VN” represents degenerate nucleotides. The reaction was treated with 3 μL of exonuclease I (Fermentas) for 15 min at 37°C followed by 2 μL of 1 M NaOH for 20 min at 98°C and neutralized with 1 μL of 2 M HCl. cDNA was run on a 10% TBE-urea gel, and products were excised and eluted from shredded gel pieces for 4 h in TE + 0.1% Tween and precipitated in ethanol overnight.

First strand cDNA was circularized for 1 h with CircLigase (Epicentre), denatured for 10 min at 80°C, and relinearized with APE I (New England Biolabs). Relinearized DNA was PCR-amplified using the Phusion Hot Start II kit, according to the manufacturer’s instructions, with the primers oNTI200 (5′-CAAGCAGAAGACGGCATA-3′) and oNTI201 (5′-AATGATACGGCGACCACCGACAGGTTCAGAGTTCTACAGTCCGACG-3′). The PCR product was run on a 10% TBE gel and eluted as in the previous step. Libraries were sequenced on an Illumina hi-Seq2000 with sequencing primer 5′-CGACAGGTTCAGAGTTCTACAGTCCGACGATC-3′.

Transfection

For PPARγ siRNA knockdown, mature differentiated 3T3-L1 adipocytes were electroporated using Amaxa Cell Line L (program A-033, Lonza). One-quarter of a 10-cm dish of cells was treated with 200 pmol of siRNA (sequences in Supplemental Table S3) and plated onto one well of a 12-well dish. Cells were harvested 24–30 h post-transfection.

ChIP-seq

ChIP was performed as described previously (Steger et al. 2008). The following antibodies were used in this study: C/EBPα (sc-61x, Santa Cruz Biotechnology), C/EBPβ (sc-150x, Santa Cruz Biotechnology), FOSL2 (sc-13017x, Santa Cruz Biotechnology), JUND (sc-74x, Santa Cruz Biotechnology), ATF2 (sc-6233x, Santa Cruz Biotechnology), MED1 (A300-793A, Bethyl Laboratories), H3K27ac (ab4729, Abcam), CBP (sc-369x, Santa Cruz Biotechnology), p300 (sc-585x, Santa Cruz Biotechnology), and NCoR (rabbit polyclonal, generated in our laboratory). Primer sequences used for ChIP-qPCR are provided in Supplemental Table S2. ChIP DNA was prepared for sequencing according to the amplification protocol provided by Illumina. Next-generation sequencing of ChIP-seq libraries was performed by the Functional Genomics Core at the University of Pennsylvania, including preprocessing and alignment of the raw tags.

GRO-seq data processing

First, adapter and poly-A sequences were trimmed off from the raw tags, and remaining tags were converted into FASTA format before alignment. Trimmed tags were aligned to the mouse genome mm8 using Bowtie (Langmead et al. 2009) with the following options: Inputs were in FASTA format (-f), three mismatches were allowed for each tag (-v 3), and only uniquely mapped tags were retained (-m 1). All of the alignments were adjusted to 50 base pairs (bp) by extending toward the 3′ end if their size was shorter than that to make libraries comparable. For the visualization of GRO-seq data, we pooled two replicates for each time point and extended each tag to 150 bp to make smooth profiles. bedGraph files were generated first using genomeCoverageBed command in the BEDTools suite (Quinlan and Hall 2010) for plus and minus strands separately and were converted to bigwig format using the bedGraphToBigWig command in BLAT suite (Kent et al. 2010). All of the microarray and high-throughput sequencing data have been deposited in Gene Expression Omnibus under the accession number {"type":"entrez-geo","attrs":{"text":"GSE56747","term_id":"56747"}}GSE56747.

ChIP-seq data processing

All of the ChIP-seq tags for C/EBPα, C/EBPβ, FOSL2, ATF2, JUND, NCoR, and MED1 were aligned to mouse genome mm8 using Bowtie with options ”-k 1 -m 1–best –strata,” and all of the redundant tags were eliminated except one before downstream analysis. For PPARγ, we used two public PPARγ ChIP-seq data, {"type":"entrez-geo","attrs":{"text":"GSM340799","term_id":"340799"}}GSM340799 and {"type":"entrez-geo","attrs":{"text":"GSM678393","term_id":"678393"}}GSM678393 (Nielsen et al. 2008; Schmidt et al. 2011) from Gene Expression Omnibus, where all of the redundant tags were eliminated in each data set, and the remaining tags were pooled into a single data set. Peak calling was performed using the findPeaks command in Homer (Heinz et al. 2010). After initial calling, all of the peaks were resized to 200 bp; the 2 reads per million (RPM) cutoff was applied for C/EBPα, C/EBPβ, FOSL2, ATF2, JUND, and PPARγ to select strong peaks; and the 1 RPM cutoff was applied for MED1, CBP, p300, and NCoR because of their indirect genomic binding.

Gene transcription analysis

To determine rosi-induced regulation of gene transcription, we used normalized tag counts within a window between +0.5 kb and 12 kb of the TSS to capture acute changes and yet minimize the bias from paused signal at promoters. When counting tags, we ignored the tags aligned onto rRNA, small nucleolar RNA (snoRNA), small nuclear RNA (snRNA), or tRNA to avoid any false contribution to gene body GRO-seq signal due to these highly abundant elements. To identify rosi-regulated genes, we performed an exact test with edgeR package (Robinson et al. 2010) for each time point—10 min, 30 min, 1 h, and 3 h—using 0 min as a control. Genes were considered rosi-regulated if the false discover rate (FDR) was <0.05 at any time point, but those with low GRO-seq signal (<0.5 RPKM [reads per kilobase per million]) or poor gene body coverage (<70%) at every time point were discarded before downstream analysis. Temporal patterns of rosi-induced regulation were graphically presented by hierarchical clustering, where GRO-seq levels were RPKM-normalized and 1 − (correlation coefficient) was used as a distance measure between genes. The initial dendrogram was first created based on Ward’s criterion using the fastcluster R package (http://danifold.net/fastcluster.html) and then subjected to optimal leaf ordering (Bar-Joseph et al. 2001). In the visualization of clustering heat maps, RPKM values were converted into log₂(fold change over 0 min) for intuitive color representation of up-regulation or down-regulation. Gene ontology analysis was performed using DAVID (Huang et al. 2008), and top ranked GO FAT terms were presented.

We thank members of the Lazar laboratory for helpful discussions and reagents, and Klaus Kaestner for critically reading the manuscript. We also thank Leighton Core, John Lis, Minna Kaikkonen, and Chris Glass for help with the GRO-seq protocol; Fang Chen and Don Baldwin for performing the microarrays; and the Functional Genomics Core of the Penn Diabetes Center (DK19525) for deep sequencing. A.P. was supported by a GEN-AU mobility grant from the Austrian Ministry of Science and Research. This work was also supported by National Institutes of Health grants R01 DK49780 (to M.A.L.), R21DK098769 (to K.-J.W.), and 5T32GM008076 (to S.E.S.) and the JPB Foundation.