DNA methylation and gene expression patterns in adipose ... - Nature

International Journal of Obesity (2016) 40, 654–661 © 2016 Macmillan Publishers Limited All rights reserved 0307-0565/16 www.nature.com/ijo

ORIGINAL ARTICLE

DNA methylation and gene expression patterns in adipose tissue differ significantly within young adult monozygotic BMI-discordant twin pairs KH Pietiläinen1,2,3,11, K Ismail4,11, E Järvinen1, S Heinonen1, M Tummers1, S Bollepalli4, R Lyle5, M Muniandy1,4, E Moilanen6, A Hakkarainen7, J Lundbom7,8,9, N Lundbom7, A Rissanen1, J Kaprio3,4,10 and M Ollikainen4 BACKGROUND: Little is known about epigenetic alterations associated with subcutaneous adipose tissue (SAT) in obesity. Our aim was to study genome-wide DNA methylation and gene expression differences in SAT in monozygotic (MZ) twin pairs who are discordant for body mass index (BMI). This design completely matches lean and obese groups for genetic background, age, gender and shared environment. METHODS: 14We analyzed DNA methylome and gene expression from SAT, together with body composition (magnetic resonance imaging/spectroscopy) and glucose tolerance test, lipids and C-reactive protein from 26 rare BMI-discordant (intrapair difference in BMI ⩾ 3 kg m − 2) MZ twin pairs identified from 10 birth cohorts of young adult Finnish twins. RESULTS: We found 17 novel obesity-associated genes that were differentially methylated across the genome between heavy and lean co-twins. Nine of them were also differentially expressed. Pathway analyses indicated that dysregulation of SAT in obesity includes a paradoxical downregulation of lipo/adipogenesis and upregulation of inflammation and extracellular matrix remodeling. Furthermore, CpG sites whose methylation correlated with metabolically harmful fat depots (intra-abdominal and liver fat) also correlated with measures of insulin resistance, dyslipidemia and low-grade inflammation, thus suggesting that epigenetic alterations in SAT are associated with the development of unhealthy obesity. CONCLUSION: This is the first study in BMI-discordant MZ twin pairs reporting genome-wide DNA methylation and expression profiles in SAT. We found a number of novel genes and pathways whose methylation and expression patterns differ within the twin pairs, suggesting that the pathological adaptation of SAT to obesity is, at least in part, epigenetically regulated. International Journal of Obesity (2016) 40, 654–661; doi:10.1038/ijo.2015.221

INTRODUCTION Adipose tissue is a plastic organ with several unique properties, enabling it to expand and contract in response to alterations in energy balance. During times of excess energy intake, adipose tissue may enlarge by increasing the size (hypertrophy) and number (hyperplasia) of adipocytes, the latter being the metabolically healthier adiposity phenotype, with high insulin sensitivity, low liver fat and absence of low-grade inflammation.1 Through mechanisms that are incompletely understood, adipogenesis, the differentiation of new adipocytes from adipose stem cells, is especially impaired in hypertrophic obesity but retained in hyperplastic obesity.2 Healthy expansion of adipose tissue requires modifications of the supporting extracellular matrix (ECM)3 and intense growth of the microvascular network.4 Inflammation and apoptosis, in contrast, have been suggested to be the hallmarks of adipose tissue dysfunction and closely connected to the development of metabolic disturbances in obesity.5 Epigenetic mechanisms such as DNA methylation and histone modifications are essential to tissue growth and remodeling, but

their role in subcutaneous adipose tissue (SAT) expansion in obesity remains largely unknown. The Infinium HumanMethylation450 BeadChip has been used to study correlations between body mass index (BMI) and DNA methylation in blood and SAT.6,7 Other studies have compared differences in methylation profiles in SAT between monozygotic (MZ) co-twins discordant for type 2 diabetes,8,9 and in visceral adipose tissue between men with and without metabolic syndrome.10 None of these previous studies addressed obesity per se. DNA methylation in SAT has been found to be moderately heritable (19%), with only a small fraction (0.2–2%) of variance being attributable to shared environmental factors.11 This finding indicates that a considerable part of DNA methylation variation is due to environmental factors and the lifestyle. Epigenetic analyses of weight-discordant MZ twin pairs represent a powerful tool to identify effects of acquired obesity that are independent of underlying genetic differences.12 In addition to the genome, MZ twins are completely matched for age, sex and family background. In the present study, we

1 Obesity Research Unit, Research Programs Unit, University of Helsinki, Helsinki, Finland; 2Endocrinology, Abdominal Center, Helsinki University Central Hospital, and University of Helsinki, Helsinki, Finland; 3Institute for Molecular Medicine FIMM, University of Helsinki, Helsinki, Finland; 4Department of Public Health, Faculty of Medicine, University of Helsinki, Helsinki, Finland; 5Oslo University Hospital and University of Oslo, Department of Medical Genetics, Oslo, Norway; 6The Immunopharmacology Research Group, University of Tampere School of Medicine and Tampere University Hospital, Tampere, Finland; 7HUS Medical Imaging Center, Radiology, University of Helsinki and Helsinki University Hospital, Helsinki, Finland; 8Institute for Clinical Diabetology, German Diabetes Center, Leibniz Center for Diabetes Research, Heinrich Heine University, Düsseldorf, Germany; 9German Center for Diabetes Research, Partner Düsseldorf, Düsseldorf, Germany and 10National Institute for Health and Welfare, Helsinki, Finland. Correspondence: Dr M Ollikainen, Department of Public Health, PO Box 41, University of Helsinki, FI-00014 Helsinki, Finland. E-mail: miina.ollikainen@helsinki.fi 11 These two authors shared first authorship. Received 4 July 2015; revised 29 August 2015; accepted 21 September 2015; accepted article preview online 26 October 2015; advance online publication, 1 December 2015

Adipose tissue methylation and gene expression in acquired obesity KH Pietiläinen et al

655 performed genome-wide DNA methylation analyses in SAT integrated with data on gene expression on a deeply phenotyped, rare sample of 26 healthy young adult BMI-discordant MZ twin pairs. We report novel genes and pathways that are differentially methylated within the BMI-discordant twin pairs, describe the relationship between DNA methylation and gene expression as well as the localization of the respective proteins in SAT. MATERIALS AND METHODS Subjects and study design This study included 26 BMI-discordant (within-pair difference in BMI (ΔBMI) ⩾ 3 kg m − 2) and 11 BMI-concordant (ΔBMI o3 kg m − 2) MZ twin pairs identified from 10 full birth cohorts of Finnish twins (n = 5417 families).13 In the discordant pairs, mean ± s.d. ΔBMI was 5.97 ± 2.78 and range was 3.06– 15.15 kg m − 2. BMI of the leaner co-twins was 25.28 ± 4.52 kg m − 2 and of the heavier co-twins was 31.25 ± 5.18 kg m − 2 (Supplementary Table 1). Percent body fat was measured by dual-energy X-ray absorptiometry, amount of SAT and intra-abdominal adipose tissue by magnetic resonance imaging, liver fat by proton magnetic resonance spectroscopy, fasting lipids, high-sensitivity C-reactive protein, leptin, adiponectin as well as glucose and insulin during a 2-h oral glucose tolerance test for the calculation of homeostatic model assessment (HOMA)-insulin resistance and Matsuda-insulin sensitivity indexes, dietary intake, smoking, alcohol use and physical activity (Supplementary Table 1) as previously described.5 SAT biopsies were available from all twins for gene expression analyses from 24 discordant and all concordant pairs for DNA methylome analyses and from 14 discordant pairs for isolation of adipocytes and subsequent gene expression analyses. The 14 pairs with adipocyte material available had the following characteristics: ΔBMI 5.90 ± 2.14, range 3.32–9.95 kg m − 2, BMI of the leaner co-twins 24.07 ± 3.27 kg m − 2 and of the heavier co-twins 29.92 ± 3.65 kg m − 2. In the 12 pairs where adipocytes were not isolated (because of limited sample material or lack of laboratory personnel during the study day), ΔBMI was 6.11 ± 3.49, range 3.06–15.15 kg m − 2, BMI of the leaner co-twins 26.70 ± 5.46 kg m − 2 and of the heavier co-twins 32.81 ± 6.35 kg m − 2, not significantly different from the pairs with adipocytes available. Only gene expression analysis was performed on the adipocytes because, unfortunately, there were not enough isolated adipocytes available for DNA extraction and for methylome analyses. All subjects provided written informed consent. The protocols were approved by the Ethics Committee of the Hospital District of Helsinki and Uusimaa.

SAT biopsies and DNA and RNA extraction SAT biopsies were collected under local anesthesia from the periumbilical area by a surgical technique, snap frozen in liquid nitrogen, and stored at − 80 ºC. High-molecular-weight DNA and total RNA were extracted using the QIAamp DNA Mini kit and the RNeasy Lipid Mini Kit (QIAGEN Nordic, Sollentuna, Sweden), respectively, according to the manufacturer’s instructions. A piece of SAT was used for the isolation of pure adipocytes that were used for calculation of adipocyte size, RNA extraction and gene expression analyses.1 RNA quality from whole SAT and isolated adipocytes was assayed by 2100 Bioanalyzer (Agilent Technologies, Espoo, Finland).

DNA methylation analysis DNA extracted from SAT was bisulfite converted using the EZ-96 DNA Methylation-Gold Kit (Zymo Research, Irvine, CA, USA) according to the manufacturer’s instructions. The co-twin samples were always converted on the same plate to minimize potential batch effects. DNA methylation status was assessed using the Infinium HumanMethylation450 BeadChip (450k) according to the manufacturer’s instructions (Illumina, San Diego, CA, USA). Preprocessing, normalization and statistical analysis of the data were performed using R-3.02(ref. 14) and Bioconductor 2.13(ref. 15) as in our previous study, where we also validated the 450k by EpiTYPER assays.16 The data were quantile normalized and adjusted for probe type using Beta-Mixture Quantile normalization (BMIQ).17 Unreliable probes (detection P-values 40.001), probes mapping to sex chromosomes or multiple locations on the genome were removed from the data along with probes with low variance, leaving 218 590 probes for analysis. All analyses were done using M values, adjusted for smoking, and β-values were used to report the outcomes. The data are available in the Gene © 2016 Macmillan Publishers Limited

Expression Omnibus database repository, www.ncbi.nlm.nih.gov/geo/ query/acc.cgi?acc = GSE68336.

Gene expression analyses and immunohistochemistry Total RNA from SAT specimens and isolated adipocytes was used to study gene expression by the Affymetrix U133 Plus 2.0 chips (Affymetrix, Vienna, Austria) that have been previously validated.1 The data were preprocessed and gcRMA normalized using limma (Bioconductor 2.13). We further validated the gene expression results by real-time quantitative PCR from commercial RNA of obese (BMI 32.1 kg m − 2) and lean individual’s SAT (RNA-T10-3 and R1234003-10, Amsbio, Abingdon, UK). Immunohistochemistry was performed on lean individual’s SAT section (Supplementary Methods).

Statistical analyses Statistical tests were conducted in R-3.02. Clinical parameters between the co-twins were analyzed using Wilcoxon’s signed rank tests. Similarity in DNA methylation between individuals was computed using Pearson’s correlation and Euclidean distance. Differential methylation and differential expression analyses in twin SAT were performed using paired moderated t-tests (limma), and P-values adjusted using the Benjamini Hochberg method. Correlations between DNA methylation, gene expression and clinical measures were computed using Spearman’s correlations, adjusted for family. Gene Set Analysis (GSA) was done pair-wise with 1000 permutations using the R package GSA. Groups of CpG sites (probe sets), representing Reactome pathways, were used as ‘gene sets’ in the algorithm. A P-value cutoff of 0.01 and a false discovery rate cutoff of 0.25 were applied to obtain the list of significant pathways. Gene expression differences in adipocytes were analyzed using paired t-tests.

RESULTS Metabolic characteristics Clinical characteristics of the twins are shown in Supplementary Table 1. In the 26 BMI-discordant MZ pairs, the heavy co-twins weighed 17.96 kg (5.97 kg m − 2) more than the lean co-twins. In addition to having higher adiposity measures, the heavy co-twins were more insulin resistant, had lower high-density lipoprotein cholesterol and adiponectin, higher low-density lipoprotein cholesterol, triglyceride and leptin levels and higher systolic and diastolic blood pressure than their lean co-twins. The discordant co-twins did not differ for energy intake or physical activity. In the 11 BMI-concordant pairs, there were subtle differences in BMI and SAT volume, but no differences in the metabolic characteristics. For these reasons, the concordant pairs were not used in any within-pair analyses. Genome-wide DNA methylation differences We first determined the reliability and consistency of the 450k, and observed that genome-wide DNA methylation profiles of SAT in MZ twin pairs were more similar (Euclidean distance = 19.59) than in unrelated twin individuals (Euclidean distance = 28.12) (Supplementary Figure 1). The subsequent analyses comparing the BMI-discordant co-twins revealed that methylation profiles were highly similar within the twin pairs (r40.985, P o2.2 10 − 16), and obesity was not associated with global hypo- or hypermethylation of SAT. There were, however, significant differences between the co-twins at individual CpG sites across the genome. Altogether, 22 CpGs were differentially methylated (false discovery rate o0.05), 12 (54.5%) of them associated with CpG islands and 17 within known genes (Table 1). Although only a single CpG site (1/17, 5.8%) was in the promoter region, the majority (12/17, 70.6%) were located in gene bodies (Supplementary Figure 2). Half of the differentially methylated CpGs were hypermethylated and most such CpGs were located in the gene bodies (8/11, 72.7%). Hypomethylation affected mostly (5/11, 45.5%) intergenic regions (Supplementary Figure 2). International Journal of Obesity (2016) 654 – 661


656 Table 1.

Differentially methylated CpG sites in subcutaneous adipose tissue of 24 BMI-discordant monozygotic twin pairs

Gene symbol CHST11 AXIN2 ZBTB16 NA NA MAD1L1 E2F5 FGFRL1 CCDC92 ACSF3 NA MAML3 BAG6 ASAP2 SORBS3 EHBP1L1 RBPMS NA TBC1D16 SERPINF1 MRPL23 NA

ENTREZ ID

Probe ID

P-value

FDR

Methylation status in heavy co-twinsa

Relation to genes

Relation to CpGi

50515 8313 7704 NA NA 8379 1875 53834 80212 197322 NA 55534 7917 8853 10174 254102 11030 NA 125058 5176 6150 NA

cg00608661 cg06405341 cg22130673 cg07611933 cg12192282 cg17187143 cg12112870 cg07727358 cg15258936 cg08551036 cg25463399 cg16218705 cg04536765 cg23902076 cg24805258 cg10142520 cg05965490 cg00874032 cg20402747 cg17141696 cg23314972 cg00960147

5.80 × 10 − 7 8.05 × 10 − 7 8.46 × 10 − 7 1.02 × 10 − 6 1.14 × 10 − 6 1.31 × 10 − 6 1.50 × 10 − 6 1.58 × 10 − 6 1.67 × 10 − 6 1.92 × 10 − 6 2.00 × 10 − 6 2.39 × 10 − 6 2.56 × 10 − 6 2.78 × 10 − 6 2.82 × 10 − 6 3.41 × 10 − 6 3.89 × 10 − 6 4.19 × 10 − 6 4.27 × 10 − 6 4.32 × 10 − 6 4.51 × 10 − 6 4.61 × 10 − 6

0.040 0.040 0.040 0.040 0.040 0.040 0.040 0.040 0.040 0.040 0.040 0.041 0.041 0.041 0.041 0.046 0.046 0.046 0.046 0.046 0.046 0.046

Hypo Hyper Hyper Hypo Hypo Hyper Hypo Hyper Hyper Hypo Hypo Hyper Hypo Hypo Hyper Hypo Hyper Hypo Hyper Hyper Hyper Hypo

Body Body Body Intergenic Intergenic Body Body Body 3′UTR Body Intergenic Body Promoter Body Body Body Body Intergenic Body 3′UTR Body Intergenic

Open Sea S Shore Open Sea S Shore Open Sea Island Open Sea S Shore Island S Shore N Shelf Open Sea S Shore Open Sea Island S Shore Open Sea Open Sea N Shore Open Sea Island Open Sea

Abbreviations: BMI, body mass index; CpGi, CpG island; CGI Shelves, 2 to 4 kb from CGI; CGI Shores, o2 kb from CGI; FDR, false discovery rate; Hyper, more methylation in the heavy co-twins; Hypo, less methylation in the heavy co-twins; Open sea, isolated CpGs outside any CGIs; UTR, untranslated region. a Methylation in heavy compared with lean co-twins.

Expression of the differentially methylated genes Next, we explored whether the DNA methylation differences we identified contribute to gene function by analyzing the expression of the differentially methylated genes. Among the 17 genes, 9 (CHST11, ZBTB16, E2F5, FGFRL1, ASAP2, EHBP1L1, RBPMS, TBC1D16 and MRPL23) were differentially expressed within the twin pairs (P o 0.05, Table 2). We then validated the twin result by real-time quantitative PCR of unrelated obese and lean individuals’ SAT RNA. In the unrelated individuals (Supplementary Figure 3) as well as in the twins (Table 2), CHST11, ASAP2, TBC1D16, EHBP1L1 and MRPL23 were upregulated, whereas ZBTB16, E2F5, FGFRL1 and RBPMS were downregulated in obese compared with lean SAT. Correlation analyses between methylation and expression in twins revealed that ASAP2 gene body hypomethylation was associated with increased expression (r = − 0.40, P = 6.3 × 10 − 4), FGFRL1 (r = − 0.59, P = 1.4 × 10 − 7) and RBPMS (r = − 0.50, P = 8.7 × 10 − 6) gene body hypermethylation with decreased expression and TBC1D16 gene body hypermethylation with increased expression (r = 0.50, P = 1.6 × 10 − 5). CHST11, ZBTB16, FGFRL1 and RBPMS are induced during adipogenesis To address the possible role of the nine differentially methylated and expressed genes in adipogenesis, we reanalyzed gene expression data from human preadipocytes undergoing differentiation, published by Mikkelsen et al.18 In these data, we observed the clearest increase in gene expression for ZBTB16 and RBPMS during the induction of adipogenesis (Supplementary Figure 4). For validation purposes, we also included PPARγ, the master regulator of adipogenesis, in the analyses and, as expected, observed that it was induced early and stayed upregulated during the maturation of adipocytes. ZBTB16 and RBPMS had similar expression patterns with PPARγ (Supplementary Figure 4). To this end, we checked the expression of PPARγ in our materials: it was downregulated in obese compared with lean SAT in twins (P = 0.003) and in International Journal of Obesity (2016) 654 – 661

Table 2.

Gene expression differences of the differentially methylated genes in subcutaneous adipose tissue of 26 BMI-discordant monozygotic twin pairs Entrez ID

Gene

Log fold change in heavy vs lean co-twins

125058 53834 8853 11030 1875 254102 6150 50515 7704 8313 197322 7917 8379 10174 55534 5176 80212

TBC1D16 FGFRL1 ASAP2 RBPMS E2F5 EHBP1L1 MRPL23 CHST11 ZBTB16 AXIN2 ACSF3 BAG6 MAD1L1 SORBS3 MAML3 SERPINF1 CCDC92

0.52 − 0.36 0.17 − 0.29 − 0.11 0.15 0.12 0.33 − 0.26 0.35 − 0.04 − 0.03 0.12 0.04 0.00 0.02 0.02

P-value

FDR

9.04 × 10 − 6 1.54 × 10 − 4 3.67 × 10 − 5 3.12 × 10 − 4 7.79 × 10 − 4 4.41 × 10 − 3 0.004 0.018 0.014 0.040 0.014 0.040 0.029 0.067 0.032 0.067 0.048 0.090 0.059 0.101 0.079 0.123 0.132 0.187 0.192 0.251 0.364 0.442 0.590 0.631 0.594 0.631 0.750 0.750

Abbreviations: BMI, body mass index; FDR, false discovery rate.

unrelated individuals (Supplementary Figure 3) and in twin adipocytes (P = 0.0305), consistently suggesting decreased adipogenesis in obesity. Differentially expressed proteins are located in multiple cell types in SAT SAT is a heterogeneous collection of cells, the proportions of which may change during the development of obesity. Therefore, we wanted, first, to explore whether the detected within-pair expression differences arise from adipocytes by analyzing gene expression data from purified adipocytes of the twins and, second, © 2016 Macmillan Publishers Limited


657

Figure 1. Genes differentially methylated between heavy and lean co-twins of the BMI-discordant pairs are present in many cell types in adipose tissue. CHST11 and ZBTB16 are located in all three cell types: adipocytes, macrophages and endothelial cells. RBPMS is present in adipocytes and endothelial cells. The different cell types were identified by cell morphological analysis by using hematoxylin and eosin staining, as well as by using known markers for adipocytes (PPARγ), macrophages (CD163) and endothelial cells (CD31). Scale bar: 10 μm.

to carefully identify the localization of the proteins encoded by our differentially methylated and expressed genes by immunohistochemistry. Of the nine genes that were differentially methylated and expressed within twin pairs in the SAT, three were also differentially expressed in isolated adipocytes (P o 0.05, Supplementary Table 2). ASAP2 was upregulated and ZBTB16 and RBPMS were downregulated in the heavy co-twins. The direction of these genes’ expression differences in the isolated adipocytes was the same as in the whole SAT (Table 2 and Supplementary Figure 3). To localize proteins of interest in the different cell types, we performed immunohistochemistry on sections of whole SAT (Figure 1). We chose the two downregulated transcription factors, ZBTB16 and RBPMS, as well as CHST11, the top differentially methylated (and upregulated) gene in the heavy compared with lean co-twins. ZBTB16 was localized sporadically in adipocyte nuclei and in endothelial cells, but was strongly present in macrophages. RBPMS staining was strong in adipocyte nuclei and endothelial cells, and it was the only studied protein not present in macrophages. CHST11 was localized in adipocyte cytoplasm and the membranes of adipocytes, macrophages and endothelial cells. Correlation of DNA methylation with clinical measures In an attempt to understand how the epigenetic variants in SAT link to clinical outcomes, we correlated the methylation values of the differentially methylated 22 CpG sites with the clinical measures in the individual twins. The results strengthened the findings of the methylation analyses of the BMI-discordant twin pairs: adiposity and its related measures were clearly associated with the differentially methylated CpG sites. All CpG sites correlated significantly with the measures of overall fatness (BMI, %body fat); the distribution of fat to the subcutaneous, intraabdominal and liver depots; adipocyte size; and leptin (Figure 2). Most CpG sites also correlated with triglycerides, insulin resistance (insulin, HOMA) and inflammation (high-sensitivity C-reactive protein). Opposite correlations were observed for beneficial (high-density lipoprotein cholesterol, Matsuda-index and adiponectin) compared with harmful clinical measures. © 2016 Macmillan Publishers Limited

Differential methylation at the pathway level In order to visualize the data in the context of biological networks, we applied GSA to the SAT methylation data in the BMI-discordant pairs. The analyses identified 26 differentially methylated pathways within the pairs (P o 0.01, false discovery rate o 0.25; Table 3 and Figure 3), from which two large networks were identified. The first, which we called lipid metabolism, included nine pathways (Table 3 and Figure 3); the gene expression of all but one of them was downregulated in the heavy co-twins (Table 3 and Supplementary Table 3). Although the downregulated pathways are mainly involved in lipogenesis, adipogenesis, fatty acid and triglyceride biosynthesis, the upregulated ABC transporter pathway also has a role in macrophage lipid metabolism.19 The second large entity was inflammation, with all nine pathways upregulated in the heavy co-twins (Table 3, Figure 3 and Supplementary Table 3). The rest of the differentially methylated pathways formed three smaller entities. Dynamic and continuous remodeling of the adipose tissue ECM20 was represented by three intertwined pathways (Figure 3) that were upregulated in the heavy co-twins (Table 3 and Supplementary Table 3). Vitamin metabolism consisted of three pathways, including two vitamin B pathways. Signaling consisted of two differentially methylated pathways within the twin pairs that included several genes involved in insulin resistance. DISCUSSION We analyzed DNA methylation of 437 180 CpG sites across the genome in SAT from a rare sample of deeply phenotyped young MZ twin pairs discordant for BMI, and integrated the data of the differentially methylated genes with their transcription and protein expression. The discordant MZ twin design is the only human model allowing analysis of epigenetic alterations independent of genetics, age and sex, the significant confounders affecting epigenetic marks. Through within-pair genome-wide DNA methylation analysis, we identified 17 differentially methylated obesity-associated genes, of which 4 were significantly downregulated and 5 upregulated in heavy compared with lean International Journal of Obesity (2016) 654 – 661


658

Figure 2.

Correlations between the differentially methylated CpG sites and clinical measures in individual twins (n = 70).

co-twins. SAT methylation of all these 17 genes correlated with BMI in a recent study.7 Pathway analyses revealed genes clustering to lipid metabolism, inflammation and ECM remodeling to be differentially methylated, with lipid metabolism showing mostly downregulation and the other two clusters showing upregulation of gene expression in heavy compared with lean co-twins. These novel results add to the previously published data on downregulation of adipogenesis21 and lipid synthesis,22 upregulation of inflammation23 and ECM remodeling20,24,25 in the pathological changes of SAT in obesity by providing a potential molecular mechanism underlying these important phenomena. Obesity,26 its related complications27 and methylation28,29 are known to be heritable. Therefore, methylation analyses in obese and lean unrelated individuals are subject to significant genetic confounding. For the same reason, MZ twins discordant for phenotype and epigenetic marks are extremely difficult to find. In the present study, we were able to identify 26 such pairs among 5417 young twin pairs, and with the detailed metabolic and molecular data, our study represents arguably the most deeply characterized modern collection of MZ twins in relation to obesity and its comorbidities. As previous methylation studies in MZ twins on blood and buccal cells,29–31 we found twins to be more similar than unrelated individuals for their genome-wide methylation levels in SAT. Despite the high within-pair similarity in DNA methylation and no significant differences in the global methylation values, we identified 22 CpG sites (within 17 genes) and 26 pathways across the genome that were different between heavy and lean MZ co-twins. In line with our finding on global methylation levels, the differentially methylated CpG sites showed the same numbers hypomethylated (n = 11) and hypermethylated (n = 11) in heavy compared with lean co-twins. Our results thus indicate that rather than associating with global methylation levels, obesity is associated with methylation variation at International Journal of Obesity (2016) 654 – 661

individual CpG sites. This is in line with methylation analyses in visceral fat between men with and without metabolic syndrome10 and in SAT and muscle in MZ twins discordant for type 2 diabetes.8,9 In unrelated individuals, global DNA methylation level in SAT correlated with waist circumference32 but not with BMI.7,32 In the present study, most of the hypomethylation events occurred in intergenic regions, whereas hypermethylation most frequently affected gene bodies. No clear relationship was observed between hypermethylation of the gene body (or hypomethylation of the gene promoter) and upregulation of gene expression. Nonetheless, our finding is in line with a recent study showing that widespread relationships with significant correlations exist between DNA methylation and gene expression.33 The paradoxical reduction in adipocyte differentiation and lipogenesis in obesity is a known phenomenon,21,22 but little is known about its molecular mechanisms. In the present study, we identified genes that are likely related to this process. Of special interest was the gene body hypermethylation of ZBTB16, FGFRL1 and RBPMS in the heavy co-twins’ SAT, along with downregulation of their gene expression. RBPMS18 and FGFR34 signaling to FGFRL1 (ref. 35) are early indicators of adipogenesis. ZBTB16 is induced during differentiation of 3T3-L1 adipocytes.18 Our analysis of the previously published gene expression data18 revealed that the expression of ZBTB16 and RBPMS, together with PPARγ, increased during adipogenesis. The decrease in their gene expression in our heavy co-twins’ SAT thus suggests that the downregulation of these genes is linked to impaired adipogenesis in obesity. In addition, consistent downregulation of lipogenesis pathways was observed in the heavy co-twins, supporting the findings of the individual genes stated above: SAT seems to have marked impairments in performing its main function—storage of excess energy—in obesity. An inability to expand and handle excess nutrients leads to susceptibility to inflammation and adipocyte death in adipose © 2016 Macmillan Publishers Limited


659 Table 3. Differentially methylated Reactome pathways from the Gene Set Analysis (GSA) of subcutaneous adipose tissue of 24 BMI-discordant monozygotic twin pairs Pathways grouped by entities

P-valuea

Lipid metabolism ChREBP activates metabolic gene expression Regulation of cholesterol biosynthesis by SREBP (SREBF) Fatty acyl-CoA biosynthesis Regulation of signaling by NODAL Acyl chain remodeling of DAG and TAG Activation of gene expression by SREBF (SREBP) Triglyceride biosynthesis Linoleic acid (LA) metabolism ABC transporters in lipid homeostasis

FDRa

Overall methylation status in heavy co-twinsb

Overall gene expression in heavy co-twinsc

o0.001 o 0.001 o0.001 o 0.001 0.001 0.105 0.001 0.105 0.001 0.105 0.001 0.105 0.002 0.144 0.003 0.182 0.004 0.219

Hyper Hyper Hyper Hyper Hyper Hyper Hyper Hyper Hyper

Down Down Down Down Down Down Down Down Up

Inflammation Downregulation of TGF-β receptor signaling o0.001 o 0.001 PTM: gamma carboxylation, hypusine formation and arylsulfatase activation o0.001 o 0.001 TGF-β receptor signaling activates SMADs 0.002 0.144 Deactivation of the β-catenin transactivating complex 0.003 0.182 FasL/CD95L signaling o0.001 o 0.001 Lysosome vesicle biogenesis o0.001 o 0.001 Interleukin-1 processing o0.001 o 0.001 Trafficking and processing of endosomal TLR o0.001 o 0.001 Hyaluronan uptake and degradation o0.001 o 0.001

Hyper Hyper Hyper Hyper Hypo Hypo Hypo Hypo Hypo

Up Up Up Up Up Up Up Up Up

Remodeling of the adipose tissue ECM Other semaphorin interactions Cell–ECM interactions Regulation of cytoskeletal remodeling and cell spreading by IPP complex components

Hyper Hyper Hyper

Up Up Up

o0.001 o 0.001 o0.001 o 0.001 0.002 0.144

Vitamin metabolism Vitamin B5 (pantothenate) metabolism FMO oxidizes nucleophiles Vitamin B6 activation to pyridoxal phosphate

0.002 0.002 0.004

0.144 0.144 0.219

Hyper Hyper Hyper

Down Down Up

Signaling SOS-mediated signaling Signal attenuation

0.001 0.003

0.105 0.182

Hyper Hyper

Up Down

Abbreviations: BMI, body mass index; DAG and TAG, di-and triacylglycerols; ECM, extracellular matrix; FDR, false discovery rate; IPP complex, ILK, PINCH and parvin; PTM: post-translational modification; TGF, transforming growth factor; TLR, Toll-like receptor. GSA was run with 1000 permutations; FDR cutoff o0.25; P-value cutoff o0.01. aP-values and FDR from GSA of methylation analyses between heavy and lean co-twins of the BMI-discordant pairs. bOverall methylation in heavy compared with lean co-twins for the pathway. cOverall gene expression in heavy compared with lean co-twins for the pathway.

tissue.36 Our pathway analyses support this view: obesity was associated with differential methylation and upregulation in the inflammation and ECM remodeling pathways. In addition, one of the lipid metabolism pathways, ABC transporters in lipid homeostasis, was upregulated in the heavy co-twins. The genes in this pathway are abundantly expressed in macrophages.19 Therefore, the overexpression signal of this pathway may relate to inflammatory processes in the SAT in obesity. ECM modifications through removal of the dying adipocytes and the remaining lipid droplets are essential in the dynamic adaptation of SAT to obesity and clearly linked to immune cell infiltration and inflammation.37,38 Inflammation and ECM modification have long been known to change pathologically in obesity, and here we provide the first evidence that they may be epigenetically regulated already before development of overt metabolic diseases. Reassuringly, similar pathways (cellular growth, lipid metabolism and inflammation) were found differentially methylated in SAT between men with and without metabolic syndrome.10 An important finding from our present study is the strong correlation of methylation values and clinical phenotypes. These results suggest that SAT methylation patterns in obesity are closely linked with early impairments in parameters of metabolic syndrome, including distribution of fat to the abdominal region © 2016 Macmillan Publishers Limited

(and to the liver), insulin resistance, high triglycerides and low high-density lipoprotein levels.39 We also found adipocyte size as well as leptin and adiponectin, secretory products specific to adipocytes, to significantly correlate with all identified differentially methylated CpG sites. As SAT contains different cell types, each having different DNA methylation profiles, it is likely that the present and previous6–8,10,11,40 results on whole SAT methylation may in part reflect the cellular composition of the tissue. On the other hand, the benefit of whole tissue analysis is that it best reflects the in vivo characteristics and cell-to-cell interactions in the tissue under examination. Immunohistochemistry experiments in our study showed that the proteins coded by the differentially methylated genes (CHST11, ZBTB16 and RBPMS) were found in at least two different cell types of the SAT. However, it is also noteworthy that the downregulation of ZBTB16 and RBPMS and the upregulation of ASAP2 gene expression were verified in isolated adipocytes. These observations indicate that, at least for these genes, the expression differences in SAT are not simply due to cell composition differences within pairs, but reflect actual differences in adipocyte expression profiles within twin pairs. A limitation of the present study is that although we were able to study transcriptomics in the adipocytes, we did not have enough adipocytes from these twins to also perform DNA methylome International Journal of Obesity (2016) 654 – 661


660

Figure 3. Significantly differentially methylated pathways between heavy and lean co-twins of the 24 BMI-discordant pairs (GSA with Reactome pathways). Three main classes of pathways were found: lipid metabolism (yellow), inflammation (orange) and dynamic and continuous remodeling of the adipose tissue extracellular matrix (ECM) (aqua). Pathways associated with vitamin metabolism (green) and signaling (purple) were also identified.

analyses. Thus far, we are not aware of any studies that would have compared methylome patterns in isolated mature adipocytes and stroma-vascular cells of SAT in humans. Such studies remain important future endeavors, as obesity is known to affect the cellular composition of SAT, with adipocyte hypertrophy1 and an overrepresentation of immune cells.37 Our pathway results (decreased lipo/adipogenesis, increased inflammation and ECM modification) suggest that epigenetic regulation at the whole tissue level is, at least in part, associated with these pathologic morphological changes. Given the cross-sectional nature of our study, the direction of causality cannot be proven. We believe that only a minority of the observed methylation and expression differences preceded obesity, but are mainly due to a complex mixture of different metabolic and clinical parameters that are related to the phenotype differences between the lean and heavy co-twins. However, we cannot preclude the possibility that some of the findings are indicative of processes that precede the onset of weight gain. In conclusion, in a rare group of young adult MZ twin pairs discordant for BMI, we show that DNA methylation of 17 genes and 26 pathways in SAT are associated with increased adiposity. The identified patterns of DNA methylation across the genome, combined with gene expression findings, uniformly demonstrate the combination of different pathogenic changes that characterize SAT in obesity. The capacity of the adipocytes for adequate enlargement, especially storage of excess energy by lipogenesis, seems to be impaired, closely connected with increased remodeling of SAT, including a significant increase in inflammation. Our study strongly suggests that the attempts of SAT to adapt to increased energy intake in obesity are epigenetically regulated. From the association between clinical parameters and DNA International Journal of Obesity (2016) 654 – 661

methylation, we further suggest that epigenetic alterations in SAT may contribute to the development of metabolic syndrome in obesity. CONFLICT OF INTEREST The authors declare no conflict of interest.

ACKNOWLEDGEMENTS We thank the participants for their invaluable contributions to the study. The Obesity Research Unit team and the staff at the Finnish Twin Cohort Study are acknowledged for their help in the collection of the data. We also thank Barbara Every, BioMedical Editor, for English language editing. This study was supported by the Academy of Finland (Grants 251316, 100499, 205585, 118555, 141054, 266286 and 272376), The Sigrid Juselius Foundation, The University of Helsinki Research Funds, Helsinki University Hospital Research Funds, grants from Novo Nordisk, the Finnish Diabetes Research Foundation, the Jalmari and Rauha Ahokas Foundation, Orion Pharmos Foundation and Emil Aaltonen Foundation, the Academy of Finland Center of Excellence in Complex Disease Genetics (Grants 213506 and 129680), the Academy of Finland Center of Excellence in Research on Mitochondria, Metabolism and Disease (272376), EPITRAIN Innovative techniques and models to understand epigenetic regulation in the pathogenesis of common diseases (EPITRAIN - FP7-PEOPLE-2012-ITN, Grant Agreement 316758) and BioSHaRE-EU (Grant Agreement HEALTH-F4-2010-261433) funded by the European Union’s Seventh Framework Programme (FP7/2007-2013).

AUTHOR CONTRIBUTIONS KHP, JK, AR and MO conceived and designed the study. JK was responsible for the twin cohort data collection from which the study pairs were recruited. KHP collected the biological samples and performed the clinical investigations of the twins. KHP and MO coordinated the study. KHP, KI, EJ and MO wrote the manuscript. JK helped in drafting the manuscript and critically commented on it. KI performed the bioinformatics and the statistical analysis of the

© 2016 Macmillan Publishers Limited


661 methylation and expression data and the clinical measurements. EJ performed the immunohistochemical staining together with MT, as well as imaging and real-time quantitative PCR validation of the array results. MM performed the gene expression data analysis. SH helped in collecting the samples, extracted the adipocytes and their RNA and determined the adipocyte cell diameters. SB analyzed the public data sets and data on unrelated individuals. RL was responsible for the generation of the methylation data. NL, JL and AH performed MR imaging and spectroscopy and analyzed the images. All authors participated in discussions related to analysis and interpretation and have read and approved the final manuscript. MO is the guarantor of this work and as such had full access to the data and takes responsibility for the integrity of the data and the accuracy of the data analysis.

REFERENCES 1 Heinonen S, Saarinen L, Naukkarinen J, Rodriguez A, Fruhbeck G, Hakkarainen A et al. Adipocyte morphology and implications for metabolic derangements in acquired obesity. Int J Obes (Lond) 2014; 38: 1423–1431. 2 Gustafson B, Hammarstedt A, Hedjazifar S, Smith U. Restricted adipogenesis in hypertrophic obesity: the role of WISP2, WNT, and BMP4. Diabetes 2013; 62: 2997–3004. 3 Bolton K, Segal D, McMillan J, Jowett J, Heilbronn L, Abberton K et al. Decorin is a secreted protein associated with obesity and type 2 diabetes. Int J Obes (Lond) 2008; 32: 1113–1121. 4 Lemoine AY, Ledoux S, Larger E. Adipose tissue angiogenesis in obesity. Thromb Haemost 2013; 110: 661–668. 5 Naukkarinen J, Heinonen S, Hakkarainen A, Lundbom J, Vuolteenaho K, Saarinen L et al. Characterising metabolically healthy obesity in weight-discordant monozygotic twins. Diabetologia 2014; 57: 167–176. 6 Dick KJ, Nelson CP, Tsaprouni L, Sandling JK, Aissi D, Wahl S et al. DNA methylation and body-mass index: a genome-wide analysis. Lancet 2014; 383: 1990–1998. 7 Ronn T, Volkov P, Gillberg L, Kokosar M, Perfilyev A, Jacobsen AL et al. Impact of age, BMI and HbA1c levels on the genome-wide DNA methylation and mRNA expression patterns in human adipose tissue and identification of epigenetic biomarkers in blood. Hum Mol Genet 2015; 24: 3792–3813. 8 Ribel-Madsen R, Fraga MF, Jacobsen S, Bork-Jensen J, Lara E, Calvanese V et al. Genome-wide analysis of DNA methylation differences in muscle and fat from monozygotic twins discordant for type 2 diabetes. PLoS One 2012; 7: e51302. 9 Nilsson E, Jansson PA, Perfilyev A, Volkov P, Pedersen M, Svensson MK et al. Altered DNA methylation and differential expression of genes influencing metabolism and inflammation in adipose tissue from subjects with type 2 diabetes. Diabetes 2014; 63: 2962–2976. 10 Guenard F, Tchernof A, Deshaies Y, Perusse L, Biron S, Lescelleur O et al. Differential methylation in visceral adipose tissue of obese men discordant for metabolic disturbances. Physiol Genomics 2014; 46: 216–222. 11 Grundberg E, Meduri E, Sandling JK, Hedman AK, Keildson S, Buil A et al. Global analysis of DNA methylation variation in adipose tissue from twins reveals links to disease-associated variants in distal regulatory elements. Am J Hum Genet 2013; 93: 876–890. 12 Locke AE, Kahali B, Berndt SI, Justice AE, Pers TH, Day FR et al. Genetic studies of body mass index yield new insights for obesity biology. Nature 2015; 518: 197–206. 13 Kaprio J. The Finnish Twin Cohort Study: an update. Twin Res Hum Genet 2013; 16: 157–162. 14 R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing: Vienna, Austria, 2013. Available from https://www.R-project.org/. 15 Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol 2004; 5: R80. 16 Ollikainen M, Ismail K, Gervin K, Kyllonen A, Hakkarainen A, Lundbom J et al. Genome-wide blood DNA methylation alterations at regulatory elements and heterochromatic regions in monozygotic twins discordant for obesity and liver fat. Clin Epigenetics 2015; 739, eCollection 2015. 17 Teschendorff AE, Marabita F, Lechner M, Bartlett T, Tegner J, Gomez-Cabrero D et al. A beta-mixture quantile normalization method for correcting probe design bias in Illumina Infinium 450 k DNA methylation data. Bioinformatics 2013; 29: 189–196.

18 Mikkelsen TS, Xu Z, Zhang X, Wang L, Gimble JM, Lander ES et al. Comparative epigenomic analysis of murine and human adipogenesis. Cell 2010; 143: 156–169. 19 Schmitz G, Kaminski WE. ABC transporters and cholesterol metabolism. Front Biosci 2001; 6: D505–D514. 20 Henegar C, Tordjman J, Achard V, Lacasa D, Cremer I, Guerre-Millo M et al. Adipose tissue transcriptomic signature highlights the pathological relevance of extracellular matrix in human obesity. Genome Biol 2008; 9: R14–2008-9-1-r14. 21 Christodoulides C, Lagathu C, Sethi JK, Vidal-Puig A. Adipogenesis and WNT signalling. Trends Endocrinol Metab 2009; 20: 16–24. 22 Lodhi IJ, Yin L, Jensen-Urstad AP, Funai K, Coleman T, Baird JH et al. Inhibiting adipose tissue lipogenesis reprograms thermogenesis and PPARgamma activation to decrease diet-induced obesity. Cell Metab 2012; 16: 189–201. 23 Hornung V, Ablasser A, Charrel-Dennis M, Bauernfeind F, Horvath G, Caffrey DR et al. AIM2 recognizes cytosolic dsDNA and forms a caspase-1-activating inflammasome with ASC. Nature 2009; 458: 514–518. 24 Catalan V, Gomez-Ambrosi J, Rodriguez A, Fruhbeck G. Role of extracellular matrix remodelling in adipose tissue pathophysiology: relevance in the development of obesity. Histol Histopathol 2012; 27: 1515–1528. 25 Kruger RP, Aurandt J, Guan KL. Semaphorins command cells to move. Nat Rev Mol Cell Biol 2005; 6: 789–800. 26 Elks CE, den Hoed M, Zhao JH, Sharp SJ, Wareham NJ, Loos RJ et al. Variability in the heritability of body mass index: a systematic review and meta-regression. Front Endocrinol (Lausanne) 2012; 3: 29. 27 Lehtovirta M, Pietilainen KH, Levalahti E, Heikkila K, Groop L, Silventoinen K et al. Evidence that BMI and type 2 diabetes share only a minor fraction of genetic variance: a follow-up study of 23,585 monozygotic and dizygotic twins from the Finnish Twin Cohort Study. Diabetologia 2010; 53: 1314–1321. 28 van Dongen J, Ehli EA, Slieker RC, Bartels M, Weber ZM, Davies GE et al. Epigenetic variation in monozygotic twins: a genome-wide analysis of DNA methylation in buccal cells. Genes (Basel) 2014; 5: 347–365. 29 Gordon L, Joo JE, Powell JE, Ollikainen M, Novakovic B, Li X et al. Neonatal DNA methylation profile in human twins is specified by a complex interplay between intrauterine environmental and genetic factors, subject to tissue-specific influence. Genome Res 2012; 22: 1395–1406. 30 Bell JT, Tsai PC, Yang TP, Pidsley R, Nisbet J, Glass D et al. Epigenome-wide scans identify differentially methylated regions for age and age-related phenotypes in a healthy ageing population. PLoS Genet 2012; 8: e1002629. 31 Wong CC, Caspi A, Williams B, Craig IW, Houts R, Ambler A et al. A longitudinal study of epigenetic variation in twins. Epigenetics 2010; 5: 516–526. 32 Keller M, Kralisch S, Rohde K, Schleinitz D, Dietrich A, Schon MR et al. Global DNA methylation levels in human adipose tissue are related to fat distribution and glucose homeostasis. Diabetologia 2014; 57: 2374–2383. 33 Wagner JR, Busche S, Ge B, Kwan T, Pastinen T, Blanchette M. The relationship between DNA methylation, genetic and expression inter-individual variation in untransformed human fibroblasts. Genome Biol 2014; 15: R37. 34 Widberg CH, Newell FS, Bachmann AW, Ramnoruth SN, Spelta MC, Whitehead JP et al. Fibroblast growth factor receptor 1 is a key regulator of early adipogenic events in human preadipocytes. Am J Physiol Endocrinol Metab 2009; 296: E121–E131. 35 Amann R, Trueb B. Evidence that the novel receptor FGFRL1 signals indirectly via FGFR1. Int J Mol Med 2013; 32: 983–988. 36 Vidal-Puig A. Adipose tissue expandability, lipotoxicity and the metabolic syndrome. Endocrinol Nutr 2013; 60(Suppl 1): 39–43. 37 Cinti S, Mitchell G, Barbatelli G, Murano I, Ceresi E, Faloia E et al. Adipocyte death defines macrophage localization and function in adipose tissue of obese mice and humans. J Lipid Res 2005; 46: 2347–2355. 38 Sun K, Kusminski CM, Scherer PE. Adipose tissue remodeling and obesity. J Clin Invest 2011; 121: 2094–2101. 39 Salmenniemi U, Ruotsalainen E, Pihlajamaki J, Vauhkonen I, Kainulainen S, Punnonen K et al. Multiple abnormalities in glucose and energy metabolism and coordinated changes in levels of adiponectin, cytokines, and adhesion molecules in subjects with metabolic syndrome. Circulation 2004; 110: 3842–3848. 40 Gehrke S, Brueckner B, Schepky A, Klein J, Iwen A, Bosch TC et al. Epigenetic regulation of depot-specific gene expression in adipose tissue. PLoS One 2013; 8: e82516.

Supplementary Information accompanies this paper on International Journal of Obesity website (http://www.nature.com/ijo)

© 2016 Macmillan Publishers Limited

International Journal of Obesity (2016) 654 – 661