Investigation of bidirectional longitudinal associations between advanced epigenetic age and peripheral biomarkers of inflammation and metabolic syndrome

Cross-lagged models: Hannum DNAm age residuals

The MetS CFA fit the data well at both T1 and T2 with all indicators loading significantly on their respective latent variables at the p ≤ 0.003 level (details available from corresponding author). At both time points, the Lipid/Obesity factor showed the strongest loading on the higher-order MetS factor (βs = 0.92 – 0.96) followed by the Blood Pressure factor (βs = 0.48 – 0.53), and the Blood Sugars factor (βs = 0.23 – 0.48). The cross-lagged panel analysis examining bidirectional longitudinal associations between DNAm age residuals and MetS severity revealed significant autoregressive effects between each variable and itself over time as well as a significant cross-lagged association. Advanced DNAm age at T1 predicted increases in MetS severity at T2 (standardized β = 0.17, p < 0.001), accounting for baseline levels of MetS (Figure 1A). T1 MetS and T1 DNAm age residuals were correlated with each other (r = 0.20, p = 0.006). Notably, the competing cross-lagged path (from T1 MetS to T2 DNAm age residuals) was not significant (standardized β = -0.05, p = 0.32). We examined potential confounds of the T1 DNAm age to worsening T2 MetS association and found that this association remained significant when additionally controlling for demographic factors (race, education), psychiatric factors (cigarette use, major depression, alcohol abuse/dependence), medication use (including psychotropic and metabolic-related medications), and time between assessments (see Supplementary Materials). Of note, in these follow-up analyses, major depressive diagnoses at T1 also predicted increasing MetS severity over time (standardized β = 0.28, p = 0.035). To further ensure no influence of additional PCs on the reported results, we also ran secondary analyses investigating the potential effects of all 20 estimated ancestry PCs (see Supplementary Materials).

Figure 1. The Figure shows the results of cross-lagged models examining longitudinal associations between Hannum DNAm age residuals and metabolic syndrome (MetS) severity factor scores (A), and Lipids/Obesity factor scores (B). Measures of each marker were residualized on age and sex (applicable to A and B). (***p < 0.005, **p < 0.01, **p < 0.05).

There were no other cross-lagged effects between any of the other peripheral biomarkers and Hannum DNAm age residuals (see Figure S1), however associations between T1 DNAm age residuals and T2 CRP (standardized β = 0.13, p = 0.053) and T2 WBCs (standardized β = 0.10, p = 0.068) just missed the threshold for statistical significance (Figure S1 A, B). When we re-analyzed these associations in follow-up models using DNAm age residuals that did not account for estimated WBCs (i.e., only age, sex, and the first two PCs were regressed out of DNAm age estimates), we found that results for measured WBCs were unchanged, but results for CRP were nominally significant, though they missed our a priori corrected p-value threshold (β = 0.15, p = 0.031; Figure S4). A significant association between T1 PTSD symptom severity and worsening CD4/CD8 t-cell profiles (e.g., controlling for baseline CD4/CD8) emerged in that model (standardized β = -0.11, p = 0.013; Figure S1C).

Given that the MetS factor scores were derived from a higher-order CFA, we wondered if the association between advanced DNAm age at T1 and increasing MetS at T2 would be evident in each of the lower-order factors that comprise MetS. Therefore, we conducted additional cross-lagged models in which each of the lower-order metabolic factor scores (Lipids/Obesity, Blood Sugars, Blood Pressure) was included in the model in place of MetS, again residualized for age and sex. Results revealed associations only between T1 DNAm age residuals and increasing Lipid/Obesity factor scores at T2 (standardized β = 0.16, p = 0.001), controlling for the significant baseline effects of the Lipid/Obesity factor scores (standardized β = 0.75, p < 0.001; Figure 1B; see Supplementary materials for results pertaining to individual indicators of this latent variable). A T1 correlation between advanced DNAm age and Lipid/Obesity factor scores was also evident (r = 0.21, p = 0.003). No cross-lagged effects emerged with the Blood Pressure or Blood Sugar factor scores (Figure S2).

Cross-lagged models: Horvath DNAm age residuals

We found no significant cross-lagged associations between Horvath DNAm age residuals and any of the peripheral biomarkers evaluated (Figure S3). There were also no significant concurrent correlations between Horvath DNAm age residuals and any of the peripheral biomarkers (Figure S3). As in the model with the Hannum DNAm age residuals, PTSD severity was associated with worsening CD4/CD8 T-cell profiles (Figure S3).

Study limitations

Results from this study should be interpreted with several limitations in mind. First, the study cohort was composed primarily of white male veterans. Future work is needed to establish that these results generalize to populations with more diverse compositions of ancestry and sex. Second, other biological variables that were not interrogated here may also play an important role. For example, unmeasured third variables (e.g. physical health diagnoses) could account for the predictive effects of either DNAm age residuals or PTSD on biological changes at T2. That said, this was a young adult cohort and individuals with neurological diseases and diabetes were excluded, which attenuates this concern. It remains to be seen if DNAm age residuals at T1 play a causative or etiological role in predicting negative health outcomes at T2, or if they are simply a marker for an underlying biological process. Third, though longitudinal, our study only covered a two-year period; our results should be interpreted with this in mind, as some associations of epigenetic aging and biological markers may not be detectable across this relatively short period. However, the unreliability in DNAm age estimates would be expected to be more impactful when examining change in DNAm age estimates over a relatively short compared to long period of time. Furthermore, DNAm age for individuals in the sample increased, on average, approximately one year for every chronological year, providing greater confidence that changes in DNAm age and DNAm age residuals over time are meaningful and not a function of error. Finally, our study was limited to analysis of DNAm and biological parameters in blood; we did not investigate tissue-specific cellular aging in the brain or other organs.

Participants

Participants were previously described in Wolf et al. [19]. They were U.S. military veterans (post-9/11 conflicts) who enrolled in the Translational Research Center for TBI and Stress Disorders (TRACTS) longitudinal study at the US Department of Veterans Affairs (VA) Rehabilitation Research and Development Traumatic Brain Injury Center of Excellence at VA Boston Healthcare System. The TRACTS longitudinal study has been described in detail previously [43]. In brief, it is an ongoing research protocol evaluating traumatic stress, traumatic brain injury (TBI), health, and neural and cognitive factors among returning veterans. Exclusion criteria for the study included the following: history of seizures unrelated to head injury, severe or unstable psychological diagnosis preventing participation, acute psychotic or bipolar disorder, neurological illness, acute homicidal and/or suicidal ideation with intent to act, and cognitive disorder due to general medical condition not related to TBI. As in Wolf et al. [19], the cohort investigated in this study was based on a subset of 179 TRACTS participants with DNAm data from two time points, Time 1 (T1) and Time 2 (T2), at the time of the second data freeze when DNA was processed. The clinical and sociodemographic characteristics of the study cohort are shown in Table 1.

Procedure

Participants provided written informed consent, and then completed a comprehensive interview and self-report-based psychological assessment. All diagnostic interviews were administered by doctoral-level psychology professionals. A team of psychologists reviewed each interview to determine consensus ratings of presence or absence of psychological diagnoses. For each time point (T1 and T2), blood was drawn for DNA extraction and metabolic assays. T1 and T2 assessments were conducted an average of 1.89 years apart (Table 1). The study was approved by the VA Boston Healthcare System IRB. All T1 and T2 samples were processed using the Illumina EPIC chip.

Measures

The Clinician Administered PTSD Scale for DSM-IV (CAPS [44],), a well-validated diagnostic interview, was used to assess PTSD status and symptom severity. The CAPS was administered by doctoral-level psychologists. Additional information regarding the administration and rating of interviews is provided in the Supplementary Materials. In this manuscript, our analyses focused on a dimensional index of current PTSD symptom severity at T1. Additional measures that were included in supplementary analyses are described in the Supplementary Materials.

DNA extraction, genotyping, and ancestry-based principal components analysis

Full details on genotyping protocols, techniques, and data cleaning procedures are detailed in Logue et al. [45], and are also summarized here. DNA extraction was performed using a Qiagen AutoPure instrument with Qiagen reagents. DNA concentrations were normalized using the Quant-iT^TM PicoGreen dsDNA fluorescent assay (Invitrogen). To determine DNA quality and quantity, TaqMan RNase P Detection assay was used (Applied Biosystems Assay, Life Technologies, Carlsbad, CA) with fluorescence detection on a 7900 Fast Real Time PCR Instrument (Applied Biosystems, Life Technologies, Carlsbad, CA). DNA was then whole-genome amplified, fragmented, precipitated, resuspended, and was then hybridized on Illumina HumanOmni2.5-8 beadchips for 20 hours at 48°C according to manufacturer’s instructions (Illumina, San Diego, CA), followed by a single-base extension and multi-layered staining process. Beadchips were imaged using the Illumina iScan System, and results were processed with the Illumina GenomeStudio v2011.1 software and the Genotyping v1.9.4 module. Genotypes were then used to develop principal components (PCs) to model ancestry: PCs were determined using 100,000 randomly chosen common (minor allele frequency >5%) single nucleotide polymorphisms (SNPs) in PLINK version 1.9 [46]; PCs were used as ancestry covariates in subsequent analyses.

Methylation

DNA was hybridized to the Infinium MethylationEPIC BeadChip per manufacturer’s instructions. T1 and T2 samples were run together on the same chip, balancing the presence of cases and controls across chips and chip positions in order to reduce systematic bias. We utilized the processing pipeline established by the Psychiatric Genetics Consortium (PGC) PTSD Epigenetics workgroup [47] for the Illumina HumanMethylation450 BeadChip updated to apply to the EPIC chip. We used GenomeStudio to derive individual-level background-corrected probe data and idat files, and cleaned DNAm data using the CpGassoc package [48] and the ChAMP package [49] in R (R Development Core Team, 2008). Individual methylation values that did not meet a detection p < 0.001 were set to missing, and probes with >10% missing values were dropped. One chip had 7 out of 8 failed samples based on a criteria of >5% missing values; data for this chip were discarded and samples were rerun on a new chip. Subsequently no samples had >5% missing data, and all were retained for analysis. No samples had intensity <50% of the experiment-wide mean or with intensity <2,000 arbitrary units (AU). Cross hybridizing probes [50] and “underperforming” EPIC probes according to Illumina Product Quality Notification PQN0223 04/19/2017 were also excluded. R v. 3.1.0 was used for data cleaning.

Epigenetic age calculation

For Hannum DNAm age estimates, data were normalized using the beta mixture quantile dilation (BMIQ) method in the wateRmelon R package [51] as previously described [5,7,19] and batch correction was performed using an empirical Bayes method implemented in COMBAT [52]. Horvath DNAm age estimates were determined following an R script based on 335 probes assessed on the EPIC chip that passed quality control (QC). We have previously shown that DNAm age estimate correlations with chronological age are similar across the EPIC and 450K chips [7,53]. The association between Hannum DNAm age estimates derived from the EPIC chip in this cohort and chronological age was r = 0.88 (p < 0.001) and r = 0.85 (p < 0.001) at T1 and T2, respectively, and for the Horvath algorithm, the association was r = 0.90 (p < 0.001) and r = 0.91 (p < 0.001) at T1 and T2, respectively [7]. Hannum and Horvath DNAm age estimates were correlated with each other at both time points (T1, r = 0.88, p < 0.001; T2, r = 0.87, p < 0.001), as were Hannum and Horvath DNAm age residuals (r = 0.44, p < 0.001, for both T1 and T2) [19].

MetS (lipid/obesity, blood pressure, sugar levels)

Height, weight, and waist-to-hip ratio were measured along with two standing and two seated blood pressure readings (taken at 1-minute intervals). Blood samples were obtained, processed immediately upon collection, and shipped the same day to a commercial laboratory (Quest Diagnostics, Cambridge, MA). This laboratory assessed HDL cholesterol, triglycerides, and glucose (fasting glucose and glycated hemoglobin A1c levels). These metabolic measures were included in an overall index of Metabolic Syndrome (MetS) severity using confirmatory factor analysis (CFA) as described in the analysis section below. Total white blood cell counts were also measured via complete blood chemistry.

Estimated white blood cell count and CD4/CD8 ratio

Specific white blood cell (WBC) type proportions at T1 and T2 were not available from the Quest metrics and instead were estimated based on the methylation data. In brief, CD4 and CD8 T-cells, natural killer cells, b-cells, and monocytes were estimated based on the methylation data using the R minifi package [54] according to procedures described by Houseman et al. [55], Jaffe and Irizarry [56], and Fortin et al. [57]. CD4 and CD8 T-cell estimates were used to calculate a ratio of CD4 to CD8 T-cells, which has previously been shown to be an index of immunosenescence (see also: Holbrook et al. [58]).

C-reactive protein (CRP) serum levels

Serum CRP was assessed in blood samples as previously described in Miller et al. [23]. Serum CRP levels were measured in a commercial laboratory (Quest Diagnostics, Cambridge, MA) using a nephelometric assay with CRP monoclonal antibodies (analytical sensitivity = 0.10 mg/dL). Laboratory assay procedures were standardized to CRP reference preparations (International Federation of Clinical Chemistry and Laboratory Medicine/Bureau Communautaire de Reference/College of American Pathologists). The sample mean at T1 was 0.19 mg/dl (SD=1.13; range: 0.09-0.69 mg/dl), and at T2 was 0.29 mg/dl (SD=0.37; range: 0.09-2.34 mg/dl). Data were log-transformed for analysis (referred to as “CRP log” due to the distribution of raw CRP values being positively-skewed as in Miller et al. [23]. Two outliers were removed (one from each time point) as their CRP estimates were 12 SDs above the mean.

Data analysis

An overall index of MetS severity was calculated using CFA of raw laboratory values and physiologic measurements as previously described in Wolf et al. [20]. CFA is ideal for measuring the common metabolic factors that underlie the covariation of various biological assays because it models their relationship to a shared latent (or common) variable. Factor scores on the latent variable can then be generated to reflect the severity of the metabolic pathology for each subject. This approach avoids concerns about arbitrary diagnostic thresholds for metabolic disease and instead models the severity and comorbidity of the metabolic components dimensionally. As in Wolf et al. [20], we developed a higher-order CFA. The lower order factors represented: (a) Blood Pressure (indicated by two seated diastolic and systolic readings); (b) Lipid/Obesity (indicated by waist-to-hip ratio, body mass index (BMI), high density lipoprotein, and triglycerides); and (c) Blood Sugars (indicated by fasting glucose and glycated hemoglobin A1c levels). The above factors were specified to load together (i.e., to be accounted by) on a higher-order factor representing overall MetS severity. The model was tested separately in the T1 and T2 data and factor scores on all latent variables were saved for subsequent analyses. Of note, insulin was not available for the majority (94%) of the subjects at T1, as this test was added later to the protocol, and thus insulin levels were not included in the MetS CFA.

We conducted cross-lagged panel models (a form of path analysis) to simultaneously evaluate bidirectional longitudinal associations between advanced epigenetic age and each peripheral biomarker of interest (Figure 2). In this analysis, each variable measured at T2 is regressed on the same variable at T1 (i.e., the auto-regressive effect) and on the competing variable at T1 (i.e., the cross-lagged effect). For example, in the model examining MetS, T2 advanced epigenetic age (as defined by DNAm age residuals) was regressed on T1 DNAm age residuals and on T1 MetS factor scores while T2 MetS was simultaneously regressed on T1 MetS and T1 DNAm age residuals. The association between T1 PTSD severity and both T2 variables was also modeled. The concurrent correlations among the predictors at T1 and the residual correlation among dependent variables at T2 were also evaluated. A significant cross-lagged path would indicate, for example, that T1 Mets predicts changes in T2 advanced epigenetic age, controlling for T1 advanced epigenetic age. DNAm age residuals at each time point were generated by regressing raw DNAm age estimates on age, sex, estimated WBCs (CD4-T, CD8-T, NK, b cells, monocytes) from the respective time point, and the top two ancestry PCs and saving the unstandardized residuals from this equation. For analyses predicting estimated CD4/CD8 ratios, DNAm age residuals were calculated by regressing raw DNAm age estimates on age, sex, and the top two ancestry PCs (but not on estimated WBCs as these were the focus of this analysis). In a similar set of analyses, we also investigated measured WBCs and CRP phenotypes using Horvath and Hannum DNAm age residuals that did not take into account the estimated WBCs from DNAm age (regressing raw Horvath and Hannum DNAm estimates on age, sex, and the top two ancestry PCs, but not on estimated WBCs).

Figure 2. The Figure shows the cross-lagged model used to examine longitudinal associations between DNAm age residuals (Hannum or Horvath) and biological variables of interest (MetS, lab-based WBC measurement, CRP levels, CD4/CD8 T-cell ratio). Measures of each biological marker were residualized on age and sex for all analyses. DNAm age residuals at each time point were generated by regressing raw DNAm age estimates on age, sex, estimated WBCs (CD4-T, CD8-T, NK, b cells, monocytes) from the respective time point, and the top two ancestry PCs and saving the unstandardized residuals from this equation. For analyses predicting estimated CD4/CD8 ratios, DNAm age residuals were calculated by regressing raw DNAm age estimates on age, sex, and the top two ancestry PCs (but not on estimated WBCs as these were the focus of this analysis).

For the sake of simplicity and consistency with the DNAm age residuals, each peripheral biomarker was also first regressed on age and sex and the residuals from this equation saved for use in the path models. This approach controls for variance in these demographic covariates by removing them from both the T1 and T2 variables.

This cross-lagged panel approach was followed, in separate analyses, for each peripheral biomarker of interest (MetS, lab-based WBC totals, CRP, and estimated CD4/CD8 ratios) and for both the Horvath and Hannum algorithms. For analyses with significant cross-lagged associations, we conducted follow-up regressions which included a series of potential confounding variables to determine if they accounted for the significant cross-lagged associations, including demographic variables (education and self-reported racial/ethnic minority), psychological variables (cigarette use, current major depressive diagnosis, current alcohol abuse or dependence diagnosis), medication variables (current use of anti-hypertensives, cholesterol lowering medication, diabetes medication, antidepressants, anti-epileptics, sedatives/hypnotics, and pain medications; see Supplementary Materials), and time between assessments.

As each analysis was executed twice (once for Horvath and once for Hannum-based indices), we took into account the correlation between the two advanced epigenetic age metrics by adjusting for 1.8 tests and set the p-value threshold for statistical significance for individual parameters of interest (e.g., the association between DNAm age residuals and the T2 biological variable controlling for the same biological variable at T1) in a model at p = 0.028. This p-value correction was derived from a permutation testing procedure as described in Miller et al. [23] and Wolf et al. [6]; based on the r = 0.49 Horvarth/Hannum DNAm age residual association in the larger cross-sectional dataset, the adjusted p-value was found to represent 1.8 tests. There was no multiple testing correction across analyses for different biological variables of interest as these analyses investigate distinct hypotheses across different families of tests. All analyses were conducted with Mplus 8.0 statistical modeling software [59]. As all models were just identified (i.e., fully saturated such that df = 0), model fit will always be perfect and thus is not reported here.