Accelerated DNA methylation age and medication use among African Americans

Participant characteristics

Phase 1 characteristics for GENOA participants are described in Table 1. Mean DNAm age estimates (ranging from 44.2 to 54.3 years) were lower than chronological age (mean = 57.0 years). DNAm age acceleration measures are also reported in Table 1 with PhenoAge acceleration (PhenoAA) being the highest (mean = 0.38 years) followed by HannumAge acceleration (HannumAA, mean = 0.15 years), HorvathAge acceleration (HorvathAA, mean = 0.14 years), and GrimAge acceleration (GrimAA, mean = 0.11 years). In the subset of the sample used for longitudinal analyses, the participants at Phase 2 were about 5.4 years older than at Phase 1 (mean age = 54.0 years at Phase 1 and 59.4 years at Phase 2).

Supplementary Table 1 shows medication use categorized by the first two digits of the Medi-Span Therapeutic Classification (MTC) code. Although diuretics, calcium channel blockers, and beta blockers are often prescribed as antihypertensive medications, they have different MTC codes than antihypertensives because they are used for other conditions as well. Drug categories with N < 30 were excluded. For the three major drug categories (antihypertensives, antihyperlipidemics, and diabetes medications), sub-drug categories were also included. The number of participants taking each drug was highest for diuretics (N=377) followed by antihypertensives (N=321), calcium channel blockers (N=221), and diabetes medications (N=164). Supplementary Table 2 shows the p-value of the associations between medication categories. Diuretic medication use was associated with most of the other medication categories, especially with the use of other antihypertensive medications (significant P range: 4.7 x 10^-11 – 0.03). Also, the use of renin-angiotensin-aldosterone system (RAAS) inhibitors was associated with diabetes medication (P< 1 x 10^-10). Both non-narcotic analgesic use and NSAID analgesic use were also significantly associated with several other medication use categories (significant P range: 9.6 x 10^-10 – 0.03).

Correlation among DNAm age estimators

The correlations between chronological age and DNAm age estimated by the four epigenetic clocks are presented in Supplementary Figure 1. As expected, chronological age was positively and significantly correlated with HorvathAge (Pearson correlation coefficient r = 0.86), HannumAge (r = 0.90), PhenoAge (r = 0.82), and GrimAge (r = 0.85). The Pearson correlation matrix of DNAm age acceleration, both before and after adjustment for white blood cell proportions, is shown in Supplementary Figure 2. The correlations range from 0.18 (GrimAA and HorvathAA) to 0.58 (HorvathAA and HannumAA). GrimAA was the least correlated with other AA measures.

Associations between medication use and DNAm age acceleration

We first fit a single multivariable model with all drug categories (N=15) to examine the unique contribution of each medication to each DNAm age acceleration measure using cross-sectional data from Phase 1 (N=1,100, Table 2). We were interested in associations that were significant at a nominal level (P<0.05) as well as at a Bonferroni-corrected significance threshold (P<.05/15 = 0.003). Model 1 adjusted for age and sex, and Model 2 also included body mass index (BMI), smoking, education, and alcohol consumption. Model 3 additionally adjusted for white blood cell proportions. Calcium channel blockers were independently associated with a 1.16-year lower HannumAA at the Bonferroni-adjusted significance threshold after adjusting for sex, BMI, smoking, education, and alcohol consumption (P=0.001, Model 2). The association remained significant after further adjustment for white blood cell proportions (Model 3).

Sensitivity analyses were performed to investigate whether the associations between medication use and DNAm age acceleration were independent of comorbidities. First, we repeated the multivariable models with all drug categories in the subgroup of participants with hypertension (N=779). The results from this analysis were not substantively different from the analysis in the full sample (Supplementary Table 3). The association between the use of calcium channel blockers and HannumAA was significant at the Bonferroni-corrected significance level in all models (Models 1-3). We further adjusted the multivariable models with all drug categories for the history of hypertension, stroke, coronary heart disease, and diabetes, as well as lipid levels to control for potential confounding effects (Table 3, Model 4). HannumAA was 1.27 years lower for those taking calcium channel blockers (P=0.001), and the association remained significant after further adjusting for white blood cell proportions at the Bonferroni-corrected significance level (P=0.0003 in Model 5).

Interaction between medication use and sex on DNAm age acceleration

We assessed whether the observed effect of medication use varied by sex by adding a medication-by-sex interaction term to the linear mixed models with one medication at a time (Supplementary Table 4). At a Bonferroni-corrected significance threshold (P=0.05/15=0.003), we observed a significant interaction between the use of RAAS inhibitors and HorvathAA after adjusting for age, sex, education, smoking, alcohol consumption, BMI, and white blood cell proportions (P=0.0015). In the sex-stratified analysis, RAAS inhibitor use was associated with 2.67-year higher HorvathAA in males after adjusting for all covariates and white blood cell proportions (P=0.001), but this association was not significant in females (Figure 1). We also included all other drug categories in the model to adjust for potential confounding effects, and the interaction remained significant (P=0.0009).

Figure 1. The plot shows the effects of RAAS inhibitors on HorvathAge acceleration in males (green) vs. females (blue). Each dot represents the least square mean of HorvathAge acceleration, and the bars represent the corresponding confidence intervals for the effect estimates. RAAS: Renin-angiotensin-aldosterone system.

Association between statin use and grimage components

Statin use was nominally associated with higher GrimAA in Models 1, 2, and 4. Since GrimAge is comprised of DNAm surrogates for seven proteins and smoking pack-years, we further assessed the association of statin use with each GrimAge component to identify the one(s) that may be driving the observed association. As shown in Supplementary Table 5, statin use was positively associated with DNAmADM, DNAmCystatinC, and DNAmTIMP1 after controlling for age, sex, education, smoking, alcohol consumption, and BMI at a nominal P-value threshold of 0.05. No associations were significant after Bonferroni correction for the number of GrimAge components (N=8). All nominally significant associations were attenuated after further adjusting for white blood cell proportions, suggesting that the effect of statins on GrimAge components may be partly due to changes in blood cell composition.

Association between medication use and longitudinal change in DNAm age acceleration

We examined whether medication use influenced the change in DNAm age acceleration from Phase 1 to Phase 2 in a subset of participants whose methylation measurements were collected at both Phases (N=266) using multivariable models. We only included drug categories with N ≥ 30 at both Phases 1 and 2 (diuretics, calcium channel blockers, antihypertensives, diabetes medications, and NSAID analgesics). We evaluated the longitudinal association between medication use and change in DNAm age acceleration after categorizing medication use as “never used,” “continued,” “stopped,” and “started.”

The associations that were significant after Bonferroni correction in the longitudinal models for the change in medication use between Phase 1 and Phase 2 are shown in Table 4. The full results are presented in Supplementary Table 6. Compared to those who never used antihypertensives, those who started to take them after Phase 1 had a 0.97-year decrease in GrimAA (P=0.006 in Model 6, Table 4), and the association remained significant after further adjusting for hypertension, stroke, coronary heart disease, diabetes, and lipid levels at the Bonferroni corrected significance level (P=0.007 in Model 7). Compared to those who never used NSAID analgesics, those who started to take them after Phase 1 had an increase in HorvathAA of about 2.61 years (P=0.0005 in Model 7). Both associations remained significant with further adjustment for white blood cell proportions (Model 8).

No other significant longitudinal effects were observed at the Bonferroni corrected significance threshold (Supplementary Table 6).

Study sample

The GENOA study is a community-based longitudinal study that was initiated to identify the genetic effects of hypertension and related target organ damage [57]. In the first phase of GENOA (Phase 1: 1996-2001), European American and African American sibships with at least 2 individuals who were clinically diagnosed with hypertension before age 60 were recruited. All other siblings were invited to participate, regardless of hypertension status. Exclusion criteria included secondary hypertension, alcoholism or drug abuse, pregnancy, insulin-dependent diabetes mellitus, active malignancy, or serum creatinine levels >2.5mg/dL [57]. In Phase 1, 1,583 European Americans (Rochester, MN) and 1,854 African Americans (Jackson, MS) were enrolled. In the first follow-up phase of GENOA (Phase 2: 2000-2004), 1,239 non-Hispanic white and 1,482 African American participants were successfully followed up, and their potential target organ damage from hypertension was measured. Demographics, medical history, clinical characteristics, information on medication use, and blood samples were obtained in each phase. In GENOA, DNA methylation levels were measured only in African American participants using blood samples collected at Phases 1 and 2. African American GENOA participants with DNA methylation available at Phase 1 were included in the current cross-sectional analysis (N=1,100). In the longitudinal analysis, those with available DNA methylation measurements at both Phases 1 and 2 and without changes in smoking status between phases were included (N=266). Written informed consent was obtained from all subjects, and approval was granted by participating institutional review boards (University of Michigan, University of Mississippi Medical Center, and Mayo Clinic).

Methylations measures

A total of 1,106 samples at Phase 1 and 304 samples at Phase 2 were assessed using the Illumina HumanMethylationEPIC BeadChip. First, raw IDAT files were imported using the Minfi R package [58]. We used the shinyMethyl R package [59] to visualize the raw intensity data and identify sex mismatches and outliers, which were removed. We also obtained the detection p-value for each sample at each probe, and individual probes with detection p-value <10^-16 were considered to be successfully detected [60]. Samples and probes with the detection rate <10% were removed [60]. Samples with incomplete bisulfite conversion identified using the QCinfo() function in the ENmix R package were removed [61]. We also checked sample identity using the 59 SNP probes implemented in the EPIC chip and removed mismatched samples. Next, Noob was used for individual background and dye-bias normalization [62]. Since two types of probes are present on the EPIC BeadChip (Infinium I and Infinium II), we used the Regression on Correlated Probes (RCP) method to adjust for probe-type bias [63]. ComBat, an empirical Bayes batch-correction method, was used sequentially to remove those batch effects (Johnson, 2007). We performed principal variance component analysis (PVCA) to quantify the variance explained by the known batch variables before and after batch adjustment to ensure that no single batch factors explained more than 3% of the variance. After exclusions, a total of 857,121 probes in 1,100 samples at Phase 1 and 294 samples at Phase 2 remained for analysis.

DNAm age calculation and blood cell counts

To calculate HorvathAge, methylation beta values were uploaded to the online Horvath epigenetic age calculator [18]. HorvathAge was estimated based on 353 CpGs [18]. This epigenetic clock captures chronological age across multiple tissue types and in different age groups. HorvathAA is the residual from regressing HorvathAge on chronological age. HannumAge was estimated based on 71 CpG sites using a single tissue (whole blood) [17], which captures chronological age. HannumAge acceleration (HannumAA) is the residual from regressing HannumAge on chronological age. PhenoAge was estimated using 513 CpG sites, and the PhenoAA measure is the residual from regressing PhenoAge on chronological age [19]. PhenoAge captures the effect of phenotypic aging as it was trained on many aging-related clinical measurements, including albumin, creatinine, glucose, C-reactive protein, alkaline phosphatase, white blood cell count, chronological age, and other metrics. GrimAge was estimated using 1,030 CpG sites, and GrimAA is the residual from regressing GrimAge on chronological age [20]. GrimAge captures the effects of the number of cigarettes smoked in pack-years as well as seven plasma proteins associated with mortality, including adrenomedullin, beta-2 microglobulin, growth differentiation factor 15, Cystatin C, leptin, plasminogen activation inhibitor 1, and tissue inhibitor metalloproteinase 1 [20]. Since GrimAge is a composite measure, we wanted to better understand the effects of pharmaceuticals on each component of GrimAge, so we additionally obtained the DNAm surrogates of all eight individual items from the online calculator. It is important to note that the EPIC array is missing 19 sites used to construct HorvathAge, and six sites that are used to calculate HannumAge; however, previous studies have demonstrated that this does not impact the performance of DNAm age predictors and acceleration measurements [64, 65]. All analyses were repeated with additional adjustment for white blood cell proportions, including CD8+ T, CD4+ T, natural killer, B cells, and granulocytes. These blood cell counts were estimated using Houseman’s method as implemented in the online epigenetic age calculator [66]. For each measure of DNAm age acceleration, outliers greater or less than 5 standard deviations from the mean were removed.

Drug classifications and other covariates

Prescription information for medications used in the past month was collected from the labels on prescriptions provided by the participants at each GENOA visit. The Medi-Span Therapeutic Classification (MTC) was utilized to codify the pharmacological class of each prescription using the first six digits. For the current analysis, we categorized medications into major classes using the first two digits of the MTC code. Major drug categories with N < 30 at Phase 1 were excluded from the analysis. The 12 major medication categories included diuretics, calcium channel blockers, beta blockers, antihypertensives, antihyperlipidemic, diabetes medication, non-narcotic analgesics, NSAID analgesics, antidepressants, antihistamines, antianxiety medication, and narcotic analgesics. Categories of antihypertensive, antihyperlipidemic, and diabetes medications were further subcategorized using the third digit of the MTC codes and were included in the analysis if N ≥ 30 at Phase 1: antihypertensives were further subcategorized into alpha blockers, sympatholytic, and RAAS inhibitors; antihyperlipidemics included statins as a sub-category; and diabetes medications were further classified into sulfonylurea and insulin. Although diuretics, calcium channel blockers, and beta blockers are often prescribed as antihypertensive medications, they have different MTC codes because they are used for other conditions as well, such as congestive heart failure, angina, and arrhythmia. As a result, a total of 15 drug categories were available at Phase 1 for the cross-sectional analysis. We only included major drug categories with N ≥ 30 at both Phases 1 and 2 (diuretics, calcium channel blockers, antihypertensives, diabetes medications, and NSAID analgesics) for the longitudinal analysis due to the limited sample size.

Self-reported smoking status was obtained at Phases 1 and 2 of GENOA. Current smokers were defined as individuals who reported smoking at the time of the exam and had smoked at least 100 cigarettes in one’s lifetime. Former smokers were individuals who had smoked at least 100 cigarettes in one’s lifetime but do not currently smoke. Never smokers were individuals who had not smoked more than 100 cigarettes in one’s lifetime. Educational attainment was collected at GENOA Phase 1 as the self-reported highest degree, as well as years of education. Participants were grouped as having less than a high school degree (<12 years of education), high school degree or equivalent (12 years or GED), and some college and above (> 12 years). Alcohol consumption was determined based on the number of drinks per week. Because the distribution of alcohol consumption was skewed, it was log-transformed prior to the analysis. Weight and height were collected at both phases, and BMI was calculated. Participants’ medical conditions were self-reported and confirmed during the medical exam. Participants were classified as having hypertension if they were taking antihypertensive medications, had a systolic blood pressure ≥ 140 mmHg, or had a diastolic blood pressure ≥ 90 mmHg. Patients taking diabetes medication or with glucose level measurement ≥ 126mg/dL were categorized as having diabetes.

Statistical analyses

We used Chi-square tests to assess independence among medication categories. We then performed cross-sectional analysis to assess the associations between medication use and the four DNAm age acceleration measures at Phase 1. For each DNAm age acceleration measure, multivariable linear mixed models were constructed with all of the drug categories as the predictors and DNAm age acceleration as the outcome, adjusting for age and sex, with batch effects and familial relationships as random effects (Model 1). We included all of the drug categories in a single model to capture the unique contribution of each medication. We also included socioeconomic and lifestyle factors that have been previously associated with DNAm age acceleration, including education, smoking status, alcohol consumption, and BMI (Model 2). Since the DNAm age clocks are correlated with white blood cell counts, we additionally adjusted for white blood cell proportions to assess whether any of the associations were due to differences in blood cell composition (Model 3). For the three major drug categories (antihypertensives, antihyperlipidemics, and diabetes medications) that had sub-categories, we included only the sub-categories into the multivariable model. Since we evaluated multiple medications, Bonferroni correction was used to account for multiple comparisons. However, the measures of medication use are not independent, and the Bonferroni approach is conservative in this setting, so we also considered a nominal significance level (P<0.05).

Since almost everyone treated with a drug has, by definition, the relevant disease to begin with, it is important to assess whether the associations of medication use with DNAm age acceleration were independent of relevant participants’ disease status or traits. To address this issue, we ran two sensitivity analyses, as follows. First, we repeated the analyses with multivariable models among a subset of participants who had hypertension since GENOA is a predominately hypertensive cohort. Second, we repeated the analyses in all participants controlling for age-related diseases and lipid levels (hypertension, stroke, coronary heart disease, diabetes, high-density lipoproteins (HDL), triglycerides, and low-density lipoproteins (LDL) (Model 4). Model 5 additionally adjusted for white blood cell proportions.

We assessed whether the observed effects of medication use on DNAm age acceleration differed by gender. For medication-by-sex interactions significant after Bonferroni correction (P<0.05/15 = 0.003), we further examined the sex-specific effects of the medication on DNAm age acceleration. Since GrimAge was comprised of DNAm surrogates of eight individual items (ADM, B2M, GDF15, CystatinC, Leptin, PAI1, TIMP1, and PACKYRS), we also assessed the observed effect with each component to identify those that may drive the association for any medication categories that were associated with GrimAA in Model 4.

Last, we conducted longitudinal analysis for drug categories with N ≥ 30 at both Phase 1 and Phase 2 on a subset of samples (N=266) whose DNA methylation was measured at both Phases. Batch effects and familial relationships were adjusted as random effects. Because smoking status substantially influences DNA methylation, we removed 27 participants whose smoking status changed between Phases. We considered the change in medication use between Phases by setting the change in DNAm age acceleration between Phases as the outcome and the change in medication use as the predictor. The change in medication use was categorized as "Never used" (reference category), "Continued," "Stopped," and "Started." Model 6 covariates included Phase 1 age, sex, education, BMI, smoking status, alcohol consumption, and DNAm age acceleration, and Model 7 additionally controlled for Phase 1 hypertension, stroke, coronary heart disease, diabetes, HDL, triglycerides, and LDL. We repeated the analyses with further adjustment for white blood cell counts to assess whether the association observed was independent of the changes in blood cell composition (Model 8).