Comparative Transcriptomic and Epigenomic Analyses Reveal New Regulators of Murine Brown Adipogenesis

Increasing energy expenditure through brown adipocyte recruitment is a promising approach to combat obesity. We report here the comprehensive profiling of the epigenome and transcriptome throughout the lineage commitment and differentiation of C3H10T1/2 mesenchymal stem cell line into brown adipocytes. Through direct comparison to datasets from differentiating white adipocytes, we systematically identify stage- and lineage-specific coding genes, lncRNAs and microRNAs. Utilizing chromatin state maps, we also define stage- and lineage-specific enhancers, including super-enhancers, and their associated transcription factor binding motifs and genes. Through these analyses, we found that in brown adipocytes, brown lineage-specific genes are pre-marked by both H3K4me1 and H3K27me3, and the removal of H3K27me3 at the late stage is necessary but not sufficient to promote brown gene expression, while the pre-deposition of H3K4me1 plays an essential role in poising the brown genes for expression in mature brown cells. Moreover, we identify SOX13 as part of a p38 MAPK dependent transcriptional response mediating early brown cell lineage commitment. We also identify and subsequently validate PIM1, SIX1 and RREB1 as novel regulators promoting brown adipogenesis. Finally, we show that SIX1 binds to adipogenic and brown marker genes and interacts with C/EBPα, C/EBPβ and EBF2, suggesting their functional cooperation during adipogenesis.


Author Summary
Obesity and its related metabolic diseases are growing problems worldwide. Brown adipose tissue (BAT) with its capability of burning off fat to generate heat is now at the center of research interest as target of therapeutic intervention for obesity treatment. In order to get a better understanding of the molecular mechanisms and transcriptional programs underlying brown adipocyte differentiation, we profiled the epigenomic and

Introduction
Obesity and its associated metabolic complications such as diabetes are increasingly responsible for significant economic and social burdens in many countries worldwide. Physiologically, obesity develops when energy intake exceeds energy expenditure, and the current treatments of obesity have been primarily focused on reducing energy intake. Unfortunately, these measures were largely inefficient in maintaining long-term weight loss [1]. The recent discovery of thermogenic adipocytes [2][3][4][5] capable of burning fat in adult humans has provided an exciting new therapeutic approach for the treatment or prevention of obesity by increasing energy expenditure [6].
Fat cells are derived from multipotent mesenchymal stem cells (MSCs), which can give rise to muscle, adipose, bone, or cartilage cells when given appropriate environmental cues. These cells can be broadly divided into fat storage cells, such as white adipocytes (WA); and fat burning cells, which include classical and inducible brown adipocytes (BA) (also known as beige or brite adipocytes) [7]. Brown cells contain high density of mitochondria and dissipate chemical energy as heat through the action of the mitochondrial protein UCP1 (uncoupling protein 1). It is evident that increased activity of thermogenic brown cells has beneficial effects on whole body metabolic homeostasis [3,5,8] and various environmental cues such as cold exposure and chemical activation of the β-adrenergic pathway can significantly up-regulate BAT activity [2,9]. Over the last couple of years, a number of protein factors as well as long non-coding RNAs (lncRNAs) and microRNAs have been identified as regulators in this process. For example, the members of the bone morphogenetic protein (BMP) family, the PPARγ co-factor PGC1α, the transcription factors (TFs) PRDM16, EBF2, KLF11, the protein deacetylase SIRT1, the secreted factors IRISIN and FGF21 as well as lncRNAs Blnc1, lncBATE1, and microRNAs miR193/365 have been shown to be essential for thermogenic fat cell recruitment [6,[10][11][12][13].
To promote thermogenic adipocyte recruitment, it is necessary to have a fundamental understanding of the gene regulation networks that control brown and white adipogenesis and identify the key differences between these two morphologically similar but functionally distinct cell types. Gene regulation networks are composed of cis-regulatory elements and transregulatory protein factors. In response to environmental stimuli, trans-factors bind to cis-elements such as enhancers or silencers to modulate gene expression. Given the large size and complexity of mammalian genomes, it has been difficult to systematically identify cis-regulatory elements at the genome-wide level. The recent discovery of signature histone modifications for these cis-elements (eg. H3K27 acetylation for active enhancers) and the advance in massive parallel DNA sequencing facilitated to comprehensively define these elements. In addition to typical enhancers, a group of so-called super-enhancers was discovered recently. These super-enhancers are large enhancer clusters containing high-density transcription factor binding and are associated with cell type specific genes [14]. They are also stronger in terms of gene activation ability and play key roles in controlling cell identity in mammals.
Trans-regulatory factors are mainly identified through differential gene expression and genetic analyses. For example, the prominent adipogenic factors PPARγ and C/EBPα are novel regulators, resulting in the identification of four activators of BA differentiation in this study.

Results
Comprehensive profiling of the transcriptome and epigenome during murine brown adipogenesis To examine the molecular control of cell fate transitioning from uncommitted progenitor cells to BAs, we used C3H10T1/2 MSCs and differentiated them into BAs, following a previously established protocol [24], where the multipotent progenitors were first committed to the brown lineage by BMP7 treatment before differentiation was triggered using a chemical cocktail ( Fig 1A). In addition, the well-established WA differentiation model 3T3-L1 was included in this study for comparison of events during differentiation. We first confirmed the efficiency and specificity of our differentiation systems by visual inspection of cell morphology, by qRT-PCR, and Western blot analyses of lineage marker gene expression (S1 Fig). The differentiation process for both lineages was highly efficient, as virtually all cells accumulated lipid droplets by day 7 of differentiation (S1A and S1B Fig). As expected, Pparg2, the master regulator of adipocyte differentiation, and other adipogenic marker genes such as Fabp4, CD36, Lpl, Adipoq and Cebpa were strongly up-regulated in both lineages in response to differentiation signals. In contrast, BA marker genes including Ucp1, Cidea, Elovl3, Ppara and Prdm16 were activated only in mature BAs, and the mitochondrial marker genes Cox7a1 and Cox8b were expressed much higher in mature BAs (S1C and S1D Fig). Corresponding expression patterns were also detected at the protein levels for PPARγ, UCP1, PPARα and CIDEA (S1E Fig). These data indicated that our differentiation processes were specific and efficient.
Next we profiled the transcriptome and epigenome during murine brown adipogenesis at five key time points: (1) day -3 (d-3, uncommitted progenitors); (2) day 0 (d0, end of brown lineage commitment by BMP7 treatment); (3) 6 hours (6h, end of epigenomic transition [23]); (4) day 2 (d2, early BA differentiation) and (5) day 7 (d7, mature BAs) (Fig 1A). For transcriptome profiling, we used RNA-seq for mRNAs and lncRNAs, and an array-based method for microRNAs. To validate our transcriptomic analysis during brown adipogenesis, a second replicate of the RNA-seq experiment was performed and the results indicated that the data were highly reproducible (S1 Table). In parallel, we also profiled the transcriptome using RNA-seq for mRNAs and lncRNAs, and microarray for microRNAs during 3T3-L1 WA differentiation. When compared with the transcriptomes of mouse adipose tissues, we found that our in vitro BA and WA systems are closely related to their corresponding in vivo tissues (S1 Table). To complement the analysis of transcriptional changes during brown adipogenesis, we also performed a comprehensive profiling of the dynamically changing chromatin landscape during BA differentiation by ChIP-seq. In this effort, we mapped a number of key chromatin marks including H3K4me1, H3K4me3, H3K9ac, H3K27ac and H3K27me3 during BA differentiation. We also performed replicates at two key time points (d0 and d7) for all histone marks and the results showed that our ChIP-seq data were highly reproducible (see S1 Table and S2A-S2C Fig). In addition, we profiled PPARγ binding using ChIP-seq in mature BAs where it is highly expressed. Examples of the epigenomic and transcriptomic landscapes as well as PPARγ binding during BA differentiation at the brown selective genes Cidea, Ucp1 and Ppara are shown in Fig 1B and S3 Fig. A corresponding epigenomic dataset for WA differentiation has been generated previously [15] and was used for subsequent comparative analyses. Prior to further analysis we validated our ChIP-seq datasets by examining the correlations between gene expression and various histone modifications. As expected, we found that highly transcribed genes were marked by active chromatin marks (H3K9ac, H3K4me1, H3K4me3, and H3K27ac) but not the repressive mark H3K27me3 at their promoters (S1 Table and S2D Fig). Together, these datasets constituted comprehensive reference maps of the epigenome and transcriptome for both BA and WA differentiation.

Stage-and lineage-specific gene expression during BA and WA differentiation
We first focused on genes that were dynamically regulated at different stages during adipogenesis, following the rationale that differentiation stage-specific expression mirrors functional roles for those genes. To this end, we systematically examined coding genes, lncRNA genes, and micro-RNA genes (Fig 2, S4 and S5 Figs, S2 Table). Using an entropy based method (See "Materials and Methods" section for details), we identified a total of 2277 (BA) and 1513 (WA) differentiation stage-specific coding genes (FPKM>5), which were 26.2% and 16.0% of the expressed genes in the respective lineages. The higher number and proportion of genes with dynamic expression in BA are in agreement with the requirement of executing additional gene programs to commit MSCs into the adipogenic lineage before differentiation. We noted a clear separation into five stages with little overlap of stage-specifically expressed genes in WAs indicating a strictly step-wise differentiation process. Specifically, 3T3-L1 cells are at the proliferation stage at d-4; at d0, these cells have been under growth arrest for 2 days [25]; after the adipogenic induction, the arrested cells re-enter the cell cycle and undergo an epigenomic transition stage at 6h [23]; while at d2, these cells are arrested again and start to differentiate [25]; finally at d7, 3T3-L1 cells are fully differentiated into mature adipocytes. In contrast, in BAs we observed a more substantial transition in gene expression between 6h and d2, whereas the time points before (d-3, d0, 6h), and after (d2, d7) showed a certain overlap of gene expression. This differentiation stage-specific gene expression pattern is likely derived from the fact that C3H10T1/2 cells continuously proliferate from d-3 to 6h without the contact inhibition and growth arrest stages observed in 3T3-L1 cells, and these cells start to accumulate fat earlier than 3T3-L1 cells (S1A and S1B Fig) at d2. To analyze our observations more systematically, we performed gene ontology (GO) analysis of stage-specific genes. The top category of enriched genes before differentiation (d-3/d-4) was "cell cycle" in both lineages, whereas the transient enrichment of "chondrocyte differentiation" was only found in BA after lineage commitment. Strikingly, the same enriched gene categories topped the list in BA at d2 and d7, i.e. "brown fat cell differentiation", "mitochondrion", and "lipid metabolic process"; but in WAs "fat cell differentiation" tops the list not before d7. This correlates with a well advanced differentiation status and accumulation of lipid droplets by d2 in BA, but not WA (see S1A and S1B Fig), and may explain similar gene expression pattern between d2 and d7 in BA.
We also examined lncRNA genes that are dynamically regulated during adipogenesis (S4B Fig and S2 Table). Using the NONCODE database [26], we found in total 1985 and 2796 expressed putative lncRNA genes during BA and WA differentiation (FPKM>0.5), respectively. Among them, 857 (BA) and 1135 (WA) lncRNAs showed a stage-specific expression pattern, which were 43.2% and 40.6% of the expressed lncRNAs, respectively. The proportions of stage-specific lncRNAs in BA and WA were therefore considerably higher than the ones for mRNAs which were 26.2% and 16.0%, suggesting lncRNA genes were regulated more dynamically during adipogenesis than coding genes. This trend was maintained even when the  comparison was limited to lowly expressed mRNAs with similar expression levels as lncRNAs, which turned out to be the least dynamic between stages (15.5% and 10.8%). In addition, we also found a previously identified lncRNA (Blnc1) that drives thermogenic adipocyte differentiation [11] to be specifically expressed in mature BAs (S4B Fig). Finally, we profiled micro-RNA gene expression along the same process, leading to the identification of known general adipogenic microRNAs (e.g. miR-378), brown lineage-specific microRNAs (miR-193), as well as several microRNAs not implicated in adipogenesis so far (S5 Fig and S2 Table).
To identify lineage-specific genes, we compared the 431 genes that were stage-specifically expressed in mature BAs (d7, Fig 2A) to those 420 genes that were specifically expressed in mature WAs (d7, Fig 2B). Among them, we found that 121 genes were robustly expressed specifically in mature BAs but not WAs, and 132 genes were only expressed in mature WAs ( Fig  2C). To further compile a list of putative lineage-specific markers for BA and WA, we selected the genes showing a similar lineage-specific expression pattern both in mouse BAT/WAT tissues and primary brown/white adipocytes, according to previously published data [27] (See S4A Fig for examples). The list for BA-specific genes contained a number of classic BA markers, such as Ucp1, Elovl3, Ppara and Cidea. In addition, Slc36a2 (also known as Pat2), recently described as a brown/beige-specific surface marker [28]; Cpt1b, the rate-controlling enzyme for long-chain fatty acid β-oxidation and several other mitochondrial protein genes (Chcd10, Sirt3, Mtfp1, Aspg and Adssl1) were also identified in this list. This observation is consistent with increased number and activity of mitochondria in BAs. Of note, we also found the gene encoding the kinase PIM1 that was specifically expressed in BAs, primary BAs, and BAT. In addition, we also noticed increased Pim1 expression upon cold exposure and chemical activation of the β-adrenergic pathway ( Fig 2D). PIM1 belongs to a group of constitutively active serine/threonine kinases and has been implicated in a number of biological functions such as apoptosis and cell cycle regulation. Recent reports suggested that PIM1 might play a role in cellular metabolism by modulating the phosphorylation status of AKT [29] and AMPK [30]. Based on this evidence, we selected PIM1 for further functional analysis as a potential regulator of BA differentiation (detailed below). For WAs, only few markers were established before, of which we rediscovered Nrip1 (also known as Rip140), a co-repressor that plays an important role in repressing a number of brown selective genes [31]. Another WA-specific gene, Trem2, was recently shown to enhance adipogenesis, promote glucose and insulin resistance, and diminish energy expenditure. Several other genes identified, such as the nuclear receptor Ear2 (Nr2f6), the ubiquitin gene Ubd (Fat10), the free fatty acid receptor Ffar2, and the insulin-like growth factor Igf1, were shown to be involved in metabolism without a clear role in white adipogenic differentiation. Finally, we also provide a list of lineage-specific lncRNAs (S2 Table,   for examples). Together, this transcriptomic dataset provides a valuable resource for the identification and further characterization of novel regulators for brown and white adipogenesis.
To ask how the lineage-specific and the commonly expressed genes in BA and WA are regulated at the chromatin level, we examined the histone modification dynamics at the gene promoters throughout both BA and WA differentiation. As expected, chromatin marks for active promoters such as H3K4me3 and H3K27ac correlate well with gene activity at the promoter regions in both lineages. Given that a recent study [32] suggested that the removal of H3K27me3 is required for brown gene expression, it was not surprising to observe a decrease of this mark at the promoters of brown specific genes in BA ( Fig 2E). Surprisingly, in WA, where these BA selective genes were not expressed, we also found a significant decrease in H3K27me3 at their promoters ( Fig 2E and S6A Fig), suggesting that the removal of H3K27me3 is not sufficient to induce the expression of these brown specific genes. Intriguingly, in contrast to WA, we found significantly higher levels of H3K4me1 at the promoters of BA specific genes throughout BA differentiation (Fig 2E and S6B Fig). This observation suggested that the pre-deposition of H3K4me1 at the early stages of brown adipogenesis was required for efficient expression of these genes at the late stage, while the removal of H3K27me3 was necessary but not sufficient to promote brown gene expression. In parallel, we found that general adipogenic genes are only marked by H3K4me1 but not H3K27me3 during both brown and white adipogenesis, suggesting their activation does not involve H3K27 demethylation.
TF binding motif analysis in stage-specific enhancers during BA and WA differentiation Enhancers are cis-regulatory elements that can activate gene expression over distance. It has been shown that enhancers are highly dynamic and play an important role in cell fate transitions [21]. To examine the dynamic regulation of enhancers during BA and WA differentiation, we identified stage-specific enhancers based on the enrichment of H3K27ac, a histone mark for active enhancers and promoters [33], and the lack of H3K4me3, a histone mark present at active promoters. We employed a similar entropy-based method as for the identification of stage-specific genes and found 24,002 and 13,429 genomic loci acting as putative stage-specific enhancers throughout BA and WA differentiation (Fig 3A and 3B). Again the higher number of dynamic enhancers in BA is in agreement with the additional commitment step in the differentiation of MSCs into the brown lineage, as compared to white adipogenesis starting from committed 3T3-L1 preadipocytes. Moreover, we observed the emergence of a distinct group of stage-specific enhancers after BMP7 treatment (compare d-3 to d0, Fig 3A), which suggested an epigenomic reprogramming during the process of brown lineage commitment. Consistent with a previous report [23], we also observed an epigenomic transition as evidenced by the formation of a new group of stage-specific enhancers within 6 hours after adipogenic induction (compare d0 to 6h, Fig 3A), while from d2 to d7, there are less changes as compared to the earlier stages ( Fig 3A). During white adipogenesis, the stage-specific enhancers at the early (d-2, d0) and late stages (d2, d7) show a certain overlap, while between d0 and d2, there is a relatively more drastic transition in enhancer formation. This pattern can be explained as at d-2 and d0, 3T3-L1 cells are under the growth arrest state [25]; after adipogenic induction, these cells go through clonal expansion between d0 and d2 [25]; beyond d2, the epigenomic reprogramming has been completed and the cells start to accumulate fat and subsequently enter the end differentiation stage at d7. To validate our analysis, we surveyed the levels of H3K4me1, another enhancer-associated epigenetic mark, at those loci and found a high concurrence between H3K27ac and H3K4me1 ( Fig 3A and 3B and S1 Table). Analysis of the genes present in the proximity of stage-specific enhancers showed that they fall into similar GO categories as the stage-specific genes analyzed earlier. This observation suggests that the stage-specific gene expression was likely regulated by the stage-specific enhancers, further confirming the role of enhancers in cell fate transitions [21].
Enhancers often serve as hubs for TFs. To identify potential TFs that bind to these stagespecific enhancers, we carried out motif analysis of enhancer associated sequences. De novo motif search led to the identification of several enriched binding motifs at each stage of BA and WA differentiation (S7 Fig). While some of the motifs are closely related to known TF binding motifs, (e.g. PPARγ in mature BA and WA) most of the motifs could not be assigned to known TFs due to our limited knowledge of the DNA binding motif for most TFs. Therefore, we examined these enhancers for the enrichment of known TF binding motifs derived from previous genome wide TF binding studies [34]. TFs with enriched motifs and robust expression at the corresponding stages are shown in Fig 3C and 3D. We found that motifs for well-known adipogenic regulators such as PPARγ, RXR, C/EBPα and FOXO1 were highly enriched in mature brown as well as white adipocytes. At the early stages, motifs for early adipogenic regulators including PBX1, KLF4 and STATs were enriched in either white or brown lineages (Fig 3C and 3D). Interestingly, the binding motif for SIX1, a homeobox transcription factor not previously implicated in the development of BAs was significantly enriched in mature BAs, but not WAs, suggesting a role for this factor in BA differentiation. To validate our finding in an in vivo setting, we also analyzed active enhancers in BAT and WAT tissue (using previously released datasets from ENCODE [35]) for enrichment of the SIX1 motif. Indeed, the SIX1 motif was found to be enriched at a much higher level in BAT than in WAT (S3 Table).

Super-enhancers mark key regulators of BA differentiation
It has been shown that key cell identity genes are often associated with super-enhancers (SEs) [36], a cluster of enhancers that are enriched for binding of TFs, mediator, and chromatin marks such as H3K27ac [14]. To search for novel regulators of BA differentiation, we sought to map the SEs and define SE-associated genes in both brown and white lineages and identify common as well as lineage specific SE genes. We employed the H3K27ac ChIP-seq data to define SEs because this allowed us to monitor SEs throughout the whole process of BA as well as WA differentiation and determine the SEs present specifically at the late (d2 and d7) but not the earlier stages (d-3 to 6h). To identify genes which are potentially regulated by these SEs, we filtered for those (1) within 100 kb of the SE and (2) whose gene expression patterns correlated with SE occurrence throughout differentiation. Through this approach, we identified 419 SEassociated genes for mature BAs and 417 SE genes for mature WAs (S4 Table), of which 324 were BA selective and 322 were WA selective (Fig 3E), respectively. As expected, well-known general adipogenic marker genes such as Cebpa, Fabp4 and Scd were associated with SEs in both lineages. And SE genes at late stages of BA differentiation included most key regulators of brown adipogenesis (Ucp1, Cidea, Fgf21 and Ppara) (Fig 3E and S8A Fig). Moreover, these

TF binding motif analysis of stage specific enhancers and identification of super-enhancer associated genes during brown adipogenesis. (A) and (B)
Differentiation stage-specific enhancers (defined as H3K27ac positive and H3K4me3 negative regions) are also marked by H3K4me1 and associate with genes reflecting their developmental stages. (C) and (D) TF binding motif analysis of stage specific enhancers with HOMER. "Site" indicates the bound TF(s), "Gene" indicates the corresponding gene coding for (one of) the TF(s). Expression of genes was determined using RNA-seq data. Sites with enriched motifs and correspondingly expressed genes were shown in the list. In differentiated BAs and WAs, the binding motifs for PPARγ, RXR and C/EBP were enriched, whereas the SIX1 motif was only enriched in BAs at late stage. (E) Genes associated with super-enhancers (SEs) in mature BAs and WAs. A subset of SE associated genes is shared between BAs and WAs, including established adipogenic markers such as Cebpa and Fabp4, whereas known brown and brite selective genes are specifically associated with SEs in BAs. (F) The Rreb1 gene was significantly up-regulated during brown adipogenesis and is associated with a SE in BAs.
doi:10.1371/journal.pgen.1006474.g003 genes tended to get transcriptionally activated (S8B Fig). Notably, Fabp3, Pdk4, Cpt1b and Cpt2, which were recently identified as putative SE associated genes in a human cell culture model of browning [10], were associated with SEs only in brown cells. Intriguingly, the TF RREB1 whose gene locus has been linked to metabolic traits like T2D susceptibility, fasting glucose levels, and body fat distribution [37][38][39] through genome-wide association studies (GWAS), was up-regulated during brown adipogenesis and also associated with a SE defined by H3K27ac enrichment (Fig 3F). This SE encompasses the entire promoter region of Rreb1. Using PPARγ binding as alternative method to define SEs in mature BAs, we found that Rreb1 was associated with one of the top SEs in BAs due to robust PPARγ binding upstream of its transcriptional start site. In addition to Rreb1, other SE genes determined by PPARγ binding signals include a whole panel of key brown cell markers such as Ucp1, Pgc1a, Cidea, Fgf21 and Ppara, and interestingly, Pim1 (S8C Fig). Based on the above observations, we selected Rreb1 for further functional analysis in BA development and function (detailed below).
BMP7 activates sox genes via a p38-dependent signaling pathway during BA lineage commitment BMP7 strongly promotes brown lineage commitment and differentiation in vitro and in vivo [24]. However, the detailed molecular mechanism underlying BMP7 function and its downstream targets during BA lineage commitment have not been thoroughly characterized. Therefore, we profiled the epigenomic and transcriptomic landscape in C3H10T1/2 cells treated with or without BMP7. To determine the molecular targets of BMP7, we compared the transcriptomic profiles between BMP7 treated and untreated C3H10T1/2 cells. On top of that, we also included the corresponding dataset from 3T3-L1 cells for comparison. To identify potential BMP7 targets, we focused on a specific group of 89 genes which showed a robust but transient induction after 3 days of BMP7 treatment (Fig 4A). Amongst those genes are modulators of WNT (Fzd9, Frzb) and TGFβ (Bambi, Scube3) signaling pathways, which were known to be involved in regulating adipogenesis (see S5 Table). In addition, the single most interesting group of genes consisted of two members of the SOX family of transcription factors, Sox8 and Sox13. Using slightly relaxed cutoff criteria, we noticed that five out of the total 20 Sox genes behave as putative BMP7 targets: Sox5, 6, 8, 9 and 13 ( Fig 4B). This observation suggested that at least part of the BMP7 response is mediated by the action of SOX proteins. Sox genes have been implicated in the regulation of embryonic development and in the determination of cell fate [40]. Sox9 expression in rat MSCs increases Cebpb expression and favors adipogenesis [41] and Sox5, 6, and 9 play important roles in chondrogenesis. Consistent with this, we observed a transient boost of chondrogenic gene expression after BMP7 treatment and enrichment of related GO categories (Fig 4C).
BMP7 activates two major signaling pathways, the SMAD, and the p38 MAPK pathways. Previous studies suggested that p38 signaling is more important for the formation of thermogenic cells and the activation of β-adrenergic pathway [24,42]. To determine if Sox genes are downstream targets of either the SMAD or the p38 pathway, we pre-incubated C3H10T1/2 cells with p38 inhibitors prior to treatment with BMP7, and surveyed their expression. As shown in Fig 4D, Sox gene activation was strictly dependent on p38 signaling as the treatment of p38 inhibitors PD169316 (PD) or SB202190 (SB) completely abolished the activation of all five Sox genes by BMP7. We also examined 27 additional genes from the list of BMP7 targets and found their expression was also dependent on p38 signaling (Fig 4E), which seemed to be the major transmitter of the BMP7 signal. Notably, at least two of those genes (Col11a2 and Col9a1) are well-known targets of Sox9. From the list of five Sox genes, we selected Sox13 for further functional studies (detailed below).
In another effort to identify relevant targets from our list of 89 candidates, we defined SEs according to the H3K27ac ChIP-seq data and generated a list containing genes associated with SEs at d0 upon BMP7 treatment. When we intersected both lists we found that 14 BMP7 induced genes were indeed associated with SEs (S9A Fig) and those genes might constitute another set of important targets of the cellular response to BMP7 activation. Amongst them, we found the fibroblast growth factor receptor Fgfr3, which is one of the receptors for FGF21 that promotes both BAT activation and subcutaneous WAT (scWAT) browning [43]. In addition to its robust induction by BMP7, we also observed significantly elevated levels of H3K4me1, H3K9ac and H3K27ac at the upstream enhancer region and increased H3K4me3 at the promoter of Fgfr3 gene (S9B Fig), underscoring its epigenetic regulation upon BMP7 treatment. Using different approaches of bioinformatics analysis, we identified a number of putative regulators for brown adipogenesis. From these candidates, we selected and validated the following four factors: (i) the kinase PIM1, found to be lineage-specifically expressed in mature BAs but not WAs (Fig 2C and 2D), and three TFs, (ii) SIX1, of which the binding motif was enriched in late stage BA enhancers (Fig 3C), (iii) RREB1, associated with a SE in BAs (Fig 3F and S8C  Fig), and finally (iv) SOX13, which is transiently induced by BMP7 during brown lineage commitment ( Fig 4B).
We first performed gain-of-function analysis of those factors during brown adipogenesis, by lenti-virally over-expressing the corresponding genes in the MSC line C3H10T1/2 before adipogenic induction without BMP7 treatment. As shown in Fig 5A, over-expression of each candidate, or EBF2, a known regulator of brown adipogenesis [16], significantly increased the differentiation efficiency of the C3H10T1/2 cells as monitored by Oil-Red-O (ORO) staining. Importantly, we also detected higher levels of brown / mitochondrial marker gene expression in the cells over-expressing the four candidate genes. These genes include the brown cell key regulators Prdm16 and Ppara, genes involved in brown cell function (Cidea and Elovl3), and genes essential for mitochondrial activity (Cox7a1 and Cox8b) ( Fig 5B). Moreover, Ucp1, the key thermogenic gene in BAs, was also up-regulated in cells expressing the four candidate genes with or without forskolin treatment (Fig 5C). And these mRNA expression changes were reflected at the protein levels as both CIDEA and PPARα proteins were up-regulated by the over-expression of the four candidates ( Fig 5D). In parallel, we also detected increased expression of general adipogenic genes such as Pparg2, Fabp4 and CD36 upon over-expression of the candidate genes, or Ebf2 (Fig 5E). This observation was consistent with increased lipid accumulation in the corresponding cells as determined by ORO staining (Fig 5A). Increased mitochondria activity leads to up-regulated oxygen consumption rate (OCR), and this is a key feature of thermogenic brown cells. To examine the effects of Pim1, Six1, Sox13, and Rreb1 over-expression on mitochondria activity, we measured the OCR in the corresponding cells 7 days after adipogenic induction. We found that both OCRs (Fig 5F) and other cellular metabolic parameters including basal respiration, proton leak, ATP production and maximal respiration ( Fig 5G) were significantly increased upon candidate over-expression. In addition, we also noticed a shift towards uncoupled respiration (Fig 5H), suggesting enhanced thermogenesis. The chemical activation of the β-adrenergic pathway by drugs such as norepinephrine can significantly stimulate BAT activity. To examine the effects of candidate over-expression on the response to β-adrenergic activation, we measured the OCR in corresponding cells after norepinephrine treatment. The results showed that cells over-expressing Pim1, Six1, Sox13, Rreb1, or Ebf2 (Fig 5J) were more susceptible to β-adrenergic activation than control cells ( Fig  5I), indicating enhanced thermogenic capability. Taken together, these results suggested that all four candidates either facilitate the commitment or the differentiation process from MSCs to functional BAs.
We showed that Pim1, Six1, Sox13 and Rreb1 were sufficient to promote brown adipogenesis. To test whether they were also necessary for BA differentiation, we performed loss-offunction analysis using Stromal Vascular Fraction (SVF) cells isolated from BAT transfected with LNA longRNA GapmeR oligonucleotides to knock down the genes of interest before adipogenic induction. As shown in Fig 6A, cells transfected with a scramble oligo (Scr) readily differentiated into mature BAs, whereas knock-down of either the candidate genes using two independent GapmeRs or Pparg led to severely reduced capabilities to differentiate as demonstrated by ORO staining. In addition, prominent brown / mitochondrial regulators and  (Fig 6B), as well as Ucp1 (Fig 6C) were down-regulated in cells with Pim1, Six1, Sox13, and Rreb1 knock-down. Consistent with the ORO staining results, adipogenic markers including Pparg2, Fabp4 and CD36 were also reduced (Fig 6D) by the knock-down of the candidate genes (Fig 6E). We also validated the function of the four candidates in SVF cells isolated from posterior scWAT. Those cells have a certain capacity to "brown" [44] and over-expression or knock-down of the four candidates resulted in similar outcomes as observed in the MSC and BAT SVF cell systems (S10 and S11 Figs).

SIX1 binds to brown marker genes and interacts with key regulators of brown adipogenesis
To gain further mechanistic insight into the mode of action for one of the identified factors, SIX1, we mapped its genomic localization via ChIP-seq in mature BAs. We found in total 7366 binding peaks for SIX1 with most of them located at intergenic regions, introns and promoters (Fig 7A), which is typical for TFs [45]. GO analysis of the SIX1 binding genes revealed that "regulation of generation of precursor metabolites and energy", "negative regulation of TGFβ receptor signaling pathway" and "brown fat cell differentiation" were amongst the most significantly enriched categories (Fig 7B). We detected SIX1 binding at the cis-regulatory regions (marked by H3K27ac) of brown markers such as Cidea and Ucp1 (Fig 7C). Moreover, we observed partial overlap of SIX1 binding to PPARγ binding at these regions. In a more quantitative analysis we measured the strength and proximity of SIX1 binding at brown-specific, white-specific and commonly expressed (i.e. white and brown) genes. We found that SIX1 bound preferentially around brown-specific and commonly expressed genes as compared to white-specific genes, suggesting a role for this factor in regulating brown selective as well as general adipogenic gene expression (Fig 7D). To decipher the molecular mechanism underlying SIX1 function, we performed a motif analysis of SIX1 bound regions. As expected, the most enriched binding motif was for SIX1 itself, which was followed by motifs for C/EBP, EBF, and NF1 TFs (Fig 7E). PPARγ and RXR binding motifs were only mildly enriched. The enrichment of C/EBP and EBF binding motifs at SIX1 binding sites suggested physical interactions between these TFs. Indeed, we verified the direct interactions between SIX1 and C/ EBPα, C/EBPβ, as well as EBF2 using co-IP assays (Fig 7F). Finally, using luciferase activity assay, we found that an upstream enhancer element of the Cidea gene harboring a SIX1 motif (S12 Fig) promotes expression in a SIX1-dependant manner (Fig 7G). Together, our findings corroborate a model in which SIX1 can be recruited to brown-specific or general adipogenic genes through either direct DNA binding (via the SIX1-binding motif) or recruitment by EBF2 and C/EBP proteins (at regions with no SIX1-binding motif).

Discussion
Promoting energy expenditure through thermogenesis is of significant interest as potential therapy for obesity and related diseases. It requires the recruitment of thermogenic fat cells such as brown and beige/brite adipocytes. Existing evidence suggests that the majority of these thermogenic cells are recruited de novo in response to environment cues [43,[46][47][48]. Therefore, to promote thermogenic adipocyte recruitment, it is essential to have a fundamental understanding of the gene regulation network that governs brown adipogenesis, especially at the lineage commitment step. In this study, we provide comprehensive profiles of the transcriptome and epigenome at five key developmental stages throughout the differentiation of murine multi-potent MSCs into mature BAs. Through in-depth bioinformatics analyses, we identified and functionally validated PIM1, SIX1, RREB1, and SOX13 as novel regulators promoting brown cell differentiation and function.
Differential gene expression analysis is a classic approach for the identification of regulators of cell type specification. A number of adipogenic and brown fat cell regulators including PPARγ, C/EBPα and PRDM16 were identified through this approach. In our study, we also used this analysis to identify brown selective genes but added additional criteria for the selection of candidates: these genes must be dynamically regulated during adipogenesis and stage-specifically expressed only in mature adipocytes (Fig 2). As the result, our list of 121 brown selective genes contains brown markers such as Cidea (#1), Elovl3 (#3), Ucp1 (#24) and Ppara (#61), as well as a number of mitochondrial genes including Cpt1b (#4). From this list, we specifically looked for factors that could potentially be involved in gene regulation or signal transduction, and we selected the kinase PIM1 (#53) for further analysis. Moreover, the Pim1 gene was later found to be associated with a SE in brown cells (S8C Fig). In our study, over-expression of Pim1 in both C3H10T1/2 cells and the scWAT SVF cells up-regulated a number of key brown cell marker genes as well as general adipogenic genes (Fig 5B-5E and S10B-S10D Fig). In addition, over-expression of this kinase also promoted the mitochondrial respiration in general and specifically uncoupled respiration, a feature of thermogenic fat cells (Fig 5F-5H and S10F and S10G Fig). In contrast, knock-down of Pim1 by GapmeRs reduced the expression of brown marker genes and adipogenesis efficiency in both primary brown cells (Fig 6) and subcutaneous white cells (S11 Fig). Therefore our analysis clearly implicates Pim1 in brown adipogenic differentiation, although future experiments will have to address if its role is solely restricted to the brown lineage. With our experimental model we cannot rule out that enhanced brown differentiation is partially caused by an increase in adipogenic differentiation in general. However, the lineage specific expression of Pim1 in BAT vs WAT, together with its increased expression in BAT upon cold exposure (Fig 2D), point towards a more specific role in BAT for this kinase. It will be interesting to investigate the functional significance and molecular mechanism of PIM1 in different tissues and under different metabolically challenging conditions such as diet induced obesity or cold exposure in vivo, especially the direct targets of this kinase.
Enhancer binding motif analysis is another powerful tool to identify novel TFs involved in specific cell differentiation processes [15,16]. In our study, we first defined stage-specific enhancers during both white and brown adipogenesis, then surveyed the enrichment of TF binding motifs at early and late stages of differentiation. While enrichment for several wellknown adipogenic factors at late stages of both brown and white adipogenesis was expected, our finding that the motif for the TF SIX1 was enriched during late brown adipogenesis was surprising (Fig 3C), since SIX1 has not been implicated in brown cell differentiation so far. Six1 belongs to the Sine Oculis Homeobox family of genes and has been reported to play a crucial role in muscle cell lineage decision as well as muscle development. Through gain-and loss-of-function assays, we confirmed that Six1 was required for the expression of brown selective and adipogenic marker genes in both the C3H10T1/2 cells and the SVF cells from scWAT and BAT. Moreover, over-expression of Six1 enhanced the mitochondrial uncoupled respiration in differentiated C3H10T1/2 cells (Fig 5H). In an attempt to decipher SIX1's mode of action, we performed genome-wide binding profiling of SIX1 in mature BAs. Strikingly, we found that SIX1 bound to brown as well as general adipogenic genes, some of whose expression were affected by the modulation of the Six1 gene, suggesting that these genes are direct targets of SIX1. Through binding motif analysis of SIX1 occupied regions and subsequent co-IP assay, we confirmed that SIX1 directly interacts and may cooperate with C/EBPα, C/EBPβ and EBF2 in regulating the transcription network during differentiation. Therefore the role of SIX1 seems not to be restricted to the regulation of BA differentiation, but it may act as a more general activator of the adipogenic transcriptional program. Again, in vivo studies to investigate the physiological role of Six1 in different adipogenic tissues will be instrumental to fully understand its biological function.
Super-enhancers are big clusters of enhancers and are often associated with genes specifying cell identity [14]. It can be defined by either mediator or TF binding or the enrichment of chromatin marks such as H3K27ac. Through analysis of super-enhancer associated genes, KLF11 was identified as an important factor promoting the browning of human mature white fat cells [10]. In our study, we profiled the stage-specific super-enhancers in a distinct process: i.e. the differentiation of brown adipocytes from multi-potent progenitor cells. That way we identified 419 genes associated with mature BA specific super-enhancers, including most brown marker genes. From this gene list, we selected the TF RREB1 for further functional analysis on the basis of its link to metabolic traits and association with one of the highest ranking SEs (#8) in BAs as defined by PPARγ binding. In our functional studies, over-expression of Rreb1 led to increased expression of brown marker genes and enhanced mitochondrial respiration in both C3H10T1/2 cells and SVF cells from scWAT. On the other hand, knock-down of Rreb1 resulted in reduced brown marker expression and impaired adipogenesis of SVF cells from BAT and scWAT. Moreover, RREB1 was very recently identified as a positive regulator of brown adipogenesis via a distinct bioinformatics approach [32] during the preparation of this manuscript. Further physiological studies on the role of Rreb1 in insulin sensitivity and energy homeostasis using corresponding gain-and loss-of-function animal models will be necessary to fully characterize the function of this positive regulator of brown adipogenesis.
BMP7 was recently established as a key factor that specifies the brown lineage from MSCs [24]. Mechanistically, BMP7 acts through either the SMAD or the p38 MAPK signaling pathway to induce its downstream target genes governing the adipogenic or thermogenic program. In order to gain an in-depth understanding of the brown lineage commitment from multipotent progenitor cells (as exemplified by BMP7 signaling), we systematically compared the gene expression profiles of C3H10T1/2 cells with or without BMP7 treatment. We found that 89 genes were transiently induced by BMP7 treatment, amongst which we identified a panel of Sox genes. We found that those genes were all p38/MAPK dependent, suggesting their involvement in promoting the thermogenic rather than the general adipogenic program. Subsequent functional analysis confirmed that one candidate gene, Sox13, promotes adipogenic differentiation, brown marker gene expression, and mitochondrial respiration. Our current results suggest that BMP7 triggers an early p38 dependent response, including the activation of Sox gene expression, important for the lineage commitment of brown cells from multi-potent progenitor cells. Clearly, future work is needed to dissect the exact contribution of the different Sox genes during this process and to define their downstream target genes.
Transcriptomic profiles and epigenomic landscapes are important resources for understanding the gene regulation network in a certain cell type or in a specific cell differentiation process. Analysis of those datasets has led to the identification of numerous novel regulators in various cellular processes. Seminal works in the field have provided valuable resources for further investigation of the molecular control of fat cell differentiation [10,15,23,[49][50][51][52][53]. Here we add a comprehensive study of the epigenomic and transcriptomic transitions at five key developmental stages throughout the process of murine brown adipogenesis. Our dataset comprises a high temporal resolution of the differentiation process as well as a broad coverage of chromatin marks. Through comparative analyses of white and brown datasets using various bioinformatics tools, we identified many potential candidates and validated four factors that promote thermogenic adipocyte differentiation in various cellular models including C3H10T1/2 cells, SVF cells from subcutaneous WAT and interscapular BAT. Moreover, through analyzing the chromatin dynamics at the promoters of lineage-specific and commonly expressed adipogenic genes in BA and WA, we found that, in addition to the mechanism proposed by a recent study [32], which suggested that the removal of H3K27me3 is required for brown gene expression, the pre-deposition of H3K4me1 at these genes during early stages of brown adipogenesis is essential for poising them for expression at a later stage. For general adipogenic genes, we found they are only marked by H3K4me1 but not H3K27me3 during both brown and white adipogenesis, suggesting their activation does not involve H3K27 demethylation. Based on these observations, we propose that the pre-deposition of H3K4me1 at brown specific genes is a critical step in the chromatin remodeling during the process of brown adipocyte lineage commitment, while the full activation of these genes is only possible once the "stop sign" (H3K27me3) is removed. Besides proposing this conceptual model for the epigenetic regulation of brown lineage specific genes, we anticipate additional factors (including lncRNAs and miRNAs), promoting or inhibiting brown cell differentiation to be identified through analyzing these datasets and further insights to be gained from these resources.

Cellular flux assay
To determine cellular oxygen consumption rates, we used the Seahorse XFe24 Extracellular Flux Analyzer. Cells were seeded on gelatinized XFe24 cell culture microplates (100777-004) at 4,000 cells/well and differentiated one day post confluency following the procedures described above. The XF Cell Mito Stress Test Kit (103015-100) was used to determine basal respiration, ATP production, proton leak, maximal respiration, and spare respiratory capacity. Concentrations of the added chemicals were: 40 μM oligomycin, 0.15 μM FCCP, 1 μM rotenone / 1 μM antimycin A. To activate the β-adrenergic pathway before measurement, cells were treated with 10 μM norepinephrine (Sigma Aldrich, A0937).

Gene expression analysis
Cells were harvested with TRIzol reagent (Ambion, 15596018) and total RNA was extracted following the manufacturer's protocol. RNA was treated with Amplification Grade DNase I (Invitrogen, 18068-015) before the generation of either RNA-seq libraries or cDNA to assess expression of individual genes. RNA-seq libraries were generated by BGI, China, following standard procedures. For cDNA generation, 500 ng of total RNA were reverse transcribed using random 9-mers and M-MLV Reverse Transcriptase (Invitrogen, 28025-021). Target gene expression was determined via quantitative Real-Time PCR analysis using Power SYBR Green PCR Master Mix from Applied Biosystems (4367659) on a 7900HT Fast Real-Time PCR machine (Applied Biosystems). All gene expression data in this study was normalized to the expression of the riboprotein gene 36B4 unless indicated differently. Primers were listed in S6 Table. ChIP-seq Chromatin immunoprecipitation was performed as described earlier [54]; antibodies used in this study were listed in S6 Table. To construct ChIP-seq libraries, we employed a method described previously [55]. In short, 5ng of ChIP DNA were used as starting material. After an end repair step, terminal addition of poly-dCs, and ligation of linkers, the DNA was amplified in a two-step PCR procedure. The final product was loaded onto a 2% agarose gel and fragments between 200bp and 500 bp were cut from the gel, purified and sequenced at BGI China.

Luciferase assay
A 222 bp fragment corresponding to the region -13,982 to -13,761 from the Cidea transcription start site (TSS) was PCR-amplified from mouse genomic DNA and cloned into pGL3-Basic plasmid (Promega) to yield pGL3-Cidea-Enh1. For the luciferase assay, HEK293 cells were transfected with pGL3-Cidea-Enh1 (or equal amounts of pGL3-Basic) and pCMV6-hSix1: myc:DDK from OriGene (RC203465) (or equal amounts of pCMV6:myc:DDK) using Lipofectamine 2000 (Invitrogen). 48h after transfection cells were lysed and Firefly luciferase activity was measured on a GloMax Multi luciferase reader using the Dual-Glo luciferase assay kit (Promega) following the manufacturer's recommendations.

ChIP-seq data analysis
All ChIP-seq datasets were aligned using Bowtie (version 2.0.4) to build version mm9 of the murine genome. Alignments were performed with the following additional parameters: -t -q -p 8 -N 1 -L 25. To visualize the ChIP-seq signals for each histone modification and PPARγ, we extended each read to 300 bp and counted the coverage of each read for each base, which was shown as the UCSC genome browser tracks. For the downstream analysis, we normalized the read counts for the ChIP samples by computing the numbers of Reads Per Kilobase of bin per Million reads sequenced (RPKM). To minimize the batch and cell type variation, the RPKM values were further normalized through Z-score transformation. MACS [56] was used to identify histone modification regions and PPARγ binding peaks by default settings.

RNA-seq data analysis
All RNA-seq datasets were aligned using Tophat (version 2.0.6) to build version mm9 of the murine genome. Alignments were performed with the following additional parameters:-p 2-solexa1.3-quals. The mapped reads were further analyzed by Cufflinks (v 2.1.1) [57] and the expression levels for each transcript were quantified as Fragments Per Kilobase of transcript per Million mapped reads (FPKM) based on refFlat database.

Identification of differentiation stage-specific coding genes, lncRNAs and microRNAs
To identify differentiation stage-specific coding genes, we used a strategy described previously based on the Shannon entropy to compute a stage-specificity index to all expressed coding genes [58][59][60] using the averages of two RNA-seq replicates. The entropy score for each gene was defined as described as follows [61]: for each gene, we defined its relative expression in a cell type i as Ri = Ei /∑E, where Ei is the FPKM value for gene expression in the cell type i; ∑E is the sum of FPKM values in all cell types; N is the total number of cell types. Then the entropy score for this gene across cell types can be defined as H = -1 Ã ∑Ri Ã logRi (1 i N), where the value of H ranges between 0 to log2(N). An entropy score close to zero indicates the expression of this gene is highly stage-specific, while an entropy score close to log2(N) indicates that this gene is expressed ubiquitously. Based on an examination of the entropy distribution genes with entropy less than 2 were selected as the candidate stage specific genes. Among these candidates, we selected genes for each stage based on the following criteria: the gene was highly expressed in this stage (FPKM>5), and high expression (FPKM>5) could not be observed in more than three additional stages. These genes were then reported in the final stage specific coding gene list and used for subsequent analyses, e.g. GO analysis (see below).
To identify stage specific lncRNAs, we used NONCODE v4 mouse annotation database [26] and removed the regions overlapping with refFlat annotation. The expressed putative lncRNAs are those with FPKM>0.5 in at least one of the five differentiation time points for each lineage. We used entropy score less than 2 to select candidates of stage specific lncRNA and applied the following criteria: the lncRNA was highly expressed in this stage (FPKM>0.5), and expression (FPKM>0.5) could not be observed in more than three additional stages for this lncRNA. These lncRNAs were then reported in the final stage specific lncRNA list.
The microRNA array signals were quantile normalized among different stages during adipogenesis and we selected those with expression varying by more than two-fold between different stages during adipogenesis. They were assigned to their specific stages based on their maximum expression.
Identification of day 7 adipocyte specific coding genes and lncRNAs Day 7 adipocyte specific coding genes were identified from the day 7 stage specific gene list obtained above. We selected the brown d7 adipocyte specific genes based on the following criteria: (1) brown adipocyte d7 expression was at least 3 fold higher than that in white d7 adipocytes; (2) this genes was not highly expressed in white d7 adipocyte (FPKM<10). To generate the list of lineage-specific genes in Fig 2C, we introduced two additional criteria: (1) expression in differentiated SVFs from BAT was at least 2.5-fold higher than in differentiated SVFs from WAT (data from [27]), and (2) expression in BAT tissue was at least 2.5-fold higher than in WAT tissue (harvested from male C57BL6 mice). The criteria for white specific genes were vice versa. Day 7 adipocyte specific lncRNAs were identified from the day 7 stage specific lncRNA list obtained above. To be defined as brown specific lncRNAs, their expression (FPKM) at day 7 had to be at least 5-fold higher in brown than that in white adipocytes. The criteria for white d7 adipocyte specific lncRNAs were vice versa.

Identification of stage specific enhancers
To identify enhancers (without overlapping promoters) we first called H3K27ac and H3K4me3 peaks for each stage by MACS. The peaks from all stages (BA: d-3, d0, 6h, d2, d7; WA: d-2, d0, d2, d7) were pooled for individual modifications and lineages. Enhancers were identified from the H3K27ac regions not overlapping with H3K4me3 regions and regions 2.5 kb up-and down-stream of Refseq transcription start sites. We calculated the RPKM within these enhancer regions and their normalized signals by Z-score normalization. Stage specific enhancers were defined based on the following criteria: the enhancer was highly active in this stage (normalized RPKM>1), and its high activity (normalized RPKM>0) could not be observed in more than three additional stages. These enhancers were then reported in the final stage specific enhancer list.

Identification of lineage specific coding genes and determination of the chromatin state at their promoters
We first extracted BA d7 and WA d7 stage specific genes as described above. Those genes were first pooled, and then separated into BA-specific, WA-specific and shared common adipogenic genes using the following criteria: a gene was defined as BA-specific if (1) brown adipocyte d7 expression was at least 3 fold higher than that in white d7 adipocytes; (2) this gene was not highly expressed in white d7 adipocyte (FPKM<10). The criteria for WA-specific genes were vice versa. The rest of the genes were identified as shared genes. Then we calculated the z-score normalized RPKM of each histone modification on promoter regions. Then all these normalized RPKM signals were mapped to its corresponding genes. The box plots show the normalized RPKM of different histone modification on different groups of genes. p-values were calculated using t-test.

Identification of super-enhancers and their potential targets
Super-enhancers of brown adipocyte and white adipocyte were identified by the software ROSE [36] based on our brown adipocyte d7 PPARγ ChIP-seq data and published 3T3-L1 d7 PPARγ ChIP-seq data. The closest gene to the super-enhancer was assigned as its potential target. Alternatively, SEs were identified using the H3K27ac datasets. To do this, we calculated Zscore normalized RPKM values on enhancer regions identified above and selected active enhancers with normalized RPKM >0 within 12.5 kb to further merge nearby signals into 'stitched' enhancers. We kept the 'stitched' enhancers larger than 12.5 kb as super-enhancer candidates. The RPKM values were calculated on these new super-enhancer candidates and ranked to find the super enhancers with slope > 0.5. To identify SEs specific to late BAs or late WAs, we used the corresponding d7 SEs and removed SEs which showed an overlap with SEs also present at early stages (i.e. d-3, d0, 6h for BA and d-2, d0 for WA). We assigned super enhancers to expressed genes (FPKM>1) within 100 kb of super-enhancer based on the Pearson correlation between expression and enhancer activity during adipogenesis. Genes with correlation coefficient higher than 0.75 were selected as the potential super-enhancer targets. If there was no gene with more than 0.75 correlation coefficiency within 100 kb of a superenhancer, we assigned the gene with the largest correlation coefficiency above 0.5 as its potential target. microRNAs within 100 kb of super-enhancer were all assigned as potential targets.

Motif analysis
To find the sequence motif enriched in the identified enhancers, we used the Homer [34] program based on mm9 genome using default settings. For motif analysis, background sequences were randomly selected from the genome, matched for GC% content to facilitate subsequent GC normalization.

Gene ontology analysis
To find enriched GO categories, we used the DAVID web-tool (using default settings) [62] with Gene Ontology database including Molecular Function terms, Biological Function terms and Cell Component terms. To find the GO categories enriched around enhancers, superenhancers and TF binding peaks, we used the GREAT tool [63] using default settings. Background calculation was based on the whole genome.
Distribution and enrichment of six1 binding peaks Six1 peaks were called by MACS as described above. The overlap between the peaks and annotated genome elements was calculated based on the mm9 genome annotation. To calculate Six1 enrichment, WA-specific, BA-specific and common genes were defined as follows: WAspecific genes: FPKM>5 in WA d7 and FPKM<5 in BA d7; BA-specific genes: FPKM>5 in BA d7 and FPKM<5 in WA d7; common genes show FPKM>5 in both WA d7 and BA d7. When comparing the Six1 signal around those groups of genes, we first defined an unbiased set of regulatory sites and potential Six1 binding site around the TSS using H3K27ac peaks from brown adipocytes (d7), H3K27ac peaks from white adipocytes (d7) or shared H3K27ac peaks from both brown and white adipocytes. Overlapping Six1 peaks were counted and enrichment was calculated as RPKM by adding the signals of Six1 peaks. Distances between genes and their closest Six1 binding sites were based on their annotated TSSs.

Accession numbers
All ChIP-seq, RNA-seq and microRNA microarray datasets were deposited in GEO under the accession number of GSE75698.

Ethics statement
All animal procedures were performed according to a protocol (IACUC#130829) approved by the Institutional Animal Care and Use Committee of the Agency for Science, Technology and Research (A Ã STAR) of Singapore. Enlarged view of the upstream region of the Cidea gene reveals the existence of an enhancer element bound by SIX1, PPARγ and marked by H3K27ac. A 222bp SIX1-bound region was PCR-amplified, cloned and tested for its regulatory potential using luciferase assay. SIX1 binding motif (TCAGGTTTC) was highlighted in bold letters. (TIF) S1