Altered intragenic DNA methylation of HOOK2 gene in adipose tissue from individuals with obesity and type 2 diabetes.

AIMS/HYPOTHESIS: Failure in glucose response to insulin is a common pathology associated with obesity. In this study, we analyzed the genome wide DNA methylation profile of visceral adipose tissue (VAT) samples in a population of individuals with obesity and assessed whether differential methylation profiles are associated with the presence of type 2 diabetes (T2D).
METHODS: More than 485,000 CpG genome sites from VAT samples from women with obesity undergoing gastric bypass (n = 18), and classified as suffering from type 2 diabetes (T2D) or not (no type 2 diabetes, NT2D), were analyzed using DNA methylation arrays.
RESULTS: We found significant differential methylation between T2D and NT2D samples in 24 CpGs that map with sixteen genes, one of which, HOOK2, demonstrated a significant correlation between differentially hypermethylated regions on the gene body and the presence of type 2 diabetes. This was validated by pyrosequencing in a population of 91 samples from both males and females with obesity. Furthermore, when these results were analyzed by gender, female T2D samples were found hypermethylated at the cg04657146-region and the cg 11738485-region of HOOK2 gene, whilst, interestingly, male samples were found hypomethylated in this latter region.
CONCLUSION: The differential methylation profile of the HOOK2 gene in individuals with T2D and obesity might be related to the attendant T2D, but further studies are required to identify the potential role of HOOK2 gene in T2D disease. The finding of gender differences in T2D methylation of HOOK2 also warrants further investigation.

Introduction

Obesity is frequently associated with the resistance of peripheral tissues (muscle and adipose) to the action of insulin. Because of a state of inflammation in adipose tissue infiltrated by inflammatory cells, secretions from these cells and fat derived adipokines promote the development of type 2 diabetes (T2D) [1-3]. Genome wide studies have led to the characterization of gene variants that are associated with an increased risk of T2D development, but their effect size is reduced when compared with traditional risk factors such as obesity, unhealthy diet or family history of diabetes [4].

Previous works have demonstrated alterations in methylation profiles of genes involved in the pathogenesis of obesity and T2D development [5]. For example, in subcutaneous adipose tissue from a group of 31 healthy men following exposure to exercise training [6], changes in methylation were found related to variations in mRNA expression in some of the T2D associated genes––HHEX (Hematopoietically-Expressed Homeobox Protein), IGF2BP2 (Insulin-like Growth Factor 2 Binding Protein 2), JAZF1 (JAZF Zinc Finger 1) and TCF7L2 (Transcription Factor 7-like 2). Ribel-Madsen et al conducted genome-wide methylation studies in insulin responsive tissues (skeletal muscle and subcutaneous adipose tissue) in elderly monozygotic (MZ) twin pairs discordant for T2D [7]. Methylation differences were found on the promoter region of CDKN2A (Cyclin-Dependent Kinase Inhibitor 2A) and HNF4A (Hepatocyte Nuclear Factor 4 Alpha) genes, although greater differences were found in genome wide repetitive DNA sequences such as LINE-1. These results were complemented by Nilsson et al, who described 1410 differentially methylated CpG sites between MZ T2D discordant twins, including KCNQ1 (Potassium Channel Voltage Gated KQT-like Subfamily Q, Member 1), NOTCH2 (Neurogenic Locus Notch Homolog Protein 2), TCF7L2 (Transcription Factor 7-like 2) and THADA (Thyroid Adenoma Associated) T2D related genes [8]. Recently, work by Chen et al has demonstrated that promoter hypermethylation of NR4A1 (Nuclear Receptor Subfamily 4 Group A Member 1) was elevated in T2D human samples as well as in T2D murine models and that this hypermethylation was associated with a decrease in mRNA levels. They also found that lack of NR4A1 (Nuclear Receptor Subfamily 4 Group A Member 1) protein is related to an increment in DNMT1 expression and the blocking of insulin signaling in patients with T2D [9]. These results were supported by the genome-wide DNA methylation analysis undertaken by Volkov et al in subcutaneous adipose tissue samples [10], whose work demonstrated an association between genetic and epigenetic mechanisms in terms of an observed relationship between single nucleotide polymorphisms and methylation at CpG sites of genes involved in metabolic patterns associated with diabetes development.

The majority of T2D cases occur in the context of a metabolic syndrome leading to a chronic excess of energetic substrates and ectopic fat storage, insulin resistance, elevated levels of inflammatory cytokines and, ultimately, a decrease in insulin secretion and the apoptosis of pancreatic β cells. The studies referred to earlier, have demonstrated the DNA methylation of T2D susceptible genes in different tissues (blood, skeletal muscle and fat) which are frequently associated with T2D development. Most of these works were carried out on subcutaneous fat, although VAT, along with skeletal muscle and liver, is in fact the most suitable candidate for determining the involvement of gene methylation in T2D development. Visceral fat not only acts as an energy storage tissue, but also as an endocrine organ, which releases hormones and adipokines, that contribute to expand macrophage population in VAT, which also liberate inflammatory cytokines. As the fat depot increases, levels of these molecules also rise and since they modulate the action of insulin in muscle and liver insulin insensitivity may result. VAT tissue has high levels of lipogenesis and lipolysis activity. This can result in hyperlipidemia and glucose intolerance [11]. Moreover, excess free fatty acids lead to ectopic lipid accumulation and lipotoxicity in muscle, where they inhibit glucose uptake, causing insulin resistance in these tissues [12, 13]

Prior research has described genome-wide DNA methylation patterns and site-specific differences in CpG methylation of candidate genes associated with T2D in VAT [14-18]. The aim of the present study is to contribute to the description of the VAT methylome in humans, and to explore the impact of epigenetics in diabetes development.

Materials and methods

Ethics statement

The study protocol (n° 68/10) was approved by the Ethical Committee at the Hospital Universitario Central de Asturias (Asturias, Spain) and all participants provided written and oral informed consent. The study was conducted in accordance with the principles of the Helsinki Declaration for human research.

Human tissue samples

Discovery cohort: a cohort of 18 (8 T2D and 10 NT2D) visceral adipose tissue (VAT) samples was obtained from patients (all female) from the Surgery and Endocrinology Departments at the Hospital Universitario Central de Asturias who underwent gastric bypass. Criteria for bariatric bypass surgery were BMI ≥35 and at least one obesity-related comorbidity (diabetes (n = 8), cardiovascular disease, CVD, (n = 9), non-alcoholic fatty liver disease (n = 4) or hypertension (n = 12)). The anthropometric and clinical characteristics of the samples are presented in .

Biological validation cohort: a cohort of 91 (55 NT2D and 36 T2D) VAT samples was obtained from patients (male and female) from the Surgery and Endocrinology Departments at the Hospital Universitario Central de Asturias (Asturias, Spain) who underwent gastric bypass. Criteria for bariatric bypass surgery were the same as the discovery cohort (data for obesity-related comorbidities: diabetes, n = 36; CVD, n = 4; non-alcoholic fatty liver disease, n = 49; and hypertension, n = 49). The anthropometric and clinical characteristics of the samples are presented in .

DNA extraction and bisulfite pyrosequencing methylation analysis

Genomic DNA was isolated from VAT samples using a phenol/chloroform protocol. DNA quality control and storage was carried out as indicated in Supplementary Methods.

Bisulfite modification and bisulfite pyrosequencing were performed as indicated in Supplementary Methods Specific pyrosequencing primers were designed using Pyromark Assay Design 2.0. Software and are described in .

Illumina 450K data methylation analysis

Methylation profiling of samples was carried out using Illumina Infinium HumanMethylation450 BeadChip Kit (Illumina Inc., USA) [19]. Data methylation analyses were performed as indicated in S1 Supplementary Methods, .

Genomic region analysis

The probes in the microarray were assigned to a genomic region according to their position relative to the transcript information extracted from the R/Bioconductor package TxDb.Hsapiens.UCSC.hg19.knownGene [20]. Genomic region analyses were conducted as indicated in Supplementary Methods.

CGI status analysis and Gene Ontology analysis

CpG island status was evaluated with respect to R. Irizarry's set of CpG island regions, included in the R/Bioconductor FDb.InfiniumMethylation.hg19 package [21, 22] and as indicated in Supplementary Methods.

To annotate the differentially methylated CpGs (dmCpGs) we used the HOMER annotation tool (http://homer.ucsd.edu/homer/ngs/annotation.html). All analyses were performed with a background reference comprising the whole set of genes present in the Illumina 450k platform. Results from HOMER (Hypergeometric Optimization of Motif EnRichment) were adjusted for multiple hypotheses using the Benjamini-Hochberg method for controlling the False Discovery Rate, FDR (threshold 0.05).

Statistical analysis

In order to identify the dmCpGs, a robust moderated t-test (R/Bioconductor package limma, version 3.24.15) was used on every probe in the array. Resulting p-values were adjusted for multiple comparisons by controlling the FDR (False Discovery Rate) using the Benjamini-Hochberg method. The difference between group means was used as a measure of effect size. In addition, the Genomic Region and CpG island status analyses were performed using a standard x2-test, and effect size was assessed by the Odds Ratio (OR). The significance threshold was fixed at 0.05 for all analyses. All work was carried out using R/Bioconductor version 3.1. A more detailed description of each of the individual types of analysis is included in Supplementary Methods.

Results

Genome wide DNA adipocyte methylation profiling and differentially methylated genes in individuals with T2D and obesity

The DNA methylation status of 485,000 methylation positions was analyzed in VAT from 8 T2D, 10 NT2D (all female) samples using Illumina Infinium Methylation Arrays. Data from probes were used to calculate a β-value between 0 and 1 (equivalent to 0%-100% methylation, respectively) (Fig 1A).

Fig 1

A. Hierarchical clustering heatmap. Showing differentially methylated CpGs of autosomal probes from women with obesity and discordant for type 2 diabetes (8 T2D and 10 NT2D). The heat map scale shows the range of methylation values, from 0 (blue) to 1 (yellow). Whether the CpG site analyzed is associated with a CpG island (CGI) or not can clearly be distinguished. Red asterisks mark the differentially methylated probes validated by pyrosequencing on HOOK2 gene. B. Strip charts showing β values of three differentially methylated CpGs (dmCpGs) located on HOOK2 gene in individual samples of the discovery cohort. The bold dotted line with the rhombus indicates the median value of each group of samples. C. Distribution of differentially methylated CpGs (dmCpGs) relative to CGIs, and relative distribution of dmCpGs across different genomic regions. Abbreviations: T2D (Type 2 Diabetes); NT2D (No Type 2 Diabetes).

We then compared the methylation (β value) of individual CpG sites in these two groups of samples. Using the robust and moderated Empirical Bayes methodology described in the Methods section, we obtained a set of 24 differentially methylated probes (adjusted p<0.05) (Data Accession number: https://doi.org/10.5281/zenodo.841654). Of these, 10 (related to seven genes) had a higher mean methylation value in T2D samples when compared to NT2D, Table 3. When we examined the CpG sites which were differentially methylated between T2D and NT2D samples, we found three positions (Illumina probes cg 11738485, cg 04657146 and cg 06417478) located in the intragenic region of chromosome 19, all of which overlay the same gene, HOOK2 (Hook Microtubule Tethering Protein 2) (Δβ ≥ 0.195) (Fig 1B). In addition, 14 CpG sites (9 genes) were hypomethylated in T2D compared with NT2D samples, Table 3.

Table 3

Differentially methylated CpG sites in T2D compared with NT2D samples.

Ilumina probe ID	CHR	probeStart	probeEnd	CpG Location	GENE SYMBOL (1)	βT2D	SD	βNT2D	SD	Δβ (2)	Adjusted pValue (3)
cg13978347	chr9	120140243	120140292	Non-CGI	ASTN2	0.8316	0.2690	0.7134	0.2072	0.1181	6.8418e-08
cg11738485	chr19	12877000	12877049	CGI	HOOK2	0.4639	0.2155	0.1676	0.2714	0.2963	1.1924e-07
cg04657146	chr19	12876947	12876996	CGI	HOOK2	0.3638	0.1407	0.1681	0.1967	0.1956	4.5798e-06
cg06417478	chr19	12876798	12876847	CGI-Shore	HOOK2	0.4211	0.1885	0.1674	0.2297	0.2537	2.6928e-06
cg01528832	chr10	65225240	65225289	CGI	JMJD1C;JMJD1C-AS1	0.1582	0.0411	0.0807	0.0436	0.0775	7.8258e-07
cg14456004	chr13	21872349	21872398	CGI	MIPEPP3	0.3598	0.0705	0.1405	0.0984	0.2192	4.0028e-06
cg07234876	chr8	600039	600088	CGI	NA	0.9501	0.0658	0.8173	0.1312	0.1327	1.0884e-05
cg20991723	chr1	152506874	152506923	Non-CGI	NA	0.9133	0.0729	0.8294	0.0399	0.0838	1.9410e-06
cg01475110	chr12	90332935	90332984	Non-CGI	NA	0.9408	0.0369	0.9071	0.0182	0.0337	4.2191e-05
cg23231268	chr3	46792462	46792511	CGI	PRSS50;PRSS45	0.9744	0.0247	0.9442	0.0478	0.0302	9.4175e-06
cg00777636	chr1	6446216	6446265	Non-CGI	ACOT7	0.8602	0.0631	0.9418	0.0242	-0.0815	1.2085e-07
cg04850148	chr17	34539744	34539793	Non-CGI	CCL4L1	0.3942	0.2744	0.6513	0.1171	-0.2570	3.4374e-05
cg01353608	chr4	15656863	15656912	CGI	FBXL5	0.0343	0.0456	0.0783	0.1052	-0.0439	1.3711e-06
cg03635532	chr19	40376835	40376884	CGI-Shore	FCGBP	0.4121	0.4372	0.6730	0.3553	-0.2609	3.6600e-06
cg05073382	chr8	2045798	2045847	CGI-Shore	MYOM2	0.5436	0.3105	0.8185	0.0330	-0.2748	3.6176e-06
cg11066601	chr1	185373438	185373487	Non-CGI	NA	0.4611	0.1220	0.6896	0.1635	-0.2285	3.6065e-08
cg06675417	chr18	77292443	77292492	CGI	NA	0.3087	0.2560	0.6583	0.1063	-0.3495	1.3992e-06
cg03070989	chr19	34311434	34311483	CGI	NA	0.1958	0.0963	0.3195	0.1256	-0.1237	8.1823e-06
cg14168080	chr7	157504135	157504184	CGI	PTPRN2	0.3583	0.4012	0.6672	0.2257	-0.3088	4.2091e-06
cg11757124	chr7	157526947	157526996	CGI-Shore	PTPRN2	0.6782	0.2053	0.8522	0.0637	-0.1740	2.8778e-06
cg17729891	chr19	51107510	51107559	CGI	SNAR-F	0.1715	0.1862	0.2931	0.1452	-0.1216	9.5658e-08
cg19825302	chr19	51107512	51107561	CGI	SNAR-F	0.1780	0.2096	0.2911	0.1554	-0.1130	2.0739e-08
cg23284931	chr11	13983273	13983322	CGI-Shore	SPON1	0.7280	0.0664	0.8348	0.0394	-0.1067	1.3765e-08
cg12587985	chr7	64295583	64295632	CGI	ZNF138	0.7722	0.0613	0.8887	0.0220	-0.1165	2.6298e-07

T2D: Type 2 Diabetes; NT2D: No Type 2 Diabetes; SD: Standard Deviation.

(1) Gene symbol “NA” represents an intergenic region to which no known genes map.

(2) A positive value for Δbeta indicates hypermethylation in T2D in comparison to NT2D samples, and a negative value indicates hypomethylation.

(3) Resulting p-values were adjusted for multiple comparisons by controlling the FDR using the Benjamini-Hochberg method.

Genomic distribution of differentially methylated CpGs in VAT samples

We compared the genomic distribution of the differentially methylated CpG sites with all the sites analyzed by the Infinium array based on the genomic location with respect to CpG islands. We did not find any significant differences, but there was an overrepresentation of CpG islands in the hypermethylated CpG sites in T2D samples with respect to NT2D samples (Pearson´s chi-squared test, adjusted p = 0.48; OR = 1.97). The same result was obtained with the hypomethylated CpGs in T2D, which are mainly located in CGI regions (CpG island, CpG shore (2 kb flanking the island), CpG shelf (2 kb flanking the shore)) (Pearson´s chi-squared test, adjusted p = 0.69; OR = 1.31) (Fig 1C).

The overall distribution of hypermethylated CpGs in relation to the nearest gene (5´UTR, 1st exon, gene body, 3´UTR) or intergenic regions in our analysis was suggestive of enrichment for the gene body (exon) in T2D compared to NT2D samples (Pearson´s chi-squared test, adjusted p = 0.08; OR = 4.69). In addition, hypomethylated CpG sites were over-represented in intergenic regions for T2D, although this difference did not reach significance (Pearson´s chi-squared test, adjusted p = 0.72; OR = 3.12) (Fig 1C).

Biological significance of the differentially methylated genes in T2D samples

To identify the biological relevance of the differentially methylated genes between T2D and NT2D samples we performed Gene Ontology analysis. We did not find any molecular functions or biological processes that were significantly associated with the differentially methylated genes (Data not shown), although a significant enrichment in specific regions on chromosome 10 was observed, S1 Table.

Methylation status of HOOK2 gene in the biological validation cohort of obese samples

To confirm the results obtained from the DNA methylation array, we developed independent DNA methylation assays using bisulfite pyrosequencing. Among the genes with differential DNA methylation levels in VAT tissue from T2D individuals compared with NT2D samples we found differences in three CpGs located within the HOOK2 gene (Δβ = 0.296; Δβ = 0.195; Δβ = 0.253).

An initial validation was made on the discovery cohort. For all CpGs analyzed, significant correlations were found between the two methods used to determine the percentage of methylation in the three probes meeting our criteria to be considered differentially methylated between T2D and NT2D (cg 04657146 probe, adjusted p = 1.75 x 10−5,R2 = 0.69); cg 06417478 probe, adjusted p = 1.95 x 10−5,R2 = 0.59; cg 11738485 probe, adjusted p = 3.31 x 10−6,R2 = 075), S2 Fig.

To further explore the involvement of HOOK2 gene methylation in T2D development, bisulfite sequencing analyses were conducted in a biological validation cohort of obese samples (n = 91, 36 T2D and 55 NT2D). The degree of DNA methylation of HOOK2 in the three regions analyzed that was observed in this way included all the differentially methylated CpGs described in the discovery cohort.

The median and the interquartile range of each CpG analyzed by pyrosequencing are shown in (Fig 2). The results demonstrate the differential methylation status of the HOOK2 gene regions between T2D and NT2D samples. The median for the cg 04657146-region (with 7 CpG nucleotides) was significantly higher in the T2D group than the NT2D with respect to the discovery cohort (adjusted p = 0.022), as well as in the biological validation cohort (adjusted p = 0.021). Interestingly, when the validation cohort was separated by gender, the median for this region in females was significantly higher in T2D samples than the NT2D (adjusted p = 0.0001334), while there was no significant difference for males. In the cg 11738485-region (formed by 5 CpG nucleotides), the pattern of higher methylation in the T2D group was repeated, but only in the discovery cohort (adjusted p = 0.020), not in the validation cohort. However, significant differences did become apparent when the validation cohort was separated by gender, with the median in female samples being significantly higher in T2D samples than the NT2D (adjusted p = 0.05), while the male T2D samples were interestingly found to be significantly hypomethylated (adjusted p = 6.25 x 10−5) compared to male NT2D samples (Fig 2). In contrast, no significant results for the cg 06417478-region (formed by 2 CpG nucleotides) were obtained for either of the cohorts analyzed, even when separated by gender.

Fig 2

Box plots illustrate the methylation values of differentially methylated CpG regions between T2D and NT2D samples validated by pyrosequencing in HOOK2 gene.

Schematic representation of the target regions studied and the methylation values from T2D and NT2D samples in the discovery cohort and the validation cohort are shown. Vertical lines represent the location of each CpG site. The analyzed region (blue line) and the CpG sites in the array (red box) are highlighted. The p-values for the comparison of the groups are also indicated (**adjusted p-value <0.05; ***adjusted p-value <0.0001). Abbreviations: T2D (Type 2 Diabetes); NT2D (No Type 2 Diabetes).

Discussion

The prevalence of T2D is increasing worldwide and while the mechanisms responsible for increased T2D risk are still poorly understood, epigenetics is thought to be key to the process [23]. This work revealed significant changes in DNA methylation throughout the entire human genome of VAT in individuals with obesity and discordant for T2D. Significant differential methylation (BH adjusted p<0.05, absolute change in methylation 20%) of 10 CpG sites was observed between the two phenotypes.

In line with previously published studies [8, 24], our results for hypermethylated CpGs associated with T2D were suggestive of methylation enrichment in the gene body, though the difference did not reach significance, possibly due to the small size of the cohorts. Gene bodies are CpG-poor regions that contain multiple repetitive and transposable elements. Methylation in these regions, may function to silence repetitive DNA elements found within the gene body [24], or be responsible for changes at intron–exon boundaries, suggesting an association with the splicing that may create chimeric transcripts with dominant–negative or neomorphic effects, or through transcriptional interference and potential RNA interference (RNAi) effects if arising from antisense transcript originating at transposon promoters [25]. Recently, the work developed by Pheiffer et al [26] has postulated that the differential methylation observed in the intergenic regions of peripheral blood samples in T2D individuals is involved in microRNA regulation and in the pathogenesis of T2D development.

The list of differentially methylated genes in VAT of individuals with obesity and T2D compared with NT2D individuals in this study comprises six genes whose Δβ is greater than 0.20. These are MIPEPP3 (Mitochondrial Intermediate Peptidase Pseudogene 3), CCL4L1 (C-C Motif Chemokine Ligand 4 Like 1), FCGBP (Fc Fragment of IgG Binding Protein), MYOM2 (Myomesin 2), PTPRN2 (Protein Tyrosine Phosphatase, Receptor Type N2) and HOOK2. Among them, only CCL4L1 could have been found to be related with insulin resistance. This gene, a CC chemokine located at chromosome 17, belongs to the same family as CCL4, which has been associated with T2D. The two proteins differ at only three amino acids and have redundant function [27]. CCL4 and its receptor CCR5 play diverse roles in the inflammatory response underlying T2D, due to their chemoattractivity towards macrophages, natural killer cells, monocytes, and immature dendritic cells, all of which have widely been described as playing a role in T2D [28]. Increased serum levels of CCL4 produced by islet cells have also been reported in type 1 diabetes and pre-diabetic condition, as well as in patients with T2D, suggesting the involvement of CCL4 in various stages of this disease [29,30]. The role of CCL4L1 in T2D has not been described, but the hypomethylation observed in VAT of individuals with T2D, is thought to perhaps contribute to the inflammatory response associated with insulin resistance in adipose tissue. PTPRN2 has been identified as an autoantigen in type 1 diabetes [31]. As a transmembrane protein, PTPRN2 shuttles between secretory vesicles and the plasma membrane, and has been implicated in insulin exocytosis [32]. However, its precise role in the secretory pathway is unknown. Differential PTPRN2 methylation has been described in placental samples of newborns with intrauterine growth restricted (IUGR) conditions [33], a condition that has been linked to the development of T2D later in life [34].

Extending the list of genes linked to T2D, this work provides evidence to support the inclusion of the gene HOOK2. We have found aberrant intragenic DNA methylation of this gene in a population with obesity and T2D when compared with a group of people with obesity but without T2D. Differential methylation status was found in all the HOOK2 gene regions validated by pyrosequencing, and, while low sample number is a potential limitation, these differences were significant with respect to two of the three regions. However, to separate the effects of HOOK2 methylation and the numerous confounding factors, we will need much larger numbers of samples and more detailed clinical information. Only through meta-analyses combining genome-wide data sets generated in different laboratories with different cohorts will we be able to reveal a more complete picture. Another potential limitation of this work is the low methylation difference we have found (close to 20%). Although previous reports indicate that methylation differences over 10% in Illumina assays reach biological significance and have a low probability of being technical artifacts [35, 36], this remains theoretical, and the work here, unfortunately, does not demonstrate the functional significance of HOOK2 methylation in patients with T2D. Another issue is the lack of quantification in our study of other C modifications, such as 5hmC (5-hydroxymethylcytosine), that have been reported to be associated with demethylation processes [37], regulation of gene expression [38], and changes in the chromatin structure [39] which may also show differences between the populations analyzed and provide another piece of the puzzle of our knowledge of T2D.

The significant hypermethylation observed on the cg04657146-region in T2D female samples indicates it may be a key region of the HOOK2 gene for DNA modification potentially involved in attendant T2D in individuals with obesity, especially females. However, additional studies are required to identify the role of DNA methylation in the cg04657146-region with respect to transcription factor recognition, chromatin remodeling and HOOK2 gene expression. In addition, and strikingly, the cg 11738485-region was found hypermethylated in T2D female samples, but hypomethylated in T2D male samples. Similar differences were observed in the work developed by Hall et al [40] who, in a genome wide DNA methylation analysis conducted on human pancreatic islets, found sex-specific differences in the DNA methylation levels of several genes located in autosomal chromosomes. The fact that these sex-biased genes involved in the functioning of pancreatic islets and insulin pathways showed differential expression could be related to the phenotypic differences observed between males and females in relation to insulin response [41, 42].

Adipose tissue function and distribution can be regulated by sexually dimorphic genes, some of which are present in VAT and controlled by sex hormones [43, 44]. From childhood on, fat distribution is different in males and females, and the difference becomes more marked during adolescence and remains this way until the menopause, when the reduction in estrogens levels means that fat distribution in women acquires an android distribution. Possibly, exposure to female sex hormones within the context of a fat and glucose rich diet, which leads to an abnormal accumulation of VAT (central adiposity), could go some way to explaining the differences observed in the methylation the HOOK2 gene in men and women with T2D. However, larger study groups, preferably with longitudinal sample collection of HOOK2 gene methylation would be needed to confirm the sexual differences observed in the methylation of this gene and its contribution to the different strategies for glycemic maintenance adopted by males and females.

The role of HOOK2 in primary cilia assembly has been described [45]. Primary cilia project out from the cellular membrane in most vertebrate cells and exert a sensor function in relation to stimuli (light, chemical and mechanical), triggering signaling pathways that maintain cell homeostasis [46]. Defects in this structure lead to a group of disorders called ciliopathies, two of which, the Bardet-Biedl and Alström syndromes, count T2D among their clinical manifestations [47]. Given that T2D and obesity are common features in ciliopathies and that HOOK2 is involved in primary cilia structure; the HOOK2 hypermethylation observed in the T2D samples in this work, would suggest that the possible role indicated here for this protein in T2D development may be mediated through alterations to primary cilia.

Another possibility is that HOOK2 protein could be involved in GLUT4 (glucose transporter type 4) traffic. The HOOK family is a group of cytoplasmic linker proteins associated with microtubules that participate in microtubule-vesicle transport and organelle support. Microtubules play a central role in GLUT4 translocation and glucose uptake [48-51]. The role of this microtubule-associated protein in insulin secretion or GLUT-4 traffic is still unclear, but HOOK2 methylation observed in T2D samples in the current work could perhaps be altering its function within the adipocyte. Benton et al [52], analyzing genome wide DNA methylation before and after bariatric surgery, suggested a similar role for the MYO1C (Myosin IC) gene, which codifies for an actin-related protein involved in GLUT4 translocation.

In conclusion, the differential methylation profiling we have described in individuals with T2D and obesity, and, more specifically, the methylation differences observed on the HOOK2 gene, points to the possible contribution of epigenetic factors, along with others previously described, to T2D predisposition. However, supplementary studies are required to identify the exact role of HOOK2 protein in the development of this pathology, as well as in terms of its potential role as a biomarker of increased risk, or in the design of gender-specific antidiabetic treatments.

Experimental workflow of Illumina 450K data methylation analysis.

(TIF)

Click here for additional data file.

Correlations between Illumina 450 K array data and pyrosequencing analysis on HOOK2 gene in the discovery cohort.

Representative data for a single CpG site for the three probes analyzed are shown. Illumina probe IDs are indicated after gene name.

(TIF)

Click here for additional data file.

Gene Ontology of hypermethylated genes in individuals with T2D compared with NT2D.

(DOC)

Click here for additional data file.

Supplementary methods.

(DOC)

Click here for additional data file.

Table 1

The anthropometric and clinical characteristics of the discovery cohort and the validation cohort.

	Discovery Cohort			Biological Validation Cohort
Characteristics	NT2D	T2D	p-value	NT2D	T2D	p-value
N(male/female)	10 (0/10)	8(0/8)		55 (14/41)	36 (16/20)
Age (years)	48.2 ± 9.16	50.3 ± 7.68	0.4562	41.5 ± 10.1	46.47 ± 11.54	0.0358
BMI(kg/m2)	50.48 ± 7.82	44.91 ± 3.8	0.0661	47.06 ± 8.92	49.06 ± 5.84	0.4482
Fasting Plasma Glucose (mg/dl)	118.20 ± 37.95	213 ± 85.44	0.0144	118.17 ± 38.43	166.5 ± 66.03	0.0002
Total-Cholesterol (mg/dl)	171.8 ± 30.8	216.33 ± 33.15	0.1682	192.73 ± 31.89	181.5 ± 39.65	0.3594
TAG (mg/dl)	87.8 ± 27.55	182 ± 46.7	0.0861	129.68 ± 58.81	158 ± 63.34	0.1250

NT2D: No Type 2 Diabetes; T2D: Type 2 Diabetes; TAG: Triacylglycerides; BMI: Body Mass Index. Groups were compared using a t test. Bold data indicate statistical significance (p<0.05)

Table 2

Set of primers used to validate by pyrosequencing.

GENE	SEQUENCE	GenBankNUMBER	POSITION CGANALYZED
HOOK2-V46
BTN-FW	5´- AGGTGGTGGGTGATATTTATATT-3´	NC_000019.10 (12765633–12766634)	12766135
RW	5´-CCCTAACCCTACTCTCTCCT-3
SEQ	5´- CCACCTAAAAAATACCA-3´
HOOK2-V78
BTN -FW	5´-TGGGAGGTGGTGGGTGGGTGAT-3´	NC_000019.10 (12765532–12766533)	12766034
RW	5´-AAAAAACCAAAAATAAACACAAA-3´
SEQ	5´-CCAAAAATAAACACAAACACTA-3´
HOOK2-V85
FW	5´-GGGAGGTGGTGGGTGATA-3´	NC_000019.10 (12765686–12766687)	12766188
BTN-RW	5´-CCTAACCCTACTCTCTCCT-3´
SEQ	5´-GTTGGTATTTTTTAGGTGG-3´