An Atlas of Promoter Chromatin Modifications and HiChIP Regulatory Interactions in Human Subcutaneous Adipose-Derived Stem Cells

The genome of human adipose-derived stem cells (ADSCs) from abdominal and gluteofemoral adipose tissue depots are maintained in depot-specific stable epigenetic conformations that influence cell-autonomous gene expression patterns and drive unique depot-specific functions. The traditional approach to explore tissue-specific transcriptional regulation has been to correlate differential gene expression to the nearest-neighbor linear-distance regulatory region defined by associated chromatin features including open chromatin status, histone modifications, and DNA methylation. This has provided important information; nonetheless, the approach is limited because of the known organization of eukaryotic chromatin into a topologically constrained three-dimensional network. This network positions distal regulatory elements in spatial proximity with gene promoters which are not predictable based on linear genomic distance. In this work, we capture long-range chromatin interactions using HiChIP to identify remote genomic regions that influence the differential regulation of depot-specific genes in ADSCs isolated from different adipose depots. By integrating these data with RNA-seq results and histone modifications identified by ChIP-seq, we uncovered distal regulatory elements that influence depot-specific gene expression in ADSCs. Interestingly, a subset of the HiChIP-defined chromatin loops also provide previously unknown connections between waist-to-hip ratio GWAS variants with genes that are known to significantly influence ADSC differentiation and adipocyte function.


Introduction
The distribution of adipose tissue throughout the body plays a significant role in predicting the health status of overweight and obese people independent of body mass index (BMI) [1,2].Excess accumulation of fat in the upper body (apple-shaped) is positively correlated with higher HbA1c, circulating triglycerides (TG), and adverse serum lipid profiles [3].In contrast, excess accumulation of fat in the lower body (pear-shaped), such as in the gluteofemoral (GF) depot, is negatively correlated with the same metabolic disease markers [4].
A recent theory to explain the effect of differential fat accumulation on metabolic health posits that the lower body adipose tissue serves as a sink for "healthy" lipid deposition and limits fat accumulation in the upper body, notably in visceral adipose tissue, the latter of which is associated with chronic inflammation and insulin resistance [5,6].Although it is well established that growth hormone, cortisol, and sex steroids influence fat distribution [7][8][9], the underlying mechanism for why people develop the apple vs. pear body shape is still not completely understood.At the cellular level, besides the obvious role of mature adipocytes to sequester lipids, there is a likely contribution from the precursor cells or adipose-derived stem cells (ADSCs) which have the capacity to differentiate into mature adipocytes to store excess lipid.In fact, previous studies have demonstrated that human primary subcutaneous abdominal (ABD) vs. GF-ADSCs, have differential adipogenic capacity in vitro [10].In addition, there are distinct transcriptional signatures, chromatin marks, and DNA methylation patterns in ABD vs. GF adipose tissue [11][12][13] that are partially maintained in isolated ADSCs following culture in vitro [14,15].Taken together, these observations provide evidence for an underlying epigenetic memory that contributes to the different patterns of gene-expression-driven phenotypes.Our working hypothesis is that these cell-autonomous epigenetic programs maintain unique ABD-and GF-AT characteristics in the cultured ADSCs and contribute to unique ABD and GF adipose tissue characteristics.These earlier studies integrated gene expression and individual epigenetic marks to help explain the sustained patterns of differential gene expression in ADSCs that were isolated from different subcutaneous adipose depots and cultured over several rounds of cell division.In the current work, we aim to extend these preliminary results and interrogate the epigenomic landscape of ABD and GF-ADSCs around the TSS of the differentially expressed genes between the two adipose tissue depots by using an extensive ChIP-seq analysis for multiple histone marks (H3K4me3, H3K4me2, H3K27me3, H3K9me3), CTCF, and RNA polymerase II along with ATAC-seq to probe chromatin openness.
In the earlier studies mentioned above, gene annotation of the regulatory regions was performed using a nearest-neighbor linear-distance approach which only showed a modest connection between differential chromatin features and differential gene expression [14,15].We recognized this was an overly simplistic approach because chromatin is highly organized and condensed with DNA packaged in a highly ordered fashion with histones and other proteins into a complex three-dimensional network [16,17].The resulting highly condensed chromatin serves to position distally located regulatory elements close to proximal gene promoters that would otherwise be located far away from each other based on a linear (2D) view of the genome.High-resolution three-dimensional methods including Hi-C and ChiAPET were developed that capture these long-range interactions after partial digestion of the DNA followed by ligation of closely juxtaposed ends that are brought into close proximity by looping out of the intervening DNA [18][19][20].
In the second part of the current work, we aimed to determine how the differential gene expression patterns in ABD and GF-ADSCs were significantly influenced by chromatin modifications and long-range chromatin interactions using the related HiChIP method which, when focused on H3K27ac, will identify loops that are anchored through an active transcriptional enhancer region [21].We next integrated the HiChIP data set with our differentially expressed RNA-seq data set that compares gene expression patterns in ABD vs. GF-derived ADSCs.These data were then compiled into an atlas that combines the different chromatin marks and active enhancer connectome related to gene expression patterns for ABD vs. GF-ADSCs isolated from apple and pear-shaped women.
This in-depth analysis of the chromatin structure and organization of ABD and GF-ADSCs provides both an initial in-depth understanding of the intrinsic genomic regulatory features that influence the functionally distinct cellular phenotypes of ABD and GF subcutaneous adipose tissue depots, and a resource for the adipose tissue research community.We also demonstrate the value of this data set as a resource by cross-referencing this information with a data set of WHR-related SNPs; together, these investigations provide significant new information for how different adipose depots contribute to differential adipose patterning and metabolic disease risk in humans.

Results
To evaluate the cell-autonomous differences and explore the molecular regulation of gene expression between ABD and GF adipose tissue depots, we isolated adipose-derived stem cells (ADSCs) from paired ABD and GF adipose tissue originating from five apple and five pear-shaped women.Principal component analysis using all the clinical parameters (the most relevant are listed in Supplementary Table S1), showed that individuals within each group are highly similar and that the two groups are well separated from each other (Figure 1A).For this reason, we analyzed the apple and pear samples separately.The ADSCs were all cultured the same way and passaged the same number of times (±1) (see Section 4 for details) prior to harvest.The overall workflow of the study is described in Figure 1B.In summary, we performed ChIP-seq for histone marks, CTCF, and RNAPII along with ATAC-seq on freshly isolated chromatin.An additional aliquot of cells was frozen and used for RNA-seq analysis.We also performed H3K27Ac-enriched HiChIP to evaluate the 3D organization of the active enhancer connectome.

Results
To evaluate the cell-autonomous differences and explore the molecular regulation of gene expression between ABD and GF adipose tissue depots, we isolated adipose-derived stem cells (ADSCs) from paired ABD and GF adipose tissue originating from five apple and five pear-shaped women.Principal component analysis using all the clinical parameters (the most relevant are listed in Supplementary Table S1), showed that individuals within each group are highly similar and that the two groups are well separated from each other (Figure 1A).For this reason, we analyzed the apple and pear samples separately.The ADSCs were all cultured the same way and passaged the same number of times (±1) (see Section 4 for details) prior to harvest.The overall workflow of the study is described in Figure 1B.In summary, we performed ChIP-seq for histone marks, CTCF, and RNAPII along with ATAC-seq on freshly isolated chromatin.An additional aliquot of cells was frozen and used for RNA-seq analysis.We also performed H3K27Ac-enriched HiChIP to evaluate the 3D organization of the active enhancer connectome.From each biopsy, the stroma vascular fraction was isolated and the ADSCs were cultured in media supplemented with serum and growth factors.Cells were harvested, their chromatin was isolated and used for RNA isolation, ChIP, ATAC, and HiChIP assays, followed by sequencing.

Transcriptomic Signatures of ABD and GF-Adipose-Derived Stem Cells
We first identified the differentially expressed genes between ABD and GF-ADSCs using RNA-seq.Setting a cut-off of a 1.7-fold change and an FDR of 0.1, the RNA-seq analysis showed a total of 599 differentially expressed genes (DEGs) between the ABD and GF-ADSC samples, of which 364 exhibited GF-enriched features and 285 were ABDenriched.This analysis revealed six clusters of DEGs (Figure 2A), stratified based on their level of expression in apple-and pear-shaped subjects.Subcutaneous adipose tissue biopsies were performed on five apple-shaped and five pear-shaped subjects.From each biopsy, the stroma vascular fraction was isolated and the ADSCs were cultured in media supplemented with serum and growth factors.Cells were harvested, their chromatin was isolated and used for RNA isolation, ChIP, ATAC, and HiChIP assays, followed by sequencing.

Transcriptomic Signatures of ABD and GF-Adipose-Derived Stem Cells
We first identified the differentially expressed genes between ABD and GF-ADSCs using RNA-seq.Setting a cut-off of a 1.7-fold change and an FDR of 0.1, the RNA-seq analysis showed a total of 599 differentially expressed genes (DEGs) between the ABD and GF-ADSC samples, of which 364 exhibited GF-enriched features and 285 were ABDenriched.This analysis revealed six clusters of DEGs (Figure 2A), stratified based on their level of expression in apple-and pear-shaped subjects.
Among the forty-two GF-enriched genes highly expressed in GF pear samples (Figure 2A, cluster 1), we found six pro-adipogenic marker genes (GPX3, PRDM1, FABP3, KLF5, WNT5B, and NRG1 [22]) which would be consistent with GF-ADSCs in pear-shaped women having the capacity to differentiate more robustly compared to GF-ADSCs from apple-shaped women.Importantly, further analysis revealed that the most significantly enriched pathway in cluster 1 contained genes involved in Wnt-β-catenin signaling (Figure 2B), which is known to play a significant role in adipocyte differentiation [23].More recently, an extensive analysis combining scRNA-seq combined with a xenograft mouse model as validation, showed that Wnt signaling preserves progenitor cell multipotency during adipose tissue development [24], which would be predicted to ensure a healthy pool of progenitor cells capable of differentiation to mature adipocytes.Cluster 2 includes 214 GF-enriched genes highly expressed in apple-shaped samples (Figure 2A), several of which are important for lipid droplet formation and others that increase in expression during adipogenesis, such as ALPL, CD44, CD36, CFD, MME, ENPP2, ELOVL2, and ABI3BP.Cluster 2 also contains fibroblastic or fibrotic marker genes (TGFBR3, LAMA3, TNFSF9, S100A4, VCAM1, CXCL12, ANPEP, COL5A) not found in cluster 1.These genes might participate in the formation of collagen, which can form a scaffold that constrains adipocyte expansion due to mechanical stress in GF apple samples [25,26].
The 285 ABD-enriched genes were divided into three clusters of roughly equal size (clusters 4, 5, and 6-Figure 2A) based on their body shape expression pattern.Cluster 4 includes genes highly expressed in ABD apple samples but with low expression in GF pear samples.This cluster includes genes activated by hypoxia (TES, STC2, DDIT4, SLC2A1), cytokine/chemokine genes (IL33, IL11, CXCL5), and PDGFA, which is known to stimulate adipose progenitor proliferation and self-renewal but also is associated with increased adipose tissue fibrosis [32,33].Cluster 4 also contains two other genes that may influence adipose tissue expansion: BAMBI, a gene known to regulate reactive oxygen levels [34] and PAWR, a suppressor of p53 [35].
Finally, cluster 6 includes genes highly expressed in ABD pear samples but with low expression in the other three samples, especially GF apple (Figure 2A).Interestingly, the pathway with the highest enrichment score in cluster 6 includes 80 genes related to epithelial-mesenchymal transition (EMT) (Figure 2B).Further analysis revealed these genes are also potential inhibitors of adipogenesis (INHBA, ITGA2, IL32, OXTR, CXCL8, MMP1, TPM1, TFPI2, MYLK, VCAN, COL7A1, PMEPA1, TGM2) [41] and this correlates well with our overall developmental/physiological hypothesis that there is reduced expansion of the ABD depot in pear subjects.
The MCODE analysis [42] revealed protein-protein interactions between the DEGs (Supplementary Figure S1).This analysis identified several proteins interacting with each other.Notably, FMN1 (up in GF depot) and EDN1 (up in ABD depot) were found at the center of two hubs, related to cell growth and inflammation, respectively.Several collagens were also interacting to a high degree.These proteins could be master regulators of the DEGs between ABD and GF-ADSCs and potential therapeutic targets to favorize lipid accumulation in the lower body adipose tissue rather than in the upper body depots.

Selective Epigenetic Hallmarks of Depot-Selective Transcription
The fact that several of the GF and ABD-ADSC-enriched genes were also enriched in the same samples from whole adipose tissue indicates the differential expression pattern is a stable feature that is likely maintained at least in part through depot-selective epigenetic regulation [14,15].To evaluate genomic signatures that might contribute to the differential gene expression patterns highlighted above, we compared the chromatin features present in ABD vs. GF-ADSCs in the same samples used for the RNA-seq analysis.We performed ATAC-seq to evaluate chromatin openness and ChIP-seq targeting different histone marks (active chromatin: H3K27Ac-active enhancers/promoters, H3K4me2-active enhancers/promoters, H3K4me3-mainly active promoters; repressed chromatin: H3K27me3-facultative heterochromatin, H3K9me3-constitutive heterochro-matin), CTCF (genome architectural protein), and the elongating form of RNAPII (phospho Ser 2 in CTD).Then, we integrated the chromatin modification data from all the samples to generate a combinatorial set of emission states using ChromHMM [43].Emission parameters were learned de novo based on genome-wide recurrent combinations of the chromatin marks studies (see above) in ADSCs (ABD and GF combined).Importantly, each emission state was defined by a specific combination of chromatin features that may be associated with distinct biological expression patterns of their linked genes.
We distinguished 10 broad classes of chromatin emission states that are labeled according to their combined predicted influence on gene activity.These include "Genomic enhancers", "Flanking Active TSS", "Active TSS", "CTCF-high", "Bivalent/poised TSS", "Repressed", "Quiescent", "Heterochromatin", "Strong Transcription", "Bivalent Enhancer" (Figure 2C), and they are independent of the tissue depot or body shape chromatin source.As a first step in understanding how these combined features influence gene activity, we displayed the average read density scores for the ten unique chromatin states around the transcription start site (TSS) of the DEGs in the two different depots (ABD and GF) and from the two different body shapes (apple and pear) in Figure 2D.Only three emission states were enriched at the TSS of differentially expressed genes.The first one "active TSS" (state 3, orange on the graphs in Figure 2D), contains all activate chromatin-associated histone marks along with CTCF, RNAPII, and the open chromatin signature (ATAC-seq peak) combined with very low levels of the repressive marks H3K27me3 and H3K9me3 (see legend Figure 2D).The active TSS state was enriched at the TSS of the GF-enriched genes belonging to clusters 1 and 3 in the GF samples (dark green arrows in Figure 2D) and at the TSS of the ABD-enriched genes belonging to cluster 5 in the ABD samples (dark orange arrows in Figure 2D).
The second enriched emission state corresponded to repressed regions (state 6, gray on the graphs in Figure 2D) which are characterized by an enrichment of the repressive H3K27me3 mark in the absence of the other features.The genes in clusters 5 and 6 are marked by a high level of the H3K27me3 repressive mark (grey arrows) at the TSS in the GF samples.This suggests that these genes are more highly expressed in ABD-ADSCs because their expression is repressed in the GF region.Taken together, these results support the concept that differential combinations of active and repressive chromatin marks at DEG's TSSs contribute to depot-specific gene expression patterns in ABD and GF-ADSCs.We also found 'bivalent domains' of histone modifications (i.e., harboring both the repressive mark H3K27me3 and the activation-associated marks) near the TSS of genes with depot-specific expression (blue lane, Figure 2D).
To obtain a better assessment of the enhancer regions in the ABD vs. GF samples, we plotted only the "Genic enhancer" state (state 1 in Figure 2C) around the promoter of the DEGs. Figure 2E displays the average read density scores for this specific state around the TSS of the DEGs for the four different groups (ABD and GF, apple and pear subjects).The genes in clusters 2 and 3 (GF-enriched genes) are marked by an increase in enhancer marks in the GF samples whereas the genes in cluster 6 (ABD-enriched genes) are marked by an increase in enhancer marks in the ABD samples.These observations suggest an active role of enhancer genomic regions in depot-specific gene regulation in human ADSCs.

Alteration of Active and Repressive Epigenetic Marks Associated with Depot-Selective Gene Expression
To evaluate the individual contributions of histone modifications in depot-selective gene expression, we separately analyzed the ChIP-seq marks within the TSS of the DEGs and compared the data between the ABD and GF-ADSCs.We focused on the individual histone marks that define the active TSS state (H3K4me3, H3K27Ac, and H3K4me2) and repressed state (H3K27me3), and calculated the differences in their respective ChIP-seq signals in ABD and GF-ADSCs around the TSS (±2 kbp) of the ABD and GF-DEGs.Those reaching statistical significance (p < 0.05) are colored (orange for ABD-enriched histone mark and green for GF-enriched histone mark) in the volcano plots of Figure 3 (apple subjects) and Supplementary Figure S2 (pear subjects).The full list of DEGs with the fold change for each histone mark is listed in Supplementary Table S2 for the apple subjects and Supplementary Table S3 for the pear subjects.For example, HOXC11 and TBX15 (GFenriched genes, cluster 3) showed an enrichment of the active histone marks and a decrease of the repressive mark H3K27me3 in the ABD samples compared to the GF samples.This provides a more detailed evaluation of the association of positive and negative histone modifications with differential gene expression patterns in ADSCs and extends what has been described in other cancer cell model systems [44,45].
Our data also suggest that histone modifications affect expression of genes involved in adipogenesis such as PRDM1, ALPL, RUNX1T1, SFRP2, and GPX3.Importantly, the newly identified GF-enriched genes in our study, ZIC1, GREM2, and IL20RA, were also associated with GF-enriched differential patterns of histone modifications at their TSSs.Among the ABD-enriched genes within the highly active TSS emission state in ABD chromatin, we detected the key developmental genes HOXA5, HOXD1, HOXD3, HOXD8, and TBX5 in cluster 5 (ABD-enriched genes in both body shape types).The differential chromatin marks were also associated with DEGs that are involved in adipose tissue expansion, such as BAMBI, PAWR, and IL8 (clusters 4 and 6; ABD-enriched genes) along with other genes known to limit adipogenesis such as INHBA (cluster 6; ABD-enriched genes highly expressed in pear samples), TIMP1, and CDKN2B (cluster 5).Interestingly, two inflammatory genes (CXCL5 and IL33) showed higher H3K27me3 levels around their TSSs in GF samples compared to their TSSs in ABD samples (cluster 4; ABD-enriched genes highly expressed in apple samples), suggesting that the expression of these two genes may be selectively repressed by H3K27me3 in GF-ADSCs.

HiChIP Regulatory Interactions in ABD and GF-ADSCs
In an earlier study comparing open chromatin regions identified by ATAC-seq with differentially expressed genes in freshly isolated adipocytes, we showed that only a small fraction of the body-shape-specific open chromatin regions were annotated to DEGs [46].In this earlier study, we used linear distance as a guide suggesting that long-range genomic interactions mediated by chromatin looping are likely involved in the differential gene expression patterns.To determine how the differential gene expression patterns in ADSCs from different adipose depots may be influenced by long-range chromatin interactions, especially by enhancer genomic regions as suggested by our data in Figure 2E, we performed H3K27ac-targeted HiChIP on chromatin from ADSCs across the ten subjects and two adipose tissue depots.
We identified 52,489 and 52,615 loops in the apple and pear samples, respectively (hichipper, FDR < 0.01).Each sample had similar levels of high-quality uniquely mapped read pairs (Supplementary Figure S3A).Principal components analysis also showed that samples from each group clustered together, and their patterns were separated based on body shape and depot source (Figure 4A).Supplementary Figure S3B shows that the overall A/B compartment score distribution across all groups was identical for chromosome 7 and this was also evident on the whole genome level as shown by the saddle plots in Supplementary Figure S3C.The median loop length was 41 kb and as expected, the number of interactions decreased with increasing distance between loop anchors (Supplementary Figure S3D).We then overlapped the anchors of the HiChIP loops with gene promoters and enhancers and the resulting loop sets were binned into three different categories: enhancer-promoter loops (18,958 in apples and 19,011 in pears), enhancer-enhancer loops (28,075 in apples and 28,138 in pears), and promoter-promoter loops (5456 in apples and 5466 in pears).
To identify differential loops between the ABD and GF-ADSC samples, we ran diffloop analysis separately for the apple and pear groups.This revealed 852 ABD-enriched loops and 493 GF-enriched loops in the apple samples with a p < 0.05 and fold change of 1.75 between the groups (orange and green symbols in Figure 4B).The number of depotenriched loops was lower for the pear groups (304 ABD-enriched and 238 GF-enriched, Supplementary Figure S3E).To validate the depot-specific regulatory loops identified by HiChIP, we overlapped the loop anchors with the read density from the independently performed H3K27ac ChIP-seq analysis on the same set of samples from Figures 2 and 3.The fold change between the ABD and GF HiChIP reads at the loop anchors highly correlated with the fold change between the ABD and GF H3K27ac signal detected by ChIP-seq at the same loop anchors.This correlation is consistent with the HiChIP pipeline used in our study, accurately identifying authentic depot-enriched loops.
This comparison resulted in a list of high-confidence H3K27ac loops in the ABD and GF-ADSC samples, including promoter and enhancer interactions that we analyzed further below.At 2.5 kb resolution, the H3K27ac HiChIP maps revealed depot-specific promoter-enhancer interactions at the promoter of HOXA genes in the ABD sample which are known ABD-enriched genes [27] (Figure 4D, left) and at the promoter of the TBX15 gene in the GF samples which is a known GF-enriched gene (Figure 4D, right).Additionally, TBX15 expression correlates with WHR and there is some evidence that it is a master transcriptional regulator in adipose tissue [47].The genomic regions that were enriched in loops (HiChIP results) also were highly enriched for the CTCF ChIP-seq peaks that were identified in the CTCF ChIP-seq data used in the ChromHMM analysis in Figure 2C.These sites co-mapped with genomic regions with a high insulation score, consistent with CTCFassociated looping organizing the 3D genomic architecture to regulate gene expression.There was also strong enrichment for H3K27ac binding along with higher RNAPII and other marks associated with gene activation in ABD chromatin at the HOXA locus (Figure 4D, left IGV snapshots and ChromHMM).Similarly, the GF-enriched TBX15 HiChIP loops were associated with more robust peaks for active histone marks and RNAPII in the GF sample (Figure 4D, right IGV snapshots and ChromHMM).

Loop Anchors Harbor DEGs and SNPs That Are Associated with Waist-Hip Ratio in Humans
To further define regulatory regions that might influence differential gene expression between ABD and GF-ADSCs, we first identified the HiChIP loop anchors that were linked to DEGs.In apples, this revealed 323 loops in the genomic regions of the DEGs between the ABD and GF samples, and these loops mapped to 64 unique DEGs.We found approximately the same results (325 loops mapping to 64 unique genes) in the pear group.From these lists we extracted the loops with at least one anchor found at the promoter region of the DEGs which corresponds to thirty-five unique DEGs (Table 1, in apples).The transcription of these genes is likely regulated by the enhancer region we identified by HiChIP (opposite loop anchor in Table 1), and some are potentially involved in fat distribution heterogeneity (HOXA, BDNF, IL33, EPHX2, IGF2BP1).In the last decade, large population studies have used GWAS to explore the genetic influences on WHR [48].Historically, human GWAS studies have been performed without regard to chromatin structure.We next asked if the location of genes important to clinical phenotypes like WHR, or genes involved in adipocyte function, might overlap with our atlas of chromatin structure in ADSCs.We identified known SNPs associated with WHR (48 studies reporting 4797 unique SNPs annotated to genes listed in Supplementary Table S4) located inside chromatin loops: all loops, differential loops in apples (described in Figure 4B), and differential loops in pears (described in Supplementary Figure S3E).This is a conservative analysis as SNPs that were nearby, but not exactly inside the loop anchors, were not included.We found 417 WHR-associated SNPs in loop anchors identified in our study (in apple and in pear samples).
To ascertain if this could be a random finding, we performed a random permutation test and found that the number of SNPs that overlap with loop anchors was significantly higher than the number expected by random (Figure 5A).
From these 417 SNPs, 39 were found in ABD-enriched loops and 7 were found in GFenriched loops (Figure 5B and Table 2).Some of the genes annotated to the depot-enriched loops, in other words having a WHR-SNP in their anchor, were also differentially expressed in one of the adipose tissue depots as depicted in the plots in Figure 5C that emphasize the differences in expression from each individual in both groups (HOXA3, MLXIP, SBF2, PPL, KCNJ6, HOXA11).Interestingly, MLXIP is a bHLH transcription factor that dimerizes with CHREB to regulate the expression of genes involved in glucose metabolism, glycogen synthesis, triglyceride synthesis, and insulin signaling [49].MLXIP expression has also been implicated in the development of metabolic diseases such as obesity, insulin resistance, and T2DM [50,51].Additionally, MLXIP has also been recently identified as a marker of a sub-population of human adipocytes that are highly responsive to insulin [52].

SNP_ID
Gene Name SNP_ID Gene Name SNP_ID Gene Name

SNP_ID
Gene Name SNP_ID Gene Name SNP_ID Gene Name

SNP_ID
Gene Name SNP_ID Gene Name SNP_ID Gene Name

SNP_ID
Gene Name SNP_ID Gene Name SNP_ID Gene Name  In the last decade, large population studies have used GWAS to explore the genetic influences on WHR [48].Historically, human GWAS studies have been performed without regard to chromatin structure.We next asked if the location of genes important to clinical phenotypes like WHR, or genes involved in adipocyte function, might overlap with our atlas of chromatin structure in ADSCs.We identified known SNPs associated with WHR (48 studies reporting 4797 unique SNPs annotated to genes listed in Supplementary Table S4) located inside chromatin loops: all loops, differential loops in apples (described in Figure 4B), and differential loops in pears (described in Supplementary Figure S3E).This is a conservative analysis as SNPs that were nearby, but not exactly inside the loop anchors, were not included.We found 417 WHR-associated SNPs in loop anchors identified in our study (in apple and in pear samples).
To ascertain if this could be a random finding, we performed a random permutation test and found that the number of SNPs that overlap with loop anchors was significantly higher than the number expected by random (Figure 5A).The overlap between WHR-SNPs and loop anchors identified as enriched in the ABD samples was more revealing than the seven WHR-SNPs found in the anchors of the GFenriched loop library (Figure 5B).Indeed, nine SNPs were found in loop anchors in the HOXA cluster on chromosome 7 (Figure 5D) and three SNPs were found in loop anchors in the PPARG gene (Figure 5E).The HOXA genes are differentially regulated in ABD vs. GF adipose tissue, preadipocytes, and adipocytes [15,27,46], whereas PPAR gamma is a master transcription factor enriched in preadipocytes and adipocytes, necessary for adipogenesis and also regulates fat and glucose metabolism [53].
Taken together, these results suggest that these SNPs may affect WHR by the regulation of ABD but not GF adipose tissue function and that this effect is driven by differential looping in the ADSCs and potentially in mature adipocytes.

Discussion
Chromatin loops can link enhancers physically close to their target genes and help to better understand the alterations of gene transcription that affect disease.Our work, described here for the first time in primary human ADSCs, provides an extensive atlas of 3D-associated regulatory interactions.To gain insight into the potential function of these long-range chromatin interactions, we integrated the HiChIP interactome with the genes differentially expressed between the ABD and GF samples, with genes known to influence adiposity and cardiometabolic traits, and with GWAS-SNPs that are associated with WHR.We also established a list of loops that describe differential 3D genomic interactions in two groups of women (apple and pear-shaped).Importantly, some of these interactions were associated with ABD and/or GF-ADSC gene expression profiles that we highlighted by RNAPII ChIP-seq analysis performed in parallel.
In our earlier study where we were limited to using linear annotation, we showed that only a small fraction of body-shape-specific open chromatin regions were annotated to DEGs [46].We proposed that long-range genomic interactions mediated by chromatin looping were likely involved in the differential gene expression patterns.Thus, in the present study, we used H3K27ac HiChIP to interrogate active enhancer-associated looping in regulating depot-enriched gene expression and this resulted in the identification of 35 unique DEGs with associated loop anchors (Table 1).The transcription of these genes is potentially regulated through the enhancer loop interaction revealed in our HiChIP data set (opposite loop anchor in Table 1) and some likely influence differential fat distribution (HOXA, BDNF, IL33, EPHX2, IGF2BP1).Importantly, our work discovered a potential new key transcription factor such as ZFP36L1, an inhibitor of adipogenesis [54,55], as a master regulator of depot-specific gene expression.We described 14 loops at its promoter (Table 1), reflecting its high potential of interaction with other distally located genomic regions.Other genes related to obesity and/or adipogenesis were identified by our HiChIP analysis as regulators of ABD vs. GF gene transcription, such as METTL15 [56] and RBM4 [57].Overall, our study supports the idea that long-range chromatin loops may affect the development or differentiation of ADSCs and could explain in part subcutaneous adipose tissue dysfunction in diseases such as T2D or PCOS.
Previously reported large population studies have relied on GWAS to link genes to WHR [48].Importantly, the gene connections have been performed relying largely on linear annotation and have not typically considered the importance of longer range chromatin interactions that are defined using more involved chromatin looping methods.Using our HiChIP data set, we connected SNPs known to be associated with WHR (48 studies reporting 4797 unique SNPs annotated to genes) with key chromatin loops: all loops (Figure 4B and Supplementary Figure S3E).It should be noted that this is a conservative estimate because we narrowly defined the SNPs to be located within the loop anchors and did not consider closely associated anchors in this analysis.Importantly, this revealed genes that were also differentially expressed in one of the adipose tissue depots (Figure 5C) that are known to influence adipose tissue function including HOXA3, MLXIP, SBF2, and PPL.Taken together, these findings demonstrate that genomic interactions play an important role in adipose depot-specific gene regulation in human ADSCs.In addition, by comparing the loops identified between the two adipose tissue depots studied (ABD vs. GF), we highlighted depot-enriched chromatin interactions that likely contribute to depot-selective 3D chromatin organization; this organization influences gene transcription and therefore the distinct functional phenotypes in ABD vs. GF-ADSCs.
We also report here for the first time in human primary ADSCs, that differential histone modifications at gene promoters influence patterns of depot-selective gene expression in ABD vs. GF depots.By studying the correlation between histone marks and differential gene expression between ABD and GF-ADSCs, our work revealed that combinations of histone marks are associated with transcriptional activity in ABD and GF-ADSCs.When the individual marks were combined to generate a combinatorial set of ChromHMM emission patterns, the data are even more supportive of the model.
However, we cannot formally conclude whether differential gene expression is the cause or consequence of differential histone modifications.Henikoff et al. showed that histone modifications were more likely the consequences than the causes of transcription, especially for H3K4me3 [58].Regardless of the direction, these histone marks provide a stable memory of recent transcriptional activity and provide a template for a robust mechanism to sustain the observed differential pattern of transcription between depots.
A limitation of our work is that we used H3K27Ac HiChIP to identify the depotspecific connectome.However, a depot-enriched loop might be identified as specific due to the fact that those regions exhibit a high depot-enriched H3K27ac signal.We cannot conclude if it is this the result of an actual architectural change or simply a difference in H3K27ac at these anchors.
We focused here on loops associated with genes that were differentially expressed in ADSCs across different depots in apple vs. pear-shaped women.It should be noted that all other key genes involved in adipose function were not differentially expressed in our study.Taken together, these and prior experiments in human ADSCs reveal a potential epigenomic mechanism by which the differential growth and function of adipose tissue depots lead to common metabolic diseases.

Participants, Tissue Collection, and Isolation of Human Adipose-Derived Stem Cells
The method of recruitment, clinical, and biochemical parameters of subjects were previously published by Divoux A. et al. [46].All procedures were performed under a research protocol approved by the AdventHealth Institutional Review Board.A subgroup of 10 healthy premenopausal, weight-stable women were used for this study.Five women displayed lower body adiposity characterized by a waist-to-hip ratio (WHR) < 0.78 (pear group; age = 34 ± 9.6 years; BMI = 29.2 ± 2.26 kg/m 2 ) and five women displayed upper body adiposity, characterized by a WHR > 0.85 (apple group; age = 38 ± 8.1 years; BMI = 28.6 ± 3.54 kg/m 2 ).Briefly, paired abdominal and gluteofemoral subcutaneous white adipose tissue samples were obtained from each participant and the stromal-vascular fractions (SVFs) were isolated by 45 min collagenase digestion (collagenase type I, Wor-thington).SVFs were plated and grown in proliferation medium containing 2.5% FBS, FGF, and EGF.Human adipose-derived stem cell (ADSC) populations were enriched as previously described [14].The cells presenting at their surface the endothelial marker CD31 (MAB2148-C, MilliporeSigma, Burlington, MA, USA) were removed by magnetic beads.

HiChIP Assay
Approximately 5 × 10 6 cells were crosslinked in 1% formaldehyde (methanol-free, dissolved in phosphate-buffered saline-PBS) for 10 min at room temperature in a 10 mL final volume.Formaldehyde was quenched with the addition of 1.5 mL 1M glycine for 5 min at room temperature.Cells were scraped and lysed in lysis buffer (1% Triton x−100, 0.1% SDS, 150 mM NaCl, 1 mM EDTA, and 20 mM Tris, pH 8.0) for 1 h in rotation in a cold room.Isolated nuclei were pelleted by centrifugation, resuspended in 100 µL 0.5% SDS, and incubated for 10 min at 65 • C. SDS was quenched by the addition of Triton-X for 15 min at 37 • C. Nuclei were incubated overnight at 37 • C in a vigorous shaker (speed-850 rpm) in the presence of MboI (375U).The following day, the samples were incubated at 65 • C for 20 min to heat the inactivate MboI.Samples were left at room temperature for 20 min to cool down.Biotin fill-in of sticky ends was performed for 1 h at 37 • C in a vigorous shaker (speed-850 rpm) followed by ligation of blunt ends at room temperature for 6 h while rotating.Nuclei were spun, resuspended in lysis buffer in the presence of 5 µg H3K27ac antibody (ab4729, Abcam, Waltham, MA, USA), and incubated overnight on a rotator at 4 • C. The next day, antibody chromatin complexes were pulled down with protein A paramagnetic beads and sequentially washed: once in wash buffer 1 (1% Triton, 0.1% SDS, 150 mM NaCl, 1 mM EDTA, 20 mM Tris, pH 8.0, and 0.1% NaDOC), twice in wash buffer 2 (1% Triton, 0.1% SDS, 500 mM NaCl, 1 mM EDTA, 20 mM Tris, pH 8.0, and 0.1% NaDOC), once in wash buffer 3 (0.25 M LiCl, 0.5% NP−40, 1 mM EDTA, 20 mM Tris, pH 8.0, 0.5% NaDOC), and once in TE-buffer.After the removal of TE-buffer, DNA was eluted from the beads in an elution buffer.DNA was quantified with the Qubit dsDNA HS kit.Approximately, 40-50 ng DNA was used for biotin pull-down with streptavidin paramagnetic beads.Sequencing libraries were constructed with the Nugen Ovation Ultralow V2 kit (Tecan, Mannedorf, Switzerland) according to the manufacturer's recommendations.Libraries were quantified with the Quibit dsDNA HS kit and subjected to bioanalyzer fragment analysis before paired-end sequencing.

HiChIP Data Processing
HiChIP data were analyzed using the default parameters of nf-core/hic (https:// zenodo.org/records/2669513,accessed on 20 December 2023; version 1.3.0).In summary, the following steps were followed: (1) mapping to the hg19 reference genome using a twostep strategy to rescue reads spanning the ligation sites (bowtie2) [60]; (2) detection for valid interaction products; (3) duplicates removal; and (4) generating raw and normalized contact maps using the ICE algorithm at various resolutions using a cooler.The quality control of the sample was included in the pipeline (HiC-Pro [61]).A/B compartments, saddle plots, and insulation scores were calculated using GENOVA [62].Representative interaction heat maps were generated using cloops2 [63] with the -corr option.We used hichipper [64] to identify chromatin loops, using the consensus H3K27ac ChIP-seq peaks per group.Significant differential looping was calculated using diffloop [65] with -nreplicates set to 3 and -nsamples set to 8.
Loop anchor positions were overlapped with BMI-adjusted waist-hip ratio SNPs (EFO:0007788) and the permutation test was performed to test the significance of overlap compared to permutated (n = 2500) locations using the regioneR package [66].

Sequencing Library Preparation
RNA-seq, ATAC-seq, ChIP-seq, and HiChIP libraries were prepared and sequenced using standard Illumina protocols for a HiSeq 2500 instrument (Illumina, San Diego, CA, USA).

RNA Sequencing and Analysis
RNA sequencing was performed as described by Divoux [15].The raw RNA-seq reads' sequencing quality was evaluated by FastQC and the reads were aligned to the hg19 reference genome using STAR (version 2.7.7a) [67].Genes were quantified using featureCounts from Rsubread (version 2.4.0)[68].The R package, edgeR with paired analysis, was used for differential gene expression analysis with cutoffs CPM > 3 in more than 4 samples.p-value < 0.05 was used to determine statistical significance for differentially expressed genes.DEGs were used for k-means clustering to create modules and visualized as a heat map.

Gene Set Enrichment and Visualization
EnrichR was used for gene set enrichment and visualization [74].The enrichment was calculated to the hallmark gene set of the Molecular Signature Database (MSigDb).Pathways with p-values < 0.05 were selected as significant.

Chromatin State Discovery with ChromHMM
Tissue-specific chromatin states were identified using the ChromHMM (version 1.21) hidden Markov model (HMM) [43].Bam files from RNAPII, CTCF, H3K27ac, H3K27me3, H3K4me2, H3K4me3 ChIP-seq, and ATAC-seq were binarized into default 200 bp bins using the function BinarizeBam from each of the 5 ABD-ADSC and 5 GF-ADSC samples in each group (apple and pear), as previously described [75].We ran ChromHMM with a range of possible states and settled on a 10 states model as it accurately captured information from higher state models and provided sufficient resolution to identify biologically meaningful patterns in a reproducible way [43].
Author Contributions: L.H.: formal analysis, visualization, writing-review and editing.A.D.: conceptualization, formal analysis, investigation, visualization, writing-original draft.K.S.: conceptualization, methodology, investigation, review and editing.E.E.: formal analysis, visualization, review and editing.B.D.: methodology, review and editing.S.R.S.: conceptualization, supervision, funding acquisition, writing-review and editing.T.F.O.: Conceptualization, methodology, super-vision, funding acquisition, writing-review and editing.All authors have read and agreed to the published version of the manuscript.

Figure 1 .
Figure 1.Sample acquisition and study design.(A) Principal component analysis (PCA) plot of the ten subjects used to isolate the ABD and GF-ADSCs based on clinical parameters.The first two principal components (PC) are plotted and colored according to body shape.PCA was performed using all clinical data collected during the clinical study.(B) Overview of the experimental workflow.Subcutaneous adipose tissue biopsies were performed on five apple-shaped and five pear-shaped subjects.From each biopsy, the stroma vascular fraction was isolated and the ADSCs were cultured in media supplemented with serum and growth factors.Cells were harvested, their chromatin was isolated and used for RNA isolation, ChIP, ATAC, and HiChIP assays, followed by sequencing.

Figure 1 .
Figure 1.Sample acquisition and study design.(A) Principal component analysis (PCA) plot of the ten subjects used to isolate the ABD and GF-ADSCs based on clinical parameters.The first two principal components (PC) are plotted and colored according to body shape.PCA was performed using all clinical data collected during the clinical study.(B) Overview of the experimental workflow.Subcutaneous adipose tissue biopsies were performed on five apple-shaped and five pear-shaped subjects.From each biopsy, the stroma vascular fraction was isolated and the ADSCs were cultured in media supplemented with serum and growth factors.Cells were harvested, their chromatin was isolated and used for RNA isolation, ChIP, ATAC, and HiChIP assays, followed by sequencing.

Figure 2 .
Figure 2. Depot-enriched gene expression and chromatin modification analysis of human ADSCs according to body shape.(A) Heat map showing the differentially expressed genes between ABD and GF-ADSCs in apple and pear-shaped subjects based on RNA-seq.The DEGs were grouped into clusters according to their level of expression in apple and pear samples.DESeq2 analysis, FDR < 0.01, FC > 0.75 Genes potentially involved in adipose tissue expansion are cited.(B) Dot plot showing the significant pathways of the DEGs in each cluster.Only the pathways with p < 0.05 are represented.(C) Legend showing the ChromHMM annotated states, with their emission values for individual chromatin marks.(D) Visualization of selected chromatin states (2, 3, 5, 6, 9, 10) around the

Figure 2 .
Figure 2. Depot-enriched gene expression and chromatin modification analysis of human ADSCs according to body shape.(A) Heat map showing the differentially expressed genes between ABD and GF-ADSCs in apple and pear-shaped subjects based on RNA-seq.The DEGs were grouped into clusters according to their level of expression in apple and pear samples.DESeq2 analysis, FDR < 0.01, FC > 0.75 Genes potentially involved in adipose tissue expansion are cited.(B) Dot plot showing the significant pathways of the DEGs in each cluster.Only the pathways with p < 0.05 are represented.(C) Legend showing the ChromHMM annotated states, with their emission values for individual chromatin marks.(D) Visualization of selected chromatin states (2, 3, 5, 6, 9, 10) around the TSS (±5 kbp) of DEGs per depot and body shape groups (rows) within gene clusters (columns).Arrows highlight when the states are visually different between ABD and GF.(E) Visualization of Genic enhancer chromatin state (state#1) around the TSS (±5 kbp) of DEGs per groups (colored) within gene clusters.Arrows highlight when the state is visually different between the ABD and GF samples.

Figure 3 .
Figure 3. Association between depot-enriched expression and depot-enriched chromatin marks at the TSS (±2 kbp) in apple samples.Volcano plots show for each gene and each histone mark studied

Figure 3 .
Figure 3. Association between depot-enriched expression and depot-enriched chromatin marks at the TSS (±2 kbp) in apple samples.Volcano plots show for each gene and each histone mark studied the average fold change of the ChIP-seq signal between ABD and GF-ADSCs at the TSS.Data are represented by cluster of DEGs (rows).Negative fold changes (green) indicate the ChIP-seq signal is significantly enriched in the GF samples, while positive fold changes (orange) indicate the ChIP-seq signal is significantly enriched in the ABD samples.

Figure 4 .
Figure 4. Mapping epigenomic landscapes in ABD and GF-ADSCs.(A) Principal component analysis (PCA) plot of normalized in-loop H3K27ac HiChIP read counts.The first two principal components (PC) are plotted and colored according to body shape.(B) Dot plot showing the correlation of read densities between ABD and GF-ADSCs in apple subjects.Differential loops are colored in

Figure 4 .
Figure 4. Mapping epigenomic landscapes in ABD and GF-ADSCs.(A) Principal component analysis (PCA) plot of normalized in-loop H3K27ac HiChIP read counts.The first two principal components (PC) are plotted and colored according to body shape.(B) Dot plot showing the correlation of read densities between ABD and GF-ADSCs in apple subjects.Differential loops are colored in yellow (ABD-enriched) and green (GF-enriched).The non-significant loops are represented in gray.p-value < 0.05 logFC > 1.75.(C) Density plot showing the correlation between differential looping (x-axis) and differential H3K27ac (y-axis) at loop anchors.The H3K27ac signal was binned into 12 groups based on the magnitude of difference in H3K27ac.Data were plotted for the apple subjects.Similar observations were made for the pear subjects.(D) Genome browser visualization of SKAP2-HOX locus (left) and TBX15 locus (right) in ABD (yellow) and in GF (green) samples.Data were derived from apple subjects.Similar observations were made with data derived from pear subjects.From top to bottom: H3K27ac HiChIP interaction matrices, domainogram of insulation score, CTCF, ATAC-seq, H3K4me3, H3K27ac, H3K4me2, RNAPII, H3K9me3, H3K27me3, ChromHMM states, H3K27ac loops, and gene annotation.Color coding for ChromHMM plots is the same as Figure 2C.

Figure 5 .
Figure 5. Integration of loop anchors and GWAS-SNPs associated with WHR.(A) Permutation test showing the overlap between loop anchors and SNPs.Green lane shows the observed overlap (n = 417) and the gray histogram shows the expected distribution of overlaps by shuffling the SNP positions 2500 times.Dotted line indicates the mean expected overlap, which was used to calculate

Table 1 .
List of DEGs with at least one loop anchor located at their promoter.
ns = not significant; NA = not applicable.

Table 2 .
List of WHR-SNPs overlapping with loop anchors found in both subcutaneous adipose tissue depots.