Identification of nuclear hormone receptor pathways causing insulin resistance by transcriptional and epigenomic analysis.

Insulin resistance is a cardinal feature of Type 2 diabetes (T2D) and a frequent complication of multiple clinical conditions, including obesity, ageing and steroid use, among others. How such a panoply of insults can result in a common phenotype is incompletely understood. Furthermore, very little is known about the transcriptional and epigenetic basis of this disorder, despite evidence that such pathways are likely to play a fundamental role. Here, we compare cell autonomous models of insulin resistance induced by the cytokine tumour necrosis factor-α or by the steroid dexamethasone to construct detailed transcriptional and epigenomic maps associated with cellular insulin resistance. These data predict that the glucocorticoid receptor and vitamin D receptor are common mediators of insulin resistance, which we validate using gain- and loss-of-function studies. These studies define a common transcriptional and epigenomic signature in cellular insulin resistance enabling the identification of pathogenic mechanisms.

Resistance to the metabolic actions of insulin is a critical factor in the pathogenesis of T2D. One puzzling feature of insulin resistance is that it arises in a wide variety of clinical scenarios, such as obesity, lipodystrophy, aging, acromegaly, pregnancy, inflammation, and steroid use, which possess few pathogenic features in common. We have previously generated transcriptional profiles from cell autonomous models of insulin resistance using cultured adipocytes treated with TNF and Dex[1]. TNF and Dex each promote insulin resistance in vivo and in vitro[2, 3], despite activating quite different molecular pathways; TNF is of course the prototypical pro-inflammatory agent while Dex is strongly anti-inflammatory. Our prior studies identified reactive oxygen species (ROS) as part of a common pathway that promotes insulin resistance, a finding corroborated by others[4, 5].

These and other studies of insulin resistance focused on membrane and cytosolic events in insulin resistance, with particular emphasis on the insulin signaling pathway, and, more recently, inflammatory pathways[6, 7]. Transcriptional and epigenetic events have not received nearly as much attention, despite a wealth of evidence that they are likely involved in the pathogenesis of insulin resistance. For example, the thiazolidinediones comprise a class of clinically effective insulin sensitizers; they work by binding and activating the transcription factor PPARγ[8]. Additionally, the phenomenon of fetal programming demonstrates that acquired metabolic traits like obesity and insulin resistance can be transmitted transgenerationally, an effect amply demonstrated in rodent studies as well as in ‘experiments’ of human history, such as the Dutch Hunger Winter[9,10]. Additionally, drugs that target chromatin-modifying enzymes, such as the histone deacetylase inhibitor valproic acid, reduce insulin sensitivity[11]. Finally, targeted ablation of genes encoding certain histone modifying enzymes (such as Jhdm2a) causes metabolic dysfunction in mice, including insulin resistance[12, 13].

We sought to understand how the epigenome is altered in states of insulin resistance, using established cellular models. We have shown previously that epigenomic information derived from such models can be used to generate testable hypotheses about the transcriptional pathways involved in a biological process of interest (for example, in the case of adipocyte differentiation)[14]. Here, we utilize comparative models of cellular insulin resistance to define a core set of epigenomic and transcriptional changes associated with this disorder, and use this information to predict and validate two causal transcriptional pathways.

Results

Establishment of cell autonomous models of insulin resistance

Murine 3T3-L1 adipocytes were treated separately with dexamethasone (Dex; 20nM) or tumor necrosis factor-α (TNF; 10ng/ml), and basal and insulin-stimulated glucose uptake were assessed at multiple time points (Supplementary Fig. 1a). These agents cause insulin resistance in rodents and humans[15-18]; in our model, they elicit an equivalent and significant reduction in insulin-stimulated glucose uptake by approximately 60% without affecting basal uptake (Fig. 1a), consistent with the degree of insulin resistance seen in T2D[19]. Time course analysis reveals a gradual diminution of insulin sensitivity that reaches a maximum after six days of treatment (Fig. 1b). Importantly, no dedifferentiation of the cells was noted, using either oil red O staining of neutral lipids or expression of genes that mark the mature state (Supplementary Figs. 1b, c, e, f). Furthermore, we compared the global transcriptional profiles induced by Dex and TNF (described below) with those associated with 3T3-L1 adipogenesis[14]. A heat map of the correlation coefficients obtained by comparing each time point in differentiation to each time point after drug treatment reveals that, while gene expression profiles change with Dex and TNF, they do not induce a state that corresponds better to undifferentiated time-points vs. control cells (Supplementary Fig. 1d). This is perhaps not surprising given that the majority of papers showing TNF-induced adipocyte de-differentiation employ higher doses than we use here[20, 21].

Figure 1

Transcriptional and epigenomic profiling of the insulin resistance model

a, Basal and insulin-stimulated glucose uptake in cells treated with Dex and TNF for 6 days (n=3 independent experiments, each point represents the mean of 6 independent replicates p<0.05 by Student’s t-test, mean±SEM). b, Adipocytes were treated with Dex or TNF for the indicated lengths of time and the development of insulin resistance was monitored by glucose uptake assay; results are expressed as percent of glucose uptake remaining compared to untreated (control) (n=3 independent experiments, each point represents the mean of 6 independent replicates p<0.05 by Student’s t-test, mean±SEM). c, Venn diagram depicting the number of up-regulated genes in response to Dex and TNF treatment. These numbers represent a ‘unity’ set in which all time points were combined. d, Fold change (log2 scale) of genes from c in the white adipose tissue of obese vs. lean mice[22], n = number of genes in each group (invariant = 160, Dex UP =357, TNF_UP = 512, Overlap_UP = 110) p<0.05 by Student’s t-test, mean±SD). e, H3K27ac peak density as a function of distance to the TSS of the nearest gene. f, Distribution of all H3K27ac peaks in untreated cells (top), and the distribution of H3K27ac peaks that change with Dex (middle) and TNF are shown (bottom). g, Sequence conservation score of DNA under distal H3K27ac peaks (>10 kb from TSS) as a function of distance from peak center, for invariant and changed peaks.

Transcriptional profiling of cellular insulin resistance

To comprehensively assess gene expression changes caused by Dex and TNF in a time-dependent manner, we profiled cells at early (2 hours), intermediate (24 hours), and late (6 days) points in the development of insulin resistance. We were able to identify groups of genes that showed Dex-selective or TNF-selective induction or repression, many of which could be predicted based on the known functions of these agents (Supplementary Fig. 1g). For example, TNF induced genes associated with inflammation while these same genes were repressed by Dex (Supplementary Tables 1 and 2). We identified 1016 genes up-regulated by Dex at any time point, and 1346 genes up-regulated by TNF; 271 of these were induced by both agents (Fig. 1c). We confirmed expression changes for a subset of these by QPCR (10/10 tested). Conversely, 1044 and 960 genes were down-regulated by Dex and TNF, respectively, with 341 in the intersecting set (Supplementary Fig. 1h; Supplementary Table 3).

We next assessed how well this system models insulin resistance associated with obesity. Comparison to publicly available gene expression data from the white adipose tissue of chow and high-fat fed C57Bl/6 and BTBR mice[22] revealed that up-regulated genes, as a whole, correlated better with these in vivo models than down-regulated genes. TNF-sensitive up-regulated genes were relatively well correlated with obesity (Supplementary Figs 2a–g), and genes that are concordantly up-regulated by both Dex and TNF were, if anything, even more highly correlated in the C57Bl/6 obesity model (Fig. 1d). The fold change of these genes in the obesity models is significant yet small, an expected finding given that the obesity data are derived from unfractionated adipose tissue containing a high percentage of inflammatory immune cells, and because Dex and TNF exert a wide array of gene expression changes, only some of which are related to insulin resistance.

Epigenomic analysis of histone modifications in insulin resistance

Because transcriptional changes are almost always associated with alterations in chromatin state, we assessed whether treatment with Dex or TNF could alter bulk histone modification using reverse-phase HPLC and mass spectrometry. The majority of histone modifications tested were not altered (Supplementary Fig. 3a), and none displayed concordant changes with both agents. Consistent with this finding, ChIP-Seq analysis using antibodies specific for histone modifications that identify promoters (H3K4me3), poised and active enhancers (H3K4me1, H3K27ac), repressed domains (H3K27me3), and transcriptional units (H3K36me3 and H3K79me2) did not reveal any significant changes in the overall number of peaks for any specific modification (Supplementary Fig. 3b).

To assess overall changes in chromatin, we annotated the epigenome in response to Dex and TNF treatment, applying a multivariate hidden Markov model that uses combinatorial patterns of histone marks to distinguish chromatin state (Supplementary Figs 3c–e)[23, 24]. We distinguished five broad classes of chromatin states, which we refer to as promoter, enhancer, transcribed, repressed and inactive states. We saw no significant changes in the total number of promoters (~13,000) or enhancers (~70,000) with either treatment. Chromatin state mapping using histone modifications has been widely used to compare different cell types[24, 25] or cells during differentiation[14], showing dynamic changes in promoter and enhancer states that are established via lineage-specific transcription factors. Here we show that promoters are largely invariant within this differentiated cell type, but both Dex and TNF induce numerous changes in enhancers (Supplementary Fig. 3f). When we compare these to the chromatin state changes computed from our preexisting adipocyte differentiation data[14] using identical state modeling parameters, we observe fewer changes at both promoters and enhancers, indicating that cellular differentiation involves a more radical restructuring of chromatin than even a prolonged pharmacological stimulus.

Focusing on H3K27ac as a marker of promoter and enhancer activity, we found that changes induced by Dex and TNF occur preferentially at distal enhancers. Across all conditions, 20% of H3K27ac peaks are found within ±2kb of the transcriptional start site (TSS) of genes and 33% are located in intergenic regions. After Dex or TNF exposure, however, only 5% of changed H3K27ac peaks localize to the TSS, while 50% localize to intergenic sites, with a disproportionate number of changes occurring >100kb from the nearest TSS (Figs. 1e, f). Increased conservation at H3K27ac peaks that change with Dex and TNF relative to invariant H3K27ac regions supports a functional role for these loci (p<10−13; Fig. 1g).

H3K27ac is a marker of active enhancers, and its degree of enrichment has been correlated to gene expression[26]. Therefore, we next asked whether induced H3K27ac peaks were located near genes that were induced by Dex and TNF. In our transcriptional analysis, we found 271 genes whose expression was up-regulated by both Dex and TNF. We delineated all up-regulated H3K27ac peaks within ±200kb from the TSS of the affected gene as candidate enhancers driving the transcriptional changes; 249/271 genes had at least one Dex and/or TNF up-regulated H3K27ac peak within this interval (Supplementary Fig. 4a; Supplementary Table 4). We identified four patterns among these gene-peak sets (Supplementary Figs. 4a–e). Some genes had either a Dex up-regulated (n=22), or a TNF up-regulated peak (n=80), but not both, within the tested region. Presumably, since the expression of these genes was coordinately regulated by both agents, the absent Dex or TNF-induced peak(s) either lie outside of the ±200kb window, or the changes were below our detection limits. One hundred and forty-seven genes, however, had both a Dex- and TNF-induced H3K27ac peak within the ±200kb interval. Of these, 67% (n=99) had differentially modified H3K27ac peaks at different genomic positions (e.g. Sbf2, Mfhas, and Ifngr1; Figs. 2b, d and Supplementary Fig. 4d), while 33% (n=48) showed overlapping peaks induced by both Dex and TNF (e.g. Tmem176a/b and Fam46b; Fig. 3a and Supplementary Fig. 4e). These two patterns suggest different modes of regulation, with the former likely the result of Dex and TNF working through different transcriptional pathways that happen to converge on the same gene from two distinct cis-regulatory elements, while the latter is consistent with the possibility of a common transcriptional mediator downstream of both Dex and TNF.

Figure 2

Motif analysis of co-induced H3K27ac enhancer peaks

a, Venn diagram depicting H3K27ac peaks up-regulated by Dex and TNF within ±200kb of the TSS of an induced gene. Also depicted are significantly enriched binding sites with TRANSFAC position weight matrix results recovered by motif analysis of ‘Dex only’, ‘TNF-only’, and ‘Dex-TNF overlapping’ peaks (with p values determine by Fisher’s Exact test). b, Tracks from an up-regulated gene (Sbf2) showing distinct Dex- and TNF-sensitive H3K27ac peaks. c, ChIP-PCR results from the Dex- and TNF-induced peaks depicted in b (data represent mean of n = 3 dishes, data from 2 additional independent experiments shown in Source Data Table), using anti-GR, anti-p65, or IgG. d, Tracks from an up-regulated gene (Mfhas1) showing ‘Dex-TNF-nonoverlapping’ H3K27ac peaks. e, ChIP-PCR results from the ‘Dex-only’ and ‘TNF-only’ peaks depicted in d (data represent mean of n = 3 dishes, data from 2 additional independent experiments shown in Source Data Table).

Figure 3

Dex and TNF cause glucocorticoid receptor to bind to predicted motifs

a, Tracks from up-regulated genes (Tmem176a/b) showing a ‘Dex-TNF overlapping’ peak. b, ChIP-PCR results from the ‘Dex-TNF overlapping’ peak depicted in a (data represent mean of n = 3 dishes, data from 2 additional independent experiments shown in Source Data Table), using anti-GR, anti-p65, or IgG in the 3T3-L1 model. c, Primary SVF cells differentiated ex vivo into adipocytes were analyzed as in b (data represent mean of n = 3 dishes from pooled cells of three individual mice, data from 1 additional independent experiment shown in Source Data Table). d, Adipose tissue from chow-fed (n=10) high fat fed (n=5) mice were harvested, pooled, and immunoprecipitated three separate times and analyzed as in b (data from 1 additional independent experiment shown in Source Data Table). e, ChIP-PCR results from a ‘Dex-TNF overlapping’ peak in the Lcn2 gene (data represent mean of n = 3 dishes, data from 2 additional independent experiments shown in Source Data Table), using anti-GR or IgG in the 3T3-L1 model. Also shown is the effect of knocking down GR with a specific shRNA. f, Primary SVF cells differentiated ex vivo into adipocytes were analyzed as in b (data represent mean of n = 3 dishes from pooled cells of three individual mice, data from 1 additional independent experiment shown in Source Data Table). g, Adipose tissue from chow-fed (n=10) high fat fed (n=5) mice were harvested, pooled, and immunoprecipitated three separate times and analyzed as in d (data from 1 additional independent experiment shown in Source Data Table). h, ChIP-PCR results from a ‘Dex-TNF overlapping’ peak in the Serpina3n gene (data represent mean of n = 3 dishes, data from 1 additional independent experiment shown in Source Data Table), using anti-GR or IgG in the 3T3-L1 model. Also shown is the effect of knocking down GR with a specific shRNA. i, Primary SVF cells differentiated ex vivo into adipocytes were analyzed as in b (data represent mean of n = 3 dishes from pooled cells of three individual mice, data from 1 additional independent experiment shown in Source Data Table). j, Adipose tissue from chow-fed (n=10) high fat fed (n=5) mice were harvested, pooled, and immunoprecipitated three separate times and analyzed as in e (data from 1 additional independent experiment shown in Source Data Table).

TNF induces glucocorticoid receptor (GR) binding to regions where H2K27ac is induced

To identify downstream transcriptional effectors of Dex and TNF, we identified overrepresented motifs in up-regulated H3K27ac peaks. Not surprisingly, the top motif identified in Dex-induced peaks was a GR binding site (p<10−52) while the top motif from TNF-induced peaks corresponded to an NF-κB site (p<10−51) (Supplementary Table 5). We next focused on the genomic space comprising ±200kb from the TSS of genes coordinately induced by Dex and TNF. Unlike the gene-centric analysis described earlier, here we quantified the number of H3K27ac peaks induced by these agents. In all, we identified 326 Dex-induced and 598 TNF-induced peaks (Fig. 2a; Supplementary Table 4) Of these, 53 peaks were co-induced by both agents. Analysis of the ‘Dex-only’ peaks recovered two motifs corresponding to binding sequences for the androgen receptor (AR) and the glucocorticoid receptor (GR) (Fig. 2a). The GR and AR are highly similar proteins, and bind nearly identical motifs[27]. The top-scoring motif for ‘TNF-only’ peaks corresponded to the binding site for NF-κB. These results suggest that Dex and TNF induce the majority of epigenomic changes at active enhancers through their classic transcriptional effector pathways, which we confirmed by ChIP-PCR (Figs. 2b–e). For Dex-TNF overlapping peaks, however, two of the most enriched motifs suggested binding sites for the GR and the vitamin D receptor (VDR). Random permutation of the peak labels (‘Dex-only’, ‘TNF-only’, and ‘Dex-TNF-overlapping’) confirmed that the GR and VDR motifs remained significant (p<0.05).

We were surprised to find the GR binding site as the most enriched motif within shared up-regulated peaks, as this suggested the counter-intuitive notion that TNF exerts some of its actions on insulin resistance via the GR. We addressed this by performing ChIP-PCR of GR in 3T3-L1 adipocytes; at ‘Dex-TNF-overlapping’ peaks, TNF and Dex both induce GR binding at the site of the GR motif (Figs. 3a, b, e, h). To establish whether this phenomenon is generalizable beyond 3T3-L1 cells, we isolated primary stromal-vascular cells from the adipose depots of C57Bl/6 mice and differentiated them into adipocytes ex vivo; Dex and TNF each induced GR binding to these locations (Figs. 3c, f, i). We also found significantly enhanced binding of the GR to these locations in the adipose tissue of obese, insulin resistant mice (Figs. 3d, g, j).

The GR mediates TNF-induced insulin resistance

To test whether the GR has a functional role in insulin resistance, we knocked it down in mature 3T3-L1 adipocytes using lentiviral delivery of four distinct shRNAs (Supplementary Fig. 5a). As expected, these cells were completely protected from Dex-induced insulin resistance (Fig. 4a). Somewhat more surprisingly, these shRNAs also conferred significant (though not total), protection against TNF-induced insulin resistance. Conversely, knockdown of p65 (encoded by Rela) did not block the ability of Dex to cause insulin resistance (Supplementary Fig. 5b, c); the effect of p65 knockdown on TNF-induced insulin resistance was not testable because cells bearing a p65 shRNA undergo apoptosis after TNF treatment. The ability of shGR to protect against both Dex- and TNF-mediated insulin resistance was confirmed in the ex vivo model (Fig. 4b).

Figure 4

TNF induces insulin resistance via ligand-independent activation of the glucocorticoid receptor

a, Several shRNAs against GR (vs. scrambled shRNA: shScr) were delivered to mature adipocytes via lentivirus, and cells were then treated with Dex or TNF and assessed for insulin-stimulated glucose uptake. Shown is the percent of insulin-stimulated glucose uptake rescued by GR knockdown (data represent mean of n = 3 dishes, data from 1 additional independent experiment shown in Source Data Table). b, Primary SVF cells differentiated ex vivo into adipocytes were analyzed as in a (data represent mean of n = 6 dishes from pooled cells of three individual mice, data from 1 additional independent experiment shown in Source Data Table). c, 3T3-L1 cells were transduced with lentivirus expressing Flag-GR and subsequently exposed to Dex or TNF for the indicated length of time. Cytosolic and nuclear cellular fractions were blotted with the indicated antibody. At right, the effect of TNF on GR is shown in whole cell extracts from the same cells. d, Insulin-stimulated glucose uptake was assessed in 3T3-L1 adipocytes expressing an shRNA against GR (vs. control; Scr) and then repleted with wild-type or mutant GR alleles (data represent mean of n = 3 dishes, data from 1 additional independent experiment shown in Source Data Table). e, Percentage restoration of drug-mediated insulin resistance from d is presented by normalizing the results of the mutant GR alleles to that of wild-type GR. The difference in percentage restoration between wild-type and each mutant (Δ) represents a ligand-dependent effect while the Δ difference between GR knock-down and each mutant GR allele represents a ligand-independent effect. (data represent mean of n = 3 dishes, data from 1 additional independent experiment shown in Source Data Table). f, RU-486, a pharmacological antagonist of GR, was administered at the indicated doses to fully differentiated adipocytes with or without Dex or TNF for 6 days and then assessed for insulin-stimulated glucose uptake (data represent mean of n = 6 dishes, data from 2 additional independent experiments shown in Source Data Table).

In the unliganded state, GR is primarily localized to the cytosol, sequestered by Hsp90 and other chaperone proteins until corticosteroid binding enables nuclear translocation and transcriptional activity[28]. We confirmed this mode of regulation in our system, with Dex causing simultaneous cytoplasmic depletion and nuclear enrichment of GR (Fig. 4c). TNF also caused GR translocation to the nucleus. This was not accompanied by depletion of cytoplasmic GR; on the contrary, TNF dramatically increased GR protein.

TNF-induced translocation of the GR raised the possibility that intracellular ligand levels could be affected. In monocytes, for example, TNF has been shown to induce the expression of 11β-hydroxysteroid dehydrogenase (11β-HSD), the enzyme responsible for the conversion of inactive to active corticosteroids[29]. This proved not to be the case in adipocytes, as Dex and TNF reduced 11β-HSD expression (Supplementary Fig. 5d). Furthermore, even high doses of the 11β-HSD inhibitor carbenoxolone did not block the ability of TNF to cause insulin resistance (Supplementary Figs. 5e, f). We next assessed the ability of mutant GR alleles that lack the ability to bind ligand to mediate TNF-induced insulin resistance. As expected, two naturally occurring human variants (V729I and D641V) and a rat mutant GR (N768) both prevented Dex from causing insulin resistance. These alleles blocked approximately half of TNF-induced insulin resistance (Figs. 4d, e), suggesting that part of the effect of TNF on the GR may be ligand-dependent. Consistent with this, the glucocorticoid antagonist RU-486 blocks some of the TNF-mediated reduction in glucose uptake (Fig. 4f). Taken together, these studies indicate that TNF causes insulin resistance via both ligand-dependent and ligand-independent GR activation.

TNF induces genome-wide binding of the GR to a subset of Dex-activated sites

We next assessed the prevalence of TNF-induced GR binding in adipocytes by performing ChIP-Seq in cells treated with Dex or TNF. While only 11 peaks were called in the absence of either agent, Dex treatment was associated with the appearance of 3,872 peaks, comparable to what others have reported in this cell type[30] (Fig. 5a). TNF induced only 403 peaks, most of which were represented in the Dex-activated set, suggesting that TNF causes the GR to associate with a subset of the sites that would normally be engaged in the presence of glucocorticoid. Consistent with this, de novo motif finding at TNF-induced GR binding sites revealed a GRE identical to the Dex-induced GRE (Fig. 5b). In general, TNF-induced GR binding sites were slightly more likely to occur at promoters (Supplementary Fig. 6a), and were more conserved than Dex-induced GR binding sites (p=0.01; Supplementary Fig. 6b). Common GR binding sites were more highly conserved than either Dex-induced (p<10−8) or TNF-induced binding sites (p<10−3). Both Dex- and TNF-induced GR binding sites were correlated with changes in the expression of the nearest gene (Fig. 5c; Supplementary Fig. 6c).

Figure 5

TNF induces GR binding at a subset of Dex-induced sites

a, Venn diagram of genome-wide GR binding sites in response to Dex and TNF. These numbers represent a unity set combining time points. b, Top de novo motifs recovered from Dex- and TNF-induced GR binding sites. c, Heat map depicting expression changes (expressed as the absolute value of the fold change) of the nearest gene to labeled sets of GR peaks. White boxes were not significant when compared to changes from a random set of H3K27ac peaks. d, GR and H3K27ac tracks from the Thada gene showing four GR binding sites (A–D).

Many loci show complex regulation by GR in these cells. One example is the GWAS-associated T2D gene Thada, where Dex induces GR binding at four separate upstream sites, some in a time-dependent manner (Fig. 5d), TNF induces GR binding and H3K27ac accumulation at two of these sites, clustered near a region that is orthologous to the site of the lead human diabetes SNP[31].

Vdr expression is induced by Dex and TNF via activation of the GR

The VDR binding motif was also enriched in H3K27ac peaks induced concordantly by Dex and TNF. To determine if VDR plays a role in Dex- and TNF-mediated insulin resistance, we first measured its expression in adipocytes treated with these agents. Vdr mRNA levels increased after treating with Dex and TNF (Fig. 6a). Examination of the Vdr locus revealed a TNF-induced H3K27ac peak approximately 17 kb upstream of the TSS, and an intronic Dex-induced peak (Fig. 6b). Interestingly, both of these peaks contain GR motifs, and ChIP-PCR confirmed GR (but not NF-κB) binding at the ‘TNF-only’ site in response to TNF (Fig. 6c). This was true in the ex vivo adipocyte model as well, as TNF and Dex induced GR binding to both of these sites (Fig. 6d). In vivo, we were able to detect GR binding to the intronic locus (Fig. 6e). Finally, knockdown of GR blocked induction of Vdr gene expression by TNF and Dex (Fig. 6f).

Figure 6

VDR is a Dex- and TNF-inducible transcription factor

a, Vdr mRNA levels were assessed in cells treated with Dex or TNF for 6 days (n=3 independent experiments, each point represents the mean of 3 independent replicates p<0.05 by Student’s t-test, mean±SEM). b, Tracks from the Vdr gene locus showing distinct Dex- and TNF-sensitive H3K27ac peaks. Both peaks contain a GR motif. c, ChIP-PCR results from the ‘Dex-only’ and ‘TNF-only’ peaks depicted in b (data represent mean of n = 3 dishes, data from 2 additional independent experiments shown in Source Data Table), using anti-GR, anti-p65, or IgG in 3T3-L1 adipocytes. d, Primary SVF cells differentiated ex vivo into adipocytes were analyzed as in b (data represent mean of n= 3 dishes from pooled cells of three individual mice, data from 1 additional independent experiment shown in Source Data Table). e, Adipose tissue from chow-fed (n=10) high fat fed (n=5) mice were harvested, pooled, and immunoprecipitated three separate times and analyzed as in d (data from 1 additional independent experiment shown in Source Data Table). f, Expression of Vdr mRNA in cells transduced with lentivirus carrying shRNA directed against VDR, p65, or GR (vs. scrambled control: shScr) and then treated with Dex or TNF for 6 days (data represent mean of n = 3 dishes, data from 2 additional independent experiments shown in Source Data Table). g, Fully differentiated adipocytes were pre-treated for 24 hrs with chemical inhibitors of various kinases, and then treated with TNF for 4 additional days. Cells were assessed for relative expression of Vdr (g, h) and insulin-stimulated glucose uptake (i, j) (data represent mean of n = 3 dishes, data from 1 additional independent experiment shown in Source Data Table).

TNF has been proposed to induce insulin resistance via several different signaling mediators, including p38, ERK, JNK, and CDKs[32-36]. In our system, p38 and CDK2 inhibitors prevented TNF-induced insulin resistance, but only p38 inhibition reduced TNF-dependent Vdr expression (Figs. 6g–j).

The VDR is a common mediator of insulin resistance

In order to demonstrate that VDR binds to the predicted genomic loci in response to Dex and TNF, we looked at several ‘Dex-TNF overlapping’ H3K27ac peaks with a VDR motif using ChIP-PCR analysis. This demonstrated Dex- and TNF-dependent VDR binding (Figs. 7a, b; Supplementary Figs. 7a–f). This was recapitulated perfectly in both the ex vivo and in vivo models (Fig. 7b). We also noted that RXR binds at the same sites as VDR, indicating that VDR binds as the classic heterodimer (Supplementary Fig. 7g).

Figure 7

VDR binds to predicted motifs and causes insulin resistance

a, Tracks at a gene (Colq) up-regulated by Dex and TNF showing a Dex-TNF-overlapping H3K27ac peak that contains a VDR motif. b, ChIP-PCR results from the ‘Dex-TNF overlapping’ peak depicted in a, using anti-VDR or IgG in the 3T3-L1 (data represent mean of n = 3 dishes, data from 2 additional independent experiments shown in Source Data Table), ex vivo (data represent mean of n = 3 dishes from pooled cells of 3 individual mice, data from 1 additional independent experiment shown in Source Data Table), and in vivo (data from 1 additional independent experiment shown in Source Data Table) models. c, Glucose uptake was assessed in cells overexpressing VDR 4 days after lentiviral transduction (data represent mean of n = 6 dishes, data from 2 additional independent experiments shown in Source Data Table). d, Several shRNAs against VDR (vs. scrambled shRNA: shScr) were delivered to mature adipocytes (Day 6) via lentivirus, and cells were then treated with Dex and TNF for 4 days and assessed for insulin-stimulated glucose uptake. Shown is the percent of insulin-stimulated glucose uptake rescued by VDR knockdown in the 3T3-L1 (left; data represent mean of n = 6 dishes, data from 2 additional independent experiments shown in Source Data Table) and ex vivo (right; data represent mean of n = 3 dishes from pooled cells of 3 individual mice, data from 1 additional independent experiment shown in Source Data Table) models. e, Vdr expression in adipose tissue from chow and high fat-fed C57Bl/6 mice (left; n=6 mice p<0.05, Student’s t-test, mean±SEM) and in adipose tissue from ob/+ and ob/ob mice in the presence and absence of rosiglitazone treatment (right; n=7 for ob/+, n=8 for ob/+ plus Rosi, n=7 for ob/ob with and without Rosi, p<0.05, Student’s t-test, mean±SEM). f, 1,25(OH)2-vitamin D3 (100 nM) or vehicle was added to fully differentiated adipocytes with Dex or TNF for 6 days and then assessed for insulin-stimulated glucose uptake (data represent mean of n = 3 dishes, data from 1 additional independent experiment shown in Source Data Table). g, Glucose uptake was assessed in cells overexpressing VDR with and without the addition of exogenous 1,25(OH)2-vitamin D3 (100 nM; data represent mean of n = 6 dishes, data from 1 additional independent experiment shown in Source Data Table). h, expression of the VDR-sensitive Txnip and Dusp10 genes in 3T3-L1 adipocytes treated with 100 nM 1,25(OH)2-vitamin D3 (data represent mean of n = 3 dishes, data from 1 additional independent experiment shown in Source Data Table).

To test the effect of VDR directly, we overexpressed it in mature adipocytes and found reduced insulin-stimulated glucose uptake to a degree comparable to the effects of Dex and TNF (Fig. 7c; Supplementary Fig. 7h). Conversely, shRNA-mediated knockdown of VDR prevented insulin resistance in both the 3T3-L1 and ex vivo models, restoring between 20–80% of the insulin sensitivity lost to Dex and TNF (Fig. 7d; Supplementary Fig. 7i).

Vdr expression in the white adipose tissue of insulin resistant high-fat fed mice and genetically obese ob/ob mice was increased (Fig. 7e). Furthermore, treatment of the latter with the insulin-sensitizing agent rosiglitazone reduced Vdr expression. In contrast, exogenous vitamin D did not prevent or worsen insulin resistance in the presence of TNF, Dex, or overexpressed VDR (Figs 7f–h).

Our prior work demonstrated that increased oxidative stress represents a pathogenic node in the development of insulin resistance[1]. In the present study, we found that the original ROS gene set was again significantly and coordinately regulated by Dex and TNF (p<0.05 for all Dex time points and two out of three TNF time points). Furthermore, ROS-related genes were significantly enriched for TNF- and Dex-sensitive H3K27ac peaks (p<0.0002). As before, chemical antioxidants block TNF- and Dex-induced insulin resistance (Supplementary Fig. 7j). Interestingly, these agents prevent maximal induction of Vdr expression, and also block the development of insulin resistance in response to VDR overexpression (Supplementary Figs. 7k, l). Thus, oxidative stress acts both upstream and downstream of VDR in the development of insulin resistance.

GR and VDR target genes represent insulin resistance genes

Finally, we sought to identify insulin resistance promoting target genes of GR and VDR. We focused on a restricted set of genes with (a) expression coordinately up-regulated by Dex and TNF, (b) expression elevated in the adipose tissue of insulin resistant rodent models (with reversibility after rosiglitazone treatment), (c) at least one shared up-regulated H3K27ac peak within ±200kb of the TSS, and (d) a VDR motif within that shared peak. Finally, we required that the gene no longer be induced by Dex or TNF if GR and VDR were knocked down (Figs. 8a, b; Supplementary Figs. 8a–d). Four genes or gene clusters met all of these criteria (Tmem176a/b, Serpina3m/n, Colq, and Lcn2). We have shown previously that Lcn2 promotes insulin resistance in 3T3-L1 cells[37], and Lcn2 null mice are at least partially protected from insulin resistance[38, 39]. Similarly, the other three genes were all individually able to reduce insulin-stimulated glucose uptake following overexpression (Fig. 8c). Taken together, these studies establish a transcriptional network that causes insulin resistance, in which the nuclear receptors GR and VDR act as integrating mediators of disparate perturbations (Fig. 8d).

Figure 8

GR and VDR target genes are elevated in obesity and promote insulin resistance

a, Expression of coordinately up-regulated genes was measured in white adipose tissue samples from ob/+ and ob/ob mice, treated with vehicle or rosiglitazone for 6 weeks (n=7 for ob/+, n=8 for ob/+ plus Rosi, n=7 for ob/ob with and without Rosi, p<0.05, Student’s t-test, mean±SEM). b, Expression of Tmem176a, Serpina3n, and Colq mRNA in cells transduced with lentivirus carrying shRNA directed against VDR or GR (vs. scrambled control: shScr) and then treated with Dex (D) or TNF (T) for 6 days (data represent mean of n = 3 dishes, data from 1 additional independent experiment shown in Source Data Table). c, Glucose uptake was assessed in cells overexpressing Tmem176a, Serpina3n, and Colq (data represent mean of n = 6 dishes, data from 1 additional independent experiment shown in Source Data Table). d, Proposed model of the transcriptional basis of Dex- and TNF-induced insulin resistance. The dashed line indicates that TNF likely operates through additional mechanisms.

Discussion

Insulin resistance is a major feature of obesity and T2D, and is associated with a wide variety of clinical and experimental conditions. Several mechanisms have been proposed to explain insulin resistance, mostly involving transient cell signaling events that include inflammation, ER stress, and ROS generation[1, 5, 7, 40]. There has been little attention paid to transcriptional and epigenetic events, despite epidemiological and molecular data that suggests that such pathways might be important. As one example, the long time course required for the development of insulin resistance in many models (including ours) is not consistent with complete dependence on signaling events that take place within seconds to minutes.

Here we utilize a strategy analogous to one we employed to identify transcriptional regulators of adipogenesis[14]. Global identification of enhancer marks and subsequent motif analysis enables the prediction of transcription factors that would be hard to identify using other means. In this study, we identify two nuclear receptor-mediated pathways that cause insulin resistance in response to both Dex and TNF, with relevance to more common forms of insulin resistance, as seen in obesity. Specifically, the first pathway involves the GR, which has been implicated in Dex-induced insulin resistance[41, 42]. Our epigenomic data suggested that TNF, a prototypical pro-inflammatory agent, might also activate the GR, a prediction that we validated.

Recent papers have suggested that TNF may enhance the ability of glucocorticoid to induce the GR to bind to its gene targets in nonadipose cells, or even to bind new targets[43, 44]. Here we demonstrate that TNF promotes GR nuclear translocation and binding to a select set of genomic elements without a requirement for exogenous or endogenous corticosteroids. How these loci are targeted by TNF is unknown, but may involve interactions with other nearby proteins. The small number of TNF-specific GR binding sites contain an indeterminate motif that does not resemble a GRE (Fig. 5b), suggesting either noncanonical direct GR binding or ‘tethering’ through other factors[45]. Consistent with the latter, TNF-specific sites tend to be weaker than overlapping sites.

NF-κB acts as the transcriptional effector for many TNF actions. Although many studies have indirectly implicated NF-κB in insulin resistance by showing the importance of the upstream kinase IKKβ[46-48] the interpretation of these studies is complicated by the fact that IKKβ can cause insulin resistance via NF-κB-independent pathways[49]. Interestingly, mice with transgenic overexpression of p65 in adipose tissue show reduced body weight and enhanced insulin sensitivity, despite a pro-inflammatory gene expression pattern; a similar phenotype is also seen in mice lacking the regulatory NF-κB subunit p50[50, 51]. This suggests that TNF promotes insulin resistance via a different transcriptional effector; we propose a role for the GR in this capacity. Additionally, the GR shares a motif with the androgen and progesterone receptors, both of which have also been implicated in metabolic pathology[52, 53]. A role for these factors as downstream mediators of TNF action should be explored as well.

Our discovery of a role for the VDR in promoting insulin resistance in adipocytes was surprising, given evidence that serum levels of vitamin D are positively associated with insulin sensitivity[54, 55]. Some have ascribed this association to sequestration of fat-soluble vitamin D in the adipose tissue of obese subjects, or to reduced sun exposure in the physically inactive. Regardless, there is precedent for ligand-receptor pairs to show different activities[56-58]; ChIP-Seq studies of VDR in human immune cells support the idea that liganded and unliganded receptor can bind to distinct genomic loci[59, 60]. Notably, in our hands, the addition of exogenous vitamin D did not affect the development of insulin resistance after Dex, TNF, or overexpression of VDR. Interestingly, transgenic mice that overexpress VDR in adipocytes show impaired glucose tolerance, albeit in the context of increased adiposity, making it hard to determine whether the effect on glucose homeostasis is direct in that system[61].

Our results demonstrate the power of epigenomic analysis as a discovery tool for biological pathways operating in physiology and disease. Specifically, these studies identified transcription factors amenable to pharmacological manipulation, and provide insight into a network of downstream genes that expand our current knowledge of the pathophysiology of this highly prevalent disorder.

Materials and Methods

Cell culture

3T3-L1 preadipocytes were obtained from ATCC and maintained and differentiated as described[1]. To generate lentivirus particles, lentiviral constructs were co-transfected with pM2DG and psPAX-expressing plasmids into 293T cells. After 48 h, virus-containing supernatant was collected, filtered through 0.45μm filters, and added to mature 3T3-L1 adipocytes for 24h along with 8μg/ml Polybrene. Transduction efficiency was determined by comparison to cells transduced with a GFP expressing lentivirus in parallel. For the ex vivo system, subcutaneous adipose tissue from wild type C57Bl/6 mice was fractionated with digestion buffer (10mg/ml collagenase D, 2.4 units Dispase II, 10mM CaCl2 in PBS). Cells from the stromal-vascular fraction (SVF) were plated in culture and differentiated as described[1]. N-acetyl-L-cysteine (NAC), Mn(III) tetrakis (4-benzoic acid) porphyrin (MnTBAP), carbenoxolone, 1,25-(OH)2-vitamin D3, RU-486, inhibitors of p38MAPK (SB203580), ERK (PD98059), JNK (SP600125), and CDK2 (GW8510) were purchased from Sigma.

Animals

For rosiglitazone studies, female mice ob/ob and ob/+ were treated starting at 8 weeks of age with rosiglitazone by gavage (10mg/kg body weight) in 0.5% carboxymethylcellulose versus 0.5% carboxymethylcellulose vehicle control daily for 6 weeks (ob/+ vehicle n=7, ob/+ rosi n=8, ob/ob vehicle n=7, ob/ob rosi n=7). For chow and high fat feeding studies, male C57Bl/6J mice were put on diet beginning at 8 wks of age and continued for three months (n=8 per dietary condition). Samples were harvested from the perigonadal fat pad. All animal work was approved by the BIDMC IACUC.

RNA extraction and quantitative PCR (qPCR)

Total RNA was extracted from cells or tissues using TRIzol reagent according to the manufacturer’s instructions. cDNA was reverse-transcribed from 1μg of RNA using the RETROscript first strand synthesis kit (Ambion). QPCR was performed with SYBR Green qPCR Master Mix (Applied Biosystems) using a 7900HT Fast Real-Time PCR System (Applied Biosystems). Primer sequences are listed in Supplementary Table 6. The relative amount of mRNA normalized to cyclophilin B was calculated by using the delta-delta method[2].

Microarrays

Transcriptional profiling was performed using Affymetrix Mouse Genome 430 2.0 arrays. Microarray. CEL files were normalized and summarized using the MAS5 algorithm [3]. Mean values, variance, paired scatter plots as well as Principal Component Analysis plots were assessed and met quality criteria. Statistically significant changes in expression were assessed using the SAM algorithm[4]; to minimize false positives we considered only the probes with maximum expression values across samples greater than 100. Genes with SAM p<0.01 and fold change >1.5 were considered as up- or down-regulated. Using complete-linkage distance calculation and K-means clustering we found that these probes could be grouped into 9 robust clusters based on their expression levels upon exposure to Dex and TNF (Fig1C, Supplementary Table 1). Annotation enrichment analysis was performed using DAVID 6.7[5]. These data were also compared to arrays generated from adipose tissue of a diet-induced obesity model (GEO ID GSE2952)[6]. To compare the gene expression profile from different platforms, we first collapsed probe ID based gene expression matrices to a gene ID based matrix. If a gene was associated with multiple probes, only the probe with maximal signal across all samples was used to represent each gene.

Plasmids

VDR, Tmem176a, Colq, and Serpina3n was subcloned into pCDH at the EcoRI (BamHI for Serpina3n) and NotI sites and hairpins targeting Vdr, Nr3c1, and Rela were subcloned at AgeI/EcoRI or purchased from Open Biosystems. Hairpin sequences are shown in Supplementary Table 6.

ChIP-Seq

3T3-L1 cells were treated with 1% formaldehyde for 10 min at 37°C. ChIP and Illumina sequencing library construction were performed as described in Mikkelsen et al.[7] and sequenced on an Illumina HiSeq 2000. Length of reads is 36nt. Aliquots containing > 107 nuclear equivalents were subjected to immunoprecipitation as described previously[7, 8]. Briefly, pellets were thawed and lysed for 20 min on ice. The chromatin was then fragmented to a size range of ~200 to 600 bp using a Branson 250 digital sonifier (sonication conditions were separately optimized for each sample). Solubilized chromatin was then diluted and incubated with 1–2 μg antibody overnight at 4°C. Immune complexes were captured with ~ 0.02 ml protein A-sepharose, washed and eluted. Enriched chromatin was then subjected to crosslink reversal and proteinase K digestion at 65°C, phenol-chloroform-isoamyl alcohol extraction and ethanol precipitation. Isolated ChIP DNA was resuspended in RNAase treated water and quantified using the Qubit assay (Invitrogen). ChIP assays were performed using the following antibodies: H3K4me1 (Abcam ab8895, lot 38311/659352, H3K4me3 (Milipore 07-473, lot DAM1623866), H3K27ac (Active Motif 39133, lot 31610003), H3K36me3 (Abcam ab9050, lot 499302, H3K27me3 (Millipore 07-449, lot DAM1514011), H3K79me2 (Cell Signaling 9757, lot 1). All antibody lots were prevalidated for specificity and efficacy by Western blots and dot blots, as described [9]. The GR antibody was kindly provided by Dr. Wally Wang (UC Berkeley).

In vitro ChIP-PCR studies

Cells were cross-linked with 1% formaldehyde for 10 min at 37°C. Genomic DNA was sheared using a sonic dismembrator model 100 (Fisher Scientific) at half maximum speed for 10 sec. IP was performed according the manufacturer’s protocol (EZ-ChIP, Millipore) using 2.5ug of anti-GR (M20, lot F2812, Santa Cruz; gift from Dr. Wang), anti-p65 (C-20, lot J0311, Sigma), anti-VDR (C-20, lot G2511, Santa Cruz), anti-RXR (SC-774, lot G2310) or normal IgG (Sigma). Primers for ChIP-PCR studies are enlisted in Supplementary Table 4. All data are corrected for multiple hypothesis testing. All data is normalized to input.

In vivo ChIP-PCR studies

Epididymal adipose tissue was collected from chow (n=15) and high fat fed (n=7) C57BL/6 mice. Pooled fat pads were minced and dounce homogenized with 10 strokes in hypotonic lysis buffer (10mM HEPES, pH7.5, 10mM KCl, 1.5mM MgCl2, 250mM Sucrose, 0.5% NP40, and protease inhibitor cocktail). Lysates were filtered through a 100um cell strainer and spun at 1500g for 5 min. Lipid and cytoplasmic fractions were removed and the nuclear pellet was resuspended in lysis buffer, cross-linked with fresh formaldehyde (1%) for 5 min at room temperature, quenched with glycine (125mM), and washed twice with PBS. Nuclei were subjected to chromatin immunoprecipitation using the EZ-CHIP kit (Millipore).

Cell fractionation and GR localization

Trypsinized cell pellets were resuspended with Buffer A (10mM HELES, pH7.9, 1.5mM MgCl2, 10mM KCl, 0.5M DTT, and protease inhibitor cocktail), kept on ice for 5 min, and dounce homogenized with 20 strokes. Dounced cells were pelleted at 228g for 5 min. The supernatant was saved as cytoplasmic fraction. The initial nuclear pellet was resuspended in buffer S1 (0.25mM Sucrose, 10mM MgCl2) overlaid with an equal volume of buffer S2 (0.88mM Sucrose, 0.5mM MgCl2) and pelleted at 2800g for 5min. Lysates from each fraction were subjected to immunoblotting.

Glucose uptake assay

Assays were performed as described previously[10]. Briefly, mature adipocytes were serum-starved for 3–6 h in DMEM and then incubated with KRH buffer (121mM NaCl, 5mM KCl, 0.33mMCaCl2, 1.2mM MgSO4, 12mM Hepes, pH 7.4) with 25nM insulin at 37°C for 20 min and 2-deoxy-D-[2,6-3H]glucose (0.33 μCi/ml) for an additional 10 min. Uptake was stopped by three rapid washes on ice with KRH containing cytochalasin B (Sigma), the cells were solubilized with KRH buffer containing 0.1 % SDS, and radioactivity was measured by liquid scintillation counting.

Mass Spec analysis of histone modifications

Histones were acid extracted from cells as previously described[11] and derivatized with propionic anhydride before and after trypsin digestion. Briefly, 100 μg of histone in 40 μl of water were mixed with 20 μl of 3:1 anhydrous isopropanol to propionic anhydride (Sigma) and immediately neutralized with ammonium hydroxide. Samples were incubated for 15 minutes at 37°C and volumes were reduced in a vacuum concentrator. This constituted one round of “propionylation”. Samples were propionylated one additional time, and samples were digested with trypsin at the ratio of 1:50 trypsin to histone for 6–8 hrs at 37°C. Following quenching with glacial acetic acid, the samples were propionylated two additional times to add propionyl groups to the N-termini of peptides. Salts were purified from histone peptides using C18 STAGE-tips constructed as previously described [12]. Peptides were separated via reversed-phase HPLC using a 75μm inner diameter column packed with 10–15 cm of 5 μm C18 (Michrom, Auburn, CA). An LTQ-Orbitrap-XL was used to detect the peptides with full scans of m/z 290–1000 and a resolution of 30,000 in the Orbitrap. Data were analyzed using a previously described algorithm [13, 14]. Briefly, the algorithm considers the following information when identifying modified histone peptides; isotopic distribution, the MS/MS fragmentation spectra, and hydrophobic elution properties relative to other modified peptides. The modified peptides identified by the algorithm were quantified by integrating the area under the peak for the extracted ion chromatograms for all observed charge states. All modified forms for a given peptide were used to determine the “relative abundance” for any given form.

High throughput sequencing

Sequencing libraries were prepared from ~1–5 ng ChIP (or input) DNA as described previously[9]. Gel electrophoresis was used to retain library fragments between 275 and 700 bp. Prior to sequencing, libraries were quantified by fluorimetry and validated by quantitative PCR analysis of expected target regions. Sequencing was performed using the Illumina Genome Analyzer 2 or HiSeq2000 according to standard operating procedures. At least 20 million 36 base reads were generated for each data set.

Sequence analysis

All sequence analyses were based on mouse mm9 reference genome and the associated annotations (http://genome.ucsc.edu). Sequence reads were compiled, postprocessed and aligned to the reference genome using MAQ[15]. A maximum of 2 mismatches were allowed and positions were randomly selected for reads with multiple hits. Duplicate reads that aligned to the same position and strand were only counted once to avoid potential PCR bias.

ChIP-Seq fragment density histograms were generated by extending each aligned read to an assumed length of 200bp and then counting the number of extended fragments that overlap every 25th nucleotide in the genome. Regions enriched for each of the histone modifications were identified as previously described[16]. To compare histone modification peaks across different treatments, significantly enriched regions from individual treatments were merged to provide a set of comparable genomic regions for which ChIP-Seq fragments were recounted for each treatment using bedtools coverage v2.16.1 program[17]. Regions that showed differences in histone modification enrichment were identified using the edgeR algorithm[18]. Regions with p value < 0.05 and minimal threshold of 20 reads in at least one sample were selected as differentially-enriched regions (DERs). We assessed by qPCR the nearest gene to all concordantly up-regulated H3K27ac DER peaks that had poor probes or were not represented on the microarray at all and added validated genes to the set of upregulated genes from the array.

Chromatin State Modeling

We applied ChromHMM[19] to ChIP-seq data for H3K4me3, H3K4me1, H3K27ac, H3K36me3, H3K27me3, and WCE to model chromatin states for each of our current TNF, Dex, and control time points, along with our previously generated ChIP-seq data sets from L1 adipocyte differentiation[16]. We varied the bin size and state number and determined 14 states with bin size of 200bp to best model enhancer and promoter states. Based on emission parameters and enrichment profiles for each chromatin state at established functional loci, we assigned likely candidate state functions, with states 1–8 associated with promoters and enhancers (Supplementary Fig. 3a). For each condition, we concatenated adjoining states from 1–8 into preliminary promoter regions calls as long as they included at least 600bp of state 3 or 4. Otherwise, we preliminarily considered all such concatenated regions with >600bp of states 1–8 as enhancers. We then merged preliminary calls for individual conditions and assigned all consistently called regions as definite promoters or enhancers. Ambiguous calls (i.e. called as a promoter in one or more conditions and as an enhancer in others) were assigned as promoters if they overlapped a TSS and as enhancers if they did not. We thus assigned 13,298 promoters, of which 12,369 (93%) overlap RefSeq TSS and 69,801 enhancers. To determine whether a promoter or enhancer was gained or lost, we required a change from no identification of a promoter/enhancer (<600bp of states 1–8) to an active state (>600 bp of states 3–4 for promoters or states 5–6 for enhancers).

Sequence Conservation

Conservation scores were acquired from UCSC genome browser (http://hgdownload.cse.ucsc.edu/goldenPath/mm9/phastCons30way/). Conservation scores were averaged in 100bp bins over the windows compared in our analysis.

Motif analysis

For each H3K27ac DER set (Supplementary Table 3) less than 5,000 bp, a 500bp central sequence was evaluated for motif instances using www.biobase-international.com match program with the TRANSFAC vertebrate non-redundant-miniFP motif database as of 2012 and default parameters other than using only high quality matrices[20, 21]. Motifs were enumerated for (1) the regions corresponding to the specific H3K27ac DERs of interest, and (2) two separate background sets of 2000 randomly selected H3K27ac enriched regions. Enrichment scores for each motif and species were then defined as the number of instances detected in (1) divided by the number of instances detected in each set from (2). Fisher’s exact test applied to calculate significant enrichment or depletion. We found no differences between comparisons using the two background sets, so enrichment scores and p-values presented are based on one (provided in Supplementary Table 2).

GR ChIP-Seq analysis

GR peak calling was performed as described for modified histones. 200bp sequences centering at the peak summit were scanned using the FIMO algorithm in the MEME package[22] with default parameters to identify transcription factor motifs defined in the Transfac database (BIOBASE, Wolfenbüttel, Germany). For de novo motifs, 200bp sequences around the peak center were input to the GLAM2 algorithm of the MEME suite using default parameters other than maximizing insertion and deletion penalties. GR target genes were defined as nearest genes to GR peaks based on the distance between the peaks to transcription start sites (TSS) of annotated RefSeq genes. At each time point, we compared absolute gene expression fold changes (Dex or TNF treatment vs. control) of GR targets to those of all other genes using a one-tailed t test. For each GR peak summit, we counted H3K27ac reads within 1000bp in each 100bp bin. The number of read number in each bin was then normalized using corresponding library sizes (total read number). A set of 2000 H3K27ac peaks was randomly selected from the H3K27ac peaks that were not up- or down-regulated at any time points and not overlapped with any GR peaks. The H3K27ac reads within 1000bp windows of this set were counted and normalized as background. GR ChIP conservation comparison was performed as for H3K27ac peak conservation.

Accession numbers

Primary: ChIP-Seq data can accessed at GEO (GSE58491). Microarray data can be accessed at GEO (GSE62635).

Secondary: Microarray data from adipose tissue of a diet-induced obesity model[22] can be accessed at GEO (GSE2952). Microarray and histone modification Chip-seq data from 3T3-L1 model of adipocyte differentiation can be accessed at GEO (GSE20752).

Statistical analysis

No statistical method was used to predetermine sample sizes for these experiments. All observed data were included without exclusion criteria. HFD and rosiglitazone treatment groups were randomized by alternating order of mice encountered. Investigators were not blinded to allocations during animal studies. Glucose uptake, ChIP-PCR, and qPCR data were normally distributed and variances were estimated to be equal between groups.