Altered pathways in methylome and transcriptome longitudinal analysis of normal weight and bariatric surgery women

DNA methylation could provide a link between environmental, genetic factors and weight control and can modify gene expression pattern. This study aimed to identify genes, which are differentially expressed and methylated depending on adiposity state by evaluating normal weight women and obese women before and after bariatric surgery (BS). We enrolled 24 normal weight (BMI: 22.5 ± 1.6 kg/m2) and 24 obese women (BMI: 43.3 ± 5.7 kg/m2) submitted to BS. Genome-wide methylation analysis was conducted using Infinium Human Methylation 450 BeadChip (threshold for significant CpG sites based on delta methylation level with a minimum value of 5%, a false discovery rate correction (FDR) of q < 0.05 was applied). Expression levels were measured using HumanHT-12v4 Expression BeadChip (cutoff of p ≤ 0.05 and fold change ≥2.0 was used to detect differentially expressed probes). The integrative analysis of both array data identified four genes (i.e. TPP2, PSMG6, ARL6IP1 and FAM49B) with higher methylation and lower expression level in pre-surgery women compared to normal weight women: and two genes (i.e. ZFP36L1 and USP32) that were differentially methylated after BS. These methylation changes were in promoter region and gene body. All genes are related to MAPK cascade, NIK/NF-kappaB signaling, cellular response to insulin stimulus, proteolysis and others. Integrating analysis of DNA methylation and gene expression evidenced that there is a set of genes relevant to obesity that changed after BS. A gene ontology analysis showed that these genes were enriched in biological functions related to adipogenesis, orexigenic, oxidative stress and insulin metabolism pathways. Also, our results suggest that although methylation plays a role in gene silencing, the majority of effects were not correlated.

entially methylated between obese women before BS and normal weight women, and most of them showed high methylation level in obese patients. Also, 666 CpGs (located in 495 genes) had methylation levels changed after RYGB. For the microarray gene expression analysis, a total of 600 genes were differently expressed (DEGs) between obese women before BS and normal weight patients. Indeed, 1,366 genes were differentially expressed between surgery times. Among these, 88 genes were up-and 512 genes were down-regulated. Figure 1 illustrates the Venn diagram showing the overlap of differently expressed genes (fold change >2, FDR p < 0.05) and differently methylated CpG sites (Δβ ± 5%, p < 0.001). When preoperative patients and normal weight women were compared, a total of six genes were, concomitantly, differentially methylated and expressed (Fig. 1A). Two of them (CCNL1 and SLC5A8) present their methylation and expression levels in the same direction: both lower in preoperative patients compared to normal weight women. However, four genes (TPP2, PSMG6, ARL6IP1 and FAM49B) presented methylation and expression levels in opposite direction: higher methylation and lower expression level in pre-surgery women compared to normal weight women. Also, these CpG sites were located at TSS200, TSS1500 and 5′UTR regions (gene promoter region).
Considering pre and postoperative times comparison, Venn diagram analysis identified a total of eight genes that were common in methylation and expression arrays (Fig. 1B). Of these, two genes (ZFP36L1 and USP32) were differentially methylated and expressed after RYGB. Methylation and expression levels of ZFP36L1 and USP32 genes were changed in opposite direction. For ZFP36L1, the methylation level decreased after surgery while its expression increased. In contrast, for USP32, methylation level increased and expression decreased after the surgery. Considering CpGs location, we observed that they were located in gene body and TSS1500 region.
The expression analysis between postoperative and normal weight individuals showed no differential gene expression, thus concomitant analysis between methylation and expression profile couldn't be performed. Table 1 summarizes the results of integrated analysis between DNA methylation and gene expression arrays, that is, gene concomitantly differently methylated and expressed between groups (obese and normal weight women) and periods (before and after BS).
To additional understanding about the biological relevance of the identified genes, GO analyses were performed. These analyses showed that these genes are related to MAPK cascade, NIK/NF-kappaB signaling, cellular response to insulin stimulus, proteolysis and others.

Discussion
We identified in the present study genes that are concomitant differentially methylated and expressed in leukocytes after BS. Relevantly, these genes were involved in molecular pathways related to obesity´s physiopathology. Also, comparing obese subjects (before BS) and normal weight women, we identified four genes that showed lower expression and higher methylation levels; which were related to adipogenesis, anti-satiety effect, oxidative stress, and insulin metabolism. Moreover, we evidenced that two genes related to proteolysis and adipogenesis, were epigenetically regulated by RYGB procedure.
Despite the fact that the impact of DNA methylation on gene expression seems to depend of cytosine methylated site and location 22,23 , results obtained here showed that promoter methylation (TSS200 and TSS1500) is associated with decreased expression; however, for specific CpG sites, intragenic methylation (gene body) also correlates with decreased gene expression. Grundberg et al. 24 , found that a large number of methylation and expression associations were positive, thus the increase in methylation levels was linked to an increase of corresponding gene expression. More interestingly, evidence suggested there an adjustment mechanism of body gene methylation levels for transcription regulation. In line of this, when promoter methylation is constant, increasing body methylation is associated with more repressed expression 25 . Thus, considering that some authors affirmed that the function of DNA methylation in intergenic and gene-body regions is less defined 24 , our results showed that regardless of the location of methylated cytosine (promoter region or body), DNA methylation promotes a reduction in gene expression.
In addition, it is important to highlight that the expression of two (CCNL1 and SLC5A8) and six genes (MED21, CD302, CLEC12B, MNAT1, CHD6 and CCNH) were not regulated by DNA methylation in obese versus normal weight and preoperative versus postoperative comparisons, respectively. Study with cancer cells showed that the presence of methylation does not always imply in gene silencing 26 . Thus, our results showed that despite methylation plays a role in gene silencing, the majority of effects were not correlated.
Another point to be discussed is that the present study is the first to identify an association between TPP2, PSMC6, ARL6IP1 and FAM49B genes with obesity. The TPP2 gene codifies a peptidase that has been related to an anti-satiety effect by degrading the cholecystokinin 8 hormone (CCK8) and adipogenesis stimulus, pointing to an important role of TPP2 in obesity 27,28 . Thus, TPP2 seems to have dual role in metabolic homoeostasis, affecting both feeding behavior and adipose tissue biology per se 28 . Moreover, TPP2 is also related to different isoforms of protein convertase gene family (PSCK) that are associated with obesity related traits (e.g PCSK1 and PCSK2) 29 .
On the other hand, ARL6IP1 and PSMC6 have been associated with insulin synthesis and pathology of diabetes, respectively 30,31 . Also, ARL6IP1 protein is associated phosphatidylethanolamine-binding family of proteins (PEBP1), which has important role in MAPK, NF-kappa B, and glycogen synthase kinase-3 (GSK-3) signaling pathways. PSMC6 gene encodes one of the ATPase subunits and its expression was inversely correlated with BMI 11 .
Lastly, FAM49B has been appointed as novel regulator of mitochondrial function and may thus be associated with oxidative stress and inflammation, however, up to now, no functional data about this protein has been published 32 . Furthermore, USP32 is a highly conserved gene and uncharacterized gene that its stable silencing caused a significant decrease in the proliferation and migration rate of cells 9 . For this reason, USP32 has been associated with growth rate of cancer cells 9,33 .
We also identified that bariatric surgery modify ZFP36L1 and USP32 methylation and expression levels in leukocytes. In culture 3T3-L1 preadipocytes, ZFP36L1 was associated with adipogenesis 34 by regulating adipogenesis rate, playing an important role in obesity development 35 . According to Tseng 36 , ZFP36L1 overexpression might repress adipogenesis at least by down-regulating PPARG2 expression. These findings are a proof that blood leukocytes are able to reflect the regulation of same adipose tissue biology-related genes as it was also demonstrated in previous reports 13,[37][38][39] . In the current research, we were unable to detect differences in methylation levels of previously identified obesity related CpG sites, probably because differences in the study design, sample size or the threshold used for selection of candidates. However, the current study adds new information to this issue. In www.nature.com/scientificreports www.nature.com/scientificreports/ fact, the identified genes were involved in adipose tissue-related pathways that were also observed in previous reports on obesity-related methylation profile 13,38,39 .
The use of DNA and RNA extracted from leukocytes is a limitation of this study because both, expression and methylation, are tissue specific molecular mechanisms. However, leukocytes samples are less invasive, more convenient and acceptable than biopsies of target tissues (i.e adipose tissue) to perform longitudinal assessments and more suitable in clinical practice. To search for obesity-associated epigenetic biomarkers for the diagnosis and management of the disease is a huge challenge as adipose tissue is inaccessible without surgery. As mentioned above, in the obesity field, relevant studies have been recently published providing evidence that blood cells can be used to identify robust and biologically relevant epigenetic variation related to BMI. Moreover, it was demonstrated that epigenetic biomarkers in blood could mirror age-related epigenetic signatures in biologically relevant target tissues such as adipose tissue. Therefore, blood leukocytes are suitable, minimally invasive biological sources to evaluate obesity signatures. On the other hand, even though the identified changes in methylation and expression pattern induced by bariatric surgery in the present study could be transient, long-term studies would indicate the degree of plasticity of the epigenomic changes associated with bariatric surgery and would establish if over a longer time, the methylation profiles returned to the obese-pre surgery levels. www.nature.com/scientificreports www.nature.com/scientificreports/ In conclusion, integrating analysis of DNA methylation and gene expression evidenced that there is a set of genes relevant to obesity that changed after BS. A gene ontology analysis showed that these genes were enriched in biological functions related to adipogenesis, orexigenic, oxidative stress and insulin metabolism pathways. Also, our results suggest that although methylation plays a role in gene silencing, the majority of effects were not correlated.

Methods
Study participants. The present study included 24 normal weight women (BMI: 22.5 ± 1.6 kg/m 2 , 36.9 ± 11.8 years) and 24 severe obese women (BMI: 43.3 ± 5.7 kg/m 2 , 36.9 ± 10.2 years) submitted to RYGB which methylation and gene expression profile had been analyzed previously 20,21 . All samples from obese women were collected at preoperative time and after 6 months of bariatric surgery. Normal weight women were evaluated once.
All obese participants underwent open performed RYGB surgery that consisted in creating a small gastric portion (30 to 50 mL) and an anastomosis of the gastric stump to the jejunum (both remaining loops measured about 100 cm). The surgical procedure was standardized at our hospital following the most usual pattern performed worldwide and well described 40 . There were not postoperative complications in the patients included in this study.
The study was conducted with the approval of the Hospital Ethics Committee and in agreement with the Declaration of Helsinki. Informed consent was obtained from all individual participants included in the study. All methods were carried out in accordance with relevant guidelines and regulations.
DnA methylation array analysis and data processing. As aforementioned 21 , genomic DNA was extracted from peripheral mononuclear blood cells (PMBC) using GE Health Care kit (Illustra blood genomic Prep Mini Spin kit) and this extracted DNA was stored at −80 °C until the next steps. DNA fragmentation or RNA contamination was analyzed by 1% agarose gel electrophoresis. Also, DNA was bisulfite converted using EZ DNA methylation kit Methylation-Gold (Zymo Research, CA, USA) according to the manufacturer's instructions and then immediately hybridized in BeadChip. Genome-wide methylation analysis was conducted using the Infinium Human Methylation 450 BeadChip (Illumina, San Diego, CA, USA). Beadchips were scanned with the Illumina iScanSQ system and DNA quality checks, bisulfite modification, hybridization, data normalization, and statistical filter were performed as described before 21 .
Analysis of methylation data were conducted using the Genome Studio software version 2011.1 (Illumina Inc.) and methylation levels were expressed in beta values (β). β values that were calculated as the intensity of the methylated channel divided by the total intensity (β = Max (SignalB, 0)/(Max (SignalA, 0) + Max (SignalB, 0) + 100), and ranged from 0 (unmethylated) to 1 (fully methylated). β values with detection p-values > 0.01 were removed from analysis because were considered to fall below the minimum intensity and threshold. Also, probes that were localized to the sex chromosomes and those CpGs that contained single nucleotide polymorphisms were filtered out.
Differences in methylation levels between preoperative and postoperative groups (Δβ) was calculated for a given CpG site by subtracting the mean beta value from the pool of pre-surgery samples (pre-surgery period) as compared to the pool of samples collected after six months of RYGB and were tested by t test. Values of p were adjusted for multiple comparisons by using the false discovery rate (FDR below 5% was considered statistically significant). Also, a threshold for the significant CpG sites based on Δβ with a minimum value of 5% (value greater than 0.05 or less than −0.05) and FDR < 0.05 was applied. We used R software (version 3.2.0) to perform these analyzes.
Gene expression array analysis and data processing. Total RNA was extracted and purified from whole blood using phenol-chloroform extraction method modified by Chomczynski & Sacchi 41 . Extracted RNA was stored at −80 °C until the next steps. After, RNA integrity number (RIN) was analyzed using Bioanalyzer (Agilent Technologies, Cedar Creek, TX, USA). Expression levels were measured using the HumanHT-12 v4 Expression BeadChip (Illumina Inc.) Expression data were visualized and analyzed using Genome Studio Software (Illumina ® ). Beadchips were also scanned with the Illumina HiScanSQ system and Genome Studio Software (Illumina ® ) was also used. For quality control and normalization, the samples with less than 6000 significantly detected probes (detection p-value < 0.01) were excluded and differentially expressed transcripts were identified with 95% confidence of no more than 1% false positive using an Mann-Whitney test. A cutoff of p ≤ 0.05 was used to detect differentially expressed probes. In addition, a cutoff of fold change ≥2.0 (symmetrical fold change ≥2.0 or ≤ −2.0) was used. p values obtained from permutations and fold change cutoff values were then used to minimize the chances of false positives. Gene expression microarray results were validated using real-time polymerase chain reaction (RT-PCR; Applied Biosystems Gene Expression Assays; Applied Biosystems, Foster City, CA, USA) 20 .
For both methylation and expression analysis, samples from obese and normal weight women were randomly scattered on each Beadchip and not one batch was one and the other was opposite.
overlapping of pooled gene expression and DnA methylation analysis. Before comparing gene expression and methylation analysis, a correction by cellular type was made in DNA methylation data (Supplementary Table 1). By crossing and comparing the differentially regulated gene (DEG) and differentially methylated genes (DMG), the list generated from the gene expression analysis 21 were overlapped with the lists from the methylation analysis 20 using Venn diagram (http://bioinfogp.cnb.csic.es/tools/venny/). Figure 2 shows a schematic diagram of integrative analysis of methylation and expression arrays.