Cross-Ancestry DNA Methylation Marks of Insulin Resistance in Pregnancy: An Integrative Epigenome-Wide Association Study

Insulin resistance is a hallmark of type 2 diabetes and is characterized by a reduced metabolic response to insulin in skeletal muscle, liver, and adipocytes (1). Insulin resistance is usually compensated by increased insulin secretion, but type 2 diabetes develops when pancreatic β-cells are unable to compensate (2). In a normal pregnancy, insulin resistance increases to secure an adequate supply of glucose to the fetus (3,4). In the third trimester, many women reach a level of insulin resistance similar to what is observed in type 2 diabetes, resulting in gestational diabetes mellitus (GDM) if the pancreatic β-cells fail to compensate (5). Although GDM usually resolves after delivery, affected women have a 10-fold increased risk for developing type 2 diabetes (6).

Currently, genetic variants explain 18% of type 2 diabetes risk (7). A large proportion of type 2 diabetes risk may be explained by the interplay between genetic and environmental factors such as obesity, physical activity, and smoking, among others (8). The relationship between environmental or genetic factors and type 2 diabetes is potentially mediated by epigenetic factors (9). Methylation at specific CpG sites in peripheral blood leukocytes has been associated with both incident and prevalent type 2 diabetes (9,10); with glucose homeostasis variables such as fasting glucose, fasting insulin, and insulin resistance in men and nonpregnant women (11–13); and with GDM in cord blood (14–16).

South Asians have an increased risk of developing GDM and type 2 diabetes in comparison with Europeans (17,18). Studies suggest that South Asians have high insulin resistance combined with reduced β-cell capacity (19,20). Cross-ancestry studies can be useful to study pathogenic mechanisms that may explain the differential risk between populations or to confirm cross-ancestry associations in diverse populations with distinct confounding structures. Still, most epigenome-wide association studies (EWAS) focus on populations of European ancestry, limiting the generalizability in other ancestries.

We have not found any published EWAS in which peripheral blood leukocytes were used to discover CpG sites associated with insulin resistance across populations with different ancestries. Furthermore, of the four studies in which associations were found between insulin resistance and epigenome-wide DNA methylation in whole blood (11,12,21,22), a replication cohort to verify findings was not used in three (11,12,22). Here, we aimed to identify cross-ancestry CpG sites in peripheral blood leukocytes associated with insulin resistance in gestational week 28, with subsequent testing for replication. We also use different data sources to leverage further knowledge about potential biological mechanisms. And we performed mediation analysis to explore the direction of effects.

The STORK Groruddalen (STORK G) study is a population-based cohort, including 823 healthy pregnant women without diabetes in gestational week 15, in the multiethnic area of Groruddalen in Oslo, Norway, 2008–2010 (23). Details about inclusion criteria are given in the Supplementary Material. The Epigenetics in Pregnancy (EPIPREG) sample include all women of European (n = 312) or South Asian (n = 168) origin participating in STORK G with available DNA collected in gestational week 28 (24).

The STORK G study, including its genetic and epigenetic data (ie, from the EPIPREG sample), is approved by the Norwegian Regional Committee for Medical Health Research Ethics South East (reference no. 2015/1035). Written informed consents were gathered from all participants.

Data were collected in gestational week 28. Venous blood (serum and plasma) was drawn into EDTA tubes. The samples were either aliquoted and stored in a biobank or subjected to routine laboratory analyses performed during the study period. Insulin resistance was estimated by the HOMA, using the HOMA2 calculator version 2.2.2 (https://www.dtu.ox.ac.uk/homacalculator) based on fasting levels of glucose (HemoCue, Angelholm, Sweden) and C-peptide (DELFIA, PerkinElmer Life Sciences, Wallac Oy, Turku, Finland) (25). The interassay coefficient of variability for C-peptide was 4–5% between 336 and 3,500 pmol/L. The 2-h glucose was measured with a point-of-care device (HemoCue). We calculated the HOMA for insulin resistance (HOMA-IR) with C-peptide instead of insulin because a recently published study from our group indicated that hepatic insulin clearance differed between European and South Asian women (26).

All women underwent a 75-g oral glucose tolerance test. GDM was diagnosed with the World Health Organization (WHO) 1999 criteria (i.e., fasting glucose ≥7.0 mmol/L or 2-h glucose ≥7.8 mmol/L). We also present a slightly modified version of the current WHO 2013 criteria (fasting glucose 5.1–6.9 mmol/L or 2-h glucose ≥8.5–11 mmol/L; 1-h glucose information was not available in our data). Fasting insulin was measured with noncompeting immunofluorometric assays (DELFIA, PerkinElmer Life Sciences). The interassay coefficient of variability for insulin was 7% between 63 and 580 pmol/L. Fasting plasma total cholesterol, LDL cholesterol, HDL cholesterol, and triglycerides levels were measured with a colorimetric method (Vitros 5.1 fs, Ortho Clinical Diagnostics, Neckargemünd, Germany).

Age of the women was calculated from date of birth reported at enrollment. With data collected in gestational week 28, we calculated BMI from measured height using a fixed stadiometer at enrolment, and measured body weight (Tanita-BC 418 MA, Tanita Corporation, Tokyo, Japan). Systolic and diastolic blood pressures were measured with the M6 Comfort HEM-7000-E device (Omron, Kyoto, Japan). Smoking status was assessed using an interviewer-administered questionnaire and dichotomized as follows: smokers (current and 3 months before pregnancy) versus nonsmokers (former smokers and never smokers).

Details about blood collection and DNA isolation and preparation can be found in the Supplementary Material. We quantified DNA methylation in maternal peripheral blood leukocytes in gestational week 28 using the Infinium MethylationEPIC BeadChip (Illumina, San Diego, CA). Quality control (QC) was performed with the meffil R package (27). We removed six outliers from the methylated to unmethylated ratio comparison (>3 SD), one outlier in control probes bisulfite 1 and bisulfite 2 (>5 SD), and one due to sex mismatch (predicted sex outliers >5 SD). We removed probes with a bead count of fewer than three and detection P value < 0.01. In total, 472 of the 480 available individuals, and 864,560 probes passed the QC. To remove technical variation, all samples were randomly distributed across chips. Furthermore, we used functional normalization implemented in the meffil R package (27), which returns normalized data standardized for batch effects of slide, row, and columns, as well as 10 principal components from the QC (25). Proportions of blood cells (namely, CD8T, CD4T, NK, β-cells, monocytes, and neutrophils) were calculated using the Houseman method (27) during the QC. After further removal of probes with infinites, missing values, Y chromosome probes, probes harboring single nucleotide polymorphisms (SNPs) and cross-reactive probes (28), 810,266 probes remained. Technical validation of the EPIC bead chip showed good agreement between pyrosequencing and the EPIC array (24).

The Illumina CoreExome chip was used for genotyping (Supplementary Material). Briefly, after QC, 300 Europeans and 138 South Asians had both genotyping and epigenetic data. Genetic ancestry was assessed through ancestry principal component analysis, which corresponded with the participants’ reported country of birth or mothers’ country of birth if the mother was born outside of Europe (24). Genotype imputation in Europeans and South Asians was done separately by using the respective imputation panel from the 1000 Genomes Project (29).

For the EWAS, we included 294 European and 162 South Asian women (n = 456) who passed the EWAS QC and did not have HOMA-IR outlying values (greater than ±3 SDs) in their respective ethnic group (Fig. 1). Of these, 287 European and 135 South Asian women had genetic data available and were included in the methylation quantitative trait loci (meQTL) analysis (Fig. 1).

The Born in Bradford (BiB) cohort (n = 430 Europeans and 449 South Asians) (30,31) included pregnant women with DNA methylation data analyzed in peripheral blood leukocytes in gestational weeks 26–28 and insulin resistance calculated with HOMA2. The Methyl Epigenome Network Association (MENA) (n = 320) and the Finnish Botnia cohort (n = 56) included nonpregnant women and men with DNA methylation in peripheral blood leukocytes and insulin resistance calculated with the standard insulin resistance formula. Additional details about the cohorts are given in the Supplementary Material and Supplementary Table 1.

Beta values were logit transformed to M values. For the main discovery analysis, we performed a cross-ancestry EWAS of insulin resistance with mixed models regression using the R packages lme4 and lmerTest (32) and the P value calculation of Satterthwaite. Insulin resistance was treated as the exposure (independent variable) and DNA methylation as the outcome (dependent variable). Blood cell composition, age, and smoking status were treated as fixed effects, and ancestry was treated as a random intercept. We treated ancestry as a random intercept to overcome potential ancestry-related differences in DNA methylation. We also performed EWAS separately on data from Europeans and South Asians to identify potential ancestry-specific CpG sites (Supplementary Material). Because only three South Asian participants were classified as smokers, a model without adjusting for smoking was also performed. For the EWAS analyses, we used a false discovery rate (FDR) threshold of 5%. We considered a CpG as replicated if it had a consistent direction of the effect with the discovery analysis and a Bonferroni corrected P value (0.05 per number of CpG sites tested for replication) in at least one replication cohort. Some of the mixed models we present had singular fit (i.e., overfitting due to the number of categories in the random intercept). Nevertheless, we do not consider this to be a major problem for the presented conclusions (Supplementary Material).

To assess if GDM influenced the identified associations, we performed a subanalysis in women without GDM; this was also done in BiB. Some of the mixed models we present had singular fit. As noted in the preceding paragraph, we do not consider this to be a major problem for the presented conclusions (Supplementary Material).

For gene ontology analysis, we included cross-ancestry CpG sites with a P value <1 × 10⁻⁴ and accepted a 5% FDR threshold (Supplementary Material).

For the association between the replicated cross-ancestry CpG sites and other cardiometabolic phenotypes, we used linear mixed models with the same covariates as in the discovery analysis. We additionally adjusted for BMI to evaluate its influence on the association between the discovered CpG sites and insulin resistance.

To identify genetic variants associated with the cross-ancestry CpG sites, we used a cross-ancestry approach. Although genetic variants can have differential effects across populations (33), it has been suggested that genetic variants associated with complex traits in more than one ancestry tend to have consistent direction of effects (34,35) and are more likely to be causal or in linkage disequilibrium with causal variants (33–35). Therefore, we expect genetic variants associated with methylation of CpG sites across ancestries to have similar direction of effects.

From the variants that passed postimputation QC in data from European (n_SNP = 5,063,268) and South Asian (n_SNP = 3,939,192) participants, we selected the variants common in both ancestries (n_SNP = 3,506,750) for the cross-ancestry meQTL analysis. To save computational power, we performed linear regression separately in data from European and South Asian participants, using the R package GEM (36), adjusted for blood cell composition, age, and smoking. The coefficients and SEs were combined using an inverse-variance weighted fixed-effects meta-analysis with custom R scripts based on METAL (37). Cis-meQTLs were defined as positioned less than ±1 Mb from the methylation site or otherwise were classified as trans-meQTLs. We report associations at P < 1.0 × 10⁻⁸ (5 × 10⁻⁸ per five meQTL analyses).

The most significant cross-ancestry meQTLs of each CpG site were examined for associations with other cardiometabolic phenotypes and mediation analysis. If some of the most significant genetic variants were located on the same gene and/or chromosome, we evaluated if they were in strong linkage disequilibrium (R² > 0.8) (Supplementary Material).

To evaluate the association between a genetic variant with other cardiometabolic phenotypes we used linear mixed models with ethnicity as a random intercept.

Mediation analysis was performed to evaluate if DNA methylation was a mediator between the cross-ancestry meQTLs and insulin resistance (Supplementary Material). The mediation analysis calculates three models: 1) the average causal mediation effect, which is the effect size that goes through the mediator (indirect effect); 2) the average direct effect, which is the effect size after adjustment for the mediator, thus removing its effect (direct effect); and 3) the total effect, which is the sum of the indirect and direct effects.

In the Skeletal Muscles, Myokines, and Glucose Metabolism (MyoGlu) study (38–40), we used Spearman analysis to test correlations between gene expression in subcutaneous adipose and muscle tissue that mapped to the replicated CpG sites and insulin sensitivity measured with the hyperinsulinemic euglycemic clamp. More details about MyoGlu are presented in the Supplementary Material.

The BIOS database (https://wiki.bbmri.nl/wiki/BIOS_start-) was consulted to identify expression quantitative trait methylation in peripheral blood leukocytes. The GTEx database (https://gtexportal.org/home/) was consulted to evaluate if the cross-ancestry meQTLs for each CpG site were also expression quantitative loci (eQTL) in peripheral blood leukocytes, subcutaneous adipose tissue, omental adipose tissue, liver, skeletal muscle, and pancreas. Data from the Atlas of GWAS Summary Statistics (https://atlas.ctglab.nl/) were extracted for cross-ancestry meQTLs to identify associations (P < 0.05) with cardiometabolic traits (details about study selection are provided in the Supplementary Material).

Because of strict regulations for genetic data and privacy protection of patients in Norway, all requests for data access are processed by the STORK G project’s steering committee. Please contact the principal investigator of STORK G ([email protected]) or the principal investigator of EPIPREG ([email protected]).

BiB data used for this submission will be made available on request to the BiB executive committee ([email protected]). The BiB data management plan (https://borninbradford.nhs.uk/research/how-to-access-data/) describes in detail the policy regarding data sharing, which is through a system of managed open access.

The clinical characteristics of the 294 European and 162 South Asian women (N = 456) included in this study are presented in Table 1.

We identified six cross-ancestry CpG sites inversely associated with insulin resistance (Table 2 and Fig. 2A), with lower levels of DNA methylation associated with higher insulin resistance. We did not find evidence of inflation (λ = 1.02) (Supplementary Fig. 1). Summary statistics for the 100 cross-ancestry CpG sites with P ≤1 × 10⁻⁴ are shown in Supplementary Table 2. In cg10318066, we detected a methylation outlier, but we decided to pursue this site because the direction of the effect was unchanged after omitting it from the analysis (effect = −0.045; SE, 0.011; P = 1.65 × 10⁻⁵).

In the ancestry-specific EWAS of insulin resistance, we identified 12 European-specific CpG sites (Supplementary Table 3), whereas no CpG sites reached the FDR threshold in the South Asian EWAS (Supplementary Table 4). The South Asian EWAS unadjusted for smoking yielded similar results in comparison with the adjusted model, and no CpG site reached the FDR threshold. The EWAS of data from European participants presented moderate inflation (λ = 1.33) (Supplementary Fig. 2), whereas no inflation seemed to be present in the EWAS of data from South Asian participants (λ = 1.03) (Supplementary Fig. 3).

All the cross-ancestry CpG sites but cg10318066 were replicated in at least one of the three independent cohorts (Table 2). None of the 12 ancestry-specific CpG sites observed in Europeans were replicated in any of the independent cohorts (Supplementary Table 5).

We performed sensitivity analyses on the replicated CpG sites by excluding individuals with GDM. When GDM cases according to WHO 1999 criteria were excluded, most of the coefficients weakened as they increased by 12–18% (Supplementary Table 6). When GDM cases according to WHO 2013 criteria were excluded, all coefficients strengthened as they decreased between 2% and 15% (Supplementary Table 6). Overall, the direction of effects did not change and remained nominally significant (P < 0.05). The same analysis was conducted in the BiB cohort, wherein similar effects were observed (Supplementary Table 6).

To identify potentially relevant biological processes, we performed gene ontology enrichment analysis in the cross-ancestry (P < 1 × 10⁻⁴) CpG sites, which annotated to 80 loci. No gene ontology terms reached an FDR < 0.5 (Supplementary Table 7).

Methylation at the replicated CpG sites were inversely associated with C-peptide, insulin, and BMI values, and additionally with at least one phenotype (Table 3 and Fig. 2B).

Because the five replicated CpG sites were associated with maternal BMI, we further adjusted the regression models for BMI. The effect sizes weakened by 8–21%, but the direction of the effects persisted and remained nominally significant (P < 0.05) (Supplementary Table 8).

From meQTL analysis in the EPIPREG sample, we obtained 478 genetic variants related to the five replicated cross-ancestry CpG sites (Supplementary Table 9). The most significant genetic variants from the cross-ancestry meQTL (Table 4) were further explored.

Because rs2297976, rs17387775, and rs34964576 associated with TNIPX sites were strongly correlated within each other (R² = 1 when using either the European or South Asian panels), only rs2297976 was explored further. In EPIPREG, the minor allele of rs2297976 (T) was associated with lower insulin resistance and C-peptide levels (Supplementary Table 10), whereas the minor allele of rs2179174 (G) and rs7006759 (C) had the opposite effect. The minor allele of rs2179174 (G) was also associated with higher fasting insulin and high systolic and diastolic blood pressures. The minor allele of rs7006759 (C) was associated with higher BMI, higher 2-h glucose levels, increased systolic and diastolic blood pressures, and GDM per WHO 1999 criteria (Supplementary Table 10).

We performed look-ups of the meQTLs in GWAS summary data for relevant outcomes (Supplementary Table 11). With nominal significance (P = 0.047 to 1.14 × 10⁻⁵), rs2297976 was associated with fasting glucose, 2-h glucose, BMI, fasting insulin, total cholesterol, and LDL cholesterol levels, and both systolic and diastolic blood pressures. The variant rs2179174 was associated with type 1 diabetes, HDL cholesterol level, and systolic blood pressure. The variant rs7006759 was not associated with any of the tested phenotypes.

The average causal mediation effects for all CpG sites were stronger than the average direct effects (Table 5), suggesting that the association between the meQTLs and insulin resistance was mediated by DNA methylation. When studying insulin resistance as a mediator between the meQTLs and DNA methylation, average direct effects were stronger than the average causal mediation effects for all the CpG sites (Table 5).

In MyoGlu, SLC7A11 expression was associated with higher insulin sensitivity (low insulin resistance) in muscle and with low insulin sensitivity in adipose tissue. ZSCAN26 expression was positively associated with high insulin sensitivity in both tissues. From expression quantitative trait methylation look-ups in data from the Biobank-based Integrative Omics Studies (BIOS) consortium higher methylation in cg19693031 in peripheral blood leukocytes was associated with lower TXNIP expression.

From eQTL look-ups in GTEx, the T allele of rs2297976 in SLC2A1 was associated with increased SLC2A1-AS1 expression levels in whole blood, pancreas, both subcutaneous and omental adipose tissue, and for SLC2A1 in pancreas (Supplementary Table 12). The C allele of the cis variant rs2179174 in ZSCAN26 was associated with increased expression of several genes from the zinc-finger family in all the tissues tested (Supplementary Table 12).

We have identified six cross-ancestry CpG sites inversely associated with insulin resistance during pregnancy in peripheral blood leukocytes, of which five were replicated in at least one of three independent cohorts. Both the CpG sites and their associated meQTLs were associated with several cardiometabolic phenotypes in the EPIPREG sample as well as in publicly available GWAS summary data. The meQTLs were also eQTLs in whole blood, subcutaneous adipose tissue, omental adipose tissue, liver, skeletal muscle, and pancreas. On the basis of mediation analysis and look-ups on gene expression, it can be inferred that the association between the meQTLs and insulin resistance was mediated by DNA methylation, presumably through modifying gene expression.

Low methylation of cg06690548 in SLC7A11 in peripheral blood leukocytes has been related to high fasting insulin levels in a nonpregnant population (13). SLC7A11 encodes a cysteine/glutamate plasma membrane transporter for the cellular uptake of cysteine in exchange for glutamate (41). Interestingly, cysteine, glutamate, and glycine conform glutathione, an antioxidant that may be deficient in type 2 diabetes, probably due to inefficient protein turnover (42). Notably, SLC7A11’s mRNA correlated positively in muscle and negatively in adipose tissue with insulin sensitivity, implying that in muscle, high cysteine/glutamate transport may be associated with low insulin resistance and with high insulin resistance in adipose tissue. Although low methylation in cg06690548 has been associated with high SLC7A11 mRNA in blood (43), the gene expression results should be interpreted with caution, because the association between SLC7A11’s mRNA and insulin resistance in peripheral blood leukocytes is currently unknown.

In accordance with our phenotype association analysis, cg06690548 has been related to BMI (44) and systolic and diastolic blood pressures (43). Notably its meQTL rs7006759 is located in PSD3, a gene that harbors genetic variants that have been associated with type 2 diabetes (45). Because rs7006759 and cg06690548 are located on different chromosomes, it is not clear how they may interact to regulate insulin resistance. Interchromosomal associations are the most common reported trans-meQTLs, and it is possible that these trans associations are in long-range linkage disequilibrium with other cis-meQTLs (46). Another possible mechanism is that trans-meQTL could target genes due to chromatin loops which could make them spatially close, thus mimicking a cis-meQTL effect (47). However, whether these mechanisms are true for rs7006759 and cg06690548 remains to be elucidated.

Regarding cg04861640 in ZSCAN26, a previous study reported that methylation of cg04861640 in peripheral blood leukocytes was decreased in patients with diabetes type 1 who had nephropathy (48), suggesting a possible role in the detrimental effects of diabetes. Notably, cg04861640 has been related to rheumatoid arthritis (49) and Sjögren syndrome (50), and its cis-meQTL, rs2179174, is associated with rheumatoid arthritis (51); both are autoimmune diseases recognized by inflammation. These and other autoimmune disorders have been associated with increased type 2 diabetes risk (52). Hence, type 2 diabetes may share mechanisms with a wide spectrum of autoimmune disorders, although the mechanisms are not well understood (52).

Because low methylation in promoters is usually associated with increased transcription (53), we can infer that lower methylation of cg04861640 might lead to higher expression of ZSCAN26, which, in turn, would increase insulin resistance. However, contrasting these possible relationships, high ZSCAN26 mRNA was associated with increased insulin sensitivity in the MyoGlu data. The effect of ZSCAN26 mRNA in peripheral blood leukocytes on insulin resistance should be elucidated to make further conclusions. Moreover, rs2179174 was an eQTL for pancreas, blood, and both omental and subcutaneous adipose tissue for several genes of the zinc-finger family, suggesting an effect on transcription of ZSCAN26 across tissues, although the implications in the context of insulin resistance and related traits are unknown.

Three of the cross-ancestry CpG sites were located in TXNIP, a gene that may induce β-cell apoptosis at high glucose levels, which may promote development of diabetes (54). Consistent with our findings, reduced methylation of cg19693031 in peripheral blood leukocytes was associated with higher insulin resistance in a general Mexican-ancestry, nonpregnant population (12). Furthermore, reduced blood methylation in cg19693031 has been observed in type 2 diabetes in a recent European-based EWAS meta-analysis (9) and in a cross-ancestry EWAS of Europeans and Indian Asians (55). Furthermore, reduced methylation of cg02988288 and cg26974062 in peripheral blood leukocytes has been associated with type 2 diabetes in Filipino individuals (56). In cord blood, reduced methylation at cg02988288 and cg26974062 was positively associated with glucose area under the curve (16) and a principal component composed of both insulin and glucose variables from an oral glucose tolerance test (57).

Low methylation in peripheral blood leukocytes of the three TXNIP CpG sites is related to high TXNIP transcription in type 2 diabetes (56), suggesting TXNIP’s mRNA levels increase in conditions with insulin resistance. TXNIP’s meQTL rs2297976 in SLC2A1 encodes the GLUT1 receptor, which is a major glucose transporter ubiquitously expressed in all tissues, especially in the blood-brain barrier and erythrocytes (58). Interestingly, SLC2A1 has a negative relationship with TXNIP’s mRNA, because SLC2A1 expression is repressed by TXNIP in the presence of glucose (55). The rs2297976 variant was also an eQTL for pancreatic tissue; hence, it is possible that rs2297976 regulates TXNIP transcription through DNA methylation in pancreas.

In line with our phenotype associations, reduced methylation levels at cg19693031 in peripheral blood leukocytes have been related to fasting glucose (13), triglyceride (59), and systolic blood pressure (60) values. Furthermore, methylation at cg19693031 is also reduced in people with metabolic syndrome (61). Thus, cg19693031 is consistently associated with a wide spectrum of cardiometabolic traits.

We replicated five cross-ancestry sites in at least one replication cohort, but only cg19693031 had a concordant direction of the effect across all three replication cohorts. This could be due to underlying differences between cohorts, or to statistical power. Most of the replicated sites were found in BiB, which had the largest sample size, an equivalent ancestry composition, and pregnant women. Furthermore, MENA and Botnia included both sexes, which may have played a role. None of the CpG sites identified in the European EWAS were replicated, suggesting that these may be false positives driven by the moderate inflation observed in this analysis.

Mediation analysis indicated that all five CpG sites mediated the association between the meQTLs and insulin resistance. That gene expression of TXNIP is associated with type 2 diabetes in peripheral blood leukocytes, whereas gene expression of SLC7A11 and ZSCAN26 in both muscle and adipose tissue is associated with insulin sensitivity, may suggest that the meQTLs alter DNA methylation, subsequently altering gene expression, which ultimately may affect insulin resistance (Fig. 3 and Supplementary Fig. 4). However, transcription levels may not necessarily directly affect insulin resistance, because other possible pathways and traits may be involved. The meQTLs were not strongly associated with insulin resistance independently of DNA methylation; hence, other environmental factors may influence the methylation of these CpG sites, which, in turn, could influence insulin resistance. Hence, Mendelian randomization studies with robust instruments should be performed to explore causality.

Major strengths of this work are the cross-ancestry approach and the successful replication in independent cohorts. Integrating data from publicly available sources enabled us to leverage more information about the CpG sites identified in this study. Important limitations include the cross-sectional design and observational nature, which make it difficult to draw a conclusion about the directionality of the associations. Because of the relatively small sample size, we are only able to detect large and medium effect sizes for the association between DNA methylation and insulin resistance. Insulin resistance was calculated with different methods across the replication cohorts, and MENA and Botnia included both men and nonpregnant women, which may partly explain why some sites were not consistently replicated across all cohorts.

In conclusion, the cross-ancestry EWAS identified five CpG sites associated with insulin resistance, which were replicated in independent cohorts. These five CpG sites showed a decrease in methylation at high insulin resistance levels and were genetically regulated. Both the CpG sites and their meQTLs showed associations with several cardiometabolic phenotypes, whereas the meQTLs were eQTLs for different tissues. Furthermore, expression of the genes in which the CpG sites are annotated showed associations with cardiometabolic traits. These relationships, along with the mediation analysis, suggest that the gene variants may affect insulin resistance levels through DNA methylation, presumably by altering gene transcription.

Acknowledgments. The authors thank the women who participated in the STORK Groruddalen study, Maria Sterner, Malin Neptin, and Gabriella Gremsperger at the Genomics Diabetes and Endocrinology CRC, Malmö, for the wet laboratory experiments of the bead chips, and Leif C. Groop, Lund University Diabetes Centre, Malmö, Sweden, for facilitating the wet laboratory experiments. The BiB study is only possible because of the enthusiasm and commitment of the Children and Parents in BiB. We are grateful to all the participants, teachers, school staff, health professionals and researchers who have made BiB happen.

Funding. EPIPREG is supported by the South-Eastern Norway Regional Health Authority (grant 2019092), and the Norwegian Diabetes Association. BiB receives core infrastructure funding from the Wellcome Trust (WT101597MA), a joint grant from the UK Medical Research Council (MRC) and Economic and Social Science Research Council (MR/N024397/1), the British Heart Foundation (CS/16/4/32482), and the National Institute for Health Research under its Applied Research Collaboration for Yorkshire and Humber and the Clinical Research Network. DNA methylation data were funded by the UK MRC via the Integrative Epidemiology Unit (MC_UU_12013/5). H.R.E. works in the MRC Integrative Epidemiology Unit at the University of Bristol, which is supported by the MRC and the University of Bristol (MC_UU_00011/5). R.P.B. received funding form EFSD/Novo Nordisk Programme for Diabetes Research in Europe, the Swedish Research Council (2021-02623), and the Crafoord Foundation. The Botnia Study has been financially supported by grants from the Sigrid Juselius Foundation, Folkhalsan Research Foundation, Nordic Center of Excellence in Disease Genetics, EU (EXGENESIS), Finnish Diabetes Research Foundation, Foundation for Life and Health in Finland, Finnish Medical Society, Helsinki University Central Hospital Research Foundation, Perklén Foundation, Ollqvist Foundation, and Narpes Health Care Foundation. The study has also been supported by the Municipal Heath Care Center and Hospital in Jakobstad and Health Care Centers in Vasa, Narpes, and Korsholm.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

The funders played no role in the study design and conduct, in these analyses, or in the decision to submit the manuscript for publication.

Author Contributions. N.F.-B., C.S., and R.B.P. contributed to the study conceptualization and design of this sub-study. C.S., E.Q., and K.I.B. conceptualized and C.S. and K.I.B. designed the EPIPREG sample. C.S., E.Q., and K.I.B. acquired genotype and DNA methylation data in STORK G. N.F.-B. conducted the statistical analyses in EPIPREG, drafted the manuscript, performed the analyses in MENA, and did the postimputation QC. H.R.E. contributed to study design and conducted the analysis in the BiB cohort. R.B.P. conducted the analysis in Botnia and realized the wet laboratory experiments regarding the methylation and genotyping chips. S.L.-Ø. curated EPIPREG data and performed the correlation analyses in the MyoGlu study. J.O.O. contributed to the data visualization and interpretation. G.-H.M. performed the QC of the genomic data and imputation in EPIPREG. K.I.B. and A.K.J. designed the STORK G project. L.S. contributed with data acquisition in STORK-G. C.A.D. and K.I.B. designed and administered the MyoGlu study. C.S. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis. All authors revised the manuscript critically and approved the final version.

Prior Presentation. This work was presented at the 57th European Association for the Study of Diabetes Annual Meeting, virtual, 28 September to 1 October 2021.