Methylation changes in the peripheral blood of Filipinos with type 2 diabetes suggest spurious transcription initiation at TXNIP

Abstract

While much work has been done in associating differentially methylated positions (DMPs) to type 2 diabetes (T2D) across different populations, not much attention has been placed on identifying its possible functional consequences. We explored methylation changes in the peripheral blood of Filipinos with T2D and identified 177 associated DMPs. Most of these DMPs were associated with genes involved in metabolism, inflammation and the cell cycle. Three of these DMPs map to the TXNIP gene body, replicating previous findings from epigenome-wide association studies (EWAS) of T2D. The TXNIP downmethylation coincided with increased transcription at the 3′ UTR, H3K36me3 histone markings and Sp1 binding, suggesting spurious transcription initiation at the TXNIP 3′ UTR as a functional consequence of T2D methylation changes. We also explored potential epigenetic determinants to increased incidence of T2D in Filipino immigrants in the USA and found three DMPs associated with the interaction of T2D and immigration. Two of these DMPs were located near MAP2K7 and PRMT1, which may point towards dysregulated stress response and inflammation as a contributing factor to T2D among Filipino immigrants.

Introduction

Type 2 diabetes (T2D) is a complex disorder characterized by chronic hyperglycemia arising from impaired insulin secretion, elevated glucagon secretion and acquired insulin resistance. Chronic hyperglycemia results in persistent oxidative stress, leading to kidney failure, retinopathy and cardiovascular complications (1, 2). T2D accounted for approximately 90% of 500 million adults living with diabetes worldwide in 2017, a number which is expected to increase based on rising prevalence in middle- to low-income countries (3, 4). In the Philippines, the latest 2013 National Nutrition Survey revealed that hyperglycemia had a prevalence of 5.6% and is projected to rise based on rapidly progressing urbanization and increasingly sedentary lifestyles (5, 6). Meanwhile, T2D prevalence is much higher among Filipino immigrants living in the USA and is even higher in contrast to a Caucasian reference (7). This has led T2D and its associated complications to become the sixth leading cause of death among Filipinos in 2016 (8).

Table 1

Characteristics of the population under study

	Control2 n = 218		Incident T2D3 n = 147
	PH n = 98	US n = 120	PH n = 100	US n = 47
Age1	48.4 ± 15.3 (18.0, 83.0)	38.6 ± 16.5 (18.0, 83.0)	61.3 ± 11.3 (33.0, 83.0)	60.3 ± 11.0 (21.0, 85.0)
Gender (%)1
Male	31 (31.6)	47 (39.2)	22 (22.0)	20 (42.6)
Female	67 (68.4)	73 (60.8)	78 (78.0)	27 (57.4)
Body mass index1	23.6 ± 4.47 (14.1, 42.0)	25.5 ± 3.66 (17.4, 38.7)	25.9 ± 5.01 (14.4, 38.7)	28.1 ± 4.20 (20.0, 38.7)
Waist circumference1	82.0 ± 11.0 (54.0, 125.2)	83.6 ± 10.7 (58.0, 114.0)	90.6 ± 12.4 (54.0, 115.0)	95.8 ± 10.1 (81.0, 116.0)

¹Reported as mean ± SD (range) or number (percentage).

²Diagnosis was made following the 2004 ADA recommendation: FBS < 100 mg/dl and OGTT < 140 mg/dl, with no history of diabetes or prediabetes.

³Diagnosis was made following the 2004 ADA recommendation: FBS ≥ 126 mg/dl, OGTT ≥200 mg/dl or a previous diagnosis of type 2 diabetes.

Genome-wide association studies (GWAS) show that genetic variants explain only a small portion (~20%) of the risk and heritability of T2D (9–11), highlighting the role of epigenetic, lifestyle and environmental risk factors. Epigenetics, in particular DNA methylation, serves as a possible molecular link for the interaction between the environment and T2D. DNA methylation involves alteration of CpG dinucleotides by the covalent transfer of a methyl group to the 5′ position of the cytosine ring, forming 5-methyl-cytosine (5mC) (12, 13). Epigenetics may contribute to the pathogenesis of complex metabolic diseases, including T2D, through the alteration of gene expression programs that predispose individuals to the diabetic phenotype (14–16). Also, altered homeostasis in T2D leads to epigenetic changes associated with the development of disease complications (15).

The development of methylation arrays has led to several studies identifying differential methylation associations to T2D in epigenome-wide association studies (EWAS) (17–22), but little has been done towards the functional characterization of these associations. In this study, we examined the peripheral blood methylation profiles of Filipinos to determine epigenetic changes associated to T2D. We also explored the contribution of methylation changes towards the higher incidence of T2D in Filipino immigrants living in the USA through comparison of methylation profiles of immigrants with Filipinos living in the Philippines. Of note, we identify diabetes-specific downmethylation at the TXNIP gene body and combined this with published RNA-seq and epigenomic data sets to provide a more detailed view of possible functional dysregulation in peripheral blood arising from T2D-associated methylation changes.

Results

Description of the study participants

A total of 365 patients were included in the study. Of those, 218 were healthy controls and 147 were cases (Table 1). Participants from USA had a lower mean age (44.7 ± 18.0) than those from the Philippines (54.9 ± 14.9; P = 9.84 × 10⁻²⁹). Most of the subjects were female (67.1%), a trend which persisted after stratification into subjects with diabetes and controls. Anthropometric measurements for obesity were higher for subjects with T2D, such as body mass index (BMI) and waist circumference (WC), which persisted after subjects were stratified by location (Supplementary Material, Table S2).

Quality control of methylation samples

Of the 365 array-processed patient blood samples, 346 (94.8% of total samples) passed quality control, with 15 (4.12%) having conflicting sex assignment and 4 (1.10%) having ambiguous sex assignment (Supplementary Material, Fig. S1A). We removed 38 365 probes (4.43% of total probes) which had a detection P value surpass the 10⁻¹⁶ threshold recommended by Lehne et al. (Supplementary Material, Fig. S1B) (23). An additional 111 574 (12.9%) were removed based on literature recommendations to mask probes containing SNPs within 5 bp of the 3′ end, mapping to the sex and mitochondrial chromosomes and exhibiting off-target specificity (Supplementary Material, Fig. S1C) (24). We were left with 715 920 (82.7%) probes for profile analysis.

Differentially methylated positions associated with T2D

We find 177 CpGs with significantly different methylation intensities corresponding to differentially methylated positions (DMPs) in patients with diabetes in contrast to healthy controls (Table 2, Supplementary Material, Table S3). The top associated DMP to T2D was downmethylation at cg19693031 (P = 4.17 × 10⁻¹⁰), located at the 3′ untranslated region (3′ UTR) of TXNIP (Table 3), which is the top T2D association in previous EWAS of blood (17–20). Also, we find an additional two DMPs located in the TXNIP gene body (Supplementary Material, Fig. S2B–D). In agreement with previous studies, we find downmethylation at cg00574958 (P = 3.92 × 10⁻³) located at the first intron of CPT1A (25). Effect sizes of the T2D-associated DMPs were small; the root mean square (RMS) of the log-transformed fold change in terms of M values was 2.80 × 10⁻¹ (Supplementary Material, Fig. S2A). Using the more intuitive β value, this results in an RMS fold change of 7.76 × 10⁻², corresponding to an average 1.05-fold increase or 0.947-fold decrease of β-methylation intensity.

Table 2

Summary of associations found in this study

	Diabetes	Diabetes* immigration1	Fasting blood glucose2	Glycated hemoglobin2	HOMA-IR2	Waist circumference
Upmethylated/direct relationship	90	0	157	470	297	93
No change/relationship	715 743	715 917	715 617	714 641	715 202	715 750
Downmethylated/indirect relationship	87	3	141	429	344	94
Sample sizes (n)
Control	204	204	204	204	204	204
Diabetic	142	142	140	45	97	142

¹Interaction experiment between T2D and US immigration.

²Some samples had missing FBG, Hb-A1C and HOMA-IR information.

Table 3

Top five methylation changes associated to type 2 diabetes

	Position1	CGI position	Gene2	Context	Effect size3	P value4
cg19693031	1:145993516	North Shelf	TXNIP	3′ UTR	−0.667	4.17 × 10−10
cg16125917	19:7903720	Island	MAP2K7	TSS-200	0.500	1.33 × 10−8
cg25544073	8: 97644077	Island	MTDH	TSS-200	0.413	3.38 × 10−5
cg04632724	15:101652191	Island	TM2D3	Gene body	−0.626	3.38 × 10−5
cg26974062	1:145994624	North Shore	TXNIP	Gene body	−0.510	3.54 × 10−4

¹Position on the GRCh38/hg38 genome. Stated as Chromosome:Base.

²Protein coding genes only.

³Reported as fold change of log₂-transformed M values.

⁴Adjusted using the Benjamini–Hochberg procedure.

DMPs associated with the interaction of T2D and immigration

To identify possible epigenetic contributions towards increased T2D incidence among Filipino immigrants, we determined the interaction of immigration with T2D and find three associated DMPs (Table 2): downmethylation at cg16125917 (P = 3.31 × 10⁻²) located at the MAP2K7 proximal promoter, cg19309676 (P = 4.20 × 10⁻²) located at the ADM5 proximal promoter and PRMT1 gene body and cg21633426 (P = 2.58 × 10⁻²) which is in proximity to two lincRNA genes (AC023115.2 and LINC01829) (Table 4). Only cg16125917 and cg19309676 are found to be also differentially methylated in T2D.

Table 4

Methylation changes associated to the interaction of type 2 diabetes and US immigration

	Position1	CGI position	Gene2	Context	Effect size3	P value4
cg21633426	2: 67213358	Open Sea	AC023115.2, LINC01829		−0.420	2.58 × 10−2
cg16125917	19:7903720	Island	MAP2K7	TSS-200	−0.453	3.31 × 10−2
cg19309676	19: 49688183	Island	PRMT1, ADM5	TSS-1500, gene body	−0.318	4.20 × 10−2

¹Position on the GRCh38/hg38 genome. Stated as Chromosome:Base.

²Protein-coding genes only.

³Reported as fold change of log₂-transformed M values.

³Adjusted using the Benjamini–Hochberg procedure.

Gene groups and pathways affected by T2D-associated methylation changes

We find 153 DMPs (86.4% of DMPs) which co-localize to 194 unique genes, either at promoters or at gene bodies. Most genes were affected with only one DMP (190, 97.9% of genes) (Fig. 1A) and were protein-coding (164, 84.5%) (Fig. 1B). Of the non-coding genes (30, 15.5%), the most prominent were antisense transcripts (13, 6.70%) and long intergenic non-coding RNAs (lincRNAs; 8, 4.12%) (Fig. 1C).

Figure 1

Genomic contexts of T2D-associated DMPs. (A) Although gene regions are 609 usually associated with several CpGs, especially with CGIs at promoter elements, most of the associated genes have only one DMP. (B) Most associated genes are protein coding, (C) while the rest are either antisense transcripts, lincRNA or pseudogenes. (D) Ontologies of gene affected by differential methylation. (E) Representative pathways affected by differentially methylated genes.

Overrepresentation analysis of gene ontology (GO) assignments by hypergeometric test shows enrichment of DMPs located at genes involved in metabolism (Fig. 1D, Supplementary Material, Table S4). We observe an enrichment of biological process (BP) terms related to metabolism, particularly of phosphorous (GO:0006793, P = 1.47 × 10⁻³; GO:0006796, P = 1.47 × 10⁻³), nitrogen compounds (GO:0006163, P = 1.28 × 10⁻²; GO:1901564, P = 4.55 × 10⁻²) and the biosynthesis of fatty acids (GO:0035338, P = 4.55 × 10⁻²; GO:0071616, P = 4.55 × 10⁻²). BP terms related to protein modification (GO:0006468, P = 1.28 × 10⁻²; GO:0036211, P = 4.55 × 10⁻²) were also enriched. Concurrent with the enrichment of metabolism-related BP terms, cellular component (CC) terms reveal that DMPs targeted genes co-localizing in the mitochondria (GO:0005739, P = 1.14 × 10⁻⁴), mitochondrial envelope (GO:0005740, P = 4.00 × 10⁻³) and cytochrome complex (GO:0070069, P = 2.83 × 10⁻²), while molecular function (MF) terms reveal enrichment of genes with kinase (GO:0004672, P = 5.11 × 10⁻³), phosphotransferase (GO:0016772, P = 7.83 × 10⁻³) and fatty acid ligase activity (GO:0015645, P = 9.37 × 10⁻³).

While GO assignments excel at gene classification, they fall short on determining implications from dysregulated gene interaction networks. We identify 134 enriched pathways containing genes dysregulated by T2D-associated DMPs (Supplementary Material, Table S5), which we simplify by hierarchical clustering to obtain 23 representatives (Fig. 1E). The top three representative pathways are cell cycle (hsa04110, P = 2.63 × 10⁻⁹), response to bacterial invasion of epithelial cells (hsa05100, P = 7.98 × 10⁻⁹) and mitogen-activated protein kinase (MAPK) signaling (hsa04010, P = 5.48 × 10⁻⁸). Pathways related to the inflammatory response were also found enriched, such as Toll-like receptor (TLR; hsa04620, P = 1.20 × 10⁻⁵), T cell receptor (TCR; hsa04660, 3.93 × 10⁻⁵), NFκB (hsa04064, P = 8.58 × 10⁻⁵) and tumor-necrosis factor (TNF; hsa04668, P = 1.34 × 10⁻⁴) signaling. We also find the pathways related to T2D pathology, such as the insulin resistance (hsa04931, P = 6.31 × 10⁻³), insulin signaling (hsa04910, P = 3.28 × 10⁻⁴) and glucagon signaling (hsa04922, P = 2.55 × 10⁻⁷).

As genes may work across multiple pathways, we took a gene-centric approach by counting the occurrence of each gene in all enriched pathways. We identify 52 unique genes with DMPs shared across the enriched pathways (Supplementary Material, Fig. S3A). ITPR1 (upmethylated at cg00448292, P = 4.01 × 10⁻²) had the highest occurrence, found in 39 pathways (20.5% of pathways). The next two were CAMK2B (upmethylated at cg08106392, P = 2.27 × 10⁻²) and CREBBP (upmethylated at cg06365577, P = 3.22 × 10⁻²), which occurred in 28 (14.7%) and 26 (13.7%) pathways, respectively. CPT1A occurs in 6 (3.16%) pathways, but TXNIP is not found in any. We examined pathways related to inflammation and find that these pathways share differentially methylated genes with the MAPK pathway (Supplementary Material, Fig. S3B).

Methylation changes associated with hyperglycemia and insulin resistance

T2D is characterized by increased insulin resistance and decreased β-cell function leading to higher blood glucose levels, so we compared associations to measures of hyperglycemia (fasting blood glucose, FBG; and glycated hemoglobin, Hb-A1C) and insulin resistance (homeostatic model assessment of insulin resistance, HOMA-IR) against T2D associations (Table 2, Supplementary Material, Table S6). As T2D incidence is also correlated with obesity, we compared WC associations as well. There is little overlap between hyperglycemia-, insulin resistance-, WC- and T2D-associated DMPs (Supplementary Material, Fig. S4A and B). Downmethylation at cg19693031 and cg02988288, both located at TXNIP, were the only common DMPs between T2D, hyperglycemia and insulin resistance (Supplementary Material, Fig. 5A), but are not found associated to WC. Downmethylation at cg19693031 was also the top association to FBG (P = 1.18 × 10⁻²⁵) (Supplementary Material, Fig. 5A and B). We find weak monotonic relationships (|Spearman’s ρ| < 0.15) between methylation intensity and its predictors (Supplementary Material, Fig. 5B—E). WC had the weakest correlation strengths among its associated DMPs (RMS Spearman’s ρ = 1.10 × 10⁻¹) (Supplementary Material, Fig. 5E).

Spurious transcription initiation at the TXNIP gene body

The three T2D-associated, TXNIP-located DMPs, cg19693031, cg26974062 and cg02988288, map to the eighth, sixth and fifth exons of the TXNIP gene body (Fig. 2). While methylation at promoter and first exon-associated CGIs strongly correlates to transcriptional silencing (26), the role CpG methylation at exons other than first is less clear. To explore the transcriptional effects to TXNIP, we utilized publicly available RNA-seq data from a study exploring interactome–transcriptome interactions in peripheral blood of diabetic and non-diabetic Han Chinese (27).

Figure 2

Genomic positions of T2D-associated DMPs located at the TXNIP body. CpG sites interrogated by the Infinium MethylationEPIC array are show together with T2D-associated DMPs. Transformed P values of associations to T2D, FBG, Hb-A1C, HOMA-IR and waist circumference is shown, with the red line signifying the significance cutoff. Mean β methylation intensities of diabetics are compared to controls. Average coverage is reported alongside variability measured as standard deviation. This is compared to average RNA-seq coverage of the gene region (from a Chinese cohort (27)) and the coverage variability. Predicted binding sites of miR-17 to the TXNIP 3′ UTR is plotted.

There is higher normalized coverage across the TXNIP gene region in patients with diabetes (Fig. 2) translating to significantly more transcripts (Supplementary Material, Fig. S6A). Exploring individual exon contributions, we find increased transcript levels for all exons other than the first (Supplementary Material, Fig. S7). Some studies suggest that gene body methylation directs isoform switching (28). TXNIP has two protein-coding isoforms differing by their N-termini and 5′ UTR (Fig. 2); however, we find no difference in the isoform ratios (Supplementary Material, Fig. S6B). Interestingly, we observe higher coverage variability at the 3′ UTR of controls (Fig. 2). Inspection of the individual coverage plots reveals that over half of controls exhibit depletion at the regions of high variability (Supplementary Material, Fig. S8), but this is not observed in any of the patients with diabetes, hinting that they may constitutively express a spurious intragenic transcript mapping to a portion of the TXNIP 3′ UTR.

Gene body methylation represses spurious intragenic transcription in concert with H3K36me3 (29). Orthogonal eFORGE analysis suggests that the three TXNIP-mapping DMPs co-locate at H3K36me3 regions, agreeing with imputed H3K36me3 regions of peripheral T cells (Supplementary Material, Fig. S9). Spurious transcription from demethylated gene regions are mediated by the Sp1 and Ets transcription factors (29); by checking DMP locations with respect to ENCODE transcription factor binding sites, we find that Sp1 colocalizes at cg02988288 in the K562 cell line (Supplementary Material, Table S7). Taken together, the data suggests the presence of a spurious transcript at the TXNIP 3′ UTR, but this will require additional experiments to establish definitively. Expressed sequence tags (ESTs) curated by ENSEMBL show two predicted non-coding transcripts mapping to TXNIP; however, these do not correspond to the regions of high transcriptional variability (Fig. 2).

Higher TXNIP transcript levels measured across exons 1–7 in diabetics suggests upregulation independent of spurious transcription events at the 3′ UTR (Supplementary Material, Fig. S6C). We hypothesized that spurious 3′ UTR transcripts contribute to increased expression by acting as dummy targets for microRNA regulation, especially since TXNIP mRNA is tightly regulated by miR-17 (30, 31). We predicted miR-17 binding to the TXNIP 3′ UTR and find colocalization at the high-variability region (Fig. 2).

Discussion

Methylation profiles in T2D, hyperglycemia and insulin resistance

Earlier studies reported differences in the peripheral blood methylation profiles in patients with T2D across different populations (17–21). Our top T2D-associated DMP agrees with those studies. We find statistically significant downmethylation at the same cg19693031 position located at TXNIP (17–20, 25). We also identify downmethylation at two other DMPs at the TXNIP gene body. Concurrently, we find a downmethylated DMP at CPT1A, but find no DMPs at other common associations, such as PHOSPHO1, ABCG1, SREBF1 and SOCS3 (25). Interestingly, differential methylation at two TXNIP CpGs is the only methylome change shared among associations to T2D, and measures of hyperglycemia (FBG and Hb-A1C) and insulin resistance (HOMA-IR). There is little overlap among associations to T2D and WC, a proxy measurement of obesity. The sparse overlap of T2D-associated DMPs with hyperglycemia, insulin resistance and obesity DMPs suggests that methylome changes based on T2D are distinct from changes based on hyperglycemic, insulin-resistant and obese states.

While some studies explored the use of peripheral blood methylation as proxies for disease-relevant tissues (32, 33), others argue that the enrichment of DMPs related to inflammatory and immunity pathways are an artifact of cell composition and do not reflect true associations to T2D (25). Our analysis reveals enrichment of inflammation-related pathways but find that this is partly explained by the overlap of immunity and inflammatory pathways with metabolism and cell-cycle pathways. We show evidence suggesting that the same group of genes play roles in cell cycle and immunity pathways. Given that GO assignments are enriched with terms related to metabolism, especially involving protein metabolism, and our ability to adjust for cell composition using complete blood counts and statistical inference, we believe that peripheral blood associations are not entirely artifacts of cell composition.

Interaction of T2D and Western immigration

T2D prevalence is higher among Filipino immigrants in the USA in contrast to Filipinos living in the Philippines and is higher in contrast to a Caucasian reference population (7). While body fat distribution accounts for increased susceptibility of Asians to T2D and metabolic diseases holding BMI constant (34), it does not fully account for the increased incidence among Filipino immigrants. We find through T2D–immigration interaction, downmethylation at MAP2K7, PRMT1 and ADM5. MAP2K7, encoding MKK7 which is part of the MAPK pathway, is involved in transducing signals in response to oxidative stress- and endoplasmic reticulum (ER) stress-induced apoptosis (35, 36), and activation of inflammatory responses via NFκB activation (37). Meanwhile, PRMT1 encodes protein arginine N-methyltransferase 1 (PRMT1), which is a protein modification enzyme that is involved in glucotoxicity-induced dysfunction of β-cells (38) and TXNIP-induced hepatic lipogenesis and inflammation (39). While methylation at these genes is significantly different in subjects with diabetes, a significant T2D–immigration interaction suggests that they may play a more consequential role in immigrants. We also find upmethylation near two lincRNAs, though these non-coding RNAs are yet to be properly functionalized.

Regulation of TXNIP by DNA methylation

Considering gene body methylation does not correlate with transcriptional repression (26), consistent T2D-, hyperglycemia- and insulin resistance-associated downmethylation of the TXNIP gene body observed in our study and in literature cannot adequately explain upregulation through transcriptional derepression. Indeed, glucose-induced expression of thioredoxin-interacting protein (TXNIP) is mediated by CRE carbohydrate response element-binding protein (ChREBP) which promotes TXNIP transcription in the presence of glucose-6-phosphate (40). Through examination of RNA-seq data and ENCODE data, we propose that gene body downmethylation of TXNIP promotes spurious intragenic transcription of the 3′ UTR. This is functionally relevant as TXNIP mRNA is heavily regulated by microRNA, notably miR-17, which we predict to bind to the suspected spurious transcript. Notably, miR-17 repression is the mechanism behind TXNIP upregulation and subsequent NLRP3 inflammasome activation in response to the unfolded protein response (31).

Conceivably, spurious 3′ UTR transcripts sequester miR-17-directed degradation by acting as a competitive endogenous RNA (ceRNA), resulting in increased TXNIP mRNA stability. TXNIP upregulation itself is important in T2D pathology by promoting hyperglycemia-induced apoptosis of pancreatic β cells, inhibiting insulin production in β cells, promoting inflammatory responses and inhibiting glucose internalization in skeletal muscle (31, 41–43). Moreover, TXNIP ablation is protective of T2D in obese mice models (44). Considering TXNIP’s ability to upregulate its own expression through positive feedback (45), initial ChREBP deregulation of TXNIP expression arising from blood glucose excursions may be exacerbated by feed-forward mechanisms and spurious transcript disruption of miR-17 targeting, leading to worsening deregulation even after the loss of initial stimulation. This may mean that TXNIP dysregulation persists in patients with history of T2D but with properly managed blood glucose, where it may still impact insulin resistance, inflammation and glucose response. We note that we were unable to directly pair our methylation data with transcriptome data based on unavailability, but given the ubiquity of TXNIP gene body downmethylation across different populations with diabetes, it is likely that some level of TXNIP methylation dysregulation is present in the methylomes of the source subjects with diabetes.

Conclusion

Through examination in the peripheral blood of Filipinos with T2D, we found methylation changes that primarily targeted genes involved with metabolism and inflammation that were distinct in contrast to changes associated with hyperglycemia, insulin resistance and obesity. We present a case for further effort in the functionalization of these DMPs through the examination of methylation changes at the gene body of TXNIP, which we hypothesize to dysregulate TXNIP expression by spurious ceRNA transcription. Given the role of TXNIP in multiple pathways leading to diabetes as found in multiple in vitro cell and animal studies, the identification of TXNIP dysregulation in the peripheral blood of Filipino patients with T2D highlights the importance of translating basic science research.

Materials and Methods

Ethics statement

Written informed consent was obtained from all participants. The research protocols used in this study were approved by the University of the Philippines Manila Research Ethics Board (UPMREB) and the Institutional Review Board (IRB) at the University of California San Francisco, USA.

Study participants

General inclusion criteria for enrolment into the study included Filipino descent and a minimum age of 18. The study investigated 198 Filipinos enrolled at the Philippine General Hospital (PGH; Manila, Philippines) and 167 Filipino immigrants enrolled at the Zuckerberg San Francisco General Hospital (ZSFG; San Francisco, California, USA). Filipino participants were identified from enrolment in another T2D study at the PGH, referral from T2D community advocates and house-to-house visits in the Pandacan and Paco communities in Manila City. Immigrant volunteers were identified through outreach at Filipino cultural and social events, advertisements in and around the ZSFG and direct contact with ZSFG outpatients.

Anthropometric measurements and fasting blood tests were taken from the prospective volunteers. Anthropometric measurements included height and weight, and waist and hip circumference. Blood testing included a complete blood count, a lipid profile and measurements of 12-h fasting blood glucose (FBG), glycated hemoglobin (Hb-A1C) and insulin. Insulin and FBG measurements were used to calculate the homeostatic model assessment for insulin resistance (HOMA-IR) (46). Prospective controls were administered an additional 2-h oral glucose tolerance test (OGTT), wherein participants drank a 75-g oral glucose beverage and abstained from taking in food or any other beverage until blood was drawn 2 h post-ingestion for blood glucose measurement.

Inclusion criteria for controls were per the 2004 American Diabetes Association (ADA) guidelines for blood glucose profiles: an FBG of less than 100 mg/dl and a 2-hr OGTT of less than 140 mg/dl (47). A past diagnosis or a family history of T2D was the exclusion criterion for controls. The inclusion criterion for T2D cases was per the 2004 ADA guidelines: an FBG of greater than or equal to 126 mg/dl or a 2-h OGTT of greater than or equal to 200 mg/dl. Participants with a self-reported diagnosis of diabetes were included as T2D cases, regardless of their FBG and OGTT measurements. Individuals with prediabetes were excluded from the analysis.

Methylation profiling

DNA was extracted from stored patient blood samples using the QIAamp® DNA Blood Mini Kit (Qiagen, Hilden, Germany) and then bisulfite-converted using the EZ DNA Methylation Kit (Zymo Research, Irvine, CA). Bisulfite-converted DNA was run using the Infinium MethylationEPIC BeadChip array (Illumina, San Diego, CA) according to the manufacturer’s protocol at the Microarray Core Laboratory, Institute of Human Genetics, National Institutes of Health, University of the Philippines Manila.

Raw intensity files were processed with R-3.5.2 (48) using the minfi-1.28.3 package (49). Background correction was performed using Illumina’s method, followed by quantile normalization as implemented by minfi. Cell type proportions were estimated from raw intensity files using the reference-based Houseman method (50) as implemented by the FlowSorted.Blood.EPIC-1.0.0 package (51). Sample quality was assessed by calculating the mean detection P value and the sex concurrence. Samples were excluded if the mean detection P value exceeded 0.05, or if the sex predicted by minfi disagreed with the clinically recorded sex. Also, samples with ambiguous sex assignment were removed, which were defined as samples whose difference between median Y intensities and median X intensities fell between −0.5 and −3.0.

Probes were filtered based on individual detection quality and on literature recommendations. Individual probes were deemed to have a poor detection quality and removed if the detection P value passed above a 10⁻¹⁶ threshold (23). Also, probes with internal SNPs occurring within 5 bp from the 3′ end, probes mapping to the sex and mitochondrial chromosomes and probes with off-target mapping were filtered out (24). The remaining probes were annotated with positions on the GRCh38/hg38 genome, CpG island contexts and GENCODE v22 genes based on annotations by Zhou et al. (24). Gene position annotations were taken from the Illumina’s MethylationEPIC v1.0 B4 manifest.

Memory intensive calculations were performed with the help of the high-performance computing facilities at the Core Facility for Bioinformatics (PGC-CFB), Philippine Genome Center, University of the Philippines System.

Association experiments

Linear regression was done with R-3.5.2 (48) using the limma-3.38.3 package (52) on the log₂-transformed M values (53). A model was constructed considering age and sex as covariates, with blocking of immigration status, array run date and cell type proportion. Within-array differences were accounted by including the first 30 principal components of the control probes (23). Continuous variables such as age, WC, hyperglycemic measures, insulin resistance and cell type compositions were modeled using natural cubic splines with three degrees of freedom as implemented by the splines-3.5.2 package. The linear models used for each association experiment found in this manuscript are described in Supplementary Material, Table S1.

Statistical inference was performed on the fitted linear models using the moderated t test as implemented by the limma package. Fold changes reported by limma were taken as the effect size for discrete variables such as phenotype. For continuous variables, monotonic relationships and effect sizes were calculated by performing Spearman correlation between methylation intensity and its predictors.

GO and pathway analysis

Gene ontologies were assigned to each probe in R-3.5.2 (48) using the missMethyl-1.16.0 package (54). Enrichment of DMPs for each ontology term was calculated using the hypergeometric test.

Pathway enrichment was done in R-3.5.2 using the pathfindR-1.3.0 package (55). Underlying protein–protein interaction networks were identified from a list of differentially methylated genes using an implementation of the greedy algorithm by the pathfindR package with a maximum search depth of 2. Enriched pathways were clustered into representative pathways using the clustering and partitioning algorithms implemented by pathfindR.

RNA-seq analysis

Raw RNA sequencing data from a 20-patient cohort of Han Chinese with 10 healthy controls and 10 patients with T2D was downloaded from the Short Read Archive (27). Paired-end sequencing reads were trimmed with trimmomatic-0.36 (56) using a 4 bp sliding window with a minimum mean window quality of 20. The trimmed reads were aligned to chromosome 1 of the GRCh38/hg38 assembly using the bwa-0.7.12 mem algorithm (57).

The coverage from the resulting BAM files was normalized into 5 bp bins using the bamCoverage program from the deepTools 2.0 software package (58). Coverage variability was computed by taking the standard deviation for each 5 bp coverage bin in the genome.

microRNA binding prediction

Binding sites for miR-17 was predicted using TargetScan by determining conserved miRNA family binding sites to the 3′ UTR of TXNIP (59). Positions in the 3′ UTR were converted to hg38 coordinates for plotting.

Adjustment for multiple-hypothesis testing

All reported P values were adjusted using the Benjamini–Hochberg procedure (60).

Acknowledgements

We would like to acknowledge the Philippine Genome Center for providing the computational resources used in the project. We would also like to acknowledge Dr Jezmhyn Marie Paelmo, Dr Zyrelle Avienn Santos, Ms Arlyn Napeñas, Mr Matthew Roces, Ms Irah Camil Rubio, Dr J. Enrico Lazaro and Dr Pia Bagamasbad.

Conflict of Interest statement. None declared.

Funding

Commission on Higher Education—Philippine–California Advanced Research Institutes (CHED-PCARI IHTM 2015-03).

Author notes

The authors wish it to be known that, in their opinion, the first two authors, Mr Albao and Dr Cutiongco-de la Paz, should be regarded as joint First Authors.

The authors wish it to be known that, in their opinion, Dr Cutiongco-de la Paz and Dr Seielstad be joint Corresponding Authors.