Table of Contents
- 1 snATAC-seq reveals distinct cell types in inguinal adipose tissue
- 2 Dynamic changes in adipocytes and identification of thermogenic beige adipocytes
- 3 Beige-specific cis– and trans-acting gene regulatory networks
- 4 Transcription factor modules associated with beige adipocyte developmental trajectories
- 5 Adipose progenitors and vascular network expansion in inguinal adipose tissue
- 6 Decreased abundance of inflammatory immune cells after cold and CL treatment
- 7 Shared and distinct mechanisms of cold exposure and CL responses in adipocytes
- 8 Cold and CL alter the composition of lipid classes in inguinal adipose tissue
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).
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).
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.
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.
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).
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.
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.
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).