Remodeling of gene regulatory networks underlying thermogenic stimuli-induced adipose beiging

snATAC-seq reveals distinct cell types in inguinal adipose tissue

In order to understand the dynamic chromatin remodeling during adipose tissue beiging, we performed snATAC-seq of inguinal WAT (iWAT) from mice exposed to cold temperature (C; 6 °C) or treated with CL at room temperature (RT; 22 °C) for 0, 1, 3, or 7 days (Fig. 1a). Histological assessment showed both cold exposure and CL induced beige adipocyte formation in iWAT depots (Supplementary Fig. 1a). Nuclei from whole adipose tissues were used to obtain cell types in entire iWAT to overcome technical bias related to isolation of cells or SVF. A total of 32,552 cells passed quality control from seven groups (day0RT, day1C, day3C, day7C, day1CL, day3CL, and day7CL) and three individual mice for each group (Supplementary Fig. 2a–s). To identify the cell types, we performed dimensionality reduction and unsupervised clustering based on top 50% variable peaks using Signac21,22. This revealed 12 highly consistent cell clusters, visualized using Uniform Manifold Approximation and Projection (UMAP) (Fig. 1b). We annotated clusters using differentially accessible genes, gene ontology (GO) analysis, cell type specific marker genes, transcription factor (TF) motif enrichment, and GREAT analysis (Fig. 1c and Supplementary Fig. 3a–e). We further confirmed that the annotation is consistent with the cell type identities predicted by publicly available adipose tissue scRNA-seq and snRNA-seq data with high correlation rates (Supplementary Fig. 3f, g). Adipocytes constitute about 90% of adipose tissue volume but comprise only 20% of total cells23. In agreement with our expectation, we found an average of 22% of adipocytes in our snATAC-seq dataset (Supplementary Fig. 4a, b). The average proportion of other non-adipocyte cell types in our dataset was about 72% for immune cells, 3% for APCs, and 2% for stromal cells (Supplementary Fig. 4a).

Fig. 1: Cell type identification in inguinal white adipose tissue.

a Schematic outline of the workflow. Two-month-old male C57BL/6 mice were exposed to cold (6 °C) or CL-316,243 (CL; 1 mg/kg/mouse/day) for 1, 3, or 7 days (n = 3 per group). Cartoons of mice and thermometers created with BioRender.com. Inguinal adipose tissues collected at day 0, 1, 3, and 7 after exposure was subjected to snATAC-seq.  b UMAP plot of 32,552 cells that passed quality control. Colors represent 12 different clusters determined by unbiased clustering. c Heatmap of normalized accessibility of cell-type marker genes. d Relative proportions of each cell type per group (n = 3 per group). *Adjusted p-value<0.05, ANOVA multiple comparisons test with Bonferroni’s post-hoc test was performed. e, f Milo analysis of cell neighborhood abundance changes in Day 0 RT vs. Day 7 cold e, or Day 7 CL f on UMAP plots. Size of points indicates the number of cells in a neighborhood; lines represent the number of cells shared between adjacent neighborhoods. Points are neighborhoods (Nhood), colored by the log fold differences at the two time points (FDR 10%). Beeswarm plots show the distribution of the log-fold in defined clusters. Macro Macrophage, Dend Dendritic cell, APC1 Adipocyte progenitor cell 1, APC2 Adipocyte progenitor cell 2, VEC Vascular endothelial cell, VSMC Vascular smooth muscle cell.

Dynamic changes in adipocytes and identification of thermogenic beige adipocytes

We analyzed the changes in cell abundance in response to cold and CL treatment in each cell type (Fig. 1d). However, the differential abundance testing of discrete clusters relies on predefined clusters and limits the ability to detect shifts in cell states within clusters. Thus, to further analyze the abundance changes caused by cold and CL treatment without relying on discrete clusters, we used Milo24 to assign cells to neighborhoods having similar cellular states. The neighborhoods around adipocytes showed the highest variation in the log2 fold changes of their UMAP coordinates after cold and CL treatment. This indicated, that adipocytes underwent especially dynamic changes in their responses to cold and CL relative to other cell types (Fig. 1e, f and Supplementary Fig. 4c–h).

Mature adipocyte populations have high accessibility at the Pparg gene and differentially accessible peaks in these cells are enriched for PPARG motifs (Fig. 2a, b). Two distinct adipocyte populations were distinctively separated from other populations on the UMAP and identified as white or adipocytes and beige (Figs. 1b and 2c). The proportion of beige adipocytes gradually increased from 1% to 5% with increased exposure time to both cold or CL treatments (Figs. 1d and 2d). Among neighborhoods identified by Milo24, those around beige adipocytes were increased after cold and CL treatment compared to day 0 at RT, confirming accumulation of beige adipocytes (Fig. 1e, f and Supplementary Fig. 4c–h).

Fig. 2: Identification of beige adipocytes within mature adipocytes.
figure 2

a UMAP plot of normalized Pparg gene accessibility. b UMAP plot of normalized PPARG motif enrichment. c UMAP plot of white and beige adipocytes among mature adipocytes (n = 7,004). d Heatmap of the relative cell abundance in white and beige adipocyte clusters for each group. e Normalized accessibilities of adipogenic or common adipocyte genes and thermogenic or beige-selective genes. f Differentially accessible genes between white and beige adipocytes. Positive and negative log2-fold changes respectively indicate increased accessibility in white and beige adipocytes. Genes with adjusted p-value < 10−100 & Abs(log2 fold-change) > 0.25 are colored and labeled. g, h GO g and KEGG pathway h analysis of genes enriched in beige adipocytes compared to white adipocytes. X-axis indicates number of genes. Color indicates adjusted p-value.

General adipocyte marker genes were highly accessible in both adipocyte populations, while beige/brown adipocyte marker genes were highly accessible in beige adipocyte population (Fig. 2e). The top five genes enriched in white adipocytes were Car3, Lrp3, Tenm4, Trabd2d, and Cmklr1 (Fig. 2f); whereas the top five genes enriched in beige adipocytes were Slc4a4, Ucp1, Pank1, Ppargc1b, and Gyk (Fig. 2f). GO analysis revealed that the genes enriched in beige adipocytes are related to FA oxidation and thermogenesis (Fig. 2g). KEGG pathway analysis showed that the genes enriched in beige adipocytes are also associated with PPAR signaling pathway and FA metabolism (Fig. 2h). Taken together, the cluster we identified as beige adipocytes showed clear thermogenic signatures compared to white adipocytes.

Beige-specific cis– and trans-acting gene regulatory networks

Chromatin accessibility analysis reveals key GRNs including TFs and cis-regulatory elements underlying dynamic genetic and epigenetic programs. To comprehensively characterize GRNs activated in beige adipocytes, we identified genomic regions with higher accessibility in beige adipocytes, which we refer to as beige-specific peaks. We firstly performed GREAT analysis25 to predict biological functions of the beige-specific peaks using annotations of the nearby genes. Predicted functions of the beige-specific regions were related to brown fat cell differentiation and FA oxidation (Fig. 3a), providing additional confidence in our original identification of beige adipocytes. We further extended these analyses by focusing on distal beige-specific peaks not residing in annotated genes. Using cicero26, we analyzed co-accessible peaks and identify putative cis-regulatory interactions and the genomes that are more likely to be enhancers of the linked genes. For instance, a set of beige-specific peaks in the upstream of Ucp1 showed high co-accessibility with each other and Ucp1 promoter (Fig. 3b). Similarly, other sets of co-accessible beige-specific peaks were linked to genes including Ppara, Pdk4, Ppargc1b, Acot11, Arhgef37, Slc4a4, Col27a1, Dio2, Elovl6, Kcnk3, and Ppargc1a. These genes themselves exhibited high accessibility in beige adipocytes (Supplementary Fig. 5a). If intergenic beige-specific peaks include functional enhancers that are active in beige adipocytes, we expected to find H3K27ac marks associated with the beige-specific peaks. To make this determination, we referred published H3K27ac ChIP-seq data collected from white and beige adipocytes27. Over 6,000 beige-specific peaks overlapped with regions carrying H3K27ac (Supplementary Fig. 5b). These overlapping regions included peaks upstream of Ucp1 that were co-accessible with the Ucp1 promoter in beige adipocytes that were largely absent from white adipocytes (Fig. 3b). This indicates that the beige-specific peaks are likely to include a large number of functional enhancers active in beige cells.

Fig. 3: Beige-specific transcription factor activity and chromatin interaction networks.
figure 3

a Biological function of beige-specific peaks by GREAT analysis. b Genome tracks showing the aggregated snATAC-seq profiles near the Ucp1 locus in adipocytes. H3K27ac ChIP-seq signal tracks (n = 2 for each)27 show the aggregated H3K27ac signals in white and beige adipocytes. snATAC-seq signal tracks show the aggregated snATAC-seq signals in white and beige adipocytes. Red boxes indicate beige-specific regions overlapping with PGC-1α binding sites39. Fragment coverage show the single cell profile (each row) of randomly selected 100 single cells from white and beige adipocytes respectively. Links show cis-coaccessibility network with multiple connections between peaks around Ucp1. c Motif enrichment for white and beige adipocytes. Positive and negative log2-fold changes respectively indicate enrichment in white and beige adipocytes. Motifs with adjusted p-value < 10−30 & Abs(log2 fold-change) > 2.0 are colored and labeled. d Most variable TF motifs are ranked. Motif sequence logos for the top variable motifs are shown. e, f TF footprints for PPARG e and ESRRA f from (white) adipocytes and beige adipocytes. Tn5 DNA sequence bias was normalized by subtracting the Tn5 bias from the footprinting signal. X-axis is distance to motif center (bp). g Venn diagram showing the number of overlaps between white-specific or beige-specific peaks and PGC-1α binding sites39. hj The number of overlapping beige-specific peaks with PGC-1α binding sites39 h, ESRRs i, and MEFs j motifs from cisbp database41. The genes closest to the overlapping peaks are sorted by the number of peaks. k Genome tracks showing the aggregate snATAC-seq profiles near the Pdk4 locus in adipocytes. H3K27ac ChIP-seq signal tracks (n = 2 merged)27 show the aggregated H3K27ac signals from white and beige adipocytes. snATAC-seq signal tracks show the aggregated snATAC-seq signals in white and beige adipocytes. Red boxes indicate beige-specific regions overlapping with PGC-1α39, ESRRα42, and MEF2α43 ChIP-seq binding sites. Links show cis-coaccessibility network with multiple connections between peaks around Pdk4.

To gain further understanding of the GRNs in beige adipocytes, we contrasted TF binding motifs enriched in beige-specific, vs white-specific peaks. Motifs enriched in beige adipocytes included members of the ESRR, MEF2, ROR, NR2, and NR4 families (Fig. 3c and Supplementary Fig. 5c, d). It was previously reported that peroxisome proliferator-activated receptor coactivator-1α (PGC-1α), a key regulator for brown adipose, promoter has two MEF2-binding sites and its transcriptional activation is mediated by MEF228. Especially, ESRRα and ESRRβ were the most highly variable motifs in adipocytes (Fig. 3d). Whereas the footprint pattern of the general adipose regulator PPARg was similar between white and beige adipocytes (Fig. 3e), beige adipocytes showed a distinct TF footprint for ESRRα (Fig. 3f). Loss of ESRR isoforms was shown to reduce thermogenic capacity of BAT upon adrenergic stimulation29. These results collectively highlight specific regulation of ESRR and its importance in beige adipocytes. In addition to ESRR and MEF2, NR2 and NR4 family members can be functionally important to induce beiging. Especially, NR4A family members were shown to bind to a regulatory region of Ucp1 and enhance Ucp1 expression in mouse and human adipocytes30. On the other hand, the most enriched TF motifs in white adipocytes is NR3C family that serves as glucocorticoid receptors (GR) (Fig. 3c). GR can function as a TF and bind to glucocorticoid response elements. It was reported that excess glucocorticoids reduce the thermogenic activity of BAT and inhibit the expression of Ucp1, suggesting that glucocorticoids are negative regulators of BAT thermogenesis31,32,33,34. In accordance with this, GR ChIP-seq data showed that NR3C1 and NR3C2 were enriched during re-warming of beige adipocytes at thermoneutrality (30 °C), suggesting that NR3Cs are involved in whitening of beige adipocytes27. Collectively, these results demonstrate that our snATAC-seq data revealed positive and negative transcriptional regulators of beige adipocyte development.

PGC-1α is a key transcriptional coactivator of white-to-brown transition and regulates pathways including mitochondrial FA oxidation35 and mitochondrial biogenesis28,36,37,38. To get better insight into how PGC-1α functionally relates to beige-specific peaks we identified, we intersected the beige-specific peaks with PGC-1α binding sites which were obtained from publicly available PGC-1α ChIP-seq data from BAT39. We found that 50%, or 465 beige-specific peaks overlapped with PGC-1α binding sites, while only 1 white-specific peaks overlapped with PGC-1α binding sites (Fig. 3g and Supplementary Fig. 5e). Notably, Ucp1 upstream has three beige-specific peaks which are PGC-1α binding sites, and overlapped with H3K27ac marks (Fig. 3b). In addition, Pdk4, Slc4a4, and Ppara also have three beige-specific peaks that have binding motifs for PGC-1α (Fig. 3h). These genes themselves exhibited high accessibility and expression level in beige adipocytes (Supplementary Fig. 5f–h). These results strongly suggest that PGC-1α, in concert with the TF motifs enriched in beige-specific peaks identified in the present study play important roles in beige adipocyte function.

To identify potential target genes of beige-specific TFs and their coactivator PGC-1α, we identified the genes satisfying four criteria: (1) Tight linkage to beige-specific peaks, (2) having binding motifs for TFs enriched in beige adipocytes obtained from JASPAR202040 or cisbp database41, (3) have validated binding sites for ESRRα42, MEF2α43, and PGC-1α39 based on published ChIP-seq data, (4) overlap with H3K27ac27 marks in beige adipocytes. Among genes satisfying all criteria, Pdk4 was identified as a target gene for beige-specific TFs and their coactivator (Fig. 3h–j and Supplementary Fig. 5i). The binding sites for PGC-1α, ESRRs, MEF2s, NR2s, and NR4s were all enriched at beige-specific peaks at the locus; accordingly, Pdk4 is a nexus for their coordinated regulation (Fig. 3h–k and Supplementary Fig. 5i). PDK4 enzyme inactivates the pyruvate dehydrogenase complex (PDH), resulting in a decrease in the formation of acetyl-CoA from glucose and a metabolic shift from glucose oxidation to FA oxidation44. The recent loss-of-function study revealed that FA oxidation was blunted in Pdk4 knockdown human adipose-derived stem cells45. Taken together, these data reveal the GRNs including TFs and cis-regulatory elements and their potential target genes underlying beiging process.

Transcription factor modules associated with beige adipocyte developmental trajectories

Our finding that the TFs function coordinately at beige-specific peaks motivated us to identify TF modules for beige adipocyte development. Three TF modules were revealed by unsupervised hierarchical clustering based on the similarity of TF activity computed by chromVAR46 across all individual adipocytes (Fig. 4a). Module 1 contained NR3C, TEAD, and STAT family motifs. It was previously shown that knockdown of TEAD4 increased expression of PPARg and C/EBPα, master regulators for adipogenesis47. Also, adipocyte specific STAT1 knockout promoted PGC-1α expression and enhanced mitochondrial function in WAT depots48. Therefore, the module 1 may represent negative regulators in beige adipogenic process. Module 3 consisted of ROR, NR2, NR4, ESRR, and MEF family motifs, which were enriched in beige adipocytes. The module 3 transcriptional regulators play critical roles in mitochondrial biogenesis and thermogenesis by interacting with PGC-1α in beige adipocytes28,36,37,38. Lastly, module 2 contains C/EBP, RXR, PPAR, and NFI family motifs. RXRs play an important role in Ucp1 induction in brown adipocytes49. RXRs can form heterodimers with PPARs50. PPARa regulates adipogenesis and FA oxidation51. C/EBPs function with PPARg to promote adipocyte development52. Also, it was reported that NFIA binds to brown fat-specific enhancers prior to brown-fat cell differentiation along with C/EBPβ in WAT and BAT53. The motifs in the module 2 have relatively higher correlations with both module 1 and module 3, suggesting intermediate states or shared features between white and beige adipocytes for adipocyte development.

Fig. 4: Uncovering transcription factor modules and beige adipocyte developmental trajectory.
figure 4

a Hierarchical clustering of top variable TF motifs (chromVAR variance > =  1.8) based on their activity scores across individual adipocytes. b Adipocyte trajectory. The smoothed arrow represents pseudotime trajectory in the UMAP embedding. c Motif enrichment dynamics along adipocyte pseudotime trajectory. d Heatmap of gene accessibility along adipocyte pseudotime trajectory. Top variable genes are labeled on the right side of the heatmap. e Gene accessibility dynamics along adipocyte pseudotime trajectory. f, g Side-by-side heatmaps of gene accessibility scores f and motif deviation scores g for TFs, for which the inferred gene accessibility is positively correlated with the chromVAR TF deviation across adipocyte pseudotime trajectory. h, i Gene accessibility of Esrra h and motif deviation of ESRRA (i) on the UMAP embedding. j, Model for the chromatin state dynamics at TF binding motifs through beige adipocyte development. NR3C in TF module 1 might bind to white adipocyte-selective promoters or enhancers to maintain chromatin states of white adipocytes. NFI, RXR, PPAR, and C/EBP in TF module 2 might promote adipocyte development in both white and beige adipocytes. In beige adipocyte, MEF and ESRR in TF module 3 with PGC-1α might promote transcription of beige-selective genes in a coordinated manner.

Given these results, and the facts that (1) individual cells can move between cell states in the same cell type, and (2) cells may exist in transition between states54, we sought to find a path covering the heterogeneous cellular states in the context of continuous beige adipocyte developmental trajectory. We constructed pseudotime trajectory of adipocytes. The inferred trajectory passes from one point where most of day 0 cells are located and to another where most of day 7 cells are located on UMAP (Fig. 4b). In line with our TF module analysis, the activity of TF motifs in module 3 was enriched in the later stage of trajectory, whereas the activity of TF motifs in module 1 was depleted (Fig. 4c). The activity of C/EBPβ and NFIA in module 2 were not significantly changed or gradually increased respectively (Fig. 4c). We ordered the genes based on changes in gene accessibility over pseudotime trajectory (Fig. 4d and Supplementary Fig. 6a, b). It showed that adipocytes lose gene accessibility of Med16, Nnat, Car3, and Fads3 at early stage (Fig. 4d, e), and there is a gradual gain of accessibility for Acly, Me1, and Acaca at the intermediate stage (Fig. 4d). As we expected, Ucp1 was highly accessible at the later stage (Fig. 4d, e). Additionally, Otop1, Cidea, Ppargc1b, Elovl6, Slc4a4, Ppara, and Pank1 were highly accessible in later pseudotime (Fig. 4d).

Finally, we identified positively correlated genes and TFs across pseudotime during beige adipocyte development. The cells at the later developmental stage gained accessibility at Ppara, Esrra, and Mef2d genes and their corresponding TF motifs (Fig. 4f, g). Specifically, the cells with high Esrra gene accessibility also had high ESRRA motif accessibility (Fig. 4h, i). Collectively, these emerging patterns of TF-encoding gene and TF motif accessibility across the adipocyte trajectory define mechanisms underlying environmentally-induced and transcriptionally-controlled adipocyte beiging (Fig. 4j).

Adipose progenitors and vascular network expansion in inguinal adipose tissue

Beige adipocytes can arise via transdifferentiation of unilocular mature adipocytes4,55 but also by de novo differentiation from various progenitors6,10,56. In our snATAC-seq data, two subtypes of APCs were identified (Fig. 5a). Pdgfra, a marker for APCs, was accessible in both APC clusters (Fig. 5b). APC1 was uniquely marked by accessibility at Ebf2, Cd34, and Ly6a (Sca1), whereas APC2 had higher accessibility at Pdgfrb and mural cell marker genes including Acta2 and Myh11 (Fig. 5b).

Fig. 5: Cell–cell interaction for vasculature development.
figure 5

a UMAP plot of adipocyte progenitor and vascular endothelial and smooth muscle cells (n = 1,861). b Heatmap of normalized gene accessibility for cell-type markers. c GO analysis using the top 100 differentially accessible genes for each cell type. d Circle plot displaying VEGF signaling network inferred by gene accessibility of ligand-receptor pairs. Circle sizes are proportional to the size of each group. The circle and line colors represent the senders. The arrows indicate the receivers that receive signal from the corresponding senders. Thicker line indicates a stronger signal interaction. e Major contributors of ligand-receptor pairs of VEGF signaling network among the senders and the receiver. APC1 Adipocyte progenitor cell 1, APC2 Adipocyte progenitor cell 2, VEC Vascular endothelial cell, VSMC Vascular smooth muscle cell.

APCs reside in a vascular niche, and thus, the development of APCs is closely related to the vasculature6,10,57. In our snATAC-seq data, two stromal cell types were identified, vascular endothelial cells (VECs) and vascular smooth muscle cells (VSMCs). VECs were uniquely marked by accessibility at endothelial cell marker genes such as Pecam1 and Chd5, while VSMCs had higher accessibility at smooth muscle cell marker genes including Mylk4, Myh4, Myog, Myf5, and Myod1 (Fig. 5b). We found that the abundance of VECs and VSMCs tended to increase following cold and CL treatment (Fig. 1d). Since they are two essential cell types in blood vessels, the expansion of VEC and VSMC populations is closely associated with vessel growth. In line with this, the genes enriched in VECs were related to vasculature development (Fig. 5c).

To gain insight into how the expansion of vasculature is regulated in iWAT, we used ligand-receptor pairs to infer cellular communication58. This analysis revealed that adipocytes are likely to interact with VECs through vascular endothelial growth factor (VEGF) signaling, mainly through Vegfra and Vegfr1 (Fig. 5d, e). VEGF is a major growth factor that promotes the proliferation of endothelial cells and the growth of new blood vessel in adipose tissue59. WAT vascular expansion is important for formation and maintenance of APC niche7. These results suggest that vascular cells and adipocytes in iWAT form a reciprocal feedback loop to promote angiogenesis in a paracrine manner and expansion of vascular cells provides a niche for APCs.

Decreased abundance of inflammatory immune cells after cold and CL treatment

Immune cells are another critical constituent of the adipose microenvironment that influence adipose tissue function and metabolic homeostasis60. In our snATAC-seq data, we found lymphocytes, including B cells, CD4+ and CD8+ T cells, and myeloid cells, including macrophages and dendritic cells (Fig. 6a, b). We observed a significant decrease in cell abundance of myeloid origin, especially dendritic cells after cold and CL treatment (Fig. 1d). Compared to other immune cells, the motifs for NF-kB family, including RELA, REL, NF-kB1, and NF-kB2, were highly enriched in dendritic cells (Fig. 6c, d). NF-kB is a key mediator of cytokine signaling and NF-kB signaling activates transcription of various proinflammatory genes61. Consistently, GO analysis showed that the less accessible genes after cold and CL treatment in dendritic cells are associated with NF-kB signaling pathway (Fig. 6e). These results suggest that dendritic cells in adipose tissue contribute to proinflammatory microenvironment, mainly through NF-kB signaling, and thus, a decrease in dendritic cell population prevents activation of inflammatory responses after cold exposure and CL treatment.

Fig. 6: Alterations in pro-inflammatory immune cell compartment.
figure 6

a UMAP of immune cell populations (n = 23,687). b Heatmap of scaled gene accessibility for cell-type markers. c UMAP of RELA motif enrichment. d Motif sequences enriched in dendritic cell population. e GO analysis of genes enriched in dendritic cells at day 0 RT. X-axis indicates number of genes. Color indicates adjusted p-value. Macro Macrophage, Dend Dendritic cell.

Shared and distinct mechanisms of cold exposure and CL responses in adipocytes

It has been suggested that beige adipocytes produced by different thermogenic stimuli may have unique genetic signatures and different metabolic properties6. Therefore, we examined similarities and differences of activated gene programs resulting from cold and CL treatments. We first identified genes that were more accessible after 1, 3, and 7 days of cold or CL treatment compared to day 0 (Fig. 7a and Supplementary Fig. 7a–i). We next identified those genes that had elevated accessibility in both treatment groups, or in just one treatment group at the 3 or 7 day time points in order to identify shared and distinct mechanisms (Fig. 7b and Supplementary Fig. 7j–l). GO analysis showed that the genes enriched after both cold and CL treatment were associated with ribose phosphate metabolic process, FA and acyl-CoA metabolic process, and adaptive thermogenesis (Fig. 7c). The genes more enriched after cold exposure were involved in glycolytic process (Fig. 7d); whereas the genes enriched in CL were related to GTPase activity and cAMP signaling (Fig. 7e). These findings are consistent with the facts that Adrb3 activation, like other G protein-coupled receptor signaling cascades, triggers cAMP signaling pathways; and in the case of Adrb3 activation, this leads to induction of thermogenesis62.

Fig. 7: Shifts in lipogenic and glycolytic genes and lipid composition in response to different thermogenic stimuli.
figure 7

a Venn diagram summary showing more accessible genes at day 1, 3, and 7 after cold exposure or CL treatment compared to Day 0 RT. The genes that are commonly more accessible after cold exposure and CL treatment from day 1 to day 7 are shown in the box. Bold font indicates the genes involved in de novo lipogenesis (DNL). b Venn diagram summary showing commonly or differentially more accessible genes after 3 or 7 days of cold exposure and CL treatment compared to Day 0 RT. c GO analysis of genes that were commonly more accessible in both cold and CL. d, e GO analysis of genes that are more accessible only in cold d or CL e. X-axis indicates number of genes. Color indicates adjusted p-value. f, g Violin and box plots showing the distribution of average gene accessibility of lipogenic genes (Acaca, Acacb, Acly, Acss2, Fasn, Scd1, Elovl6) f and glycolytic genes (Hk1, Gpi1, Pfkl, Pfkp, Gapdh, Pgk1, Pgm1, Eno1, Pkm, Ldha) (g) of fold changes after cold exposure and CL treatment compared to Day0 RT. The lower and upper hinges indicate the 25th and 75th percentiles. The horizontal line in the middle of box plot denotes the median. The mean for each day for cold and CL are connected by a line. *p-value<0.05, Welch’s two sample t-test was performed. h Relative concentration of quantified lipid species after cold exposure and CL treatment (n = 3 per group; More lipid species in Supplementary Fig. 8c). i Ratio of desaturation and elongation from lipidomics data (n = 3 per group). j, Relative proportion of saturated FA (SFA), monounsaturated FA (MUFA), and polyunsaturated FA (PUFA) by lipidomics analysis (n = 3 per group). k Ratio of unsaturated FA (UFA) to SFA (n = 3 per group). The slope for each cold and CL from day 1 to day 7 is calculated by a linear regression function. The p-value for the comparison of the slopes between cold and CL is indicated. l, Model for the contribution of cold temperature and CL treatment to FA metabolic process through functional accessibility changes at lipidomic genes and changes in lipid profile. *Adjusted p-value<0.05, ANOVA multiple comparisons test with Bonferroni’s post-hoc test was performed for hj.

The genes that are commonly more accessible by cold and CL treatment for 1, 3, and 7 days were involved in FA metabolism, specifically de novo lipogenesis (DNL) (Fig. 7a). Interestingly, we found some kinetic differences in accessibility of lipogenic genes between cold and CL across time course (Fig. 7f). In cold, the accessibility of lipogenic genes gradually increased from day 1 to day 7, whereas the accessibility of lipogenic genes was highly increased at day 3 in CL. In contrast, the accessibility of glycolytic genes gradually increased after cold exposure, while no significant changes were observed after CL treatment (Fig. 7g). Therefore, cold and CL share common pathways during beiging, yet differ in specific signaling cascades and in timing of changes in gene accessibility.

Cold and CL alter the composition of lipid classes in inguinal adipose tissue

The genes encoding enzymes involved in DNL and FA elongation such as Acaca, Acly, Fasn, and Elovl6 underwent profound changes in adipocytes after cold exposure and CL treatment (Fig. 7a and Supplementary Fig. 8a, b). FA synthetase (FAS), encoded by Fasn, converts malonyl-CoA into palmitic acid, a major product of DNL63. Then, FAs can be further elongated by FA elongases (ELOVLs)63. If any of these changes were functionally important, we expect that FAs extracted from adipose tissue would exhibit different profiles from RT after cold or CL treatment, or that profiles found in cold and CL may differ as well. Therefore, we measured a total of 21 lipid species from the same iWAT used for snATAC-seq by gas-chromatography mass spectrometry (GC-MS)-based lipid profiling. The most abundant free FAs were oleic acid (18:1n-9), linoleic acid (18:2n-6), and palmitic acid (16:0), making up around 82% of the FA content of the depot (Fig. 7h and Supplementary Fig. 8c). Palmitoleic acid (16:1n-7), stearic acid (18:0), and vaccenic acid (18:1n-7) were the next most abundant free FAs, which account for 10% (Fig. 7h and Supplementary Fig. 8c).

Based on the composition of major FAs and its elongation products, we calculated the metabolic activity of elongation and found that both cold and CL significantly increased the ratio of elongation (Fig. 7i). Consistently, our snATAC-seq analysis showed changes in FA elongation and accessibility of Elovl6, a gene encoding a key enzyme that catalyzes the elongation of FAs (Figs. 2h and 7a), providing evidence for the correlation between chromatin accessibility and their functional modulation of lipid profiles. Next, we characterized changes in lipid composition after cold and CL treatment. Despite the increased elongation by both cold and CL, different types of FAs were increased by cold and CL. Strikingly, CL treatment significantly increased the proportion of oleic acid, a major component of monounsaturated FA (MUFA), not cold treatment (Fig. 7h, j). In addition, the proportion of some lipid species have opposite dynamics between cold and CL. For example, cold exposure tended to decrease palmitic acid gradually over time course, while CL treatment dramatically decreased palmitic acid at day 1 then restored to the day 0 level by day 7 of CL treatment (Fig. 7h). Accordingly, the ratio of unsaturated FA (UFA) to saturated FA (SFA) had opposite trends between cold and CL over the time course of our study (Fig. 7k). Overall, these results reveal that both cold and CL treatment increase abundance of longer chain FAs in adipose, but they increase different types of UFA in adipose with the different dynamics (Fig. 7l).