Extensive enhancer crosstalk controls PPARG2 activation during adipogenesis

Domain: biology
Published: 2026-03-15
Editorial outcome: minor_revision

Abstract
Transcriptional master regulators drive cell fate transitions. Peroxisome proliferator-activated receptor γ (PPARγ) is the master regulator of adipogenesis, and its expression must therefore be tightly regulated and efficiently induced in response to adipogenic cues. Here we decipher the regulatory mechanisms of the highly connected enhancer community driving activation of the PPARG locus during adipocyte differentiation of human mesenchymal stem cells. By systematically deleting 9 individual enhancers, spanning upstream, promoter-proximal, and downstream super-enhancer constituents, we demonstrate elaborate enhancer crosstalk in cis involving stabilization of C/EBPβ recruitment prior to chromatin remodeling. We show that the super-enhancer constituent E+102 plays a dual role in cis crosstalk and feedback activation and is obligate for activation of PPARG expression. Non-coding genetic variants associated with cardiometabolic traits and predicted to regulate PPARG expression map E+102 and other essential enhancers in the community, thereby supporting the importance of these enhancers in human physiology and disease.

Introduction
Adipocytes are highly specialized cells that play an important role in metabolic energy storage and homeostasis. Genetic or acquired conditions resulting in decreased lipid storage and endocrine function, so-called partial lipodystrophies, are closely associated with ectopic lipid deposition and cardiometabolic diseases1. Similarly, mounting evidence points to compromised adipocyte function as a key driver of cardiometabolic co-morbidities in obesity2. Adipogenesis is the process by which stromal progenitor cells of the mesenchymal lineage develop into lipid-laden adipocytes3,4. This process is of major biomedical importance, as it is required for the development of adipose tissue, as well as for healthy expansion of adipose tissue5. The importance of this process, and the fact that it can be efficiently induced in vitro, has made it one of the most well-studied differentiation model systems.

Studies of many individual genes, as well as genome-wide transcriptional and epigenomic studies, have documented that the transcriptional networks regulating in vitro adipogenesis consist of at least two consecutive waves of transcription factors (TFs). The first wave is induced directly by the adipogenic inducers and drives the remodelling of the chromatin template as well as the activation of the master regulators peroxisome proliferator-activated receptor γ (PPARγ) and CCAAT/enhancer-binding protein α (C/EBPα). These master regulators are key drivers of the second wave, which activates the adipocyte gene program. The most important of the two master regulators is PPARγ, which is indispensable for adipocyte differentiation in vitro as well as in vivo and sufficient to drive trans-differentiation of other mesenchymal lineages into adipocytes3,6,7. The master regulator attributes of PPARγ rely at least in part on its biochemical properties, as suggested by comparative studies of lipodystrophy mutants of PPARγ8 as well as comparative studies with other members of the PPAR family9,10. In addition, the ability to dramatically increase PPARG expression in response to adipogenic cues may play a role. Thus, the expression of PPARG must be tightly controlled; initially repressed in the progenitor state, and then rapidly induced in response to the transcription factors of the first wave of adipogenesis. Several transcription factors, most notably C/EBPβ and C/EBPδ, have been shown to be directly involved in the activation of the PPARG locus11,12. Surprisingly, however, the mechanisms by which the PPARG locus is activated and the transcription of PPARG is switched on during adipogenesis remain poorly understood.

Adipogenesis is associated with dramatic chromatin remodeling and de novo enhancer activation12-14 and dynamic rewiring of promoter-enhancer interactions (Siersbæk et al., 2017). Using enhancer capture Hi-C (ECHi-C), we identified extensive 3D interconnectivity of enhancers in multipotent human stromal progenitors (hMSC-TERAT4), forming highly interconnected enhancer (HICE) communities that are dynamically regulated during adipogenesis15. Interestingly, lineage-selective promoters are generally associated with the most highly connected HICE communities, and these communities are typically established already at the stem cell state, indicating that HICE communities may preprogram lineage fates in progenitor cells. In response to the adipogenic cocktail, the connectivity of the adipogenic HICE communities increased correlating with increased expression of the associated genes15. This suggests that enhancer connectivity plays an important role in lineage determination possibly by mediating cooperativity between enhancers.

Several other studies have reported cooperativity between enhancers; however, the molecular mechanisms underlying this cooperativity is unclear. Depending on the context, enhancers have been reported to display redundancy, additive or synergistic cooperativity16-23. This emphasizes the importance of studying the cooperativity and crosstalk between enhancers in their natural context using detailed molecular analyses of endogenous loci. Currently, such studies are primarily carried out in embryonic stem cells22 and erythroid cells16,19.

In this work we have dissected the enhancer community that drives the induction of the most important locus in adipogenesis, that of PPARG, in hMSC-TERT4 cells. Using genome editing, we identify key regulatory enhancer constituents in the PPARG community and demonstrate extensive enhancer crosstalk and interdependencies. Importantly, we show that the key enhancers overlap with cardiometabolic disease-associated genetic variants, predicted to regulate PPARG2 transcription, thereby underscoring the in vivo relevance of these enhancers. Our findings offer a conceptual framework for how enhancer crosstalk governs cell fate decision, which may be disrupted in a disease context.

Results
The PPARG locus is primed and connected in undifferentiated human MSC
Taking advantage of our previously published Hi-C data from undifferentiated human bone marrow-derived mesenchymal stem cells (hMSC-TERT4)15, we inspected the overall chromatin topology around the human PPARG locus. This showed that PPARG is isolated in a topologically associating domain (TAD) of 350 kb (chr3: 12,174,868 - 12,526,858) (Fig. 1a) at the stem cell stage prior to chromatin remodeling and induction of PPARG expression13. Notably, this TAD also contains TIMP4 and the 3’end of SYN2 with no known functions in adipogenesis. Consistent with the TAD structure, ChIP-seq revealed strong CTCF binding with divergent motifs at the boundaries of the PPARG TAD (Fig. 1a). Analysis of publicly available Hi-C data from other stem cells24 indicated that the locus is also isolated into a partially defined TAD in embryonic (ESC), mesendodermal cells and mesenchymal stem cells (MSC), but not in neuronal progenitor cells, indicating that the locus is organized very early during development but loses its connectivity in some lineages (Supplementary Fig. 1a).

Inspection of the PPARG locus connectivity at higher resolution using our previously published ECHi-C data from hMSC-TERT4 cells15 confirmed that the locus is engaged in an isolated chromatin network in the undifferentiated state and becomes even more connected and confined during differentiation (Fig. 1b). ChIP-seq analysis revealed that the locus is already primed by H3K4me1 in the undifferentiated state but only gains the active enhancer mark H3K27ac by day 3 of differentiation, coinciding with the activation of PPARG2 expression (Fig. 1a, Fig. 1c).

The PPARG locus is transcribed from two major transcription start sites (TSS), PPARG1 and PPARG2, where PPARG2 is the adipocyte-specific TSS25. The PPARG2 transcript encodes PPARγ2, which is characterized by an additional 28 amino acids in the N-terminal compared with PPARγ1. PPARγ2 has been shown to be the most adipogenic isoform26,27 and to be required for efficient adipocyte differentiation in vitro and in vivo in mice28-30. Notably, in hMSC-TERT4 cells, the PPARG1 promoter is inactive (Supplementary Fig. 1b-c), thereby enabling us to specifically study the activation of the PPARG2 promoter without interference from the PPARG1 promoter. Similarly to PPARG2, the two other genes located in the TAD, SYN2 and TIMP4, are selectively induced during adipocyte differentiation, whereas the expression of genes outside the TAD borders, TAMM41 and TSEN2, are not affected by differentiation (Fig. 1c).

Collectively, these results show that the PPARG2 enhancer community is demarcated prior to cellular commitment and then activated in response to the adipogenic cocktail. Importantly, the regulatory functions of these enhancers appear to be contained within the TAD.

The connectivity and activity of the PPARG enhancer community is increased during early adipogenesis
Integration of ECHi-C, MED1 ChIP-seq, and DNase-seq data obtained during adipogenesis and osteogenesis of hMSC-TERT4 cells13,15 revealed marked activation of several putative enhancers in the PPARG locus between 1 and 3 days of adipocyte differentiation (Fig. 2a) and showed that these enhancers are part of a HICE community (Fig. 2a top panel). To investigate the mechanistic workings of the PPARG enhancer community during adipogenesis, we selected 9 putative enhancers that exhibit three different temporal activation patterns during differentiation, yet all are engaged in the PPARG enhancer community already in the undifferentiated state. Enhancers were named according to their position (in kb) relative to the PPARG2 TSS. Two distal upstream enhancers (E-130, E-123) are already active in undifferentiated cells and remain active throughout both adipocyte and osteoblast differentiation. In contrast, the two promoter-proximal enhancers (E-1, E+2) are activated de novo by day 1 of adipogenesis, and the 5 distal downstream enhancers (E+89, E+102, E+107, E+111, E+125), collectively classified as a super-enhancer by the ROSE algorithm31,32 (Supplementary Fig. 2b), are activated de novo by day 1-3 (Fig. 2a, Supplementary Fig. 2a). During the time course of adipogenesis, none of the nine enhancers bind CTCF, except for E+107, where CTCF binding is observed from day 1 of adipogenesis (Fig. 2a).

To specifically investigate the temporal dynamics of chromatin interactions in the PPARG enhancer community during differentiation, we performed Capture-C in the progenitor state and on day 10 of adipogenesis using the nine putative enhancers and the PPARG2 TSS as bait. We focused on pairwise interactions within the specific region of the PPARG locus (chr3: 12.1-12.6 Mb), excluding interactions occurring within ± 2kb of a viewpoint. To compare interactions across samples, we normalized the unique interaction counts within the region (cis interactions) to 50,000 counts per sample and viewpoint (Supplementary Fig. 3). In line with the dynamics of the ECHi-C-based interactions seen in the PPARG2 HICE community (Fig. 2a), Capture-C showed that the main long-distance interaction in the undifferentiated state is between the upstream enhancer E-123 and the super-enhancer constituent E+125, which thereby “pinches” the locus (Fig. 2b-d, Supplementary Fig. 4a). In addition, there are strong interactions between the two upstream enhancers E-130 and E-123 (Fig. 2c).

Interestingly, we also observed strong connections between the super-enhancer constituents, which are primed but not yet remodeled or active (Fig. 2c). The most dramatic changes in connectivity during adipogenesis are the strong inductions of interactions between the PPARG2 promoter-proximal regions and the downstream super-enhancer constituents (Fig. 2b-c).

Taken together, the connectivity of the PPARG enhancer community is markedly increased during adipogenesis, concomitantly with the sequential activation of putative enhancers near the TSS as well as in the downstream super-enhancer region (Fig. 2d).

PPARG2 activation requires cooperativity between many non-redundant enhancers
To investigate the role of individual enhancers in chromatin topology and activation of PPARG2 transcription, we generated a hMSC-TERT4-based clonal cell line with inducible Cas9 (iCas9; hMSC-TERT4/iCas9-4) (Fig. 3a, Supplementary Fig. 5a). hMSC-TERT4 cells were transduced with Lenti-iCas9-Neo33 encoding a flag-tagged human codon-optimized S. pyogenes Cas9 coupled to green fluorescent protein (GFP) by a Porcine teschovirus-1 2A (P2A) self-cleavage peptide (Supplementary Fig. 5b-c). In this system, the expression of GFP-coupled Cas9 can be transiently induced by doxycycline (Supplementary Fig. 5b-d), thereby minimizing off-target effects from constitutive Cas9 expression.

The differentiation potential of the selected hMSC-TERT4/iCas9-4 clone is comparable to wild-type hMSC-TERT4 cells (Supplementary Fig. 5e). The hMSC-TERT4/iCas9-4 cell line was transfected with pairs of gRNAs designed based on the position and width of the corresponding DNase-seq peak to excise each of the putative enhancers (Fig. 3a, Fig. S9, Supplementary Table 6). For each deletion, we selected three independent clones. As controls, we used three independent clones transfected with control gRNA (‘Ctr gRNA’) targeting intron 1 of the PPP1R12C gene, which is suggested to be a “safe harbor” for incorporation of adeno-associated virus and thus genetic modifications34. For each clone, the zygosity and the size of the deletion were confirmed by Sanger sequencing of regions flanking the deleted sites (Supplementary Table 6, Supplementary Table 3).

The control and enhancer deletion clones were induced to undergo adipocyte differentiation, RNA was harvested on day 10 and used for quantification of transcripts by qPCR. Interestingly, deletion of several individual enhancers, i.e., the promoter-proximal enhancer E-1 and E+2, as well as the super-enhancer constituents E+102, E+107, and E+111, results in a significantly decreased expression of PPARG2 and SYN2, with the deletion of E+102 almost completely preventing the induction of PPARG2 (Fig. 3b). Deletion of the distal upstream E-123 is associated with a slight decrease in PPARG2 expression, whereas deletion of the distal upstream E-130 and distal downstream E+125 enhancers do not affect PPARG2 (Fig. 3b). These results indicate an important and non-redundant role of the promoter-proximal enhancers E-1 and E+2, as well as the super-enhancer constituents E+102, E+107, and E+111. Inspection of epigenetic profiles of these enhancers shows that they are all remodeled and activated (Fig. 2a) and become much more interconnected (Fig. 2c) during adipogenesis. Furthermore, during adipogenesis, all super-enhancer constituents, E+102 in particular, become more connected to the PPARG2 TSS (Fig. 2b).

Notably, none of the enhancer deletions affect the expression of TAMM41 and TSEN2, which are located outside the TAD, or the expression of TIMP4 in the TAD (Fig. 3b). By contrast, the expressions of SYN2 and PPARG2 are similarly affected by the individual enhancer deletions (Fig. 3b). The canonical SYN2 promoter region is outside the TAD and does not seem to interact with the PPARG community enhancers on day 10 of adipogenesis (Fig. 1b), suggesting SYN2 might be transcribed from an alternative TSS. Moreover, both controls, ‘No gRNA’ and ‘Ctr gRNA’, exhibit comparable levels of expression of PPARG2 and neighboring genes (Fig. 3b) and therefore, for subsequent experiments, we proceeded with ‘No gRNA’ as the control, unless otherwise stated.

Consistent with PPARγ being a master regulator of adipogenesis and a direct activator of most genes involved in lipid handling and metabolism in adipocytes6, there is a strong correlation between PPARG2 expression and lipid accumulation in the different enhancer-deleted clones at day 9 of adipocyte differentiation (Pearson R=0.99, Fig. 3c, Supplementary Fig. 5f). Moreover, deletions of enhancers that are important for the expression of PPARG2 also dramatically reduced the expression of late adipocyte genes, such the PPARγ target genes FABP4, CEBPA, but not the early adipogenic transcription factor, CEBPB, acting upstream of PPARγ (Supplementary Fig. 5g). In summary, at least six out of nine putative enhancers in the PPARG HICE community are required for the full induction of PPARG2 expression and ultimately adipocyte differentiation.

Deletion of key enhancers affects the connectivity in the PPARG locus
To explore the role of individual enhancers in the development of local topology in the PPARG locus, we performed Capture-C in controls and enhancer-deleted clones before induction of differentiation (day 0) and following exposure to the adipogenic cocktail for 10 days. Results from ‘No gRNA’ and ‘Ctr gRNA’ controls correlated well (Supplementary Fig. 4b) and were therefore considered as combined controls (n=4) for all statistical tests. Interactions where one of the interacting elements was deleted were disregarded. In the progenitor cell state, where neither the enhancers in the promoter-proximal region nor the super-enhancer constituents are accessible, as assessed by DNase-seq (Fig. 2a), or active, as assessed by H3K27ac and MED1 ChIP-seq (Fig. 2a, Supplementary Fig. 2a), only a few enhancer interactions are affected by the deletion of other enhancers (Fig. 4a). Interestingly, however, the deletion of E+2, E+107 or E+125 leads to decreased interaction between E-1 and E+89, while that between E+102 and the PPARG2 TSS is increased. This indicates that the E+102-TSS connectivity may be limited by other super-enhancer constituents and E+2, possibly by some competitive mechanisms. Overall, these data indicate that the connectivity of the PPARG HICE community at the progenitor stage is relatively resistant to the deletion of constituents. Moreover, the architectural changes that occur upon enhancer deletion at this stage do not correlate well with the ability of the enhancers to induce PPARG2.

At day 10 of adipogenesis, when the community is more connected, many enhancer interactions are affected by the deletion of other enhancers. The most notable changes in enhancer interaction occur in response to the deletion of enhancers that are important for the induction of PPARG2 expression (Fig. 4a), in particular E+102 (Fig. 4b) and E+2 (Fig. 4c). Deletion of E+102 inherently ablates the strongest and most differentiation-dependent interactions between the super-enhancer and the TSS; however, it also leads to a general destabilization of the connectivity between super-enhancer constituents (e.g., the interactions between E+89 and E+107; and between E+107 and E+111) as well as a decrease in the interactions between other super enhancer constituents and the TSS and promoter proximal enhancers (e.g., interactions between E+89 and TSS and E-1; and between E+107 and E+2) (Fig. 4b). Deletion of E+2 leads to a dramatic decrease in the interaction between the super-enhancer constituents and the TSS as well as a decreased connectivity of E+111 with E+107 and with promoter-proximal E-1, respectively (Fig. 4b). Interestingly, however, the deletion of E+2 also increases the interaction frequency between E+102 and E-123 and between E+102 and E-107 (Fig. 4c). This indicates that E+2 may help direct the super-enhancer to interact with the TSS and in doing so may compete with other super-enhancer constituents for interactions with E+102.

Taken together, while many genomic interactions in the PPARG locus are unaffected by the deletion of individual enhancers, deletion of E+102 or E+2 interferes with adipogenesis-associated stabilization of interactions in the PPARG HICE community, specifically between the PPARG2 proximal regions and the distal downstream super-enhancer (Fig. 4d).

Enhancers in the PPARG locus display wide-spread functional crosstalk
To investigate whether the enhancers within the PPARG enhancer community affect the activity of each other, we assessed MED1 occupancy in individual clones at day 10 of adipocyte differentiation, where all enhancers are active in control cells (Fig. 2a). Interestingly, deletion of E+102 leads to decreased MED1 recruitment to other enhancers activated during adipogenesis, i.e. E-1, E+2, E+89, E+107, E+111 and E+125, as assessed by ChIP-PCR (Fig. 5a, Supplementary Fig. 6a) as well as ChIP-seq (Supplementary Fig. 6b). Similar locus-wide effects, but with a smaller magnitude, are observed in response to deletion of other enhancers that are important for PPARG2 expression (Fig. 3b, Fig. 5a). Thus, at day 10 of adipogenesis, MED1 occupancy is not only lost at the deleted enhancer or its closest neighbors but is also lost at other enhancers in the locus. The summed MED1 occupancy at the 9 putative enhancers correlated well (Pearson R=0.83, Fig. 5b) with PPARG2 expression, suggesting that transcription of the locus may be determined by the sum of enhancer activity.

We also assessed MED1 ChIP-seq at day 3 of adipogenesis, the earliest time point with detectable MED1 occupancy at the super-enhancer (Fig. 2a). Similar to MED1 ChIP-qPCR data at day 10 of adipogenesis (Fig. 5a, Fig. 6b), this showed that the most pronounced effect of E+102 deletion is a lower MED1 signal at the other super-enhancer constituents, E+89, E+107, E+111, and E+125 (Fig. 5c), whereas the most pronounced effect of E+2 deletion is a decreased MED1 occupancy at the neighboring promoter-proximal E-1 and the distal E+125 (Fig. 5c). Importantly, the effects of both E+102 and E+2 deletions are restricted to the PPARG enhancer community (Supplementary Fig. 6c-d).

Collectively, these data indicate profound crosstalk between enhancers in the PPARG locus, in particular between those that interact in undifferentiated cells (Fig. 5d).

PPARG2 enhancer network is activated through extensive enhancer crosstalk
We have previously shown that at a genome-wide level, there is a strong positive correlation between the occupancy of a TF at the bait and the target enhancers of two interacting enhancers bound by the same TF. We hypothesized that increased TF occupancy at one enhancer positively affects the occupancy at interacting enhancers and termed the phenomenon ‘cross-interaction stabilization’ (CIST)15. The PPARG enhancer deletion mutants provide a unique opportunity to challenge this hypothesis and determine whether CIST may be involved in cis cooperativity between enhancers in the PPARG locus. C/EBPβ is known to be a key driver of early adipogenesis and has been proposed to directly activate Pparg gene expression in mouse preadipocytes35. Consistent with that, C/EBPβ ChIP-seq analysis showed rapid recruitment of C/EBPβ in response to the adipogenic cocktail (Supplementary Fig. 7a). In the PPARG locus, the distal upstream enhancers, E-130 and E-123, are bound by C/EBPβ in the progenitor state, while the promoter-proximal enhancers, E-1 and E+2 and super-enhancer constituents E+102, E+107, E+111 recruit C/EBPβ at day 1 of differentiation (Fig. 6a). Notably, the binding of C/EBPβ to the super-enhancer at day 1 (Fig. 6a) precedes chromatin remodeling (Fig. 2a), suggesting that C/EBPβ might serve as a pioneering factor in the super-enhancer region. Similarly, C/EBPβ has been shown to bind to chromatin prior to remodeling in the 3T3-L1 preadipocyte cell line14.

To explore whether early C/EBPβ occupancy in the PPARG HICE community is affected by the deletion of other C/EBPβ bound interacting enhancers, we performed C/EBPβ ChIP-seq in control and E+2 and E+102 deletion cell lines at day 1 of adipogenesis (Fig. 6b). Interestingly, and consistent with CIST, deletion of E+102 leads to significantly decreased C/EBPβ occupancy at super-enhancer constituents E+107 and E+125, which are connected to E+102 at the progenitor state (Fig. 6b, Fig. 2c). By contrast, C/EBPβ binding at the distal upstream enhancers, E-130 and E-123, and at the promoter-proximal E-1 and E+2 are not significantly affected (Fig. 6b). Deletion of E+2 decreases C/EBPβ occupancy at distal super-enhancer constituent E+125, the only enhancer that interacts with E+2 at the progenitor state (Fig. 6b, Fig. 2c). Importantly, these effects are restricted to the PPARG enhancer community (Supplementary Fig. 7b), as the C/EBPβ ChIP-seq signal in the E+2 and E+102 deletion cell lines is comparable to control lines outside the PPARG locus (Supplementary Fig. 7b) and at a genome-wide level (Supplementary Fig. 7c). Collectively, this provides experimental support for the CIST model and indicates that C/EBPβ bound enhancers cooperate in cis through interactions and stabilization of C/EBPβ binding (Fig. 6c).

PPARγ has been proposed to positively feedback regulate its own expression through binding to promoter-proximal regulatory elements in the PPARG2 in mice and humans36 as well as through induction of other transcription factors, including C/EBPα, which regulate the PPARG2 promoter35,37. The importance of this feedback regulation is supported by early loss-of-function studies38. In hMSC-TERT4 cells, PPARγ binds to E-1 and E+102 at day 3 of adipogenesis (Fig. 6d), both of which have moderately conserved PPAR response elements (PPREs) (Supplementary Fig. 7d). At later stages of adipogenesis, PPARγ occupancy is also gained at E+107 and E+2, although no consensus PPRE is found at E+2. In summary, three out of the four most important enhancers for regulating PPARG2 expression, MED1 occupancy, and lipid accumulation exhibit direct binding by PPARγ itself. Moreover, C/EBPα ChIP-seq at day 10 of differentiation showed that C/EBPα binds to the promoter-proximal E-1, E+2, and to super-enhancer constituents E+102, E+107, E+111 by day 10 of adipogenesis, indicating that C/EBPα is also involved in a positive feedback loop (Fig. 6d). Thus, a strong enhancer in the PPARG locus may affect the activity of other enhancers not just in cis but also in trans by contributing to increased expression of PPARγ which then either directly activates other PPARG enhancers or the expression of other adipogenic transcription factors (Fig. 6e). Interestingly, we also noticed a strong binding of PPARγ by day 3 near the TSS of a shorter predicted isoform of SYN2 (Fig. 6d), suggesting that the co-regulation of SYN2 and PPARG2 during adipogenesis might be mediated, at least in part, through PPARγ activation of promoter-proximal enhancers in the SYN2 locus.

Taken together, the enhancers in the PPARG HICE community show clear evidence of cis cooperativity between interacting enhancers, including increased C/EBPβ binding at early stages of differentiation prior to expression of PPARG. Once PPARγ is induced, it binds to both promoter-proximal enhancers as well as some super-enhancer constituents thereby potentially enhancing its own expression.

Key PPARG2 enhancers overlap with non-coding genetic variants associated with metabolic disease traits
Given the importance of well-functioning adipocytes for metabolic health and the obligate role of PPARγ in adipocyte differentiation and function2,6, we rationalized that if the enhancers identified in vitro are also important in vivo, they would overlap with metabolic disease-associated genetic variants. We queried the PPARG TAD (chr3: 12,174,868 - 12,526,858) and extracted 796 haplotypes associated with a series of cardiometabolic traits comprising in total 3596 variants in high linkage disequilibrium (LD) (r2> 0.8, CEU) from GWAS in the public domain (Fig. 7a). The vast majority of these genetic variants associated with common complex metabolic traits map to non-coding elements in the PPARG locus. Because any given variant in a haplotype could be functional and, therefore, mediate the association with metabolic disease, we used the Enformer model39, a state-of-the-art sequence-based deep learning approach for the prediction of non-coding variant effects on gene expression and chromatin states. We predicted the variant effects on PPARG2 transcript expression as well as PPARG enhancer regions and ranked the variants based on their positive or negative effects on the target regions (Fig. 7b, Supplementary Data 1, Supplementary Data 2). The analysis revealed that 1759 out of 3596 variants have predicted significant effects on PPARG2 gene expression.

Notably, 11 out of the 40 top-scoring variants with predicted (positive and negative) effects on PPARG2 signal and/or enhancer accessibility, overlap with enhancers identified in the hMSC-TERT4 model system (Fig. 7b-c, Supplementary Data 1, Supplementary Data 2). In particular, many of these genetic variants are located in the three most important enhancers, E+102, E+2, and E-1 (Fig. 7b-c, Supplementary Data 1, Supplementary Data 2), indicating that these enhancers are also important for the regulation of PPARG expression in human physiology. As part of this analysis, we found that 14 of these Enformer top-scored variants are fine-mapped to phenotypes (FinnGen, PIP>0.01) including rs12486170, rs7649970, and rs7647481 located in E-1 (Supplementary Data 3).

These variants may change transcription factor affinity for the enhancer sequence, thereby impacting transcription factors binding, and enhancer activity. One of the top non-coding variants, rs7647481 A/G, located in E-1, is a common variant associated with decreased risk of developing type 2 diabetes (rs7647481 A/G, p=2.07E-77, Supplementary Data 1). Using label-free quantitative LC–MS/MS proteomics, this variant has previously been shown to enhance the binding of YY1 leading to increased activation of the PPARG2 promoter40. Moreover, 8 top-scoring common variants from the Enformer analysis, including rs7647481 and rs7649970 located in E-1, colocalize with eQTLs for PPARG expression (Fig. 7b, Supplementary Data 4), indicating a physiological relevance of these cis-acting elements in the PPARG locus.

Another interesting, low-frequency variant, rs181025382 G/A (MAF=0.0001, CEU) associated with lipid traits including decreased levels of total cholesterol and non-HDL cholesterol resides in E+102 and is predicted to negatively impact both PPARG2 expression and enhancer accessibility level (Fig. 7c). Interestingly, this variant is located in a PPRE in E+102, which is the enhancer first and most strongly bound by PPARγ from day 3 of differentiation (Fig. 6e). In the minor allele variant, the A in the spacer between the two half-sites of the direct repeat is replaced with a less optimal G (Fig. 7d). To estimate the effect of a G in this position on the binding affinity of the PPARγ:RXR heterodimer, we took advantage of the recently published nuclear receptor sequence specific binding preferences for binding to naked DNA in the presence or absence of agonists based on high throughput SELEX data41. In E+102, the minor variant with a G in the spacer showed an approximately 3-fold lower affinity for PPARγ:RXR compared with the major variant, which has an A in that position (Fig. 7d). The synthetic agonist rosiglitazone slightly increased the affinity of PPARγ:RXR for both PPRE variants (Fig. 7d), consistent with previous ChIP-seq data42. To investigate the biological relevance of a 3-fold difference in binding affinity for the PPARγ:RXR heterodimer, we compared PPARγ:RXR affinity (± rosiglitazone) between the 1,000 top PPARγ-bound enhancers and 1000 PPARγ-unbound enhancers from our datasets derived from hMSC-TERT4 (Supplementary Fig. 8a). Importantly, PPARγ generally exhibits about 3-fold higher affinity for bound sites compared to unbound enhancers, with affinity being higher for rosiglitazone-activated PPARγ (Supplementary Fig. 7a). This indicates that a 3-fold difference in naked DNA affinity may have major functional implications for binding of PPARγ to sites in chromatin.

To further interrogate PPARγ binding to chromatin in hMSC-TERT4 cells, we applied BPNet43, a deep learning model to predict base-resolution binding syntax from DNA sequences within PPARγ binding peaks in hMSC-TERT4 (Fig. 7e). We used in silico mutagenesis to investigate the effect of rs181025382 on PPARγ:RXR binding in E+102 and found that substitution of the major allele A with minor allele G leads to -0.12 log fold change (Fig. 7e).

Further, we used TF-MoDIScO for motif discovery within identified seqlets in PPARγ-bound sequences. By inferring the data with previously annotated motifs from JASPAR database, we identified PPARγ:RXR as a top de novo motif matching with PPARγ:RXR (MA0065.2) motif (FDR= 1.3e-10, Supplementary Table 7) (Fig. 7f). To ensure the predicted decrease in binding upon A to G substitution was not specific to E+102, we identified PPARγ binding sites genome-wide that contain the PPARγ motif with an A at the 9th position. At all these sites, we mutated the A to a G and predicted the difference in binding using the trained BPNet model. This showed -0.09 log fold change in PPARγ:RXR binding affinity genome-wide upon substitution of A with G, similar to results for rs181025282 (Supplementary Fig. 8b-c). This suggests that the G variant of rs181025382, is likely to interfere with the ability of PPARγ to activate this enhancer, thereby reducing the positive feedback activation of PPARγ on PPARG2 expression.

Our results show that PPARG enhancers that are critical for in vitro activation of PPARG2 in response to adipogenic inducers overlap with functional non-coding genetic variants associated with metabolic disease. The genetic variants in these enhancers are among the top-scoring genetic variants predicted to affect PPARG2 expression. These results indicate that the enhancers identified in our hMSC-TERT4 model system are important for human physiology. Furthermore, the finding that one of these notably low-frequency, non-coding genetic variants decreases the binding affinity of PPARγ:RXR for the most important enhancer constituent E+102 in the super-enhancer further supports the notion that positive feedback activation of the locus plays an important role in the activation of PPARG2 expression during adipogenesis, a function relevant to cardiometabolic disease in humans.

Methods
Cell culture and differentiation
Telomerase-immortalized human bone marrow-derived stromal progenitor cells (hMSC-TERT4) were kindly provided by Moustapha Kassem49,50, and maintained in α-MEM with 10 % fetal bovine serum (FBS) and 1% penicillin-streptomycin (P/S) (growth medium) at 37°C, 95% humidity and 5% carbon dioxide. Osteo- or adipogenic differentiation was induced two days post confluence (day 0) by changing the growth medium to the differentiation medium. Osteogenic (α-MEM with 10% FBS, 1% P/S and 10 nM freshly added vitamin D3 in 96% ethanol, 50 μg/mL L-ascorbic acid in MilliQ water, 10 mM β-glycerophosphate in MilliQ water and 10 nM dexamethasone) or adipogenic differentiation medium (DMEM with 10% FBS, 1% P/S and 100nM freshly added dexamethasone in 96% ethanol, 0.5 mM 3-isobutyl-1-methylxanthine in 0.1 M KOH, 2 μM human insulin, and 1 μM rosiglitazone in DMSO) was renewed every second day until cells were harvested. HEK293T cells were maintained in DMEM with 10% FBS and 1% P/S (HEK growth medium) at 37°C, 95% humidity and 5% carbon dioxide.

Production of the hMSC-TERT4/iCas9-4 cell line
The Lenti-iCas9-Neo plasmid (Addgene #85400, gift from Qin Yan33) was used to transform competent Escherichia coli DH5α and purified in parallel with psPAX2 (Addgene #12260) and pMD2.G (Addgene #12259). Integrity was evaluated by linearization (using appropriate restriction enzymes for each plasmid) followed by gel electrophoresis.

Lenti-iCas9-neo virus was produced in HEK293T cells and used to transduce hMSC-TERT4. Medium was changed to growth medium with 0.5 μg/mL doxycycline 34 hours after transduction to induce GFP and Cas9 expression. After 24 hours, GFP-positive cells were selected by fluorescence-activated cell sorting (FACS) using FACS Aria III (BD Biosciences) and sorted into and maintained in 37°C growth medium. Data analysis was performed with BD FACSDiva software. The initial GFP-positive population was called hMSC-TERT4/iCas9. Single-cell clones of this population were generated and the line hMSC-TERT4/iCas9-4 was selected for further studies.

Quantification of incorporated virus copy number
Genomic DNA was used for qPCR-based determination of lentiviral transfer sequence copy number in hMSC-TERT4/iCas9-4 cells. Cells were harvested in cold lysis buffer (0.1% SDS, 1% Triton X-100, 0.15 M NaCl, 1 mM EDTA, 20 mM Tris pH 8) and stored at -80 °C prior to sonication. Entire cell lysates were sonicated in 130 μL microTUBE AFA fiber Pre-Slit Snap-Cap 6x16 mm tubes (Covaris) on an ME220 Focused-ultrasonicator (Covaris) for 20 s, peak power 75, duty% factor 26.67, cycles/burst 200. Shreded DNA was purified and used as a template for each qPCR reaction with forward and reverse primer targeting the long terminal repeat (LTR) of the transfer sequence and a genomic region (LTR and NoGene2 primer, respectively, Supplementary Information Table 2) and FastStart Essential DNA Green Master (Roche) on a Lightcycler 480 II (Roche). Incorporated virus copy number was calculated as:
2−Ct LTR
Copy# = 2−Ct NoGene2
Since both primers amplify two target regions per viral transfer sequence and per genome, respectively.

Genome editing in hMSC-TERT4/iCas9-4
To design spacer sequences for guide RNAs (gRNAs), enhancer sequences were extracted from the UCSC Genome Browser51 (genome.ucsc.edu, hg19 assembly) based on the location and width of peaks from DNase-sequencing of hMSC-TERT4 cells during osteogenic and adipogenic differentiation13. Spacer sequences were found using the online tool sgRNA Designer (portals.broadinstitute.org/gpp/public/analysis-tools/sgrna-design)52,53, selected by on-target and off-target rank (higher rank is predicted to be more active and specific) and ordered as default Alt-R CRISPR-Cas9 crRNAs (Integrated DNA Technologies). For spacer sequences see Supplementary Table 4. tracrRNA for forming tracrRNA:crRNA duplexes, i.e. gRNAs, was also ordered from Integrated DNA Technologies.

hMSC-TERT4/iCas9-4 cells were induced with 36 μg doxycycline in MilliQ water per 25 mL medium at 80-90% confluence. After 6 hours, cells were reverse transfected in a 12-well plate with a total of 12 pmol crRNA:tracrRNA duplexes (i.e. gRNA), 2 μL DharmaFECT 1 (Thermo Scientific) and 25,000cells/cm2 per well. Cells were kept in transfection medium (1/3 Opti-MEM, 2/3 α-MEM with 10 % FBS) for 2 days at 37 °C, 95% humidity, and 5% carbon dioxide, and seeded as single cells in 96-well plates (0.5 cells/100 μL medium). Clonal enhancer-deletion cell lines were then maintained in growth medium. Two controls were included: ‘No gRNA’, which received the same treatment as all other cell lines, except the transfection complexes contained no gRNA; ‘Ctr gRNA’, which received a pair of gRNAs mediating a deletion in intron 1 of PPP1R12C (adeno-associated virus “safe harbor” site)34.

For screening of enhancer-deleted clones, genomic DNA was harvested in DNeasy Blood and Tissue Kit for DNA isolation (Qiagen), and DNA was extracted according to the manufacturer’s instructions. Target regions in the deleted regions were amplified by end-point PCR using Phusion High-Fidelity DNA polymerase (NEB) (primers in Supplementary Table 3), and the presence of PCR product was evaluated on a Fragment Analyzer (Advanced Analytical Technologies). For each enhancer deletion, three independent clones lacking the wild-type PCR product were used to assess the size of the deleted region by PCR analysis of genomic region, followed by purification of PCR products with QIAquick PCR Purification Kit (Qiagen), and Sanger sequencing (primers in Supplementary Table 3, deletion size in Supplementary Table 6, and Supplementary Fig. 9).

Oil-Red-O and Alizarin Red S staining
Cells were fixed in 4% formaldehyde in PBS for 15 min. For Oil-red-O staining of lipids, cells were incubated with 0.3% Oil-red-O in 60% isopropanol (filtered) at room temperature in the dark for 30 min. For Alizarin Red S staining of calcified matrix, cells were incubated with 1% Alizarin Red S in MilliQ water, pH adjusted to 4-4.2 with HCl and NH4OH. Cells were washed five times with water and subjected to macroscopic and microscopic photography.

Western blotting
Cell lysates were harvested from hMSC-TERT4 after 0 or 7 days of adipogenic differentiation, or hMSC-TERT4/iCas9-4 cells treated with 0.5 μg/mL doxycycline for 6 hours, in lysis buffer (50 mM Tris-HCl pH 6.8, 10% glycerol, 2.5% SDS, 10 mM β-glycerophosphate, 10 mM NaF, 0.1 mM sodium orthovanadate, 1 mM PMSF (added fresh), 1x Complete (Roche, added fresh), boiled and treated with benzonase. Protein concentration was determined using the Pierce BCA Protein Assay Kit (Thermo Scientific) in accordance with the manufacturer’s instructions. Dithioeythritol (DTE) was then added to 0.01 M. For Western blotting, 10-40 μg protein sample was used per well in a 10% gel (separation: 10% acrylamide/bis-acrylamide solution, 0.4 M Tris-HCl pH 8.8, 0.1% SDS, 0.1% ammonium persulfate in MilliQ (APS), 0.001% tetramethylethylenediamine (TEMED); stacking: 6% acrylamide/bis-acrylamide solution, 0.125 M Tris-HCl pH 6.8, 0.1% SDS, 0.1% APS, 0.001% TEMED). For specific protein detection, the wet-blotted poly-vinyl-difluoride membrane was incubated for 1 hour at room temperature or overnight at 4°C with primary antibody: rabbit-anti-Lamin A (1:10,000, L1293, Sigma-Aldrich), mouse-anti-FLAG (1:500, F1804-1MG, Sigma-Aldrich), rabbit-anti-PPARG (1:1000, 2443S, Cell Signaling) in PBS with 1.5% BSA, and then 1 hour at room temperature with secondary antibody: goat-anti-rabbit (1:1,000, P0448, DAKO) or goat-anti-mouse (1:2,000, P0447, DAKO) in PBS with 1.5% BSA. Enhanced chemiluminescence detection was performed using Immobilon Classico Western HRP substrate (Merck millipore), and either exposing the membrane to an X-ray film, or imaging using the Amersham Imager 680.

Quantification of lipid droplet area
Lipid accumulation was imaged using Nikon Eclipse Ts2-FL with differential interference contrast and 40x total magnification. Lipid droplet area was quantified using ImageJ54, by thresholding on the background subtracted greyscale image (using identical thresholds for all images) and quantifying the area in pixels that reached the threshold.

Analysis of gene expression using qPCR
hMSC-TERT4/iCas9-4 enhancer deletion and control cell lines (three clones for each condition) were harvested using Isol-RNA Lysis Reagent (5PRIME) after 10 days of adipogenic differentiation. RNA was purified using phenol-chloroform extraction and EconoSpin columns (Epoch Life Sciences). For cDNA synthesis, up to 1000 ng RNA (always equal amounts in replicates) was treated with DNase I and reverse-transcribed using Moloney murine leukemia virus reverse transcriptase (Invitrogen). cDNA was used for qPCR with FastStart Essential DNA Green Master (Roche) on a Lightcycler 480 II (Roche) (primers Supplementary Table 2). Data for each sample were normalized to the level of transcript of TATA-box protein II B (TBP).

Gene expression analysis using RNA-seq, and enhancer activity analysis using DNase-seq and MED1 ChIP-seq
Our previously published RNA-seq data from hMSC-TERT4 at day 0, 4 hours, day 1, 3, 7, and 14 of adipogenic and osteogenic differentiation (GEO; GSE113253) were analyzed as described13. Previously published DNase-seq and MED1 ChIP-seq data from hMSC-TERT4 at day 0, and day 1, 3, 7, and 14 of adipogenic and day 0, and day 1, 3, and 7 of osteogenic differentiation (GEO; GSE113253) were analyzed as described13, and visualized using using the pyGenomeTracks software55.

Analysis of putative enhancer-anchored interactions using ECHi-C
Our previously published ECHi-C data from hMSC-TERT4 at day 0 and day 10 of adipogenic differentiation (GEO; GSE140782) was analyzed as described15. Briefly, raw sequencing reads were mapped to human genome (hg19) and quality controlled using HiCUP56. CHiCAGO57 pipeline was used for generating interaction counts and calculating p-values for interactions. Interactions anchored in putative enhancers were visualized using the WashU Epigenome Browser58.

Analysis of chromatin conformation in PPARG locus in different stem cell states using public Hi-C data
Hi-C data set for hMSC-TERT4 at progenitor state (GEO; GSE140782)15 were quality controlled using HiCUP56, mapped to human genome (hg19), filtered and processed using runHiC59 using default parameters. HiC matrixes for ESC, mesodermal cells, MSC, neuronal progenitor cells (GEO; GSE52457) were directly accessed and visualized using Genome Browser from YUE lab (http://3dgenome.fsm.northwestern.edu/) at 40 kb resolution.

Analysis of super-enhancers in hMSC-TERT4
To identify and rank super-enhancers, we used Rank Ordering of Super-Enhancers (ROSE)31,32 algorithm, using MED1 ChIP-seq files at day 14 of adipogenesis in hMSC-TERT4, and using default parameters: stitching distance 12,5 kb and TSS exclusion 0 bp.

ChIP
ChIP was performed as described13. Briefly, only for MED1 ChIP cells were crosslinked in 2 mM disuccinimidyl glutarate (DSG) in PBS (100 µL of 0.5 M DSG (Proteochem) in DMSO per 25 ml of PBS) for 20 minutes. For H3K4me1, CTCF, C/EBPβ, C/EBPα and PPARγ ChIP, this step was omitted. Then, cells were crosslinked in 1% formaldehyde in PBS for 10 minutes followed by quenching with 0.125 M glycine for 10 minutes. Cells were harvested in cold lysis buffer (0.1% SDS, 1% Triton X-100, 0.15 M NaCl, 1 mM EDTA, 20 mM Tris pH 8) and stored at -80°C prior to sonication. Samples were sonicated using ME220 Focused-ultrasonicator (Covaris) for 15 min, peak power 75, duty% factor 5, cycles/burst 1000 (for MED1 ChIP-qPCR, μL microTUBE AFA fiber Pre-Slit Snap-Cap 6x16 mm tubes (Covaris)), and for 10 min, peak power 75, duty% factor 26.66, cycles/burst 500 (for ChIP samples followed by sequencing, miliTUBE 1 ml AFA Fiber(100) (Covaris)). Antibodies directed against MED1 (A3E00-793A, Bethyl Laboratories, 1 μl/IP), H3K4me1 (ab8895, abcam, 3 μl/IP), CTCF (61311, Active Motif, 3 μl/IP), C/EBPβ (sc-7962, Santa Cruz, 1 μl/IP), C/EBPα (8178S, Cell Signaling, 3 μl/IP) and PPARγ (2443S, Cell Signaling, 5μl/IP) were used for immunoprecipitation.

For ChIP-qPCR precipitated DNA was used as a template for qPCR with FastStart Essential DNA Green Master (Roche) on a Lightcycler 480 II (Roche) (primers Supplementary Table 2). MED1 occupancy was determined as % recovery of input sample, and the signal in each sample was normalized to signal at the promoter of the house-keeping gene ZRANB3 to account for ChIP efficiency (see Supplementary Fig. 6 for MED1 % recovery in different genomic regions where signal was expected to be similar across samples).

For ChIP-seq, precipitated DNA was prepared for sequencing following the manufacturer's protocol (Illumina). ChIP-seq was performed on minimum two independent biological experiments. ChIP libraries were sequenced on the Illumina NovaSeq 6000, and the quality of sequenced reads was estimated using FastQC (https://qubeshub.org/resources/fastqc). Sequencing reads were mapped to the human genome (hg19) using STAR60. Picard (http://broadinstitute.github.io/picard/) was used to remove duplicates from aligned reads. Peaks were identified using MACS261, and reproducible peaks were identified using BEDTools Intersect62. Bigwig files were created using bamCoverage with RPKM normalization and averaged between replicates using bigwigAverage using deepTools63, and then visualized using pyGenomeTracks software55. For MED1 ChIP-seq at day 3 of adipogenesis in control ‘No gRNA’ and in E+2 deletion and in E+102 deletion in hMSC-TERT4-iCas9/4 , and for C/EBPβ ChIP-seq at day 1 of adipogenesis in control ‘No gRNA’ and in E+2 deletion and in E+102 deletion in hMSC-TERT4-iCas9/4, size factors were estimate using DESeqDataSetFromMatrix and estimateSizeFactors functions in DESeq2 v.1.42.064. Then size factors were used to create normalized bigwig files with BamCoverage from deepTools63, and bigwig replicate files were averaged with bigwigAverage from deepTools63. For MED1 ChIP-seq at day 3 of adipogenesis, and C/EBPβ ChIP-seq at day 1 of adipogenesis in control ‘No gRNA’ and in E+2 and E+102 deletion in hMSC-TERT4-iCas9/4, differential binding was performed using DiffBind v.3.12.0 with edgeR (unadjusted P-value ≤ 0.05)65. Plot profiles, and heatmaps were created using DeepTools63, with commands computeMatrix, plotHeatmap or plotProfile using normalized and averaged bigwig files.

Capture-C
hMSC-TERT4/iCas9-4 enhancer deletion and control cell lines (one clone for each condition) were harvested at day 0 and 10 of adipogenic differentiation. Capture-C was performed using a customized version of a previously published protocol66. Briefly, cells were cross-linked for 10 minutes in 2% formaldehyde in PBS, quenched with 0.125 M glycine, harvested in fresh cold lysis buffer (10 mM Tris-HCl pH 8, 10 mM NaCl, 0.2% Igepal CA-630 (Sigma), 1x Complete (Roche)), incubated 20 minutes on ice and snap-frozen in liquid nitrogen before storage at -80°C. 3C libraries were generated by in-nucleus double digestion with 3x 150 U of each DdeI and DraI restriction enzymes (New England Biolabs) added over approximately 24 hours (final buffer conditions after last addition of enzyme: 1x enzyme buffer, 0.06% sodium dodecyl sulfate, 1.5% Triton X-100). Restriction enzymes were chosen based on the distribution of restriction sites in the genomic region surrounding the PPARG locus, aiming for the regulatory elements of interest located on restriction fragments (RFs) of approximately 120-250 bp. Digestion overhangs were filled in with dNTP and blunt-end ligated. De-crosslinked chromatin was purified. Approximately 5 μg of purified 3C library was sheared to an average size of 450 bp by sonication on an ME220 Focused-ultrasonicator (Covaris) and cleaned up with 0.85:1 AMPure XP beads (Beckman Coulter) sample ratio to remove fragments <200 bp. For sequencing library preparation, sheared DNA was end-repaired, adenine-tailed and ligated to NEBNext hairpin adaptors (New England Biolabs), followed by U excision by USER Enzyme (New England Biolabs). Libraries were amplified using Q5 high-fidelity DNA polymerase (New England Biolabs), universal primer and 8-bp-index primers (New England Biolabs) in 7 cycles of PCR.

For capture, 5’ biotinylated 116-120 bp DNA oligos were designed to target the middle of single RFs spanning, or located as close as possible to, the regulatory elements of interest, avoiding repetitive regions (Sigma standard DNA oligos in tubes (single-stranded), Supplementary Table 5). From each sample, 1650 ng of amplified 3C library was used to generate one multiplexed library for capture (biological replicates were kept separate). Per sample in the multiplexed library, 13 fmol in total of pooled biotinylated DNA oligos targeting all 10 regulatory elements of interest was used for capture, i.e. 1.3 fmol of each oligo, and 130 fmol total oligo pool if the multiplexed library consisted of 10 samples. Sequencing adapters and indexes were blocked using a total of 2 nmol custom blockers (Integrated DNA Technologies, PAGE purified DNA oligos) per sample in the multiplexed library, and two rounds of capture were performed as previously described66, using 10 cycles of amplification after each capture. According to the manufacturer's instructions, multiplexed and double-captured libraries were subjected to 150 bp paired-end sequencing (Illumina).

The resulting raw sequencing reads were mapped to hg19 and quality-filtered using the CCseqBasicS pipeline67 in a version customized to handle reads produced with the restriction enzymes DraI and DdeI (Supplementary Fig. 3). Resulting unique interaction counts from chromosome 3 were normalized to 50,000 counts per sample and viewpoint. As expected, in samples where a captured viewpoint was deleted, the raw interaction count was generally very low compared to samples where the viewpoint was intact (Supplementary Fig. 4a). However, for the viewpoints RE-130, E-123 and E+89, the capture oligo and deleted sequence only displayed partial overlap, resulting in some capture of the deleted region. To simplify downstream analysis, counts in all samples, where a captured viewpoint was deleted, were set to 0 for that viewpoint. The analysis was then restricted to the region chr3: 12.1-12.6 Mb (region of interest, ROI) containing the PPARG locus, and interactions to regions ±2 kb from the viewpoint were excluded from further analysis. Normalized interaction counts were averaged over the RFs overlapping bins of 1 kb or 5 kb, giving a count for each 1 kb or 5 kb bin along the ROI seen from each viewpoint in each sample. The expected count for each viewpoint in each sample (at either RF or binned resolution), as a function of distance from the viewpoint, was then predicted by fitting a cubic polynomial model to the observed sample count per 1 kb bin (for element specific interactions, see below) or 5 kb bin (for ROI-wide interactions), given by the geometric mean count per bin across all intact viewpoints in that sample. P-values describing the chance of an observed interaction count at a given distance from a viewpoint randomly being higher than the predicted count for that distance were calculated using Student’s t-test and adjusted for multiple testing using the Benjamini-Hochberg method. Significant interactions (FDR ≤ 0.01, 208 ROI-wide unique interactions out of 1000 possible across all viewpoints) identified in any condition were kept for further analysis. Differential interactions were identified using Student’s t-test but not adjusted for multiple testing to avoid inflating the number of false negatives in this relatively small dataset. Since interaction counts in the ‘No gRNA’ and ‘Ctr gRNA’ controls correlated well, compared to the remaining conditions (Supplementary Fig. 4b), these conditions were considered as four replicates of one control for all statistical tests (instead of two replicates of two controls). To specifically analyze interactions between the selected PPARG community enhancers, element-specific observed and predicted interaction counts were calculated based on un-binned counts as the mean interaction count seen from one viewpoint (enhancer) to all non-captured RFs overlapping another enhancer. This resulted in two interaction counts for each element-specific interaction, one count seen from each viewpoint, which were summed to obtain a single total observed or predicted count for that interaction. Significant interactions (FDR ≤ 0.01, 34 element-specific interactions across all viewpoints out of 45 possible) were identified and subjected to differential analysis as described above. Element-specific interactions were visualized using the circlize package for R68.

Analysis of genetic variants within PPARG locus and their predicted effect size using Enfomer
To identify genetic variants within the PPARG TAD (chr3:12174868–12526858, hg19) associated with cardiometabolic traits, we integrated evidence from genome-wide association studies (GWAS) for common variants and from whole-exome/whole-genome sequencing (WES/WGS) for rare variants; (a) the Common Metabolic Diseases Knowledge Portal69, (b) the CARDIoGRAMplusC4D (Coronary ARtery Disease Genome wide Replication and Meta-analysis - CARDIoGRAM), plus The Coronary Artery Disease (C4D) Genetics) [The data on coronary artery disease contributed by the CARDIoGRAMplusC4D and UK Biobank CardioMetabolic Consortium CHD working group who used the UK Biobank Resource (application number 9922). Data have been downloaded from www.CARDIOGRAMPLUSC4D.ORG, (c) the DIAGRAM consortium [https://diagram-consortium.org/]70-73 and (d) the variant lists from the GWAS Catalog associated with BMI, T2D, RT2DadjBMI, WHR, and WHRadjBMI. For the Common variants (GWAS), the significant variants (p-value < 5 x 10-8) were extracted and submitted to the HaploReg server74 for identifying their haplotypes in LD > 0.8. This LD expansion captures the full haplotype structure of known PPARG association signals. As expected from LD architecture and power constraints, individual p-values of proxies span several orders of magnitude; such variation does not indicate lack of biological relevance but rather the statistical resolution of the locus. For the rare and low-frequency variants (WES/WGS), due to low allele counts, such variants cannot reach genome-wide significance, but they represent bona fide PPARG-locus associations identified in large sequencing studies. Rare variants were not LD-expanded and are reported individually. Across these sources, we assembled a combined set of 3,596 variants (common + rare) within the PPARG locus.

To predict variant effects on the PPARG2 and our enhancer regions we used the Enformer model39, and we evaluated variant effects on the genomic coordinates of the enhancer regions and the two transcripts of PPARG2 (ENST00000287820 [chr3:12351472-12434356] and ENST00000683699 [chr3:12351500-12434162]). In this analysis, we considered the effects of alternate alleles vs. the reference allele of the variants and their haplotypes to evaluate the effects on PPARG2 enhancer activity (using H3K27ac as proxy) and its expression (CAGE) in adipocyte-related tissues (Supplementary Table 1). The maximum absolute effects of the queried variants in these signals were used to represent variant effects on PPARG2 enhancer and CAGE activities. To define a background effect and significance level on the enhancer and CAGE activities, we selected 4000 variants at random from the GWAS catalogue list of significant hits (p-value < 5 x 10-8). Using Enformer, we calculated the effect of these randomly selected variants on the region of their closest transcripts (the transcript regions were denied using GENCODE v43). The absolute maximum effects of the variants on the same list of adipocyte-related signals were used to define a background effect. The 90 percentiles of the effects from these randomly selected variants were then used as the cutoffs for identifying significant effects of our of variants [Enhancer (90%); negative cutoff: -0.0090, positive cutoff: 0.0088, CAGE (90%): negative cutoff: -0.1147, positive cutoff: 0.1585]. The variants with maximum Enformer effects are reported in the (Supplementary Data 1, Supplementary Data 2).

To further annotate the variants and their corresponding haplotypes, we queried their associated PheWAS from the Common Metabolic Diseases Knowledge Portal API - HuGeAMP69, and aggregated traits related to T2D with an association p < 0.05 as a nominal threshold for secondary trait annotation only. We note that this threshold was not used for variant inclusion, which was determined exclusively by the GWAS and WES/WGS criteria described above. These phenotypes or associated traits and their corresponding p-values are reported in the (Supplementary Data 1, Supplementary Data 2). These tables reports the global and ancestry-specific allelic frequencies of the variants which were directly queried from the NIH Allele Frequency Aggregator (AALFA) [https://www.ncbi.nlm.nih.gov/snp/docs/gsr/alfa/].

We further investigated these variations using FinnGen data75 to identify their significant eQTL, pQTL, and GWAS associations, as well as phenotypic trait fine-mapping PIP. A subset of these findings is presented in Supplementary Datas 1 and 2, with the comprehensive reports accessible in Supplementary Data 3.

Analysis of binding affinity using HT-SELEX datasets
Analysis was performed following the study41. Briefly, effect of the genetic variant; rs181025382 on DNA binding was estimated as log fold change in MinSeqs using two HT-SELEX data sets for PPARγ: without ligand treatment, and for rosiglitazone-treated PPARγ. These datasets were used to determine the enrichment of sequence flanking ± 20 bp from rs181025382 (chr3:12494623-12494663/hg19) for both the reference allele (A/hg19) and the rs181025382-altered allele (G/hg19).

As controls, we selected the top 1000 PPARγ-bound enhancers in hMSC-TERT4 (identified based on a reproducible signal from DNase- and MED1-ChIP seq at day 14 of adipogenesis and PPARγ signal at day 10 of adipogenesis in hMSC-TERT4), and 1000 randomly selected enhancers that were not bound by PPARγ (based on the same signals). MinSeq sequences analysis from two PPARγ datasets mentioned earlier was then performed on 1000 bp sequences these two control enhancers sets.

The log fold change in the enrichment was calculated with log ( E(alt.allele)+ η / E(Ref.allele)+ η ), where E(alt.allele) and E(Ref.allele) are enrichment values calculated by scanning MinSeq sequences within the queried reference and altered sequences. η is a parameter to facilitate unintended issues caused by division by small numbers.

BPNet Model analysis
A BPNet43 model was trained on PPARγ ChIP at day 10 of adipogenesis in hMSC-TERT4 data. Replicates were merged using samtools v. 1.1376. 5’ ends of reads were extracted using bedtools v.2.30.062 using subcommand genomecov for each strand, and base pair resolution bigwig files were constructed using bedGraphToBigWig. Low threshold peaks were called with macs2 v.2.2.9.1 using subcommand callpeak using the parameters -p 0.1 --gsize 2701495761 --keep-dup all --shift -75 --extsize 150. Peak regions were blacklist filtered77 with bedtools intersect using the ENCODE hg19 blacklist. The model was trained using the BPNet-lite v.0.8.1 implementation of BPNet (https://github.com/jmschrei/bpnet-lite). We used an input sequence length of 2114 bp and an output profile width of 1000 bp. We used 64 convolutional filters and an 8 layers deep tower of dilated convolutions. Augmentation was used on the training sequences, using a max jitter (random lengthwise translation) of 128 bp and reverse complementing with a probability of 0.5. Chromosomes 8 and 10 were used for validation, while the remaining autosomal chromosomes plus chromosome X were used for training. The model with the lowest loss on the validation chromosomes after training for 50 epochs was selected.

Motif discovery and attributions
To identify sequence motifs driving the ChIP-seq signal we utilized TF-MoDISco78, running the updated implementation tfmodisco-lite v.2.2.0 (https://github.com/jmschrei/tfmodisco-lite). The identified motifs were matched against the JASPAR motif database79 using tomtom80 from meme v.5.5.5.

SNP analysis
To analyze the effect of individual SNP at or near the PPARG locus we used tangermeme v.0.2.1 (https://github.com/jmschrei/tangermeme) and performed in-silico saturation mutagenesis (ISM) using the count predictions.

Genome-wide PPARG motif mutation analysis
We extracted the PPARG:RXR dimer motif (ID: MA0065.2) from the JASPAR database. We then ran FIMO using tangermeme v.0.2.1 (https://github.com/jmschrei/tangermeme) to identify motif matches across all PPARγ peak region PPARγ ChIP-seq from day 10 of adipogenesis in hMSC-TERT4, using a p-value threshold of 0.01. The matches were then narrowed down to unique and high scoring matches, selecting only those with a score above 10. Finally, we selected those matches with an A at position 9 of the motif match, so that we can induce the A > G mutation as for rs181025382. We then calculated the count prediction for sequences centered on the variant both with and without the A > G mutation (Supplementary Fig. 8b-c).

Statistics and reproducibility
For each enhancer deletion or hMSC-TERT4/iCas9-4 control, mRNA expression and MED1 ChIP-qPCR were assessed based on two independent differentiation experiments of the three independent clones (n=2 independent biological replicates). For ChIP-seq in enhancer deletion and controls two clonal cell lines per enhancer deletion was used in two independent differentiation experiments (n=2 independent biological replicates, n=2 clonal cell line per condition). For Capture-C, only one clonal cell line per enhancer deletion was used in two independent differentiation experiments (n=2 independent biological replicates). For publicly available sequencing data from hMSC-TERT4, data were used from at least two independent differentiation experiments (n=2 independent biological replicates, n=3 for RNA-sequencing). Lipid accumulation in response to enhancer deletion was assayed based on one differentiation experiment using 3 clonal cell lines per enhancer deletion. For statistics not performed within published packages or methods, normality was evaluated visually before testing using the qqnorm() function in R. Equal variances were either not assumed or tested using var.test() in R. Student’s t-test (for data with equal variances) and Welch’s t-test (for data with non-equal variances) were performed in R using the t.test() function in a two-sided manner. P-values were adjusted for multiple testing in R using the p.adjust() function with the Benjamini-Hochberg method where indicated. Polynomial regression was performed in R using the lm() function. All bar plots show individual measurements and/or mean across a group (defined in Figure legends), and error bars represent the standard deviation (SD), standard error of the mean (SEM) or range between two experiments (indicated in Figure legends). Line plots of Capture-C interactions show the mean of two independent experiments with error bars representing a range.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This is a rigorous, well-written, and highly impactful manuscript that provides novel insights about the complex regulatory mechanisms of PPARG2 activation during adipogenesis. The authors identify an extensive, locus-wide \"highly interconnected enhancer (HICE) community\" that drives PPARG2 activation. This community is primed in undifferentiated cells and undergoes dynamic remodeling and increased connectivity during adipogenesis. Through carefully designed inducible CRISPR/Cas9 genome editing experiments, the authors show that a subset of enhancers are essential for efficient PPARG2 induction and adipocyte differentiation. Strikingly, the deletion of the super-enhancer E+102 almost completely prevented PPARG2 activation and was highly correlated with reduction in lipid accumulation, linking enhancer molecular biology to cellular physiology. Overall, this work represents a strong advance in understanding PPARG2 regulation and its relevance to metabolic health. I only have a few minor suggestions that can hopefully further strengthen the work.\n1. Instead of “race” the authors should use “ancestry” as think that’s what they mean.\n2. More information about the selected GWAS SNPs would be helpful to better understand the results. In supp tables 1 and 2, please include nominal p-values for the listed SNPs. The methods state that SNPs with p<0.05 were selected. This is quite permissive. If the SNP p-values do not surpass genome-wide significance, then they are unlikely to be the functional SNP in the GWAS signal.\n3. The enhancer variant to target gene prediction approaches are good, but it would be more convincing if there’s also eQTL data to support the relationships. Is that the case? It’s not critical to show this and there may be missing eQTL data or low-powered studies unable to detect relevant eQTL, but some comment on this would be helpful.\n4. More information should be provided about the SNPs that are the focal point of discussion in the enhancer communities. Based on the GWAS data, are any of the SNPs in the credible set? If so, what is the PIP? Such information would help to establish a stronger argument that the relevant SNPs are good candidates for causal SNPs in the signal.", "reviewer_2": "This manuscript presents a comprehensive and well-executed investigation into the enhancer regulation of the PPARG locus during adipogenesis. By integrating chromatin conformation analyses, CRISPR/Cas9-mediated enhancer deletions, and ChIP-seq experiments, the authors demonstrate that multiple enhancers act in a cooperative yet non-redundant manner to control PPARG2 expression. The study highlights the particular importance of enhancer E+102 in establishing enhancer connectivity and transcriptional output, while also implicating C/EBPβ and PPARγ as key regulators of enhancer activity. Furthermore, by linking these mechanistic insights to GWAS-identified variants associated with metabolic disease, the work provides an important translational perspective. Overall, the findings are novel, convincing, and of broad interest to the fields of adipocyte biology, enhancer function, and transcriptional regulation.\nMinor Comments\n1. Overall, the data convincingly support the idea that E+102 plays an important role in the enhancer network regulating PPARG2. However, as exemplified by the statement on page 10 (lines 244–247), “during adipogenesis, all super-enhancer constituents, E+102 in particular, become more connected to the PPARG2 TSS (Fig. 2b),” the evidence presented in Figure 2b does not appear entirely clear-cut to justify such a strong emphasis. Nevertheless, given that E+102 emerges as a key enhancer from multiple lines of analysis, it would greatly strengthen the manuscript if the authors could discuss possible mechanistic explanations for why E+102, among all super-enhancer constituents, is activated earlier and more robustly. Such a discussion would help to reconcile the somewhat modest support from connectivity data with the overall conclusion of its essential role as a “seed enhancer.”\n2. In Figure 6b, the ChIP-seq data suggest that deletion of E+2 may also reduce C/EBPβ occupancy at several enhancers, even if the effect is less pronounced compared to E+102 deletion. This makes it somewhat unclear on what basis the authors describe certain changes as “significant.” Since the term appears in several parts of the manuscript, it would be helpful if the authors could clarify how “significant” is defined in their analyses—for example, whether based on statistical testing, fold-change thresholds, or other criteria. A brief explanation would make the interpretation of the data more transparent and consistent.\n3. The overlap between enhancers and GWAS-identified variants shown in Figure 7 is highly interesting and adds an important translational dimension to the study. However, it would strengthen the manuscript if the Discussion addressed more explicitly how these variants are thought to act. Specifically, do the authors consider that the variants alter enhancer function itself, for example by affecting chromatin accessibility or enhancer–promoter interactions? Or are they more likely to primarily influence transcription factor binding, such as that of PPARγ or C/EBPβ? Alternatively, is the precise mechanism still unclear? Even a brief speculative discussion of these possibilities would enhance the interpretation and clinical relevance of the findings.\n4. Several enhancer–enhancer interaction labels appear to be missing in Fig. 4a. Please revise to ensure all relevant interactions are clearly annotated.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "Cetnarowska et al. investigate how PPARG2 is activated during adipocyte differentiation using human stromal progenitors (hMSC-TERT4) as model. Using Hi-C, Capture-C, ChIP-seq, and CRISPR/Cas9 enhancer deletions, the authors show that the PPARG locus contains a highly interconnected enhancer community. They demonstrate that multiple enhancers act in a non-redundant and cooperative manner, with promoter-proximal and super-enhancer constituents being essential for PPARG2 induction. Feedback regulation by PPARγ and C/EBPα further amplifies enhancer activity. Finally, they integrate GWAS data and deep learning predictions, showing that key enhancers overlap with non-coding variants linked to metabolic disease.\nOverall, this is an excellent manuscript that provides important mechanistic insight into the regulation of PPARG2, a master regulator of adipogenesis. The study is highly relevant to both basic biology and metabolic disease. The manuscript has only a few limitations. E.g. it relies on immortalized hMSC-TERT4 cells, which may not fully capture in vivo adipose biology. Functional readouts are limited to mRNA expression and lipid accumulation; additional assays such as insulin sensitivity, adipokine secretion, or metabolic flux analyses could further support the conclusion that enhancer deletions affect adipocyte function. The integration of GWAS data was not further experimentally validated. Also, the work focuses on only the PPARG2 locus, thereby limiting its potential broader conclusion.\nDespite these caveats, the manuscript is very strong, methodologically rigorous, and conceptually innovative. I believe it is already highly mature and suitable for publication after only minor revisions.\nSome typos:\nline 75: communitiesis -> communities\nline 220: three independent clones not transfected > three independent clones transfected\nline 574: studes -> studies\nline 575: distance-depednent > distance-dependent\nline 581: essays > assays"}, "author_response": "Response: We thank the reviewer for her/his positive feedback and for recognizing the importance of our studies.\n1. Instead of “race” the authors should use “ancestry” as think that’s what they mean.\nResponse: We thank the reviewer for recognizing this mistake, which we have now corrected in the manuscript.\n2. More information about the selected GWAS SNPs would be helpful to better understand the results. In supp tables 1 and 2, please include nominal p-values for the listed SNPs. The methods state that SNPs with p<0.05 were selected. This is quite permissive. If the SNP p-values do not surpass genome-wide significance, then they are unlikely to be the functional SNP in the GWAS signal.\nResponse: We thank the reviewer for this helpful comment. We agree that the original wording conflated two different analytical steps and did not adequately separate (i) variant inclusion criteria from (ii) secondary trait annotation criteria in the CMD-Knowledge Portal. We have now modified the Methods (pages 33-35) and relevant Results sentences to clarify the following three points:\n(1) Variant inclusion is based on genome-wide significant (GWS) GWAS signals for common variants and large WES/WGS consortia for rare variants. Common variants were obtained by identifying all genome-wide significant cardiometabolic GWAS index SNPs (p < 5×10⁻⁸) within the PPARG TAD across CMD-KP, DIAGRAM, CARDIoGRAMplusC4D, and the GWAS Catalog. These index variants were then expanded to high-LD proxies (r² > 0.8, CEU), as is standard at established loci.\nRare and low-frequency variants were obtained from CARDIoGRAMplusC4D WES/WGS and UK Biobank CardioMetabolic Consortium CHD analyses. As expected for variants with very low allele counts, such variants cannot achieve genome-wide significance due to limited power, but they constitute bona fide PPARG-locus associations from large sequencing resources. Rare variants were not LD-expanded.\n(2) We agree with the reviewer that variants whose p-values are far weaker than their LD index SNP have a lower statistical prior for causality. We therefore now report all variant-specific p-values in Supplementary Tables 1 and 2.\n(3) The previously mentioned “p < 0.05” applied only to PheWAS trait annotation, not variant selection. This threshold was used exclusively to identify secondary phenotypes nominally associated with variants already present in PPARG GWAS or sequencing-based loci. It did not influence which variants were included in the analysis. This is now clarified explicitly in the Methods.\nAs requested, Supplementary Tables 1 and 2 now include each variant’s association p-value, enabling readers to evaluate the relative strength of evidence for each variant.\nFinally, we emphasize that the goal of this section is not to declare locus-level associations, these are already well established for PPARG, but rather to integrate:\n(i) Variants across the frequency distribution with prior support for cardiometabolic trait associations from large-scale variant association studies, (ii) Enformer-predicted regulatory effects, (iii) enhancer maps and ChIP-seq, (iv) and, thanks to this reviewer, now fine-mapping and QTL evidence to prioritize candidate functional variants within the PPARG regulatory landscape.\n3. The enhancer variant to target gene prediction approaches are good, but it would be more convincing if there’s also eQTL data to support the relationships. Is that the case? It’s not critical to show this and there may be missing eQTL data or low-powered studies unable to detect relevant eQTL, but some comment on this would be helpful.\nResponse: We appreciate reviewer’s feedback, and now we have evaluated QTLs for the top-scoring Enformer predicted regulatory variants (Supplementary Table 11). The data revealed that we indeed find colocalized eQTLs for eight common variants, including rs7647481 and rs7649970 located in enhancer E-1. We have addressed this important point in more detail in the results section (page 19). As expected, due to limited power we did not detect a colocalized eQTL for the rare variant rs181025382 (MAF=0.0001, CEU) in enhancer E-102.\n4. More information should be provided about the SNPs that are the focal point of discussion in the enhancer communities. Based on the GWAS data, are any of the SNPs in the credible set? If so, what is the PIP? Such information would help to establish a stronger argument that the relevant SNPs are good candidates for causal SNPs in the signal.\nResponse: We appreciate this excellent point and have now queried the FinnGen data-server to retrieve available QTL (p < 5e-8) and fine-mapping information of the variants. Since the posterior inclusion probabilities (PIP) of variants depends on the number of variants in the corresponding LD blocks we used FinnGen’s default threshold of 0.01 to identify fine-mapping associations. In a new Supplementary Table 10, we have included variants fine-mapped to phenotypes and their PIP from the top-scoring Enformer predicted regulatory variants. This information is now also included in the results section (page 18).\nResponse: We thank the reviewer for her/his positive feedback and for recognizing the novelty and broad interest to the field.\nMinor Comments\n1. Overall, the data convincingly support the idea that E+102 plays an important role in the enhancer network regulating PPARG2. However, as exemplified by the statement on page 10 (lines 244–247), “during adipogenesis, all super-enhancer constituents, E+102 in particular, become more connected to the PPARG2 TSS (Fig. 2b),” the evidence presented in Figure 2b does not appear entirely clear-cut to justify such a strong emphasis. Nevertheless, given that E+102 emerges as a key enhancer from multiple lines of analysis, it would greatly strengthen the manuscript if the authors could discuss possible mechanistic explanations for why E+102, among all super-enhancer constituents, is activated earlier and more robustly. Such a discussion would help to reconcile the somewhat modest support from connectivity data with the overall conclusion of its essential role as a “seed enhancer.”\nResponse: We agree with the reviewer that the important role of E+102 is intriguing and cannot be explained by connectivity alone. In the Discussion in the original submission (page 25), we stated “The difference in importance of super-enhancer constituents cannot be explained solely by connectivity, as all constituents interact with the promoter-proximal enhancers and the PPARG2 TSS. However, the importance of these enhancers appears to be closely linked to their effect on the activity of other enhancers in the locus.“ Later in the Discussion (page 24), we suggested that “the obligate role of E+102 may be due to its dual functions, acting as a seed enhancer in the super-enhancer and as a key target enhancer for feedback action by PPARγ (and possibly other adipogenic TFs)”.\nWe have now further elaborated on this discussion to propose possible molecular mechanisms leading to E+102 being the first and most strongly activated enhancer constituent. First, the affinity of the E+102 sequence for early proadipogenic transcription factors is likely to be important. At day 1 of adipogenesis, where the super-enhancer constituents are not yet DNase accessible, the early proadipogenic transcription factor C/EBPβ appears to bind with similar strength to E+102, E+107 and E+111 (Figure 6a), indicating that early C/EBPβ occupancy cannot explain the difference on occupancy. The combination of other early proadipogenic transcription factors binding to E+102 compared with other super-enhancer constituents are likely to be important for the special role of E+102.\n2. In Figure 6b, the ChIP-seq data suggest that deletion of E+2 may also reduce C/EBPβ occupancy at several enhancers, even if the effect is less pronounced compared to E+102 deletion. This makes it somewhat unclear on what basis the authors describe certain changes as “significant.” Since the term appears in several parts of the manuscript, it would be helpful if the authors could clarify how “significant” is defined in their analyses—for example, whether based on statistical testing, fold-change thresholds, or other criteria. A brief explanation would make the interpretation of the data more transparent and consistent.\nResponse: We thank the reviewer for pointing this out. In all ChIP-seq analyses of binding to the 9 enhancers, “significant” refers to statistical significance (EdgeR differential binding analysis (DiffBind)) between binding in cells with deletion of enhancer E+2 or E+102 compared to control cells (p-value <0.05). We did not apply multiple testing corrections, given the small, hypothesis-driven set of tests. Cell lines with deletion of enhancer E+2 show fewer sites with a significant decrease of C/EBPβ occupancy at day 1 compared to cell lines with deletion of enhancer E+102. Similarly, we assessed MED1 differential binding at day 3 of adipogenesis in cell lines with deletion of enhancers: E+2, and E+102 compared to control cells (Fig. 5c). The statistical testing is indicated in the relevant figure legends as well as in the Methods section.\n3. The overlap between enhancers and GWAS-identified variants shown in Figure 7 is highly interesting and adds an important translational dimension to the study. However, it would strengthen the manuscript if the Discussion addressed more explicitly how these variants are thought to act. Specifically, do the authors consider that the variants alter enhancer function itself, for example by affecting chromatin accessibility or enhancer–promoter interactions? Or are they more likely to primarily influence transcription factor binding, such as that of PPARγ or C/EBPβ? Alternatively, is the precise mechanism still unclear? Even a brief speculative discussion of these possibilities would enhance the interpretation and clinical relevance of the findings.\nResponse: We thank the reviewer for raising this point. The listed variants are predicted to affect enhancer accessibility, based on the H3H27ac signal. We have now explained more explicitly (page 18) that we expect causal variants to affect the affinity of transcription factor binding to their target enhancers, thereby affecting enhancer activity as assessed by H3K27ac and affecting PPARG2 promoter activity.\n4. Several enhancer–enhancer interaction labels appear to be missing in Fig. 4a. Please revise to ensure all relevant interactions are clearly annotated.\nResponse: We thank the reviewer for alerting us to Figure 4a, which we agree was difficult to understand. We had only indicated the interactions that are significantly affected by one or more enhancer deletions on day 0 or day 10 of adipogenesis. We have now revised Figure 4a such that all interactions detected between the enhancers are indicated on the figure, the interactions which are significantly affected by enhancer deletions are indicated with bold font, whereas interactions that are unaffected by deletions are indicated in grey color. Moreover, we replaced crosses with ‘NA’, to indicate interactions that cannot be assessed as they involve the deleted enhancers.\nResponse: We thank the reviewer for her/his positive feedback, and for noting its importance for both basic biology and metabolic disease.\nThe manuscript has only a few limitations. E.g. it relies on immortalized hMSC-TERT4 cells, which may not fully capture in vivo adipose biology. Functional readouts are limited to mRNA expression and lipid accumulation; additional assays such as insulin sensitivity, adipokine secretion, or metabolic flux analyses could further support the conclusion that enhancer deletions affect adipocyte function. The integration of GWAS data was not further experimentally validated. Also, the work focuses on only the PPARG2 locus, thereby limiting its potential broader conclusion.\nResponse: We appreciate the reviewer’s thoughtful comments and agree that the study has a few limitations.\nModel system. Our work relies on adipogenic differentiation of immortalized hMSC-TERT4 cells, which does not fully recapitulate in vivo adipose biology. We now state this explicitly in the Discussion.\nFunctional readouts. The aim of this study was to investigate transcriptional regulation at the PPARG locus, not the downstream effects of limiting PPARG expression. We have focused on mapping chromatin accessibility and interactions, transcription factor occupancy, PPARG2 expression, and lipid accumulation. Assays such as insulin sensitivity, adipokine secretion, or metabolic flux could be valuable follow-ups in model systems (ideally in vivo) that reach a more mature adipocyte state than is the case for hMSC-TERT4. To assess the downstream effects of modulating PPARG expression, we have measured expression of a few adipogenic marker genes: CEBPB (which is upstream of PPARγ), FABP4, CEBPA (which are PPARγ target genes) at day 10 of differentiation in control and enhancer-deleted hMSC-TERT4-iCas9-4 lines. These data are now provided as a new Supplementary Fig. 5g.\nGWAS integration/allelic effects. We acknowledge this limitation. However, an experimental validation of GWAS associations for variants located within studied enhancers and predicted to affect PPARG2 expression is beyond the scope of this study. Our primary aim was to focus on the mechanistic workings of the enhancer community in the PPARG locus and to characterize the contribution of individual enhancers to PPARG2 activation and adipogenesis.\nLocus focus. In our hMSC-TERT4 model system, differentiation leads to induced expression of only the PPARG2 isoform (Supplementary Figure 1b). We acknowledge that this limits generalization beyond PPARG2; however, as stated on page 4, it also provides a clean framework to investigate cooperativity within the enhancer network regulating the adipocyte-specific PPARG2 promoter without interference from the PPARG1 promoter.\nDespite these caveats, the manuscript is very strong, methodologically rigorous, and conceptually innovative. I believe it is already highly mature and suitable for publication after only minor revisions.\nResponse: We thank the reviewer for her/his very positive feedback.\nSome typos:\nline 75: communitiesis -> communities\nline 220: three independent clones not transfected > three independent clones transfected\nline 574: studes -> studies\nline 575: distance-depednent > distance-dependent\nline 581: essays > assays\nResponse: We thank the reviewer for noticing these typos, which we have now corrected."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed all comment. Congratulations on this important study!", "reviewer_2": "I consider that the authors have adequately addressed all points I previously raised. I have no further essential comments.", "reviewer_3": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.", "reviewer_4": "The authors made considerable efforts to revise the manuscript and have satisfyingly addressed the reviewers' comments. In my view, the quality of the manuscript has further improved and is of excellent quality, sufficient for publication."}, "author_response": "Response: We thank the reviewer for her/his positive feedback and for recognizing the importance of our studies.\nResponse: We thank the reviewer for her/his positive feedback.\nResponse: We thank the reviewer for co-revieweing the manuscript, and the Nature Communications initiative that supports peer-review training and recognition for Early Career Researchers.\nResponse: We thank the reviewer for her/his positive feedback and for recognizing the importance of our studies."}]

Decision letter