Methylome and transcriptome maps of human visceral and subcutaneous adipocytes reveal key epigenetic differences at developmental genes

Adipocytes support key metabolic and endocrine functions of adipose tissue. Lipid is stored in two major classes of depots, namely visceral adipose (VA) and subcutaneous adipose (SA) depots. Increased visceral adiposity is associated with adverse health outcomes, whereas the impact of SA tissue is relatively metabolically benign. The precise molecular features associated with the functional differences between the adipose depots are still not well understood. Here, we characterised transcriptomes and methylomes of isolated adipocytes from matched SA and VA tissues of individuals with normal BMI to identify epigenetic differences and their contribution to cell type and depot-specific function. We found that DNA methylomes were notably distinct between different adipocyte depots and were associated with differential gene expression within pathways fundamental to adipocyte function. Most striking differential methylation was found at transcription factor and developmental genes. Our findings highlight the importance of developmental origins in the function of different fat depots.


Differential gene expression between SA and VA
The transcriptomes of SA and VA were generated, using both polyA and whole transcriptome RNA-seq and aligned against a set of 57,818 genes. Differentially-expressed (DE) genes were determined in each dataset, (Supplementary Tables S20 and S21 respectively). Levels of differential expression were highly concordant, R 2 =0.9144 (Supplementary Fig. S1a), so a combined set of DE genes was used for subsequent analysis.
As well as genes discussed in the main text, a number of genes that were highly differentiallyexpressed warrant consideration in relation to the differential biology of SA and VA. Those most differentially expressed in VA include: • cell surface proteins mesothelin (MSLN), uroplakin 3B (UPK3B) and mucin16 (MUC16), • cytokeratin 19 (KRT19) • LDL receptor related protein 2 (LRP2) that binds a variety of ligands, including lipoprotein lipase and has been shown to be involved in leptin signalling (Byun et al. 2014) • PKHD1L1, potentially associated with immune function (Hogan et al. 2003) and • arachidonate lipoxygenase 15 (ALOX15) for which linoleic acid also a substrate (Cole et al. 2012). • ISL1, that has been shown to be a key early driver of adipogenesis in 3T3-L cells ) • Basonuclin (BNC1) a zinc finger transcription not previously associated with a function in adipose tissue.
The most highly differentially-expressed gene in SA is SIM1 (Drosophila single-minded (sim) gene homolog) that has been associated with obesity, but where studies to date have been focused toward neural effects and hyperphagia (Xi et al. 2012). Other genes much more highly expressed in SA are Tubulin-beta-2A (TUBB2A), the homeobox transcription factor IRX5 and the divergentlytranscribed lncRNA CRNDE that is involved in insulin signalling (Ellis et al. 2014), Cyclin D1, CCND1, required for pre-adipocyte proliferation (Marquez et al. 2017), and its neighbour MYEOV.
Notably most of these genes that show highly differential expression between isolated visceral and subcutaneous adipocytes have been shown to be differentially expressed in studies of whole tissue, significant overlap (Fisher pval = 4.41e-07) was found between our list of differentially expressed genes and those identified by Wolfs and colleagues in a study examining whole tissue SAT and VAT from 75 obese adults (Wolfs et al. 2010). Of the 810 genes we found to be up in SA 275 of these were also found to be up in SAT vs VAT, while 333 of the 1395 genes down regulated in SA were also found to be down between SAT and VAT in their data.

DNA methylation analyses
As well as WGBS of the core three subjects, the same DNA samples, as well as samples from two normal weight male subjects (Supplementary Table S1) were analysed using Illumina Infinium Bead Chip Human450 arrays (450K arrays). The MDS plot of 450K array data for the 5 individuals showed very similar clustering as WGBS data (Supplementary Fig. S1b), and for CpG sites on the 450K arrays, WGBS and 450K data were highly correlated, with a Pearson's correlation coefficient >0.94 for all sample comparisons (Supplementary Fig. S1d) and. Comparison of 450K data from our samples with publically available data on SAT and VAT showed a broad clustering by depot type (Supplementary Fig. S1e) 450K data was used to examine variance of methylation levels at each CpG site. Median variance was about 1.5% (Supplementary Fig. S2a). Lowest variance was found promoter and first exon regions, corresponding with regions of very low methylation including CpG islands. Higher variance was seen at 3' UTRs, intergenic and gene body regions (Supplementary Fig. S2b).
Variance at individual CpG sites was well correlated between SA and VA. The minority of CpG sites that show high variance in methylation commonly show such variance in SA, VA and PBL of the same individual, suggesting a genotype influence rather than being tissue-specific effects ( Supplementary Fig. S2c, d).
Supplementary Figure S2. Concordance of DNA methylation between individuals.
(a) Box plots of the within tissue variance of DNA methylation for SA, VA, VAT and PBL, as measured across all probes on 450k arrays. Y-axis = standard deviation, x-axis = annotation category. (b) Box plots of the variance of DNA methylation of combined tissues across genomic annotations as measured by 450k array. Y-axis = standard deviation, x-axis = annotation category.
(c) Plot of correlation of within tissue variance of DNA methylation as measured by 450k array between SA and VA. (d)

Differentially Methylated Regions
Comparing VA and SA 450K array data for the 5 subjects we identified 2,393 DMRs (1,493 up in SA and 900 down in SA) (Supplementary Table S9). We found significant overlap in differential methylation between our 450K and WGBS data. Using GenomicRanges (Lawrence et al. 2013) to directly determine overlaps of DMRs, 1,728 450K DMRs (72%) overlapped or were contained within 1,910 WGBS DMRs (Supplementary Fig. S3d). A strong overlap of differential methylation is also evident at the single CpG level. Using the normally distributed test statistics for differential methylation, we tested concordance of the differential signal between 450K and WGBS, for CpGs with probes on the 450K array. There was an overall concordance value of r = 0.63 (Supplementary Fig. S3e).
Recently, Macartney-Coxson et al. (Macartney-Coxson et al. 2017) used 450K arrays to compare the DNA methylomes of subcutaneous and visceral adipose tissue of 15 morbidly obese subjects before and after weight reduction following bariatric surgery. They identified 703 DMRs associated with 385 genes in common between before and after surgery SA/VA comparisons. We used GenomicRanges to identify overlaps between our 450K DMRs identified on isolated SA and VA and these 703 DMRs. 438 (62%) of the DMRs were shared with ours. This strong overlap is despite our comparison being of isolated adipocytes of normal weight subjects and theirs of whole adipose tissue from obese subjects. Notably, highly differential methylation in regions of important developmental genes was observed in both datasets.

Association of DMRs with body fat distribution and obesity-related SNPs
IntervalStats (Chikina and Troyanskaya 2012) was used to determine the probability of overlaps or proximity of SA/VA DMRs (67,048) with two sets of SNPs, 97 obesity-related SNPs  and 50 body fat distribution SNPs (Shungin et al. 2015 The data shows no relationship between BMI risk SNPs and SA/VA DMRs. However, a suggestive association is seen between SA/VA DMRs and Body fat Distribution. When the analysis is broken down to the highest ranked SNP/DMR pairs for each query direction (

Ontology of genes associated with DMRs
Supplementary Figure

DMRs associated with genes encoding Transcription Factors
Given the enrichment of terms relating to DNA binding in the GO analysis we further examined the relationship of DMRs with transcription factor genes. There are many instances in which multiple transcription factor genes from the same family have DMRs within their promoter or gene body. These include -HOXs, PITXs, CEBPs, DMRTs, GATAs, IRXs, and TBXs. Examples are shown in Supplementary Figs. 5 and 6. Within most of these families both elevated and reduced methylation in SA compared to VA was observed (Supplementary Table S12). However, it is interesting to note that there are a number of gene families where most members have reduced methylation in SA. These include: the key regulator of adipocyte differentiation, PPARG, and its family members PPARα and PPARδ and binding partner RXRA; the HIF family, master regulators of cellular response to hypoxia (HIF1a, ARNT, EPAS1, ARNT2, HIF3a) along with other related genes encoding basic helix-loop-helix genes ARNTL, ARNTL2 and SIM1; and the SMAD family, modulators of numerous signalling pathways (SMAD1,2,3,4,7,9).
Supplementary Figure S5. Methylation profiles across HOX gene clusters. Plots of DNA methylation profiles using the Integrative Genomics Viewer, IGV (Thorvaldsdottir et al. 2013). Complete plot of HOX clusters. SA (orange) and VA (purple) DNA methylation levels from 0% to 100%. Chromosomal location top, RefSeq diagrammatic representation of genes bottom tracks.

Relationship of DNA methylation and DMRs with gene expression
To broadly characterize the relationship between DNA methylation and gene expression in SA and VA we examined the DNA methylation profile around the TSS of all annotated gene promoters. As has been noted for other tissues (Roadmap Epigenomics et al. 2015) decreasing DNA methylation at the TSS was associated with progressively increased gene expression ( Supplementary Fig. S7a). This association is particularly evident for promoters with medium and high CpG content.
Notably, lowest methylation levels were observed about 200 bp after the TSS. Also notable for the most highly expressed genes was the increased levels of methylation both upstream and downstream of the TSS (for more detail with respect to CpG density and gene expression see Supplementary Material and Supplementary Fig. S7b). Interestingly, promoters with low CpG density show a more heterogeneous relationship with gene expression and methylation relative to promoters with higher CpG density.
Supplementary Figure S7. Relationship of Promoter DNA methylation with gene expression. (a) Genes were grouped based on VA expression level into non-expressing (blue), and those in the low (yellow), medium (green) and high (red) thirds. Loess regression was used to fit a curve representing overall methylation levels in VA 3kb either side of the TSS of these grouped sets of genes. (b) Gene sets were further split by CpG density of promoter region into Low (<70 CpGs across 4kb region), medium (70 to 230 CpGs across 4 kb ) and high (>230 CpG across 4 kb) CpG density groups (Weber et al. 2007). Plots show overlayed curves of VA DNA methylation across individual promoters in each group. Colour-coding by gene expression level as in (a). (c) As for Panel (a), but with promoters grouped by CpG density, low (blue), medium (orange) and high (green).

Features and Functions of DE-DM Genes
The visceral fat depot releases more free fatty acid than subcutaneous fat due to lower sensitivity to insulin to suppress lipolysis (Arner 2005). We found that genes associated with lipid droplet and lipid mobilization such as CIDEC, monoglyceride lipase (MGLL), triglyceride lipase (PNPLA3) and PLIN4 all have higher promoter methylation and reduced expression in VA. Other key genes in metabolic regulation such as CPT1A and PDK4, exhibit both elevated expression and promoter methylation in VA. Numerous developmental genes associated with fat distribution and waist-hip ratio were DE-DM ( Fig. 3d and Supplementary Table S14), these included HOXA11 and TBX15 that had both higher expression and DNA methylation in SA (Shungin et al. 2015). Fat distribution is highly regulated by the presence of sex hormones and it has been shown that increased estrogen (17β-estradiol) in abdominal SAT can potentiate lipolysis (Gavin et al. 2013). Interestingly we found that both estrogen receptors 1 and 2 show higher expression and lower methylation in SA. Differences were also observed for mesenteric estrogen-dependent adipose gene, MEDAG, which has been identified as a novel gene of visceral adiposity (Zhang et al. 2012). Angiogenesis and extra cellular matrix (ECM) remodeling are crucial processes for proper formation of adipose tissue, and vascularisation and angiogenic capacity are known to differ between different fat depots (Gealekman et al. 2011). A number of DE-DM genes in our data set are involved in angiogenesis, hypoxia and ECM remodeling ( Fig. 3d and Supplementary Table S14). One of these genes is hypoxia inducible factor HIF3A, for which elevated methylation in blood and subcutaneous fat was found to be associated with obesity/BMI (Dick et al. 2014).  (14). DNA methylation cliffs (outlined in blue) associated with ALX1, DMRT3, EN1, IRX3, TBX18 and SLFN12L. SA (orange) and VA (purple) DNA methylation levels from 0% to 100%. Chromosomal location top, RefSeq diagrammatic representation of genes bottom tracks. Side bars denote direction of differential expression, orange is higher in SA and purple is higher in VA.

Regulatory elements identified within DNA methylomes of purified SA and VA
Putative regulatory elements (UMRs and LMRs) as well as DMVs were identified from WGBS data separately for each sample (for each cell type from each individual). The numbers of regions identified are shown in Supplementary Fig. S10a. Venn diagrams showing the overlaps of scored UMRs and LMRs between the three individuals are shown in Supplementary Fig. S10b. All individual samples had approximately 20,000 UMRs with a strong overall concordance within each sample type (>87% in common). Between-sample calling of LMRs was more variable than UMRs with about 50% being called in a given sample type for all three individuals and a further ~18-22% being called in 2 of 3 individuals. For comparisons between cell types, regions being scored in at least two individuals are included in Supplementary Tables S15-S20.
The averaged DNA methylation profiles across UMRs, LMRs and DMVs for each of SA, VA and PBL are very similar (Supplementary Fig. S10c).

Clustering of TF binding sites in D-LMRs
To identify potential regulatory regions distinguishing SA and VA, we chose LMRs that were in the lowest 25% of methylation and specific to either SA or VA. We further overlapped these with VA-SA DMRs. ENCODE ChIP-seq binding sites (162 TFs) and subcutaneous adipocyte PPARγ ChIPseq binding sites were then mapped onto these DMRs. Binding sites for individual TFs within SA and VA-specific DMR/LMRs are shown in Supplementary Tables S26 and S27, and unsupervised clustering of TF binding and individual regions in Figure 5c and Supplementary Fig. S14a.
Clusters of regions containing binding sites for common set of TFs are indicated in the figures. Clusters containing similar sets of TFs in SA and VA are indicated by letters A-F. Figure S14. TF binding to D-LMRs.

Supplementary
(a) Unsupervised cluster plot of the presence (yellow) or absence (red) of ChIP-seq peak for 162 TFs (vertical axis) within 862 VA D-LMRs (horizontal axis). TF clusters A-F represent clusters of regions containing binding sites for common set of TFs. See table immediately following for detail. (b) UCSC browser image of LRP5. LMRs, WGBS and RNA-seq profiles shown as orange, purple and green respectively for SA, VA and PBL. Transcription factor ChIP-seq peaks shown as grey bars.
The TFs in each cluster are shown in the table below, with those in common between SA and VA shown in bold. TFs contributing to these clusters are discussed further below. For both SA and VA this cluster includes MYBL2, HNF4A, HNF4G and RXRA (the heterodimer binding partner of PPARγ) and binding is also commonly associated with PPARγ and/or CEBPδ. In SA other TFs in this cluster are NFIC, SP1, CEBPD and TEAD4. In post-natal life, members of the FOXA family control glucose metabolism through the regulation of multiple target genes in the liver, pancreas, and adipose tissue. FOXA proteins function as 'pioneer factors' whose binding to promoters and enhancers enable chromatin access for other tissue-specific transcription factors (Friedman and Kaestner 2006). NFIC has recently been shown to regulate adipocyte differentiation through the canonical wnt pathway (Zhou et al. 2017) and to bind to and activate the LRP5 promoter in mice. Among SA LMR/DMRs, we identified a candidate regulatory region in the first intron of the LRP5 gene (about 1.2 kb downstream of the TSS) that contains ChIP binding sites for NFIC, RXRA, HNF4G and FOXA1 and that could contribute to the higher expression of LRP5 in SA compared with VA (Supplementary Fig. S14b).

Cluster C
This includes the chromatin structural proteins CTCF, RAD21 and SMC3 along with ZNF143 and YY1, indicative of regions that are involved in connecting promoters with distal regulatory elements. ZNF143 provides sequence specificity to secure chromatin interactions at gene promoters (Bailey et al. 2015). CTCF, RAD21 and SMC3 are components of the cohesion complex and orchestrate the mitotic chromatin interaction landscape (Kakui and Uhlmann 2017).

Cluster D
While segregating as a separate cluster in SA, the TFs in this cluster are commonly associated in SA LMR/DMRs with Cluster A TFs. TCF7L2 is a major wnt pathway TF and repressor of adipogenesis, while AP2 (TFAP2A and TFAP2C) represses adiponectin expression and has roles in glucose uptake. In VA many of these TFs also form part of an extension of Cluster A; for example TCF7L2 binding sites are found in 114 VA LMR/DMRs of which 75 also contain EP300 ChIP sites (Supplementary Table S27).

Cluster E
Cluster E also contains TFs important in adipocyte function. IRF4 is important in control of lipid handling (Eguchi et al. 2011) while ATF2 is critically involved in early steps of adipogenesis (Maekawa et al. 2010). EBF1 regulates adipocyte morphology, lipolysis and inflammatory pathways (Griffin et al. 2013;Gao et al. 2014). STAT5A also has a significant function in promoting adipogenesis and regulates expression of a number of genes, including adiponectin (Able et al. 2017).
Cluster F This cluster contains four genes in common between SA and VA LMR/DMRs, CBX3, TRIM28, NR2F2 and TAL1. TRIM28 is an epigenetic regulator for which haploinsufficiency has been implicated in obesity (Dalgaard et al. 2016) and that is also associated with CBX3 as a key regulator of nephrogenesis in mice (Dihazi et al. 2015). NR2F2 (COUPTFII) has been previously described as an essential regulator of adipogenesis (Xu et al. 2008). In addition of Cluster F TFs, D-SA LMRs that contain NR2F2 sites show high frequency of binding sites for a number of sequence-specific TFs including TEAD4, GATA2, STAT5A as well as CEBPβ and PPARγ.

CONTACT FOR REAGENT AND RESOURCE SHARING
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Susan Clark (s.clark@garvan.org.au)

Blood collection
Venous blood was collected on the day of elective surgery from the antecubital vein following an overnight fast from midnight. Samples were collected for glucose, lipids, CRP and circulating leukocytes were isolated for DNA methylation measurements. For the latter, 10 mL blood was collected into K2EDTA vacutainers. Following centrifugation at 1,300 rcf for 10 minutes at room temperature, the top plasma layer was removed and the buffy coat layer collected. Cells were washed three times, with vigorous resuspension, in 10 mL Tris-EDTA buffer (100 mM Tris-HCl, 0.1 mM EDTA), with pellets collected after centrifugation at 10,000 rcf for 10 minutes. The final pellet was resuspended in 500 µL of buffer and stored at -80ºC. Blood biochemistry was performed by Sydpath -the Pathology Service of St Vincent's Hospital Sydney, Australia.

Adipose tissue collection
Visceral adipose tissue (VAT) was collected from the greater omental region; subcutaneous adipose tissue (SAT) from the periumbilical site of surgical incision. Adipose samples were separated into two samples: one was snap frozen for whole tissue analysis; the second was prepared for adipocyte isolation -see below.

Adipocyte isolation
Adipocyte isolation was adapted from (Rodbell 1964;Reynisdottir et al. 1994). In brief, freshly isolated VAT and SAT samples were weighed and minced with surgical scissors. For each gram of minced tissue, 1.25 mg/mL of Collagenase I (Sigma Aldrich, CAT#C6885-1G) was added to 2 mL of HEPES buffer containing 7.9 g/L NaCl, 323 mg/L CaCl2.2H2O, 308 mg/L MgSO4.7H2O, 154.3 mg/L Na2HPO4 (anhydrous), 5 mL/L of 1M Hepes (Sigma Aldrich), 450 mg/L D-glucose, pH 7.2 -7.4, 20 g/L BSA. The tissue was incubated for 30-45 minutes at 37ºC before being diluted 1:9 with additional HEPES buffer. The mixture was passed through a 250 micron mesh before being centrifuged at 300 rcf for 6 minutes at room temperature. The top lipid layer was removed and the adipocytes collected for snap-freezing in liquid nitrogen.

RNA isolation
RNA was extracted from VA and SA sample types following the RNeasy Lipid Tissue Mini Kit (Qiagen, Cat#74804), as we have described before. The protocol was followed exactly except that homogenisation used the probe tissue ruptor, and on-column DNase digestion was carried out. RNA quality and quantity was determined on the Agilent Bioanalyser.

DNA isolation
DNA from the VA and SA sample types was extracted using an in-house method. For every 200 -350 µL of frozen sample, 1 mL of TP lysis buffer (Tris 50mM pH 7.5, NaCl 0.1M, SDS 0.5%, EDTA 5mM) and 100 µL of proteinase K (20 mg/ mL, Promega, Cat#V3021) was applied. This solution was then incubated overnight at 55ºC with constant shaking (700 rpm). The next day, after centrifugation at 2,000 rpm, the lipid layer was discarded and an equal volume of phenol: chloroform: isoamyl alcohol was applied. After centrifugation the aqueous phase was collected and precipitated in 0.3M sodium acetate, 66% v/v ethanol and 100 mg/ mL Glycoblue (Life Tech, cat#AM9515) while incubating in dry ice for 30 minutes. The DNA was pelleted and washed in 75% ethanol before being resuspended in 100 µL milliQ water.
VAT DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Cat#69504) following an amended protocol for total DNA from animal tissue. Briefly, 30 mg pieces of tissue were homogenised in 360 µL of ALT buffer using the bead beater for 2 minutes set at full speed, and a 5 mm stainless steel bead (Qiagen, Cat#69989). 40 µL of Proteinase K was added and the solution incubated at 56ºC overnight with constant shaking (500 -700 rpm). 400 µL of Buffer AL was mixed with the sample the following day and this mixture was passed through the DNeasy Mini spin column. The adherent DNA was then washed with Buffers AW1 and AW2 and then eluted in 200 µL of Buffer AE.
DNA was extracted from PBL sample types following the Gentra Puregene Blood Cell Kit (Qiagen, Cat#158445).

Whole Genome Bisulfite Sequencing (WGBS)
Whole genome bisulfite sequencing libraries were prepared following Illumina's "Whole-Genome Bisulfite sequencing for Methylation Analysis" protocol. Briefly, 1 µg of genomic DNA was spiked with 0.5% unmethylated lambda DNA and sonicated to generate fragments of size between 150 to 300bp. Library preparation was performed for using the Illumina's Paired-end DNA Sample Prep Kit (discontinued, Illumina, CA, USA) according to the manufacturer's protocol. The size selected libraries were then subjected to bisulfite conversion as previously described (Clark et al. 2006). Adaptor-ligated bisulfite treated DNA was enriched by 10 cycles of PCR amplification using the PfuTurbo Cx Hotstart DNA Polymerase (Stratagene). Qualitative and quantitative checks of the libraries were performed using Agilent's High sensitivity DNA kit (Agilent) and KAPA Library quantification kit (KAPA Biosystems). Three lanes of paired end 100bp sequencing was performed for each of the library on the Illumina HiSeq2500 platform using the TruSeq v3 cluster kits and SBS kits to achieve coverage ranging between 25-30x.

Poly-A RNA-seq
RNA-seq libraries were prepared using Illumina TruSeq Stranded poly-A RNA Library Prep Kit by the University of Western Sydney Next Generation Sequencing Facility. 500 ng of RNA was used for each sample. Libraries were sequenced using 100bp paired-end HiSeq2500 chemistry3 by UWS, yielding at least 30 million reads per sample, Supplementary Table S1.
Whole RNA-seq RNA-seq libraries were prepared from 500 ng RNA, following ribosomal RNA reduction using Illumina TruSeq Stranded Total RNA Library Prep Kit by the Australian Genome Research Facility (AGRF). Libraries were sequenced using 100bp paired-end HiSeq2500 chemistry4 by AGRF, yielding 35 to 40 million aligned reads per sample, Supplementary Table S1.

Illumina 450K analysis
Briefly normalization was performed using the dasen method from watermelon (Pidsley et al. 2013). Some probes were excluded from the analysis. In particular these were, the union of probes on the sex chromosomes (n = 11,648), probes targeting CpGs located two or fewer nucleotides from a known single nucleotide polymorphism (SNP) with a minor allele frequency >0.05 (n = 29,476), and known cross-hybridizing probes (n = 30,969) published by Chen et al (Chen et al. 2013). Probes were also excluded if they failed quality control metrics in one or more samples, based on a detection P value >0.05.Probe M-values were used to detect differentially methylated probes via Limma (Ritchie et al. 2015). Differentially methylated region calling was performed with DMRcate (Peters et al. 2015), using default parameters.

Calling methylation changes between groups of samples
DMRs were called using the bsseq Bioconductor package (Hansen et al. 2012). Smoothing was performed using the parameters "ns=70, h=1000" and only CpG sites with minimum 2x coverage in at least 2 samples in each group being compared was used. T-stats were called using the parameters "estimate.var='paired', local.correct=TRUE" and the DMRs called using a t-stat cutoff of 4 and filtered for minimum 3 CpG sites and an absolute change in average methylation of 10%. Given our sequencing depth across 3 biological replicates 10% DMRs should have a true positive rate (TPR) of >80% and false discovery rate (FDR) < 10%, while 20% DMRs TPR >90% and FRD <5% (Ziller et al. 2015).

Plotting methylation data
Aggregated DNA methylation profiles were plotted using Loess regression. Default settings except span = 0.2 was used.

Alignment/Mapping & Counting:
Adapters from the sequencing reads were trimmed using Trimmomatic (v0.33) (Bolger et al. 2014). Fastq sequences were then quality checked using FastQC bioinformatics tool (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). Reads were then mapped to the human genome (hg19) using custom bioinformatics pipeline NGSANE (Buske et al. 2014), which utilizes TopHat (Trapnell et al. 2013) (TopHat v2.0.6) to perform the alignment in a strand specific manner. It is to be noted that TopHat was run using default parameters with options set as "fr-firststand" to specify library type and "-transcriptome-index" set to known transcriptome built from RefSeq.
Following the mapping process, htseqcount feature from the bioinformatics tool htseq (Anders et al. 2015) was employed to count the number of uniquely mapped reads of genomic features such as genes and transcripts guided by the GENCODE GTF (v19) (Harrow et al. 2012) annotation with the following parameter (HTSEQCOUNT_UNIQUE = 1 and GTF = "gencode.v19.annotation.gtf").

Differential Expression:
Differential expression analysis for subcutaneous adipocytes vs visceral adipocytes was performed using edgeR bioconductor package McCarthy et al. 2012). Previous studies have suggested that edgeR is better suited for evaluating differential expression between biological conditions utilizing read count of genomic regions of interest, such as genes and transcripts McCarthy et al. 2012). The generalized linear model (GLM) algorithm of edgeR is able to separate out technical variability and extract feature-wise variation between condition even only using a few biological replicates (50,51). The topTag() function of edgeR was used to identify differentially expressed (DE) genes/transcripts with a Benjamin-Hochberg (BH) FDR correction method. Only genes/transcripts with FDR value <0.05 have been considered as DE genes/transcripts.

GO Analysis:
Typically Gene Ontology analysis (GO) is used to understand the underlying biological phenomenon/functions. For this study we have used goseq (Young et al. 2010) to perform gene functional enrichment analysis as it has been shown to discount for gene length bias that may have affected any differential expression call. Goseq typically identifies the over-represented GO terms that are statistically significant having a p-value <0.05. We then have further applied BH method to account for multiple testing corrections and only have considered the GO terms with adjusted FDR < 0.05 for further investigation and biological implication. In order to aid biological interpretation of terms only GO terms with between 15-500 genes in the category were used for generation of bar plots.
For GO analysis on methylation data GO-seq was modified to account for CpG density of regions tested. Only DMRs within the promoters (2kb either side of TSS) or genebody (-2kb of TSS up until transcriptional end site) were used in GO analyses. To objectively summarise DMR GO terms tree plots were generated using REViGO (13). The top 100 GO terms were used as input for this process. Table  Number Title  S1 Subject characteristics and genomic data S2