Epigenetic Link Between Statin Therapy and Type 2 Diabetes
Abstract
OBJECTIVE To investigate the role of epigenetics in statins’ diabetogenic effect comparing DNA methylation (DNAm) between statin users and nonusers in an epigenome-wide association study in blood.
RESEARCH DESIGN AND METHODS Five cohort studies’ participants (n = 8,270) were classified as statin users when they were on statin therapy at the time of DNAm assessment with Illumina 450K or EPIC array or noncurrent users otherwise. Associations of DNAm with various outcomes like incident type 2 diabetes, plasma glucose, insulin, and insulin resistance (HOMA of insulin resistance [HOMA-IR]) as well as with gene expression were investigated.
RESULTS Discovery (n = 6,820) and replication (n = 1,450) phases associated five DNAm sites with statin use: cg17901584 (1.12 × 10−25 [DHCR24]), cg10177197 (3.94 × 10−08 [DHCR24]), cg06500161 (2.67 × 10−23 [ABCG1]), cg27243685 (6.01 × 10−09 [ABCG1]), and cg05119988 (7.26 × 10−12 [SC4MOL]). Two sites were associated with at least one glycemic trait or type 2 diabetes. Higher cg06500161 methylation was associated with higher fasting glucose, insulin, HOMA-IR, and type 2 diabetes (odds ratio 1.34 [95% CI 1.22, 1.47]). Mediation analyses suggested that ABCG1 methylation partially mediates the effect of statins on high insulin and HOMA-IR. Gene expression analyses showed that statin exposure and ABCG1 methylation were associated with ABCG1 downregulation, suggesting epigenetic regulation of ABCG1 expression. Further, outcomes insulin and HOMA-IR were significantly associated with ABCG1 expression.
CONCLUSIONS This study sheds light on potential mechanisms linking statins with type 2 diabetes risk, providing evidence on DNAm partially mediating statins’ effects on insulin traits. Further efforts shall disentangle the molecular mechanisms through which statins may induce DNAm changes, potentially leading to ABCG1 epigenetic regulation.
Introduction
Statins effectively reduce the risk of cardiovascular disease (1). However, clinical trials and observational studies show that statins lead to insulin resistance and type 2 diabetes (2,3). The underlying mechanisms remain unclear.
Statins are associated with epigenetic changes, including histone acetylation, miRNA regulation (4), and DNA methylation (DNAm), particularly at genes related to lipid and insulin metabolism (5). DNAm is linked to type 2 diabetes pathophysiology (6); thus, it may be a potential mechanism contributing to the increased risk of type 2 diabetes observed in statin therapy. Nevertheless, this hypothesis has not been investigated (4).
We conducted an epigenome-wide association study (EWAS) in blood investigating the association between statin use and changes in DNAm at sites in the genome called CpGs. Further, we sought to replicate the findings and study the associations of the statin-related CpGs with gene expression, glycemic traits, and type 2 diabetes. Finally, we explored the potential role of the statin-related CpGs mediating the association of statins with glycemic traits and type 2 diabetes.
Research Design and Methods
Study Design and Population
Overall, 8,270 participants from five prospective cohort studies were included. In EWAS analyses, the discovery panel included 6,820 Caucasian individuals from Epidemiologische Studie zu Chancen der Verhütung, Früherkennung und optimierten Therapie chronischer Erkrankungen in der älteren Bevölkerung (ESTHER) (7), a population-based study recruited in Saarland, Germany, with participants aged 50–75 years; Kooperative Gesundheitsforschung in der Region Augsburg-F4 (KORA-F4) (8), a follow-up study of the KORA-S4 cohort, established in 1996 in Augsburg, Germany; Study of Health in Pomerania (SHIP)-Trend (9), a population-based study of participants from northeast Germany; and a population-based cohort of South Asians residing in London, U.K., aged 35–75 years: the London Life Sciences Prospective Population (LOLIPOP) study (10). Our study included 998 individuals from ESTHER, 1,726 from KORA-F4, 247 from SHIP-Trend, and 3,848 from LOLIPOP. Hence, the discovery panel was composed of 56% South Asians and 44% Caucasians.
The replication panel was selected based on ethnic homogeneity and availability of data for post-EWAS analyses. It included 1,450 Caucasians from the Rotterdam Study (RS) (subcohorts RS-III-1 [first visit of the third RS cohort] and RS-BIOS [comprising the third visit of the second cohort (RS-II-3) and the second visit of the third RS cohort (RS-III-2)]) (11), a prospective population-based cohort study of participants living in Ommoord, Rotterdam, the Netherlands, aged 45 years and older. Only RS had information on former statin users (n = 74). Additional association analyses between DNAm and glycemic traits (n = 1,165) and gene expression (n = 731) included RS data, whereas association analyses with incident type 2 diabetes included RS, KORA, and LOLIPOP (n = 4,422). See Fig. 1 for study design.
Study design. An EWAS was conducted to identify DNAm changes at CpG sites related to current statin use compared with noncurrent statin use. Four cohorts were used as a discovery panel, and findings were replicated in an independent cohort. Association analyses of the identified statin-related CpGs with glycemic traits, incident type 2 diabetes (T2D), and gene expression were performed. The large white arrows show the proceeding steps in the study design. The dotted black arrows represent a potential indirect effect of statins through changes in cg06500161 methylation. The black thin arrows indicate the direct effect of statins on glycemic traits (insulin and HOMA-IR).
Statin Use
Information on medication was obtained from pharmacy records in RS; self-declaration in KORA-F4, LOLIPOP, and SHIP-Trend; and general practitioner records in ESTHER.
Statin use status was categorized in current users and noncurrent users. Users were defined as current users were if the statins prescription occurred at the time of blood drawing for DNAm assessment or as noncurrent users otherwise. RS additionally identified former statin users as participants who had previously used statins but were no longer users on the blood draw date.
DNAm Data
DNA was extracted from whole peripheral blood (stored in EDTA tubes) by standardized salting out methods. Genome-wide DNAm was measured in bisulfite-converted genomic DNA using the Illumina Infinium HumanMethylation450 BeadChip or Illumina Infinium MethylationEPIC BeadChip according to the manufacturer’s protocols. Samples showing inadequate hybridization, incomplete bisulfite treatment, and gender swaps were excluded. Details on cohort-specific methods are provided in Supplementary Table 1. CpG methylation proportion was reported as a normalized β-value ranging from 0 to 1, where 1 represents 100% methylation. For avoidance of sex bias, CpGs annotated to genes in sex chromosomes were excluded. After normalization and quality control, overlapping CpGs across studies were meta-analyzed. CpGs missing in one cohort were also included, while CpGs missing in two or more cohorts were excluded, leaving 332,357 CpGs for final analysis.
mRNA Expression Data
mRNA data included a subset of the replication panel (n = 731 [RS-III-1]). Total RNA was isolated (PAXgene Blood RNA kits; QIAGEN) from whole blood (PAXgene Tubes; Becton Dickinson). Samples were analyzed using the LabChip GX (Caliper Life Sciences) according to the manufacturer’s instructions. Samples with an RNA quality score >7 were amplified, labeled (TotalPrep RNA; Ambion), and hybridized to the Illumina HumanHT-12 v4 Expression BeadChips according to the manufacturer’s protocol. Expression data were quantile normalized to the median distribution and log (2) transformed. Probe and sample means were centered to zero. Genes were considered significantly expressed when detection P values calculated by Genome Studio were <0.05 in >10% of all discovery samples. Thus, 21,238 probes were available to study. Quality control used eQTL (expression quantitative trait loci) mapping pipeline (https://github.com/molgenis/systemsgenetics/tree/master/eqtl-mapping-pipeline) (12). Only probes that uniquely mapped to the human genome build 37 and represented gene mRNA expression were analyzed (13). Samples were processed at the Genetic Laboratory of Internal Medicine, Erasmus Medical Center. Expression data are available for 881 samples at the Gene Expression Omnibus (GEO) public repository (accession no. GSE338828).
Covariates, Glycemic Traits, and Type 2 Diabetes
Covariates were selected based on previous literature on DNAm studies, including age, sex, BMI, smoking status, history of coronary heart disease (CHD), antihypertensive medication, systolic blood pressure (SBP), total cholesterol, HDL, triglycerides, and LDL. BMI was calculated by dividing weight in kilograms by the square of height in meters. Smoking history was self-reported and classified as current, former, and never smokers. SBP was the mean of two consecutive measurements at the right brachial artery with a random-zero sphygmomanometer with the participant in a sitting position. Antihypertensive medication use was defined as the use of diuretics, β-blockers, ACE inhibitors, and calcium channel blockers. CHD was defined as the occurrence of myocardial infarction, revascularization, coronary artery bypass graft surgery, or percutaneous coronary intervention. Triglycerides, HDL, and total cholesterol concentrations were measured in blood using an automated enzymatic method.
Glycemic traits examined were plasma fasting glucose, insulin, HOMA of insulin resistance (HOMA-IR), and type 2 diabetes incidence. Type 2 diabetes cases were defined as the presence of either serum glucose level ≥7.0 mmol/L or glucose-lowering medication use in RS; self-report in KORA; physician diagnosis or use of glucose-lowering medication in ESTHER; diabetes medication use, HbA1c ≥6.5%, serum glucose level (nonfasting) ≥11.1 mmol/L, or self-report in SHIP-Trend; and diabetes history or fasting glucose ≥7.0 mmol/L or HbA1c >6.5% in LOLIPOP. The interassay coefficients of variations are insulin <8%, glucose <1.4%, and lipids <2.1%. LDL was calculated using the Friedewald formula (total cholesterol − HDL − triglycerides/5).
Statistical Analyses
EWAS
The association between CpG methylation and current versus noncurrent statin users was investigated in a two-step EWAS approach: discovery and replication. Normalized DNAm β-values at each CpG were the dependent variable, and current statin use (yes = 1, no = 0) was the predictor. Per individual CpG, participants with methylation levels higher than three times the interquartiles range were excluded. The regression analysis used the lme4 package in R, version 3.4.2 (https://cran.r-project.org/web/packages/lme4/index.html), to perform linear mixed effects models. The random part of the model included adjustment for array number and position on array as random terms to account for batch effects in the DNAm measurement. The fixed part of the model included leukocyte proportions (to account for cell mixture), sex, age, smoking status, BMI, SBP, antihypertensive medication, and prevalent CHD and type 2 diabetes. Leukocyte proportions (B cells, CD4 T cells, CD8 T cells, granulocytes, monocytes, and NK cells) were estimated as previously described by Houseman (14), using the R minfi package (15). Statistical significance for discovery was determined by a Bonferroni-corrected threshold (1.0 × 10−7). Heterogeneity (I2) of effect estimates was assessed to account for differences between cohorts. Results across the discovery cohorts were combined in a meta-analysis by inverse variance–weighted method in METAL, version 2011-03-25.
The identified CpGs were replicated in an independent panel using identical models. Bonferroni correction for replication was 0.05 divided by the number of discovery findings. Results of discovery and replication panels were combined in a meta-analysis as aforementioned.
Pearson test was used to assess the correlation among the replicated CpGs in the replication panel (RS).
Associations Between Statin-Related CpGs and Glycemic Traits and Incident Type 2 Diabetes
Standardized methylation values at the replicated CpGs as predictors of fasting glucose, insulin, and HOMA-IR (RS [n = 1,165]) were tested in linear mixed effects models. For exploration of reverse causation, prevalent type 2 diabetes case subjects and former statin users were excluded. Logistic regression was used to assess the association with incident type 2 diabetes (RS, KORA, and LOLIPOP [n = 4,422, including 1,173 cases]). Results across cohorts were meta-analyzed as aforementioned. The adjusting covariates included leukocyte proportions, batch effects, sex, age, smoking status, SBP, antihypertensive medication, CHD, and statin use. Bonferroni correction was applied (<0.0025). Variables with right-skewed distributions were transformed to the natural logarithmic scale (glucose, insulin, HOMA-IR, and HDL).
Causal Mediation Analysis
To investigate whether DNAm changes are part of the pathway though which statins exert their diabetogenic effect, we explored the significant associations in the previous stage in a nonparametric causal mediation analysis using the R mediation package (16) as previously described (17). The mediation analysis dissects the total effect of the exposure (statin treatment) on the outcome (glycemic traits and incident type 2 diabetes) into the 1) indirect effect acting via the mediator of interest (DNAm) and 2) direct effect acting directly or via a mediator other than the one under study. The mediation package compares the effect estimates of each association between statin use and traits in subjects with different methylation levels. Nonparametric bootstrap with 20,000 resamples was performed. Data from the RS were used (n = 1,180), with exclusion of prevalent type 2 diabetes case subjects and former statin users. The normalized CpG methylation β-value was modeled as mediator and the CpG-associated trait as outcome. Supplementary Fig. 1 shows a schematic representation. Sex, age, and cohort were considered confounders. For investigation of potential unmeasured confounding, the sequential ignorability assumption (no confounder assumption) was tested in a sensitivity analysis using the correlated residuals method (18), which estimates the correlation that quantifies how strong the confounder would have to be to change the conclusion.
Gene Expression-Association Analysis
For investigation of possible biological pathways, relations between the identified CpGs and transcriptome-wide gene expression were explored (RS) (n = 731). Statin-related CpGs were regressed against methylation batch variables as random terms; a similar approach was applied for 21,238 transcription probes and the expression batch variables. Next, linear regressions between the residuals from the CpGs regression (predictors) and the expression probe residuals (dependent variable) were performed. The Bonferroni correction was 0.05/number of probes × number of CpGs tested.
Additionally, the significant expression probes in the aforementioned analysis were tested in association with statin use exposure (n = 631) and as predictors of the significant glycemic traits in the mediation analyses (n = 616). Type 2 diabetes case subjects and former statin users were excluded.
Dose Effects
In 303 current statin users from SHIP-Trend and RS with available data on statin dose, we examined the association between the replicated CpGs and the defined daily dose (DDD) (WHO ATC/DDD classification [https://www.whocc.no/atc_ddd_index/]), using the EWAS model. Next, the association between gene expression and statin dose was studied in the RS.
Confounding and Ethnicity
To investigate whether the associations of statin use with the replicated CpGs were independent of blood lipids, we additionally adjusted the EWAS for serum total cholesterol or for individual blood lipids instead (HDL, LDL, and triglycerides) in both discovery and replication.
In an effort to rule out confounding by indication, we restricted the EWAS to subjects with LDL ≥70 mg/dL or ≥100 mg/dL—cutoffs used to advise statin prescription (19,20). Further, we reran the EWAS excluding prevalent type 2 diabetes and prediabetes case subjects.
To account for ethnicity heterogeneity in the discovery panel (South Asians and Caucasians), we reran the EWAS, taking Caucasians (n = 4,349) for discovery and South Asians (n = 3,849) for replication.
Post Hoc Power Analysis
To calculate the range of power for replication of the findings, we used the GPower 3.1 tool (21) with six predictors representing all discovered CpGs except for cg03467813, which was not included due to its high heterogeneity (I2 = 99%). CpGs’ smallest β-estimate (0.0039 [cg10177197]) and probability error 0.008 or the largest estimate (0.0219 [cg17901584]) was used.
Results
Association Between Statin Use and DNA Methylation
The mean (SD) age in the discovery panel ranged from 51.5 (13.8) years (SHIP-Trend) to 62.1 (6.5) years (KORA-F4). Of the total discovery panel, 2,969 (43.5%) participants were women. Supplementary Table 2A and B shows other baseline characteristics.
In the discovery panel, seven CpGs passed the Bonferroni threshold for significance (P < 1.0 × 10−7), being differentially methylated in current statins users compared with noncurrent users, of which current users had lower methylation at four CpGs, annotated to genes DHCR24 (cg17901584, located at TSS1500), FAM50B (cg03467813, located at the gene body region), SC4MOL (cg05119988, located at 5′ untranslated region [5′UTR]) and AHRR (cg05575921, located at gene body). Three CpGs annotated to ABCG1 (cg06500161 and cg27243685, located at gene body) and DHCR24 (cg10177197, located at gene body) showed higher methylation among current users (Table 1 and Supplementary Fig. 1). All associations were independently replicated, except for two CpGs annotated to AHRR and FAM50B, which failed to reach statistical significance after Bonferroni correction for replication (P < 7.14 × 10−3). The post hoc power analysis showed a power ranging from 39% to 99% for our replication study. In the replication panel, correlations among CpGs ranged from r = −0.011 between cg06500161 (ABCG1) and cg17901584 (DHCR24) to r = 0.477 for cg06500161 and cg27243685 (both annotated to ABCG1) (Supplementary Table 3).
DNAm sites associated with current statin use (compared with noncurrent use)
The meta-analysis across the discovery and replication panels revealed one new CpGs at DHCR24 (cg17475467). The heterogeneity in the associations between statin use and the replicated CpGs across cohorts varied from I2 = 0 for cg10177197 (DHCR24) to I2 = 74.8 for cg17901584 (DHCR24) (Table 1). Unreplicated CpGs were not further investigated.
The sensitivity analysis using RS data with never statin users as the reference group (former users excluded) revealed similar results, replicating the five CpGs (Supplementary Table 4).
Associations of the Statin-Related CpGs With Glycemic Traits and Type 2 Diabetes
After a Bonferroni correction (2.5 × 10−3) and comprehensive assessment of potential confounders, increase in 1 SD of methylation at cg06500161 (ABCG1) was associated with incident type 2 diabetes (odds ratio 1.34 [95% CI 1.22, 1.47]), fasting glucose (1.03 × 10−3), insulin (4.63 × 10−8), and HOMA-IR (1.05 × 10−8). Increase in 1 SD of methylation at cg27243685 (ABCG1) was associated with incident type 2 diabetes (1.24 [1.13, 1.38]) and approaching borderline statistical significance with augmented insulin (4.62 × 10−3) but not statistically significantly associated with glucose or HOMA-IR. Increase in 1 SD of methylation at cg17901584 was approaching borderline statistical significance with incident type 2 diabetes (6.2 × 10−3) but was not significantly associated with glucose, insulin, or HOMA-IR. No significant findings were observed for cg05119988 (SC4MOL) or cg10177197 (DHCR24) (Table 2).
Association of DNAm of the replicated CpGs associated with statin use with glycemic traits and type 2 diabetes
Causal Mediation
Significant results were obtained for the models with outcome fasting insulin and HOMA-IR. For both cases, the analyses showed that 1) current statin use has a significant overall effect of β = 0.275 (0.190, 0.360) on insulin and β = 0.291 (0.199, 0.380) on HOMA-IR and 2) part of that effect goes directly or via mediator(s) other than cg06500161, called average direct effect (ADE). We found an ADE of β = 0.233 (0.149, 0.320) on insulin and β = 0.242 (0.151, 0.330) on HOMA-IR. 3) Another part of such effect may operate via an indirect path (indirect effect), possibly through cg06500161 methylation, with an average causal mediator effect of β = 0.043 (0.024, 0.060) on insulin and β = 0.049 (0.028, 0.070) on HOMA-IR; 4) consequently, of the total effect of statins on insulin and HOMA-IR, 15.5% (0.086, 0.260) and 16.8% (0.095, 0.270) are suggested to act via cg06500161 methylation, respectively. According to these results, cg06500161 methylation does not explain 100% of the effect of statins on insulin/HOMA-IR; thus, we may call cg06500161 methylation a partial mediator instead of a total mediator.
No significant findings were observed for the mediation model with outcome glucose (Table 3). We also investigated the potential mediation of cg27243685 on incident type 2 diabetes, with no significant findings. The sensitivity analysis showed that residual correlation due to unmeasured confounding of r ≥ 0.2 would be needed to violate the sequential ignorability assumption.
Causal mediation analysis on the significant associations between the statin-related CpGs (cg06500161 and cg27243685) as mediators and glycemic traits and incident type 2 diabetes as outcomes
Further mediation models were performed to account for potential confounders. In addition to age, sex, and cohort, we further adjusted for 1) BMI, 2) BMI and HDL, and 3) smoking status, SBP, antihypertensive medication, and presence of CHD. Results did not materially change with these further adjustments, although the effect sizes attenuated (Supplementary Tables 5–7).
Gene Expression
For exploration of potential epigenetic regulations of gene expression by DNAm, the identified CpGs (predictors) were tested in association with transcriptome-wide gene expression in the RS. cg06500161 was inversely associated with two expression probes at ABCG1 (ILMN_1794782 and ILMN_2329927) after Bonferroni correction (4.71 × 10−7). An inverse association was observed for cg27243685 and one probe (ILMN_1794782). These results indicate a lower expression of ABCG1 with increasing cg06500161 and cg27243685 methylation. Furthermore, cg17901584 methylation was directly associated with ABCG1 expression (ILMN_1794782) and suggestively associated with another probe at the same gene (ILMN_2329927). There were no significant findings for the other two CpGs after Bonferroni correction (Supplementary Table 8). Moreover, exposure to statin treatment was associated with lower ABCG1 expression (Supplementary Table 9). When ABCG1 expression was studied as a predictor of insulin traits, we found an inverse association between ILMN_1794782 and fasting insulin and HOMA-IR (Supplementary Table 10).
Dose Effect
Statin dose exposure was significantly associated with the five identified CpGs after Bonferroni correction (P < 0.01) (Supplementary Table 11). Exposure to increasing statin dose was associated nominally with lower levels of ABCG1 expression at ILMN_1794782 and suggestively with ILMN_2329927 (Supplementary Table 12).
Sensitivity Analyses
Residual confounding and confounding by indication were as follows: 1) additional adjustment of the EWAS for serum total cholesterol or blood lipids HDL, LDL, and triglycerides instead (Supplementary Table 13); 2) exclusion of prevalent type 2 diabetes and pre–type 2 diabetes case subjects (Supplementary Table 14); or 3) restriction of the EWAS to participants with LDL ≥70 mg/dL or ≥100 mg/dL (Supplementary Table 15) did not change the associations of the five CpGs. Further, exclusion of type 2 diabetes case subjects from the other analyses did not affect the results (data not shown).
With regard to transethnic replication, the new EWAS taking Caucasians only for discovery and South Asians for replication showed that all five CpGs passed the Bonferroni correction threshold (Supplementary Table 16).
Conclusions
The current study sheds light on potential mechanisms linking statin use and type 2 risk diabetes by, firstly, identifying and replicating associations between statin therapy and methylation at five CpGs (cg06500161, cg27243685, cg17901584, cg10177197, and cg05119988) independent of blood lipids and, secondly, providing evidence on the partial mediation via ABCG1 methylation in the effect of statins on increased levels of insulin and insulin resistance under the sequential ignorability assumption.
Our results on the association between statin use and DNAm agree with a report from the Framingham Study, which found epigenome-wide significance at only two sites (cg17901584 and cg06500161) (5). The Framingham Study examined this association in a smaller sample (n = 1,545 vs. 6,820 for our discovery panel), and contrary to our investigation, they did not proceed with replication or examine associations among gene expression, metabolic markers, and type 2 diabetes. Moreover, our work reports novel associations and adds to current knowledge that statin dosage might be implicated in the degree of methylation and ABCG1 expression. The latter finding goes in line with an experimental study where macrophages treated with various types of statins showed lower ABCG1 expression as dose increased (22).
Sites cg06500161 and cg27243685 are annotated to ABCG1, cg17901584 and cg10177197 to DHCR24, and site cg05119988 is annotated to SC4MOL. ABCG1 (ATP-binding cassette member-1 subfamily-G) encodes a protein that mediates the transport of different molecules, such as cholesterol efflux to the HDL, oxysterols, and phospholipid transport in macrophages (23,24). It is also involved in insulin secretion and sensitivity (25). ABCG1 expression was reduced in statin users compared with nonusers (26) and in patients with type 2 diabetes (27,28). Genes DHCR24 (24-dehydrocholesterol reductase) and SC4MOL (sterol-C4-methyl oxidase-like) code enzymes catalyzing different steps during cholesterol biosynthesis. DHCR24 mutations are related to desmosterolosis (29,30) and Alzheimer disease (31). Deficiency of the protein coded by SC4MOL produces congenital cataracts, microcephaly, growth delay, skin conditions, and immune dysfunction (32). An observational study found that DHCR24 and SC4MOL were upregulated among statin users (26), while SC4MOL upregulation increased type 2 diabetes risk (27,28).
Methylation at four of the identified CpGs (cg06500161, cg27243685, cg17901584, and cg05119988) was previously associated with glycemic traits, type 2 diabetes, and blood lipids, while cg06500161 (ABCG1) was associated with increased glucose, insulin, HbA1c, HDL, and triglyceride levels (33–35). Our causal mediation analyses provided evidence of cg06500161 methylation partially mediating the association between statin use and higher fasting insulin and HOMA-IR. However, given the cross-sectional nature of these analyses, our results must be interpreted with caution. Since the mediator was measured at the same time as the outcome, reverse causation cannot be completely ruled out. However, in an effort to overcome this issue, we excluded type 2 diabetes case subjects. The mediation effect observed may be the result of more complex interactions of combined effects across the epigenome. Moreover, the lack of significant findings in mediation by the other CpG at ABCG1 (cg27243685) could reflect a power issue. Hence, future studies shall explore how the CpGs within the same gene interact and mediate the effects.
Based on our findings and the available evidence, we hypothesize that cg06500161 hypermethylation may be a consequence of statin use, possibly inducing a decrease in ABCG1 transcription in blood. Further research is needed to investigate to what extent this downregulation could in turn compromise downstream signals, resulting in impaired insulin metabolism. In this line, impaired insulin sensitivity and secretion as a consequence of statin treatment were recently observed in a longitudinal study (36). Furthermore, a functional study suggested an epigenetic regulation of ABCG1 mediated by methylation-dependent transcription factor binding (37). Our hypothesized model is displayed in Fig. 2. Nevertheless, further investigations should test this hypothesis and assess the effect of statins in target tissues and to what extent these findings are generalizable to populations with ethnic backgrounds other than Caucasians and South Asians.
Scheme of the hypothesized mechanism linking statin therapy and risk of type 2 diabetes. The solid black arrow represents the effect of statins on plasmatic insulin and HOMA-IR levels that operates directly or through a pathway different from the mediator analyzed in the current study (methylation at cg06500161). The dotted black arrows represent the suggested alternative pathway, where an indirect effect of statins on insulin and HOMA-IR is mediated by cg06500161 methylation. The black thin arrows indicate decreased or increased level of the parameter; in this figure they are representing decrease in gene expression of ABCG1 and increase in serum levels of fasting insulin and HOMA-IR.
We highlight the large sample size with DNAm data in our study, the use of a replication panel, and the comprehensive assessment of confounding factors. Addition of complementary data like gene expression, transethnic replication, and use of a causal inference method as mediation analysis also add value to our research.
One limitation is the lack of data on DNAm and gene expression from specific tissues of interest for type 2 diabetes and drug metabolism (pancreas, liver, adipose tissue), as both DNAm and gene expression may be tissue specific (38). Instead, we used whole peripheral blood, allowing for cell mixture. Consequently, we included adjustments for estimated blood leukocytes proportions.
Most cohorts had limited information on former statin users; thus, it is likely that they were included in the reference group. However, the replication panel showed that the same results are valid for the group of current statin users only, although we had limited numbers to investigate such associations. Furthermore, it was not possible to explore potential effects of cessation time or treatment duration on DNAm. A Mendelian randomization approach based on genetic instrumental variables to explore a potential causal effect of statins on DNAm was not possible due to the lack of adequate genetic variants that may mimic the randomization of the hypothesized adverse effect of statins on type 2 diabetes risk. The use of variants at the gene encoding 3-hydroxy-3-methylglutaryl CoA reductase (39), statins’ target protein, would mimic the main effect of statins, rather than an alternative pathway leading to the hypothesized effect. Moreover, these variants were associated with LDL (39), reflecting cholesterol metabolism. This violates the Mendelian randomization assumption that the instrumental variable is associated with the outcome only through the exposure, since high cholesterol is a risk factor for type 2 diabetes.
All cohorts used Illumina 450K array to measure DNAm, except one, which used Illumina EPIC array (n = 247). Despite differences between the two, EPIC array covers ∼90% of the same sites on the 450K chip; both use the same technology and showed overall per-sample correlations of r > 0.99 (40), thus encouraging the application of combined data. Moreover, cohorts used different methods to diagnose type 2 diabetes, which might introduce error measurement.
The heterogeneity observed in the meta-analysis for some of the CpGs might be a consequence of differences in statin use prevalence (ranged from 9.2% to 18.1%) and types and dose of statin used across cohorts. Additionally, differences in ethnicity could have influenced heterogeneity (the association of statin use and cg06500161 in the overall meta-analysis showed I2 = 52.30, while I2 = 5.20 across Caucasians only), as could different methods applied to collect data on medication and disease (pharmacy and general practitioner records or self-report).
Finally, some analyses were performed using data from the replication panel only, which may introduce a type II error.
Conclusion
We identified and replicated five DNAm sites associated with statin therapy, of which cg06500161 may be partially mediating the effect of statins on increasing insulin levels and HOMA-IR. This could be one potential mechanism linking statin therapy and type 2 diabetes onset observed in clinical trials and observational studies. Nevertheless, other biological pathways should not be discarded. From our study, it is not clear what might be the mechanism by which statins may induce DNAm changes; this shall be a challenge for future research.
Meanwhile, close monitoring of risk factors for type 2 diabetes in these patients is imperative to prevent further complications.
Article Information
Funding. C.O.-R. reports receiving grant support from Agencia Nacional de Investigación y Desarrollo (ANID)/Scholarship Program/Doctorado Becas Chile/2016 - 72170524. P.P.M. and T.L. reported receiving financial support from the Academy of Finland [grants 286284, 134309 (Eye), 126925, 121584, 124282, 129378 (Salve), 117787 (Gendi), and 41071 (Skidi)], the Social Insurance Institution of Finland, Competitive State Research Financing of the Expert Responsibility area of Kuopio, Tampere and Turku University Hospitals (grant X51001), Juho Vainio Foundation, Paavo Nurmi Foundation, Finnish Foundation for Cardiovascular Research, Finnish Cultural Foundation, Tampere Tuberculosis Foundation, Emil Aaltonen Foundation, Yrjö Jahnsson Foundation, Signe and Ane Gyllenberg Foundation, Diabetes Research Foundation of Finnish Diabetes Association, and EU Horizon 2020 (grant 755320 for TAXINOMISIS Project). H.J.G., S.B.F., and A.T. report receiving financial support from the Federal Ministry of Education and Research (grants 01ZZ9603, 01ZZ0103, and 01ZZ0403), the Ministry of Cultural Affairs as well as the Social Ministry of the Federal State of Mecklenburg-West Pomerania, and the network “Greifswald Approach to Individualized Medicine (GANI_MED)” funded by the Federal Ministry of Education and Research (grant 03IS2061A) and the DZHK (grant 81X34000104). J.C. is supported by the Singapore Ministry of Health’s National Medical Research Council under its Singapore Translational Research Investigator (STaR) Award (NMRC/STaR/0028/2017). KORA was initiated and financed by Helmholtz Zentrum München (German Research Center for Environmental Health), which is funded by the German Federal Ministry of Education and Research (BMBF) and by the State of Bavaria. Furthermore, KORA research has been supported within the Munich Center of Health Sciences (MC-Health), Ludwig-Maximilians-Universität, as part of LMUinnovativ. This work was supported by a grant (WA 4081/1-1) from the German Research Foundation.
Duality of Interest. H.J.G. has received travel grants and speakers honoraria from Fresenius Medical Care, Neuraxpharm, Servier, and Janssen Cilag, as well as research funding from Fresenius Medical Care. O.H.F. reports receiving grants or research support from Metagenics, Inc., and Nestlé Nutrition (Nestec Ltd.). No other potential conflicts of interest relevant to this article were reported.
Author Contributions. T.M. and C.O.-R. contributed to the study concept and design. C.O.-R., R.W., B.L., P.P.M., X.G., M.L., D.J.-Q., J.B.J.v.M., W.Z., J.S.K., H.J.G., S.B.F., B.S., Y.Z., C.G., M.M.-N., M.H., A.P., and A.T. contributed to data generation and analyses. C.O.-R., T.M., R.W., B.L., P.P.M., X.G., D.J.-Q., M.L., E.P.-F., J.N., O.H.F., M.G., O.L.R.-O., A.T., H.B., M.W., M.A.I., J.B.J.v.M., T.V., J.C., and B.H.S. contributed to drafting of the manuscript. O.H.F., M.G., O.L.R.-O., T.L., A.T., H.B., M.W., M.A.I., J.B.J.v.M., T.V., J.C., B.H.S., and T.M. contributed to critical revision of the manuscript for important intellectual content. All authors approved the final version of the manuscript. C.O.-R. 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.
Prior Presentation. Parts of this study were presented in abstract form at the 78th Scientific Sessions of the American Diabetes Association, Orlando, FL, 22–26 June 2018, and at EuroPrevent 2019, Lisbon, Portugal, 11–13 April 2019.
Footnotes
This article contains Supplementary Data online at https://care.diabetesjournals.org/lookup/suppl/doi:10.2337/dc19-1828/-/DC1.
- Received September 13, 2019.
- Accepted January 14, 2020.
- © 2020 by the American Diabetes Association.
Readers may use this article as long as the work is properly cited, the use is educational and not for profit, and the work is not altered. More information is available at https://www.diabetesjournals.org/content/license.
References
