Epigenetic changes in the metabolically healthy obese: A case-control versus a prospective study.

INTRODUCTION

Obesity has become a worldwide health problem, in which latest statistics reported an increased trend over the last years. Accordingly, long‐term obesity is involved in an increased risk of developing several metabolic diseases, such as cardiovascular diseases, insulin resistance or type 2 diabetes. Although it is common to find obesity‐related disorders in the majority of individuals with obesity, some of them do not show metabolic complications. Obese people without these metabolic disorders are defined as metabolically healthy obesity (MHO) individuals, whereas those who present metabolic complications are known as metabolically unhealthy obesity (MUO) subjects. Importantly, the prevalence of this healthy phenotype varies widely, in which a correct diagnosis is critical for a good prognosis. ,  Nevertheless, the variability of this prevalence may be due to the nature of related‐study designs and the studied populations or the confounding variables considered. Therefore, conducting appropriately designed studies to evaluate important clinical research questions is crucial, to elucidate the impact of MHO in the progression of detrimental obesity‐related phenotypes. In this case, epigenetic arises as a valuable tool for clinical applications. Epigenetic mechanisms might be implicated in the regulation of several metabolic disorders. However, until now, only a few studies have been conducted to understand the epigenetic changes in individuals with MHO and MUO. , ,

Therefore, we hypothesized that epigenetic changes may be involved in the development and the transition from MHO to MUO phenotype. Consequently, the aim of this study was to analyze the epigenome‐wide DNA methylation study in peripheral blood mononuclear cells (PMBCs) from participants with MHO and MUO, by using two study designs, such as case‐control and prospective approaches.

MATERIALS AND METHODS

Study design

The participants of the current study were selected from the Pizarra cohort study (Malaga, Spain). This study was a population‐based prospective study from 1051 participants started in 1995–1997, in which 547 subjects completed a follow‐up of 11 years. A total of 17 participants were selected from this prospective cohort. We performed two approaches: (a) case‐control analysis: we compared 9 participants with MHO (control) and 8 participants with MUO at 11 years of follow‐up; (b) prospective analysis: we included 8 participants with MHO (control) at the baseline who transitioned to MUO phenotype after 11 years of follow‐up. The MHO and MUO were classified following the National Cholesterol Program Adult Treatment Panel criteria (NCEP‐ATPIII) (Table 1). They were considered as MHO if had abdominal obesity (waist circumference >102 cm in men and >88 cm in women) and <2 of the NCEP ATPIII metabolic syndrome criteria was present: systolic blood pressure ≥135 mmHg or diastolic blood pressure ≥85 mmHg; fasting plasma glucose concentration ≥100 mg/dl; HDL‐C concentration <40 mg/dl in men and<50 mg/dl in women; fasting plasma TG concentration ≥150 mg/dl; or treatment with antihypertensive, lipid lowering, or glucose‐lowering medication. Serum samples were used for measuring fasting glucose, HDL‐c and triglycerides by oxidase (Accu‐Chek, Roche Diagnostics, Barcelona, Spain) and enzymatic methods respectively. All participants gave their informed consent. The study was reviewed and approved by the Ethics Committees of “Carlos Haya” Regional University Hospital (Málaga, Spain).

TABLE 1

Anthropometric and biochemical characteristics in the population study

Case‐control design
Prospective design
Variables	MHO (n = 8) baseline	MUO (n = 8) 11 years of follow‐up	MHO (n = 9) 11 years of follow‐up
Age	53.2 ± 8.9	62.42 ± 8.7	55.9 ± 11.2
Sex (Male/Female)	2/7		3/6
Fasting glucose (mg/dl)	107.9 ± 10.1	109.25 ± 16.6	90 ± 4.9b
BMI (kg/m2)	29.02 ± 4.3	31.16 ± 4.3a	29.9 ± 3.4
Triglycerides (mg/dl)	92.6 ± 37.8	100.3 ± 53.0	82 ± 25.2
HDL‐cholesterol	51.4 ± 9.5	53.5 ± 8.4	60.7 ± 6.3
DBP (mm Hg)	88.6 ± 16.8	90.5 ± 11.4	75.7 ± 9.3b
SBP (mm Hg)	138 ± 26	153.18 ± 23	126.3 ± 19.3b
HTA treatment (%)	22.2	55.6a	NAb

Data are expressed as the mean ± standard deviation, or as percentages. Statistical significat were calculated using the Kruskal‐Walli test, and for categoriacal data they were calculated using the chi‐square test or Fisher's exact test. Statistical significant were defined as p < 0.05 a) prospective study b) case‐control study. They were considered as MHO if had abdominal obesity (waist circumference >102 cm in men and >88 cm in women) and <2 of the NCEP ATPIII metabolic syndrome criteria was present: systolic blood pressure ≥135 mmHg or diastolic blood pressure ≥85 mmHg; fasting plasma glucose concentration ≥100 mg/dl; HDL‐C concentration <40 mg/dl in men and<50 mg/dl in women; fasting plasma TG concentration ≥150 mg/dl; or treatment with antihypertensive, lipid lowering, or glucose‐lowering medications. Abbreviations: BMI, body mass index; HDL cholesterol, high density liporptein cholesterol; DBP, diastolic blood pressure; SBP, systolic blood pressure; HTA treatment, arterial hypertension treatment.

DNA isolation, bisulfite reaction and epigenome‐wide DNA methylation analysis

DNA was extracted from peripheral blood using the QIAmp DNA Blood Mini Kit (Qiagen, Hilden, Germany). DNA concentration was quantified with a Qubit 3.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) using Qubit dsDNA HS Assay Kit Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). 500 ng of genomic DNA was bisulfite‐treated using Zymo Research Infinitum HD FFPE Methylation kit (Zymo Research Corp, Irvine, CA, USA) and purified using DNA‐Clean‐Up kit (Zymo Research Corp, Irvine, CA, USA).

Genome wide DNA methylation profiling

Over 850.000 CpG were used in the 850 k Infinium Methylation EPIC Bead Chip Kit (Illumina, San Diego, CA, USA) following the Infinium HD Assay Methylation instructions. Methylation analysis workflow was based on a previous work. Raw data extraction and downstream processing were performed with R software (version 4.0.0) using bioconductor packages. Minfi package was used to DNA quality checks, bisulfite modification, hybridization, data normalization, statistical filtering, and value calculation were performed as described on Maksimovic work (2017). The annotation of the CpGs were analysed using IlluminaHumanMethylationEPICmanifest (version 3.0) and IlluminaHumanMethylationEPICanno.ilm10b4.hg19 (version 6.0). Finally, probes with single nucleotide polymorphism (SNP) or single‐base extension sites and sex chromosomes were removed.

Methylation data analysis

The differentially methylated CpGs (DMPs) for case‐control and prospective analysis were obtained on the matrix of β‐values in limma and minfi package using M‐values (p < 0.05) For case‐control study, we adjusted data by age sex, fasting glucose and arterial hypertension. For prospective analysis, we adjusted data by age and sex. We obtained differentially methylated regions (DMRs) in both case‐control and prospective models, by using the dmrcate package and β‐values. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) were obtained using the gometh package. A network analysis was performed using the STRING analysis (version 11.00) (https://string‐db.org/). Clustering was conducted by using web‐platform cytoscape (https://cytoscape.org/).

Differences in anthoprometric and biochemical variables with the chi‐square test for categorical data and Kruskall–Wallis test for continous data. Statiscally significant values were set at p < 0.05.

RESULTS

Clinical characteristic of study participants

Principal clinical variables are shown in Table 1. In prospective design, body mass index (BMI) and arterial hypertension (HTA) treatment were significant higher in MUO at 11 years of follow‐up (p < 0.05). In case‐control at 11‐year follow‐up desing, fasting glucose, diastolic and systolic blood pressure values were significantly higher in MUO group (p < 0.05).

Differential Methylation Patterns in the case‐control and prospective studies

DNA methylation patterns from both case‐control and prospective models clearly distinguished two clusters (individuals with MHO and MUO) by the methylation levels, as is shown in the principal component analysis (PCA) (Figure 1). Additionally, by studying differential methylation analysis, we found significant DMPs around the genome, as Manhattan plot display (Figure 1B,E). Moreover, Volcano plot shows differentially hypo‐ and hypermethylated CpGs (Figure S1A,C), in both approaches (p < 0.05). Overall, in the case‐control approach, we observed that the average global DNA methylation was decreased in the MUO group (p < 0.001) (Figure 1C), whereas in the prospective design, individuals with MUO at 11‐years follow‐up showed increased average global DNA methylation in comparison with subjects at baseline (MHO) (p < 0.001).

FIGURE 1

Differential methylation patterns: (A) PCA plot of Case‐control analysis, (B) CpGs chromosome distribution in Manhattan plot in case‐control (C) Boxplot of global β‐methylation of significant DMPs case‐control; (D) PCA plot of prospective analysis (E) CpGs chromosome distribution in Manhattan plot in prospective analysis

DMPs from both case‐control and prospective analyses were mainly distributed in open sea regions, body genes and intergenic regions (Figure 2A,B). At the promoter region, the majority of DMPs were located at the transcription start site (TSS)1500, TSS200 and 5´UTR (untranslated region) (Figure 2B). The chromosomal distribution of DMPs are shown in the Figure 2C. The illustration of hypo‐ and hypermethylated probes according to the chromosomal location and the genomic regions are summarized in the Figure S2.

FIGURE 2

Localization of Differentially methylated positions (A) % of DMPs distribution according to CpG island position (OpenSea, Island, South‐shore, North‐shore, North‐Shelf, South‐Shelf) in Case‐control (Grey) and prospective (black) analysis. (B) % of DMPs distribution according to gene localization (TSS1500, TSS200, 3’‐UTR, 5’‐UTR, Exon Boundary, Body, 3’UTR and IGR) in Case‐control (grey) and prospective (black) approaches. (C) % of DMPs according to chromosome in Case‐control (grey) and prospective (black) analysis D) The most enriched DMPs in case control and in prospective Model (logFC>|30|, p‐value < 0.001)

Differential Methylation Positions and Regions in the case‐control and prospective analyses

The logistic regression analysis identified 46,035 DMPs in the case‐control model (25,768 hypermethylated probes and 20,267 hypomethylated probes), and 50,464 DMPs in the prospective model (18,092 hypermethylated probes and 32,372 hypomethylated probes) (p < 0.05). Table S1 and Table S2 show the most 20 differentially methylated probes. We found 2824 common DMPs between the case‐control and prospective models. We further filtered them using those DMPs located at the promoter regions and within island regions. We also selected those genes that contain ≥2 DMPs. Thus, we found 14 enriched genes, such as AC008806.2, AC018464.3, C5orf15, CPNE2, GNG12, HYAL1, JAK1, L3MBTL4, LZTS1, NKAIN1, OR10H1, RP11‐522N14.1, TAS2R7 and VWC2.

Additionally, to find additional enriched genes that are strongly associated with MUO phenotype, we selected those genes that contain more than 20 DMPs. A total of 19 genes (U6, Y_RNA, TBCD, PI16, snoU13, PRKCZ, SDK1, C7orf50, FOXP1, ATP11A, INPP5A, HDAC4, RPTOR, B3GNTL1, AKAP13, DAB1, DDR1 and RBFOX3) were found in the case‐control model, and 17 genes (U6, Y_RNA, snoU13, RP11‐73O6.4, 7SK, C7orf50, TBCD, RPTOR, PI16, PRDM8, RBFOX3, MSLN, FOXP1, PRKCZ, SDK1, SPRN and SHANK2) from the prospective model (Figure 3A). Interestingly, U6, TBCD, RPTOR, PI16 or RBFOX3 genes were common in both analysis.

FIGURE 3

(A) The most enriched gene in Case‐control and prospective approaches. Choice the commun DMPs from case‐control and prospective. Gene filtered by those located at promoter area and islands, with more than 2 DMPs. (B) String net

In addition, to identify potential epigenetically regulated genes in blood samples and with a strong association with the MUO phenotype, we first filtered DMPs from both case‐control and prospective models, by using a cut‐off value of p < 0.001 and fold change >|30%|. Interestingly, 9 DMPs were found in the case‐control analysis and only cg24147543 was selected from the prospective model (Figure 2D).

Finally, to identify highly correlated CpG regional clusters, we studied the DMRs. Our DMR analysis identified 2973 DMRs from the case control model and 2491 DMRs from the prospective model. Table S3 and Table S4 show the top 20 DMRs. Importantly, ISO2 and DPYS genes were commun in both models, being related to metabolic pathways. Furthermore, a correlation analysis was performed to test whether these DMRs were associated with the clinical variables of the study population. DMR3 and DMR4 (overlapping DPYS gene) were correlated with blood pressure and fasting glucose variables in both case‐control and prospective designs (Figure [Link], [Link]).

Functional enrichment analysis

To perform the functional analysis, we selected a total of 50 genes. First of all, we selected 14 genes which were in common between case‐control and prospective analysis, also they were located at promoter and islands area and they had more than 2 DMPs per gene. In addition, we looked for those genes in case‐control a prospective with more than 20 DMPs. Then, we obtained 19 genes in case‐control and 17 in the prospective study (Figure 3A). We used STRING to perform a network analysis with the 50 genes selected (Figure 3B). The network obtained showed two differentiated clusters. One of them related to insulin and mammalian target of rapamycin (mTOR) signalling pathways, whereas the other cluster was related to JAK‐STAT signalling pathway. These results were according to the GO and KEGG analyses (Figure [Link], [Link]). Interestingly, KEGG also provided important pathways, such as pancreatic and insulin secretion, as well as chemokine and cGMP‐PKG signaling pathways or carbohydrate digestion and absorption.

DISCUSSION

This study presented an epigenome‐wide DNA methylation analysis for the epigenetic changes in the MUO phenotype, using both the case‐control and prospective designs, with a follow‐up of 11 years. After analyzing both approaches, this study demonstrated that the DNA methylation landscape in blood samples between MHO and MUO phenotypes was significantly different. In this context, our findings offer several methylated candidate genes, which may be further considered in the epigenetic regulation processes in metabolic related‐pathways.

Previous studies evaluating the role of epigenetic in the transition from MHO to MUO are relatively scarce, despite the potential role of epigenetic in metabolic disorders. Accordingly, a case‐control study conducted by Desiderio and colleagues (2019) found that the methylation of the ANKRD26 gene was significantly increased in individuals with MUO when compared to subjects with MHO (p = 0.009). A prospective study carried out by Gallardo‐Escribano et al. (2020) reported that several metabolic genes (LPL, SCD, SREBF1, LXR) were epigenetically altered in the MHO group after 12 months of follow‐up. Recently, our group identified in a prospective study two differentially methylated positions located in the ZFPM2 and CYP2E1 genes, closely associated with the transition from MHO to MUO state. A recent meta‐analysis conducted by Andrade et al. (2021) found more than 2000 DMRs in the adipose tissue of individuals with MUO compared to individuals with MHO, as our study reported.  There are several genes that overlapped our DMRs, such as FLRT, FRMPD4 PSMD5, NKX or LRRC, suggesting a similar effect in obesity. , However, all of these studies did not provide potential candidates to better predict the pathways implicated in the transition from MHO to MUO.

In our study, we found a list of potential candidates for the transition from MHO to MUO phenotype. These genes were strongly related to insulin and mTOR signalling pathways, suggesting that these pathways are epigenetically associated with the progression of metabolic disorders. Interestingly, we found potential DMRs overlapped genes, in which two genes were common in both approaches, such as ISCO2 and DPYS. In this way, a study reported that ISOC2 was able to inhibit p16 (INK4a) in a dose‐dependent manner, in which p16 may enhance glucose uptake and mitochondrial activity. Similarly, DPYS gene provides instructions for making an enzyme called dihydropyrimidinase. This enzyme leads to the production of beta‐aminoisobutyric acid as metabolite of enzymatic reaction, in which beta‐aminoisobutyric acid increases the production and release of leptin, an important hormone in obesity and insulin secretion.

In spite of the sample size limitation, the selected cohort analyses was very homogeneous. Moreover, by using both approaches provide enough data to avoid possible biases in the analysis, to establish strong relationship between epigenetic and MUO phenotype. In addition, we also adjusted our analyses by potential confounding variables, such as glucose, sex or age, which can strongly influence the DNA methylation pattern of the participants.

CONCLUSION

In conclusion, our study displayed the potential of DNA methylation profile to differentiate between MHO and MUO phenotypes. We found potential candidate genes and significant related‐pathways, which suggest that DNA methylation changes may have a strong effect in the regulation of insulin metabolism. However, future studies are needed to clarify the implication of these pathways in the transition from MHO to MUO phenotype.

CONFLICT OF INTEREST

The authors have declared no conflict of interest.

AUTHORS’ CONTRIBUTION

FJT and SM designed research; TMLP and HB conducted research; FJT, GRM, and SV were the coordinators of subject recruitment at the outpatient clinics; TMLP, EAL and HB analyzed the data; TMLP, HB, and SM interpreted statistical analysis and data; HB, TMLP, SM and MMG drafted the paper; FJT and SM supervised the study. All authors revised the manuscript for important intellectual content, read and approved the final manuscript.

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Fig S3

Click here for additional data file.

Fig S4

Click here for additional data file.